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