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 : 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 : end subroutine get_atm_tau_base
216 :
217 :
218 132 : 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 : 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 : 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 : 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 0 : 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_table_base
441 :
442 : type (star_info), pointer :: s
443 : logical, intent(in) :: skip_partials
444 : real(dp), intent(in) :: L, R, M, cgrav
445 : real(dp), intent(in) :: Teff
446 : real(dp), intent(out) :: lnT
447 : real(dp), intent(out) :: dlnT_dL
448 : real(dp), intent(out) :: dlnT_dlnR
449 : real(dp), intent(out) :: dlnT_dlnM
450 : real(dp), intent(out) :: dlnT_dlnkap
451 : real(dp), intent(out) :: lnP
452 : real(dp), intent(out) :: dlnP_dL
453 : real(dp), intent(out) :: dlnP_dlnR
454 : real(dp), intent(out) :: dlnP_dlnM
455 : real(dp), intent(out) :: dlnP_dlnkap
456 : integer, intent(out) :: ierr
457 :
458 : real(dp) :: Z
459 : real(dp) :: tau_base
460 : integer :: table_id
461 : real(dp) :: alfa
462 : real(dp) :: beta
463 : real(dp) :: lnT_a
464 : real(dp) :: dlnT_dL_a
465 : real(dp) :: dlnT_dlnR_a
466 : real(dp) :: dlnT_dlnM_a
467 : real(dp) :: dlnT_dlnkap_a
468 : real(dp) :: lnP_a
469 : real(dp) :: dlnP_dL_a
470 : real(dp) :: dlnP_dlnR_a
471 : real(dp) :: dlnP_dlnM_a
472 : real(dp) :: dlnP_dlnkap_a
473 : real(dp) :: lnT_b
474 : real(dp) :: dlnT_dL_b
475 : real(dp) :: dlnT_dlnR_b
476 : real(dp) :: dlnT_dlnM_b
477 : real(dp) :: dlnT_dlnkap_b
478 : real(dp) :: lnP_b
479 : real(dp) :: dlnP_dL_b
480 : real(dp) :: dlnP_dlnR_b
481 : real(dp) :: dlnP_dlnM_b
482 : real(dp) :: dlnP_dlnkap_b
483 : real(dp) :: kap_a
484 :
485 : include 'formats'
486 :
487 : ! Check that tau_factor is correct
488 :
489 0 : if (s% tau_factor /= 1._dp) then
490 0 : write(*,*) 'Cannot use atm_option == ''table'' with tau_factor /= 1'
491 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
492 : end if
493 :
494 0 : Z = s% Z(1)
495 :
496 : ! Sanity check on L
497 :
498 0 : if (L < 0._dp) then
499 0 : s% retry_message = 'atm get_table: L < 0'
500 0 : if (s% report_ierr) then
501 0 : write(*,2) 'atm get_table: L < 0', s% model_number, L
502 : end if
503 0 : ierr = -1
504 0 : return
505 : end if
506 :
507 : ! Get the table id
508 :
509 0 : call get_table_id(s% atm_table, table_id, ierr)
510 0 : if (ierr /= 0) then
511 0 : s% retry_message = 'get_table_id failed in get_table'
512 0 : if (s% report_ierr) write(*, *) s% retry_message
513 0 : return
514 : end if
515 :
516 : ! Set up alfa and beta for doing table blending with off-table
517 : ! option
518 :
519 : call atm_get_table_alfa_beta( &
520 0 : L, Teff, R, M, cgrav, table_id, alfa, beta, ierr)
521 0 : if (ierr /= 0) then
522 0 : s% retry_message = 'Call to atm_get_table_alfa_beta failed in get_table'
523 0 : if (s% report_ierr) write(*, *) s% retry_message
524 0 : return
525 : end if
526 :
527 : ! If completely off the table, use the same
528 : ! attachment depth for the table and backup atmosphere.
529 0 : call atm_get_table_base(table_id, tau_base, ierr)
530 0 : if (ierr /= 0) then
531 0 : if (s% report_ierr) then
532 0 : write(*,*) 'Call to atm_get_table_base failed in get_table'
533 : end if
534 0 : return
535 : end if
536 0 : s% tau_base = tau_base
537 0 : if (beta == 0._dp) then
538 0 : write(*,*) "outside range for atm_table,"
539 0 : write(*,*) "setting T_tau value for tau_base =", tau_base
540 : end if
541 :
542 : ! Evaluate temperature and pressure from the table
543 :
544 0 : if (beta /= 0._dp) then
545 :
546 : call atm_eval_table( &
547 : L, R, M, cgrav, table_id, Z, skip_partials, &
548 : Teff, &
549 : lnT_b, dlnT_dL_b, dlnT_dlnR_b, dlnT_dlnM_b, dlnT_dlnkap_b, &
550 : lnP_b, dlnP_dL_b, dlnP_dlnR_b, dlnP_dlnM_b, dlnP_dlnkap_b, &
551 0 : ierr)
552 0 : if (ierr /= 0) then
553 0 : s% retry_message = 'Call to atm_eval_table failed in get_table'
554 0 : if (s% report_ierr) write(*, *) s% retry_message
555 0 : return
556 : end if
557 :
558 : else
559 :
560 0 : lnT_b = 0._dp
561 0 : dlnT_dL_b = 0._dp
562 0 : dlnT_dlnR_b = 0._dp
563 0 : dlnT_dlnM_b = 0._dp
564 0 : dlnT_dlnkap_b = 0._dp
565 :
566 0 : lnP_b = 0._dp
567 0 : dlnP_dL_b = 0._dp
568 0 : dlnP_dlnR_b = 0._dp
569 0 : dlnP_dlnM_b = 0._dp
570 0 : dlnP_dlnkap_b = 0._dp
571 :
572 : end if
573 :
574 : ! Evaluate temperature and pressure from the backup atmosphere
575 : ! option
576 :
577 0 : if (alfa /= 0._dp) then
578 :
579 0 : select case (s% atm_off_table_option)
580 :
581 : case ('T_tau')
582 :
583 : call get_T_tau( &
584 : s, s% tau_base, L, R, M, cgrav, &
585 : s% atm_T_tau_relation, s% atm_T_tau_opacity, skip_partials, &
586 : Teff, kap_a, &
587 : lnT_a, dlnT_dL_a, dlnT_dlnR_a, dlnT_dlnM_a, dlnT_dlnkap_a, &
588 : lnP_a, dlnP_dL_a, dlnP_dlnR_a, dlnP_dlnM_a, dlnP_dlnkap_a, &
589 0 : ierr)
590 0 : if (ierr /= 0) then
591 0 : s% retry_message = 'Call to get_T_tau failed in get_table'
592 0 : if (s% report_ierr) write(*, *) s% retry_message
593 0 : return
594 : end if
595 :
596 : case ('')
597 :
598 0 : write(*,*) 'Attempt to interpolate outside atmosphere table'
599 0 : call mesa_error(__FILE__,__LINE__,'Try setting the which_off_table_option control in your inlist file')
600 0 : stop
601 :
602 : case default
603 :
604 0 : write(*,*) 'Unknown value for atm_off_table_option: ' // trim(s% atm_off_table_option)
605 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
606 :
607 : end select
608 :
609 : else
610 :
611 0 : lnT_a = 0._dp
612 0 : dlnT_dL_a = 0._dp
613 0 : dlnT_dlnR_a = 0._dp
614 0 : dlnT_dlnM_a = 0._dp
615 0 : dlnT_dlnkap_a = 0._dp
616 :
617 0 : lnP_a = 0._dp
618 0 : dlnP_dL_a = 0._dp
619 0 : dlnP_dlnR_a = 0._dp
620 0 : dlnP_dlnM_a = 0._dp
621 0 : dlnP_dlnkap_a = 0._dp
622 :
623 : end if
624 :
625 : ! Blend the results together
626 :
627 0 : lnT = alfa*lnT_a + beta*lnT_b
628 0 : lnP = alfa*lnP_a + beta*lnP_b
629 :
630 0 : if (.not. skip_partials) then
631 :
632 0 : dlnT_dL = alfa*dlnT_dL_a + beta*dlnT_dL_b
633 0 : dlnT_dlnR = alfa*dlnT_dlnR_a + beta*dlnT_dlnR_b
634 0 : dlnT_dlnM = alfa*dlnT_dlnM_a + beta*dlnT_dlnM_b
635 0 : dlnT_dlnkap = alfa*dlnT_dlnkap_a + beta*dlnT_dlnkap_b
636 :
637 0 : dlnP_dL = alfa*dlnP_dL_a + beta*dlnP_dL_b
638 0 : dlnP_dlnR = alfa*dlnP_dlnR_a + beta*dlnP_dlnR_b
639 0 : dlnP_dlnM = alfa*dlnP_dlnM_a + beta*dlnP_dlnM_b
640 0 : dlnP_dlnkap = alfa*dlnP_dlnkap_a + beta*dlnP_dlnkap_b
641 :
642 : end if
643 :
644 : ! Set the effective temperature
645 :
646 : ! if (alfa /= 0._dp) then
647 : ! if (beta /= 0._dp) then
648 : ! if (Teff_a /= Teff_b) then
649 : ! ierr = -1
650 : ! s% retry_message = 'Mismatch between Teff values in get_tables'
651 : ! write(*,*) 'Mismatch between Teff values in get_tables: ', Teff_a, Teff_b
652 : ! !call mesa_error(__FILE__,__LINE__)
653 : ! end if
654 : ! end if
655 : ! Teff = Teff_a
656 : ! else
657 : ! Teff = Teff_b
658 : ! end if
659 :
660 : return
661 :
662 : end subroutine get_table
663 :
664 :
665 0 : subroutine get_table_id (table_name, table_id, ierr)
666 :
667 : use atm_def, only: &
668 : ATM_TABLE_TAU_100, &
669 : ATM_TABLE_TAU_10, &
670 : ATM_TABLE_TAU_1, &
671 : ATM_TABLE_TAU_1M1, &
672 : ATM_TABLE_PHOTOSPHERE, &
673 : ATM_TABLE_WD_TAU_25, &
674 : ATM_TABLE_DB_WD_TAU_25
675 :
676 : character(*), intent(in) :: table_name
677 : integer, intent(out) :: table_id
678 : integer, intent(out) :: ierr
679 :
680 0 : ierr = 0
681 :
682 : ! Get the table id
683 :
684 0 : select case (table_name)
685 : case ('tau_100')
686 0 : table_id = ATM_TABLE_TAU_100
687 : case ('tau_10')
688 0 : table_id = ATM_TABLE_TAU_10
689 : case ('tau_1')
690 0 : table_id = ATM_TABLE_TAU_1
691 : case ('tau_1m1')
692 0 : table_id = ATM_TABLE_TAU_1M1
693 : case ('photosphere')
694 0 : table_id = ATM_TABLE_PHOTOSPHERE
695 : case ('WD_tau_25')
696 0 : table_id = ATM_TABLE_WD_TAU_25
697 : case ('DB_WD_tau_25')
698 0 : table_id = ATM_TABLE_DB_WD_TAU_25
699 : case default
700 0 : write(*,*) 'Unknown value for atm_table: ' // trim(table_name)
701 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
702 : end select
703 :
704 0 : return
705 :
706 : end subroutine get_table_id
707 :
708 :
709 0 : subroutine get_irradiated( &
710 : s, irradiated_opacity, skip_partials, L, R, M, cgrav, &
711 : Teff, &
712 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
713 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
714 : ierr)
715 :
716 : use atm_lib, only: atm_eval_irradiated
717 : use kap_support, only: prepare_kap
718 :
719 : type (star_info), pointer :: s
720 : character(*), intent(in) :: irradiated_opacity
721 : logical, intent(in) :: skip_partials
722 : real(dp), intent(in) :: L, R, M, cgrav
723 : real(dp), intent(in) :: Teff
724 : real(dp), intent(out) :: lnT_surf
725 : real(dp), intent(out) :: dlnT_dL
726 : real(dp), intent(out) :: dlnT_dlnR
727 : real(dp), intent(out) :: dlnT_dlnM
728 : real(dp), intent(out) :: dlnT_dlnkap
729 : real(dp), intent(out) :: lnP_surf
730 : real(dp), intent(out) :: dlnP_dL
731 : real(dp), intent(out) :: dlnP_dlnR
732 : real(dp), intent(out) :: dlnP_dlnM
733 : real(dp), intent(out) :: dlnP_dlnkap
734 : integer, intent(out) :: ierr
735 :
736 : real(dp) :: kap_for_atm
737 : real(dp) :: kap
738 : real(dp) :: tau_surf
739 :
740 : include 'formats'
741 :
742 0 : if (s% solver_iter > 0) then
743 0 : kap_for_atm = s% opacity_start(1)
744 : else
745 0 : kap_for_atm = s% opacity(1)
746 : end if
747 :
748 : ! Sanity check on L
749 :
750 0 : if (L < 0._dp) then
751 0 : s% retry_message = 'get_irradiated: L <= 0'
752 0 : if (s% report_ierr) then
753 0 : write(*,2) 'get_irradiated: L <= 0', s% model_number, L
754 : !call mesa_error(__FILE__,__LINE__)
755 : end if
756 0 : ierr = -1
757 0 : return
758 : end if
759 :
760 : ! Evaluate temperature and pressure by dispatching to the
761 : ! appropriate atm_lib routine
762 :
763 0 : select case (irradiated_opacity)
764 :
765 : case ('fixed')
766 :
767 : call atm_eval_irradiated( &
768 : L, R, M, cgrav, s% atm_irradiated_T_eq, s% atm_irradiated_P_surf, &
769 : kap_for_atm, s% atm_irradiated_kap_v, s% atm_irradiated_kap_v_div_kap_th, &
770 : eos_proc_for_get_irradiated, kap_proc_for_get_irradiated, &
771 : s% atm_irradiated_errtol, 0, skip_partials, &
772 : Teff, kap, tau_surf, &
773 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
774 0 : ierr)
775 :
776 : case ('iterated')
777 :
778 0 : call prepare_kap(s, ierr)
779 0 : if (ierr /= 0) then
780 0 : s% retry_message = 'Failed in call to prepare_kap'
781 0 : if (s% report_ierr) write(*, *) s% retry_message
782 0 : return
783 : end if
784 :
785 : call atm_eval_irradiated( &
786 : L, R, M, cgrav, s% atm_irradiated_T_eq, s% atm_irradiated_P_surf, &
787 : kap_for_atm, s% atm_irradiated_kap_v, s% atm_irradiated_kap_v_div_kap_th, &
788 : eos_proc_for_get_irradiated, kap_proc_for_get_irradiated, &
789 : s% atm_irradiated_errtol, s% atm_irradiated_max_iters, skip_partials, &
790 : Teff, kap, tau_surf, &
791 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
792 0 : ierr)
793 :
794 : case default
795 :
796 0 : write(*,*) 'Unknown value for atm_irradiated_opacity: ' // trim(irradiated_opacity)
797 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
798 :
799 : end select
800 :
801 : ! Set up remaining values
802 :
803 0 : lnP_surf = log(s% atm_irradiated_P_surf)
804 :
805 0 : dlnP_dL = 0._dp
806 0 : dlnP_dlnR = 0._dp
807 0 : dlnP_dlnM = 0._dp
808 0 : dlnP_dlnkap = 0._dp
809 :
810 : ! Update tau_factor
811 :
812 0 : s% tau_factor = tau_surf/s% tau_base
813 :
814 0 : return
815 :
816 : contains
817 :
818 0 : subroutine eos_proc_for_get_irradiated( &
819 : lnP, lnT, &
820 0 : lnRho, res, dres_dlnRho, dres_dlnT, &
821 : ierr)
822 :
823 : real(dp), intent(in) :: lnP
824 : real(dp), intent(in) :: lnT
825 : real(dp), intent(out) :: lnRho
826 : real(dp), intent(out) :: res(:)
827 : real(dp), intent(out) :: dres_dlnRho(:)
828 : real(dp), intent(out) :: dres_dlnT(:)
829 : integer, intent(out) :: ierr
830 :
831 : call eos_proc( &
832 : s, lnP, lnT, &
833 : lnRho, res, dres_dlnRho, dres_dlnT, &
834 0 : ierr)
835 :
836 0 : end subroutine eos_proc_for_get_irradiated
837 :
838 :
839 0 : subroutine kap_proc_for_get_irradiated( &
840 0 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
841 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
842 : ierr)
843 :
844 : real(dp), intent(in) :: lnRho
845 : real(dp), intent(in) :: lnT
846 : real(dp), intent(in) :: res(:)
847 : real(dp), intent(in) :: dres_dlnRho(:)
848 : real(dp), intent(in) :: dres_dlnT(:)
849 : real(dp), intent(out) :: kap
850 : real(dp), intent(out) :: dlnkap_dlnRho
851 : real(dp), intent(out) :: dlnkap_dlnT
852 : integer, intent(out) :: ierr
853 :
854 : call kap_proc( &
855 : s, lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
856 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
857 0 : ierr)
858 :
859 0 : end subroutine kap_proc_for_get_irradiated
860 :
861 : end subroutine get_irradiated
862 :
863 :
864 0 : subroutine get_legacy (s, ierr)
865 :
866 : type(star_info), pointer :: s
867 : integer, intent(out) :: ierr
868 :
869 : ! Handles legacy / invalid choices for atm_option. This is
870 : ! to guide users toward the newer atmosphere controls; in the long
871 : ! term, it can be removed
872 :
873 0 : select case (s% atm_option)
874 :
875 : case ('simple_photosphere')
876 :
877 0 : write(*,*) 'The ''simple_photosphere'' choice for atm_option is deprecated'
878 0 : write(*,*) 'To achieve the same results, set:'
879 0 : write(*,*) ' atm_option = ''T_tau'''
880 0 : write(*,*) ' atm_T_tau_relation = ''Eddington'''
881 0 : write(*,*) ' atm_T_tau_opacity = ''fixed'''
882 :
883 : case ('grey_and_kap')
884 :
885 0 : write(*,*) 'The ''grey_and_kap'' choice for atm_option is deprecated'
886 0 : write(*,*) 'To achieve the same results, set:'
887 0 : write(*,*) ' atm_option = ''T_tau'''
888 0 : write(*,*) ' atm_T_tau_relation = ''Eddington'''
889 0 : write(*,*) ' atm_T_tau_opacity = ''iterated'''
890 :
891 : case ('Eddington_grey')
892 :
893 0 : write(*,*) 'The ''Eddington_grey'' choice for atm_option is deprecated'
894 0 : write(*,*) 'To achieve the same results, set:'
895 0 : write(*,*) ' atm_option = ''T_tau'''
896 0 : write(*,*) ' atm_T_tau_relation = ''Eddington'''
897 0 : write(*,*) ' atm_T_tau_opacity = ''varying'''
898 :
899 : case ('solar_Hopf_grey')
900 :
901 0 : write(*,*) 'The ''solar_Hopf_grey'' choice for atm_option is deprecated'
902 0 : write(*,*) 'To achieve the same results, set:'
903 0 : write(*,*) ' atm_option = ''T_tau'''
904 0 : write(*,*) ' atm_T_tau_relation = ''solar_Hopf'''
905 0 : write(*,*) ' atm_T_tau_opacity = ''varying'''
906 :
907 : case ('Krishna_Swamy')
908 :
909 0 : write(*,*) 'The ''Krishna_Swamy'' choice for atm_option is deprecated'
910 0 : write(*,*) 'To achieve the same results, set:'
911 0 : write(*,*) ' atm_option = ''T_tau'''
912 0 : write(*,*) ' atm_T_tau_relation = ''Krishna_Swamy'''
913 0 : write(*,*) ' atm_T_tau_opacity = ''varying'''
914 :
915 : case ('tau_100_tables')
916 :
917 0 : write(*,*) 'The ''tau_100_tables'' choice for atm_option is deprecated'
918 0 : write(*,*) 'To achieve the same results, set:'
919 0 : write(*,*) ' atm_option = ''table'''
920 0 : write(*,*) ' atm_table = ''tau_100'''
921 :
922 : case ('tau_10_tables')
923 :
924 0 : write(*,*) 'The ''tau_10_tables'' choice for atm_option is deprecated'
925 0 : write(*,*) 'To achieve the same results, set:'
926 0 : write(*,*) ' atm_option = ''table'''
927 0 : write(*,*) ' atm_table = ''tau_10'''
928 :
929 : case ('tau_1_tables')
930 :
931 0 : write(*,*) 'The ''tau_1_tables'' choice for atm_option is deprecated'
932 0 : write(*,*) 'To achieve the same results, set:'
933 0 : write(*,*) ' atm_option = ''table'''
934 0 : write(*,*) ' atm_table = ''tau_1'''
935 :
936 : case ('tau_1m1_tables')
937 :
938 0 : write(*,*) 'The ''tau_1m1_tables'' choice for atm_option is deprecated'
939 0 : write(*,*) 'To achieve the same results, set:'
940 0 : write(*,*) ' atm_option = ''table'''
941 0 : write(*,*) ' atm_table = ''tau_1m1'''
942 :
943 : case ('photosphere_tables')
944 :
945 0 : write(*,*) 'The ''tau_10_tables'' choice for atm_option is deprecated'
946 0 : write(*,*) 'To achieve the same results, set:'
947 0 : write(*,*) ' atm_option = ''table'''
948 0 : write(*,*) ' atm_table = ''photosphere'''
949 :
950 : case ('grey_irradiated')
951 :
952 0 : write(*,*) 'The ''grey_irradiated'' choice for atm_option is deprecated'
953 0 : write(*,*) 'To achieve the same results, set:'
954 0 : write(*,*) ' atm_option = ''irradiated_grey'''
955 0 : write(*,*) ' atm_irradiated_opacity = ''fixed'' (if atm_grey_irradiated_simple_kap_th = .true.)'
956 0 : write(*,*) ' atm_irradiated_opacity = ''varying'' (if atm_grey_irradiated_simple_kap_th = .false.)'
957 :
958 : case default
959 :
960 0 : write(*,*) 'Unknown value for atm_option: ' // trim(s% atm_option)
961 :
962 : end select
963 :
964 : ! Stop because we can't continue
965 :
966 0 : ierr = -1 ! ifort complains if this isn't set
967 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
968 :
969 : return
970 :
971 : end subroutine get_legacy
972 :
973 :
974 0 : subroutine build_atm( &
975 : s, L, R, Teff, M, cgrav, ierr)
976 :
977 : type(star_info), pointer :: s
978 : real(dp), intent(in) :: L, R, Teff, M, cgrav
979 : integer, intent(out) :: ierr
980 :
981 : ! Create the atmosphere structure by dispatching to the
982 : ! appropriate internal routine
983 :
984 0 : select case (s% atm_option)
985 :
986 : case ('T_tau')
987 :
988 : call build_T_tau( &
989 : s, s% tau_factor*s% tau_base, L, R, Teff, M, cgrav, &
990 : s% atm_T_tau_relation, s% atm_T_tau_opacity, &
991 : s% atm_structure_num_pts, s% atm_structure, &
992 0 : ierr)
993 :
994 : case default
995 :
996 0 : write(*,*) 'Cannot create atm structure for atm_option: ' // TRIM(s% atm_option)
997 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
998 :
999 : end select
1000 :
1001 0 : end subroutine build_atm
1002 :
1003 :
1004 0 : subroutine build_T_tau( &
1005 : s, tau_surf, L, R, Teff, M, cgrav, T_tau_relation, T_tau_opacity, &
1006 : atm_structure_num_pts, atm_structure, &
1007 : ierr)
1008 :
1009 : use atm_lib, only: &
1010 : atm_build_T_tau_uniform, &
1011 : atm_build_T_tau_varying
1012 :
1013 : type(star_info), pointer :: s
1014 : real(dp), intent(in) :: tau_surf, L, R, Teff, M, cgrav
1015 : character(*), intent(in) :: T_tau_relation
1016 : character(*), intent(in) :: T_tau_opacity
1017 : integer, intent(out) :: atm_structure_num_pts
1018 : real(dp), pointer :: atm_structure(:,:)
1019 : integer, intent(out) :: ierr
1020 :
1021 : real(dp) :: kap
1022 : real(dp) :: lnT_surf
1023 : real(dp) :: dlnT_dL
1024 : real(dp) :: dlnT_dlnR
1025 : real(dp) :: dlnT_dlnM
1026 : real(dp) :: dlnT_dlnkap
1027 : real(dp) :: lnP_surf
1028 : real(dp) :: dlnP_dL
1029 : real(dp) :: dlnP_dlnR
1030 : real(dp) :: dlnP_dlnM
1031 : real(dp) :: dlnP_dlnkap
1032 : integer :: T_tau_id
1033 :
1034 : ierr = 0
1035 :
1036 : ! First evaluate the temperature, pressure and opacity
1037 :
1038 : call get_T_tau( &
1039 : s, tau_surf, L, R, M, cgrav, T_tau_relation, T_tau_opacity, .TRUE., &
1040 : Teff, kap, &
1041 : lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
1042 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
1043 0 : ierr)
1044 0 : if (ierr /= 0) then
1045 0 : s% retry_message = 'Call to get_T_tau failed in build_T_tau'
1046 0 : if (s% report_ierr) write(*, *) s% retry_message
1047 0 : return
1048 : end if
1049 :
1050 : ! Get the T-tau id
1051 :
1052 0 : call get_T_tau_id(T_tau_relation, T_tau_id, ierr)
1053 0 : if (ierr /= 0) then
1054 0 : s% retry_message = 'Call to get_T_tau_id failed in build_T_tau'
1055 0 : if (s% report_ierr) write(*, *) s% retry_message
1056 0 : return
1057 : end if
1058 :
1059 : ! Build the atmosphere structure by dispatching to the
1060 : ! appropriate atm_lib routine, with the supplied T-tau relation
1061 :
1062 0 : select case (T_tau_opacity)
1063 :
1064 : case ('fixed', 'iterated')
1065 :
1066 : call atm_build_T_tau_uniform( &
1067 : tau_surf, L, R, Teff, M, cgrav, kap, s% Pextra_factor, s% atm_build_tau_outer, &
1068 : T_tau_id, eos_proc_for_build_T_tau, kap_proc_for_build_T_tau, &
1069 : s% atm_build_errtol, s% atm_build_dlogtau, &
1070 : atm_structure_num_pts, atm_structure, &
1071 0 : ierr)
1072 0 : if (ierr /= 0) then
1073 0 : s% retry_message = 'Call to atm_build_T_tau_uniform failed in build_T_tau'
1074 0 : if (s% report_ierr) write(*, *) s% retry_message
1075 0 : return
1076 : end if
1077 :
1078 : case ('varying')
1079 :
1080 : call atm_build_T_tau_varying( &
1081 : tau_surf, L, R, Teff, M, cgrav, lnP_surf, s% atm_build_tau_outer, &
1082 : T_tau_id, eos_proc_for_build_T_tau, kap_proc_for_build_T_tau, &
1083 : s% atm_build_errtol, s% atm_build_dlogtau, &
1084 : atm_structure_num_pts, atm_structure, &
1085 0 : ierr)
1086 0 : if (ierr /= 0) then
1087 0 : s% retry_message = 'Call to atm_build_T_tau_varying failed in build_T_tau'
1088 0 : if (s% report_ierr) write(*, *) s% retry_message
1089 0 : return
1090 : end if
1091 :
1092 : case default
1093 :
1094 0 : write(*,*) 'Unknown value for atm_T_tau_opacity: ' // trim(T_tau_opacity)
1095 0 : call mesa_error(__FILE__,__LINE__,'Please amend your inlist file to correct this problem')
1096 :
1097 : end select
1098 :
1099 : return
1100 :
1101 : contains
1102 :
1103 0 : subroutine eos_proc_for_build_T_tau( &
1104 : lnP, lnT, &
1105 0 : lnRho, res, dres_dlnRho, dres_dlnT, &
1106 : ierr)
1107 :
1108 : real(dp), intent(in) :: lnP
1109 : real(dp), intent(in) :: lnT
1110 : real(dp), intent(out) :: lnRho
1111 : real(dp), intent(out) :: res(:)
1112 : real(dp), intent(out) :: dres_dlnRho(:)
1113 : real(dp), intent(out) :: dres_dlnT(:)
1114 : integer, intent(out) :: ierr
1115 :
1116 : call eos_proc( &
1117 : s, lnP, lnT, &
1118 : lnRho, res, dres_dlnRho, dres_dlnT, &
1119 0 : ierr)
1120 :
1121 0 : end subroutine eos_proc_for_build_T_tau
1122 :
1123 :
1124 0 : subroutine kap_proc_for_build_T_tau( &
1125 0 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
1126 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
1127 : ierr)
1128 :
1129 : real(dp), intent(in) :: lnRho
1130 : real(dp), intent(in) :: lnT
1131 : real(dp), intent(in) :: res(:)
1132 : real(dp), intent(in) :: dres_dlnRho(:)
1133 : real(dp), intent(in) :: dres_dlnT(:)
1134 : real(dp), intent(out) :: kap
1135 : real(dp), intent(out) :: dlnkap_dlnRho
1136 : real(dp), intent(out) :: dlnkap_dlnT
1137 : integer, intent(out) :: ierr
1138 :
1139 : call kap_proc( &
1140 : s, lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
1141 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
1142 0 : ierr)
1143 :
1144 0 : end subroutine kap_proc_for_build_T_tau
1145 :
1146 : end subroutine build_T_tau
1147 :
1148 :
1149 0 : subroutine eos_proc( &
1150 : s, lnP, lnT, &
1151 : lnRho, res, dres_dlnRho, dres_dlnT, &
1152 : ierr)
1153 :
1154 : use eos_def, only: num_eos_basic_results
1155 : use eos_lib, only: radiation_pressure, eos_gamma_PT_get_rho_energy
1156 : use eos_support, only: solve_eos_given_PgasT
1157 :
1158 : type(star_info), pointer :: s
1159 : real(dp), intent(in) :: lnP
1160 : real(dp), intent(in) :: lnT
1161 : real(dp), intent(out) :: lnRho
1162 : real(dp), intent(out) :: res(num_eos_basic_results)
1163 : real(dp), intent(out) :: dres_dlnRho(num_eos_basic_results)
1164 : real(dp), intent(out) :: dres_dlnT(num_eos_basic_results)
1165 : integer, intent(out) :: ierr
1166 :
1167 : real(dp), parameter :: LOGRHO_TOL = 1E-11_dp
1168 : real(dp), parameter :: LOGPGAS_TOL = 1E-11_dp
1169 :
1170 : real(dp) :: T
1171 : real(dp) :: P
1172 : real(dp) :: Prad
1173 : real(dp) :: Pgas
1174 : real(dp) :: rho, gamma, energy
1175 : real(dp) :: logRho, logRho_guess
1176 0 : real(dp) :: dres_dxa(num_eos_basic_results,s% species)
1177 :
1178 0 : T = exp(lnT)
1179 0 : P = exp(lnP)
1180 :
1181 0 : Prad = radiation_pressure(T)
1182 0 : Pgas = MAX(1.E-99_dp, P - Prad)
1183 :
1184 0 : gamma = 5d0/3d0
1185 : call eos_gamma_PT_get_rho_energy( &
1186 0 : s% abar(1), P, T, gamma, rho, energy, ierr)
1187 0 : if (ierr /= 0) then
1188 0 : s% retry_message = 'Call to eos_gamma_PT_get_rho_energy failed in eos_proc'
1189 0 : if (s% report_ierr) write(*, *) trim(s% retry_message)
1190 : end if
1191 :
1192 0 : logRho_guess = log10(rho)
1193 :
1194 : call solve_eos_given_PgasT( &
1195 : s, 0, s% xa(:,1), &
1196 : lnT/ln10, log10(Pgas), logRho_guess, LOGRHO_TOL, LOGPGAS_TOL, &
1197 : logRho, res, dres_dlnRho, dres_dlnT, dres_dxa, &
1198 0 : ierr)
1199 0 : if (ierr /= 0) then
1200 0 : s% retry_message = 'Call to solve_eos_given_PgasT failed in eos_proc'
1201 0 : if (s% report_ierr) write(*, *) trim(s% retry_message)
1202 : end if
1203 :
1204 0 : lnRho = logRho*ln10
1205 :
1206 :
1207 0 : return
1208 :
1209 : end subroutine eos_proc
1210 :
1211 :
1212 0 : subroutine kap_proc( &
1213 0 : s, lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
1214 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
1215 : ierr)
1216 :
1217 : use eos_def, only: i_lnfree_e, i_eta
1218 : use kap_def, only: num_kap_fracs
1219 : use kap_support, only: get_kap
1220 :
1221 : type(star_info), pointer :: s
1222 : real(dp), intent(in) :: lnRho
1223 : real(dp), intent(in) :: lnT
1224 : real(dp), intent(in) :: res(:)
1225 : real(dp), intent(in) :: dres_dlnRho(:)
1226 : real(dp), intent(in) :: dres_dlnT(:)
1227 : real(dp), intent(out) :: kap
1228 : real(dp), intent(out) :: dlnkap_dlnRho
1229 : real(dp), intent(out) :: dlnkap_dlnT
1230 : integer, intent(out) :: ierr
1231 :
1232 : real(dp) :: kap_fracs(num_kap_fracs)
1233 :
1234 : include 'formats'
1235 :
1236 : call get_kap( &
1237 : s, 0, s% zbar(1), s% xa(:,1), lnRho/ln10, lnT/ln10, &
1238 : res(i_lnfree_e), dres_dlnRho(i_lnfree_e), dres_dlnT(i_lnfree_e), &
1239 : res(i_eta), dres_dlnRho(i_eta), dres_dlnT(i_eta), &
1240 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, &
1241 0 : ierr)
1242 0 : if (ierr /= 0) then
1243 0 : s% retry_message = 'Call to get_kap failed in kap_proc'
1244 0 : if (s% report_ierr) write(*, *) s% retry_message
1245 : end if
1246 :
1247 0 : if (kap <= 0d0 .or. is_bad(kap)) then
1248 0 : write(*,1) 'bad kap', kap
1249 0 : write(*,1) 's% zbar(1)', s% zbar(1)
1250 0 : write(*,1) 'lnRho/ln10', lnRho/ln10
1251 0 : write(*,1) 'lnT/ln10', lnT/ln10
1252 0 : write(*,1) 'res(i_eta)', res(i_eta)
1253 0 : ierr = -1
1254 0 : return
1255 : call mesa_error(__FILE__,__LINE__,'atm kap_proc')
1256 : end if
1257 :
1258 : return
1259 :
1260 : end subroutine kap_proc
1261 :
1262 : end module atm_support
|
