! ::::::::::::::::
!
! run_star_extras including an experimental wind routine written by Z. Keszthelyi and J. Puls. 
!
! ::::::::::::::::

            module run_star_extras

      use star_lib
      use star_def
      use const_def
      use star_private_def
      use chem_def, only: ih1, ihe4, ic12, ic13, in14, io16
      
      implicit none 

      ! these routines are called by the standard run_star check_model

      contains
      
      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

         ! this is the place to set any procedure pointers you want to change
         ! e.g., other_wind, other_mixing, other_energy  (see star_data.inc)    

         s% other_wind => experimental  ! experimental wind routine used by Keszthelyi, Puls, Wade (2017) A&A 598, A4 
       
         ! Uncomment these lines if you wish to use the functions in this file
         ! otherwise we use a null_ version which does nothing.
        
        ! 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

         ! Once you have set the function pointers you want,
         ! then uncomment this (or set it in your star_job inlist)
         ! to disable the printed warning message,
          s% job% warn_run_star_extras =.false.

      end subroutine extras_controls

     

!.............................................................................

! The subroutine below is the experimental wind routine used by
! Keszthelyi, Puls, Wade (2017) A&A, 598, A4 (KPW2017, hereafter).
! The mass loss calculation is based on the Wind-momentum Luminosity
! Relation [WLR; Kudritzki, Lennon, Puls (1995) svlt.conf, 246;
! Puls et al. (1996) A&A, 305, 171; Kudritzki & Puls (2000), ARA&A, 38, 613],
! and allows for testing different scenarios around the debated bi-stability jumps
! related to Fe recombination in the wind [Pauldrach & Puls (1990) A&A, 237, 409;
! Vink, de Koter, Lamers (1999) A&A, 350, 181; Vink, de Koter, Lamers (2000)
! A&A, 362, 295; Vink, de Koter, Lamers (2001) A&A, 369, 574].
! This routine is applicable for hot (10 - 50 kK) massive stars.
!
! Since this routine requires some experience/understanding of stellar
! winds and related problems, this routine is *NOT* a default
! wind scheme that could be adopted in 'standard runs'.  

!.............................................................................

subroutine experimental(id, L_phot,M_phot, R_phot, T_phot, w, ierr)   

  type (star_info), pointer :: s
  real(dp), intent(in) :: L_phot, M_phot, R_phot, T_phot ! photospheric values (cgs)
  integer, intent(in) :: id
  integer, intent(out) :: ierr
  integer :: h1, he4, nz
  real(dp) :: wind, alfa1, alfa2, I_he, &   ! alfa1 and alfa2 are used only for interpolation
  X, Y, Z, L1, M1, R1, T1, &
  YHe, &                         ! Helium abundance by number, N_He/N_H  
  w1, w2, w3, logvesc, vesc, &   ! w1, w2, w3 are mass-loss rates being partitioned
  Gamma_e, sigma_e, logvinf      ! Eddington factor for pure electron scattering, electron scattering
                                 ! cross section, and the log of the terminal wind velocity   
  real(dp) :: surface_h1, surface_he4     
  logical, parameter :: dbg = .false.
  real(dp), intent(out) :: w
  real(dp) ::  logMdot, logDmom, logD0_2, logD0_3     ! log mass-loss rates, log wind momentum,
                                                      ! and constants to increase the wind momentum over the bi-stability 

! MESA constants are defined in const_def.f in mesa/const/public, pi4 is 4 * pi,
! gravitational constant is 'standard_cgrav', lightspeed is 'clight', seconds in a year is 'secyer'.

!=== INPUT PARAMETER VALUES ===

! --- Values related to the WLR ---

  real(dp), parameter :: alphaprime = 0.543d0
! alphaprime = alpha - delta (force-multiplier parameters), corresponding to the inverse slope
! of the logarithmic WLR, Eq. 13 by KPW2017.
! See Pauldrach, Puls, & Kudritzki, (1986) A&A, 164, 86;
! and Puls, Springmann, Lennon (2000) A&A 141, 23. 
! Literature values are listed in Table 3 by KPW2017.  

  real(dp), parameter :: logD0_1 = 18.40d0  
! The offset of the logarithmic WLR before (hot side) the first bi-stability,
! Eq. 12 and Table 3 by KPW2017.
! The value of this offset is changed after the bi-stability regions. Units refer to cgs.
! Clumping-corrected values can be different, 
! e.g. 20.23d0 by Mokiem et al. (2007) A&A 473, 603.

  real(dp), parameter :: deltaD0_1 = 0.35d0
! The size (=increase of log Dmom) over the first bi-stability jump, Eq. 21 by KPW2017. 

  real(dp), parameter :: deltaD0_2 = 0.d0
! The size of the second bi-stability jump, assumed to be zero here

  real(dp), parameter :: fscal = 1.d0
! Global scaling factor 


! --- Values related to the bi-stability ---

! The calculations are partitioned into three Teff regimes. 

  real(dp), parameter :: Teff_jump1 = 20.d3  
! FIRST bi-stability jump temperature (from observations)
  
  real(dp), parameter :: Teff_jump2 = 10.d3  
