this is for holding javascript data
Matteo Cantiello edited MLT.f.tex
over 2 years ago
Commit id: 61d97251bbff5320b5f182215eb8399eaa8325e8
deletions | additions
diff --git a/MLT.f.tex b/MLT.f.tex
index 7020895..e69de29 100644
--- a/MLT.f.tex
+++ b/MLT.f.tex
...
\section{Mlt.f}
Below the routine used in MESA to calculate the quantities described above.
It goes beyond, as it also includes thermohaline mixing etc...
\begin{lstlisting}
! ***********************************************************************
!
! Copyright (C) 2012 Bill Paxton, Mike Montgomery, Pascale Garaud
!
! MESA is free software; you can use it and/or modify
! it under the combined terms and restrictions of the MESA MANIFESTO
! and 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.
!
! You should have received a copy of the MESA MANIFESTO along with
! this software; if not, it is available at the mesa website:
! http://mesa.sourceforge.net/
!
! 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 mlt
use mlt_def
use const_def
use crlibm_lib
use utils_lib, only: is_bad_num
implicit none
#ifdef offload
!dir$ options /offload_attribute_target=mic
#endif
! largely based on standard mixing length theory
! as, for example, described in Cox & Giuli, Chapter 14. (C&G)
! with extensive additions for semiconvection and thermohaline and more.....
integer, parameter :: nvbs = num_mlt_partials
contains
subroutine do_mlt_eval( &
have_L_need_gradT, cgrav, m, r, T, rho, L, dg, P, &
chiRho, chiT, Cp, Cv, csound, X, opacity, grada, dPrad_dm, &
gradr_factor, gradL_composition_term, &
alpha_semiconvection, semiconvection_option, &
thermohaline_coeff, thermohaline_option, &
dominant_iso_for_thermohaline, &
mixing_length_alpha, alt_scale_height, remove_small_D_limit, &
MLT_option, Henyey_y_param, Henyey_nu_param, &
prev_conv_vel, max_conv_vel, dt, tau, MLT_dbg, &
mixing_type, mlt_basics, mlt_partials1, ierr)
logical, intent(in) :: have_L_need_gradT
real(dp), intent(in) :: cgrav, m, r, T, Rho, L, dg, P
real(dp), intent(in) :: chiRho, chiT, Cp, Cv, csound, X, &
opacity, grada, dPrad_dm, gradr_factor
real(dp), intent(in) :: gradL_composition_term
real(dp), intent(in) :: alpha_semiconvection, thermohaline_coeff
real(dp), intent(in) :: mixing_length_alpha
logical, intent(in) :: alt_scale_height
character (len=*), intent(in) :: thermohaline_option, MLT_option, semiconvection_option
integer, intent(in) :: dominant_iso_for_thermohaline
real(dp), intent(in) :: Henyey_y_param, Henyey_nu_param, &
prev_conv_vel, max_conv_vel, dt, tau, remove_small_D_limit
logical, intent(in) :: MLT_dbg
integer, intent(out) :: mixing_type
real(dp), intent(out) :: mlt_basics(:) ! (num_mlt_results)
real(dp), intent(out), pointer :: mlt_partials1(:) ! =(num_mlt_partials, num_mlt_results)
! e.g., mlt_partials(mlt_dlnT,mlt_gradT) has partial wrt lnT of gradT
integer, intent(out) :: ierr
real(dp), pointer :: mlt_partials(:,:)
mlt_partials(1:num_mlt_partials,1:num_mlt_results) => &
mlt_partials1(1:num_mlt_partials*num_mlt_results)
if (MLT_dbg) write(*,*) 'start do_mlt_eval'
ierr = 0
mlt_basics = 0
mlt_partials = 0
call Get_results( &
have_L_need_gradT, cgrav, m, r, T, Rho, L, dg, P, &
chiRho, chiT, Cp, Cv, csound, X, opacity, grada, dPrad_dm, &
gradr_factor, gradL_composition_term, &
alpha_semiconvection, semiconvection_option, &
thermohaline_coeff, thermohaline_option, &
dominant_iso_for_thermohaline, &
mixing_length_alpha, alt_scale_height, remove_small_D_limit, &
MLT_option, Henyey_y_param, Henyey_nu_param, &
prev_conv_vel, max_conv_vel, dt, tau, MLT_dbg, mixing_type, &
mlt_basics(mlt_gradT), mlt_partials(:,mlt_gradT), &
mlt_basics(mlt_L_conv), mlt_partials(:,mlt_L_conv), &
mlt_basics(mlt_L_total), mlt_partials(:,mlt_L_total), &
mlt_basics(mlt_gradr), mlt_partials(:,mlt_gradr), &
mlt_basics(mlt_gradL), mlt_partials(:,mlt_gradL), &
mlt_basics(mlt_scale_height), mlt_partials(:,mlt_scale_height), &
mlt_basics(mlt_Lambda), mlt_partials(:,mlt_Lambda), &
mlt_basics(mlt_convection_velocity), mlt_partials(:,mlt_convection_velocity), &
mlt_basics(mlt_D), mlt_partials(:,mlt_D), &
mlt_basics(mlt_gamma), mlt_partials(:,mlt_gamma), &
mlt_basics(mlt_conv_dP_term), mlt_partials(:,mlt_conv_dP_term), &
ierr)
if (mlt_basics(mlt_D) < 0) then
stop 'do_mlt_eval: mlt_basics(mlt_D) < 0'
end if
if (MLT_dbg) write(*,*) 'done do_mlt_eval'
end subroutine do_mlt_eval
subroutine Get_results( &
have_L_need_gradT, cgrav, m, r, T, Rho, L, dg, P, &
chiRho, chiT, Cp, Cv, csound, xh, opacity, grada, dPrad_dm, &
gradr_factor, gradL_composition_term, &
alpha_semiconvection, semiconvection_option, &
thermohaline_coeff, thermohaline_option, &
dominant_iso_for_thermohaline, &
mixing_length_alpha, alt_scale_height, remove_small_D_limit, &
MLT_option, Henyey_y_param, Henyey_nu_param, &
prev_conv_vel, max_conv_vel, dt, tau, MLT_dbg, mixing_type, &
gradT, d_gradT_dvb, &
L_conv, d_L_conv_dvb, &
L_total, d_L_total_dvb, &
gradr, d_gradr_dvb, &
gradL, d_gradL_dvb, &
scale_height, d_scale_height_dvb, &
Lambda, d_Lambda_dvb, &
conv_vel, d_conv_vel_dvb, & ! convection velocity
D, d_D_dvb, &
Gamma, d_Gamma_dvb, &
conv_P, d_conv_P_dvb, &
ierr)
use utils_lib, only: is_bad_num
use chem_def, only: chem_isos
use atm_lib, only: eval_Paczynski_gradr
logical, intent(in) :: have_L_need_gradT
real(dp), intent(in) :: cgrav, m, r, T, Rho, L, dg, P
real(dp), intent(in) :: chiRho, chiT, Cp, Cv, csound, &
xh, opacity, dPrad_dm, grada, gradr_factor
real(dp), intent(in) :: gradL_composition_term
real(dp), intent(in) :: alpha_semiconvection, thermohaline_coeff
real(dp), intent(in) :: mixing_length_alpha
logical, intent(in) :: alt_scale_height
character (len=*), intent(in) :: &
MLT_option, thermohaline_option, semiconvection_option
integer, intent(in) :: dominant_iso_for_thermohaline
real(dp), intent(in) :: Henyey_y_param, Henyey_nu_param, &
prev_conv_vel, max_conv_vel, dt, tau, remove_small_D_limit
logical, intent(in) :: MLT_dbg
integer, intent(out) :: mixing_type
real(dp), intent(out) :: gradT, d_gradT_dvb(:)
real(dp), intent(out) :: L_conv, d_L_conv_dvb(:)
real(dp), intent(out) :: L_total, d_L_total_dvb(:)
real(dp), intent(out) :: gradr, d_gradr_dvb(:)
real(dp), intent(out) :: gradL, d_gradL_dvb(:)
real(dp), intent(out) :: scale_height, d_scale_height_dvb(:)
real(dp), intent(out) :: Lambda, d_Lambda_dvb(:)
real(dp), intent(out) :: conv_vel, d_conv_vel_dvb(:)
real(dp), intent(out) :: D, d_D_dvb(:)
real(dp), intent(out) :: Gamma, d_Gamma_dvb(:) ! convective efficiency
real(dp), intent(out) :: conv_P, d_conv_P_dvb(:)
integer, intent(out) :: ierr
real(dp) :: scale_height1, scale_height2
real(dp) :: Pg, Pr, dP_dvb(nvbs), dPg_dvb(nvbs), dPr_dvb(nvbs), dRho_dvb(nvbs)
real(dp) :: dT_dvb(nvbs), dL_dvb(nvbs), alpha, phi, dgrad, denom, tmp
real(dp) :: grav, d_grav_dvb(nvbs)
real(dp) :: diff_grads, d_diff_grads_dvb(nvbs)
real(dp) :: convective_conductivity, d_cc_dvb(nvbs)
real(dp) :: radiative_conductivity, d_rc_dvb(nvbs)
real(dp) :: surf, dsurf_dvb(nvbs)
real(dp) :: beta, d_beta_dvb(nvbs)
real(dp) :: chi, d_chi_dvb(nvbs)
real(dp) :: D_div_B, d_D_div_B_dvb(nvbs)
real(dp) :: Q, dQ_dvb(nvbs)
real(dp) :: A, dA_dvb(nvbs)
real(dp) :: Bcubed, d_Bcubed_dvb(nvbs)
real(dp) :: Zeta, d_Zeta_dvb(nvbs)
real(dp) :: d_Cp_dvb(nvbs)
real(dp) :: dR_dvb(nvbs)
real(dp) :: d_opacity_dvb(nvbs)
real(dp) :: d_grada_dvb(nvbs)
real(dp) :: Dconv, d_Dconv_dvb(nvbs)
real(dp) :: delta, d_delta_dvb(nvbs)
real(dp) :: f, f0, d_f0_dvb(nvbs)
real(dp) :: f1, d_f1_dvb(nvbs)
real(dp) :: f2, d_f2_dvb(nvbs)
real(dp) :: x, d_x_dvb(nvbs)
real(dp) :: d_chiT_dvb(nvbs), d_chiRho_dvb(nvbs)
real(dp) :: a0, omega, theta
real(dp) :: d_omega_dvb(nvbs), d_a0_dvb(nvbs), d_theta_dvb(nvbs)
integer :: i
real(dp), parameter :: tiny = 1d-30, min_D_th = 1d-3
character (len=256) :: message
logical :: quit
real(dp) :: diff_grad, K, gamma0, L_ratio, frac, s, sqrt_x, &
dilution_factor, conv_tau, init_conv_vel, area
real(dp) :: K_T, K_mu, nu_rad, nu_mol, nu, grad_mu, R0, r_th, H_P
logical, parameter :: debug = .false.
include 'formats.dek'
ierr = 0
call set_no_mixing
if (MLT_option == 'none') return
if (max_conv_vel == 0d0) return
if (alpha_semiconvection > 1d0) then
write(*,1) 'MLT alpha_semiconvection must be <= 1', alpha_semiconvection
ierr = -1
return
end if
d_grada_dvb = 0
d_grada_dvb(mlt_dgrada) = 1
grav = cgrav*m / (r*r)
d_grav_dvb = 0
d_grav_dvb(mlt_dlnR) = -2*grav
d_grav_dvb(mlt_dlnm) = grav
if (grav < 0) then
write(*,1) 'grav', grav
write(*,1) 'cgrav', cgrav
write(*,1) 'm', m
write(*,1) 'r', r
stop 'mlt'
end if
dP_dvb = 0
dP_dvb(mlt_dP) = 1
dRho_dvb = 0
dRho_dvb(mlt_dlnd) = Rho
dT_dvb = 0
dT_dvb(mlt_dlnT) = T
dL_dvb = 0
dL_dvb(mlt_dL) = 1
dR_dvb = 0
dR_dvb(mlt_dlnR) = r
d_chiT_dvb = 0
d_chiT_dvb(mlt_dchiT) = 1
d_chiRho_dvb = 0
d_chiRho_dvb(mlt_dchiRho) = 1
d_Cp_dvb = 0
d_Cp_dvb(mlt_dCp) = 1
d_opacity_dvb = 0
d_opacity_dvb(mlt_dopacity) = 1
surf = pi4*r*r
if (debug) write(*,1) 'surf', surf
dsurf_dvb = 8*pi*r*dR_dvb
Pr = one_third*crad*T*T*T*T
if (debug) write(*,1) 'Pr', Pr
dPr_dvb = 0
dPr_dvb(mlt_dlnT) = 4*Pr
! Ledoux temperature gradient (same as Schwarzschild if composition term = 0)
gradL = grada + gradL_composition_term
d_gradL_dvb = d_grada_dvb ! ignore partials of composition term
if (have_L_need_gradT) then
gradr = eval_Paczynski_gradr(P,opacity,L,m,cgrav,Pr,tau,T,r,rho)
gradr = gradr*gradr_factor
d_gradr_dvb(mlt_dP) = gradr/P
d_gradr_dvb(mlt_dopacity) = gradr/opacity
d_gradr_dvb(mlt_dL) = gradr_factor*P*opacity/(16*pi*clight*m*cgrav*Pr)
d_gradr_dvb(mlt_dlnm) = -gradr
d_gradr_dvb(mlt_dlnT) = -4*gradr
diff_grads = gradr - gradL ! convective if this is > 0
d_diff_grads_dvb = d_gradr_dvb - d_gradL_dvb
else
! dg = gradT - grada
! diff_grads = gradT - gradL = dg + grada - gradL
gradT = grada + dg
diff_grads = gradT - gradL ! convective if this is > 0
d_diff_grads_dvb = -d_gradL_dvb
d_diff_grads_dvb(mlt_dgradT) = 1
end if
! if (diff_grads > 0 .and. tiny_dg > 0) then
! ! smooth the cutoff as diff_grads goes to 0
! x = exp(-tiny_dg/diff_grads)
! d_diff_grads_dvb = x*d_diff_grads_dvb*(1 + tiny_dg/diff_grads)
! diff_grads = x*diff_grads
! end if
Pg = P - Pr
if (debug) write(*,1) 'Pg', Pg
if (Pg < tiny .or. r < tiny .or. opacity < tiny .or. rho < tiny .or. Cp < tiny) then
call set_no_mixing
return
end if
dPg_dvb = dP_dvb - dPr_dvb
beta = Pg/P
if (debug) write(*,1) 'beta', beta
d_beta_dvb = beta*(dPg_dvb/Pg - dP_dvb/P)
scale_height = P / (grav*rho)
d_scale_height_dvb = scale_height*(dP_dvb/P - d_grav_dvb/grav - dRho_dvb/Rho)
if (alt_scale_height) then
! consider sound speed*hydro time scale as an alternative scale height
! (this comes from Eggleton's code.)
scale_height2 = sqrt(P / cgrav) / rho
if (scale_height2 < scale_height) then
scale_height = scale_height2
d_scale_height_dvb = scale_height*(0.5d0*dP_dvb/P - dRho_dvb/Rho)
end if
end if
H_P = scale_height
if (debug) write(*,1) 'scale_height', scale_height
! mixing length, Lambda
Lambda = mixing_length_alpha*scale_height
if (debug) write(*,1) 'Lambda', Lambda
d_Lambda_dvb = mixing_length_alpha*d_scale_height_dvb
if (mixing_length_alpha <= 0) then
call set_no_mixing
return
end if
! 'Q' param C&G 14.24
Q = chiT/chiRho
dQ_dvb = Q*( d_chiT_dvb/chiT - d_chiRho_dvb/chiRho )
if (Q <= 0) then
call set_no_mixing
return
end if
radiative_conductivity = (4*crad*clight / 3)*T*T*T / (opacity*rho) ! erg / (K cm sec)
if (debug) write(*,1) 'radiative_conductivity', radiative_conductivity
d_rc_dvb = 0
d_rc_dvb(mlt_dlnd) = -radiative_conductivity
d_rc_dvb(mlt_dlnT) = 3*radiative_conductivity
d_rc_dvb(mlt_dopacity) = -radiative_conductivity/opacity
if (diff_grads <= 0d0) then ! not convective (Ledoux stable)
call set_no_mixing
gradr = gradT
if (gradL_composition_term < 0) then ! composition unstable
call set_thermohaline
else if (gradr > grada) then ! Schw unstable
call set_semiconvection
else
call set_no_mixing
end if
if (D < remove_small_D_limit) call set_no_mixing
if (.not. have_L_need_gradT) then ! L_total = L_rad
! calculate L_rad using dPrad/dm instead of gradT
! to avoid assumption of hydrostatic balance
area = 4*pi*r*r
L_total = -dPrad_dm*clight*area*area/opacity
d_L_total_dvb(mlt_d_dPrad_dm) = -clight*area*area/opacity
d_L_total_dvb(mlt_dlnR) = 4*L_total
d_L_total_dvb(mlt_dopacity) = -L_total/opacity
end if
call set_conv_P
return
end if
call set_convective_mixing
call set_conv_P
if (quit) return
D = conv_vel*Lambda/3 ! diffusion coefficient [cm^2/sec]
if (debug) write(*,1) 'D', D
d_D_dvb = (d_conv_vel_dvb*Lambda + conv_vel*d_Lambda_dvb)/3
mixing_type = convective_mixing
if (debug .or. D < 0) then
write(*,*) 'get_gradT: convective_mixing'
write(*,1) 'D', D
write(*,1) 'conv_vel', conv_vel
write(*,1) 'Pg/P', Pg/P
write(*,1) 'H_P', H_P
write(*,1) 'scale_height', scale_height
write(*,1) 'scale_height/H_P', scale_height/H_P
write(*,1) 'm/Msun', m/Msun
write(*,1) 'r/Rsun', r/Rsun
write(*,1) 'T', T
write(*,1) 'rho', rho
write(*,1) 'grada', grada
write(*,1) 'chiT', chiT
write(*,1) 'chiRho', chiRho
write(*,2) 'mixing_type', mixing_type
write(*,*)
stop 'MLT: get_gradT'
end if
if (D < remove_small_D_limit) call set_no_mixing
contains
#ifdef offload
!dir$ attributes offload: mic :: set_conv_P
#endif
subroutine set_conv_P
if (conv_vel == 0d0) then
conv_P = 0d0
d_conv_P_dvb = 0d0
else
conv_P = rho*conv_vel*conv_vel/(3*P) ! see C&G (14.69) for Pturb/P
d_conv_P_dvb = rho*2*d_conv_vel_dvb*conv_vel/(3*P) + &
dRho_dvb*conv_vel*conv_vel/(3*P) - conv_P*dP_dvb/P
end if
end subroutine set_conv_P
#ifdef offload
!dir$ attributes offload: mic :: set_thermohaline
#endif
subroutine set_thermohaline
real(dp) :: kipp_D, tau, Pr, BGS_C, Nu_mu, l2, lmda, phi, r_guess, &
sqrt_1_plus_phi, sqrt_Pr
logical, parameter :: dbg = .false.
include 'formats'
if (dbg) write(*,*) 'set_thermohaline ' // trim(thermohaline_option)
diff_grad = max(1d-40, grada - gradr) ! positive since Schwarzschild stable
K = 4*crad*clight*T*T*T/(3*opacity*rho) ! thermal conductivity
if (thermohaline_option == 'Kippenhahn') then
! Kippenhahn, R., Ruschenplatt, G., & Thomas, H.-C. 1980, A&A, 91, 175
D = -thermohaline_coeff*3*K/(2*rho*cp)*gradL_composition_term/diff_grad
else if (thermohaline_option == 'Brown_Garaud_Stellmach_13' .or. &
thermohaline_option == 'Traxler_Garaud_Stellmach_11') then
call get_diff_coeffs(K_T,K_mu,nu)
R0 = (gradr - grada)/gradL_composition_term
Pr = nu/K_T
tau = K_mu/K_T
r_th = (R0 - 1d0)/(1d0/tau - 1d0)
if (r_th >= 1d0) then ! stable if R0 >= 1/tau
D = 0d0
else if (thermohaline_option == 'Traxler_Garaud_Stellmach_11') then
! Traxler, Garaud, & Stellmach, ApJ Letters, 728:L29 (2011).
! also see Denissenkov. ApJ 723:563–579, 2010.
D = 101d0*sqrt(K_mu*nu)*exp_cr(-3.6d0*r_th)*pow_cr(1d0 - r_th,1.1d0) ! eqn 24
else
! if (thermohaline_option == 'Brown_Garaud_Stellmach_13') then
D = K_mu*Numu(R0,r_th,pr,tau)
endif
D = thermohaline_coeff*D
else
D = 0
ierr = -1
write(*,*) 'unknown value for MLT thermohaline_option' // trim(thermohaline_option)
return
end if
if (D < min_D_th) then
call set_no_mixing
return
end if
d_D_dvb = 0
conv_vel = min(max_conv_vel, 3*D/Lambda)
d_conv_vel_dvb = 0
mixing_type = thermohaline_mixing
return
end subroutine set_thermohaline
#ifdef offload
!dir$ attributes offload: mic :: get_diff_coeffs
#endif
subroutine get_diff_coeffs(kt,kmu,vis)
use chem_def, only: chem_isos
use crlibm_lib
real(dp) :: kt,kmu,vis,qe4
real(dp) :: loglambdah,loglambdacx,loglambdacy,ccx,ccy
real(dp) :: Bcoeff
real(dp) :: chemA,chemZ,acx,acy
real(dp), parameter :: sqrt5 = sqrt(5d0)
kt = K/(Cp*rho) ! thermal diffusivity (assumes radiatively dominated)
qe4=qe*qe*qe*qe
! Log Lambda for pure H (equation 10 from Proffitt Michaud 93)
loglambdah = -19.26 - 0.5*log_cr(rho) + 1.5*log_cr(T) - 0.5*log_cr(1. + 0.5*(1+xh))
nu_rad = 4*crad*T*T*T*T/(15*clight*opacity*rho*rho) ! radiative viscosity
nu_mol = 0.406*sqrt(amu)*pow_cr(boltzm*T,2.5d0)/(qe4*loglambdah*rho)
! From Spitzer "Physics of Fully Ionized Gases equation 5-54
! Assumes pure H. Still trying to work out what it would be for a mixture.
vis = nu_mol + nu_rad ! total viscosity
! The following is from Proffitt & Michaud, 1993.
! Their constant B (equation 15)
Bcoeff = (15./16.)*sqrt(2.*amu/(5*pi))*pow_cr(boltzm,2.5d0)/qe4
! Extract what species drives the thermohaline concvection
chemA = chem_isos%Z_plus_N(dominant_iso_for_thermohaline)
chemZ = chem_isos%Z(dominant_iso_for_thermohaline)
if(chemZ.gt.2) then
! This is if the driving chemical is NOT He.
! Log Lambda for H-dominant chem mixture (equation 10)
loglambdacx = loglambdah - log_cr(chemz)
! Log Lambda for He-dominant chem mixture (equation 10)
loglambdacy = loglambdah - log_cr(2.*chemz)
! Calculation of C_ij coeffs (equation 12)
ccx = log_cr(exp_cr(1.2*loglambdacx)+1.)/1.2
ccy = log_cr(exp_cr(1.2*loglambdacy)+1.)/1.2
! Reduced masses (I had to guess, from Bahcall & Loeb 1990), with H and He
acx = (1.*chemA)/(1.+chemA)
acy = 4*chemA/(4.+chemA)
! My formula (see notes) based on Proffitt and Michaud 1993
kmu = 2*Bcoeff*pow_cr(T,2.5d0)/(sqrt5*rho*chemZ*chemZ)/ &
(xh*sqrt(acx)*ccx + (1-xh)*sqrt(acy)*ccy)
else
! Log Lambda for H-He mixture (equation 10)
loglambdah = -19.26 - log_cr(2.d0) - 0.5*log_cr(rho) + &
1.5*log_cr(T) - 0.5*log_cr(1. + 0.5*(1+xh))
! Calculation of C_ij coeffs (equation 12)
ccy = log_cr(exp_cr(1.2*loglambdah)+1.)/1.2
! My formula (see notes) based on Proffitt and Michaud 1993
kmu = (Bcoeff*pow_cr(T,2.5d0)/(rho*ccy))*(3+xh)/((1+xh)*(3+5*xh)*(0.7+0.3*xh))
endif
! write(57,*) kt,kmu,vis,chemZ
end subroutine get_diff_coeffs
#ifdef offload
!dir$ attributes offload: mic :: numu
#endif
double precision function numu(R0,r_th,prandtl,diffratio)
!Function calculates Nu_mu from input parameters, following Brown et al. 2013.
!Written by P. Garaud (2013). Please email [email protected] for troubleshooting.
real(dp), intent(in) :: R0,r_th,prandtl,diffratio
real(dp) :: maxl2,maxl,lambdamax
real(dp) :: myvars(2)
integer :: ierr, iter
! Initialize guess using estimates from Brown et al. 2013
call analytical_estimate_th(maxl,lambdamax,r_th,prandtl,diffratio)
myvars(1) = maxl
myvars(2) = lambdamax
!Call Newton relaxation algorithm
call NR(myvars,prandtl,diffratio,R0,ierr)
!If the growth rate is negative, then try another set of parameters as first guess.
!Repeat as many times as necessary until convergence is obtained.
iter = 1
do while((myvars(2)<0).or.(ierr /= 0))
!write(*,*) 'Alternative', r_th,prandtl,diffratio,iter
!Reset guess values
myvars(1) = maxl
myvars(2) = lambdamax
!Call relaxation for slightly different Pr, tau, R0.
call NR(myvars,prandtl*(1.+iter*1.d-2),diffratio,R0/(1.+iter*1.d-2),ierr)
!If it converged this time, call NR for the real parameters.
if(ierr.eq.0) call NR(myvars,prandtl,diffratio,R0,ierr)
!write(*,*) prandtl,diffratio,R0,myvars(1),myvars(2),ierr
!Otherwise, increase counter and try again.
iter = iter + 1
enddo
!Plug solution into "l^2" and lambda.
maxl2 = myvars(1)*myvars(1)
lambdamax = myvars(2)
!write(*,*) prandtl,diffratio,r_th,maxl2,lambdamax
!Calculate Nu_mu using Formula (33) from Brown et al, with C = 7.
numu = 1. + 49.*lambdamax*lambdamax/(diffratio*maxl2*(lambdamax+diffratio*maxl2))
return
end function numu
#ifdef offload
!dir$ attributes offload: mic :: thermohaline_rhs
#endif
subroutine thermohaline_rhs(myx,myf,myj,prandtl,diffratio,R0)
! This routine is needed for the NR solver.
! Inputs the two following equations for lambda and maxl2:
! lambda^3 + a_2 lambda^2 + a_1 lambda + a_0 = 0 (eq. 19 of Brown et al.)
! b_2 lambda^2 + b_1 lambda + b_0 = 0 (eq. 20 of Brown et al.)
! Inputs f, the equations, and j, their jacobian.
! Written by P. Garaud (2013). Please email [email protected] for troubleshooting.
real(dp), intent(in) :: myx(2), prandtl, diffratio, R0
real(dp), intent(out) :: myf(2), myj(2,2)
real(dp) :: a_2,a_1,a_0,b_2,b_1,b_0,myterm,myx1_2,myx1_3,myx1_4
!This inputs the coefficients.
b_2 = 1.+prandtl+diffratio
myx1_2 = myx(1)*myx(1)
myx1_3 = myx1_2*myx(1)
myx1_4 = myx1_3*myx(1)
a_2 = myx1_2*b_2
myterm = diffratio*prandtl+prandtl+diffratio
b_1 = 2*myx1_2*myterm
a_1 = myx1_4*myterm + prandtl*(1. - (1./R0))
b_0 = 3.*myx1_4*diffratio*prandtl + prandtl*(diffratio - (1./R0))
a_0 = myx1_4*myx1_2*diffratio*prandtl + myx1_2*prandtl*(diffratio - (1./R0))
! write(*,*) a_2,a_1,a_0,b_2,b_1,b_0
!These are equations 19 and 20
myf(1) = ((myx(2) + a_2)*myx(2) + a_1)*myx(2) + a_0
myf(2) = b_2*myx(2)*myx(2) + b_1*myx(2) + b_0
!These are their Jacobians for the NR relaxation.
myj(1,1) = 2*myx(1)*b_2*myx(2)*myx(2) + &
4*myx1_3*myterm*myx(2) + 6*myx1_4*myx(1)*diffratio*prandtl &
+ 2*myx(1)*prandtl*(diffratio - (1./R0))
myj(1,2) = 3*myx(2)*myx(2) + 2*a_2*myx(2) + a_1
myj(2,1) = 4*myx(1)*myterm*myx(2) + 12.*myx1_3*diffratio*prandtl
myj(2,2) = 2*b_2*myx(2) + b_1
return
end subroutine thermohaline_rhs
#ifdef offload
!dir$ attributes offload: mic :: analytical_estimate_th
#endif
subroutine analytical_estimate_th(maxl,lambdamax,r_th,prandtl,diffratio)
!Inputs analytical estimates for l and lambda from Brown et al. 2013.
real(dp) :: prandtl, diffratio, maxl, lambdamax, r_th, phi, maxl4, maxl6
phi = diffratio/prandtl
if(r_th.lt.min(0.5*phi/(1.+phi),0.5*(1.+phi)*(1.+phi)*(1.+phi)/(4.*phi))) then
if(r_th.gt.prandtl) then
maxl = pow_cr((1./(1.+phi)) - 2.*dsqrt(r_th*phi)/pow_cr(1.+phi,2.5d0),0.25d0)
! Equation (B14)
maxl4 = maxl*maxl*maxl*maxl
maxl6 = maxl4*maxl*maxl
lambdamax = 2*prandtl*phi*maxl6/(1.-(1.+phi)*maxl4) ! Equation (B11)
else
maxl = dsqrt(dsqrt(1./(1.+phi)) - dsqrt(prandtl)*(1.+phi/((1.+phi)*(1.+phi))))
! Equation (B5)
lambdamax = dsqrt(prandtl) - prandtl*dsqrt(1.+phi) !Equation (B5)
endif
else
maxl = pow_cr((1./3.)*(1.-r_th) + (1-r_th)*(1-r_th)*(5-4*phi)/27.,0.25d0)
! Equation (B19) carried to next order (doesn't work well otherwise)
maxl4 = maxl*maxl*maxl*maxl
maxl6 = maxl4*maxl*maxl
lambdamax = 2*prandtl*phi*maxl6/(1.-(1.+phi)*maxl4) ! Equation (B11)
endif
if(lambdamax<0) then ! shouldn't be needed, but just as precaution
maxl = 0.5d0
lambdamax = 0.5d0
endif
return
end subroutine analytical_estimate_th
#ifdef offload
!dir$ attributes offload: mic :: NR
#endif
subroutine NR(xrk,prandtl,diffratio,R0,ierr)
! Newton Relaxation routine used to solve cubic & quadratic in thermohaline case.
! Written by P. Garaud (2013). Please email [email protected] for troubleshooting.
real(dp), parameter :: acy = 1.d-13 ! accuracy of NR solution.
integer, parameter :: niter = 20 ! max number of iterations allowed before giving up.
integer, parameter :: & !array dimension input parameters for dgesvx
n = 2, &
nrhs = 1, &
lda = n, &
ldaf = n, &
ldb = n, &
ldx = n
integer :: iter,ierr
real(dp) :: xrk(2), f(2) ! Functions f
real(dp) :: j(2,2) ! Jacobian
real(dp) :: err,errold ! Error at each iteration
real(dp) :: x1_sav,x2_sav
real(dp) :: prandtl, diffratio, R0
real(dp) :: A(lda,n), AF(ldaf,n), R(n), C(n), B(ldb,nrhs), X(ldx,nrhs), &
rcond, ferr(nrhs), berr(nrhs), work(4*n)
character :: fact, trans, equed
integer :: ipiv(n), iwork(n)
include 'formats'
!Initialize flags and other counters.
ierr = 0
iter = 0
err = 0d0
errold = 0d0
!Save input guess (probably not necessary here but useful in other routines)
x1_sav = xrk(1)
x2_sav = xrk(2)
!While error is too large .and. decreasing, iterate.
do while ((err.gt.acy).and.(ierr.eq.0).and.(iter.lt.niter))
call thermohaline_rhs(xrk,f,j,prandtl,diffratio,R0)
fact = 'E'
trans = 'N'
equed = ''
A = j
B(1,1) = f(1)
B(2,1) = f(2)
#ifdef offload
!dir$ attributes offload: mic :: dgesvx
#endif
call dgesvx( fact, trans, n, nrhs, A, lda, AF, ldaf, ipiv, &
equed, r, c, B, ldb, x, ldx, rcond, ferr, berr, &
work, iwork, ierr )
if (ierr /= 0) then
!write(*,*) 'dgesvx failed in thermohaline routine', iter
!write(*,2) j(1,1),j(1,2)
!write(*,2) j(2,1),j(2,2)
else
iter = iter + 1
f(1) = X(1,1)
f(2) = X(2,1)
err = dsqrt(f(1)*f(1)+f(2)*f(2)) ! Calculate the new error
! If, after a while, the error is still not decreasing, give up and exit NR.
! Otherwise, continue.
if((iter.gt.5).and.(err.gt.errold)) then
! Write(*,2) 'Error not decreasing at iter', iter, err, errold
ierr = 1
! Reset xs and exit loop.
xrk(1) = x1_sav
xrk(2) = x2_sav
else
xrk = xrk - f ! The solution is now in f, so update x
errold = err
endif
endif
enddo
!if(iter.eq.niter) write(*,2) 'Failed to converge'
return
end subroutine NR
#ifdef offload
!dir$ attributes offload: mic :: set_convective_mixing
#endif
subroutine set_convective_mixing
! need to set gradT, d_gradT_dvb, conv_vel, d_conv_vel_dvb
include 'formats.dek'
real(dp) ff1, ff2, ff3, ff4, ff5, aa, bb, y0, xres, a1, a2
real(dp) :: A_0, A_1, A_2, A_numerator, A_denom
real(dp), dimension(nvbs) :: &
dA_0_dvb, dA_1_dvb, dA_2_dvb, dA_numerator_dvb, dA_denom_dvb
real(dp), parameter :: two_13 = 1.2599210498948730d0 ! = pow_cr(2d0,1d0/3d0)
real(dp), parameter :: four_13 = 1.5874010519681994d0 ! = pow_cr(4d0,1d0/3d0)
! options for MLT_option are:
! 'ML1' Bohm-Vitense 1958 MLT
! 'ML2' Bohm and Cassinelli 1971 MLT
! 'Mihalas' Mihalas 1978, Kurucz 1979 MLT
! 'Henyey' Henyey, Rardya, and Bodenheimer 1965 MLT
! Values of the f1..f4 coefficients are taken from Table 1 of Ludwig et al. 1999, A&A, 346, 111
! with the following exception: their value of f3 for Henyey convection is f4/8 when it should be
! 8*f4, i.e., f3=32*pi**2/3 and f4=4*pi**2/3. f3 and f4 are related to the henyey y parameter, so
! for the 'Henyey' case they are set based on the value of Henyey_y_param. The f1..f4 parameters
! have been renamed with a double ff, i.e., ff1..ff4, to avoid a collision of variable names with
! f1 in the cubic root solver.
quit = .false.
x = Q*Rho / (2*P)
d_x_dvb = x*(drho_dvb/rho + dQ_dvb/Q - dP_dvb/P)
convective_conductivity = Cp*grav*Lambda*Lambda*Rho*(sqrt(x)) / 9 ! erg / (K cm sec)
if (convective_conductivity < 0) then
write(*,1) 'convective_conductivity', convective_conductivity
write(*,1) 'Cp', Cp
write(*,1) 'grav', grav
write(*,1) 'Lambda', Lambda
write(*,1) 'Rho', Rho
write(*,1) 'x', x
stop 'mlt'
end if
if (debug) write(*,1) 'convective_conductivity', convective_conductivity
d_cc_dvb = convective_conductivity* &
(d_Cp_dvb/Cp + d_grav_dvb/grav + &
2*d_Lambda_dvb/Lambda + dRho_dvb/rho + d_x_dvb/(2*x))
if (MLT_option == 'Cox') then ! this assumes optically thick
a0 = 9d0/4d0
! 'A' param is ratio of convective to radiative conductivities C&G 14.98
A = convective_conductivity / radiative_conductivity ! unitless.
if (debug) write(*,1) 'A', A
dA_dvb = (d_cc_dvb - d_rc_dvb*A) / radiative_conductivity
if (A < 0) then
write(*,1) 'A', A
write(*,1) 'convective_conductivity', convective_conductivity
write(*,1) 'radiative_conductivity', radiative_conductivity
stop 'mlt'
end if
else
select case(trim(MLT_option))
case ('Henyey')
ff1=1./Henyey_nu_param
ff2=1./2.
! popular values for y are 1/3 or 3/(4*pi**2)
ff3=8./Henyey_y_param
ff4=1./Henyey_y_param
case ('ML1')
ff1=1./8.
ff2=1./2.
ff3=24.0
ff4=0.0
case ('ML2')
ff1=1.
ff2=2.
ff3=16.0
ff4=0.0
case ('Mihalas')
ff1=1./8.
ff2=1./2.
ff3=16.0
ff4=2.0
case default
write(*,'(3a)') 'Error: ',trim(MLT_option), &
' is not an allowed MLT version for convection'
write(*,*)
return
end select
omega = Lambda*Rho*opacity !dimensionless
d_omega_dvb = omega*( d_Lambda_dvb/Lambda + dRho_dvb/Rho + d_opacity_dvb/opacity)
! the variable theta in no longer needed
! theta = omega / ( 1d0 + Henyey_y_param*omega**2 )
! d_theta_dvb = d_omega_dvb*(1d0 - Henyey_y_param*omega**2 ) /
! ( ( 1d0 + Henyey_y_param*omega**2 )**2 )
! a0 = 0.75d0*omega*theta
! d_a0_dvb = a0*( d_omega_dvb/omega + d_theta_dvb/theta )
a0 = (3./16.)*ff2*ff3/(1.+ff4/(omega*omega))
d_a0_dvb = a0*2*ff4*d_omega_dvb/(ff4 + omega*omega)/omega
! A = sqrt(P*Q*rho/Henyey_nu_param)*(Cp*mixing_length_alpha)/
! (2*crad*clight*T**3*theta)
A_0 = sqrt(ff1*P*Q*rho)
dA_0_dvb = ff1*(dP_dvb*Q*rho + P*dQ_dvb*rho + P*Q*drho_dvb)/(2*A_0)
A_1 = 4*A_0*Cp
dA_1_dvb = 4*(dA_0_dvb*Cp + A_0*d_Cp_dvb)
A_2 = mixing_length_alpha*omega*(1.+ff4/(omega*omega))
dA_2_dvb = mixing_length_alpha*(1-ff4/(omega*omega))*d_omega_dvb
A_numerator = A_1*A_2
dA_numerator_dvb = dA_1_dvb*A_2 + A_1*dA_2_dvb
A_denom = ff3*crad*clight*T*T*T
dA_denom_dvb = A_denom*3*dT_dvb/T
A = A_numerator/A_denom
dA_dvb = dA_numerator_dvb/A_denom - A_numerator*dA_denom_dvb/(A_denom*A_denom)
if (A < 0) then
write(*,1) 'A', A
write(*,1) 'A_numerator', A_numerator
write(*,1) 'A_denom', A_denom
write(*,1) 'A_1', A_1
write(*,1) 'ff3', ff3
write(*,1) 'A_0', A_0
stop 'mlt'
end if
end if
if (have_L_need_gradT) then ! C&G, section 14.7
! 'B' param C&G 14.81
Bcubed = (A*A / a0)*diff_grads
d_Bcubed_dvb = (A*A / a0)*d_diff_grads_dvb + (2*A*dA_dvb / a0)*diff_grads
if (debug) write(*,1) 'Bcubed', Bcubed
! now solve cubic equation for convective efficiency, Gamma
! a0*Gamma^3 + Gamma^2 + Gamma - a0*Bcubed == 0 C&G 14.82,
! rewritten in terms of Gamma
! leave it to Mathematica to find an expression for the root we want (with a0 = 9/4)
delta = a0*Bcubed
d_delta_dvb = a0*d_Bcubed_dvb !+ Bcubed*d_a0_dvb
if (debug) write(*,1) 'a0', a0
if (debug) write(*,1) 'delta', delta
f = -2 + 9*a0 + 27*a0*a0*delta !ignoring derivative wrt a0
if (f > 1d100) then
f0 = f
d_f0_dvb = 27*a0*a0*d_delta_dvb
else
f0 = sqrt(f*f + 4*(-1 + 3*a0)*(-1 + 3*a0)*(-1 + 3*a0))
d_f0_dvb = 27*a0*a0*f*d_delta_dvb / f0
end if
if (debug) write(*,1) 'f0', f0
f1 = -2 + 9*a0 + 27*a0*a0*delta + f0
if (f1 < 0 .or. is_bad_num(f1)) then
call set_no_mixing
quit = .true.
return
end if
f1 = pow_cr(f1,one_third)
d_f1_dvb = (27*a0*a0*d_delta_dvb + d_f0_dvb) / (3*f1*f1)
f2 = 2*two_13*(1 - 3*a0) / f1
d_f2_dvb = -f2*d_f1_dvb / f1
Gamma = (four_13*f1 + f2 - 2) / (6*a0)
d_Gamma_dvb = (four_13*d_f1_dvb + d_f2_dvb) / (6*a0)
if (is_bad_num(Gamma) .or. Gamma < 0) then
call set_no_mixing
quit = .true.
return
end if
! average convection velocity, vbar C&G 14.86b
! vbar = vsound*Sqrt(Q)*alpha*Gamma / (2*Sqrt(2*Gamma1)*A)
! vsound = Sqrt(Gamma1*P / rho), so
! vbar = Sqrt(Q*P / (8*rho))*alpha*Gamma / A
x = Q*P / (8*rho)
sqrt_x = sqrt(x)
conv_vel = sqrt_x*mixing_length_alpha*Gamma / A
if (conv_vel > max_conv_vel) then
conv_vel = max_conv_vel
d_conv_vel_dvb = 0
else
d_conv_vel_dvb = 0.5d0*conv_vel* &
(-2*dA_dvb / A + 2*d_Gamma_dvb / Gamma + &
dP_dvb / P + dQ_dvb / Q - drho_dvb / rho)
end if
call limit_conv_vel
Zeta = Gamma*Gamma*Gamma / Bcubed ! C&G 14.80
d_Zeta_dvb = Zeta*(3*d_Gamma_dvb / Gamma - d_Bcubed_dvb / Bcubed)
! Zeta must be >= 0 and < 1
if (Zeta < 0d0) then
Zeta = 0
d_Zeta_dvb = 0
else if (Zeta >= 1d0) then
Zeta = 1d0
d_Zeta_dvb = 0
end if
gradT = (1 - Zeta)*gradr + Zeta*gradL ! C&G 14.79
if (debug) then
write(*,1) 'Bcubed', Bcubed
write(*,1) 'Gamma', Gamma
write(*,1) 'Zeta', Zeta
write(*,1) 'gradT', gradT
end if
d_gradT_dvb = (1 - Zeta)*d_gradr_dvb + Zeta*d_gradL_dvb + &
(gradL - gradr)*d_Zeta_dvb
if (is_bad_num(gradT)) then
call set_no_mixing
quit = .true.
return
end if
else
x = 4*A*A*diff_grads
d_x_dvb = 4d0*A*(2d0*dA_dvb*diff_grads + A*d_diff_grads_dvb)
tmp = sqrt(1+x)
Gamma = 0.5d0*(tmp-1) ! C&G 14.106
if (debug) write(*,1) 'Gamma', Gamma
if (Gamma <= 0) then
call set_no_mixing
return
end if
d_Gamma_dvb = 0.25d0*d_x_dvb/tmp
conv_vel = sqrt(Q*P/(8d0*rho))*mixing_length_alpha*Gamma/A ! C&G 14.110
if (debug) write(*,1) 'conv_vel', conv_vel
if (conv_vel <= 0) then
call set_no_mixing
return
end if
if (conv_vel > max_conv_vel) then
conv_vel = max_conv_vel
d_conv_vel_dvb = 0
else
d_conv_vel_dvb = conv_vel* &
(0.5d0*(dP_dvb/P + dQ_dvb/Q - dRho_dvb/Rho) &
+ d_Gamma_dvb/Gamma - dA_dvb/A)
end if
call limit_conv_vel
L_conv = 0.5d0*surf*rho*conv_vel*cp*T*mixing_length_alpha*diff_grads ! C&G 14.18
if (debug) write(*,1) 'L_conv', L_conv
d_L_conv_dvb = L_conv* &
(dsurf_dvb/surf + dRho_dvb/Rho + d_conv_vel_dvb/conv_vel &
+ d_cp_dvb/cp + dT_dvb/T + d_diff_grads_dvb/diff_grads)
x = a0*Gamma*Gamma*Gamma/(A*A)
d_x_dvb = x*(3*d_Gamma_dvb/Gamma - 2*dA_dvb/A)
gradr = gradT + x ! C&G 14.107
if (debug) write(*,1) 'gradr', gradr
d_gradr_dvb = d_x_dvb
d_gradr_dvb(mlt_dgradT) = 1
L_total = L_conv*gradr/x ! C&G 14.109
if (debug) write(*,1) 'L_total', L_total
d_L_total_dvb = L_total*(d_L_conv_dvb/L_conv + d_gradr_dvb/gradr - d_x_dvb/x)
end if
end subroutine set_convective_mixing
#ifdef offload
!dir$ attributes offload: mic :: limit_conv_vel
#endif
subroutine limit_conv_vel
real(dp) :: inv_sqrt_x, inv_conv_tau
real(dp), dimension(nvbs) :: &
d_inv_sqrt_x_dvb, d_inv_conv_tau_dvb
include 'formats'
conv_tau = 1d99
init_conv_vel = conv_vel
if (dt > 0 .and. conv_vel > prev_conv_vel) then
! time limit increase in conv_vel; no infinite accelerations please!
conv_tau = scale_height/(conv_vel + prev_conv_vel)
! Arnett, W.D., 1969, Ap. and Space Sci, 5, 180.
! Wood, P.R., ApJ, 190:609-630, 1974. (Appendix V, eqns 1-3)
! modified to use scale_height instead of 0.5*Lambda
! to avoid cases of alpha set tiny for tiny convective region.
if (dt < conv_tau) then
inv_conv_tau = (conv_vel + prev_conv_vel)/scale_height
d_inv_conv_tau_dvb = ( &
d_conv_vel_dvb - inv_conv_tau*d_scale_height_dvb)/scale_height
conv_vel = prev_conv_vel + (conv_vel - prev_conv_vel)*inv_conv_tau*dt
d_conv_vel_dvb = ( &
d_conv_vel_dvb*inv_conv_tau + &
(conv_vel - prev_conv_vel)*d_inv_conv_tau_dvb)*dt
if (conv_vel > max_conv_vel) then
if (MLT_dbg) write(*,1) 'conv_vel > max_conv_vel', conv_vel, max_conv_vel
conv_vel = max_conv_vel
d_conv_vel_dvb = 0
end if
if (have_L_need_gradT) then
! recalculate Gamma to match modified conv_vel
if (A <= 1d-99 .or. sqrt_x <= 1d-99) then
Gamma = 1d25
d_Gamma_dvb = 0
if (MLT_dbg) write(*,1) 'A or sqrt_x too small', A, sqrt_x
else
inv_sqrt_x = 1d0/sqrt_x
d_inv_sqrt_x_dvb = &
(8d0*x*(drho_dvb*P - rho*dP_dvb) - P*P*dQ_dvb)/ &
(16d0*x*sqrt_x*P*rho)
Gamma = conv_vel*A*inv_sqrt_x/mixing_length_alpha
d_Gamma_dvb = ( &
d_conv_vel_dvb*A*inv_sqrt_x + &
conv_vel*dA_dvb*inv_sqrt_x + &
conv_vel*A*d_inv_sqrt_x_dvb)/mixing_length_alpha
end if
end if
end if
end if
if (MLT_dbg) &
write(*,1) 'prev/init init/final conv_vel', &
prev_conv_vel/init_conv_vel, &
prev_conv_vel, init_conv_vel, conv_vel
if (debug) write(*,1) 'conv_vel', conv_vel
if (conv_vel < 0) then
!$omp critical
write(*,1) 'conv_vel', conv_vel
write(*,1) 'mixing_length_alpha', mixing_length_alpha
write(*,1) 'x', x
write(*,1) 'A', A
write(*,1) 'Gamma', Gamma
stop 'MLT: set_convective_mixing'
!$omp end critical
end if
end subroutine limit_conv_vel
#ifdef offload
!dir$ attributes offload: mic :: set_semiconvection
#endif
subroutine set_semiconvection ! Langer 1983 & 1985
real(dp) :: alpha, bc, LG, &
a0, a1, a2, a3, a4, a5, a6, a, &
b1, b2, b3, b4, b5, b6, b7, b, div, bsq
real(dp), dimension(nvbs) :: &
d_bc_dvb, d_LG_dvb, d_a0_dvb, d_a1_dvb, d_a2_dvb, d_a3_dvb, d_a4_dvb, &
d_a5_dvb, d_a6_dvb, d_a_dvb, d_b1_dvb, d_b2_dvb, d_b3_dvb, d_b4_dvb, &
d_b5_dvb, d_b6_dvb, d_b7_dvb, d_b_dvb, d_div_dvb
include 'formats.dek'
if (MLT_dbg) write(*,*) 'check for semiconvection'
call set_no_mixing ! sets gradT = gradr
D = alpha_semiconvection*radiative_conductivity/(6*Cp*rho) &
*(gradr - grada)/(gradL - gradr)
if (D <= 0) then
if (MLT_dbg) then
write(*,1) 'set_no_mixing D', D
write(*,1) 'alpha_semiconvection', alpha_semiconvection
write(*,1) 'radiative_conductivity', radiative_conductivity
write(*,1) 'gradr - grada', gradr - grada
write(*,1) 'gradL - gradr', gradL - gradr
stop
end if
call set_no_mixing
return
end if
d_D_dvb = 0 ! not used, so skip for now.
conv_vel = 3*D/Lambda
d_conv_vel_dvb = 0
mixing_type = semiconvective_mixing
if (MLT_dbg) write(*,2) 'mixing_type', mixing_type
if (semiconvection_option == 'Langer_85 mixing; gradT = gradr') return
if (semiconvection_option /= 'Langer_85') then
write(*,*) 'MLT: unknown values for semiconvection_option ' // &
trim(semiconvection_option)
ierr = -1
return
end if
! Solve[{
! L/Lrad - Lsc/Lrad - 1 == 0,
! Lrad == grad LG,
! gradMu == (4 - 3*beta)/beta*gradL_composition_term,
! Lsc/Lrad == alpha (grad - gradA)/(2 grad (gradL - grad))
! (grad - gradA - (beta (8 - 3 beta))/bc gradMu)},
! grad, {Lsc, Lrad, gradMu}] // Simplify
alpha = alpha_semiconvection
bc = 32d0 - 24d0*beta - beta*beta
d_bc_dvb = - 24d0*d_beta_dvb - 2d0*d_beta_dvb*beta
LG = (16d0/3d0*pi*clight*m*cgrav*crad*T*T*T*T)/(P*opacity)
d_LG_dvb = 0
d_LG_dvb(mlt_dP) = -LG/P
d_LG_dvb(mlt_dopacity) = -LG/opacity
d_LG_dvb(mlt_dlnm) = LG
d_LG_dvb(mlt_dlnT) = 4d0*LG
a0 = alpha*gradL_composition_term*LG
d_a0_dvb = alpha*gradL_composition_term*d_LG_dvb
a1 = -2d0*bc*L
d_a1_dvb = -2d0*L*d_bc_dvb
d_a1_dvb(mlt_dL) = d_a1_dvb(mlt_dL) - 2d0*bc
a2 = 2d0*alpha*bc*grada*LG
d_a2_dvb = 2d0*alpha*(d_bc_dvb*grada*LG + bc*d_grada_dvb*LG + bc*grada*d_LG_dvb)
a3 = -2d0*bc*gradL*LG
d_a3_dvb = -2d0*(d_bc_dvb*gradL*LG + bc*d_gradL_dvb*LG + bc*gradL*d_LG_dvb)
a4 = 32d0*a0
d_a4_dvb = 32d0*d_a0_dvb
a5 = -36d0*beta*a0
d_a5_dvb = -36d0*(d_beta_dvb*a0 + beta*d_a0_dvb)
a6 = 9d0*beta*beta*a0
d_a6_dvb = 9d0*(2d0*beta*d_beta_dvb*a0 + beta*beta*d_a0_dvb)
a = a1 + a2 + a3 + a4 + a5 + a6
d_a_dvb = d_a1_dvb + d_a2_dvb + d_a3_dvb + d_a4_dvb + d_a5_dvb + d_a6_dvb
b1 = 32d0 - 36d0*beta + 9d0*beta*beta
d_b1_dvb = - 36d0*d_beta_dvb + 18d0*beta*d_beta_dvb
b2 = b1*a0
d_b2_dvb = d_b1_dvb*a0 + b1*d_a0_dvb
b3 = -2d0*gradL*L + alpha*grada*grada*LG
d_b3_dvb = -2d0*d_gradL_dvb*L + alpha*(2d0*grada*d_grada_dvb*LG + grada*grada*d_LG_dvb)
d_b3_dvb(mlt_dL) = d_b3_dvb(mlt_dL) - 2d0*gradL
b4 = (-alpha*gradA + gradL)*LG
d_b4_dvb = (-alpha*d_grada_dvb + d_gradL_dvb)*LG + (-alpha*gradA + gradL)*d_LG_dvb
b5 = -b2 + 2d0*bc*(L + b4)
d_b5_dvb = -d_b2_dvb + 2d0*d_bc_dvb*(L + b4) + 2d0*bc*d_b4_dvb
d_b5_dvb(mlt_dL) = d_b5_dvb(mlt_dL) + 2d0*bc
b6 = b2*grada + bc*b3
d_b6_dvb = d_b2_dvb*grada + b2*d_grada_dvb + d_bc_dvb*b3 + bc*d_b3_dvb
b7 = -4d0*(-2d0 + alpha)*bc*LG*b6
d_b7_dvb = -4d0*(-2d0 + alpha)*(d_bc_dvb*LG*b6 + bc*d_LG_dvb*b6 + bc*LG*d_b6_dvb)
b = b7 + b5*b5
d_b_dvb = d_b7_dvb + 2*b5*d_b5_dvb
div = 2d0*(-2d0 + alpha)*bc*LG
d_div_dvb = 2d0*(-2d0 + alpha)*(d_bc_dvb*LG + bc*d_LG_dvb)
bsq = sqrt(b)
gradT = (a + bsq)/div
d_gradT_dvb = -gradT*d_div_dvb/div + d_a_dvb/div + 0.5d0*d_b_dvb/(div*bsq)
end subroutine set_semiconvection
#ifdef offload
!dir$ attributes offload: mic :: show_partials
#endif
subroutine show_partials(str, val, partials)
character (len=*), intent(in) :: str
real(dp), intent(in) :: val, partials(nvbs)
integer :: i
include 'formats'
write(*,1) trim(str), val
i = mlt_dlnm
write(*,1) 'dlnm', partials(i)
i = mlt_dlnR
write(*,1) 'dlnR', partials(i)
i = mlt_dlnT
write(*,1) 'dlnT', partials(i)
i = mlt_dlnd
write(*,1) 'dlnd', partials(i)
i = mlt_dL
write(*,1) 'dL', partials(i)
i = mlt_dP
write(*,1) 'dP', partials(i)
i = mlt_dchiRho
write(*,1) 'dchiRho', partials(i)
i = mlt_dchiT
write(*,1) 'dchiT', partials(i)
i = mlt_dCp
write(*,1) 'dCp', partials(i)
i = mlt_dopacity
write(*,1) 'dopacity', partials(i)
i = mlt_dgrada
write(*,1) 'dgrada', partials(i)
write(*,*)
end subroutine show_partials
#ifdef offload
!dir$ attributes offload: mic :: set_no_mixing
#endif
subroutine set_no_mixing
mixing_type = no_mixing
if (have_L_need_gradT) then
gradT = gradr
d_gradT_dvb = d_gradr_dvb
else
L_conv = 0
d_L_conv_dvb = 0
end if
conv_vel = 0
d_conv_vel_dvb = 0
D = 0
d_D_dvb = 0
conv_P = 0
d_conv_P_dvb = 0
end subroutine set_no_mixing
#ifdef offload
!dir$ attributes offload: mic :: show_args
#endif
subroutine show_args
1 format(a30,1pe26.16)
write(*,1) 'cgrav = ', cgrav
write(*,1) 'm = ', m
write(*,1) 'r = ', r
write(*,1) 'T = ', T
write(*,1) 'Rho = ', Rho
write(*,1) 'L = ', L
write(*,1) 'P = ', P
write(*,1) 'chiRho = ', chiRho
write(*,1) 'chiT = ', chiT
write(*,1) 'Cp = ', Cp
write(*,1) 'Cv = ', Cv
write(*,1) 'csound = ', csound
write(*,1) 'xh = ', xh
write(*,1) 'opacity = ', opacity
write(*,1) 'grada = ', grada
write(*,1) 'mixing_length_alpha = ', mixing_length_alpha
end subroutine show_args
end subroutine Get_results
#ifdef offload
!dir$ end options
#endif
end module mlt
\end{lstlisting}
