! ***********************************************************************
!
!   Copyright (C) 2010  Bill Paxton
!
!   this file is part of mesa.
!
!   mesa is free software; you can redistribute it and/or modify
!   it under the terms of the gnu general library public license as published
!   by the free software foundation; either version 2 of the license, or
!   (at your option) any later version.
!
!   mesa is distributed in the hope that it will be useful, 
!   but without any warranty; without even the implied warranty of
!   merchantability or fitness for a particular purpose.  see the
!   gnu library general public license for more details.
!
!   you should have received a copy of the gnu library general public license
!   along with this software; if not, write to the free software
!   foundation, inc., 59 temple place, suite 330, boston, ma 02111-1307 usa
!
! ***********************************************************************
 !search debug to turn on/off heat 

      module run_star_extras

      use star_lib
      use star_def
      use const_def
      use crlibm_lib
      
      implicit none
      real(dp) :: accum_ener1, accum_ener2, accum_ener3, accum_ener4, accum_ener5, E_heat_accum
      real(dp) :: delta_ener1 ,delta_ener2 ,delta_ener3 ,delta_ener4 ,delta_ener5, E_heat_delta

      ! these routines are called by the standard run_star check_model
      contains
      
      !include 'standard_run_star_extras.inc'

      subroutine extras_controls(id, ierr)
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
       
         s% extras_startup => extras_startup
         s% extras_check_model => extras_check_model
         s% extras_finish_step => extras_finish_step
         s% extras_after_evolve => extras_after_evolve
         s% how_many_extra_history_columns => how_many_extra_history_columns
         s% data_for_extra_history_columns => data_for_extra_history_columns
         s% how_many_extra_profile_columns => how_many_extra_profile_columns
         s% data_for_extra_profile_columns => data_for_extra_profile_columns  
         s% job% warn_run_star_extras = .false.     


	 !New routines to add

	 s% other_energy => wave_energy
	 !s% other_torque => wave_torque


            
      end subroutine extras_controls




      subroutine wave_energy(id, ierr)
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         integer :: k,nz,k_max,j,ki,kshock,l
		 real(dp) :: conv_freq,heatfrac,L_heat,M_heat,vg,tprop
		 real(dp) :: conv_freq_max,L_waves_max,L_conv_max
		 real(dp) :: om,omsq,lambl1,lambl2,lambl3,dr,krl1,krl2,krl3,intkrl1,intkrl2,intkrl3,Tl1,Tl2,Tl3
		 real(dp) :: L_waves,L_max
		 real(dp) :: diff,diff_t,brunt,vel, gamphot
		 real(dp) :: f_damp_int,f_damp,gamnu,f_fac,x_fac
		 real(dp) :: maxbrunt,k_min,heatsmooth
		 real(dp) :: Tl1b,intkrl1b,conv_edge_flag,intkrl1min,Tl1min,minflag,intmin,krl1min,kmin,bruntmin
		 real(dp) :: M_damp_rad,M_damp_shock1,M_damp_shock2,M_damp_shock
		 real(dp) :: flagp,shockflag,shockcon,ri,rhoi,csi,Lmaxi,ui,kri,integrand,rshock,rhoshock
		 real(dp) :: krxir, krxir_max, max_nl_rad, lmax,shellrsize,scaleval, Lfracpt1a, Lfracpt2a, Lfrac_a, &
		 			Lfracpt1b, Lfracpt2b, Lfrac_b, Lfrac_int, Ldenom, Lfrac,Lheatl,Lfractestval
		 real(dp) :: temp, scaleheight!
		 integer :: i, n1, largest_shell, n2, swap_index,itemp
		 real(dp), allocatable :: conv_qs(:), ups_fin(:), lows_fin(:), shellwidth(:),temp1(:),temp2(:),temp4(:), &
		 							Lheatl_array(:), Lheatltest(:), delta_ener_array(:)
		 integer, allocatable :: conv_locs(:), ups(:), lows(:), temp_lows(:),upper_ms(:),lower_ms(:), &
		 						 temp_ups(:), shellindex(:),temp3(:),temp5(:)
		 logical, allocatable :: large_shells_mask(:)
		 real(dp), allocatable :: L_conv_max_array(:), conv_freq_max_array(:), L_waves_max_array(:), &
	  	                          om_array(:), omsq_array(:), heatfrac_array(:),L_waves_array(:), &
	  	                          lmax_array(:),scale_array(:) 
	  	 integer, allocatable :: k_max_array(:)


         ierr = 0
         call star_ptr(id,s,ierr)
         if (ierr /= 0) return

      conv_locs = PACK(s% conv_bdy_loc(:), s% conv_bdy_loc(:)/=0)

      ups=PACK(conv_locs,s% top_conv_bdy(1:SIZE(conv_locs)))
      lows=PACK(conv_locs,s% top_conv_bdy(1:SIZE(conv_locs)) .eqv. .false.)

      n1 = 1

      if (MOD(SIZE(s% top_conv_bdy(1:SIZE(conv_locs))),2) /= 0) n1=2

      allocate(upper_ms(SIZE(ups)),lower_ms(SIZE(ups)))
      upper_ms(:) = 0
      lower_ms(:) = 0

      if (n1==2) then 
      	allocate(temp_lows(SIZE(lows)+1))
      	temp_lows = (/s% nz, lows /)
      	deallocate(lows)
      	allocate(lows(SIZE(temp_lows)))
      	lows=temp_lows
      	deallocate(temp_lows)
      endif

      if (SIZE(upper_ms)>0) then
		upper_ms(1) = ups(1)
		lower_ms(1) = lows(1)
      endif
	  do i=2,SIZE(ups)

	  	if(((s% m(lows(i))/Msun-s% m(ups(i-1))/Msun) < 1.) .and. & 
	  		!in larger mass models, sometimes very far away shells have similar abundance
		    ((ALL(ABS(PACK(s% xa(:,lows(i)),s% xa(:,lows(i))>1d-3 .or. s% xa(:,ups(i-1))>1d-3 )- &
	  		& PACK(s% xa(:,ups(i-1)),s% xa(:,lows(i))>1d-3 .or. s% xa(:,ups(i-1))>1d-3 ))/MAX(1d-3, &
	  		& PACK(s% xa(:,lows(i)),s% xa(:,lows(i))>1d-3 .or. s% xa(:,ups(i-1))>1d-3 ), &
	  		& PACK(s% xa(:,ups(i-1)),s% xa(:,lows(i))>1d-3 .or. s% xa(:,ups(i-1))>1d-3 )) < 0.1)))) then!)) then
	  		upper_ms(i)=ups(i)
	  		lower_ms(i)=lower_ms(i-1)
	  	else
	  		upper_ms(i)=ups(i)
	  		lower_ms(i)=lows(i)
	  	endif
	  enddo 
	
	allocate(large_shells_mask(SIZE(upper_ms)))

	large_shells_mask(:) = .true.
	do i = SIZE(upper_ms)-1, 1, -1
		large_shells_mask(i) = .not.(ANY(lower_ms(i+1:)==lower_ms(i)))
	enddo

	ups_fin=s% m(PACK(upper_ms,large_shells_mask))/Msun
	lows_fin=s% m(PACK(lower_ms,large_shells_mask))/Msun

	allocate(shellwidth(SIZE(ups_fin)),shellindex(SIZE(ups_fin)),temp_ups(SIZE(ups_fin)),temp_lows(SIZE(ups_fin)))
	temp_ups=PACK(upper_ms,large_shells_mask)
	temp_lows=PACK(lower_ms,large_shells_mask)

	shellwidth=ups_fin-lows_fin
	shellindex=(/(i, i=1,SIZE(shellwidth), 1)/)


	!if (SIZE(shellwidth)>2) then
		do i=1, SIZE(shellwidth)-1
			largest_shell=MAXLOC(shellwidth(i:SIZE(shellwidth)),dim=1)
			
			swap_index=largest_shell+i-1
			if (swap_index /= i) then
				temp = shellwidth(i)
				shellwidth(i)=shellwidth(swap_index)
				shellwidth(swap_index)=temp
				itemp=shellindex(i)
				shellindex(i)=shellindex(swap_index)
				shellindex(swap_index)=itemp
			endif
		enddo
		allocate(temp4(SIZE(PACK(shellwidth,shellwidth>0.015*s% x_ctrl(1) )))) !0.5/17, 0.3/15 ~ 0.025*M
		temp4=PACK(shellindex,shellwidth>0.015*s% x_ctrl(1))
		deallocate(shellindex)
		allocate(shellindex(SIZE(temp4)))
		shellindex=temp4
		deallocate(temp4)

		allocate(temp1(SIZE(shellindex)),temp2(SIZE(shellindex)),temp3(SIZE(shellindex)),temp5(SIZE(shellindex)))
		temp1=0.
		temp2=0.
		temp3=0
		temp5=0
		do i=1,SIZE(shellindex)
			temp1(i)=ups_fin(shellindex(i))
			temp2(i)=lows_fin(shellindex(i))
			temp3(i)=temp_ups(shellindex(i))
			temp5(i)=temp_lows(shellindex(i))
		enddo
		ups_fin=temp1
		lows_fin=temp2 !these are now sorted from largest to smallest shell
		temp_ups=temp3
		temp_lows=temp5
	!endif

	shellindex=(/(i, i=1,SIZE(shellindex), 1)/)

	do i = 1, MIN(SIZE(ups_fin),5)-1
		!sort shells from largest radial coord to smallest
		largest_shell=MAXLOC(s% r(temp_ups(i:MIN(SIZE(ups_fin),5))),dim=1)
			
		swap_index=largest_shell+i-1
		if (swap_index /= i) then
			temp = temp_ups(i)
			temp_ups(i)=temp_ups(swap_index)
			temp_ups(swap_index)=temp

			temp = temp_lows(i)
			temp_lows(i)=temp_lows(swap_index)
			temp_lows(swap_index)=temp
			itemp=shellindex(i)
			shellindex(i)=shellindex(swap_index)
			shellindex(swap_index)=itemp
		endif
	enddo
	do i=1,SIZE(shellindex)
		temp1(i)=ups_fin(shellindex(i))
		temp2(i)=lows_fin(shellindex(i))
		!temp3(i)=temp_ups(shellindex(i))
	enddo
	ups_fin=temp1
	lows_fin=temp2 !these are now sorted from largest to smallest radial coord
	!temp_ups=temp3

  	!print*, 'ups_fin', ups_fin
	!print*, 'lows_fin', lows_fin


	  deallocate(upper_ms,lower_ms,large_shells_mask,temp1,temp2,temp3,temp5)

	  !!!!!!!!!!!!!!!!!!!!!!!!!!Calculate wave frequency and luminosity!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	  om = 1d-3
	  omsq = 1d-3
	  k_max = s% nz
	  
	  allocate(k_max_array(6),  conv_freq_max_array(6), L_waves_max_array(6), &
	  	& om_array(6), omsq_array(6), lmax_array(6), scale_array(6))

	  k_max_array(:)=0. !actually the location of excitation point 
	  conv_freq_max_array(:)=0.  !actually this is a mass weighted average
	  L_waves_max_array(:)=0. !actually this is a mass weighted average
	  lmax_array(:)=0.
	  scale_array(:)=0
	  om_array(:)=0. 
	  omsq_array(:)=0.

	  do i = 1, MIN(SIZE(ups_fin),6)
	  	  shellrsize=s% r(temp_ups(i))-s% r(temp_lows(i))
	  	  conv_freq_max= 0. !1d-99 
	      L_waves_max = 0. !1d-99
		  do k=1,s% nz
			if ( (s% m(k)/Msun > lows_fin(i)) .and. (s% m(k)/Msun < ups_fin(i)) ) then 
				L_waves_max = L_waves_max + s% dm(k) * (s% conv_vel(k)/s% csound(k))*s% L_conv(k)
				conv_freq_max = conv_freq_max + s% dm(k) * s% conv_vel(k)/(4.*s% scale_height(k))
			end if		
		  end do
		  k_max = temp_ups(i) !location of excitation is top of shell
		  conv_freq_max = conv_freq_max/((ups_fin(i)-lows_fin(i))*Msun) !mass-weighted average
		  L_waves_max = L_waves_max/((ups_fin(i)-lows_fin(i))*Msun) !mass weighted average
          scaleheight=s% P(k_max)/(s% rho(k_max) * s% grav(k_max))
		  scaleval=MIN(shellrsize, scaleheight) !^pressure scale height
		  lmax = s% r(k_max)/scaleval !maximum angular wavenumber: 
		  k_max_array(i)=k_max
		  conv_freq_max_array(i)=conv_freq_max
		  L_waves_max_array(i)=L_waves_max
		  lmax_array(i) = lmax
		  scale_array(i)=scaleval/Rsun 
		  om_array(i) = 2.*pi*conv_freq_max
		  omsq_array(i) = pow2(2.*pi*conv_freq_max)
	  end do

	  heatfrac = 1d-99 ! repeat all scalar quantities for shell
	  L_waves = 1d-99

	  allocate(Lheatl_array(6),Lheatltest(6),heatfrac_array(6),L_waves_array(6))
	  Lheatl_array(:) = 0.
	  Lheatltest(:) = 0.
	  L_waves_array(:)=0.
	  heatfrac_array(:)=0.
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Calculate transmission coefficient!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	  do i = 1, MIN(SIZE(ups_fin),6)
	  	  k_max=k_max_array(i)
	  	  omsq=omsq_array(i)
	  	  lmax=lmax_array(i)
	  	  L_waves=L_waves_max_array(i)

		  intmin = 0.
		  

		  Lheatl=0.
		  Lfractestval=0.
        do l=1, 10 !loop over l=1 to l=10
	      intkrl1min = 0.
		  do k=1,k_max !this is for all zones outside wave excitation point 

			minflag = 0.

	        ! vel = 0.
	        ! if (abs(s% v(k)) < 1d11) vel = abs(s% v(k)) !Required if v turned off
	        !print tau(k), vel
	        !if ( (s% tau(k) > 1d4) .and. (vel < s% csound(k)) ) then !Avoid outflowing surface regions

			if (s% brunt_N2(k) < 0.) brunt = 0. !Avoid artificial drops in Brunt
			if (s% brunt_N2(k) > 0.) brunt = s% brunt_N2(k) !Avoid artificial drops in Brunt
	        
			lambl1 = l*(l+1)*pow2(s% csound(k)/(s% r(k)))

			krl1 = (brunt - omsq)*(lambl1 - omsq)/(omsq*pow2(s% csound(k)))

			dr = s% dm(k)/(4.*pi*s% rho(k)*pow2(s% r(k)))

			!Find smallest transmission coefficient 
			if (krl1 < 0.) then
				minflag = 1.
				intmin = intmin + dr*sqrt(-krl1)
			
			end if
			if (minflag < 0.5) then
				if (intmin > intkrl1min) then 
					intkrl1min = intmin
				end if
				intmin = 0.
			end if
			!end if
		  end do !loop over k

		  Tl1min = exp_cr(-2.*intkrl1min) !it says l1 but it's for whatever l value in loop

		  tprop = 0.
		  f_damp_int = 0.
		  x_fac = 0.
		  f_fac = 0.
		  do k=k_max,1,-1
			lambl1 = l*(l+1)*pow2(s% csound(k)/(s% r(k)))
			if (s% brunt_N2(k) < 0.) brunt = 0. 
			if (s% brunt_N2(k) > 0.) brunt = s% brunt_N2(k) 
			krl1 = (brunt - omsq)*(lambl1 - omsq)/(omsq*pow2(s% csound(k)))
			dr = s% dm(k)/(4.*pi*s% rho(k)*pow2(s% r(k)))
			if (krl1 > 0.) then
				!Damp energy via radiative diffusion
				diff_t = (16d0*boltz_sigma*pow3(s% T(k)))/(3d0*pow2(s% rho(k))*s% cp(k)*s% opacity(k))
				gamphot = krl1*diff_t
				if (omsq > brunt) then
					vg = s% csound(k)
					tprop = tprop + dr/vg
				endif
				if (omsq < brunt) then
					vg = omsq*s% r(k)/sqrt(l*(l+1)*brunt) !l(l+1) = 2 for l=1
					tprop = tprop + dr/vg
					gamnu = -pow2(s% grav(k)*s% gamma1(k)*s% grada(k))
					gamnu = gamnu*(s% d_epsnuc_dlnT(k) - s% d_nonnucneu_dlnT(k))/(brunt*pow2(s% csound(k)*s% csound(k)))
					if (gamnu>0.) then
						f_damp_int = f_damp_int + 2.*gamnu*dr/vg
					end if
					if (gamphot>0.) then
						f_damp_int = f_damp_int + 2.*gamphot*dr/vg
					end if

				end if
			end if
		  end do
	  	  f_damp = exp_cr(f_damp_int)
		  !f_fac = 1.-1./f_damp
		  f_fac = f_damp-1
		  x_fac = f_damp_int
		  !heatfrac = (1.)/(1. + x_fac/Tl1min)
		  heatfrac = (1.)/(1. + f_fac/Tl1min)

		  Lfracpt1a = -2.*exp_cr(-pow2(DFLOAT(l))/pow2(lmax))* &
		  					(pow2(DFLOAT(l))+DFLOAT(l)*lmax+pow2(lmax))/pow2(lmax)
		  Lfracpt2a = sqrt(pi)*Erf(DFLOAT(l)/lmax)
		  Lfrac_a = (Lfracpt1a+Lfracpt2a)/4.

		  Lfracpt1b = -2.*exp_cr(-pow2(DFLOAT(l)-1)/pow2(lmax))* &
		  					(pow2(DFLOAT(l)-1)+(DFLOAT(l)-1)*lmax+pow2(lmax))/pow2(lmax)
		  Lfracpt2b = sqrt(pi)*Erf((DFLOAT(l)-1)/lmax)
		  Lfrac_b = (Lfracpt1b+Lfracpt2b)/4.

		  Lfrac_int = Lfrac_a - Lfrac_b
		  Ldenom= (2.+sqrt(pi))/4.

		  Lfrac = Lfrac_int/Ldenom
		  ! heatfrac=1./(1.+1./Tl1min)
		 
		  Lheatl = Lheatl + L_waves*Lfrac*heatfrac ! add up the heating rate over l=1 to 10 for this shell
		  Lfractestval=Lfractestval+Lfrac

		  if (l == 1) then
		  	!Tl1min_array(i)=Tl1min
		  	heatfrac_array(i)=heatfrac
		  	L_waves_array(i)=heatfrac*Lfrac*L_waves
		  end if

	    end do !loop over l

	    if (s% center_he4 < 1d-3) then !.or. s% x_ctrl(2) == 1.) then
	! 		L_waves_array(i) = L_waves
			Lheatltest(i) = Lfractestval
	    	Lheatl_array(i) = Lheatl
		endif

	  end do !loop over shells

	  allocate(temp1(6))
	  temp1(:) = 0.
	  temp1(1:MIN(SIZE(ups_fin),6))=lows_fin(1:MIN(SIZE(ups_fin),6))
	  lows_fin = temp1
	  deallocate(temp1)
	
	if (s% x_ctrl(2) == 1.)  then !to make it work for helium stars 
		allocate(temp1(SIZE(L_waves_max_array)+1))
		temp1(1)=0
		temp1(2:)=L_waves_max_array
		L_waves_max_array=temp1

		temp1(2:)=lmax_array
		lmax_array=temp1

		temp1(2:)=lows_fin
		lows_fin=temp1

		temp1(2:)=omsq_array
		omsq_array=temp1

		deallocate(temp1)

		allocate(temp3(SIZE(k_max_array)+1))
		temp3(1)=s% nz
		temp3(2:)=k_max_array
		k_max_array=temp3
		deallocate(temp3)
	endif

        !print*, 'k_max_array', k_max_array, ups_fin
	
		s% xtra1 = L_waves_max_array(2)/Lsun
		s% xtra2 = lmax_array(2)
		s% xtra3 = omsq_array(2)
		s% xtra4 = lows_fin(2)
		s% ixtra1=k_max_array(2)

		s% xtra5 = L_waves_max_array(3)/Lsun
		s% xtra6 = lmax_array(3)
		s% xtra7 = omsq_array(3)
	  	s% xtra8 = lows_fin(3)
	  	s% ixtra2 = k_max_array(3)

		s% xtra9 = L_waves_max_array(4)/Lsun
		s% xtra10 = lmax_array(4)
		s% xtra11 = omsq_array(4)
		s% xtra12 = lows_fin(4)
		s% ixtra3 =k_max_array(4)

		s% xtra13 = L_waves_max_array(5)/Lsun
		s% xtra14 = lmax_array(5)
		s% xtra15 = omsq_array(5)
		s% xtra16 = lows_fin(5)
		s% ixtra4 = k_max_array(5)
	
		s% xtra17 = L_waves_max_array(6)/Lsun
		s% xtra18 = lmax_array(6)
		s% xtra19 = omsq_array(6)
		s% xtra20 = lows_fin(6)
		s% ixtra5 = k_max_array(6)

		allocate(delta_ener_array(5))
		delta_ener_array(:)=0.
		if (SIZE(ups_fin)>1) then
			do i = 2, MIN(SIZE(ups_fin),6)
				if (Lheatl_array(i) > 1.e3*Lsun) then
						delta_ener_array(i-1) = s% dt * Lheatl_array(i)
				endif
			end do
		endif
		delta_ener1=delta_ener_array(1)
		delta_ener2=delta_ener_array(2)
		delta_ener3=delta_ener_array(3)
		delta_ener4=delta_ener_array(4)
		delta_ener5=delta_ener_array(5)

		s% extra_heat(:) = 1d-99

		deallocate(k_max_array, conv_freq_max_array, L_waves_max_array, lmax_array, temp_ups, temp_lows, &
		&  om_array, omsq_array, shellwidth,shellindex, heatfrac_array, scale_array, &  
		& L_waves_array, Lheatltest, Lheatl_array )

      end subroutine wave_energy









      integer function how_many_extra_profile_columns(id, id_extra)
         use star_def, only: star_info
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_profile_columns = 38
      end function how_many_extra_profile_columns
      
      
      subroutine data_for_extra_profile_columns(id, id_extra, n, nz, names, vals, ierr)
         use star_def, only: star_info, maxlen_profile_column_name
         use const_def, only: dp
         integer, intent(in) :: id, id_extra, n, nz
         character (len=maxlen_profile_column_name) :: names(n)
         real(dp) :: vals(nz,n)
         real(dp) :: om_con,diff,diff_t,diff_c,vg,epsheat,ebind,debind,brunts,lamb1,lamb2,lamb3,Lwavemax,Lwaveprop, vel
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         integer :: k
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         
         !note: do NOT add the extra names to profile_columns.list
         ! the profile_columns.list is only for the built-in profile column options.
         ! it must not include the new column names you are adding here.

         !if (n /= 1) stop 'data_for_extra_profile_columns'
         names(1) = 'L_waves'
	 names(2) = 'log_L_waves'
	 names(3) = 'conv_freq'
	 names(4) = 'log_conv_freq'
	 names(5) = 'M_damp'
	 names(6) = 'log_M_damp'
	 names(7) = 't_therm'
	 names(8) = 'log_t_therm'
	 names(9) = 't_diff'
	 names(10) = 'log_t_diff'
	 names(11) = 't_heat'
	 names(12) = 'log_t_heat'
	 names(13) = 't_cross'
	 names(14) = 'log_t_cross'
	 names(15) = 't_cool'
	 names(16) = 'log_t_cool'
	 names(17) = 't_dyn'
	 names(18) = 'log_t_dyn'	 
	 names(19) = 'log_L_conv'
	 names(20) = 'Mdot'
	 names(21) = 'Ebind'
	 names(22) = 'logEbind'
	 names(23) = 'log_N'
	 names(24) = 'log_L1'
	 names(25) = 'log_L2'
	 names(26) = 'log_L3'
	 names(27) = 'log_f_con'

         names(28) = 'L_maximum'
	 names(29) = 'log_L_max'
         names(30) = 'L_acmax'
	 names(31) = 'log_L_acmax'
	 names(32) = 'M_shock'
	 names(33) = 'log_M_shock'
	 names(34) = 'L_wave_prop'
	 names(35) = 'log_L_wave_prop'
	 names(36) = 't_adv'
	 names(37) = 'log_t_adv'
	 names(38) = 'diff_therm'

	 !names(28) = 'log_alpha_RTI'
	 !names(28) = 'dlnepsnucdlnT'
	 !names(29) = 'dlnnonnucdlnT'
	 !names(30) = 'depstot'


	 ebind = 0.
	 Lwavemax = s% xtra3!s% xtra1*s% xtra4
	 Lwaveprop = 0.

         vel = 0.

         do k = 1, nz

            if (abs(s% v(k)) < 1d11) vel = abs(s% v(k)) !Required if v turned off

	    if (s% brunt_N2(k) < 1d-14) brunts = 1d-14 
	    if (s% brunt_N2(k) > 1d-14) brunts = s% brunt_N2(k) 

	    !Wave generation power
            vals(k,1) = (s% conv_vel(k)/s% csound(k))*s% L_conv(k)/Lsun !Lwave calculation
	    vals(k,2) = safe_log10_cr(vals(k,1))

	    !Linear convective turnover frequency
	    vals(k,3) = (s% conv_vel(k)/(4.*s% scale_height(k)))
	    vals(k,4) = safe_log10_cr(vals(k,3))

	    !Calculate damping mass in solar masses
	    om_con = 2.*pi*s% xtra6
	    diff_t = 16.*boltz_sigma*pow3(s% T(k))/(3.*pow2(s% rho(k))*s% cp(k)*s% opacity(k))
	    diff_c = 0.!s% conv_vel(k)*s% scale_height(k)
	    diff = diff_t + diff_c
	    vals(k,5) = 4.*pi*s% rho(k)*pow2(s% rmid(k))*pow3(s% csound(k))/(pow2(om_con)*diff*Msun)
	    vals(k,6) = safe_log10_cr(vals(k,5))

	    !Thermal time scale in years
	    vals(k,7) = 4.*pi*s% rho(k)*pow2(s% rmid(k))*pow2(s% csound(k))*s% scale_height(k)/(s% L(k)*3.15d7)
	    vals(k,8) = safe_log10_cr(vals(k,7))

	    !Photon diffusion time scale in years
	    vals(k,9) = (s% r(1) - s% r(k))*s% tau(k)/(clight*3.15d7)
	    vals(k,10) = safe_log10_cr(vals(k,9))

	    !Wave heating time scale in years
	    vals(k,11) = pow2(s% csound(k))/(s% extra_heat(k)*3.15d7)
	    vals(k,12) = safe_log10_cr(vals(k,11))

	    !Gravity wave crossing time scale across a scale height, in years
	    vg = pow2(om_con)*s% r(k)/sqrt(2.*brunts)
	    vals(k,13) = s% scale_height(k)/(vg*3.15d7)
	    vals(k,14) = safe_log10_cr(vals(k,13))

	    !Neutrino cooling time scale, in years
	    epsheat = s% eps_nuc(k) - s% non_nuc_neu(k)
	    vals(k,15) = (s% cp(k)*s% T(k))/(-10.*epsheat*3.15d7)
	    vals(k,16) = safe_log10_cr(vals(k,15))

	    !Local dynamical time, in years
	    vals(k,17) = (s% scale_height(k))/(s% csound(k)*3.15d7)
	    vals(k,18) = safe_log10_cr(vals(k,17))

	    !Convective luminosity
	    vals(k,19) = safe_log10_cr(s% L_conv(k)/Lsun)

	    !Local mass outflow rate, in Msun/year
	    vals(k,20) = 4.*pi*s% rho(k)*pow2(s% r(k))*vel*(3.156d7/Msun)