! SECOND bi-stability jump temperature (from observations) 

  real(dp), parameter :: dT1 = 3.d3
! Interpolation width of the first bi-stability jump. 
! Note that a small dT results in a steep change
! in Mdot while a larger dT yields a more gradual change. 

  real(dp), parameter :: dT2 = 1.5d3
! Interpolation width of the second bi-stability jump. 


! --- Values related to terminal wind velocity

  real(dp), parameter :: fvinf1 = 2.6d0  ! ratio of terminal to escape velocity before Teff_jump_1
  real(dp), parameter :: fvinf2 = 1.3d0  ! ratio of terminal to escape velocity after Teff_jump_1 and before Teff_jump_2
  real(dp), parameter :: fvinf3 = 0.7d0  ! ratio of terminal to escape velocity after Teff_jump_2

! Values from Kudritzki & Puls (2000), ARA&A, 38, 613; and Markova & Puls (2008), A&A, 478, 823

! Note: As of yet, a steep change is adopted, however observations point to a gradual change
! in v_inf / v_esc (KPW2017). 


! === calculate WLR offset for after the bi-stabilities

  logD0_2 = logD0_1 + deltaD0_1 
! The offset of the WLR after the first bi-stability (and before the second bi-stability). 

  logD0_3 = logD0_2 + deltaD0_2 
! The offset of the WLR after the second bi-stability. 


!=== end of input values ===

  call get_star_ptr(id,s,ierr)            

!--- elements ---- 

  nz = s% nz
  h1 = s% net_iso(ih1)
  he4 = s% net_iso(ihe4)
  surface_h1 = s% xa(h1,1)
  surface_he4 = s% xa(he4,1)
  X = surface_h1            ! mass fraction of H
  Y = surface_he4           ! mass fraction of He
  Z = 1.d0 - (X + Y)        ! surface metallicity 

  YHe = Y/X / 4.d0  
! surface Helium content by number, 
! n_He/n_H; approximation: factor 4 should be m_He/m_H


!=== photospheric parameters in solar units ===

  L1 = L_phot
  M1 = M_phot
  T1 = T_phot
  R1 = s% photosphere_r 

!=== Stop if parameter input needs adjustment ===

 if (Teff_jump1 < Teff_jump2) stop 'adjust bi-stability jump temperatures in subr. experimental '

 if (dT1 < 0.d0 .or. dT2 < 0.d0) stop 'adjust the interpolation width of the bi-stability jump in subr. experimental'

 if (Teff_jump1 - dT1 .le. Teff_jump2 + dT2) stop 'adjust jump temperature and interpolation width in subr. experimental'

 if (T1 < Teff_jump2 - dT2) stop 'This Teff is below the applicability of this wind scheme' 


!-------------------------------------

! === Calculation of Mdot for Teff > Teff_jump1 + dT  and for jump region (hot) ===

! --- PATH 1 --- 

! === Calculation of Mdot for T > Teff_jump1 - dT1 ===
!     corresponding to hot region and first jump (hot solution)

 if (T1 .ge. Teff_jump1 - dT1) then 

    I_he = 2.d0  ! the number of free electrons per He nuclei (hot regime)
    sigma_e = 0.398d0 * (1.d0 + YHe * I_he) / (1.d0 + 4.d0 * YHe)  ! Eq. 66 by Kudritzki et al. (1989) A&A, 219, 205 
    Gamma_e =  sigma_e * L1/(M1 * pi4 * standard_cgrav * clight) 
   
! escape velocity corrected for electron scattering 

    vesc = sqrt((1.d0 - Gamma_e) * 2.d0 * standard_cgrav * M1 / R1 )

! Eq. 15 KPW2017
    
    logvinf = log10_cr(fvinf1) + log10_cr(vesc)  ! log of the terminal velocity for this Teff range

! === The WLR, Eq. 12 KPW2017 ===  

    logDmom = log10_cr(L1/Lsun)/alphaprime + logD0_1  
    
! Eq. 19 KPW2017

    logMdot = logDmom - logvinf - 0.5d0 * log10_cr(R1/Rsun) + log10_cr(secyer) - log10_cr(Msun)

! Since the first three terms on the r.h.s. are in cgs, but the output mass-loss rate needs to
! be in solar units, it is scaled by the last two terms to get units of Msun/yr.
    
    w1 =  exp10_cr(logMdot)
    w = w1 
 endif

! --- PATH 2 --- 

! === Calculation of Mdot for Teff_jump1 + dT1  > T > Teff_jump2 - dT1  ===
!     corrsponding to first jump (cool solution), cool region and 2nd jump (warm solution)  

 if (T1 .le. Teff_jump1 + dT1 .and. T1 .ge. Teff_jump2 - dT2) then 

    I_he = 1.d0     ! cooler side of jump, Helium recombined to HeII 
    sigma_e = 0.398d0 * (1.d0 + YHe * I_he) / (1.d0 + 4.d0 * YHe)  
    Gamma_e =  sigma_e * L1/(M1 * pi4 * standard_cgrav * clight) 
    vesc = sqrt(2.d0 * standard_cgrav * M1 / R1 * (1.d0-Gamma_e))
    logvinf = log10_cr(fvinf2) + log10_cr(vesc) ! log of the terminal velocity for this Teff range
    logDmom = log10_cr(L1/Lsun)/alphaprime + logD0_2 
    logMdot = logDmom - logvinf - 0.5d0 *log10_cr(R1/Rsun) + log10_cr(secyer) - log10_cr(Msun)  
    w2 = exp10_cr(logMdot)
    w = w2
 endif

