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_support
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10
24 : use utils_lib, only: mesa_error, is_bad
25 :
26 : implicit none
27 :
28 : private
29 : public :: get_atm_PT
30 : public :: get_atm_PT_legacy_grey_and_kap
31 : public :: get_atm_tau_base
32 : public :: get_T_tau_id
33 : public :: build_atm
34 :
35 : contains
36 :
37 44 : subroutine get_atm_PT( &
38 : s, tau_surf, L, R, M, cgrav, skip_partials, &
39 : Teff, &
40 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
41 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
42 : ierr)
43 :
44 : type (star_info), pointer :: s
45 : real(dp), intent(in) :: tau_surf, L, R, M, cgrav
46 : logical, intent(in) :: skip_partials
47 : real(dp), intent(in) :: Teff
48 : real(dp), intent(out) :: lnT_surf
49 : real(dp), intent(out) :: dlnT_dL
50 : real(dp), intent(out) :: dlnT_dlnR
51 : real(dp), intent(out) :: dlnT_dlnM
52 : real(dp), intent(out) :: dlnT_dlnkap
53 : real(dp), intent(out) :: lnP_surf
54 : real(dp), intent(out) :: dlnP_dL
55 : real(dp), intent(out) :: dlnP_dlnR
56 : real(dp), intent(out) :: dlnP_dlnM
57 : real(dp), intent(out) :: dlnP_dlnkap
58 : integer, intent(out) :: ierr
59 :
60 44 : real(dp) :: kap
61 :
62 44 : if (is_bad(tau_surf)) then
63 0 : write(*,*) 'tau_surf', tau_surf
64 0 : ierr = -1
65 0 : return
66 : call mesa_error(__FILE__,__LINE__,'bad tau_surf arg for get_atm_PT')
67 : end if
68 :
69 : ! Evaluate surface temperature and pressure by dispatching to the
70 : ! appropriate internal routine
71 :
72 88 : select case (s% atm_option)
73 :
74 : case ('T_tau')
75 :
76 : call get_T_tau( &
77 : s, tau_surf, L, R, M, cgrav, &
78 : s% atm_T_tau_relation, s% atm_T_tau_opacity, skip_partials, &
79 : Teff, kap, &
80 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
81 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
82 44 : ierr)
83 :
84 : case ('table')
85 :
86 : call get_table( &
87 : s, skip_partials, L, R, M, cgrav, &
88 : Teff, &
89 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
90 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
91 0 : ierr)
92 :
93 : case ('irradiated_grey')
94 :
95 : call get_irradiated( &
96 : s, s% atm_irradiated_opacity, skip_partials, L, R, M, cgrav, &
97 : Teff, &
98 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
99 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
100 0 : ierr)
101 :
102 : case default
103 :
104 44 : call get_legacy(s, ierr)
105 :
106 : end select
107 :
108 : end subroutine get_atm_PT
109 :
110 :
111 0 : subroutine get_atm_PT_legacy_grey_and_kap( &
112 : s, tau_surf, L, R, M, cgrav, skip_partials, &
113 : Teff, kap, &
114 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
115 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
116 : ierr)
117 :
118 : type(star_info), pointer :: s
119 : real(dp), intent(in) :: tau_surf, L, R, M, cgrav
120 : logical, intent(in) :: skip_partials
121 : real(dp), intent(in) :: Teff
122 : real(dp), intent(out) :: kap
123 : real(dp), intent(out) :: lnT_surf
124 : real(dp), intent(out) :: dlnT_dL
125 : real(dp), intent(out) :: dlnT_dlnR
126 : real(dp), intent(out) :: dlnT_dlnM
127 : real(dp), intent(out) :: dlnT_dlnkap
128 : real(dp), intent(out) :: lnP_surf
129 : real(dp), intent(out) :: dlnP_dL
130 : real(dp), intent(out) :: dlnP_dlnR
131 : real(dp), intent(out) :: dlnP_dlnM
132 : real(dp), intent(out) :: dlnP_dlnkap
133 : integer, intent(out) :: ierr
134 :
135 : ! This routine is for legacy callers who require the old
136 : ! grey_and_kap behavior
137 :
138 : call get_T_tau( &
139 : s, tau_surf, L, R, M, cgrav, 'Eddington', 'iterated', skip_partials, &
140 : Teff, kap, &
141 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
142 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
143 0 : ierr)
144 :
145 0 : return
146 :
147 : end subroutine get_atm_PT_legacy_grey_and_kap
148 :
149 :
150 1 : subroutine get_atm_tau_base(s, tau_base, ierr)
151 :
152 : use atm_lib, only: &
153 : atm_get_T_tau_base, &
154 : atm_get_table_base
155 :
156 : type(star_info), pointer :: s
157 : real(dp), intent(out) :: tau_base
158 : integer, intent(out) :: ierr
159 :
160 : integer :: T_tau_id
161 : integer :: table_id
162 :
163 1 : ierr = 0
164 :
165 : ! Get the base optical depth
166 :
167 2 : select case (s% atm_option)
168 :
169 : case ('T_tau')
170 :
171 1 : call get_T_tau_id(s% atm_T_tau_relation, T_tau_id, ierr)
172 1 : if (ierr /= 0) then
173 0 : s% retry_message = 'Call to get_T_tau_id failed in get_atm_tau_base'
174 0 : if (s% report_ierr) write(*, *) s% retry_message
175 0 : return
176 : end if
177 :
178 1 : call atm_get_T_tau_base(T_tau_id, tau_base, ierr)
179 1 : if (ierr /= 0) then
180 0 : s% retry_message = 'Call to atm_get_T_tau_base failed in get_atm_tau_base'
181 0 : if (s% report_ierr) write(*, *) s% retry_message
182 0 : return
183 : end if
184 :
185 : case ('table')
186 :
187 0 : call get_table_id(s% atm_table, table_id, ierr)
188 0 : if (ierr /= 0) then
189 0 : s% retry_message = 'Call to get_table_id failed in get_atm_tau_base'
190 0 : if (s% report_ierr) write(*, *) s% retry_message
191 0 : return
192 : end if
193 :
194 0 : call atm_get_table_base(table_id, tau_base, ierr)
195 0 : if (ierr /= 0) then
196 0 : s% retry_message = 'Call to atm_get_table_base failed in get_atm_tau_base'
197 0 : if (s% report_ierr) write(*, *) s% retry_message
198 0 : return
199 : end if
200 :
201 : case ('irradiated_grey')
202 :
203 0 : tau_base = 2._dp/3._dp
204 :
205 : case default
206 :
207 : ! All other options use this
208 :
209 2 : tau_base = 2._dp/3._dp
210 :
211 : end select
212 :
213 : return
214 :
215 1 : end subroutine get_atm_tau_base
216 :
217 :
218 44 : subroutine get_T_tau( &
219 : s, tau_surf, L, R, M, cgrav, T_tau_relation, T_tau_opacity, skip_partials, &
220 : Teff, kap, &
221 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
222 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
223 : ierr)
224 :
225 1 : use atm_lib, only: &
226 : atm_eval_T_tau_uniform, &
227 : atm_eval_T_tau_varying
228 : use kap_support, only : prepare_kap
229 :
230 : type (star_info), pointer :: s
231 : real(dp), intent(in) :: tau_surf, L, R, M, cgrav
232 : character(*), intent(in) :: T_tau_relation
233 : character(*), intent(in) :: T_tau_opacity
234 : logical, intent(in) :: skip_partials
235 : real(dp), intent(in) :: Teff
236 : real(dp), intent(out) :: kap
237 : real(dp), intent(out) :: lnT_surf
238 : real(dp), intent(out) :: dlnT_dL
239 : real(dp), intent(out) :: dlnT_dlnR
240 : real(dp), intent(out) :: dlnT_dlnM
241 : real(dp), intent(out) :: dlnT_dlnkap
242 : real(dp), intent(out) :: lnP_surf
243 : real(dp), intent(out) :: dlnP_dL
244 : real(dp), intent(out) :: dlnP_dlnR
245 : real(dp), intent(out) :: dlnP_dlnM
246 : real(dp), intent(out) :: dlnP_dlnkap
247 : integer, intent(out) :: ierr
248 :
249 : integer :: T_tau_id
250 44 : real(dp) :: kap_guess
251 :
252 : include 'formats'
253 :
254 : ! Sanity check on L
255 :
256 44 : if (L < 0._dp) then
257 0 : s% retry_message = 'get_T_tau -- L <= 0'
258 0 : if (s% report_ierr) then
259 0 : write(*,2) 'get_T_tau: L <= 0', s% model_number, L
260 : !call mesa_error(__FILE__,__LINE__)
261 : end if
262 0 : ierr = -1
263 0 : return
264 : end if
265 :
266 : ! Get the T-tau id
267 :
268 44 : call get_T_tau_id(T_tau_relation, T_tau_id, ierr)
269 44 : if (ierr /= 0) then
270 0 : s% retry_message = 'Call to get_T_tau_id failed in get_T_tau'
271 0 : if (s% report_ierr) write(*, *) s% retry_message
272 0 : return
273 : end if
274 :
275 : ! Evaluate temperature and pressure by dispatching to the
276 : ! appropriate atm_lib routine, with the supplied T-tau relation
277 :
278 44 : select case (T_tau_opacity)
279 :
280 : case ('fixed')
281 :
282 : ! ok to use s% opacity(1) for fixed
283 : call atm_eval_T_tau_uniform( &
284 : tau_surf, L, R, M, cgrav, s% opacity(1), s% Pextra_factor, &
285 : T_tau_id, eos_proc_for_get_T_tau, kap_proc_for_get_T_tau, &
286 : s%atm_T_tau_errtol, 0, skip_partials, &
287 : Teff, kap, &
288 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
289 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
290 44 : ierr)
291 :
292 : case ('iterated')
293 :
294 0 : call prepare_kap(s, ierr)
295 0 : if (ierr /= 0) then
296 0 : s% retry_message = 'Call to prepare_kap failed in get_T_tau'
297 0 : if (s% report_ierr) write(*, *) s% retry_message
298 0 : return
299 : end if
300 :
301 : ! need to start iterations from same kap each time, so use opacity_start
302 0 : if (s% solver_iter > 0) then
303 0 : kap_guess = s% opacity_start(1)
304 : else
305 0 : kap_guess = s% opacity(1)
306 : end if
307 : call atm_eval_T_tau_uniform( &
308 : tau_surf, L, R, M, cgrav, kap_guess, s% Pextra_factor, &
309 : T_tau_id, eos_proc_for_get_T_tau, kap_proc_for_get_T_tau, &
310 : s%atm_T_tau_errtol, s%atm_T_tau_max_iters, skip_partials, &
311 : Teff, kap, &
312 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
313 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
314 0 : ierr)
315 :
316 : case ('varying')
317 :
318 0 : call prepare_kap(s, ierr)
319 0 : if (ierr /= 0) then
320 0 : s% retry_message = 'Call to prepare_kap failed in get_T_tau'
321 0 : if (s% report_ierr) write(*, *) s% retry_message
322 0 : return
323 : end if
324 :
325 : call atm_eval_T_tau_varying( &
326 : tau_surf, L, R, M, cgrav, &
327 : T_tau_id, eos_proc_for_get_T_tau, kap_proc_for_get_T_tau, &
328 : s% atm_T_tau_errtol, s% atm_T_tau_max_steps, skip_partials, &
329 : Teff, &
330 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
331 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
332 0 : ierr)
333 :
334 0 : kap = 0._dp ! This value is not used
335 :
336 : case default
337 :
338 0 : write(*,*) 'Unknown value for atm_T_tau_opacity: ' // trim(T_tau_opacity)
339 44 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
340 :
341 : end select
342 :
343 44 : return
344 :
345 : contains
346 :
347 0 : subroutine eos_proc_for_get_T_tau( &
348 : lnP, lnT, &
349 0 : lnRho, res, dres_dlnRho, dres_dlnT, &
350 : ierr)
351 :
352 : real(dp), intent(in) :: lnP
353 : real(dp), intent(in) :: lnT
354 : real(dp), intent(out) :: lnRho
355 : real(dp), intent(out) :: res(:)
356 : real(dp), intent(out) :: dres_dlnRho(:)
357 : real(dp), intent(out) :: dres_dlnT(:)
358 : integer, intent(out) :: ierr
359 : include 'formats'
360 :
361 : call eos_proc( &
362 : s, lnP, lnT, &
363 : lnRho, res, dres_dlnRho, dres_dlnT, &
364 0 : ierr)
365 :
366 44 : end subroutine eos_proc_for_get_T_tau
367 :
368 :
369 0 : subroutine kap_proc_for_get_T_tau( &
370 0 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
371 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
372 : ierr)
373 :
374 : real(dp), intent(in) :: lnRho
375 : real(dp), intent(in) :: lnT
376 : real(dp), intent(in) :: res(:)
377 : real(dp), intent(in) :: dres_dlnRho(:)
378 : real(dp), intent(in) :: dres_dlnT(:)
379 : real(dp), intent(out) :: kap
380 : real(dp), intent(out) :: dlnkap_dlnRho
381 : real(dp), intent(out) :: dlnkap_dlnT
382 : integer, intent(out) :: ierr
383 : include 'formats'
384 :
385 : call kap_proc( &
386 : s, lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
387 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
388 0 : ierr)
389 :
390 0 : end subroutine kap_proc_for_get_T_tau
391 :
392 : end subroutine get_T_tau
393 :
394 :
395 45 : subroutine get_T_tau_id (T_tau_relation, T_tau_id, ierr)
396 :
397 : use atm_def, only: &
398 : ATM_T_TAU_EDDINGTON, &
399 : ATM_T_TAU_SOLAR_HOPF, &
400 : ATM_T_TAU_KRISHNA_SWAMY, &
401 : ATM_T_TAU_TRAMPEDACH_SOLAR
402 :
403 : character(*), intent(in) :: T_tau_relation
404 : integer, intent(out) :: T_tau_id
405 : integer, intent(out) :: ierr
406 :
407 45 : ierr = 0
408 :
409 : ! Get the T-tau id
410 :
411 45 : select case (T_tau_relation)
412 : case ('Eddington')
413 45 : T_tau_id = ATM_T_TAU_EDDINGTON
414 : case ('solar_Hopf')
415 0 : T_tau_id = ATM_T_TAU_SOLAR_HOPF
416 : case ('Krishna_Swamy')
417 0 : T_tau_id = ATM_T_TAU_KRISHNA_SWAMY
418 : case ('Trampedach_solar')
419 0 : T_tau_id = ATM_T_TAU_TRAMPEDACH_SOLAR
420 : case default
421 0 : write(*,*) 'Unknown value for atm_T_tau_relation: ' // trim(T_tau_relation)
422 45 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
423 : end select
424 :
425 45 : return
426 :
427 : end subroutine get_T_tau_id
428 :
429 :
430 0 : subroutine get_table( &
431 : s, skip_partials, L, R, M, cgrav, &
432 : Teff, &
433 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
434 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
435 : ierr)
436 :
437 : use atm_lib, only: &
438 : atm_eval_table, &
439 : atm_get_table_alfa_beta, &
440 : atm_get_T_tau_base, &
441 : atm_get_table_base
442 :
443 : type (star_info), pointer :: s
444 : logical, intent(in) :: skip_partials
445 : real(dp), intent(in) :: L, R, M, cgrav
446 : real(dp), intent(in) :: Teff
447 : real(dp), intent(out) :: lnT
448 : real(dp), intent(out) :: dlnT_dL
449 : real(dp), intent(out) :: dlnT_dlnR
450 : real(dp), intent(out) :: dlnT_dlnM
451 : real(dp), intent(out) :: dlnT_dlnkap
452 : real(dp), intent(out) :: lnP
453 : real(dp), intent(out) :: dlnP_dL
454 : real(dp), intent(out) :: dlnP_dlnR
455 : real(dp), intent(out) :: dlnP_dlnM
456 : real(dp), intent(out) :: dlnP_dlnkap
457 : integer, intent(out) :: ierr
458 :
459 0 : real(dp) :: Z
460 0 : real(dp) :: tau_base
461 : integer :: T_tau_id
462 : integer :: table_id
463 0 : real(dp) :: alfa
464 0 : real(dp) :: beta
465 0 : real(dp) :: lnT_a
466 0 : real(dp) :: dlnT_dL_a
467 0 : real(dp) :: dlnT_dlnR_a
468 0 : real(dp) :: dlnT_dlnM_a
469 0 : real(dp) :: dlnT_dlnkap_a
470 0 : real(dp) :: lnP_a
471 0 : real(dp) :: dlnP_dL_a
472 0 : real(dp) :: dlnP_dlnR_a
473 0 : real(dp) :: dlnP_dlnM_a
474 0 : real(dp) :: dlnP_dlnkap_a
475 0 : real(dp) :: lnT_b
476 0 : real(dp) :: dlnT_dL_b
477 0 : real(dp) :: dlnT_dlnR_b
478 0 : real(dp) :: dlnT_dlnM_b
479 0 : real(dp) :: dlnT_dlnkap_b
480 0 : real(dp) :: lnP_b
481 0 : real(dp) :: dlnP_dL_b
482 0 : real(dp) :: dlnP_dlnR_b
483 0 : real(dp) :: dlnP_dlnM_b
484 0 : real(dp) :: dlnP_dlnkap_b
485 0 : real(dp) :: kap_a
486 :
487 : include 'formats'
488 :
489 : ! Check that tau_factor is correct
490 :
491 0 : if (s% tau_factor /= 1._dp) then
492 0 : write(*,*) 'Cannot use atm_option == ''table'' with tau_factor /= 1'
493 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
494 : end if
495 :
496 0 : Z = s% Z(1)
497 :
498 : ! Sanity check on L
499 :
500 0 : if (L < 0._dp) then
501 0 : s% retry_message = 'atm get_table: L < 0'
502 0 : if (s% report_ierr) then
503 0 : write(*,2) 'atm get_table: L < 0', s% model_number, L
504 : end if
505 0 : ierr = -1
506 0 : return
507 : end if
508 :
509 : ! Get the table id
510 :
511 0 : call get_table_id(s% atm_table, table_id, ierr)
512 0 : if (ierr /= 0) then
513 0 : s% retry_message = 'get_table_id failed in get_table'
514 0 : if (s% report_ierr) write(*, *) s% retry_message
515 0 : return
516 : end if
517 :
518 : ! Set up alfa and beta for doing table blending with off-table
519 : ! option
520 :
521 : call atm_get_table_alfa_beta( &
522 0 : L, Teff, R, M, cgrav, table_id, alfa, beta, ierr)
523 0 : if (ierr /= 0) then
524 0 : s% retry_message = 'Call to atm_get_table_alfa_beta failed in get_table'
525 0 : if (s% report_ierr) write(*, *) s% retry_message
526 0 : return
527 : end if
528 :
529 : ! If completely off the table, may need to reset tau_base to the
530 : ! T_Tau value to get the expected off-table behavior.
531 0 : if(beta == 0._dp) then
532 0 : call get_T_tau_id(s% atm_T_tau_relation, T_tau_id, ierr)
533 0 : if (ierr /= 0) then
534 0 : if (s% report_ierr) then
535 0 : write(*,*) 'Call to get_T_tau_id failed in get_table'
536 : end if
537 0 : return
538 : end if
539 :
540 0 : call atm_get_T_tau_base(T_tau_id, tau_base, ierr)
541 0 : if (ierr /= 0) then
542 0 : if (s% report_ierr) then
543 0 : write(*,*) 'Call to atm_get_T_tau_base failed in get_table'
544 : end if
545 0 : return
546 : end if
547 :
548 0 : if(s% tau_base /= tau_base) then
549 0 : write(*,*) "outside range for atm_table,"
550 0 : write(*,*) "setting T_tau value for tau_base =", tau_base
551 0 : s% tau_base = tau_base
552 : end if
553 : else ! check to see if need to switch back to table value for tau_base
554 0 : call atm_get_table_base(table_id, tau_base, ierr)
555 0 : if (ierr /= 0) then
556 0 : if (s% report_ierr) then
557 0 : write(*,*) 'Call to atm_get_table_base failed in get_table'
558 : end if
559 0 : return
560 : end if
561 :
562 0 : if(s% tau_base /= tau_base) then
563 0 : write(*,*) "back on atm_table,"
564 0 : write(*,*) "resetting to atm_table tau_base =", tau_base
565 0 : s% tau_base = tau_base
566 : end if
567 : end if
568 :
569 : ! Evaluate temperature and pressure from the table
570 :
571 0 : if (beta /= 0._dp) then
572 :
573 : call atm_eval_table( &
574 : L, R, M, cgrav, table_id, Z, skip_partials, &
575 : Teff, &
576 : lnT_b, dlnT_dL_b, dlnT_dlnR_b, dlnT_dlnM_b, dlnT_dlnkap_b, &
577 : lnP_b, dlnP_dL_b, dlnP_dlnR_b, dlnP_dlnM_b, dlnP_dlnkap_b, &
578 0 : ierr)
579 0 : if (ierr /= 0) then
580 0 : s% retry_message = 'Call to atm_eval_table failed in get_table'
581 0 : if (s% report_ierr) write(*, *) s% retry_message
582 0 : return
583 : end if
584 :
585 : else
586 :
587 0 : lnT_b = 0._dp
588 0 : dlnT_dL_b = 0._dp
589 0 : dlnT_dlnR_b = 0._dp
590 0 : dlnT_dlnM_b = 0._dp
591 0 : dlnT_dlnkap_b = 0._dp
592 :
593 0 : lnP_b = 0._dp
594 0 : dlnP_dL_b = 0._dp
595 0 : dlnP_dlnR_b = 0._dp
596 0 : dlnP_dlnM_b = 0._dp
597 0 : dlnP_dlnkap_b = 0._dp
598 :
599 : end if
600 :
601 : ! Evaluate temperature and pressure from the backup atmosphere
602 : ! option
603 :
604 0 : if (alfa /= 0._dp) then
605 :
606 0 : select case (s% atm_off_table_option)
607 :
608 : case ('T_tau')
609 :
610 : call get_T_tau( &
611 : s, s% tau_base, L, R, M, cgrav, &
612 : s% atm_T_tau_relation, s% atm_T_tau_opacity, skip_partials, &
613 : Teff, kap_a, &
614 : lnT_a, dlnT_dL_a, dlnT_dlnR_a, dlnT_dlnM_a, dlnT_dlnkap_a, &
615 : lnP_a, dlnP_dL_a, dlnP_dlnR_a, dlnP_dlnM_a, dlnP_dlnkap_a, &
616 0 : ierr)
617 0 : if (ierr /= 0) then
618 0 : s% retry_message = 'Call to get_T_tau failed in get_table'
619 0 : if (s% report_ierr) write(*, *) s% retry_message
620 0 : return
621 : end if
622 :
623 : case ('')
624 :
625 0 : write(*,*) 'Attempt to interpolate outside atmosphere table'
626 0 : call mesa_error(__FILE__,__LINE__,'Try setting the which_off_table_option control in your inlist file')
627 0 : stop
628 :
629 : case default
630 :
631 0 : write(*,*) 'Unknown value for atm_off_table_option: ' // trim(s% atm_off_table_option)
632 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
633 :
634 : end select
635 :
636 : else
637 :
638 0 : lnT_a = 0._dp
639 0 : dlnT_dL_a = 0._dp
640 0 : dlnT_dlnR_a = 0._dp
641 0 : dlnT_dlnM_a = 0._dp
642 0 : dlnT_dlnkap_a = 0._dp
643 :
644 0 : lnP_a = 0._dp
645 0 : dlnP_dL_a = 0._dp
646 0 : dlnP_dlnR_a = 0._dp
647 0 : dlnP_dlnM_a = 0._dp
648 0 : dlnP_dlnkap_a = 0._dp
649 :
650 : end if
651 :
652 : ! Blend the results together
653 :
654 0 : lnT = alfa*lnT_a + beta*lnT_b
655 0 : lnP = alfa*lnP_a + beta*lnP_b
656 :
657 0 : if (.not. skip_partials) then
658 :
659 0 : dlnT_dL = alfa*dlnT_dL_a + beta*dlnT_dL_b
660 0 : dlnT_dlnR = alfa*dlnT_dlnR_a + beta*dlnT_dlnR_b
661 0 : dlnT_dlnM = alfa*dlnT_dlnM_a + beta*dlnT_dlnM_b
662 0 : dlnT_dlnkap = alfa*dlnT_dlnkap_a + beta*dlnT_dlnkap_b
663 :
664 0 : dlnP_dL = alfa*dlnP_dL_a + beta*dlnP_dL_b
665 0 : dlnP_dlnR = alfa*dlnP_dlnR_a + beta*dlnP_dlnR_b
666 0 : dlnP_dlnM = alfa*dlnP_dlnM_a + beta*dlnP_dlnM_b
667 0 : dlnP_dlnkap = alfa*dlnP_dlnkap_a + beta*dlnP_dlnkap_b
668 :
669 : end if
670 :
671 : ! Set the effective temperature
672 :
673 : ! if (alfa /= 0._dp) then
674 : ! if (beta /= 0._dp) then
675 : ! if (Teff_a /= Teff_b) then
676 : ! ierr = -1
677 : ! s% retry_message = 'Mismatch between Teff values in get_tables'
678 : ! write(*,*) 'Mismatch between Teff values in get_tables: ', Teff_a, Teff_b
679 : ! !call mesa_error(__FILE__,__LINE__)
680 : ! end if
681 : ! end if
682 : ! Teff = Teff_a
683 : ! else
684 : ! Teff = Teff_b
685 : ! end if
686 :
687 : return
688 :
689 0 : end subroutine get_table
690 :
691 :
692 0 : subroutine get_table_id (table_name, table_id, ierr)
693 :
694 0 : use atm_def, only: &
695 : ATM_TABLE_TAU_100, &
696 : ATM_TABLE_TAU_10, &
697 : ATM_TABLE_TAU_1, &
698 : ATM_TABLE_TAU_1M1, &
699 : ATM_TABLE_PHOTOSPHERE, &
700 : ATM_TABLE_WD_TAU_25, &
701 : ATM_TABLE_DB_WD_TAU_25
702 :
703 : character(*), intent(in) :: table_name
704 : integer, intent(out) :: table_id
705 : integer, intent(out) :: ierr
706 :
707 0 : ierr = 0
708 :
709 : ! Get the table id
710 :
711 0 : select case (table_name)
712 : case ('tau_100')
713 0 : table_id = ATM_TABLE_TAU_100
714 : case ('tau_10')
715 0 : table_id = ATM_TABLE_TAU_10
716 : case ('tau_1')
717 0 : table_id = ATM_TABLE_TAU_1
718 : case ('tau_1m1')
719 0 : table_id = ATM_TABLE_TAU_1M1
720 : case ('photosphere')
721 0 : table_id = ATM_TABLE_PHOTOSPHERE
722 : case ('WD_tau_25')
723 0 : table_id = ATM_TABLE_WD_TAU_25
724 : case ('DB_WD_tau_25')
725 0 : table_id = ATM_TABLE_DB_WD_TAU_25
726 : case default
727 0 : write(*,*) 'Unknown value for atm_table: ' // trim(table_name)
728 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
729 : end select
730 :
731 0 : return
732 :
733 : end subroutine get_table_id
734 :
735 :
736 0 : subroutine get_irradiated( &
737 : s, irradiated_opacity, skip_partials, L, R, M, cgrav, &
738 : Teff, &
739 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
740 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
741 : ierr)
742 :
743 : use atm_lib, only: atm_eval_irradiated
744 : use kap_support, only: prepare_kap
745 :
746 : type (star_info), pointer :: s
747 : character(*), intent(in) :: irradiated_opacity
748 : logical, intent(in) :: skip_partials
749 : real(dp), intent(in) :: L, R, M, cgrav
750 : real(dp), intent(in) :: Teff
751 : real(dp), intent(out) :: lnT_surf
752 : real(dp), intent(out) :: dlnT_dL
753 : real(dp), intent(out) :: dlnT_dlnR
754 : real(dp), intent(out) :: dlnT_dlnM
755 : real(dp), intent(out) :: dlnT_dlnkap
756 : real(dp), intent(out) :: lnP_surf
757 : real(dp), intent(out) :: dlnP_dL
758 : real(dp), intent(out) :: dlnP_dlnR
759 : real(dp), intent(out) :: dlnP_dlnM
760 : real(dp), intent(out) :: dlnP_dlnkap
761 : integer, intent(out) :: ierr
762 :
763 0 : real(dp) :: kap_for_atm
764 0 : real(dp) :: kap
765 0 : real(dp) :: tau_surf
766 :
767 : include 'formats'
768 :
769 0 : if (s% solver_iter > 0) then
770 0 : kap_for_atm = s% opacity_start(1)
771 : else
772 0 : kap_for_atm = s% opacity(1)
773 : end if
774 :
775 : ! Sanity check on L
776 :
777 0 : if (L < 0._dp) then
778 0 : s% retry_message = 'get_irradiated: L <= 0'
779 0 : if (s% report_ierr) then
780 0 : write(*,2) 'get_irradiated: L <= 0', s% model_number, L
781 : !call mesa_error(__FILE__,__LINE__)
782 : end if
783 0 : ierr = -1
784 0 : return
785 : end if
786 :
787 : ! Evaluate temperature and pressure by dispatching to the
788 : ! appropriate atm_lib routine
789 :
790 0 : select case (irradiated_opacity)
791 :
792 : case ('fixed')
793 :
794 : call atm_eval_irradiated( &
795 : L, R, M, cgrav, s% atm_irradiated_T_eq, s% atm_irradiated_P_surf, &
796 : kap_for_atm, s% atm_irradiated_kap_v, s% atm_irradiated_kap_v_div_kap_th, &
797 : eos_proc_for_get_irradiated, kap_proc_for_get_irradiated, &
798 : s% atm_irradiated_errtol, 0, skip_partials, &
799 : Teff, kap, tau_surf, &
800 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
801 0 : ierr)
802 :
803 : case ('iterated')
804 :
805 0 : call prepare_kap(s, ierr)
806 0 : if (ierr /= 0) then
807 0 : s% retry_message = 'Failed in call to prepare_kap'
808 0 : if (s% report_ierr) write(*, *) s% retry_message
809 0 : return
810 : end if
811 :
812 : call atm_eval_irradiated( &
813 : L, R, M, cgrav, s% atm_irradiated_T_eq, s% atm_irradiated_P_surf, &
814 : kap_for_atm, s% atm_irradiated_kap_v, s% atm_irradiated_kap_v_div_kap_th, &
815 : eos_proc_for_get_irradiated, kap_proc_for_get_irradiated, &
816 : s% atm_irradiated_errtol, s% atm_irradiated_max_iters, skip_partials, &
817 : Teff, kap, tau_surf, &
818 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
819 0 : ierr)
820 :
821 : case default
822 :
823 0 : write(*,*) 'Unknown value for atm_irradiated_opacity: ' // trim(irradiated_opacity)
824 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
825 :
826 : end select
827 :
828 : ! Set up remaining values
829 :
830 0 : lnP_surf = log(s% atm_irradiated_P_surf)
831 :
832 0 : dlnP_dL = 0._dp
833 0 : dlnP_dlnR = 0._dp
834 0 : dlnP_dlnM = 0._dp
835 0 : dlnP_dlnkap = 0._dp
836 :
837 : ! Update tau_factor
838 :
839 0 : s% tau_factor = tau_surf/s% tau_base
840 :
841 0 : return
842 :
843 : contains
844 :
845 0 : subroutine eos_proc_for_get_irradiated( &
846 : lnP, lnT, &
847 0 : lnRho, res, dres_dlnRho, dres_dlnT, &
848 : ierr)
849 :
850 : real(dp), intent(in) :: lnP
851 : real(dp), intent(in) :: lnT
852 : real(dp), intent(out) :: lnRho
853 : real(dp), intent(out) :: res(:)
854 : real(dp), intent(out) :: dres_dlnRho(:)
855 : real(dp), intent(out) :: dres_dlnT(:)
856 : integer, intent(out) :: ierr
857 :
858 : call eos_proc( &
859 : s, lnP, lnT, &
860 : lnRho, res, dres_dlnRho, dres_dlnT, &
861 0 : ierr)
862 :
863 0 : end subroutine eos_proc_for_get_irradiated
864 :
865 :
866 0 : subroutine kap_proc_for_get_irradiated( &
867 0 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
868 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
869 : ierr)
870 :
871 : real(dp), intent(in) :: lnRho
872 : real(dp), intent(in) :: lnT
873 : real(dp), intent(in) :: res(:)
874 : real(dp), intent(in) :: dres_dlnRho(:)
875 : real(dp), intent(in) :: dres_dlnT(:)
876 : real(dp), intent(out) :: kap
877 : real(dp), intent(out) :: dlnkap_dlnRho
878 : real(dp), intent(out) :: dlnkap_dlnT
879 : integer, intent(out) :: ierr
880 :
881 : call kap_proc( &
882 : s, lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
883 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
884 0 : ierr)
885 :
886 0 : end subroutine kap_proc_for_get_irradiated
887 :
888 : end subroutine get_irradiated
889 :
890 :
891 0 : subroutine get_legacy (s, ierr)
892 :
893 : type(star_info), pointer :: s
894 : integer, intent(out) :: ierr
895 :
896 : ! Handles legacy / invalid choices for atm_option. This is
897 : ! to guide users toward the newer atmosphere controls; in the long
898 : ! term, it can be removed
899 :
900 0 : select case (s% atm_option)
901 :
902 : case ('simple_photosphere')
903 :
904 0 : write(*,*) 'The ''simple_photosphere'' choice for atm_option is deprecated'
905 0 : write(*,*) 'To achieve the same results, set:'
906 0 : write(*,*) ' atm_option = ''T_tau'''
907 0 : write(*,*) ' atm_T_tau_relation = ''Eddington'''
908 0 : write(*,*) ' atm_T_tau_opacity = ''fixed'''
909 :
910 : case ('grey_and_kap')
911 :
912 0 : write(*,*) 'The ''grey_and_kap'' choice for atm_option is deprecated'
913 0 : write(*,*) 'To achieve the same results, set:'
914 0 : write(*,*) ' atm_option = ''T_tau'''
915 0 : write(*,*) ' atm_T_tau_relation = ''Eddington'''
916 0 : write(*,*) ' atm_T_tau_opacity = ''iterated'''
917 :
918 : case ('Eddington_grey')
919 :
920 0 : write(*,*) 'The ''Eddington_grey'' choice for atm_option is deprecated'
921 0 : write(*,*) 'To achieve the same results, set:'
922 0 : write(*,*) ' atm_option = ''T_tau'''
923 0 : write(*,*) ' atm_T_tau_relation = ''Eddington'''
924 0 : write(*,*) ' atm_T_tau_opacity = ''varying'''
925 :
926 : case ('solar_Hopf_grey')
927 :
928 0 : write(*,*) 'The ''solar_Hopf_grey'' choice for atm_option is deprecated'
929 0 : write(*,*) 'To achieve the same results, set:'
930 0 : write(*,*) ' atm_option = ''T_tau'''
931 0 : write(*,*) ' atm_T_tau_relation = ''solar_Hopf'''
932 0 : write(*,*) ' atm_T_tau_opacity = ''varying'''
933 :
934 : case ('Krishna_Swamy')
935 :
936 0 : write(*,*) 'The ''Krishna_Swamy'' choice for atm_option is deprecated'
937 0 : write(*,*) 'To achieve the same results, set:'
938 0 : write(*,*) ' atm_option = ''T_tau'''
939 0 : write(*,*) ' atm_T_tau_relation = ''Krishna_Swamy'''
940 0 : write(*,*) ' atm_T_tau_opacity = ''varying'''
941 :
942 : case ('tau_100_tables')
943 :
944 0 : write(*,*) 'The ''tau_100_tables'' choice for atm_option is deprecated'
945 0 : write(*,*) 'To achieve the same results, set:'
946 0 : write(*,*) ' atm_option = ''table'''
947 0 : write(*,*) ' atm_table = ''tau_100'''
948 :
949 : case ('tau_10_tables')
950 :
951 0 : write(*,*) 'The ''tau_10_tables'' choice for atm_option is deprecated'
952 0 : write(*,*) 'To achieve the same results, set:'
953 0 : write(*,*) ' atm_option = ''table'''
954 0 : write(*,*) ' atm_table = ''tau_10'''
955 :
956 : case ('tau_1_tables')
957 :
958 0 : write(*,*) 'The ''tau_1_tables'' choice for atm_option is deprecated'
959 0 : write(*,*) 'To achieve the same results, set:'
960 0 : write(*,*) ' atm_option = ''table'''
961 0 : write(*,*) ' atm_table = ''tau_1'''
962 :
963 : case ('tau_1m1_tables')
964 :
965 0 : write(*,*) 'The ''tau_1m1_tables'' choice for atm_option is deprecated'
966 0 : write(*,*) 'To achieve the same results, set:'
967 0 : write(*,*) ' atm_option = ''table'''
968 0 : write(*,*) ' atm_table = ''tau_1m1'''
969 :
970 : case ('photosphere_tables')
971 :
972 0 : write(*,*) 'The ''tau_10_tables'' choice for atm_option is deprecated'
973 0 : write(*,*) 'To achieve the same results, set:'
974 0 : write(*,*) ' atm_option = ''table'''
975 0 : write(*,*) ' atm_table = ''photosphere'''
976 :
977 : case ('grey_irradiated')
978 :
979 0 : write(*,*) 'The ''grey_irradiated'' choice for atm_option is deprecated'
980 0 : write(*,*) 'To achieve the same results, set:'
981 0 : write(*,*) ' atm_option = ''irradiated_grey'''
982 0 : write(*,*) ' atm_irradiated_opacity = ''fixed'' (if atm_grey_irradiated_simple_kap_th = .true.)'
983 0 : write(*,*) ' atm_irradiated_opacity = ''varying'' (if atm_grey_irradiated_simple_kap_th = .false.)'
984 :
985 : case default
986 :
987 0 : write(*,*) 'Unknown value for atm_option: ' // trim(s% atm_option)
988 :
989 : end select
990 :
991 : ! Stop because we can't continue
992 :
993 0 : ierr = -1 ! ifort complains if this isn't set
994 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
995 :
996 0 : return
997 :
998 : end subroutine get_legacy
999 :
1000 :
1001 0 : subroutine build_atm( &
1002 : s, L, R, Teff, M, cgrav, ierr)
1003 :
1004 : type(star_info), pointer :: s
1005 : real(dp), intent(in) :: L, R, Teff, M, cgrav
1006 : integer, intent(out) :: ierr
1007 :
1008 : ! Create the atmosphere structure by dispatching to the
1009 : ! appropriate internal routine
1010 :
1011 0 : select case (s% atm_option)
1012 :
1013 : case ('T_tau')
1014 :
1015 : call build_T_tau( &
1016 : s, s% tau_factor*s% tau_base, L, R, Teff, M, cgrav, &
1017 : s% atm_T_tau_relation, s% atm_T_tau_opacity, &
1018 : s% atm_structure_num_pts, s% atm_structure, &
1019 0 : ierr)
1020 :
1021 : case default
1022 :
1023 0 : write(*,*) 'Cannot create atm structure for atm_option: ' // TRIM(s% atm_option)
1024 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
1025 :
1026 : end select
1027 :
1028 0 : end subroutine build_atm
1029 :
1030 :
1031 0 : subroutine build_T_tau( &
1032 : s, tau_surf, L, R, Teff, M, cgrav, T_tau_relation, T_tau_opacity, &
1033 : atm_structure_num_pts, atm_structure, &
1034 : ierr)
1035 :
1036 : use atm_lib, only: &
1037 : atm_build_T_tau_uniform, &
1038 : atm_build_T_tau_varying
1039 :
1040 : type(star_info), pointer :: s
1041 : real(dp), intent(in) :: tau_surf, L, R, Teff, M, cgrav
1042 : character(*), intent(in) :: T_tau_relation
1043 : character(*), intent(in) :: T_tau_opacity
1044 : integer, intent(out) :: atm_structure_num_pts
1045 : real(dp), pointer :: atm_structure(:,:)
1046 : integer, intent(out) :: ierr
1047 :
1048 : real(dp) :: kap
1049 : real(dp) :: lnT_surf
1050 : real(dp) :: dlnT_dL
1051 : real(dp) :: dlnT_dlnR
1052 : real(dp) :: dlnT_dlnM
1053 : real(dp) :: dlnT_dlnkap
1054 : real(dp) :: lnP_surf
1055 : real(dp) :: dlnP_dL
1056 : real(dp) :: dlnP_dlnR
1057 : real(dp) :: dlnP_dlnM
1058 : real(dp) :: dlnP_dlnkap
1059 : integer :: T_tau_id
1060 :
1061 : ierr = 0
1062 :
1063 : ! First evaluate the temperature, pressure and opacity
1064 :
1065 : call get_T_tau( &
1066 : s, tau_surf, L, R, M, cgrav, T_tau_relation, T_tau_opacity, .TRUE., &
1067 : Teff, kap, &
1068 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
1069 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
1070 0 : ierr)
1071 0 : if (ierr /= 0) then
1072 0 : s% retry_message = 'Call to get_T_tau failed in build_T_tau'
1073 0 : if (s% report_ierr) write(*, *) s% retry_message
1074 0 : return
1075 : end if
1076 :
1077 : ! Get the T-tau id
1078 :
1079 0 : call get_T_tau_id(T_tau_relation, T_tau_id, ierr)
1080 0 : if (ierr /= 0) then
1081 0 : s% retry_message = 'Call to get_T_tau_id failed in build_T_tau'
1082 0 : if (s% report_ierr) write(*, *) s% retry_message
1083 0 : return
1084 : end if
1085 :
1086 : ! Build the atmosphere structure by dispatching to the
1087 : ! appropriate atm_lib routine, with the supplied T-tau relation
1088 :
1089 0 : select case (T_tau_opacity)
1090 :
1091 : case ('fixed', 'iterated')
1092 :
1093 : call atm_build_T_tau_uniform( &
1094 : tau_surf, L, R, Teff, M, cgrav, kap, s% Pextra_factor, s% atm_build_tau_outer, &
1095 : T_tau_id, eos_proc_for_build_T_tau, kap_proc_for_build_T_tau, &
1096 : s% atm_build_errtol, s% atm_build_dlogtau, &
1097 : atm_structure_num_pts, atm_structure, &
1098 0 : ierr)
1099 0 : if (ierr /= 0) then
1100 0 : s% retry_message = 'Call to atm_build_T_tau_uniform failed in build_T_tau'
1101 0 : if (s% report_ierr) write(*, *) s% retry_message
1102 0 : return
1103 : end if
1104 :
1105 : case ('varying')
1106 :
1107 : call atm_build_T_tau_varying( &
1108 : tau_surf, L, R, Teff, M, cgrav, lnP_surf, s% atm_build_tau_outer, &
1109 : T_tau_id, eos_proc_for_build_T_tau, kap_proc_for_build_T_tau, &
1110 : s% atm_build_errtol, s% atm_build_dlogtau, &
1111 : atm_structure_num_pts, atm_structure, &
1112 0 : ierr)
1113 0 : if (ierr /= 0) then
1114 0 : s% retry_message = 'Call to atm_build_T_tau_varying failed in build_T_tau'
1115 0 : if (s% report_ierr) write(*, *) s% retry_message
1116 0 : return
1117 : end if
1118 :
1119 : case default
1120 :
1121 0 : write(*,*) 'Unknown value for atm_T_tau_opacity: ' // trim(T_tau_opacity)
1122 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
1123 :
1124 : end select
1125 :
1126 0 : return
1127 :
1128 : contains
1129 :
1130 0 : subroutine eos_proc_for_build_T_tau( &
1131 : lnP, lnT, &
1132 0 : lnRho, res, dres_dlnRho, dres_dlnT, &
1133 : ierr)
1134 :
1135 : real(dp), intent(in) :: lnP
1136 : real(dp), intent(in) :: lnT
1137 : real(dp), intent(out) :: lnRho
1138 : real(dp), intent(out) :: res(:)
1139 : real(dp), intent(out) :: dres_dlnRho(:)
1140 : real(dp), intent(out) :: dres_dlnT(:)
1141 : integer, intent(out) :: ierr
1142 :
1143 : call eos_proc( &
1144 : s, lnP, lnT, &
1145 : lnRho, res, dres_dlnRho, dres_dlnT, &
1146 0 : ierr)
1147 :
1148 0 : end subroutine eos_proc_for_build_T_tau
1149 :
1150 :
1151 0 : subroutine kap_proc_for_build_T_tau( &
1152 0 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
1153 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
1154 : ierr)
1155 :
1156 : real(dp), intent(in) :: lnRho
1157 : real(dp), intent(in) :: lnT
1158 : real(dp), intent(in) :: res(:)
1159 : real(dp), intent(in) :: dres_dlnRho(:)
1160 : real(dp), intent(in) :: dres_dlnT(:)
1161 : real(dp), intent(out) :: kap
1162 : real(dp), intent(out) :: dlnkap_dlnRho
1163 : real(dp), intent(out) :: dlnkap_dlnT
1164 : integer, intent(out) :: ierr
1165 :
1166 : call kap_proc( &
1167 : s, lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
1168 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
1169 0 : ierr)
1170 :
1171 0 : end subroutine kap_proc_for_build_T_tau
1172 :
1173 : end subroutine build_T_tau
1174 :
1175 :
1176 0 : subroutine eos_proc( &
1177 : s, lnP, lnT, &
1178 : lnRho, res, dres_dlnRho, dres_dlnT, &
1179 : ierr)
1180 :
1181 : use eos_def, only: num_eos_basic_results
1182 : use eos_lib, only: radiation_pressure, eos_gamma_PT_get_rho_energy
1183 : use eos_support, only: solve_eos_given_PgasT
1184 :
1185 : type(star_info), pointer :: s
1186 : real(dp), intent(in) :: lnP
1187 : real(dp), intent(in) :: lnT
1188 : real(dp), intent(out) :: lnRho
1189 : real(dp), intent(out) :: res(num_eos_basic_results)
1190 : real(dp), intent(out) :: dres_dlnRho(num_eos_basic_results)
1191 : real(dp), intent(out) :: dres_dlnT(num_eos_basic_results)
1192 : integer, intent(out) :: ierr
1193 :
1194 : real(dp), parameter :: LOGRHO_TOL = 1E-11_dp
1195 : real(dp), parameter :: LOGPGAS_TOL = 1E-11_dp
1196 :
1197 0 : real(dp) :: T
1198 0 : real(dp) :: P
1199 0 : real(dp) :: Prad
1200 0 : real(dp) :: Pgas
1201 0 : real(dp) :: rho, gamma, energy
1202 0 : real(dp) :: logRho, logRho_guess
1203 0 : real(dp) :: dres_dxa(num_eos_basic_results,s% species)
1204 :
1205 0 : T = exp(lnT)
1206 0 : P = exp(lnP)
1207 :
1208 0 : Prad = radiation_pressure(T)
1209 0 : Pgas = MAX(1.E-99_dp, P - Prad)
1210 :
1211 0 : gamma = 5d0/3d0
1212 : call eos_gamma_PT_get_rho_energy( &
1213 0 : s% abar(1), P, T, gamma, rho, energy, ierr)
1214 0 : if (ierr /= 0) then
1215 0 : s% retry_message = 'Call to eos_gamma_PT_get_rho_energy failed in eos_proc'
1216 0 : if (s% report_ierr) write(*, *) trim(s% retry_message)
1217 : end if
1218 :
1219 0 : logRho_guess = log10(rho)
1220 :
1221 : call solve_eos_given_PgasT( &
1222 : s, 0, s% xa(:,1), &
1223 : lnT/ln10, log10(Pgas), logRho_guess, LOGRHO_TOL, LOGPGAS_TOL, &
1224 : logRho, res, dres_dlnRho, dres_dlnT, dres_dxa, &
1225 0 : ierr)
1226 0 : if (ierr /= 0) then
1227 0 : s% retry_message = 'Call to solve_eos_given_PgasT failed in eos_proc'
1228 0 : if (s% report_ierr) write(*, *) trim(s% retry_message)
1229 : end if
1230 :
1231 0 : lnRho = logRho*ln10
1232 :
1233 :
1234 0 : return
1235 :
1236 0 : end subroutine eos_proc
1237 :
1238 :
1239 0 : subroutine kap_proc( &
1240 0 : s, lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
1241 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
1242 : ierr)
1243 :
1244 0 : use eos_def, only: i_lnfree_e, i_eta
1245 : use kap_def, only: num_kap_fracs
1246 : use kap_support, only: get_kap
1247 :
1248 : type(star_info), pointer :: s
1249 : real(dp), intent(in) :: lnRho
1250 : real(dp), intent(in) :: lnT
1251 : real(dp), intent(in) :: res(:)
1252 : real(dp), intent(in) :: dres_dlnRho(:)
1253 : real(dp), intent(in) :: dres_dlnT(:)
1254 : real(dp), intent(out) :: kap
1255 : real(dp), intent(out) :: dlnkap_dlnRho
1256 : real(dp), intent(out) :: dlnkap_dlnT
1257 : integer, intent(out) :: ierr
1258 :
1259 0 : real(dp) :: kap_fracs(num_kap_fracs)
1260 :
1261 : include 'formats'
1262 :
1263 : call get_kap( &
1264 : s, 0, s% zbar(1), s% xa(:,1), lnRho/ln10, lnT/ln10, &
1265 : res(i_lnfree_e), dres_dlnRho(i_lnfree_e), dres_dlnT(i_lnfree_e), &
1266 : res(i_eta), dres_dlnRho(i_eta), dres_dlnT(i_eta), &
1267 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, &
1268 0 : ierr)
1269 0 : if (ierr /= 0) then
1270 0 : s% retry_message = 'Call to get_kap failed in kap_proc'
1271 0 : if (s% report_ierr) write(*, *) s% retry_message
1272 : end if
1273 :
1274 0 : if (kap <= 0d0 .or. is_bad(kap)) then
1275 0 : write(*,1) 'bad kap', kap
1276 0 : write(*,1) 's% zbar(1)', s% zbar(1)
1277 0 : write(*,1) 'lnRho/ln10', lnRho/ln10
1278 0 : write(*,1) 'lnT/ln10', lnT/ln10
1279 0 : write(*,1) 'res(i_eta)', res(i_eta)
1280 0 : ierr = -1
1281 0 : return
1282 : call mesa_error(__FILE__,__LINE__,'atm kap_proc')
1283 : end if
1284 :
1285 : return
1286 :
1287 0 : end subroutine kap_proc
1288 :
1289 : end module atm_support
|