!s% v(k)*(3.156d7/Msun)

	    !Calculate total binding energy of overlying envelope
	    debind = 6.67d-8*s% m(k)*s% dm(k)/s% r(k)
	    ebind = ebind + debind
	    vals(k,21) = ebind
	    vals(k,22) = safe_log10_cr(vals(k,21))

	    vals(k,23) = safe_log10_cr(sqrt(brunts)/(2.*pi))

	    !Calculate Lamb frequencies
	    lamb1 = 2.*pow2(s% csound(k)/(s% r(k)))
	    lamb2 = 6.*pow2(s% csound(k)/(s% r(k)))
	    lamb2 = 12.*pow2(s% csound(k)/(s% r(k)))

	    vals(k,24) = safe_log10_cr(sqrt(lamb1)/(2.*pi))
	    vals(k,25) = safe_log10_cr(sqrt(lamb2)/(2.*pi))
	    vals(k,26) = safe_log10_cr(sqrt(lamb3)/(2.*pi))

	    !Wave frequency
	    vals(k,27) = safe_log10_cr(om_con/(2.*pi))

	    !L_max
	    vals(k,28) = 2.*pi*s% rho(k)*pow2(s% r(k))*pow3(s% csound(k))
	    vals(k,29) = safe_log10_cr(vals(k,28)/Lsun)

	    !Wave action max, accounts for background flow
	    vals(k,30) = vals(k,28)*(1d0+vel/s% csound(k))