! --- INTERPOLATION for the FIRST bi-stability jump ---  

 if (Teff_jump1 + dT1 > T1 .and. T1 > Teff_jump1 - dT1) then
    alfa1 = (T1 - (Teff_jump1 - dT1)) / (2.d0 * dT1)
    w = alfa1 * w1 + (1.d0 - alfa1) * w2
 endif

! --- PATH 3 ---

! === Calculation of Mdot for Teff < Teff_jump2 + dT2 ===
!     corresponding to 2nd jump (cool solution) and coolest region

 if (T1 .le. Teff_jump2 + dT2) then
     
    sigma_e = 0.398d0 * 1.d0 / (1.d0 + 4.d0 * YHe) ! since I_he = 0 for T1 < 10 kK. 
    Gamma_e =  sigma_e * L1/(M1 * pi4 * standard_cgrav * clight) 
    vesc = sqrt(2.d0 * standard_cgrav * M1 / R1 * (1.d0-Gamma_e))   
    logvinf = log10_cr(fvinf3) + log10_cr(vesc) ! log of the terminal velocity for this Teff range
    logDmom = log10_cr(L1/Lsun)/alphaprime + logD0_3 
    logMdot = logDmom - logvinf - 0.5d0 *log10_cr(R1/Rsun) + log10_cr(secyer) - log10_cr(Msun)  
    w3 = exp10_cr(logMdot)
    w = w3
 endif

! --- INTERPOLATION for the SECOND bi-stability jump ---    

 if (Teff_jump2 + dT2 > T1 .and. T1 > Teff_jump2 - dT2) then
    alfa2 = (T1 - (Teff_jump2 - dT2)) / (2.d0 * dT2)  
    w = alfa2 * w2 + (1.d0 - alfa2) * w3
 endif

!:::::::: FINAL MDOT :::::::::::::::::::::: 

  w = w * fscal  ! fscal is a global scaling factor 

if (dbg) write(*,*) 'wind', w
    
end subroutine experimental


! This is now the rest from 'standard_run_star_extras.inc’ 
  
      integer function extras_startup(s, id, restart, ierr)
         type (star_info), pointer :: s
         integer, intent(in) :: id
         logical, intent(in) :: restart
         integer, intent(out) :: ierr
         ierr = 0
         extras_startup = 0
         if (.not. restart) then
            call alloc_extra_info(s)
         else ! it is a restart
            call unpack_extra_info(s)
         end if
      end function extras_startup
      

      ! returns either keep_going, retry, backup, or terminate.
      integer function extras_check_model(s, id, id_extra)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         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


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 = 0      ! don't forget to synchronize this.
end function how_many_extra_history_columns
      
      
      
  subroutine data_for_extra_history_columns(id, id_extra, n, names, vals, ierr)

   ! MESA provides routines for taking logarithms which will handle
   ! cases where you pass them zero, etc.  Take advantage of this!
   ! use crlibm_lib, only: safe_log10_cr

    use const_def, only: dp
    use star_def, only: maxlen_history_column_name

   integer, intent(in) :: id, id_extra, n
   character (len=maxlen_history_column_name) :: names(n)
   real(dp) :: vals(n)
   integer, intent(out) :: ierr
   type (star_info), pointer :: s
    
   ierr = 0
   call star_ptr(id, s, ierr)
   if (ierr /= 0) return


  ! extra column 1
  !   names(1) = ""
  !   vals(1) = 

  ! extra column 2
  !   names(2) = ""
  !   vals(2) = 

      end subroutine data_for_extra_history_columns

      
      integer function how_many_extra_profile_columns(s, id, id_extra)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         how_many_extra_profile_columns = 0
      end function how_many_extra_profile_columns
      
      
      subroutine data_for_extra_profile_columns(s, id, id_extra, n, nz, names, vals, ierr)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra, n, nz
         character (len=maxlen_profile_column_name) :: names(n)
         real(dp) :: vals(nz,n)
         integer, intent(out) :: ierr
         integer :: k
         ierr = 0
         
         !note: do NOT add these 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.

         ! here is an example for adding a profile column
         !if (n /= 1) stop 'data_for_extra_profile_columns'
         !names(1) = 'beta'
         !do k = 1, nz
         !   vals(k,1) = s% Pgas(k)/s% P(k)
         !end do
         
      end subroutine data_for_extra_profile_columns
      

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

         ! 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(s, id, id_extra, ierr)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         integer, intent(out) :: ierr
         ierr = 0
      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       
         
         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)
            case (extra_info_put)
               s% extra_work(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