Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2021 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module MLT
21 :
22 : implicit none
23 :
24 : private
25 : public :: calc_MLT
26 :
27 : contains
28 :
29 : !> Calculates the outputs of convective mixing length theory.
30 : !!
31 : !! @param MLT_option A string specifying which MLT option to use. Currently supported are Cox, Henyey, ML1, ML2, Mihalas. Note that 'TDC' is also a valid input and will return the Cox result. This is for use when falling back from TDC -> MLT, as Cox is the most-similar prescription to TDC.
32 : !! @param mixing_length_alpha The mixing length parameter.
33 : !! @param Henyey_MLT_nu_param The nu parameter in Henyey's MLT prescription.
34 : !! @param Henyey_MLT_y_param The y parameter in Henyey's MLT prescription.
35 : !! @param chiT dlnP/dlnT|rho
36 : !! @param chiRho dlnP/dlnRho|T
37 : !! @param Cp Specific heat at constant pressure (erg/g/K).
38 : !! @param grav The acceleration due to gravity (cm/s^2).
39 : !! @param Lambda The mixing length (cm).
40 : !! @param rho density (g/cm^3).
41 : !! @param T temperature (K).
42 : !! @param opacity opacity (cm^2/g)
43 : !! @param gradr The radiative temperature gradient dlnT/dlnP_{rad}
44 : !! @param grada The adiabatic temperature gradient dlnT/dlnP|s
45 : !! @param gradL The Ledoux temperature gradient dlnT/dlnP
46 : !! @param Gamma The convective efficiency parameter (output).
47 : !! @param gradT The temperature gradient dlnT/dlnP (output).
48 : !! @param Y_face The superadiabaticity (dlnT/dlnP - grada, output).
49 : !! @param conv_vel The convection speed (cm/s).
50 : !! @param D The chemical diffusion coefficient (cm^2/s).
51 : !! @param mixing_type Set to convective if convection operates (output).
52 : !! @param ierr Tracks errors (output).
53 67974 : subroutine calc_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, &
54 : Henyey_MLT_y_param, chiT, chiRho, Cp, grav, Lambda, rho, &
55 : P, T, opacity, gradr, grada, gradL, Gamma, gradT, Y_face, &
56 : conv_vel, D, mixing_type, max_conv_vel, ierr)
57 : use const_def, only: dp, clight, convective_mixing, crad, two_13, four_13, one_third
58 : use num_lib
59 : use utils_lib
60 : use auto_diff
61 : type(auto_diff_real_star_order1), intent(in) :: chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, gradr, grada, gradL
62 : character(len=*), intent(in) :: MLT_option
63 : real(dp), intent(in) :: mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel
64 : type(auto_diff_real_star_order1), intent(out) :: Gamma, gradT, Y_face, conv_vel, D
65 : integer, intent(out) :: mixing_type, ierr
66 :
67 67974 : real(dp) :: ff1, ff2, ff3, ff4
68 : type(auto_diff_real_star_order1) :: &
69 : Q, omega, a0, ff4_omega2_plus_1, A_1, A_2, &
70 : A_numerator, A_denom, A, Bcubed, delta, Zeta, &
71 : f, f0, f1, f2, radiative_conductivity, convective_conductivity, csound
72 : include 'formats'
73 67974 : if (gradr > gradL) then
74 : ! Convection zone
75 :
76 11713 : Q = chiT/chiRho ! 'Q' param C&G 14.24
77 11713 : if (MLT_option == 'Cox' .or. MLT_option == 'TDC') then ! this assumes optically thick
78 11713 : a0 = 9d0/4d0
79 : convective_conductivity = &
80 11713 : Cp*grav*pow2(Lambda)*rho*(sqrt(Q*rho/(2d0*P)))/9d0 ! erg / (K cm sec)
81 : radiative_conductivity = &
82 11713 : (4d0/3d0*crad*clight)*pow3(T)/(opacity*rho) ! erg / (K cm sec)
83 11713 : A = convective_conductivity / radiative_conductivity ! unitless.
84 : else
85 0 : select case(trim(MLT_option))
86 : case ('Henyey')
87 0 : ff1=1.0d0/Henyey_MLT_nu_param
88 0 : ff2=0.5d0
89 0 : ff3=8.0d0/Henyey_MLT_y_param
90 0 : ff4=1.0d0/Henyey_MLT_y_param
91 : case ('ML1')
92 0 : ff1=0.125d0
93 0 : ff2=0.5d0
94 0 : ff3=24.0d0
95 0 : ff4=0.0d0
96 : case ('ML2')
97 0 : ff1=1.0d0
98 0 : ff2=2.0d0
99 0 : ff3=16.0d0
100 0 : ff4=0.0d0
101 : case ('Mihalas')
102 0 : ff1=0.125d0
103 0 : ff2=0.5d0
104 0 : ff3=16.0d0
105 0 : ff4=2.0d0
106 : case default
107 0 : write(*,'(3a)') 'Error: ', trim(MLT_option), ' is not an allowed MLT option'
108 0 : call mesa_error(__FILE__,__LINE__)
109 : end select
110 :
111 0 : omega = Lambda*rho*opacity
112 0 : ff4_omega2_plus_1 = ff4/pow2(omega) + 1d0
113 0 : a0 = (3d0/16d0)*ff2*ff3/ff4_omega2_plus_1
114 0 : A_1 = 4d0*Cp*sqrt(ff1*P*Q*rho)
115 0 : A_2 = mixing_length_alpha*omega*ff4_omega2_plus_1
116 0 : A_numerator = A_1*A_2
117 0 : A_denom = ff3*crad*clight*pow3(T)
118 0 : A = A_numerator/A_denom
119 : end if
120 :
121 : ! 'B' param C&G 14.81
122 11713 : Bcubed = (pow2(A)/a0)*(gradr - gradL)
123 :
124 : ! now solve cubic equation for convective efficiency, Gamma
125 : ! a0*Gamma^3 + Gamma^2 + Gamma - a0*Bcubed == 0 C&G 14.82,
126 : ! leave it to Mathematica to find an expression for the root we want
127 11713 : delta = a0*Bcubed
128 11713 : f = -2d0 + 9d0*a0 + 27d0*a0*a0*delta
129 11713 : if (f > 1d100) then
130 0 : f0 = f
131 : else
132 11713 : f0 = pow2(f) + 4d0*(-1d0 + 3d0*a0)*(-1d0 + 3d0*a0)*(-1d0 + 3d0*a0)
133 11713 : if (f0 <= 0d0) then
134 0 : f0 = f
135 : else
136 11713 : f0 = sqrt(f0)
137 : end if
138 : end if
139 11713 : f1 = -2d0 + 9d0*a0 + 27d0*a0*a0*delta + f0
140 11713 : if (f1 <= 0d0) return
141 11713 : f1 = pow(f1,one_third)
142 11713 : f2 = 2d0*two_13*(1d0 - 3d0*a0) / f1
143 11713 : Gamma = (four_13*f1 + f2 - 2d0) / (6d0*a0)
144 :
145 11713 : if (Gamma <= 0d0) return
146 :
147 : ! average convection velocity C&G 14.86b
148 11713 : conv_vel = mixing_length_alpha*sqrt(Q*P/(8d0*rho))*Gamma / A
149 :
150 : ! convective velocity limiter
151 11713 : if (conv_vel%val > max_conv_vel) conv_vel%val = max_conv_vel
152 :
153 11713 : D = conv_vel*Lambda/3d0 ! diffusion coefficient [cm^2/sec]
154 :
155 : !Zeta = pow3(Gamma)/Bcubed ! C&G 14.80
156 11713 : Zeta = exp(3d0*log(Gamma) - log(Bcubed)) ! write it this way to avoid overflow problems
157 : else
158 : ! Radiative zone, because this means that gradr < gradL
159 56261 : Gamma = -1d99
160 56261 : Zeta = 0d0
161 56261 : conv_vel = 0d0
162 56261 : D = 0d0
163 : end if
164 :
165 : ! Zeta must be >= 0 and <= 1.
166 : ! By construction (above) it cannot be less than zero,
167 : ! so we just check that it is a valid number and is not greater
168 : ! than one.
169 67974 : if (is_bad(Zeta%val)) return
170 67974 : if (Zeta > 1d0) then
171 0 : Zeta = 1d0
172 : end if
173 :
174 67974 : gradT = (1d0 - Zeta)*gradr + Zeta*gradL ! C&G 14.79
175 67974 : Y_face = gradT - gradL
176 :
177 67974 : if (Y_face > 0d0) then
178 11713 : mixing_type = convective_mixing
179 : end if
180 :
181 67974 : end subroutine calc_MLT
182 :
183 : end module MLT
|