!+s% v(k)/s% csound(k))
	    vals(k,31) = safe_log10_cr(vals(k,30)/Lsun)

	    !Damping mass due to weak shock dissipation
	    vals(k,32) = (3d0*pi/(s% gamma1(k) + 1.))*(vals(k,28)/(om_con*pow2(s% csound(k))))/Msun
	    vals(k,32) = vals(k,32)*sqrt(vals(k,28)/Lwavemax)*(1d0+vel/s% csound(k))**(5./2.)
!s% v(k)/s% csound(k))**(5./2.)
	    vals(k,33) = safe_log10_cr(vals(k,32))

	    !Propagating wave flux
	    Lwaveprop = Lwaveprop + s% dm(k)*s% extra_heat(k)
	    vals(k,34) = Lwaveprop/Lsun
	    !if (k > s% xtra6) then
            !    Lwaveprop = 1d-99
	    !endif
	    vals(k,35) = safe_log10_cr(vals(k,34))

	    !Advective time scale
	    vals(k,36) = (s% scale_height(k)/vel)/3.156d7
!s% v(k))/3.156d7
	    vals(k,37) = safe_log10_cr(vals(k,36))
            
            !Thermal diffusivity
	    vals(k,38) = 16.*boltz_sigma*pow3(s% T(k))/(3.*pow2(s% rho(k))*s% cp(k)*s% opacity(k))

	    !vals(k,28) = safe_log10_cr(s% alpha_RTI(k))
	    !vals(k,28) = s% d_epsnuc_dlnT(k)/(s% eps_nuc(k))
	    !vals(k,29) = s% d_nonnucneu_dlnT(k)/(s% non_nuc_neu(k))
	    !vals(k,30) = s% d_epsnuc_dlnT(k) - s% d_nonnucneu_dlnT(k)

         end do

      end subroutine data_for_extra_profile_columns





      integer function how_many_extra_history_columns(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_history_columns = 165
      end function how_many_extra_history_columns
      
      
      subroutine data_for_extra_history_columns(id, id_extra, n, names, vals, ierr)
      	 use chem_def, only: ipp, icno, i3alf, i_burn_ne, i_burn_si, icc, ico, ioo
         integer, intent(in) :: id, id_extra, n
         character (len=maxlen_history_column_name) :: names(n)
         real(dp) :: vals(n)
         real(dp) :: lmax, krxir, krxir_max, max_nl_rad, omsq, heatfrac, L_waves, krl1, lambl1, xirl1, &
         			 krxir_maxofell, ell_maxnl,max_nl_radofell, ell_maxL
         real(dp) :: burn_h, burn_he, burn_c, burn_o, burn_ne, burn_si, dr, brunt, r, intmin, intkrl1min
         real(dp) :: Lheatl, dLheatl, tprop, f_damp, f_damp_int, x_fac,f_fac, gamnu, Ldenom, Lfrac_int
         real(dp) :: Lfrac, Lfrac_a, Lfrac_b, Lfracpt1a, Lfracpt2a, Lfracpt1b, Lfracpt2b, dLheatl_max
         real(dp) :: Tl1min, vg, gamphot, diff_t
         real(dp), allocatable :: krxir_max_array(:),max_nl_rad_array(:), omsq_array(:), L_waves_max_array(:), &
         						  heatfrac_array(:), L_waves_array(:), upper_ms(:), lower_ms(:),  &
         						  burntypes_array(:),  maxburn_array(:), lmax_array(:), Lheatl_array(:), & 
         						  ell_maxnl_array(:), dLheatl_max_array(:), ell_maxL_array(:)
         real(dp) :: dLheatl_array(5,10), fesc_array(5,10) 
         integer, allocatable :: k_max_array(:)
         integer :: k_max, length, i, k, maxburn, minflag, l, count
         character(2) :: lchar
         character(1) :: ichar
         character(9) :: fescname
         character(11) :: dLheatname
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         !allocate all the necessary arrays and make the arrays

  		 allocate(k_max_array(5),lower_ms(5),upper_ms(5),omsq_array(5),L_waves_max_array(5),lmax_array(5))
	 	 k_max_array=(/ s% ixtra1, s% ixtra2, s% ixtra3, s% ixtra4, s% ixtra5 /)
	 	 L_waves_max_array = (/ s% xtra1, s% xtra5, s% xtra9, s% xtra13, s% xtra17 /)
         lower_ms=(/ s% xtra4, s% xtra8, s% xtra12, s% xtra16, s% xtra20 /)
         omsq_array=(/ s% xtra3, s% xtra7, s% xtra11, s% xtra15, s% xtra19 /)
         lmax_array=(/ s% xtra2, s% xtra6, s% xtra10, s% xtra14, s% xtra18 /)

         length=count(k_max_array/=0)
  !        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Calculate maximum nonlinearity!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 	  	allocate(burntypes_array(5),maxburn_array(5),krxir_max_array(5),max_nl_rad_array(5),Lheatl_array(5))
 	  	allocate(ell_maxnl_array(5),dLheatl_max_array(5),ell_maxL_array(5))
	  	Lheatl_array(:) = 0.
	  	dLheatl_max_array(:) = 0.
	  	ell_maxL_array(:) = 0.
	  	upper_ms(:)=0.
	 	krxir_max_array(:) = 0.
	  	max_nl_rad_array(:) = 0.
	  	ell_maxnl_array(:)=0.
	  	burntypes_array(:) =  0. !possible burn labels
	  	maxburn_array = (/ 1, 4, 12, 16, 20, 28 /) !burn labels per shell by atomic number (/h, he, c, o, ne, si /)
	!   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Calculate transmission coefficient!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		dLheatl_array(:,:) = 0.
		fesc_array(:,:) = 0.
		do i = 1, length
			k_max=k_max_array(i)
			omsq=omsq_array(i)
			lmax=lmax_array(i)
			L_waves=L_waves_max_array(i)*Lsun
			krxir_max=0 !the max nonlinearity
			max_nl_rad=0 !radius of the maximum nonlinearity
			krxir_maxofell=0
			ell_maxnl = 0
			max_nl_radofell = 0
			intmin = 0.
			
			Lheatl=0.
			dLheatl_max=0.
			ell_maxL=0.

		    do l=1, 10 !loop over l=1 to l=10
				intkrl1min = 0.
				
				do k=1,k_max !this is for all zones outside wave excitation point 

					minflag = 0.

					! vel = 0.
					! if (abs(s% v(k)) < 1d11) vel = abs(s% v(k)) !Required if v turned off
					!print tau(k), vel
					!if ( (s% tau(k) > 1d4) .and. (vel < s% csound(k)) ) then !Avoid outflowing surface regions
					    !is this generally true for entire mass range

					if (s% brunt_N2(k) < 0.) brunt = 0. !Avoid artificial drops in Brunt
					if (s% brunt_N2(k) > 0.) brunt = s% brunt_N2(k) !Avoid artificial drops in Brunt

					lambl1 = l*(l+1)*pow2(s% csound(k)/(s% r(k)))

					krl1 = (brunt - omsq)*(lambl1 - omsq)/(omsq*pow2(s% csound(k)))

					dr = s% dm(k)/(4.*pi*s% rho(k)*pow2(s% r(k)))

					!Find smallest transmission coefficient and find f_damp_int
					if (krl1 < 0.) then
						minflag = 1.
						intmin = intmin + dr*sqrt(-krl1)
					end if
					if (minflag < 0.5) then
						if (intmin > intkrl1min) then 
							intkrl1min = intmin
						end if
						intmin = 0.
					end if
				!end if
				end do !loop over k

				Tl1min = exp_cr(-2.*intkrl1min) !it says l1 but it's for whatever l value in loop

				tprop = 0.
				f_damp_int = 0.
				x_fac = 0.
				f_fac = 0.
				do k=k_max,1,-1
					lambl1 = l*(l+1)*pow2(s% csound(k)/(s% r(k)))
					if (s% brunt_N2(k) < 0.) brunt = 0. 
					if (s% brunt_N2(k) > 0.) brunt = s% brunt_N2(k) 
					krl1 = (brunt - omsq)*(lambl1 - omsq)/(omsq*pow2(s% csound(k)))
					dr = s% dm(k)/(4.*pi*s% rho(k)*pow2(s% r(k)))
					if (krl1 > 0.) then
						!Damp energy via radiative diffusion
						diff_t = (16d0*boltz_sigma*pow3(s% T(k)))/(3d0*pow2(s% rho(k))*s% cp(k)*s% opacity(k))
						gamphot = krl1*diff_t
						if (omsq > brunt) then
							vg = s% csound(k)
							tprop = tprop + dr/vg
						endif
						if (omsq < brunt) then
							vg = omsq*s% r(k)/sqrt(l*(l+1)*brunt) !l(l+1) = 2 for l=1
							tprop = tprop + dr/vg
							gamnu = -pow2(s% grav(k)*s% gamma1(k)*s% grada(k))
							gamnu = gamnu*(s% d_epsnuc_dlnT(k) - s% d_nonnucneu_dlnT(k))/&
									(brunt*pow2(s% csound(k)*s% csound(k)))
							if (gamnu>0.) then
								f_damp_int = f_damp_int + 2.*gamnu*dr/vg
							end if
							if (gamphot>0.) then
								f_damp_int = f_damp_int + 2.*gamphot*dr/vg
							end if
						end if
					end if
				end do
				f_damp = exp_cr(f_damp_int)
				!f_fac = 1.-1./f_damp
				f_fac = f_damp - 1.
				x_fac = f_damp_int
				heatfrac = (1.)/(1. + f_fac/Tl1min)

				Lfracpt1a = -2.*exp_cr(-pow2(DFLOAT(l))/pow2(lmax))* &
									(pow2(DFLOAT(l))+DFLOAT(l)*lmax+pow2(lmax))/pow2(lmax)
				Lfracpt2a = sqrt(pi)*Erf(DFLOAT(l)/lmax)
				Lfrac_a = (Lfracpt1a+Lfracpt2a)/4.

				Lfracpt1b = -2.*exp_cr(-pow2(DFLOAT(l)-1)/pow2(lmax))* &
									(pow2(DFLOAT(l)-1)+(DFLOAT(l)-1)*lmax+pow2(lmax))/pow2(lmax)
				Lfracpt2b = sqrt(pi)*Erf((DFLOAT(l)-1)/lmax)
				Lfrac_b = (Lfracpt1b+Lfracpt2b)/4.

				Lfrac_int = Lfrac_a - Lfrac_b
				Ldenom= (2.+sqrt(pi))/4.

				Lfrac = Lfrac_int/Ldenom
				! heatfrac=1./(1.+1./Tl1min)
				dLheatl=L_waves*Lfrac*heatfrac
				Lheatl = Lheatl + dLheatl  ! add up the heating rate over l=1 to 10 for this shell

				if (dLheatl > dLheatl_max)  then
					dLheatl_max = dLheatl
					ell_maxL = FLOAT(l)
				end if
				!for saving in columns
				dLheatl_array(i,l) = dLheatl
				fesc_array(i,l) = heatfrac 
					!calculate the nonlinearity for each l value
				do k=1,k_max !this is for all zones outside wave excitation point

					lambl1 = l*(l+1)*pow2(s% csound(k)/(s% r(k)))

					if (s% brunt_N2(k) < 0.) brunt = 0. !Avoid artificial drops in Brunt
					if (s% brunt_N2(k) > 0.) brunt = s% brunt_N2(k) !Avoid artificial drops in Brunt

					krl1 = (brunt - omsq)*(lambl1 - omsq)/(omsq*pow2(s% csound(k))) !square of radial wavenumber

					dr = s% dm(k)/(4.*pi*s% rho(k)*pow2(s% r(k)))

					if ((omsq < brunt) .and. (omsq < lambl1)) then !finding the g-mode cavity
						!calculate the nonlinearity
						xirl1=(2/heatfrac)*(dLheatl*sqrt(DFLOAT(l)*(DFLOAT(l)+1)))/&
								(4.*pi*s% rho(k)*pow_cr(s% r(k),3.d0)*sqrt(brunt)*omsq)
						!xirl1 is radial displacement squared
						krxir=sqrt(krl1*xirl1)
						if ((krxir > krxir_max) .and. sqrt(brunt/lambl1)<5d0)  then
							krxir_max = krxir 
							max_nl_rad = s% r(k)
						end if
					end if
					
				end do !loop over k to find nonlinearity

				if (krxir_max > krxir_maxofell)  then
					krxir_maxofell = krxir_max
					max_nl_radofell = max_nl_rad
					ell_maxnl = FLOAT(l)
				end if
		    end do !loop over l

		    if (s% center_he4 < 1d-3) then 
		    	Lheatl_array(i) = Lheatl
		    	dLheatl_max_array(i) = dLheatl_max
		    	ell_maxL_array(i) = ell_maxL
			    krxir_max_array(i) = krxir_maxofell
			    max_nl_rad_array(i) = max_nl_radofell
			    ell_maxnl_array(i) = ell_maxnl
			endif
		end do !loop over shells

		do i = 1, length
			k_max=k_max_array(i)
			upper_ms(i) = s% m(k_max)/Msun

			burn_h = 0 !power in burning
			burn_he = 0
			burn_c = 0
			burn_o = 0
			burn_ne = 0
			burn_si = 0

			do k=1,s% nz !first sum up all the eps_nuc (erg/s/g) within the shell 
				if ( (s% m(k)/Msun > lower_ms(i)) .and. (s% m(k)/Msun < upper_ms(i)) ) then 
					burn_h = burn_h + (s% eps_nuc_categories(ipp,k) + s% eps_nuc_categories(icno,k))
					burn_he = burn_he + s% eps_nuc_categories(i3alf,k)
					burn_c = burn_c + (s% eps_nuc_categories(icc,k) + s% eps_nuc_categories(ico,k))
					burn_o = burn_o + (s% eps_nuc_categories(ico,k) + s% eps_nuc_categories(ioo,k))
					burn_ne = burn_ne + s% eps_nuc_categories(i_burn_ne,k)
					burn_si = burn_si + s% eps_nuc_categories(i_burn_si,k)

				end if		
			end do
			
			maxburn = MAXLOC( (/burn_h, burn_he, burn_c, burn_o, burn_ne, burn_si/), 1)
			!choose label based on the maximum eps_nuc
			burntypes_array(i) = maxburn_array(maxburn)

		end do

		!add L_waves_max; max_nl_rad; add dLheat and escape fraction in each \ell value

		lchar = '  '
		ichar = ' '
		fescname = 'fescl00si'
		dLheatname = 'dLheatl00si'
		k = 0
		do i = 1, 5
			do l = 1, 10
				lchar = '  '
				ichar = ' '
				k = 10*(i-1)+l

				write(lchar,'(I0.2)') l
				write(ichar,'(I1)') i

				dLheatname(8:9) = lchar
				fescname(6:7) = lchar
				dLheatname(11:11) = ichar
				fescname(9:9) = ichar

				names(k) = dLheatname
				vals(k) = dLheatl_array(i,l)
				names(k+50) = fescname
				vals(k+50) = fesc_array(i,l)

			end do
		end do


		names(101) = 'Lheatl1'
		vals(101) = Lheatl_array(1)/Lsun
		names(102)= 'Lheatl2'
		vals(102) = Lheatl_array(2)/Lsun
		names(103)= 'Lheatl3'
		vals(103) = Lheatl_array(3)/Lsun
		names(104)= 'Lheatl4'
		vals(104) = Lheatl_array(4)/Lsun
		names(105)= 'Lheatl5'
		vals(105) = Lheatl_array(5)/Lsun

		names(106) = 'burntype1'
		vals(106) = burntypes_array(1)
		names(107) = 'burntype2'
		vals(107) = burntypes_array(2)
		names(108) = 'burntype3'
		vals(108) = burntypes_array(3)
		names(109) = 'burntype4'
		vals(109) = burntypes_array(4)
		names(110) = 'burntype5'
		vals(110) = burntypes_array(5)

		names(111) = 'omsq1'
		vals(111) = omsq_array(1)
		names(112) = 'omsq2'
		vals(112) = omsq_array(2)
		names(113) = 'omsq3'
		vals(113) = omsq_array(3)
		names(114) = 'omsq4'
		vals(114) = omsq_array(4)
		names(115) = 'omsq5'
		vals(115) = omsq_array(5)

		names(116) = 'lmax1'
		vals(116) = lmax_array(1)
		names(117) = 'lmax2'
		vals(117) = lmax_array(2)
		names(118) = 'lmax3'
		vals(118) = lmax_array(3)
		names(119) = 'lmax4'
		vals(119) = lmax_array(4)
		names(120) = 'lmax5'
		vals(120) = lmax_array(5)

		names(121) = 'top1'
		vals(121) = upper_ms(1)
		names(122) = 'top2'
		vals(122) = upper_ms(2)
		names(123) = 'top3'
		vals(123) = upper_ms(3)
		names(124) = 'top4'
		vals(124) = upper_ms(4)
		names(125) = 'top5'
		vals(125) = upper_ms(5)

		names(126) = 'bottom1'
		vals(126) = lower_ms(1)
		names(127) = 'bottom2'
		vals(127) = lower_ms(2)
		names(128) = 'bottom3'
		vals(128) = lower_ms(3)
		names(129) = 'bottom4'
		vals(129) = lower_ms(4)
		names(130) = 'bottom5'
		vals(130) = lower_ms(5)

         names(131)='E_waves1'
         names(132)='E_waves2'
         names(133)='E_waves3'
         names(134)='E_waves4'
         names(135)='E_waves5'
 		 vals(131) = accum_ener1 
		 vals(132) = accum_ener2 
		 vals(133) = accum_ener3 
		 vals(134) = accum_ener4 
		 vals(135) = accum_ener5 

 		names(136) = 'krxir_max1'
		vals(136) = krxir_max_array(1)
		names(137) = 'krxir_max2'
		vals(137) = krxir_max_array(2)
		names(138) = 'krxir_max3'
		vals(138) = krxir_max_array(3)
		names(139) = 'krxir_max4'
		vals(139) = krxir_max_array(4)
		names(140) = 'krxir_max5'
		vals(140) = krxir_max_array(5)

 		names(141) = 'ell_maxnl1'
		vals(141) = ell_maxnl_array(1)
		names(142) = 'ell_maxnl2'
		vals(142) = ell_maxnl_array(2)
		names(143) = 'ell_maxnl3'
		vals(143) = ell_maxnl_array(3)
		names(144) = 'ell_maxnl4'
		vals(144) = ell_maxnl_array(4)
		names(145) = 'ell_maxnl5'
		vals(145) = ell_maxnl_array(5)

 		names(146) = 'max_nl_rad1'
		vals(146) = max_nl_rad_array(1)
		names(147) = 'max_nl_rad2'
		vals(147) = max_nl_rad_array(2)
		names(148) = 'max_nl_rad3'
		vals(148) = max_nl_rad_array(3)
		names(149) = 'max_nl_rad4'
		vals(149) = max_nl_rad_array(4)
		names(150) = 'max_nl_rad5'
		vals(150) = max_nl_rad_array(5)

 		names(151) = 'dLheatl_max1'
		vals(151) = dLheatl_max_array(1)
		names(152) = 'dLheatl_max2'
		vals(152) = dLheatl_max_array(2)
		names(153) = 'dLheatl_max3'
		vals(153) = dLheatl_max_array(3)
		names(154) = 'dLheatl_max4'
		vals(154) = dLheatl_max_array(4)
		names(155) = 'dLheatl_max5'
		vals(155) = dLheatl_max_array(5)

 		names(156) = 'ell_maxL1'
		vals(156) = ell_maxL_array(1)
		names(157) = 'ell_maxL2'
		vals(157) = ell_maxL_array(2)
		names(158) = 'ell_maxL3'
		vals(158) = ell_maxL_array(3)
		names(159) = 'ell_maxL4'
		vals(159) = ell_maxL_array(4)
		names(160) = 'ell_maxL5'
		vals(160) = ell_maxL_array(5)

 		names(161) = 'L_waves_max1'
		vals(161) = L_waves_max_array(1)
		names(162) = 'L_waves_max2'
		vals(162) = L_waves_max_array(2)
		names(163) = 'L_waves_max3'
		vals(163) = L_waves_max_array(3)
		names(164) = 'L_waves_max4'
		vals(164) = L_waves_max_array(4)
		names(165) = 'L_waves_max5'
		vals(165) = L_waves_max_array(5)

         deallocate(k_max_array,lower_ms,upper_ms,omsq_array,L_waves_max_array,lmax_array,ell_maxnl_array)
         deallocate(burntypes_array,maxburn_array,krxir_max_array,Lheatl_array,dLheatl_max_array)
		 deallocate(ell_maxL_array,max_nl_rad_array)
         !note: do NOT add the extras names to history_columns.list
         ! the history_columns.list is only for the built-in log column options.
         ! it must not include the new column names you are adding here.

      end subroutine data_for_extra_history_columns



      
      ! None of the following functions are called unless you set their
      ! function point in extras_control.
      
      
      integer function extras_startup(id, restart, ierr)
         integer, intent(in) :: id
         logical, intent(in) :: restart
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_startup = 0
         if (.not. restart) then
            call alloc_extra_info(s)
            accum_ener1=0 !for energy accumulation
	        accum_ener2=0
	        accum_ener3=0
	        accum_ener4=0
	        accum_ener5=0
	        E_heat_accum=0
         else ! it is a restart
            call unpack_extra_info(s)
         end if

         	delta_ener1=0
	        delta_ener2=0
	        delta_ener3=0
	        delta_ener4=0
	        delta_ener5=0
	        E_heat_delta=0
      end function extras_startup
      

      ! returns either keep_going, retry, backup, or terminate.
      integer function extras_check_model(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_check_model = keep_going         
         if (.false. .and. s% star_mass_h1 < 0.35d0) then
            ! stop when star hydrogen mass drops to specified level
            extras_check_model = terminate
            write(*, *) 'have reached desired hydrogen mass'
            return
         end if


         ! if you want to check multiple conditions, it can be useful
         ! to set a different termination code depending on which
         ! condition was triggered.  MESA provides 9 customizeable
         ! termination codes, named t_xtra1 .. t_xtra9.  You can
         ! customize the messages that will be printed upon exit by
         ! setting the corresponding termination_code_str value.
         ! termination_code_str(t_xtra1) = 'my termination condition'

         ! by default, indicate where (in the code) MESA terminated
         if (extras_check_model == terminate) s% termination_code = t_extras_check_model
      end function extras_check_model


      

      ! returns either keep_going or terminate.
      ! note: cannot request retry or backup; extras_check_model can do that.
      integer function extras_finish_step(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_finish_step = keep_going
         call store_extra_info(s)

         accum_ener1 = accum_ener1 + delta_ener1
         accum_ener2 = accum_ener2 + delta_ener2
         accum_ener3 = accum_ener3 + delta_ener3
         accum_ener4 = accum_ener4 + delta_ener4
         accum_ener5 = accum_ener5 + delta_ener5
         E_heat_accum = E_heat_accum + E_heat_delta


         	delta_ener1=0
	        delta_ener2=0
	        delta_ener3=0
	        delta_ener4=0
	        delta_ener5=0
	        E_heat_delta=0
         ! to save a profile, 
            ! s% need_to_save_profiles_now = .true.
         ! to update the star log,
            ! s% need_to_update_history_now = .true.

         ! see extras_check_model for information about custom termination codes
         ! by default, indicate where (in the code) MESA terminated
         if (extras_finish_step == terminate) s% termination_code = t_extras_finish_step
      end function extras_finish_step
      
      
      subroutine extras_after_evolve(id, id_extra, ierr)
         integer, intent(in) :: id, id_extra
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
      end subroutine extras_after_evolve
      
      
      ! routines for saving and restoring extra data so can do restarts
         
         ! put these defs at the top and delete from the following routines
         !integer, parameter :: extra_info_alloc = 1
         !integer, parameter :: extra_info_get = 2
         !integer, parameter :: extra_info_put = 3
      
      
      subroutine alloc_extra_info(s)
         integer, parameter :: extra_info_alloc = 1
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_alloc)
      end subroutine alloc_extra_info
      
      
      subroutine unpack_extra_info(s)
         integer, parameter :: extra_info_get = 2
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_get)
      end subroutine unpack_extra_info
      
      
      subroutine store_extra_info(s)
         integer, parameter :: extra_info_put = 3
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_put)
      end subroutine store_extra_info
      
      
      subroutine move_extra_info(s,op)
         integer, parameter :: extra_info_alloc = 1
         integer, parameter :: extra_info_get = 2
         integer, parameter :: extra_info_put = 3
         type (star_info), pointer :: s
         integer, intent(in) :: op
         
         integer :: i, j, num_ints, num_dbls, ierr
         
         i = 0
         ! call move_int or move_flg    
         num_ints = i
         
         i = 0
         ! call move_dbl       
         
         call move_dbl(accum_ener1)
         call move_dbl(accum_ener2)
         call move_dbl(accum_ener3)
         call move_dbl(accum_ener4)
         call move_dbl(accum_ener5)
         call move_dbl(E_heat_accum)

         call move_dbl(delta_ener1)
         call move_dbl(delta_ener2)
         call move_dbl(delta_ener3)
         call move_dbl(delta_ener4)
         call move_dbl(delta_ener5)
         ! call move_dbl(E_heat_delta)

         num_dbls = i
         
         if (op /= extra_info_alloc) return
         if (num_ints == 0 .and. num_dbls == 0) return
         
         ierr = 0
         call star_alloc_extras(s% id, num_ints, num_dbls, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in star_alloc_extras'
            write(*,*) 'alloc_extras num_ints', num_ints
            write(*,*) 'alloc_extras num_dbls', num_dbls
            stop 1
         end if
         
         contains
         
         subroutine move_dbl(dbl)
            real(dp) :: dbl
            i = i+1
            select case (op)
            case (extra_info_get)
               dbl = s% extra_work(i)
               !dbl = s% xtra1_array(i)
            case (extra_info_put)
               s% extra_work(i) = dbl
	       !s% xtra1_array(i) = dbl
            end select
         end subroutine move_dbl
         
         subroutine move_int(int)
            integer :: int
            i = i+1
            select case (op)
            case (extra_info_get)
               int = s% extra_iwork(i)
            case (extra_info_put)
               s% extra_iwork(i) = int
            end select
         end subroutine move_int
         
         subroutine move_flg(flg)
            logical :: flg
            i = i+1
            select case (op)
            case (extra_info_get)
               flg = (s% extra_iwork(i) /= 0)
            case (extra_info_put)
               if (flg) then
                  s% extra_iwork(i) = 1
               else
                  s% extra_iwork(i) = 0
               end if
            end select
         end subroutine move_flg
      
      end subroutine move_extra_info

      
      










      end module run_star_extras
      
