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_irradiated
21 :
22 : use const_def, only: dp
23 : use math_lib
24 :
25 : implicit none
26 :
27 : private
28 :
29 : public :: eval_irradiated
30 :
31 : contains
32 :
33 : ! Evaluate irradiated atmosphere data with fixed, uniform opacity
34 :
35 18 : subroutine eval_irradiated( &
36 : L, R, M, cgrav, T_eq, P_surf, kap_guess, kap_v, gamma, &
37 : eos_proc, kap_proc, errtol, max_iters, skip_partials, &
38 : Teff, kap, tau, &
39 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
40 : ierr)
41 :
42 : use atm_def, only: atm_eos_iface, atm_kap_iface
43 : use atm_utils, only: eval_Teff_g
44 : use eos_def, only: num_eos_basic_results, i_chiRho, i_chiT
45 :
46 : real(dp), intent(in) :: L
47 : real(dp), intent(in) :: R
48 : real(dp), intent(in) :: M
49 : real(dp), intent(in) :: cgrav
50 : real(dp), intent(in) :: T_eq
51 : real(dp), intent(in) :: P_surf
52 : real(dp), intent(in) :: kap_guess
53 : real(dp), intent(in) :: kap_v
54 : real(dp), intent(in) :: gamma
55 : procedure(atm_eos_iface) :: eos_proc
56 : procedure(atm_kap_iface) :: kap_proc
57 : real(dp), intent(in) :: errtol
58 : integer, intent(in) :: max_iters
59 : logical, intent(in) :: skip_partials
60 : real(dp), intent(in) :: Teff
61 : real(dp), intent(out) :: kap
62 : real(dp), intent(out) :: tau
63 : real(dp), intent(out) :: lnT
64 : real(dp), intent(out) :: dlnT_dL
65 : real(dp), intent(out) :: dlnT_dlnR
66 : real(dp), intent(out) :: dlnT_dlnM
67 : real(dp), intent(out) :: dlnT_dlnkap
68 : integer, intent(out) :: ierr
69 :
70 : real(dp) :: T_int
71 2 : real(dp) :: g
72 2 : real(dp) :: lnP
73 : integer :: iters
74 2 : real(dp) :: lnRho
75 54 : real(dp) :: res(num_eos_basic_results)
76 54 : real(dp) :: dres_dlnRho(num_eos_basic_results)
77 54 : real(dp) :: dres_dlnT(num_eos_basic_results)
78 2 : real(dp) :: dlnkap_dlnRho
79 2 : real(dp) :: dlnkap_dlnT
80 2 : real(dp) :: kap_prev
81 2 : real(dp) :: err
82 2 : real(dp) :: chiRho
83 2 : real(dp) :: chiT
84 2 : real(dp) :: dlnkap_dlnT_P
85 :
86 : include 'formats'
87 :
88 2 : ierr = 0
89 :
90 : ! Sanity checks
91 :
92 2 : if (L <= 0._dp .OR. R <= 0._dp .OR. M <= 0._dp) then
93 0 : ierr = -1
94 0 : return
95 : end if
96 :
97 : ! Evaluate the 'interior' temperature & gravity
98 :
99 2 : call eval_Teff_g(L, R, M, cgrav, T_int, g)
100 :
101 : ! Evaluate atmosphere data using kap_guess as the opacity
102 :
103 2 : kap = kap_guess
104 :
105 : call eval_data( &
106 : T_int, g, L, T_eq, P_surf, kap, kap_v, gamma, skip_partials, &
107 : tau, lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
108 2 : ierr)
109 :
110 : ! Iterate to find a consistent opacity
111 :
112 2 : lnP = log(P_surf)
113 :
114 16 : iterate_loop : do iters = 1, max_iters
115 :
116 : ! Calculate the density & eos results
117 :
118 : call eos_proc( &
119 : lnP, lnT, &
120 : lnRho, res, dres_dlnRho, dres_dlnT, &
121 16 : ierr)
122 16 : if (ierr /= 0) then
123 0 : write(*,*) 'Call to eos_proc failed in eval_T_tau_uniform'
124 0 : return
125 : end if
126 :
127 : ! Update the opacity
128 :
129 16 : kap_prev = kap
130 :
131 : call kap_proc( &
132 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
133 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
134 16 : ierr)
135 16 : if (ierr /= 0) then
136 0 : write(*,*) 'Call to kap_proc failed in eval_T_tau_uniform'
137 0 : return
138 : end if
139 :
140 : ! Check for convergence
141 :
142 16 : err = abs(kap_prev - kap)/(errtol + errtol*kap)
143 :
144 16 : if (err < 1._dp) exit iterate_loop
145 :
146 14 : kap = kap_prev + 0.5_dp*(kap - kap_prev) ! under correct
147 :
148 : ! Re-evaluate atmosphere data
149 :
150 : call eval_data( &
151 : T_int, g, L, T_eq, P_surf, kap, kap_v, gamma, skip_partials, &
152 : tau, lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
153 16 : ierr)
154 :
155 : end do iterate_loop
156 :
157 2 : if (max_iters > 0 .AND. iters > max_iters) then
158 0 : write(*,*) 'Exceeded max_iters iterations in eval_irradiated'
159 0 : ierr = -1
160 0 : return
161 : end if
162 :
163 : ! If necessary, fix up the partials to account for the implicit
164 : ! dependence of the opacity on the final T
165 :
166 2 : if (max_iters > 0 .AND. .NOT. skip_partials) then
167 :
168 0 : chiRho = res(i_chiRho)
169 0 : chiT = res(i_chiT)
170 :
171 0 : dlnkap_dlnT_P = dlnkap_dlnT - dlnkap_dlnRho*chiT/chiRho
172 :
173 0 : dlnT_dL = dlnT_dL/(1._dp - dlnkap_dlnT_P*dlnT_dlnkap)
174 0 : dlnT_dlnR = dlnT_dlnR/(1._dp - dlnkap_dlnT_P*dlnT_dlnkap)
175 0 : dlnT_dlnM = dlnT_dlnM/(1._dp - dlnkap_dlnT_P*dlnT_dlnkap)
176 :
177 0 : dlnT_dlnkap = 0._dp
178 :
179 : end if
180 :
181 : ! Set the effective temperature. This is equal to T_int, because
182 : ! irradiation has no effect on the *net* flux emerging from the
183 : ! atmosphere
184 :
185 : ! Teff = T_int
186 :
187 : return
188 :
189 2 : end subroutine eval_irradiated
190 :
191 :
192 : ! Evaluate atmosphere data
193 :
194 16 : subroutine eval_data( &
195 : T_int, g, L, T_eq, P_surf, kap, kap_v, gamma, skip_partials, &
196 : tau, lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
197 : ierr)
198 :
199 2 : use atm_utils, only: eval_E2
200 :
201 : real(dp), intent(in) :: T_int
202 : real(dp), intent(in) :: g
203 : real(dp), intent(in) :: L
204 : real(dp), intent(in) :: T_eq
205 : real(dp), intent(in) :: P_surf
206 : real(dp), intent(in) :: kap
207 : real(dp), intent(in) :: kap_v
208 : real(dp), intent(in) :: gamma
209 : logical, intent(in) :: skip_partials
210 : real(dp), intent(out) :: tau
211 : real(dp), intent(out) :: lnT
212 : real(dp), intent(out) :: dlnT_dL
213 : real(dp), intent(out) :: dlnT_dlnR
214 : real(dp), intent(out) :: dlnT_dlnM
215 : real(dp), intent(out) :: dlnT_dlnkap
216 : integer, intent(out) :: ierr
217 :
218 16 : real(dp) :: gamma_eff
219 : real(dp) :: x
220 : real(dp) :: E2
221 : real(dp) :: dE2_dx
222 16 : real(dp) :: f_1
223 16 : real(dp) :: f_2
224 16 : real(dp) :: T4_int
225 16 : real(dp) :: T4_eq
226 16 : real(dp) :: T4
227 16 : real(dp) :: df_1_dtau
228 16 : real(dp) :: df_2_dtau
229 16 : real(dp) :: df_1_dx
230 16 : real(dp) :: df_2_dx
231 16 : real(dp) :: dlnT_dlnT_int
232 16 : real(dp) :: dlnT_dlntau
233 16 : real(dp) :: dlnT_dlnx
234 16 : real(dp) :: dlnT_int_dlnL
235 16 : real(dp) :: dlnT_int_dlnR
236 16 : real(dp) :: dlnT_int_dlnM
237 16 : real(dp) :: dlnT_int_dlnkap
238 16 : real(dp) :: dlntau_dlnL
239 16 : real(dp) :: dlntau_dlnR
240 16 : real(dp) :: dlntau_dlnM
241 16 : real(dp) :: dlntau_dlnkap
242 16 : real(dp) :: dlnx_dlnL
243 16 : real(dp) :: dlnx_dlnR
244 16 : real(dp) :: dlnx_dlnM
245 16 : real(dp) :: dlnx_dlnkap
246 :
247 : ! Calculate the optical depth corresponding to P_surf
248 :
249 16 : tau = P_surf*kap/g
250 :
251 : ! Evaluate irradiation terms [cf. eq. 6 of Guillot & Havel (2011,
252 : ! A&A 527, A20)]
253 :
254 16 : if (gamma > 0._dp) then
255 : gamma_eff = gamma
256 : else
257 16 : gamma_eff = kap_v/kap
258 : end if
259 :
260 16 : x = gamma_eff*tau
261 :
262 16 : call eval_E2(x, E2, dE2_dx, ierr)
263 16 : if (ierr /= 0) return
264 :
265 16 : f_1 = 2._dp*(1._dp + (x/2._dp - 1._dp)*exp(-x))*tau/(3._dp*x)
266 16 : f_2 = 2._dp*x*(1 - tau*tau/2._dp)*E2/(3._dp*tau)
267 :
268 : ! Evaluate the temperature
269 :
270 16 : T4_int = T_int*T_int*T_int*T_int
271 16 : T4_eq = T_eq*T_eq*T_eq*T_eq
272 :
273 16 : T4 = 0.75_dp*(T4_int*(tau + 2._dp/3._dp) + T4_eq*(f_1 + f_2 + 2._dp/3._dp))
274 :
275 16 : lnT = 0.25_dp*log(T4)
276 :
277 : ! Set up partials
278 :
279 16 : if (.NOT. skip_partials) then
280 :
281 0 : df_1_dtau = (2._dp + (x - 2._dp)*exp(-x))/(3._dp*x)
282 0 : df_2_dtau = -x*(2._dp + tau*tau)*E2/(3._dp*tau*tau)
283 :
284 0 : df_1_dx = -(2._dp - (2._dp + (2._dp - x)*x)*exp(-x))*tau/(3._dp*x*x)
285 0 : df_2_dx = -(tau*tau - 2._dp)*(E2 + x*dE2_dx)/(3._dp*tau)
286 :
287 0 : dlnT_dlnT_int = 0.75_dp*(tau + 2._dp/3._dp)*T4_int/T4
288 0 : dlnT_dlntau = 0.75_dp*(T4_int + T4_eq*(df_1_dtau + df_2_dtau))*tau/(4._dp*T4)
289 0 : dlnT_dlnx = 0.75_dp*T4_eq*(df_1_dx + df_2_dx)*x/(4._dp*T4)
290 :
291 0 : dlnT_int_dlnL = 0.25_dp
292 0 : dlnT_int_dlnR = -0.5_dp
293 0 : dlnT_int_dlnM = 0._dp
294 0 : dlnT_int_dlnkap = 0._dp
295 :
296 0 : dlntau_dlnL = 0._dp
297 0 : dlntau_dlnR = 2._dp
298 0 : dlntau_dlnM = -1._dp
299 0 : dlntau_dlnkap = 1._dp
300 :
301 0 : dlnx_dlnL = 0._dp
302 0 : dlnx_dlnR = 2._dp
303 0 : dlnx_dlnM = -1._dp
304 0 : if (gamma > 0._dp) then
305 : dlnx_dlnkap = 1._dp
306 : else
307 0 : dlnx_dlnkap = 0._dp
308 : end if
309 :
310 : dlnT_dL = (dlnT_dlnT_int*dlnT_int_dlnL + &
311 : dlnT_dlntau*dlntau_dlnL + &
312 0 : dlnT_dlnx*dlnx_dlnL)/L
313 : dlnT_dlnR = dlnT_dlnT_int*dlnT_int_dlnR + &
314 : dlnT_dlntau*dlntau_dlnR + &
315 0 : dlnT_dlnx*dlnx_dlnR
316 : dlnT_dlnM = dlnT_dlnT_int*dlnT_int_dlnM + &
317 : dlnT_dlntau*dlntau_dlnM + &
318 0 : dlnT_dlnx*dlnx_dlnM
319 : dlnT_dlnkap = dlnT_dlnT_int*dlnT_int_dlnkap + &
320 : dlnT_dlntau*dlntau_dlnkap + &
321 0 : dlnT_dlnx*dlnx_dlnkap
322 :
323 : else
324 :
325 16 : dlnT_dL = 0._dp
326 16 : dlnT_dlnR = 0._dp
327 16 : dlnT_dlnM = 0._dp
328 16 : dlnT_dlnkap = 0._dp
329 :
330 : end if
331 :
332 16 : end subroutine eval_data
333 :
334 : end module atm_irradiated
|
