! *********************************************************************** ! ! 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 ! ! *********************************************************************** module run_star_extras use star_lib use star_def use const_def use crlibm_lib use chem_def use rates_def, only: i_rate use interp_1d_def, only: pm_work_size use interp_1d_lib, only: interp_pm, interp_values, interp_value implicit none integer :: time0, time1, clock_rate, redo_count logical :: dbg = .false. logical :: pgstar_flag 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 s% extras_startup => extras_startup s% extras_start_step => extras_start_step 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 pgstar_flag = s% job% pgstar_flag ! this is optional if(s%use_other_wind) then if(s%x_integer_ctrl(2)==1) then s% other_wind => brott_wind ! Hamen else s% other_wind=>Dutch_tanh_intrp_wind ! Tramper and N&L end if end if s% other_adjust_mdot => my_adjust_mdot s% other_before_struct_burn_mix => my_before_struct_burn_mix if(s%use_other_neu) then s% other_neu => my_neu end if end subroutine extras_controls subroutine brott_wind(id, Lsurf, Msurf, Rsurf, Tsurf, w, ierr) use star_def integer, intent(in) :: id real(dp), intent(in) :: Lsurf, Msurf, Rsurf, Tsurf ! surface values (cgs) ! NOTE: surface is outermost cell. not necessarily at photosphere. ! NOTE: don't assume that vars are set at this point. ! so if you want values other than those given as args, ! you should use values from s% xh(:,:) and s% xa(:,:) only. ! rather than things like s% Teff or s% lnT(:) which have not been set yet. real(dp), intent(out) :: w ! wind in units of Msun/year (value is >= 0) integer, intent(out) :: ierr integer :: h1, he4 real(dp) :: Xs, Ys, Z_div_Z_solar, Teff_jump, alfa, L1, M1, R1, T1, & vink_wind, nieu_wind, hamann_wind type (star_info), pointer :: s ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return L1 = Lsurf M1 = Msurf R1 = Rsurf T1 = Tsurf h1 = s% net_iso(ih1) he4 = s% net_iso(ihe4) Xs = s% xa(h1,1) Ys = s% xa(he4,1) ! Z=0.0142 is Z from Asplund et al. 2009 if(s%x_integer_ctrl(4)==1) then Z_div_Z_solar = 1d-3/0.0142d0 else Z_div_Z_solar = s% Zbase/0.0142d0 end if ! use Vink et al 2001, eqns 14 and 15 to set "jump" temperature Teff_jump = 1d3*(61.2d0 + 2.59d0*(-13.636d0 + 0.889d0*log10_cr(Z_div_Z_solar))) vink_wind = 0d0 nieu_wind = 0d0 hamann_wind = 0d0 w = 0 call eval_Vink_wind(vink_wind) call eval_Nieuwenhuijzen_wind(nieu_wind) call eval_Hamann_wind(hamann_wind) vink_wind = vink_wind*s%x_ctrl(24) nieu_wind = nieu_wind*s%x_ctrl(24) hamann_wind = hamann_wind*s%x_ctrl(24) ! use 1/10 hamann hamann_wind = hamann_wind/10d0 if (T1 < Teff_jump) then ! low T wind w = max(vink_wind, nieu_wind) else ! high T wind alfa = 0d0 if (Xs > 0.7d0) then alfa = 1d0 else if (Xs > 0.4d0 .and. Xs < 0.7d0) then alfa = (Xs - 0.4d0)/0.3d0 end if w = alfa * vink_wind + (1d0-alfa) * hamann_wind end if ierr = 0 contains subroutine eval_Vink_wind(w) real(dp), intent(inout) :: w real(dp) :: alfa, w1, w2, logMdot, dT, vinf_div_vesc ! alfa = 1 for hot side, = 0 for cool side if (T1 > 27500d0) then alfa = 1 else if (T1 < 22500d0) then alfa = 0 else dT = 100d0 if (T1 > Teff_jump + dT) then alfa = 1 else if (T1 < Teff_jump - dT) then alfa = 0 else alfa = (T1 - (Teff_jump - dT)) / (2*dT) end if end if if (alfa > 0) then ! eval hot side wind (eqn 24) vinf_div_vesc = 2.6d0 ! this is the hot side galactic value vinf_div_vesc = vinf_div_vesc*pow_cr(Z_div_Z_solar,0.13d0) ! corrected for Z logMdot = & - 6.697d0 & + 2.194d0*log10_cr(L1/Lsun/1d5) & - 1.313d0*log10_cr(M1/Msun/30) & - 1.226d0*log10_cr(vinf_div_vesc/2d0) & + 0.933d0*log10_cr(T1/4d4) & - 10.92d0*pow2(log10_cr(T1/4d4)) & + 0.85d0*log10_cr(Z_div_Z_solar) w1 = exp10_cr(logMdot) else w1 = 0 end if if (alfa < 1) then ! eval cool side wind (eqn 25) vinf_div_vesc = 1.3d0 ! this is the cool side galactic value vinf_div_vesc = vinf_div_vesc*pow_cr(Z_div_Z_solar,0.13d0) ! corrected for Z logMdot = & - 6.688d0 & + 2.210d0*log10_cr(L1/Lsun/1d5) & - 1.339d0*log10_cr(M1/Msun/30) & - 1.601d0*log10_cr(vinf_div_vesc/2d0) & + 1.07d0*log10_cr(T1/2d4) & + 0.85d0*log10_cr(Z_div_Z_solar) w2 = exp10_cr(logMdot) else w2 = 0 end if w = alfa*w1 + (1 - alfa)*w2 end subroutine eval_Vink_wind subroutine eval_Nieuwenhuijzen_wind(w) ! Nieuwenhuijzen, H.; de Jager, C. 1990, A&A, 231, 134 (eqn 2) real(dp), intent(out) :: w real(dp) :: log10w include 'formats' log10w = -14.02d0 & +1.24d0*log10_cr(L1/Lsun) & +0.16d0*log10_cr(M1/Msun) & +0.81d0*log10_cr(R1/Rsun) & +0.85d0*log10_cr(Z_div_Z_solar) w = exp10_cr(log10w) end subroutine eval_Nieuwenhuijzen_wind subroutine eval_Hamann_wind(w) ! Hamann, W.-R.; Koesterke, L.; Wessolowski, U. 1995, A&A, 299, 151 real(dp), intent(out) :: w real(dp) :: log10w include 'formats' log10w = -11.95d0 & +1.5d0*log10_cr(L1/Lsun) & -2.85d0*Xs & + 0.85d0*log10_cr(Z_div_Z_solar) w = exp10_cr(log10w) end subroutine eval_Hamann_wind end subroutine brott_wind subroutine Dutch_tanh_intrp_wind(id, Lsurf, Msurf, Rsurf, Tsurf, w, ierr) use crlibm_lib type (star_info), pointer :: s integer, intent(in) :: id real(dp), intent(in) :: Lsurf, Msurf, Rsurf, Tsurf ! surface values (cgs) ! NOTE: surface is outermost cell. not necessarily at photosphere. ! NOTE: don't assume that vars are set at this point. ! so if you want values other than those given as args, ! you should use values from s% xh(:,:) and s% xa(:,:) only. ! rather than things like s% Teff or s% lnT(:) which have not been set yet. real(dp), intent(out) :: w ! wind in units of Msun/year (value is >= 0) integer, intent(out) :: ierr real(dp) :: L1, M1, R1, T1, xsurf, etaWR, etaHOT, etaCOOL, Zsolar, Y, Z real(dp) :: log10w, w1, w2, T_high, T_low, alfa, fe call get_star_ptr(id,s,ierr) w = 0 ierr = 0 L1 = Lsurf M1 = Msurf R1 = Rsurf T1 = Tsurf Zsolar = 0.0142d0 !I use the same factor everywhere etaHOT = s% x_ctrl(24) etaCOOL = s% x_ctrl(24) etaWR = s% x_ctrl(24) !Assumes the use of approx21.net !xsurf = s%xa(s%net_iso(iprot), 1)+ s%xa(s%net_iso(ih1), 1) !H mass fraction of the outermost cell, species 1 being neutrons ! Assumes to NOT use approx* family of networks xsurf = s%xa(s%net_iso(ih1), 1) Y = s% xa(s%net_iso(ihe3),1) + s% xa(s%net_iso(ihe4),1) Z = s% Zbase T_high = 11000d0 T_low = 10000d0 if (xsurf .GE. 0.4d0) then if (T1 <= T_low) then print *, 'using: de_Jager, eta=', etaCOOL call eval_de_Jager_wind(w) w = w * etaCOOL else if (T1 >= T_high) then print *, 'using: Vink_tanh, eta=', etaHOT call eval_tanh_intrp_vink_wind(w) w = w * etaHOT else ! transition print *, 'interpolating de Jager and Vink_tanh, eta=', etaHOT call eval_de_Jager_wind(w1) call eval_tanh_intrp_vink_wind(w2) alfa = (T1 - T_low)/(T_high - T_low) w = (1-alfa)*w1 + alfa*w2 w = w * etaHOT end if else !means it's a WR star if(s%x_integer_ctrl(3)==1) then print *, 'using: N&L, eta=', etaWR w = 1d-11 * pow_cr(L1/Lsun,1.29d0) * pow_cr(Y,1.7d0) * sqrt(s% Zbase) w = w * etaWR else if(s% x_integer_ctrl(3) ==2) then ! H free WR from tramper 2016 if(s%net_iso(ife56)==0) then fe = s% Zbase/zsolar ! When we are using basic.net we dont have fe in the net else fe = s% xa(s%net_iso(ife56),1)/(1.1686000000000001D-003) end if w = 10.d0**(-9.2d0) * pow_cr(L1/Lsun,0.85d0) * pow_cr(Y,0.44d0) * & pow_cr(fe,1.3d0) * etaWR end if end if contains subroutine eval_de_Jager_wind(w) ! de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259. real(dp), intent(out) :: w real(dp) :: log10w include 'formats' log10w = 1.769d0*log10_cr(L1/Lsun) - 1.676d0*log10_cr(T1) - 8.158d0 w = exp10_cr(log10w) end subroutine eval_de_Jager_wind subroutine eval_tanh_intrp_vink_wind(w) real(dp), intent(inout) :: w real(dp) :: alfa, w1, w2, w_c, w_h, Teff_jump, logMdot, dT, vinf_div_vesc w = 0 if (T1 > 27500d0) then vinf_div_vesc = 2.6d0 ! this is the hot side galactic value vinf_div_vesc = vinf_div_vesc*pow_cr(Z/Zsolar,0.13d0) ! corrected for Z logMdot = & - 6.697d0 & + 2.194d0*log10_cr(L1/Lsun/1d5) & - 1.313d0*log10_cr(M1/Msun/30) & - 1.226d0*log10_cr(vinf_div_vesc/2d0) & + 0.933d0*log10_cr(T1/4d4) & - 10.92d0*pow2(log10_cr(T1/4d4)) & + 0.85d0*log10_cr(Z/Zsolar) w1 = exp10_cr(logMdot) w = w1 else if (T1 < 22500d0) then vinf_div_vesc = 1.3d0 ! this is the cool side galactic value vinf_div_vesc = vinf_div_vesc*pow_cr(Z/Zsolar,0.13d0) ! corrected for Z logMdot = & - 6.688d0 & + 2.210d0*log10_cr(L1/Lsun/1d5) & - 1.339d0*log10_cr(M1/Msun/30) & - 1.601d0*log10_cr(vinf_div_vesc/2d0) & + 1.07d0*log10_cr(T1/2d4) & + 0.85d0*log10_cr(Z/Zsolar) w2 = exp10_cr(logMdot) w = w2 else ! use Vink et al 2001, eqns 14 and 15 to set "jump" temperature Teff_jump = 1d3*(61.2d0 + 2.59d0*(-13.636d0 + 0.889d0*log10_cr(Z/Zsolar))) dT = 100d0 ! print *, 'Teff_jump[K]', Teff_jump ! first evaluate the mass loss rate one would have if T = 27500, all rest held constant vinf_div_vesc = 2.6d0 ! this is the hot side galactic value vinf_div_vesc = vinf_div_vesc*pow_cr(Z/Zsolar,0.13d0) ! corrected for Z w_h = & - 6.697d0 & + 2.194d0*log10_cr(L1/Lsun/1d5) & - 1.313d0*log10_cr(M1/Msun/30) & - 1.226d0*log10_cr(vinf_div_vesc/2d0) & + 0.933d0*log10_cr(27500/4d4) & - 10.92d0*pow2(log10_cr(27500/4d4)) & + 0.85d0*log10_cr(Z/Zsolar) w_h = exp10_cr(w_h) !then evaluate the mass loss rate one would have if T= 22500, all rest held constant vinf_div_vesc = 1.3d0 ! this is the cool side galactic value vinf_div_vesc = vinf_div_vesc*pow_cr(Z/Zsolar,0.13d0) ! corrected for Z w_c = & - 6.688d0 & + 2.210d0*log10_cr(L1/Lsun/1d5) & - 1.339d0*log10_cr(M1/Msun/30) & - 1.601d0*log10_cr(vinf_div_vesc/2d0) & + 1.07d0*log10_cr(22500/2d4) & + 0.85d0*log10_cr(Z/Zsolar) w_c = exp10_cr(w_c) !then interpolate between this two values w = ((w_h - w_c)/2)*tanh_cr((T1-Teff_jump)/(1d0*dT)) & +(w_c+w_h)/2 end if end subroutine eval_tanh_intrp_vink_wind end subroutine Dutch_tanh_intrp_wind subroutine my_adjust_mdot(id, ierr) use star_def integer, intent(in) :: id integer, intent(out) :: ierr type (star_info), pointer :: s real(dp) :: Lrad_div_Ledd ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return if (s% generations > 2) then write(*,*) "check mdots", s% mstar_dot, s% mstar_dot_old if (abs(s% mstar_dot) > 1.05*abs(s% mstar_dot_old)) then s% mstar_dot = 1.05*s% mstar_dot_old else if (abs(s% mstar_dot) < 0.95*abs(s% mstar_dot_old)) then s% mstar_dot = 0.95*s% mstar_dot_old end if end if end subroutine my_adjust_mdot integer function extras_startup(id, restart, ierr) integer, intent(in) :: id logical, intent(in) :: restart integer, intent(out) :: ierr type (star_info), pointer :: s include 'formats' ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return extras_startup = 0 call system_clock(time0,clock_rate) if (.not. restart) then call alloc_extra_info(s) ! lxtra1 is true while hydro is on s% lxtra1 = .false. ! lxtra2 is true after hydro is turned on for the first time s% lxtra2 = .false. ! lxtra4 is true, then use fixed vsurf boundary conditions s% lxtra4 = .false. ! lxtra5 is true, after a profile is saved when gamma_integral=0 ! One profile is saved during each hydro phase, this is turned ! back to false after a relax s% lxtra5 = .false. ! same as before, but it's kept as .true. after the first time ! gamma_integral=0. Used to count time from first pulse. s% lxtra6 = .false. ! lxtra7 indicates if BB BCs are being used s% lxtra7 = .false. ! lxtra8 is true, after a profile is saved when gamma_integral=0.005 ! Done one time per simulation s% lxtra8 = .false. ! xtra 1 contains the time at which a pulse starts s% xtra1 = -1d0 ! xtra2 contains the estimate for the dynamical time when that happens s% xtra2 = -1d0 ! xtra3 contains gamma_integral in bound regions s% xtra3 = -1d0 ! xtra4 contains time in seconds since first time gamma_integral=0 s% xtra4 = 0d0 ! This is used to assess fallback (see extra history columns) s% xtra10 = 0d0 s% xtra11 = 0d0 s% xtra12 = 0d0 ! Surface T for fixed Tsurf atm option ! this is used when layers at large radii are removed s% xtra13 = 0d0 ! Store tau_factor which is set by remove_surface s% xtra14 = 0d0 ! Store velocity the first time remove_surface is used ! after that point we fix the surface velocity to this s% xtra15 = 0d0 ! ixtra1 counts the steps during which the conditions to relax a model ! are met s% ixtra1 = 0 ! ixtra2 counts the number of times the star has relaxed s% ixtra2 = 0 ! ixtra3 counts the number of steps since last relax s% ixtra3 = 0 ! ixtra4 counts the number of steps since hydro was turned on s% ixtra4 = 0 else ! it is a restart call unpack_extra_info(s) if (s% lxtra1) then call star_read_controls(id, 'inlist_hydro_on', ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return call star_set_u_flag(id, .true., ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return else if (s% lxtra2) then call star_read_controls(id, 'inlist_hydro_off', ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return end if if (s% lxtra2) then call star_read_controls(id, 'inlist_after_first_pulse', ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return end if end if end function extras_startup subroutine extras_after_evolve(id, id_extra, ierr) integer, intent(in) :: id, id_extra integer, intent(out) :: ierr type (star_info), pointer :: s double precision :: dt character (len=strlen) :: test ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return call system_clock(time1,clock_rate) dt = dble(time1 - time0) / clock_rate / 60 write(*,'(/,a50,f12.2,99i10/)') 'runtime (minutes), retries, backups, steps', & dt, s% num_retries, s% num_backups, s% model_number ierr = 0 end subroutine extras_after_evolve ! returns either keep_going, retry, backup, or terminate. integer function extras_check_model(id, id_extra) integer, intent(in) :: id, id_extra integer :: ierr, k real(dp) :: conv_fraction type (star_info), pointer :: s include 'formats' ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return extras_check_model = keep_going if (.not. s% lxtra7) then if (s% u_flag .and. s% dt > 1d0) then if (abs(s% xh(s% i_u, 1) - s% xh_old(s% i_u, 1))> 2d7) then extras_check_model = retry end if end if end if 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 = 23 end function how_many_extra_history_columns subroutine data_for_extra_history_columns(id, id_extra, n, names, vals, ierr) integer, intent(in) :: id, id_extra, n character (len=maxlen_history_column_name) :: names(n) real(dp) :: vals(n), v_esc integer, intent(out) :: ierr type (star_info), pointer :: s integer :: k, k0 ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return names(1) = "relax_count" vals(1) = s% ixtra2 names(2) = "log_R_098" ! Radius of 98% of the mass of the star names(3) = "log_R_095" ! of 95% names(4) = "log_R_090" ! and 90% names(5) = "helium_core_mass" names(6) = "co_core_mass" names(7) = "log_R_100" ! Radius of outermost layer names(8) = "log_R_vesc" ! Radius of layers just below the escape velocity names(9) = "log_R_vesc_098" ! Radius of 98% of the mass just below the esc. vel. names(10) = "log_R_vesc_095" ! of 95% names(11) = "log_R_vesc_090" ! and 90% names(12) = "M_below_vesc" ! Mass of layers below vesc names(13) = "u_flag" ! 1 if hydro is on names(14) = "Mbound_below_10Rsun" ! When layers are above the escape velocity, this indicates which of ! the mass that remains bound has been pushed beyond 10Rsun names(15) = "Mbound_below_20Rsun" ! Same for 20Rsun names(16) = "Mbound_below_30Rsun" ! Same for 30Rsun names(17) = "U_ejecta" ! internal energy of layers above vesc names(18) = "T_ejecta" ! kinetic energy of layers above vesc names(19) = "Omega_ejecta" ! gravitational energy of layers above vesc names(20) = "gamma_integral" vals(20) = s% xtra3 names(21) = "max_velocity" vals(21) = maxval(s% u(1:s% nz)) names(22) = "min_velocity" vals(22) = minval(s% u(1:s% nz)) names(23) = "yr_since_coll" vals(23) = s% xtra4/secyer vals(2) = -100 vals(3) = -100 vals(4) = -100 vals(5) = -100 vals(6) = 0 vals(7) = log10_cr(s% r(1)/Rsun) vals(9) = -100 vals(10) = -100 vals(11) = -100 vals(14) = s% xtra10 vals(15) = s% xtra11 vals(16) = s% xtra12 vals(17) = 0d0 vals(18) = 0d0 vals(19) = 0d0 do k=1, s% nz if (s% q(k) < 0.98 .and. vals(2) == -100) then vals(2) = log10_cr(s% r(k)/Rsun) else if (s% q(k) < 0.95 .and. vals(3) == -100) then vals(3) = log10_cr(s% r(k)/Rsun) else if (s% q(k) < 0.9 .and. vals(4) == -100) then vals(4) = log10_cr(s% r(k)/Rsun) end if if (vals(5) < 0d0 .and. s% X(k) < 0.01) then vals(5) = s% m(k)/Msun end if if (vals(6) <= 0d0 .and. s% Y(k) < 0.01) then vals(6) = s% m(k)/Msun end if end do if (s% u_flag) then do k = s% nz, 1, -1 v_esc = s% x_ctrl(1)*sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k))) if (s% u(k) > v_esc) then exit end if end do if (k==0) then k=1 end if else k=1 end if vals(8) = log10_cr(s% r(k)/Rsun) vals(12) = s% m(k)/Msun ! get internal radii do k0=k, s% nz if (s% q(k0)/s% q(k) < 0.98 .and. vals(9) == -100) then vals(9) = log10_cr(s% r(k0)/Rsun) else if (s% q(k0)/s% q(k) < 0.95 .and. vals(10) == -100) then vals(10) = log10_cr(s% r(k0)/Rsun) else if (s% q(k0)/s% q(k) < 0.9 .and. vals(11) == -100) then vals(11) = log10_cr(s% r(k0)/Rsun) end if end do ! get energies if (k>1) then do k0 = 1, k vals(17) = vals(17) + s% dm(k0)*s% energy(k0) if (s% u_flag) then vals(18) = vals(18) + 0.5d0*s% dm(k0)*s% u(k0)*s% u(k0) end if vals(19) = vals(19) - s% dm(k0)*s% cgrav(k0)*s% m(k0)/s% r(k0) end do end if if (s% u_flag) then vals(13) = 1d0 else vals(13) = 0d0 end if end subroutine data_for_extra_history_columns 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 = 9 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), alpha integer, intent(out) :: ierr type (star_info), pointer :: s integer :: k ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return names(1) = "vesc" names(2) = "v_div_vesc" names(3) = "specific_grav_e" names(4) = "specific_kin_e" names(5) = "specific_thermal_e" names(6) = "total_specific_e" names(7) = "mlt_vc" names(8) = "gradr_start" names(9) = "mlt_vc_start" do k = s% nz, 1, -1 vals(k,1) = sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k))) vals(k,2) = s% u(k)/vals(k,1) vals(k,3) = -s% cgrav(k)*s% m(k)/s% r(k) vals(k,4) = 0.5d0*s% u(k)**2 vals(k,5) = two_thirds*avo*kerg*s% T(k)/(2*s% mu(k)*s% rho(k)) & + crad*pow4(s% T(k))/s% rho(k) vals(k,6) = vals(k,3) + vals(k,4) + vals(k,5) vals(k,7) = s% mlt_vc(k) vals(k,8) = s% gradr_start(k) vals(k,9) = s% mlt_vc_start(k) end do end subroutine data_for_extra_profile_columns integer function extras_start_step(id, id_extra) integer, intent(in) :: id, id_extra integer :: ierr type (star_info), pointer :: s integer :: k, k0, k1, num_pts, species, model_number, num_trace_history_values real(dp) :: v_esc, time, gamma1_integral, integral_norm, tdyn, & max_center_cell_dq, avg_v_div_vesc, energy_removed_layers, dt_next, dt, & total_specific_energy, max_years_for_timestep real(dp), pointer :: q(:), xq(:), xa(:,:), entropy(:) real(dp) :: core_mass, rmax, alfa, log10_r, conv_vel_temp real(dp), pointer :: interp_work(:), conv_vel_interp(:) real(dp) :: cumulative_acoustic_L, cumulative_acoustic_L_center, & cumulative_eps_grav, cumulative_delta_total_energy, cumulative_energy_error, & cumulative_L_center, cumulative_L_surf, cumulative_extra_heating, & cumulative_irradiation_heating, cumulative_Ne22_sedimentation_heating, & cumulative_nuclear_heating, cumulative_non_nuc_neu_cooling, & cumulative_sources_and_sinks logical :: just_did_relax character (len=200) :: fname include 'formats' extras_start_step = terminate ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return !this is used to ensure BCs are right s% use_other_before_struct_burn_mix = .true. ! for some reason setting this in the inlist crashes the initial model s% max_conv_vel_div_csound = 1d0 ! Ignore energy checks before first time hydro is turned on ! otherwise need small steps during core helium burning and it ! slows down the test_suite if (.not. s% lxtra2) then s% cumulative_eps_grav = 0d0 s% cumulative_delta_total_energy = 0d0 s% cumulative_energy_error = 0d0 s% cumulative_L_center = 0d0 s% cumulative_L_surf = 0d0 s% cumulative_extra_heating = 0d0 s% cumulative_irradiation_heating = 0d0 s% cumulative_Ne22_sedimentation_heating = 0d0 s% cumulative_nuclear_heating = 0d0 s% cumulative_non_nuc_neu_cooling = 0d0 s% cumulative_sources_and_sinks = 0d0 write(*,*) & "Setting energy conservation error to zero until hydro is turned on for the first time" end if if (pgstar_flag) then s% job% pgstar_flag = .true. end if just_did_relax = .false. ! get point where v v_esc) then exit end if end do if (k>0 .and. k < s% nz) k = k+1 ! s% xtra3 stores the gamma integral ! during pulses integrate gamma only inside bound regions ! otherwise compute it for the whole star s% xtra3 = -1d0 ! can be adjusted below if nearing breakout s% profile_interval = 10000 if (s% u_flag .and. k > 0 .and. s% xtra1 > 0d0) then !Compute mass that has gone beyond certain radii ("fallback") ! do this for 10, 50 and 100 Rsun if (s% r(k) > 10*Rsun) then do k0 = k, s% nz-1 if (s% r(k0+1) < 10*Rsun) then s% xtra10 = max(s% xtra10, (s% m(k) - s% m(k0))/Msun) exit end if end do end if if (s% r(k) > 20*Rsun) then do k0 = k, s% nz-1 if (s% r(k0+1) < 20*Rsun) then s% xtra11 = max(s% xtra11, (s% m(k) - s% m(k0))/Msun) exit end if end do end if if (s% r(k) > 30*Rsun) then do k0 = k, s% nz-1 if (s% r(k0+1) < 30*Rsun) then s% xtra12 = max(s% xtra12, (s% m(k) - s% m(k0))/Msun) exit end if end do end if ! check energy and average escape velocity of outer layers avg_v_div_vesc = 0d0 energy_removed_layers = 0d0 do k0 = 1, k avg_v_div_vesc = avg_v_div_vesc + s% dm(k0)*s% u(k0)/sqrt(2*s% cgrav(k0)*s% m(k0)/(s% r(k0))) energy_removed_layers = energy_removed_layers + & 0.5d0*s% dm(k0)*s% u(k0)*s% u(k0) - s% dm(k0)*s% cgrav(k0)*s% m(k0)/s% r(k0) end do avg_v_div_vesc = avg_v_div_vesc/(s% m(1) - s% m(k)) ! compute gamma integral in what will remain of the star integral_norm = 0d0 gamma1_integral = 0d0 do k0=k,s% nz integral_norm = integral_norm + s% P(k0)*s% dm(k0)/s% rho(k0) gamma1_integral = gamma1_integral + & (s% gamma1(k0)-4.d0/3.d0)*s% P(k0)*s% dm(k0)/s% rho(k0) end do gamma1_integral = gamma1_integral/max(1d-99,integral_norm) s% xtra3 = gamma1_integral do k1 = 1, s% nz if (s% q(k1) < 0.9d0) then !increase profile density near breakout, check that at q=0.9 velocities are ! above 500 km/s, and that at the surface they're below that. ! for the first pulse, increase profile density from the onset of instability ! to breakout if ((s% u(k1)>5d7 .and. s% u(1)<5d7) & .or. (s% ixtra2 == 0 .and. gamma1_integral < 0d0 .and. s% u(1)<5d7)) then s% profile_interval = 10000 else s% profile_interval = 10000 end if exit end if end do ! To relax the star after a pulse, we check for a series of conditions that ! must apply in layers below q=x_ctrl(7) do k0 = k, s% nz if (s% q(k0) < s% x_ctrl(7)*s% q(k)) then exit end if end do write(*,*) 'Layers above q=', s% q(k), 'will be removed' write(*,*) 'checking for conditions inside q=', s% x_ctrl(7), 'of material that will remain' write(*,*) 'check time left', s% xtra1 + s% xtra2*s% x_ctrl(6) - s% time write(*,*) 'max vel inside fraction of cutoff', maxval(abs(s% u(k0:s% nz)))/1e5 write(*,*) 'max c/cs inside fraction of cutoff', maxval(abs(s% u(k0:s% nz)/s% csound(k0:s% nz))) write(*,*) 'average v/vesc outside cutoff', avg_v_div_vesc write(*,*) 'Kinetic plus potential energy outside cutoff', energy_removed_layers write(*,*) 'mass inside cutoff', k, s% m(k)/Msun write(*,*) 'relax counter', s% ixtra1 write(*,*) 'log max v [km/s]=', safe_log10_cr(maxval(s% u(1:s% nz))/1e5) ! If enough time has passed after a pulse, then turn off Riemann hydro ! also, wait for the inner s% x_ctrl(7) of what would remain of the star to be below a threshold ! in v/cs, and below a threshold in v ! Verify also that the conditions to turn on Riemann hydro are not still satisfied ! For details on all these options check inlist_project ! ignore if s% q(k) < 1d-3, in that case it's very likely a PISN if (s% q(k) > 1d-3) then if (s% lxtra1 .and. s% xtra1 > 0d0 .and. s% time > s% xtra1 + s% xtra2*s% x_ctrl(6)& .and. maxval(abs(s% u(k0:s% nz)))/1e5 < s% x_ctrl(8) & .and. maxval(abs(s% u(k0:s% nz)/s% csound(k0:s% nz))) < s% x_ctrl(9) & .and. log10_cr(s% T(s% nz)) < s% x_ctrl(10) & .and. (.not. (log10_cr(s% T(s% nz)) > s% x_ctrl(2) & .and. gamma1_integral < s% x_ctrl(3))) & .and. (log10_cr(s% power_neutrinos) < s% x_ctrl(16)) & .and. (log10_cr(s% power_nuc_burn) < s% x_ctrl(20))) then s% ixtra1 = s% ixtra1 + 1 else s% ixtra1 = 0 end if else do k0 = k, 1, -1 v_esc = s% x_ctrl(1)*sqrt(2*s% cgrav(k0)*s% m(k0)/(s% r(k0))) if (s% u(k0) < v_esc) then write(*,*) "Likely PISN", s% q(k0), s% u(k0)/v_esc exit end if end do if (k0 == 0) then extras_start_step = terminate write(fname, fmt="(a19)") 'LOGS/pisn_prof.data' call star_write_profile_info(id, fname, id_extra, ierr) write(*,*) "Entire star is expanding above the escape velocity, PISN!" write(*,*) "termination" return end if end if if(s% ixtra1 >= s% x_integer_ctrl(1) .and. s% q(k) > 1d-3) then write(*,*) "Relaxing model to lower mass!" s% ixtra2 = s% ixtra2 + 1 !save a profile just before relaxation write(fname, fmt="(a18,i0.3,a5)") 'LOGS/prerelax_prof', s% ixtra2, '.data' call star_write_profile_info(id, fname, id_extra, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return time = s% time model_number = s% model_number num_trace_history_values = s% num_trace_history_values cumulative_eps_grav = s% cumulative_eps_grav cumulative_delta_total_energy = s% cumulative_delta_total_energy cumulative_energy_error = s% cumulative_energy_error cumulative_L_center = s% cumulative_L_center cumulative_L_surf = s% cumulative_L_surf cumulative_extra_heating = s% cumulative_extra_heating cumulative_irradiation_heating = s% cumulative_irradiation_heating cumulative_Ne22_sedimentation_heating = s% cumulative_Ne22_sedimentation_heating cumulative_nuclear_heating = s% cumulative_nuclear_heating cumulative_non_nuc_neu_cooling = s% cumulative_non_nuc_neu_cooling cumulative_sources_and_sinks = s% cumulative_sources_and_sinks s% lxtra1 = .false. s% lxtra4 = .false. s% lxtra5 = .false. s% lxtra7 = .false. ! Reset masses of "fallback" s% xtra10 = 0d0 s% xtra11 = 0d0 s% xtra12 = 0d0 s% ixtra3 = 0 write(*,*) "removing cells", k, s% m(k)/Msun species = s% species num_pts = s% nz - k + 1 allocate(q(s% nz - k + 1), xq(s% nz - k + 1), xa(species, s% nz - k + 1), & entropy(s% nz - k + 1)) !need to compute cell-centered q for remaining mass xq(1) = s% dq(k)/2/s% q(k) do k0 = 1, s% nz - k xq(1+k0) = xq(1+k0-1) + (s% dq(k+k0) + s% dq(k+k0-1))/s% q(k)/2 end do ! create interpolant for conv_vel, notice that conv_vel is edge valued ! we also interpolate luminosities to get a consistent gradr right after ! relax allocate(interp_work((s% nz - k + 1)*pm_work_size), & conv_vel_interp(4*(s% nz - k + 1)), stat=ierr) do k0 = 1, num_pts conv_vel_interp(4*k0-3) = s% conv_vel(k0+k-1) q(k0) = s% q(k0+k-1)/s% q(k) end do call interp_pm(q, num_pts, conv_vel_interp,& pm_work_size, interp_work, 'conv_vel interpolant', ierr) xa(:,:) = s% xa(:,k:s% nz) entropy(:) = s% entropy(k:s% nz)*avo*kerg max_center_cell_dq = s% max_center_cell_dq s% max_center_cell_dq = s% dq(s% nz) dt = s% dt dt_next = s% dt_next call star_read_controls(id, 'inlist_hydro_off', ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return s% job% pgstar_flag = .false. call star_set_conv_vel_flag(id, .false., ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return call star_set_u_flag(id, .false., ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return call my_before_struct_burn_mix(s% id, s% dt, extras_start_step) s% use_atm_PT_at_center_of_surface_cell = .false. s% have_previous_conv_vel = .false. s% initial_z = 0.02 s% initial_y = 0.28 s% initial_mass = s% m(k)/Msun ! avoid making photos s% photo_interval=100000 call star_change_to_new_net(id, .true., 'basic.net', ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return call star_load_zams(id, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return call star_change_to_new_net(id, .true., s% job% new_net_name, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return max_years_for_timestep = s% max_years_for_timestep s% max_years_for_timestep = 1d0 s% delta_lgL_nuc_limit = -1d0 s% delta_lgL_nuc_hard_limit = -1d0 s% dxdt_nuc_factor = 0d0 s% non_nuc_neu_factor = 0d0 s% eps_nuc_factor = 0d0 s% use_other_before_struct_burn_mix = .false. call star_set_conv_vel_flag(id, .true., ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return s% num_trace_history_values = 0 call star_relax_composition( & id, s% job% num_steps_to_relax_composition, num_pts, species, xa, xq, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return call star_relax_entropy( & id, s% job% max_steps_to_relax_entropy, num_pts, entropy, xq, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return s% dxdt_nuc_factor = 1d0 s% non_nuc_neu_factor = 1d0 s% eps_nuc_factor = 1d0 s% max_years_for_timestep = max_years_for_timestep s% photo_interval=100000 s% dt_next = min(1d5, dt_next) s% dt = min(1d5, dt) s% time = time s% model_number = model_number ! this avoids the history file from being rewritten s% doing_first_model_of_run = .false. ! reload relax counter s% ixtra1 = 0 !take care of convective velocities do k0=1, s% nz call interp_value(q, num_pts, conv_vel_interp, s% q(k0), s% conv_vel(k0), ierr) !avoid extending regions with non-zero conv vel do k=2, num_pts-1 if(s% q(k0) < q(k) .and. s% q(k0) > q(k+1) & .and. (conv_vel_interp(4*k-3)<1d-5 .or. conv_vel_interp(4*(k+1)-3)<1d-5)) then s% conv_vel(k0) = 0d0 exit end if end do s% xh(s% i_ln_cvpv0, k0) = log_cr(s% conv_vel(k0)+s% conv_vel_v0) end do do k0=1, s% nz_old call interp_value(q, num_pts, conv_vel_interp, s% q_old(k0), conv_vel_temp, ierr) !avoid extending regions with non-zero conv vel do k=2, num_pts-1 if(s% q_old(k0) < q(k) .and. s% q_old(k0) > q(k+1) & .and. (conv_vel_interp(4*k-3)<1d-5 .or. conv_vel_interp(4*(k+1)-3)<1d-5)) then conv_vel_temp = 0d0 exit end if end do s% xh_old(s% i_ln_cvpv0, k0) = log_cr(conv_vel_temp+s% conv_vel_v0) end do s% generations = 1 s% max_center_cell_dq = max_center_cell_dq just_did_relax = .true. !save a profile and a model just after relaxation write(fname, fmt="(a19,i0.3,a5)") 'LOGS/postrelax_prof', s% ixtra2, '.data' call star_write_profile_info(id, fname, id_extra, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return write(fname, fmt="(a20,i0.3,a4)") 'LOGS/postrelax_model', s% ixtra2, '.mod' call star_write_model(id, fname, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return s% cumulative_eps_grav = cumulative_eps_grav s% cumulative_delta_total_energy = cumulative_delta_total_energy s% cumulative_energy_error = cumulative_energy_error s% cumulative_L_center = cumulative_L_center s% cumulative_L_surf = cumulative_L_surf s% cumulative_extra_heating = cumulative_extra_heating s% cumulative_irradiation_heating = cumulative_irradiation_heating s% cumulative_Ne22_sedimentation_heating = cumulative_Ne22_sedimentation_heating s% cumulative_nuclear_heating = cumulative_nuclear_heating s% cumulative_non_nuc_neu_cooling = cumulative_non_nuc_neu_cooling s% cumulative_sources_and_sinks = cumulative_sources_and_sinks s% num_trace_history_values = num_trace_history_values deallocate(q, xq, xa, entropy, conv_vel_interp, interp_work) s% use_other_before_struct_burn_mix = .true. end if end if ! turn on Riemman hydro if close to instability integral_norm = 0d0 gamma1_integral = 0d0 do k=1,s% nz integral_norm = integral_norm + s% P(k)*s% dm(k)/s% rho(k) gamma1_integral = gamma1_integral + & (s% gamma1(k)-4.d0/3.d0)*s% P(k)*s% dm(k)/s% rho(k) end do gamma1_integral = gamma1_integral/max(1d-99,integral_norm) if (s% xtra3 == -1d0) then ! if xtra3 is different from -1d0 it means it was computed for the bound material, ! use that instead s% xtra3 = gamma1_integral end if ! Save profile if gamma1_integral becomes negative if (gamma1_integral < 0d0 .and. .not. s% lxtra5) then s% lxtra5 = .true. s% lxtra6 = .true. write(fname, fmt="(a16,i0.3,a5)") 'LOGS/gamma1_prof', s% ixtra2+1, '.data' call star_write_profile_info(id, fname, id_extra, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return end if ! Save profile shortly before pulses begin if (gamma1_integral < 0.005d0 .and. .not. s% lxtra8) then s% lxtra8 = .true. write(fname, fmt="(a15)") 'LOGS/preSN.data' call star_write_profile_info(id, fname, id_extra, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return end if write(*,*) "check gamma integral", gamma1_integral if (s% ixtra3 > 10 .and. .not. just_did_relax .and. .not. s% lxtra1 & .and. ((log10_cr(s% T(s% nz)) > s% x_ctrl(2) .and. & gamma1_integral < s% x_ctrl(3)) .or. log10_cr(s% T(s% nz)) > s% x_ctrl(10) & .or. log10_cr(s% power_neutrinos) > s% x_ctrl(11) & .or. (s% lxtra2 .and. log10_cr(s% power_neutrinos) > s% x_ctrl(17)))) then write(*,*) "Turning on Riemann hydro!" call star_read_controls(id, 'inlist_hydro_on', ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return call star_read_controls(id, 'inlist_after_first_pulse', ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return !! sponge fixes velocities to zero just after hydro is turned on !! this is removed in extras_finish_step s% dt_next = min(1d0,s% dt_next) s% dt = min(1d0,s% dt) write(fname, fmt="(a18,i0.3,a5)") 'LOGS/prehydro_prof', s% ixtra2+1, '.data' call star_write_profile_info(id, fname, id_extra, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return write(fname, fmt="(a19,i0.3,a4)") 'LOGS/prehydro_model', s% ixtra2+1, '.mod' call star_write_model(id, fname, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return s% job% pgstar_flag = .false. call star_set_u_flag(id, .true., ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return s% lxtra1 = .true. s% lxtra2 = .true. s% xtra1 = -1 s% ixtra4 = 0 !force initial velocities to zero to prevent issues in outer layers s% xh(s% i_u,:) = 0d0 s% xh_old(s% i_u,:) = 0d0 s% xh_older(s% i_u,:) = 0d0 s% generations = 1 end if ! check time to relax model ! the star is considered to be in a pulse if the velocity is high enough ! only check velocity in x_ctrl(7) of the mass of the star, this avoids ! considering as a pulse cases where the surface goes a bit nuts do k0 = 1, s% nz if (s% q(k0) < s% x_ctrl(7)) then exit end if end do if (s% lxtra1 .and. s% xtra1 < 0d0) then if ((maxval(abs(s% xh(s% i_u, k0:s% nz)))/1e5 > s% x_ctrl(4) .or. gamma1_integral < 1d-3)) then core_mass = s% m(1) if(s% o_core_mass > 0d0) then core_mass = s% o_core_mass*Msun else if(s% c_core_mass > 0d0) then core_mass = s% c_core_mass*Msun else if (s% he_core_mass > 0d0) then core_mass = s% he_core_mass*Msun end if do k=1, s% nz if (s% m(k)/core_mass < s% x_ctrl(5)) then exit end if end do tdyn = 1/sqrt(s% m(k)/(4d0/3d0*pi*s% r(k)**3)*standard_cgrav) s% xtra1 = s% time s% xtra2 = tdyn write(*,*) "reached high velocities", maxval(abs(s% xh(s% i_u, k0:s% nz)))/1e5, & "stopping in at least", s% x_ctrl(6)*tdyn, "seconds" end if end if !Always call this at the end to ensure we are using the BCs we want call my_before_struct_burn_mix(s% id, s% dt, extras_start_step) extras_start_step = keep_going end function extras_start_step subroutine my_before_struct_burn_mix(id, dt, res) use const_def, only: dp use star_def integer, intent(in) :: id real(dp), intent(in) :: dt integer, intent(out) :: res ! keep_going, redo, retry, backup, terminate real(dp) :: power_photo, v_esc integer :: ierr, k type (star_info), pointer :: s include 'formats' ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return !read proper options when hydro is on/off !do this to ensure proper behaviour of retries/backups if(s% u_flag) then call star_read_controls(id, 'inlist_hydro_on', ierr) if (s% xtra1 > 0d0) then s% v_drag = 1d5*s% x_ctrl(12) s% v_drag_factor = 1d0 if (.not. s% lxtra7) then s% max_timestep = s% x_ctrl(23) else s% max_timestep = 1d99 end if s% max_q_for_convection_with_hydro_on = 1d10 do k = s% nz, 1, -1 v_esc = s% x_ctrl(1)*sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k))) if (s% u(k) > 4*v_esc) then exit end if end do if (k > 1) then s% max_q_for_convection_with_hydro_on = s% q(k) end if else !this prevents surface from going mad right when hydro is turned on !right after a relax s% v_drag = 0d0 s% v_drag_factor = 1d0 s% max_timestep = 1d99 s% max_q_for_convection_with_hydro_on = 0.999d0 end if do k=1, s% nz s% xh(s% i_u,k) = min(s% xh(s% i_u,k), 1d5*s% x_ctrl(12)) s% u(k) = s% xh(s% i_u,k) end do ! use fixed_vsurf if surface v remains too high if (s% xh(s% i_u,1) >= 1d5*s% x_ctrl(12)) then s% use_fixed_vsurf_outer_BC = .true. s% use_momentum_outer_BC = .false. s% fixed_vsurf = s% x_ctrl(12)*1d5 else if (s% lxtra7) then! .and. s% xh(s% i_u,1) <= 5d2*1d5) then s% use_fixed_vsurf_outer_BC = .true. s% use_momentum_outer_BC = .false. !s% fixed_vsurf = 5d2*1d5 s% fixed_vsurf = s% xtra15 else s% use_fixed_vsurf_outer_BC = .false. s% use_momentum_outer_BC = .true. end if end if if (s% xtra1 < 0d0) then !this allows the outermost layers to merge right after hydro is turned on !it helps prevent a tiny very outermost layer from going mad s% merge_amr_max_abs_du_div_cs = 1d99 else s% merge_amr_max_abs_du_div_cs = 0.1d0 end if else s% max_timestep = 1d99 call star_read_controls(id, 'inlist_hydro_off', ierr) end if !ensure we are using the proper BCs and options every time before struct_burn_mix if (s% lxtra7) then s% use_atm_PT_at_center_of_surface_cell = .true. s% which_atm_option = "fixed_Tsurf" s% atm_fixed_Tsurf = s% xtra13 s% tau_factor = s% xtra14 s% force_tau_factor = s% xtra14 else if (s% ixtra2 > 0 .and. (s% xh(s% i_lnR, 1)-log_cr(Rsun))/log_cr(10d0) > 2d0) then s% use_atm_PT_at_center_of_surface_cell = .true. else s% use_atm_PT_at_center_of_surface_cell = .false. end if s% which_atm_option = "simple_photosphere" s% tau_factor = 1d0 s% force_tau_factor = 1d0 end if !ignore L_nuc limit if L_phot is too high or if we just did a relax !also ignore if timesteps are too small !(ixtra3 is set to zero right after a relax) power_photo = dot_product(s% dm(1:s% nz), s% eps_nuc_categories(iphoto,1:s% nz))/Lsun if (s% ixtra3 == 0 .or. log10_cr(abs(power_photo)) > s% x_ctrl(22) .or. s% dt < 10d0) then s% delta_lgL_nuc_limit = -1d0 s% delta_lgL_nuc_hard_limit = -1d0 else s% delta_lgL_nuc_limit = s% x_ctrl(21) s% delta_lgL_nuc_hard_limit = 2d0*s% x_ctrl(21) end if if (s% use_fixed_vsurf_outer_BC) then write(*,*) "Using fixed_vsurf", s% fixed_vsurf/1e5 end if !ignore winds if neutrino luminosity is too high or for a few steps after !a relax if(s% ixtra3 < 50 .or. log10_cr(s% power_neutrinos) > s% x_ctrl(19) .or. s% u_flag) then s% use_other_wind = .false. else s% use_other_wind = .true. end if ! reading inlists can turn this flag off for some reason s% use_other_before_struct_burn_mix = .true. res = keep_going end subroutine my_before_struct_burn_mix subroutine null_binary_controls(id, binary_id, ierr) integer, intent(in) :: id, binary_id integer, intent(out) :: ierr ierr = 0 end subroutine null_binary_controls ! returns either keep_going or terminate. integer function extras_finish_step(id, id_extra) use run_star_support integer, intent(in) :: id, id_extra integer :: ierr,k real(dp) :: max_vel_inside, v_esc type (star_info), pointer :: s include 'formats' ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return extras_finish_step = keep_going if(abs(s% xtra25_old - s% center_h1) > s% x_ctrl(18)) then s% dt_next = min(s% dt_next, s% dt * s% min_timestep_factor) write(*,*) "reducing dt due to large change in central hydrogen" else if(abs(s% xtra26_old - s% center_he4) > s% x_ctrl(18) .and. s% center_c12 > 1d-4) then s% dt_next = min(s% dt_next, s% dt * s% min_timestep_factor) write(*,*) "reducing dt due to large change in central helium" else if(abs(s% xtra27_old - s% center_c12) > s% x_ctrl(18)) then s% dt_next = min(s% dt_next, s% dt * s% min_timestep_factor) write(*,*) "reducing dt due to large change in central carbon" else if(abs(s% xtra28_old - s% center_o16) > s% x_ctrl(18)) then s% dt_next = min(s% dt_next, s% dt * s% min_timestep_factor) write(*,*) "reducing dt due to large change in central oxygen" end if s% xtra25 = s% center_h1 s% xtra26 = s% center_he4 s% xtra27 = s% center_c12 s% xtra28 = s% center_o16 if (s% u_flag) then s% xtra29 = s% u(1) else s% xtra29 = 0d0 end if !count time since first collapse if (s% lxtra6) then s% xtra4 = s% xtra4 + s% dt end if do k=1, s% nz if (s% conv_vel(k) < 1d-5 .and. s% mlt_vc(k) <= 0d0) then s% xh(s% i_ln_cvpv0, k) = 0d0 s% conv_vel(k) = 0d0 end if end do if (s% x_logical_ctrl(1) .and. s% u_flag & .and. s% log_surface_radius > 4d0) then do k = 1, s% nz if (log10_cr(s% r(k)/Rsun) < 4d0) then exit end if end do s% which_atm_option = "fixed_Tsurf" s% atm_fixed_Tsurf = s% T(k) s% xtra13 = s% T(k) v_esc = sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k))) write(*,*) "Check v_esc limit", 10*v_esc/1e5 s% xtra15 = max(s% u(k), 10*v_esc) s% lxtra7 = .true. write(*,*) "Removing surface layers", k, s% m(k)/Msun call star_remove_surface_at_cell_k(s% id, k, ierr) if (dbg) write(*,*) "check ierr", ierr if (ierr /= 0) return s% xtra14 = s% tau_factor end if s% ixtra3 = s% ixtra3 + 1 s% ixtra4 = s% ixtra4 + 1 if (s% ixtra2 == 1 .and. s% xtra3 < 0d0 .and. s% x_logical_ctrl(2)) then !for the test_suite, terminate at the onset of the second pulse extras_finish_step = terminate write(*,*) "Terminate due to onset of second pulse" end if if (extras_finish_step == terminate) s% termination_code = t_extras_finish_step end function extras_finish_step 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 call move_int(redo_count) num_ints = i i = 0 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) double precision :: 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 subroutine my_neu( & id, k, T, log10_T, Rho, log10_Rho, abar, zbar, z2bar, log10_Tlim, flags, & loss, sources, ierr) !use neu_lib, only: neu_get use neu_def integer, intent(in) :: id ! id for star integer, intent(in) :: k ! cell number or 0 if not for a particular cell real(dp), intent(in) :: T ! temperature real(dp), intent(in) :: log10_T ! log10 of temperature real(dp), intent(in) :: Rho ! density real(dp), intent(in) :: log10_Rho ! log10 of density real(dp), intent(in) :: abar ! mean atomic weight real(dp), intent(in) :: zbar ! mean charge real(dp), intent(in) :: z2bar ! mean charge squared real(dp), intent(in) :: log10_Tlim logical, intent(inout) :: flags(num_neu_types) ! true if should include the type of loss real(dp), intent(inout) :: loss(num_neu_rvs) ! total from all sources real(dp), intent(inout) :: sources(num_neu_types, num_neu_rvs) integer, intent(out) :: ierr integer :: i real(dp) :: scale(num_neu_types) type(star_info),pointer :: s ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return call my_neutrinos( & T, log10_T, Rho, log10_Rho, abar, zbar, z2bar, log10_Tlim, flags, & loss, sources, s%x_ctrl(26), s% x_ctrl(27), ierr) scale=(/0.1d0,0.05d0,0.01d0,0.1d0,0.1d0/) scale = 1d0 - scale * s%x_ctrl(25) loss=loss*scale do i=1,num_neu_rvs sources(:,i) = sources(:,i) * scale end do end subroutine my_neu ! These are copied from mod_neu.f90 as we can't change the weinberg angle from ! the public interface real(dp) function ifermi12(f) !..this routine applies a rational function expansion to get the inverse !..fermi-dirac integral of order 1/2 when it is equal to f. !..maximum error is 4.19d-9. reference: antia apjs 84,101 1993 !..declare integer i,m1,k1,m2,k2 real(dp) f,an,a1(12),b1(12),a2(12),b2(12),rn,den,ff !..load the coefficients of the expansion data an,m1,k1,m2,k2 /0.5d0, 4, 3, 6, 5/ data (a1(i),i=1,5)/ 1.999266880833d4, 5.702479099336d3, & 6.610132843877d2, 3.818838129486d1, & 1.0d0/ data (b1(i),i=1,4)/ 1.771804140488d4, -2.014785161019d3, & 9.130355392717d1, -1.670718177489d0/ data (a2(i),i=1,7)/-1.277060388085d-2, 7.187946804945d-2, & -4.262314235106d-1, 4.997559426872d-1, & -1.285579118012d0, -3.930805454272d-1, & 1.0d0/ data (b2(i),i=1,6)/-9.745794806288d-3, 5.485432756838d-2, & -3.299466243260d-1, 4.077841975923d-1, & -1.145531476975d0, -6.067091689181d-2/ if (f .lt. 4.0d0) then rn = f + a1(m1) do i=m1-1,1,-1 rn = rn*f + a1(i) enddo den = b1(k1+1) do i=k1,1,-1 den = den*f + b1(i) enddo ifermi12 = log_cr(f * rn/den) else ff = 1.0d0/pow_cr(f,1.0d0/(1.0d0 + an)) rn = ff + a2(m2) do i=m2-1,1,-1 rn = rn*ff + a2(i) enddo den = b2(k2+1) do i=k2,1,-1 den = den*ff + b2(i) enddo ifermi12 = rn/(den*ff) end if return end function ifermi12 real(dp) function zfermim12(x) !..this routine applies a rational function expansion to get the fermi-dirac !..integral of order -1/2 evaluated at x. maximum error is 1.23d-12. !..reference: antia apjs 84,101 1993 !..declare integer i,m1,k1,m2,k2 real(dp) x,an,a1(12),b1(12),a2(12),b2(12),rn,den,xx !..load the coefficients of the expansion data an,m1,k1,m2,k2 /-0.5d0, 7, 7, 11, 11/ data (a1(i),i=1,8)/ 1.71446374704454d7, 3.88148302324068d7, & 3.16743385304962d7, 1.14587609192151d7, & 1.83696370756153d6, 1.14980998186874d5, & 1.98276889924768d3, 1.0d0/ data (b1(i),i=1,8)/ 9.67282587452899d6, 2.87386436731785d7, & 3.26070130734158d7, 1.77657027846367d7, & 4.81648022267831d6, 6.13709569333207d5, & 3.13595854332114d4, 4.35061725080755d2/ data (a2(i),i=1,12)/-4.46620341924942d-15, -1.58654991146236d-12, & -4.44467627042232d-10, -6.84738791621745d-8, & -6.64932238528105d-6, -3.69976170193942d-4, & -1.12295393687006d-2, -1.60926102124442d-1, & -8.52408612877447d-1, -7.45519953763928d-1, & 2.98435207466372d0, 1.0d0/ data (b2(i),i=1,12)/-2.23310170962369d-15, -7.94193282071464d-13, & -2.22564376956228d-10, -3.43299431079845d-8, & -3.33919612678907d-6, -1.86432212187088d-4, & -5.69764436880529d-3, -8.34904593067194d-2, & -4.78770844009440d-1, -4.99759250374148d-1, & 1.86795964993052d0, 4.16485970495288d-1/ if (x .lt. 2.0d0) then xx = exp_cr(x) rn = xx + a1(m1) do i=m1-1,1,-1 rn = rn*xx + a1(i) enddo den = b1(k1+1) do i=k1,1,-1 den = den*xx + b1(i) enddo zfermim12 = xx * rn/den !.. else xx = 1.0d0/(x*x) rn = xx + a2(m2) do i=m2-1,1,-1 rn = rn*xx + a2(i) enddo den = b2(k2+1) do i=k2,1,-1 den = den*xx + b2(i) enddo zfermim12 = sqrt(x)*rn/den end if return end function zfermim12 subroutine my_neutrinos(T, logT, Rho, logRho, abar, zbar, z2bar, log10_Tlim, & flags, loss, sources, xnufam, theta, info) use utils_lib, only: is_bad, mesa_error use neu_def !..this routine computes neutrino losses from the analytic fits of !..itoh et al. apjs 102, 411, 1996, and also returns their derivatives. ! provide T or logT or both (the code needs both, so pass 'em if you've got 'em!) ! same for Rho and logRho real(dp), intent(in) :: T ! temperature real(dp), intent(in) :: logT ! log10 of temperature real(dp), intent(in) :: Rho ! density real(dp), intent(in) :: logRho ! log10 of density real(dp), intent(in) :: abar ! mean atomic weight real(dp), intent(in) :: zbar ! mean charge real(dp), intent(in) :: z2bar ! mean charge squared real(dp), intent(in) :: log10_Tlim ! start to cutoff at this temperature logical, intent(in) :: flags(num_neu_types) ! true if should include the type ! in most cases for stellar evolution, you may want to include brem, plas, pair, and phot ! but skip reco. it is fairly expensive to compute and typically makes only a small contribution. real(dp), intent(inout) :: loss(num_neu_rvs) ! total from all sources real(dp), intent(inout) :: sources(num_neu_types, num_neu_rvs) integer, intent(out) :: info ! 0 means AOK. !..various numerical constants real(dp), parameter :: con1 = 1.0d0/5.9302d0 real(dp), parameter :: sixth = 1.0d0/6.0d0 real(dp), parameter :: iln10 = 4.342944819032518d-1 !..theta is sin**2(theta_weinberg) = 0.2319 plus/minus 0.00005 (1996) !..xnufam is the number of neutrino flavors = 3.02 plus/minus 0.005 (1998) !..change theta and xnufam if need be, and the changes will automatically !..propagate through the routines. cv and ca are the vector and axial currents. real(dp),intent(in) :: theta != 0.2319d0 real(dp),intent(in) :: xnufam != 3.0d0 real(dp) :: cv != 0.5d0 + 2.0d0 * theta real(dp) :: cvp != 1.0d0 - cv real(dp) :: ca != 0.5d0 real(dp) :: cap != 1.0d0 - ca real(dp) :: tfac1 != cv*cv + ca*ca + (xnufam-1.0d0) * (cvp*cvp+cap*cap) real(dp) :: tfac2 != cv*cv - ca*ca + (xnufam-1.0d0) * (cvp*cvp - cap*cap) real(dp) :: tfac3 != tfac2/tfac1 real(dp) :: tfac4 != 0.5d0 * tfac1 real(dp) :: tfac5 != 0.5d0 * tfac2 real(dp) :: tfac6 != cv*cv + 1.5d0*ca*ca + (xnufam - 1.0d0)*(cvp*cvp + 1.5d0*cap*cap) !..local variables real(dp) :: temp,logtemp,den,logden,tcutoff_factor,fac1,fac2,fac3 real(dp) :: ye,deni,tempi,abari,zbari real(dp) :: t9,xl,xldt,xlp5,xl2,xl3,xl4,xl5,xl6,xl7,xl8,xl9,xlmp5,xlm1,xlm2,xlm3,xlm4 real(dp) :: snu, dsnudt, dsnudd, dsnuda, dsnudz real(dp) spair,spairdt,spairdd,spairda,spairdz, & splas,splasdt,splasdd,splasda,splasdz, & sphot,sphotdt,sphotdd,sphotda,sphotdz, & sbrem,sbremdt,sbremdd,sbremda,sbremdz, & sreco,srecodt,srecodd,srecoda,srecodz real(dp) :: a0,a1,a2,a3,b1,b2,c00,c01,c02,c03,c04,c05,c06,xlnt,cc,den6,tfermi, & c10,c11,c12,c13,c14,c15,c16,c20,c21,c22,c23,c24, & c25,c26,dd00,dd01,dd02,dd03,dd04,dd05,dd11,dd12, & dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,b,c,d,f0,f1,f2,f3,z, & dum,dumdt,dumdd,dumda,dumdz, & gum,gumdt,gumdd,gumda,gumdz real(dp) :: rm,rmdd,rmda,rmdz,rmi,gl,gldt, & zeta,zetadt,zetadd,zetada,zetadz,zeta2,zeta3, & xnum,xnumdt,xnumdd,xnumda,xnumdz, & xden,xdendt,xdendd,xdenda,xdendz !..photo real(dp) tau,taudt,cos1,cos2,cos3,cos4,cos5,sin1,sin2, & sin3,sin4,sin5,last,xast, & fphot,fphotdt,fphotdd,fphotda,fphotdz, & qphot,qphotdt,qphotdd,qphotda,qphotdz !..brem real(dp) t8,t812,t832,t82,t83,t85,t86,t8m1,t8m2,t8m3,t8m5,t8m6, & eta,etadt,etadd,etada,etadz,etam1,etam2,etam3, & fbrem,fbremdt,fbremdd,fbremda,fbremdz, & gbrem,gbremdt,gbremdd,gbremda,gbremdz, & u,gm1,gm2,gm13,gm23,gm43,gm53,v,w,fb,gt,gb, & fliq,fliqdt,fliqdd,fliqda,fliqdz, & gliq,gliqdt,gliqdd,gliqda,gliqdz !..plasma real(dp) gl2,gl2dt,gl2dd,gl2da,gl2dz,gl12,gl32,gl72,gl6, & ft,ftdt,ftdd,ftda,ftdz,fl,fldt,fldd,flda,fldz, & fxy,fxydt,fxydd,fxyda,fxydz ! pair real(dp) :: fpair,fpairdt,fpairdd,fpairda,fpairdz, & qpair,qpairdt,qpairdd,qpairda,qpairdz ! recombination real(dp) :: nu,nudt,nudd,nuda,nudz, & nu2,nu3,bigj,bigjdt,bigjdd,bigjda,bigjdz cv = 0.5d0 + 2.0d0 * theta cvp = 1.0d0 - cv ca = 0.5d0 cap = 1.0d0 - ca tfac1 = cv*cv + ca*ca + (xnufam-1.0d0) * (cvp*cvp+cap*cap) tfac2 = cv*cv - ca*ca + (xnufam-1.0d0) * (cvp*cvp - cap*cap) tfac3 = tfac2/tfac1 tfac4 = 0.5d0 * tfac1 tfac5 = 0.5d0 * tfac2 tfac6 = cv*cv + 1.5d0*ca*ca + (xnufam - 1.0d0)*(cvp*cvp + 1.5d0*cap*cap) info = 0 if ((T /= arg_not_provided .and. T .le. Tmin_neu) .or. & (logT /= arg_not_provided .and. logT .le. log10Tmin_neu)) then loss = 0d0 sources(1:num_neu_types, 1:num_neu_rvs) = 0d0 return end if if (T == arg_not_provided) then if (logT == arg_not_provided) then info = -2 return end if temp = exp10_cr(logT) else temp = T end if if (T <= 0) then info = -1 return end if if (logT == arg_not_provided) then logtemp = log10_cr(T) else logtemp = logT end if if (logtemp > 20) then info = -1 return end if if (Rho == arg_not_provided) then if (logRho == arg_not_provided) then info = -3 return end if den = exp10_cr(logRho) else den = Rho end if if (Rho <= 0) then info = -1 return end if if (logRho == arg_not_provided) then logden = log10_cr(Rho) else logden = logRho end if if (logden > 20) then info = -1 return end if fac1 = 5.0d0 * pi / 3.0d0 fac2 = 10.0d0 * pi fac3 = pi / 5.0d0 !..initialize spair = 0.0d0 spairdt = 0.0d0 spairdd = 0.0d0 spairda = 0.0d0 spairdz = 0.0d0 splas = 0.0d0 splasdt = 0.0d0 splasdd = 0.0d0 splasda = 0.0d0 splasdz = 0.0d0 sphot = 0.0d0 sphotdt = 0.0d0 sphotdd = 0.0d0 sphotda = 0.0d0 sphotdz = 0.0d0 sbrem = 0.0d0 sbremdt = 0.0d0 sbremdd = 0.0d0 sbremda = 0.0d0 sbremdz = 0.0d0 sreco = 0.0d0 srecodt = 0.0d0 srecodd = 0.0d0 srecoda = 0.0d0 srecodz = 0.0d0 snu = 0.0d0 dsnudt = 0.0d0 dsnudd = 0.0d0 dsnuda = 0.0d0 dsnudz = 0.0d0 !..to avoid lots of divisions deni = 1.0d0 / den tempi = 1.0d0 / temp abari = 1.0d0 / abar zbari = 1.0d0 / zbar !..some composition variables ye = zbar * abari !xmue = abar * zbari !..some frequent factors t9 = temp * 1.0d-9 xl = t9 * con1 xldt = 1.0d-9 * con1 xlp5 = sqrt(xl) xl2 = xl * xl xl3 = xl2 * xl xl4 = xl3 * xl xl5 = xl4 * xl xl6 = xl5 * xl xl7 = xl6 * xl xl8 = xl7 * xl xl9 = xl8 * xl xlmp5 = 1.0d0 / xlp5 xlm1 = 1.0d0 / xl xlm2 = xlm1*xlm1 xlm3 = xlm1*xlm2 xlm4 = xlm1*xlm3 rm = den*ye rmdd = ye rmda = -rm*abari rmdz = den*abari rmi = 1.0d0/rm a0 = rm * 1.0d-9 a1 = pow_cr(a0,one_third) zeta = a1 * xlm1 zetadt = -a1 * xlm2 * xldt a2 = one_third * a1*rmi * xlm1 zetadd = a2 * rmdd zetada = a2 * rmda zetadz = a2 * rmdz zeta2 = zeta * zeta zeta3 = zeta2 * zeta !..do the requested types if (flags(pair_neu_type)) call pair_neu if (flags(plas_neu_type)) call plas_neu if (flags(phot_neu_type)) call phot_neu if (flags(brem_neu_type)) call brem_neu if (flags(reco_neu_type)) call reco_neu !..convert from erg/cm^3/s to erg/g/s !..comment these out to duplicate the itoh et al plots spair = spair * deni spairdt = spairdt * deni spairdd = spairdd * deni - spair * deni spairda = spairda * deni spairdz = spairdz * deni splas = splas * deni splasdt = splasdt * deni splasdd = splasdd * deni - splas * deni splasda = splasda * deni splasdz = splasdz * deni sphot = sphot * deni sphotdt = sphotdt * deni sphotdd = sphotdd * deni - sphot * deni sphotda = sphotda * deni sphotdz = sphotdz * deni sbrem = sbrem * deni sbremdt = sbremdt * deni sbremdd = sbremdd * deni - sbrem * deni sbremda = sbremda * deni sbremdz = sbremdz * deni sreco = sreco * deni srecodt = srecodt * deni srecodd = srecodd * deni - sreco * deni srecoda = srecoda * deni srecodz = srecodz * deni !..calculate temperature cutoff factor if (logtemp <= log10_Tlim .and. log10_Tlim > log10Tmin_neu) then tcutoff_factor = 0.5d0* & (1 - cospi_cr((logtemp - log10Tmin_neu)/(log10_Tlim - log10Tmin_neu))) splas = tcutoff_factor * splas splasdt = tcutoff_factor * splasdt splasdd = tcutoff_factor * splasdd splasda = tcutoff_factor * splasda splasdz = tcutoff_factor * splasdz spair = tcutoff_factor * spair spairdt = tcutoff_factor * spairdt spairdd = tcutoff_factor * spairdd spairda = tcutoff_factor * spairda spairdz = tcutoff_factor * spairdz sphot = tcutoff_factor * sphot sphotdt = tcutoff_factor * sphotdt sphotdd = tcutoff_factor * sphotdd sphotda = tcutoff_factor * sphotda sphotdz = tcutoff_factor * sphotdz sbrem = tcutoff_factor * sbrem sbremdt = tcutoff_factor * sbremdt sbremdd = tcutoff_factor * sbremdd sbremda = tcutoff_factor * sbremda sbremdz = tcutoff_factor * sbremdz sreco = tcutoff_factor * sreco srecodt = tcutoff_factor * srecodt srecodd = tcutoff_factor * srecodd srecoda = tcutoff_factor * srecoda srecodz = tcutoff_factor * srecodz end if !..the total neutrino loss rate snu = splas + spair + sphot + sbrem + sreco dsnudt = splasdt + spairdt + sphotdt + sbremdt + srecodt dsnudd = splasdd + spairdd + sphotdd + sbremdd + srecodd dsnuda = splasda + spairda + sphotda + sbremda + srecoda dsnudz = splasdz + spairdz + sphotdz + sbremdz + srecodz if (is_bad(snu)) then info = -1 return end if !..packup the results loss(1) = snu loss(2) = dsnudt loss(3) = dsnudd loss(4) = dsnuda loss(5) = dsnudz call store(pair_neu_type, spair, spairdt, spairdd, spairda, spairdz) call store(plas_neu_type, splas, splasdt, splasdd, splasda, splasdz) call store(phot_neu_type, sphot, sphotdt, sphotdd, sphotda, sphotdz) call store(brem_neu_type, sbrem, sbremdt, sbremdd, sbremda, sbremdz) call store(reco_neu_type, sreco, srecodt, srecodd, srecoda, srecodz) contains subroutine store(neu_type, s, sdt, sdd, sda, sdz) integer, intent(in) :: neu_type real(dp), intent(in) :: s, sdt, sdd, sda, sdz sources(neu_type,1) = s sources(neu_type,2) = sdt sources(neu_type,3) = sdd sources(neu_type,4) = sda sources(neu_type,5) = sdz end subroutine store subroutine phot_neu !..photoneutrino process section !..for reactions like e- + gamma => e- + nu_e + nubar_e !.. e+ + gamma => e+ + nu_e + nubar_e !..equation 3.8 for tau, equation 3.6 for cc, !..and table 2 written out for speed if (temp .ge. 1.0d7 .and. temp .lt. 1.0d8) then tau = logtemp - 7d0 cc = 0.5654d0 + tau c00 = 1.008d11 c01 = 0.0d0 c02 = 0.0d0 c03 = 0.0d0 c04 = 0.0d0 c05 = 0.0d0 c06 = 0.0d0 c10 = 8.156d10 c11 = 9.728d8 c12 = -3.806d9 c13 = -4.384d9 c14 = -5.774d9 c15 = -5.249d9 c16 = -5.153d9 c20 = 1.067d11 c21 = -9.782d9 c22 = -7.193d9 c23 = -6.936d9 c24 = -6.893d9 c25 = -7.041d9 c26 = -7.193d9 dd01 = 0.0d0 dd02 = 0.0d0 dd03 = 0.0d0 dd04 = 0.0d0 dd05 = 0.0d0 dd11 = -1.879d10 dd12 = -9.667d9 dd13 = -5.602d9 dd14 = -3.370d9 dd15 = -1.825d9 dd21 = -2.919d10 dd22 = -1.185d10 dd23 = -7.270d9 dd24 = -4.222d9 dd25 = -1.560d9 else if (temp .ge. 1.0d8 .and. temp .lt. 1.0d9) then tau = logtemp - 8d0 cc = 1.5654d0 c00 = 9.889d10 c01 = -4.524d8 c02 = -6.088d6 c03 = 4.269d7 c04 = 5.172d7 c05 = 4.910d7 c06 = 4.388d7 c10 = 1.813d11 c11 = -7.556d9 c12 = -3.304d9 c13 = -1.031d9 c14 = -1.764d9 c15 = -1.851d9 c16 = -1.928d9 c20 = 9.750d10 c21 = 3.484d10 c22 = 5.199d9 c23 = -1.695d9 c24 = -2.865d9 c25 = -3.395d9 c26 = -3.418d9 dd01 = -1.135d8 dd02 = 1.256d8 dd03 = 5.149d7 dd04 = 3.436d7 dd05 = 1.005d7 dd11 = 1.652d9 dd12 = -3.119d9 dd13 = -1.839d9 dd14 = -1.458d9 dd15 = -8.956d8 dd21 = -1.549d10 dd22 = -9.338d9 dd23 = -5.899d9 dd24 = -3.035d9 dd25 = -1.598d9 else if (temp .ge. 1.0d9) then tau = logtemp - 9d0 cc = 1.5654d0 c00 = 9.581d10 c01 = 4.107d8 c02 = 2.305d8 c03 = 2.236d8 c04 = 1.580d8 c05 = 2.165d8 c06 = 1.721d8 c10 = 1.459d12 c11 = 1.314d11 c12 = -1.169d11 c13 = -1.765d11 c14 = -1.867d11 c15 = -1.983d11 c16 = -1.896d11 c20 = 2.424d11 c21 = -3.669d9 c22 = -8.691d9 c23 = -7.967d9 c24 = -7.932d9 c25 = -7.987d9 c26 = -8.333d9 dd01 = 4.724d8 dd02 = 2.976d8 dd03 = 2.242d8 dd04 = 7.937d7 dd05 = 4.859d7 dd11 = -7.094d11 dd12 = -3.697d11 dd13 = -2.189d11 dd14 = -1.273d11 dd15 = -5.705d10 dd21 = -2.254d10 dd22 = -1.551d10 dd23 = -7.793d9 dd24 = -4.489d9 dd25 = -2.185d9 end if taudt = iln10*tempi !..equation 3.7, compute the expensive trig functions only one time cos1 = cos_cr(fac1*tau) cos2 = cos_cr(fac1*2.0d0*tau) cos3 = cos_cr(fac1*3.0d0*tau) cos4 = cos_cr(fac1*4.0d0*tau) cos5 = cos_cr(fac1*5.0d0*tau) last = cos_cr(fac2*tau) sin1 = sin_cr(fac1*tau) sin2 = sin_cr(fac1*2.0d0*tau) sin3 = sin_cr(fac1*3.0d0*tau) sin4 = sin_cr(fac1*4.0d0*tau) sin5 = sin_cr(fac1*5.0d0*tau) xast = sin_cr(fac2*tau) a0 = 0.5d0*c00 & + c01*cos1 + dd01*sin1 + c02*cos2 + dd02*sin2 & + c03*cos3 + dd03*sin3 + c04*cos4 + dd04*sin4 & + c05*cos5 + dd05*sin5 + 0.5d0*c06*last f0 = taudt*fac1*(-c01*sin1 + dd01*cos1 - c02*sin2*2.0d0 & + dd02*cos2*2.0d0 - c03*sin3*3.0d0 + dd03*cos3*3.0d0 & - c04*sin4*4.0d0 + dd04*cos4*4.0d0 & - c05*sin5*5.0d0 + dd05*cos5*5.0d0) & - 0.5d0*c06*xast*fac2*taudt a1 = 0.5d0*c10 & + c11*cos1 + dd11*sin1 + c12*cos2 + dd12*sin2 & + c13*cos3 + dd13*sin3 + c14*cos4 + dd14*sin4 & + c15*cos5 + dd15*sin5 + 0.5d0*c16*last f1 = taudt*fac1*(-c11*sin1 + dd11*cos1 - c12*sin2*2.0d0 & + dd12*cos2*2.0d0 - c13*sin3*3.0d0 + dd13*cos3*3.0d0 & - c14*sin4*4.0d0 + dd14*cos4*4.0d0 - c15*sin5*5.0d0 & + dd15*cos5*5.0d0) - 0.5d0*c16*xast*fac2*taudt a2 = 0.5d0*c20 & + c21*cos1 + dd21*sin1 + c22*cos2 + dd22*sin2 & + c23*cos3 + dd23*sin3 + c24*cos4 + dd24*sin4 & + c25*cos5 + dd25*sin5 + 0.5d0*c26*last f2 = taudt*fac1*(-c21*sin1 + dd21*cos1 - c22*sin2*2.0d0 & + dd22*cos2*2.0d0 - c23*sin3*3.0d0 + dd23*cos3*3.0d0 & - c24*sin4*4.0d0 + dd24*cos4*4.0d0 - c25*sin5*5.0d0 & + dd25*cos5*5.0d0) - 0.5d0*c26*xast*fac2*taudt !..equation 3.4 dum = a0 + a1*zeta + a2*zeta2 dumdt = f0 + f1*zeta + a1*zetadt + f2*zeta2 + 2.0d0*a2*zeta*zetadt dumdd = a1*zetadd + 2.0d0*a2*zeta*zetadd dumda = a1*zetada + 2.0d0*a2*zeta*zetada dumdz = a1*zetadz + 2.0d0*a2*zeta*zetadz z = exp_cr(-cc*zeta) xnum = dum*z xnumdt = dumdt*z - dum*z*cc*zetadt xnumdd = dumdd*z - dum*z*cc*zetadd xnumda = dumda*z - dum*z*cc*zetada xnumdz = dumdz*z - dum*z*cc*zetadz xden = zeta3 + 6.290d-3*xlm1 + 7.483d-3*xlm2 + 3.061d-4*xlm3 dum = 3.0d0*zeta2 xdendt = dum*zetadt - xldt*(6.290d-3*xlm2 & + 2.0d0*7.483d-3*xlm3 + 3.0d0*3.061d-4*xlm4) xdendd = dum*zetadd xdenda = dum*zetada xdendz = dum*zetadz dum = 1.0d0/xden fphot = xnum*dum fphotdt = (xnumdt - fphot*xdendt)*dum fphotdd = (xnumdd - fphot*xdendd)*dum fphotda = (xnumda - fphot*xdenda)*dum fphotdz = (xnumdz - fphot*xdendz)*dum !..equation 3.3 a0 = 1.0d0 + 2.045d0 * xl xnum = 0.666d0*pow_cr(a0,-2.066d0) xnumdt = -2.066d0*xnum/a0 * 2.045d0*xldt dum = 1.875d8*xl + 1.653d8*xl2 + 8.449d8*xl3 - 1.604d8*xl4 dumdt = xldt*(1.875d8 + 2.0d0*1.653d8*xl + 3.0d0*8.449d8*xl2 & - 4.0d0*1.604d8*xl3) z = 1.0d0/dum xden = 1.0d0 + rm*z xdendt = -rm*z*z*dumdt xdendd = rmdd*z xdenda = rmda*z xdendz = rmdz*z z = 1.0d0/xden qphot = xnum*z qphotdt = (xnumdt - qphot*xdendt)*z dum = -qphot*z qphotdd = dum*xdendd qphotda = dum*xdenda qphotdz = dum*xdendz !..equation 3.2 sphot = xl5 * fphot sphotdt = 5.0d0*xl4*xldt*fphot + xl5*fphotdt sphotdd = xl5*fphotdd sphotda = xl5*fphotda sphotdz = xl5*fphotdz a1 = sphot sphot = rm*a1 sphotdt = rm*sphotdt sphotdd = rm*sphotdd + rmdd*a1 sphotda = rm*sphotda + rmda*a1 sphotdz = rm*sphotdz + rmdz*a1 a1 = tfac4*(1.0d0 - tfac3 * qphot) a2 = -tfac4*tfac3 a3 = sphot sphot = a1*a3 sphotdt = a1*sphotdt + a2*qphotdt*a3 sphotdd = a1*sphotdd + a2*qphotdd*a3 sphotda = a1*sphotda + a2*qphotda*a3 sphotdz = a1*sphotdz + a2*qphotdz*a3 if (.false.) then write(*,*) 'logT', logtemp write(*,*) 'logRho', logden write(*,*) 'sphot', sphot write(*,*) end if if (sphot .le. 0.0) then sphot = 0.0d0 sphotdt = 0.0d0 sphotdd = 0.0d0 sphotda = 0.0d0 sphotdz = 0.0d0 end if end subroutine phot_neu subroutine brem_neu_weak_degen real(dp) :: p !..equation 5.3 dum = 7.05d6 * t832 + 5.12d4 * t83 dumdt = (1.5d0*7.05d6*t812 + 3.0d0*5.12d4*t82)*1.0d-8 z = 1.0d0/dum eta = rm*z etadt = -rm*z*z*dumdt etadd = rmdd*z etada = rmda*z etadz = rmdz*z etam1 = 1.0d0/eta etam2 = etam1 * etam1 etam3 = etam2 * etam1 !..equation 5.2 a0 = 23.5d0 + 6.83d4*t8m2 + 7.81d8*t8m5 f0 = (-2.0d0*6.83d4*t8m3 - 5.0d0*7.81d8*t8m6)*1.0d-8 xnum = 1.0d0/a0 dum = 1.0d0 + 1.47d0*etam1 + 3.29d-2*etam2 z = -1.47d0*etam2 - 2.0d0*3.29d-2*etam3 dumdt = z*etadt dumdd = z*etadd dumda = z*etada dumdz = z*etadz c00 = 1.26d0 * (1.0d0+etam1) z = -1.26d0*etam2 c01 = z*etadt c02 = z*etadd c03 = z*etada c04 = z*etadz z = 1.0d0/dum xden = c00*z xdendt = (c01 - xden*dumdt)*z xdendd = (c02 - xden*dumdd)*z xdenda = (c03 - xden*dumda)*z xdendz = (c04 - xden*dumdz)*z fbrem = xnum + xden fbremdt = -xnum*xnum*f0 + xdendt fbremdd = xdendd fbremda = xdenda fbremdz = xdendz !..equation 5.9 a0 = 230.0d0 + 6.7d5*t8m2 + 7.66d9*t8m5 f0 = (-2.0d0*6.7d5*t8m3 - 5.0d0*7.66d9*t8m6)*1.0d-8 z = 1.0d0 + rm*1.0d-9 dum = a0*z dumdt = f0*z z = a0*1.0d-9 dumdd = z*rmdd dumda = z*rmda dumdz = z*rmdz xnum = 1.0d0/dum z = -xnum*xnum xnumdt = z*dumdt xnumdd = z*dumdd xnumda = z*dumda xnumdz = z*dumdz p = pow_cr(t8,3.85d0) c00 = 7.75d5*t832 + 247.0d0*p dd00 = (1.5d0*7.75d5*t812 + 3.85d0*247.0d0*p/T8)*1.0d-8 p = pow_cr(t8,1.4d0) c01 = 4.07d0 + 0.0240d0 * p dd01 = 1.4d0*0.0240d0*(p/T8)*1.0d-8 p = pow_cr(t8,-0.110d0) c02 = 4.59d-5 * p dd02 = -0.11d0*4.59d-5*(p/T8)*1.0d-8 z = pow_cr(den,0.656d0) dum = c00*rmi + c01 + c02*z dumdt = dd00*rmi + dd01 + dd02*z z = -c00*rmi*rmi dumdd = z*rmdd + 0.656d0*c02*pow_cr(den,-0.454d0) dumda = z*rmda dumdz = z*rmdz xden = 1.0d0/dum z = -xden*xden xdendt = z*dumdt xdendd = z*dumdd xdenda = z*dumda xdendz = z*dumdz gbrem = xnum + xden gbremdt = xnumdt + xdendt gbremdd = xnumdd + xdendd gbremda = xnumda + xdenda gbremdz = xnumdz + xdendz !..equation 5.1 dum = 0.5738d0*zbar*ye*t86*den dumdt = 0.5738d0*zbar*ye*6.0d0*t85*den*1.0d-8 dumdd = 0.5738d0*zbar*ye*t86 dumda = -dum*abari dumdz = 0.5738d0*2.0d0*ye*t86*den z = tfac4*fbrem - tfac5*gbrem sbrem = dum * z sbremdt = dumdt*z + dum*(tfac4*fbremdt - tfac5*gbremdt) sbremdd = dumdd*z + dum*(tfac4*fbremdd - tfac5*gbremdd) sbremda = dumda*z + dum*(tfac4*fbremda - tfac5*gbremda) sbremdz = dumdz*z + dum*(tfac4*fbremdz - tfac5*gbremdz) end subroutine brem_neu_weak_degen subroutine brem_neu_liquid_metal !..liquid metal with c12 parameters (not too different for other elements) !..equation 5.18 and 5.16 u = fac3 * (logden - 3.0d0) a0 = iln10*fac3 * deni !..compute the expensive trig functions of equation 5.21 only once cos1 = cos_cr(u) cos2 = cos_cr(2.0d0*u) cos3 = cos_cr(3.0d0*u) cos4 = cos_cr(4.0d0*u) cos5 = cos_cr(5.0d0*u) sin1 = sin_cr(u) sin2 = sin_cr(2.0d0*u) sin3 = sin_cr(3.0d0*u) sin4 = sin_cr(4.0d0*u) sin5 = sin_cr(5.0d0*u) !..equation 5.21 fb = 0.5d0 * 0.17946d0 + 0.00945d0*u + 0.34529d0 & - 0.05821d0*cos1 - 0.04969d0*sin1 & - 0.01089d0*cos2 - 0.01584d0*sin2 & - 0.01147d0*cos3 - 0.00504d0*sin3 & - 0.00656d0*cos4 - 0.00281d0*sin4 & - 0.00519d0*cos5 c00 = a0*(0.00945d0 & + 0.05821d0*sin1 - 0.04969d0*cos1 & + 0.01089d0*sin2*2.0d0 - 0.01584d0*cos2*2.0d0 & + 0.01147d0*sin3*3.0d0 - 0.00504d0*cos3*3.0d0 & + 0.00656d0*sin4*4.0d0 - 0.00281d0*cos4*4.0d0 & + 0.00519d0*sin5*5.0d0) !..equation 5.22 ft = 0.5d0 * 0.06781d0 - 0.02342d0*u + 0.24819d0 & - 0.00944d0*cos1 - 0.02213d0*sin1 & - 0.01289d0*cos2 - 0.01136d0*sin2 & - 0.00589d0*cos3 - 0.00467d0*sin3 & - 0.00404d0*cos4 - 0.00131d0*sin4 & - 0.00330d0*cos5 c01 = a0*(-0.02342d0 & + 0.00944d0*sin1 - 0.02213d0*cos1 & + 0.01289d0*sin2*2.0d0 - 0.01136d0*cos2*2.0d0 & + 0.00589d0*sin3*3.0d0 - 0.00467d0*cos3*3.0d0 & + 0.00404d0*sin4*4.0d0 - 0.00131d0*cos4*4.0d0 & + 0.00330d0*sin5*5.0d0) !..equation 5.23 gb = 0.5d0 * 0.00766d0 - 0.01259d0*u + 0.07917d0 & - 0.00710d0*cos1 + 0.02300d0*sin1 & - 0.00028d0*cos2 - 0.01078d0*sin2 & + 0.00232d0*cos3 + 0.00118d0*sin3 & + 0.00044d0*cos4 - 0.00089d0*sin4 & + 0.00158d0*cos5 c02 = a0*(-0.01259d0 & + 0.00710d0*sin1 + 0.02300d0*cos1 & + 0.00028d0*sin2*2.0d0 - 0.01078d0*cos2*2.0d0 & - 0.00232d0*sin3*3.0d0 + 0.00118d0*cos3*3.0d0 & - 0.00044d0*sin4*4.0d0 - 0.00089d0*cos4*4.0d0 & - 0.00158d0*sin5*5.0d0) !..equation 5.24 gt = -0.5d0 * 0.00769d0 - 0.00829d0*u + 0.05211d0 & + 0.00356d0*cos1 + 0.01052d0*sin1 & - 0.00184d0*cos2 - 0.00354d0*sin2 & + 0.00146d0*cos3 - 0.00014d0*sin3 & + 0.00031d0*cos4 - 0.00018d0*sin4 & + 0.00069d0*cos5 c03 = a0*(-0.00829d0 & - 0.00356d0*sin1 + 0.01052d0*cos1 & + 0.00184d0*sin2*2.0d0 - 0.00354d0*cos2*2.0d0 & - 0.00146d0*sin3*3.0d0 - 0.00014d0*cos3*3.0d0 & - 0.00031d0*sin4*4.0d0 - 0.00018d0*cos4*4.0d0 & - 0.00069d0*sin5*5.0d0) dum = 2.275d-1 * zbar * zbar*t8m1 * pow_cr(den6*abari, one_third) dumdt = -dum*tempi dumdd = one_third*dum * deni dumda = -one_third*dum*abari dumdz = 2.0d0*dum*zbari gm1 = 1.0d0/dum gm2 = gm1*gm1 gm13 = pow_cr(gm1,one_third) gm23 = gm13 * gm13 gm43 = gm13*gm1 gm53 = gm23*gm1 !..equation 5.25 and 5.26 v = -0.05483d0 - 0.01946d0*gm13 + 1.86310d0*gm23 - 0.78873d0*gm1 a0 = one_third*0.01946d0*gm43 - two_thirds*1.86310d0*gm53 + 0.78873d0*gm2 w = -0.06711d0 + 0.06859d0*gm13 + 1.74360d0*gm23 - 0.74498d0*gm1 a1 = -one_third*0.06859d0*gm43 - two_thirds*1.74360d0*gm53 + 0.74498d0*gm2 !..equation 5.19 and 5.20 fliq = v*fb + (1.0d0 - v)*ft fliqdt = a0*dumdt*(fb - ft) fliqdd = a0*dumdd*(fb - ft) + v*c00 + (1.0d0 - v)*c01 fliqda = a0*dumda*(fb - ft) fliqdz = a0*dumdz*(fb - ft) gliq = w*gb + (1.0d0 - w)*gt gliqdt = a1*dumdt*(gb - gt) gliqdd = a1*dumdd*(gb - gt) + w*c02 + (1.0d0 - w)*c03 gliqda = a1*dumda*(gb - gt) gliqdz = a1*dumdz*(gb - gt) !..equation 5.17 dum = 0.5738d0*zbar*ye*t86*den dumdt = 0.5738d0*zbar*ye*6.0d0*t85*den*1.0d-8 dumdd = 0.5738d0*zbar*ye*t86 dumda = -dum*abari dumdz = 0.5738d0*2.0d0*ye*t86*den z = tfac4*fliq - tfac5*gliq sbrem = dum * z sbremdt = dumdt*z + dum*(tfac4*fliqdt - tfac5*gliqdt) sbremdd = dumdd*z + dum*(tfac4*fliqdd - tfac5*gliqdd) sbremda = dumda*z + dum*(tfac4*fliqda - tfac5*gliqda) sbremdz = dumdz*z + dum*(tfac4*fliqdz - tfac5*gliqdz) end subroutine brem_neu_liquid_metal subroutine brem_neu !..bremsstrahlung neutrino section !..for reactions like e- + (z,a) => e- + (z,a) + nu + nubar !.. n + n => n + n + nu + nubar !.. n + p => n + p + nu + nubar real(dp) :: alfa, beta, sb, sbdt, sbdd, sbda, sbdz real(dp), parameter :: tflo = 0.05d0, tfhi = 0.3d0 !..equation 4.3 den6 = den * 1.0d-6 t8 = temp * 1.0d-8 t812 = sqrt(t8) t832 = t8 * t812 t82 = t8*t8 t83 = t82*t8 t85 = t82*t83 t86 = t85*t8 t8m1 = 1.0d0/t8 t8m2 = t8m1*t8m1 t8m3 = t8m2*t8m1 t8m5 = t8m3*t8m2 t8m6 = t8m5*t8m1 tfermi = 5.9302d9*(sqrt(1.0d0+1.018d0* pow_cr(den6*ye,two_thirds) )-1.0d0) if (temp .ge. tfhi * tfermi) then call brem_neu_weak_degen() else if (temp .le. tflo * tfermi) then call brem_neu_liquid_metal() else ! blend call brem_neu_weak_degen() sb = sbrem sbdt = sbremdt sbdd = sbremdd sbda = sbremda sbdz = sbremdz call brem_neu_liquid_metal() alfa = (temp / tfermi - tflo) / (tfhi - tflo) alfa = 0.5d0 * (1d0 - cospi_cr(alfa)) beta = 1d0 - alfa sbrem = alfa * sb + beta * sbrem sbremdt = alfa * sbdt + beta * sbremdt sbremdd = alfa * sbdd + beta * sbremdd sbremda = alfa * sbda + beta * sbremda sbremdz = alfa * sbdz + beta * sbremdz end if end subroutine brem_neu subroutine reco_neu !..recombination neutrino section !..for reactions like e- (continuum) => e- (bound) + nu_e + nubar_e !..equation 6.11 solved for nu xnum = 1.10520d8 * den * ye /(temp*sqrt(temp)) xnumdt = -1.50d0*xnum*tempi xnumdd = xnum * deni xnumda = -xnum*abari xnumdz = xnum*zbari !..the chemical potential nu = ifermi12(xnum) !..a0 is d(nu)/d(xnum) a0 = 1.0d0/(0.5d0*zfermim12(nu)) nudt = a0*xnumdt nudd = a0*xnumdd nuda = a0*xnumda nudz = a0*xnumdz nu2 = nu * nu nu3 = nu2 * nu !..table 12 if (nu .ge. -20.0 .and. nu .lt. 0.0) then a1 = 1.51d-2 a2 = 2.42d-1 a3 = 1.21d0 b = 3.71d-2 c = 9.06e-1 d = 9.28d-1 f1 = 0.0d0 f2 = 0.0d0 f3 = 0.0d0 else if (nu .ge. 0.0 .and. nu .le. 10.0) then a1 = 1.23d-2 a2 = 2.66d-1 a3 = 1.30d0 b = 1.17d-1 c = 8.97e-1 d = 1.77d-1 f1 = -1.20d-2 f2 = 2.29d-2 f3 = -1.04d-3 end if !..equation 6.7, 6.13 and 6.14 if (nu .ge. -20.0 .and. nu .le. 10.0) then zeta = 1.579d5*zbar*zbar*tempi zetadt = -zeta*tempi zetadd = 0.0d0 zetada = 0.0d0 zetadz = 2.0d0*zeta*zbari c00 = 1.0d0/(1.0d0 + f1*nu + f2*nu2 + f3*nu3) c01 = f1 + f2*2.0d0*nu + f3*3.0d0*nu2 dum = zeta*c00 dumdt = zetadt*c00 + zeta*c01*nudt dumdd = zeta*c01*nudd dumda = zeta*c01*nuda dumdz = zetadz*c00 + zeta*c01*nudz z = 1.0d0/dum dd00 = pow_cr(dum,-2.25d0) dd01 = pow_cr(dum,-4.55d0) c00 = a1*z + a2*dd00 + a3*dd01 c01 = -(a1*z + 2.25*a2*dd00 + 4.55*a3*dd01)*z z = exp_cr(c*nu) dd00 = b*z*(1.0d0 + d*dum) gum = 1.0d0 + dd00 gumdt = dd00*c*nudt + b*z*d*dumdt gumdd = dd00*c*nudd + b*z*d*dumdd gumda = dd00*c*nuda + b*z*d*dumda gumdz = dd00*c*nudz + b*z*d*dumdz z = exp_cr(nu) a1 = 1.0d0/gum bigj = c00 * z * a1 bigjdt = c01*dumdt*z*a1 + c00*z*nudt*a1 - c00*z*a1*a1 * gumdt bigjdd = c01*dumdd*z*a1 + c00*z*nudd*a1 - c00*z*a1*a1 * gumdd bigjda = c01*dumda*z*a1 + c00*z*nuda*a1 - c00*z*a1*a1 * gumda bigjdz = c01*dumdz*z*a1 + c00*z*nudz*a1 - c00*z*a1*a1 * gumdz !..equation 6.5 z = exp_cr(zeta + nu) dum = 1.0d0 + z a1 = 1.0d0/dum a2 = 1.0d0/bigj sreco = tfac6 * 2.649d-18 * ye * powi_cr(zbar,13) * den * bigj*a1 srecodt = sreco*(bigjdt*a2 - z*(zetadt + nudt)*a1) srecodd = sreco*(1.0d0 * deni + bigjdd*a2 - z*(zetadd + nudd)*a1) srecoda = sreco*(-1.0d0*abari + bigjda*a2 - z*(zetada+nuda)*a1) srecodz = sreco*(14.0d0*zbari + bigjdz*a2 - z*(zetadz+nudz)*a1) end if end subroutine reco_neu subroutine plas_neu !..plasma neutrino section !..for collective reactions like gamma_plasmon => nu_e + nubar_e !..equation 4.6 a1 = 1.019d-6*rm a2 = pow_cr(a1,two_thirds) a3 = two_thirds*a2/a1 b1 = sqrt(1.0d0 + a2) b2 = 1.0d0/b1 c00 = 1.0d0/(temp*temp*b1) gl2 = 1.1095d11 * rm * c00 gl2dt = -2.0d0*gl2*tempi d = rm*c00*b2*0.5d0*b2*a3*1.019d-6 gl2dd = 1.1095d11 * (rmdd*c00 - d*rmdd) gl2da = 1.1095d11 * (rmda*c00 - d*rmda) gl2dz = 1.1095d11 * (rmdz*c00 - d*rmdz) gl = sqrt(gl2) gl12 = sqrt(gl) gl32 = gl * gl12 gl72 = gl2 * gl32 gl6 = gl2 * gl2 * gl2 !..equation 4.7 ft = 2.4d0 + 0.6d0*gl12 + 0.51d0*gl + 1.25d0*gl32 gum = 1.0d0/gl2 a1 =(0.25d0*0.6d0*gl12 +0.5d0*0.51d0*gl +0.75d0*1.25d0*gl32)*gum ftdt = a1*gl2dt ftdd = a1*gl2dd ftda = a1*gl2da ftdz = a1*gl2dz !..equation 4.8 a1 = 8.6d0*gl2 + 1.35d0*gl72 a2 = 8.6d0 + 1.75d0*1.35d0*gl72*gum b1 = 225.0d0 - 17.0d0*gl + gl2 b2 = -0.5d0*17.0d0*gl*gum + 1.0d0 c = 1.0d0/b1 fl = a1*c d = (a2 - fl*b2)*c fldt = d*gl2dt fldd = d*gl2dd flda = d*gl2da fldz = d*gl2dz !..equation 4.9 and 4.10 cc = log10_cr(2.0d0*rm) xlnt = logtemp xnum = sixth * (17.5d0 + cc - 3.0d0*xlnt) xnumdt = -iln10*0.5d0*tempi a2 = iln10*sixth*rmi xnumdd = a2*rmdd xnumda = a2*rmda xnumdz = a2*rmdz xden = sixth * (-24.5d0 + cc + 3.0d0*xlnt) xdendt = iln10*0.5d0*tempi xdendd = a2*rmdd xdenda = a2*rmda xdendz = a2*rmdz !..equation 4.11 if (abs(xnum) .gt. 0.7d0 .or. xden .lt. 0.0d0) then fxy = 1.0d0 fxydt = 0.0d0 fxydd = 0.0d0 fxydz = 0.0d0 fxyda = 0.0d0 else a1 = 0.39d0 - 1.25d0*xnum - 0.35d0*sin_cr(4.5d0*xnum) a2 = -1.25d0 - 4.5d0*0.35d0*cos_cr(4.5d0*xnum) b1 = 0.3d0 * exp_cr(-1.0d0*(4.5d0*xnum + 0.9d0)**2) b2 = -b1*2.0d0*(4.5d0*xnum + 0.9d0)*4.5d0 c = min(0.0d0, xden - 1.6d0 + 1.25d0*xnum) if (c .eq. 0.0) then dumdt = 0.0d0 dumdd = 0.0d0 dumda = 0.0d0 dumdz = 0.0d0 else dumdt = xdendt + 1.25d0*xnumdt dumdd = xdendd + 1.25d0*xnumdd dumda = xdenda + 1.25d0*xnumda dumdz = xdendz + 1.25d0*xnumdz end if d = 0.57d0 - 0.25d0*xnum a3 = c/d c00 = exp_cr(-1.0d0*a3**2) f1 = -c00*2.0d0*a3/d c01 = f1*(dumdt + a3*0.25d0*xnumdt) c02 = f1*(dumdd + a3*0.25d0*xnumdd) c03 = f1*(dumda + a3*0.25d0*xnumda) c04 = f1*(dumdz + a3*0.25d0*xnumdz) fxy = 1.05d0 + (a1 - b1)*c00 fxydt = (a2*xnumdt - b2*xnumdt)*c00 + (a1-b1)*c01 fxydd = (a2*xnumdd - b2*xnumdd)*c00 + (a1-b1)*c02 fxyda = (a2*xnumda - b2*xnumda)*c00 + (a1-b1)*c03 fxydz = (a2*xnumdz - b2*xnumdz)*c00 + (a1-b1)*c04 end if !..equation 4.1 and 4.5 splas = (ft + fl) * fxy splasdt = (ftdt + fldt)*fxy + (ft+fl)*fxydt splasdd = (ftdd + fldd)*fxy + (ft+fl)*fxydd splasda = (ftda + flda)*fxy + (ft+fl)*fxyda splasdz = (ftdz + fldz)*fxy + (ft+fl)*fxydz a2 = exp_cr(-gl) a3 = -0.5d0*a2*gl*gum a1 = splas splas = a2*a1 splasdt = a2*splasdt + a3*gl2dt*a1 splasdd = a2*splasdd + a3*gl2dd*a1 splasda = a2*splasda + a3*gl2da*a1 splasdz = a2*splasdz + a3*gl2dz*a1 a2 = gl6 a3 = 3.0d0*gl6*gum a1 = splas splas = a2*a1 splasdt = a2*splasdt + a3*gl2dt*a1 splasdd = a2*splasdd + a3*gl2dd*a1 splasda = a2*splasda + a3*gl2da*a1 splasdz = a2*splasdz + a3*gl2dz*a1 a2 = 0.93153d0 * 3.0d21 * xl9 a3 = 0.93153d0 * 3.0d21 * 9.0d0*xl8*xldt a1 = splas splas = a2*a1 splasdt = a2*splasdt + a3*a1 splasdd = a2*splasdd splasda = a2*splasda splasdz = a2*splasdz end subroutine plas_neu subroutine pair_neu !..pair neutrino section !..for reactions like e+ + e- => nu_e + nubar_e !..equation 2.8 gl = 1.0d0 - 13.04d0*xl2 +133.5d0*xl4 +1534.0d0*xl6 +918.6d0*xl8 gldt = xldt*(-26.08d0*xl +534.0d0*xl3 +9204.0d0*xl5 +7348.8d0*xl7) !..equation 2.7 a1 = 6.002d19 + 2.084d20*zeta + 1.872d21*zeta2 a2 = 2.084d20 + 2.0d0*1.872d21*zeta if (t9 .lt. 10.0) then b1 = exp_cr(-5.5924d0*zeta) b2 = -b1*5.5924d0 else b1 = exp_cr(-4.9924d0*zeta) b2 = -b1*4.9924d0 end if xnum = a1 * b1 c = a2*b1 + a1*b2 xnumdt = c*zetadt xnumdd = c*zetadd xnumda = c*zetada xnumdz = c*zetadz if (t9 .lt. 10.0) then a1 = 9.383d-1*xlm1 - 4.141d-1*xlm2 + 5.829d-2*xlm3 a2 = -9.383d-1*xlm2 + 2.0d0*4.141d-1*xlm3 - 3.0d0*5.829d-2*xlm4 else a1 = 1.2383d0*xlm1 - 8.141d-1*xlm2 a2 = -1.2383d0*xlm2 + 2.0d0*8.141d-1*xlm3 end if b1 = 3.0d0*zeta2 xden = zeta3 + a1 xdendt = b1*zetadt + a2*xldt xdendd = b1*zetadd xdenda = b1*zetada xdendz = b1*zetadz a1 = 1.0d0/xden fpair = xnum*a1 fpairdt = (xnumdt - fpair*xdendt)*a1 fpairdd = (xnumdd - fpair*xdendd)*a1 fpairda = (xnumda - fpair*xdenda)*a1 fpairdz = (xnumdz - fpair*xdendz)*a1 !..equation 2.6 a1 = 10.7480d0*xl2 + 0.3967d0*xlp5 + 1.005d0 a2 = xldt*(2.0d0*10.7480d0*xl + 0.5d0*0.3967d0*xlmp5) xnum = 1.0d0/a1 xnumdt = -xnum*xnum*a2 a1 = 7.692d7*xl3 + 9.715d6*xlp5 a2 = xldt*(3.0d0*7.692d7*xl2 + 0.5d0*9.715d6*xlmp5) c = 1.0d0/a1 b1 = 1.0d0 + rm*c xden = pow_cr(b1,-0.3d0) d = -0.3d0*xden/b1 xdendt = -d*rm*c*c*a2 xdendd = d*rmdd*c xdenda = d*rmda*c xdendz = d*rmdz*c qpair = xnum*xden qpairdt = xnumdt*xden + xnum*xdendt qpairdd = xnum*xdendd qpairda = xnum*xdenda qpairdz = xnum*xdendz include 'formats.dek' !..equation 2.5 a1 = exp_cr(-2.0d0*xlm1) a2 = a1*2.0d0*xlm2*xldt spair = a1*fpair spairdt = a2*fpair + a1*fpairdt spairdd = a1*fpairdd spairda = a1*fpairda spairdz = a1*fpairdz a1 = spair spair = gl*a1 spairdt = gl*spairdt + gldt*a1 spairdd = gl*spairdd spairda = gl*spairda spairdz = gl*spairdz a1 = tfac4*(1.0d0 + tfac3 * qpair) a2 = tfac4*tfac3 a3 = spair spair = a1*a3 spairdt = a1*spairdt + a2*qpairdt*a3 spairdd = a1*spairdd + a2*qpairdd*a3 spairda = a1*spairda + a2*qpairda*a3 spairdz = a1*spairdz + a2*qpairdz*a3 end subroutine pair_neu end subroutine my_neutrinos end module run_star_extras