Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 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 atm_utils
21 :
22 : use const_def, only: dp, pi, boltz_sigma, ln10, clight, crad
23 : use math_lib
24 :
25 : implicit none
26 :
27 : integer, parameter :: E2_NPAIRS = 571
28 :
29 : real(dp), target, save :: E2_x(E2_NPAIRS)
30 : real(dp), save :: E2_pairs(2*E2_NPAIRS)
31 : real(dp), target, save :: E2_f_ary(4*E2_NPAIRS)
32 : real(dp), pointer, save :: E2_f1(:), E2_f(:,:)
33 : logical, save :: have_E2_interpolant = .false.
34 :
35 : private
36 : public :: init
37 : public :: shutdown
38 : public :: eval_Teff_g
39 : public :: eval_Paczynski_gradr
40 : public :: eval_E2
41 :
42 : contains
43 :
44 5 : subroutine init(use_cache, ierr)
45 :
46 : use table_atm, only: table_atm_init
47 :
48 : logical, intent(in) :: use_cache
49 : integer, intent(out) :: ierr
50 :
51 : ierr = 0
52 :
53 3 : E2_f1 => E2_f_ary
54 3 : E2_f(1:4,1:E2_NPAIRS) => E2_f1(1:4*E2_NPAIRS)
55 :
56 3 : call table_atm_init(use_cache, ierr)
57 :
58 : include 'e2_pairs.dek'
59 :
60 3 : end subroutine init
61 :
62 :
63 1 : subroutine shutdown()
64 :
65 3 : use table_atm, only : table_atm_shutdown
66 :
67 1 : call table_atm_shutdown()
68 :
69 1 : end subroutine shutdown
70 :
71 :
72 2 : subroutine eval_Teff_g(L, R, M, cgrav, Teff, g)
73 :
74 : real(dp), intent(in) :: L
75 : real(dp), intent(in) :: R
76 : real(dp), intent(in) :: M
77 : real(dp), intent(in) :: cgrav
78 : real(dp), intent(out) :: Teff
79 : real(dp), intent(out) :: g
80 :
81 : ! Evaluate the effective temperature and surface gravity
82 :
83 2 : Teff = pow(L/(4._dp*pi*R*R*boltz_sigma), 0.25_dp)
84 :
85 2 : g = cgrav * M / (R*R)
86 :
87 1 : end subroutine eval_Teff_g
88 :
89 :
90 0 : function eval_Paczynski_gradr( &
91 0 : T, P, rho, tau, kap, L, M, R, cgrav) result (gradr)
92 :
93 : use eos_lib, only: radiation_pressure
94 :
95 : real(dp), intent(in) :: T
96 : real(dp), intent(in) :: P
97 : real(dp), intent(in) :: rho
98 : real(dp), intent(in) :: tau
99 : real(dp), intent(in) :: kap
100 : real(dp), intent(in) :: L
101 : real(dp), intent(in) :: R
102 : real(dp), intent(in) :: M
103 : real(dp), intent(in) :: cgrav
104 : real(dp) :: gradr
105 :
106 0 : real(dp) :: Prad
107 0 : real(dp) :: dilution_factor
108 0 : real(dp) :: s
109 0 : real(dp) :: f
110 :
111 : ! Evaluate the radiative temperature gradient, using expressions
112 : ! from Paczynski (1969, Act Ast, 19, 1)
113 :
114 0 : Prad = radiation_pressure(T)
115 :
116 0 : gradr = P*kap*L / (16._dp*pi*clight*M*cgrav*Prad)
117 :
118 0 : if (tau < 2._dp/3._dp) then ! Eqn. 8
119 :
120 : ! Eqn. 15
121 :
122 0 : s = (2._dp*crad*T*T*T*SQRT(R))/(3._dp*cgrav*M*rho)*pow(L/(8._dp*pi*boltz_sigma), 0.25_dp)
123 :
124 : ! Eqn. 8
125 :
126 0 : f = 1._dp - 1.5_dp*tau
127 :
128 0 : dilution_factor = (1._dp + f*s*(4._dp*pi*cgrav*clight*M)/(kap*L))/(1._dp + f*s)
129 :
130 0 : gradr = gradr*dilution_factor
131 :
132 : end if
133 :
134 : return
135 :
136 0 : end function eval_Paczynski_gradr
137 :
138 16 : subroutine eval_E2(x, E2, dE2_dx, ierr)
139 :
140 0 : use interp_1d_lib
141 : use interp_1d_def
142 :
143 : real(dp), intent(in) :: x
144 : real(dp), intent(out) :: E2
145 : real(dp), intent(out) :: dE2_dx
146 : integer, intent(out) :: ierr
147 :
148 16 : real(dp) :: val
149 16 : real(dp) :: slope
150 :
151 : ierr = 0
152 :
153 : ! Evaluate the E2 exponential integral
154 :
155 16 : if (.not. have_E2_interpolant) then
156 2 : call create_E2_interpolant(ierr)
157 2 : if (ierr /= 0) return
158 : end if
159 :
160 16 : call interp_value_and_slope(E2_x, E2_NPAIRS, E2_f1, x, val, slope, ierr)
161 16 : if (ierr /= 0) then
162 0 : write(*,*) 'call to interp_value_and_slope failed in eval_E2'
163 0 : return
164 : end if
165 :
166 : ! val = log10[E2]
167 16 : E2 = exp10(val)
168 16 : dE2_dx = slope*ln10*E2
169 :
170 16 : end subroutine eval_E2
171 :
172 2 : subroutine create_E2_interpolant(ierr)
173 :
174 16 : use interp_1d_lib
175 : use interp_1d_def
176 :
177 : use atm_def
178 :
179 : integer, intent(out) :: ierr
180 :
181 : integer, parameter :: NWORK = pm_work_size
182 :
183 3428 : real(dp), target :: work_ary(E2_NPAIRS*NWORK)
184 : real(dp), pointer :: work(:)
185 : integer :: i
186 :
187 2 : ierr = 0
188 :
189 : ! Set up an interpolating spline for the E2 exponential integral
190 :
191 2 : work => work_ary
192 :
193 1144 : do i = 1, E2_NPAIRS
194 1142 : E2_x(i) = E2_pairs(2*i-1)
195 1144 : E2_f(1,i) = E2_pairs(2*i)
196 : end do
197 :
198 2 : call interp_pm(E2_x, E2_NPAIRS, E2_f1, nwork, work, 'atm_utils', ierr)
199 2 : if (ierr /= 0) then
200 0 : write(*,*) 'call to interp_pm failed in create_E2_interpolant'
201 : end if
202 :
203 2 : have_E2_interpolant = .true.
204 :
205 2 : end subroutine create_E2_interpolant
206 :
207 : end module atm_utils
|
