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 :
21 : module semiconvection
22 :
23 : implicit none
24 :
25 : private
26 : public :: calc_semiconvection
27 :
28 : contains
29 :
30 : !> Calculates the outputs of semiconvective mixing theory.
31 : !!
32 : !! @param L Luminosity across a face (erg/s).
33 : !! @param Lambda The mixing length (cm).
34 : !! @param m Mass inside the face (g).
35 : !! @param T temperature (K).
36 : !! @param P pressure (erg/cm^3).
37 : !! @param Pr radiation pressure (erg/cm^3).
38 : !! @param beta ratio of gas pressure to radiation pressure.
39 : !! @param opacity opacity (cm^2/g).
40 : !! @param rho density (g/cm^3).
41 : !! @param alpha_semiconvection The semiconvective alpha parameter.
42 : !! @param semiconvection_option A string specifying which semiconvection theory to use. Currently supported are 'Langer_85 mixing; gradT = gradr' and 'Langer_85'.
43 : !! @param cgrav gravitational constant (erg*cm/g^2).
44 : !! @param Cp Specific heat at constant pressure (erg/g/K).
45 : !! @param gradr The radiative temperature gradient dlnT/dlnP_{rad}
46 : !! @param grada The adiabatic temperature gradient dlnT/dlnP|s
47 : !! @param gradL The Ledoux temperature gradient dlnT/dlnP
48 : !! @param gradL_composition_term The contribution of composition gradients to the Ledoux temperature gradient.
49 : !! @param gradT The temperature gradient dlnT/dlnP (output).
50 : !! @param Y_face The superadiabaticity (dlnT/dlnP - grada, output).
51 : !! @param conv_vel The convection speed (cm/s).
52 : !! @param D The chemical diffusion coefficient (cm^2/s).
53 : !! @param mixing_type Set to semiconvective if convection operates (output).
54 : !! @param ierr Tracks errors (output).
55 0 : subroutine calc_semiconvection(L, Lambda, m, T, P, Pr, beta, opacity, rho, alpha_semiconvection, &
56 : semiconvection_option, cgrav, Cp, gradr, grada, gradL, &
57 : gradL_composition_term, &
58 : gradT, Y_face, conv_vel, D, mixing_type, ierr) ! Langer 1983 & 1985
59 : use const_def, only: dp, pi, clight, crad, semiconvective_mixing
60 : use num_lib
61 : use utils_lib
62 : use auto_diff
63 : type(auto_diff_real_star_order1), intent(in) :: L, Lambda, T, P, Pr, beta, opacity, rho
64 : type(auto_diff_real_star_order1), intent(in) :: Cp, gradr, grada, gradL
65 : character(len=*), intent(in) :: semiconvection_option
66 : real(dp), intent(in) :: alpha_semiconvection, cgrav, gradL_composition_term, m
67 : type(auto_diff_real_star_order1), intent(out) :: gradT, Y_face, conv_vel, D
68 : integer, intent(out) :: mixing_type, ierr
69 :
70 : type(auto_diff_real_star_order1) :: bc, LG, &
71 : radiative_conductivity, a0, a1, a2, a3, a4, a5, a6, a, &
72 : b1, b2, b3, b4, b5, b6, b7, b, div, bsq
73 0 : real(dp) :: alpha
74 : include 'formats'
75 :
76 : ! Pre-compute common pieces
77 : radiative_conductivity = &
78 0 : (4d0/3d0*crad*clight)*pow3(T)/(opacity*rho) ! erg / (K cm sec)
79 : D = alpha_semiconvection*radiative_conductivity/(6d0*Cp*rho) &
80 0 : *(gradr - grada)/(gradL - gradr)
81 0 : if (D%val <= 0) return
82 0 : conv_vel = 3d0*D/Lambda
83 :
84 0 : if (semiconvection_option == 'Langer_85 mixing; gradT = gradr') then
85 0 : gradT = gradr
86 0 : Y_face = gradT - grada
87 0 : mixing_type = semiconvective_mixing
88 0 : return
89 0 : else if (semiconvection_option == 'Langer_85') then
90 : ! Solve[{
91 : ! L/Lrad - Lsc/Lrad - 1 == 0,
92 : ! Lrad == grad LG,
93 : ! gradMu == (4 - 3*beta)/beta*gradL_composition_term,
94 : ! Lsc/Lrad == alpha (grad - gradA)/(2 grad (gradL - grad))
95 : ! (grad - gradA - (beta (8 - 3 beta))/bc gradMu)},
96 : ! grad, {Lsc, Lrad, gradMu}] // Simplify
97 0 : alpha = min(1d0, alpha_semiconvection)
98 0 : bc = 32d0 - 24d0*beta - beta*beta
99 0 : LG = (16d0*pi*clight*m*cgrav*Pr)/(P*opacity)
100 0 : a0 = alpha*gradL_composition_term*LG
101 0 : a1 = -2d0*bc*L
102 0 : a2 = 2d0*alpha*bc*grada*LG
103 0 : a3 = -2d0*bc*gradL*LG
104 0 : a4 = 32d0*a0
105 0 : a5 = -36d0*beta*a0
106 0 : a6 = 9d0*beta*beta*a0
107 0 : a = a1 + a2 + a3 + a4 + a5 + a6
108 0 : b1 = 32d0 - 36d0*beta + 9d0*beta*beta
109 0 : b2 = b1*a0
110 0 : b3 = -2d0*gradL*L + alpha*grada*grada*LG
111 0 : b4 = (-alpha*gradA + gradL)*LG
112 0 : b5 = -b2 + 2d0*bc*(L + b4)
113 0 : b6 = b2*grada + bc*b3
114 0 : b7 = -4d0*(alpha - 2d0)*bc*LG*b6
115 0 : b = b7 + b5*b5
116 0 : div = 2d0*(alpha - 2d0)*bc*LG
117 0 : bsq = sqrt(b)
118 0 : gradT = (a + bsq)/div
119 0 : Y_face = gradT - grada
120 0 : conv_vel = 3d0*D/Lambda
121 0 : mixing_type = semiconvective_mixing
122 : else
123 : write(*,*) 'turb: unknown values for semiconvection_option ' // &
124 0 : trim(semiconvection_option)
125 0 : ierr = -1
126 0 : return
127 : end if
128 :
129 0 : end subroutine calc_semiconvection
130 :
131 : end module semiconvection
|
