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_T_tau_uniform
21 :
22 : use const_def, only: dp, pi, ln10, clight, crad
23 : use math_lib
24 : use utils_lib, only: mesa_error
25 : use utils_lib, only: is_bad
26 :
27 : implicit none
28 :
29 : private
30 : public :: eval_T_tau_uniform
31 : public :: build_T_tau_uniform
32 :
33 : contains
34 :
35 : ! Evaluate atmosphere data from T-tau relation with uniform opacity
36 :
37 116 : subroutine eval_T_tau_uniform( &
38 : tau_surf, L, R, M, cgrav, kap_guess, Pextra_factor, &
39 : T_tau_id, eos_proc, kap_proc, errtol, max_iters, skip_partials, &
40 : Teff, kap, &
41 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
42 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
43 : ierr)
44 :
45 : use atm_def, only: atm_eos_iface, atm_kap_iface
46 : use eos_def, only: num_eos_basic_results, i_chiRho, i_chiT
47 :
48 : real(dp), intent(in) :: tau_surf
49 : real(dp), intent(in) :: L
50 : real(dp), intent(in) :: R
51 : real(dp), intent(in) :: M
52 : real(dp), intent(in) :: cgrav
53 : real(dp), intent(in) :: kap_guess
54 : real(dp), intent(in) :: Pextra_factor
55 : integer, intent(in) :: T_tau_id
56 : procedure(atm_eos_iface) :: eos_proc
57 : procedure(atm_kap_iface) :: kap_proc
58 : real(dp), intent(in) :: errtol
59 : integer, intent(in) :: max_iters
60 : logical, intent(in) :: skip_partials
61 : real(dp), intent(in) :: Teff
62 : real(dp), intent(out) :: kap
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 : real(dp), intent(out) :: lnP
69 : real(dp), intent(out) :: dlnP_dL
70 : real(dp), intent(out) :: dlnP_dlnR
71 : real(dp), intent(out) :: dlnP_dlnM
72 : real(dp), intent(out) :: dlnP_dlnkap
73 : integer, intent(out) :: ierr
74 :
75 : real(dp) :: g
76 : integer :: iters
77 48 : real(dp) :: lnRho
78 1296 : real(dp) :: res(num_eos_basic_results)
79 1296 : real(dp) :: dres_dlnRho(num_eos_basic_results)
80 1296 : real(dp) :: dres_dlnT(num_eos_basic_results)
81 48 : real(dp) :: dlnkap_dlnRho
82 48 : real(dp) :: dlnkap_dlnT
83 48 : real(dp) :: kap_prev
84 48 : real(dp) :: err
85 48 : real(dp) :: chiRho
86 48 : real(dp) :: chiT
87 48 : real(dp) :: dlnkap_dlnP_T
88 48 : real(dp) :: dlnkap_dlnT_P
89 48 : real(dp) :: dlnP_dL_
90 48 : real(dp) :: dlnT_dL_
91 48 : real(dp) :: dlnP_dlnR_
92 48 : real(dp) :: dlnT_dlnR_
93 48 : real(dp) :: dlnP_dlnM_
94 48 : real(dp) :: dlnT_dlnM_
95 :
96 : include 'formats'
97 :
98 : ierr = 0
99 :
100 : ! Sanity checks
101 :
102 48 : if (L <= 0._dp .OR. R <= 0._dp .OR. M <= 0._dp) then
103 0 : ierr = -1
104 0 : write(*,*) 'atm: eval_T_tau_uniform: L, R, or M bad', L, R, M
105 0 : return
106 : end if
107 :
108 : ! Evaluate the gravity
109 :
110 48 : g = cgrav*M/(R*R)
111 :
112 : ! Evaluate atmosphere data at optical depth tau_surf,
113 : ! using kap_guess as the opacity
114 :
115 48 : kap = kap_guess
116 :
117 : call eval_data( &
118 : tau_surf, Teff, g, L, M, cgrav, &
119 : kap, Pextra_factor, T_tau_id, skip_partials, &
120 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
121 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
122 48 : ierr)
123 48 : if (ierr /= 0) then
124 0 : write(*,*) 'atm: Call to eval_data failed in eval_T_tau_uniform'
125 0 : return
126 : end if
127 48 : if (is_bad(lnT)) then
128 0 : ierr = -1
129 0 : write(*,*) 'atm: eval_T_tau_uniform logT from eval_data', lnT/ln10
130 0 : return
131 : end if
132 48 : if (is_bad(lnP)) then
133 0 : ierr = -1
134 0 : write(*,*) 'atm: eval_T_tau_uniform logP from eval_data', lnP/ln10
135 0 : return
136 : end if
137 :
138 : ! Iterate to find a consistent opacity
139 :
140 68 : iterate_loop : do iters = 1, max_iters
141 :
142 : ! Calculate the density & eos results
143 :
144 : call eos_proc( &
145 : lnP, lnT, &
146 : lnRho, res, dres_dlnRho, dres_dlnT, &
147 22 : ierr)
148 22 : if (ierr /= 0) then
149 0 : write(*,*) 'atm: call to eos_proc failed in eval_T_tau_uniform logP logT logPrad', &
150 0 : lnP/ln10, lnT/ln10, log10(crad*exp(4d0*lnT)/3d0)
151 0 : return
152 : end if
153 :
154 : ! Update the opacity
155 :
156 22 : kap_prev = kap
157 :
158 : call kap_proc( &
159 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
160 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
161 22 : ierr)
162 22 : if (ierr /= 0) then
163 0 : write(*,*) 'atm: Call to kap_proc failed in eval_T_tau_uniform logRho logT', lnRho/ln10, lnT/ln10
164 0 : return
165 : end if
166 :
167 : ! Check for convergence
168 :
169 22 : err = abs(kap_prev - kap)/(errtol + errtol*kap)
170 :
171 22 : if (err < 1._dp) exit iterate_loop
172 :
173 20 : kap = kap_prev + 0.5_dp*(kap - kap_prev) ! under correct
174 :
175 : ! Re-evaluate atmosphere data
176 :
177 : call eval_data( &
178 : tau_surf, Teff, g, L, M, cgrav, &
179 : kap, Pextra_factor, T_tau_id, skip_partials, &
180 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
181 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
182 20 : ierr)
183 68 : if (ierr /= 0) then
184 0 : write(*,*) 'Call to eval_data failed in eval_T_tau_uniform'
185 0 : return
186 : end if
187 :
188 : end do iterate_loop
189 :
190 48 : if (max_iters > 0 .AND. iters > max_iters) then
191 0 : write(*,*) 'atm: Exceeded max_iters iterations in eval_T_tau_uniform'
192 0 : ierr = -1
193 0 : return
194 : end if
195 :
196 : ! If necessary, fix up the partials to account for the implicit
197 : ! dependence of the opacity on the final P and T
198 :
199 48 : if (max_iters > 0 .AND. .NOT. skip_partials ) then
200 :
201 0 : chiRho = res(i_chiRho)
202 0 : chiT = res(i_chiT)
203 :
204 0 : dlnkap_dlnP_T = dlnkap_dlnRho/chiRho
205 0 : dlnkap_dlnT_P = dlnkap_dlnT - dlnkap_dlnRho*chiT/chiRho
206 :
207 : dlnP_dL_ = (dlnP_dL + dlnkap_dlnT_P*(dlnP_dlnkap*dlnT_dL - dlnP_dL*dlnT_dlnkap))/&
208 0 : (1._dp - dlnkap_dlnP_T*dlnP_dlnkap - dlnkap_dlnT_P*dlnT_dlnkap)
209 : dlnT_dL_ = (dlnT_dL + dlnkap_dlnP_T*(dlnT_dlnkap*dlnP_dL - dlnT_dL*dlnP_dlnkap))/&
210 0 : (1._dp - dlnkap_dlnP_T*dlnP_dlnkap - dlnkap_dlnT_P*dlnT_dlnkap)
211 :
212 : dlnP_dlnR_ = (dlnP_dlnR + dlnkap_dlnT_P*(dlnP_dlnkap*dlnT_dlnR - dlnP_dlnR*dlnT_dlnkap))/ &
213 0 : (1._dp - dlnkap_dlnP_T*dlnP_dlnkap - dlnkap_dlnT_P*dlnT_dlnkap)
214 : dlnT_dlnR_ = (dlnT_dlnR + dlnkap_dlnP_T*(dlnT_dlnkap*dlnP_dlnR - dlnT_dlnR*dlnP_dlnkap))/ &
215 0 : (1._dp - dlnkap_dlnP_T*dlnP_dlnkap - dlnkap_dlnT_P*dlnT_dlnkap)
216 :
217 : dlnP_dlnM_ = (dlnP_dlnM + dlnkap_dlnT_P*(dlnP_dlnkap*dlnT_dlnM - dlnP_dlnM*dlnT_dlnkap))/ &
218 0 : (1._dp - dlnkap_dlnP_T*dlnP_dlnkap - dlnkap_dlnT_P*dlnT_dlnkap)
219 : dlnT_dlnM_ = (dlnT_dlnM + dlnkap_dlnP_T*(dlnT_dlnkap*dlnP_dlnM - dlnT_dlnM*dlnP_dlnkap))/ &
220 0 : (1._dp - dlnkap_dlnP_T*dlnP_dlnkap - dlnkap_dlnT_P*dlnT_dlnkap)
221 :
222 0 : dlnP_dL = dlnP_dL_
223 0 : dlnT_dL = dlnT_dL_
224 :
225 0 : dlnP_dlnR = dlnP_dlnR_
226 0 : dlnT_dlnR = dlnT_dlnR_
227 :
228 0 : dlnP_dlnM = dlnP_dlnM_
229 0 : dlnT_dlnM = dlnT_dlnM_
230 :
231 0 : dlnP_dlnkap = 0._dp
232 0 : dlnT_dlnkap = 0._dp
233 :
234 : end if
235 :
236 : return
237 :
238 48 : end subroutine eval_T_tau_uniform
239 :
240 :
241 : ! Build atmosphere structure data from T-tau relation with uniform
242 : ! opacity
243 :
244 0 : subroutine build_T_tau_uniform( &
245 : tau_surf, L, R, Teff, M, cgrav, kap, Pextra_factor, tau_outer, &
246 : T_tau_id, eos_proc, kap_proc, errtol, dlogtau, &
247 : atm_structure_num_pts, atm_structure, &
248 : ierr)
249 :
250 48 : use atm_def, only: atm_eos_iface, atm_kap_iface, num_results_for_build_atm
251 : use atm_utils, only: eval_Teff_g
252 : use num_lib, only: dopri5_work_sizes, dopri5
253 :
254 : real(dp), intent(in) :: tau_surf
255 : real(dp), intent(in) :: L
256 : real(dp), intent(in) :: R
257 : real(dp), intent(in) :: Teff
258 : real(dp), intent(in) :: M
259 : real(dp), intent(in) :: cgrav
260 : real(dp), intent(in) :: kap
261 : real(dp), intent(in) :: Pextra_factor
262 : real(dp), intent(in) :: tau_outer
263 : integer, intent(in) :: T_tau_id
264 : procedure(atm_eos_iface) :: eos_proc
265 : procedure(atm_kap_iface) :: kap_proc
266 : real(dp), intent(in) :: errtol
267 : real(dp), intent(in) :: dlogtau
268 : integer, intent(out) :: atm_structure_num_pts
269 : real(dp), pointer :: atm_structure(:,:)
270 : integer, intent(out) :: ierr
271 :
272 : integer, parameter :: INIT_NUM_PTS = 100
273 : integer, parameter :: NUM_VARS = 1
274 : integer, parameter :: NRDENS = 0
275 : integer, parameter :: MAX_STEPS = 0
276 : integer, parameter :: LRPAR = 0
277 : integer, parameter :: LIPAR = 0
278 : integer, parameter :: IOUT = 1
279 : integer, parameter :: LOUT = 0
280 :
281 0 : real(dp) :: g
282 : integer :: liwork
283 : integer :: lwork
284 0 : real(dp), pointer :: work(:)
285 0 : integer, pointer :: iwork(:)
286 : real(dp), target :: rpar_ary(LRPAR)
287 : integer, target :: ipar_ary(LIPAR)
288 0 : real(dp), target :: y_ary(NUM_VARS)
289 0 : real(dp), pointer :: rpar(:)
290 0 : integer, pointer :: ipar(:)
291 0 : real(dp), pointer :: y(:)
292 0 : real(dp) :: lntau_surf
293 0 : real(dp) :: lntau_outer
294 0 : real(dp) :: rtol(NUM_VARS)
295 0 : real(dp) :: atol(NUM_VARS)
296 0 : real(dp) :: dlntau
297 0 : real(dp) :: dlntau_max
298 : integer :: idid
299 :
300 0 : ierr = 0
301 :
302 : ! Sanity check
303 :
304 0 : if (dlogtau <= 0._dp) then
305 0 : write(*,*) 'atm: Invalid dlogtau in build_T_tau_uniform:', dlogtau
306 0 : call mesa_error(__FILE__,__LINE__)
307 : end if
308 :
309 : ! Evaluate the gravity
310 :
311 0 : g = cgrav*M/(R*R)
312 :
313 : ! Allocate atm_structure at its initial size
314 :
315 0 : allocate(atm_structure(num_results_for_build_atm,INIT_NUM_PTS))
316 :
317 0 : atm_structure_num_pts = 0
318 :
319 0 : if (tau_outer > tau_surf) return
320 :
321 : ! Allocate work rrays for the integrator
322 :
323 0 : call dopri5_work_sizes(NUM_VARS, NRDENS, liwork, lwork)
324 0 : allocate(work(lwork), iwork(liwork), stat=ierr)
325 0 : if (ierr /= 0) then
326 0 : write(*,*) 'atm: allocate failed in build_T_tau_uniform'
327 0 : deallocate(atm_structure)
328 0 : return
329 : end if
330 :
331 0 : work = 0._dp
332 0 : iwork = 0
333 :
334 : ! Set pointers (simply because dopri5 wants pointer args)
335 :
336 0 : rpar => rpar_ary
337 0 : ipar => ipar_ary
338 :
339 0 : y => y_ary
340 :
341 : ! Set starting values for the integrator (y(1) = delta_r)
342 :
343 0 : y(1) = 0._dp
344 :
345 : ! Integrate from the atmosphere base outward
346 :
347 0 : lntau_outer = log(tau_outer)
348 0 : lntau_surf = log(tau_surf)
349 :
350 0 : dlntau_max = -dlogtau*ln10
351 0 : dlntau = 0.5_dp*dlntau_max
352 :
353 0 : rtol = errtol
354 0 : atol = errtol
355 :
356 : call dopri5( &
357 : NUM_VARS, build_fcn, lntau_surf, y, lntau_outer, &
358 : dlntau, dlntau_max, MAX_STEPS, &
359 : rtol, atol, 1, &
360 : build_solout, IOUT, &
361 : work, lwork, iwork, liwork, &
362 : LRPAR, rpar, LIPAR, ipar, &
363 0 : LOUT, idid)
364 0 : if (idid < 0) then
365 0 : write(*,*) 'atm: Call to dopri5 failed in build_T_tau_uniform: idid=', idid
366 0 : ierr = -1
367 : end if
368 :
369 : ! Reverse the data ordering (since by convention data in atm_structure
370 : ! should be ordered outward-in)
371 :
372 : atm_structure(:,:atm_structure_num_pts) = &
373 0 : atm_structure(:,atm_structure_num_pts:1:-1)
374 :
375 : ! Deallocate arrays
376 :
377 0 : deallocate(work, iwork)
378 :
379 : contains
380 :
381 0 : subroutine build_fcn(n, x, h, y, f, lr, rpar, li, ipar, ierr)
382 :
383 0 : use atm_def
384 :
385 : integer, intent(in) :: n, lr, li
386 : real(dp), intent(in) :: x, h
387 : real(dp), intent(inout) :: y(:)
388 : real(dp), intent(inout) :: f(:)
389 : integer, intent(inout), pointer :: ipar(:)
390 : real(dp), intent(inout), pointer :: rpar(:)
391 : integer, intent(out) :: ierr
392 :
393 0 : real(dp) :: atm_structure_sgl(num_results_for_build_atm)
394 0 : real(dp) :: tau
395 0 : real(dp) :: rho
396 :
397 : ierr = 0
398 :
399 : ! Evaluate structure data
400 :
401 0 : call build_data(x, y(1), atm_structure_sgl, ierr)
402 0 : if (ierr /= 0) then
403 : ! Non-zero ierr signals we must use a smaller stepsize;
404 : ! whereas in fact we want to terminate the integration.
405 : ! Needs fixing!
406 0 : return
407 : end if
408 :
409 : ! Set up the rhs for the optical depth equation
410 : ! dr/dlntau = -tau/(kappa*rho)
411 :
412 0 : tau = exp(x)
413 0 : rho = exp(atm_structure_sgl(atm_lnd))
414 :
415 0 : f(1) = -tau/(kap*rho)
416 :
417 : end subroutine build_fcn
418 :
419 :
420 0 : subroutine build_solout( &
421 0 : nr, xold, x, n, y, rwork_y, iwork_y, interp_y, lrpar, rpar, lipar, ipar, irtrn)
422 :
423 : use utils_lib, only: realloc_double2
424 :
425 : integer, intent(in) :: nr, n, lrpar, lipar
426 : real(dp), intent(in) :: xold, x
427 : real(dp), intent(inout) :: y(:)
428 : real(dp), intent(inout), target :: rwork_y(*)
429 : integer, intent(inout), target :: iwork_y(*)
430 : integer, intent(inout), pointer :: ipar(:)
431 : real(dp), intent(inout), pointer :: rpar(:)
432 : interface
433 : real(dp) function interp_y(i, s, rwork_y, iwork_y, ierr)
434 : use const_def, only: dp
435 : implicit none
436 : integer, intent(in) :: i
437 : real(dp), intent(in) :: s
438 : real(dp), intent(inout), target :: rwork_y(*)
439 : integer, intent(inout), target :: iwork_y(*)
440 : integer, intent(out) :: ierr
441 : end function interp_y
442 : end interface
443 : integer, intent(out) :: irtrn
444 :
445 : integer :: ierr
446 : integer :: sz
447 0 : real(dp) :: atm_structure_sgl(num_results_for_build_atm)
448 :
449 : ierr = 0
450 :
451 : ! Evaluate structure data
452 :
453 0 : call build_data(x, y(1), atm_structure_sgl, ierr)
454 0 : if (ierr /= 0) then
455 0 : irtrn = -1
456 0 : return
457 : end if
458 :
459 : ! If necessary, expand arrays
460 :
461 0 : atm_structure_num_pts = atm_structure_num_pts + 1
462 0 : sz = size(atm_structure, dim=2)
463 :
464 0 : if (atm_structure_num_pts > sz) then
465 0 : sz = 2*sz + 100
466 : call realloc_double2( &
467 0 : atm_structure,num_results_for_build_atm,sz,ierr)
468 : end if
469 :
470 : ! Store data
471 :
472 0 : atm_structure(:,atm_structure_num_pts) = atm_structure_sgl
473 :
474 : return
475 :
476 : end subroutine build_solout
477 :
478 :
479 0 : subroutine build_data(lntau, delta_r, atm_structure_sgl, ierr)
480 :
481 : use atm_def
482 : use atm_utils, only: eval_Paczynski_gradr
483 : use eos_def
484 :
485 : real(dp), intent(in) :: lntau
486 : real(dp), intent(in) :: delta_r
487 : real(dp), intent(out) :: atm_structure_sgl(:)
488 : integer, intent(out) :: ierr
489 :
490 0 : real(dp) :: tau
491 0 : real(dp) :: lnT
492 : real(dp) :: dlnT_dL
493 : real(dp) :: dlnT_dlnR
494 : real(dp) :: dlnT_dlnM
495 : real(dp) :: dlnT_dlnkap
496 0 : real(dp) :: lnP
497 : real(dp) :: dlnP_dL
498 : real(dp) :: dlnP_dlnR
499 : real(dp) :: dlnP_dlnM
500 : real(dp) :: dlnP_dlnkap
501 0 : real(dp) :: lnRho
502 0 : real(dp) :: res(num_eos_basic_results)
503 0 : real(dp) :: dres_dlnRho(num_eos_basic_results)
504 0 : real(dp) :: dres_dlnT(num_eos_basic_results)
505 0 : real(dp) :: gradr
506 :
507 0 : ierr = 0
508 :
509 : ! Evaluate temperature and pressure at optical depth tau
510 :
511 0 : tau = exp(lntau)
512 :
513 : call eval_data( &
514 : tau, Teff, g, L, M, cgrav, &
515 : kap, Pextra_factor, T_tau_id, .TRUE., &
516 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
517 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
518 0 : ierr)
519 0 : if (ierr /= 0) then
520 0 : write(*,*) 'atm: Call to eval_data failed in build_data'
521 0 : return
522 : end if
523 :
524 : ! Calculate the density & eos results
525 :
526 : call eos_proc( &
527 : lnP, lnT, &
528 : lnRho, res, dres_dlnRho, dres_dlnT, &
529 0 : ierr)
530 0 : if (ierr /= 0) then
531 0 : write(*,*) 'atm: Call to eos_proc failed in build_data'
532 0 : return
533 : end if
534 :
535 : ! Evaluate radiative temperature gradient
536 :
537 0 : gradr = eval_Paczynski_gradr(exp(lnT), exp(lnP), exp(lnRho), tau, kap, L, M, R, cgrav)
538 :
539 : ! Store data
540 :
541 0 : atm_structure_sgl(atm_xm) = 0._dp ! We assume negligible mass in the atmosphere
542 0 : atm_structure_sgl(atm_delta_r) = delta_r
543 0 : atm_structure_sgl(atm_lnP) = lnP
544 0 : atm_structure_sgl(atm_lnd) = lnRho
545 0 : atm_structure_sgl(atm_lnT) = lnT
546 0 : atm_structure_sgl(atm_gradT) = gradr ! by assumption, atm is radiative
547 0 : atm_structure_sgl(atm_kap) = kap
548 0 : atm_structure_sgl(atm_gamma1) = res(i_gamma1)
549 0 : atm_structure_sgl(atm_grada) = res(i_grad_ad)
550 0 : atm_structure_sgl(atm_chiT) = res(i_chiT)
551 0 : atm_structure_sgl(atm_chiRho) = res(i_chiRho)
552 0 : atm_structure_sgl(atm_cp) = res(i_Cp)
553 0 : atm_structure_sgl(atm_cv) = res(i_Cv)
554 0 : atm_structure_sgl(atm_tau) = tau
555 0 : atm_structure_sgl(atm_lnfree_e) = res(i_lnfree_e)
556 0 : atm_structure_sgl(atm_dlnkap_dlnT) = 0._dp
557 0 : atm_structure_sgl(atm_dlnkap_dlnd) = 0._dp
558 0 : atm_structure_sgl(atm_lnPgas) = res(i_lnPgas)
559 0 : atm_structure_sgl(atm_gradr) = gradr
560 :
561 0 : end subroutine build_data
562 :
563 : end subroutine build_T_tau_uniform
564 :
565 :
566 : ! Evaluate atmosphere data
567 :
568 68 : subroutine eval_data( &
569 : tau, Teff, g, L, M, cgrav, &
570 : kap, Pextra_factor, T_tau_id, skip_partials, &
571 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
572 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
573 : ierr)
574 :
575 : use atm_t_tau_relations, only: eval_T_tau
576 :
577 : real(dp), intent(in) :: tau
578 : real(dp), intent(in) :: Teff
579 : real(dp), intent(in) :: g
580 : real(dp), intent(in) :: L
581 : real(dp), intent(in) :: M
582 : real(dp), intent(in) :: cgrav
583 : real(dp), intent(in) :: kap
584 : real(dp), intent(in) :: Pextra_factor
585 : integer, intent(in) :: T_tau_id
586 : logical, intent(in) :: skip_partials
587 : real(dp), intent(out) :: lnT
588 : real(dp), intent(out) :: dlnT_dL
589 : real(dp), intent(out) :: dlnT_dlnR
590 : real(dp), intent(out) :: dlnT_dlnM
591 : real(dp), intent(out) :: dlnT_dlnkap
592 : real(dp), intent(out) :: lnP
593 : real(dp), intent(out) :: dlnP_dL
594 : real(dp), intent(out) :: dlnP_dlnR
595 : real(dp), intent(out) :: dlnP_dlnM
596 : real(dp), intent(out) :: dlnP_dlnkap
597 : integer, intent(out) :: ierr
598 :
599 : real(dp) :: P0
600 : real(dp) :: Pextra
601 68 : real(dp) :: Pfactor
602 : real(dp) :: P
603 68 : real(dp) :: dlogg_dlnR
604 68 : real(dp) :: dlogg_dlnM
605 68 : real(dp) :: dP0_dlnR
606 68 : real(dp) :: dP0_dlnkap
607 68 : real(dp) :: dP0_dL
608 68 : real(dp) :: dP0_dlnM
609 68 : real(dp) :: dPfactor_dlnR
610 68 : real(dp) :: dPfactor_dlnkap
611 68 : real(dp) :: dPfactor_dL
612 68 : real(dp) :: dPfactor_dlnM
613 68 : real(dp) :: dlnTeff_dL
614 68 : real(dp) :: dlnTeff_dlnR
615 68 : real(dp) :: dlnT_dlnTeff
616 :
617 : include 'formats'
618 :
619 68 : ierr = 0
620 :
621 : ! The analytic P(tau) relation is
622 : ! P = (tau*g/kap)*[1 + (kap/tau)*(L/M)/(6*pi*c*G)]
623 : ! The factor in square brackets comes from including nonzero Prad at tau=0
624 : ! see, e.g., Cox & Giuli, Section 20.1
625 :
626 68 : P0 = tau*g/kap
627 :
628 68 : Pextra = Pextra_factor*(kap/tau)*(L/M)/(6._dp*pi*clight*cgrav)
629 :
630 68 : Pfactor = 1._dp + Pextra
631 68 : P = P0*Pfactor
632 68 : lnP = log(P)
633 68 : if (is_bad(lnP)) then
634 0 : ierr = -1
635 0 : write(*,1) 'bad logP in atm_t_tau_uniform eval_data', lnP/ln10
636 0 : write(*,1) 'tau', tau
637 0 : write(*,1) 'g', g
638 0 : write(*,1) 'kap', kap
639 0 : write(*,1) 'P0', P0
640 0 : write(*,1) 'Pextra', Pextra
641 0 : write(*,1) 'P', P
642 : !call mesa_error(__FILE__,__LINE__,'atm_t_tau_uniform')
643 0 : return
644 : end if
645 :
646 68 : call eval_T_tau(T_tau_id, tau, Teff, lnT, ierr)
647 68 : if (ierr /= 0) then
648 0 : write(*,*) 'atm: Call to eval_T_tau failed in eval_data'
649 0 : return
650 : end if
651 :
652 : ! Set up partials
653 :
654 68 : if (.NOT. skip_partials) then
655 :
656 44 : dlogg_dlnR = -2._dp
657 44 : dlogg_dlnM = 1._dp
658 :
659 44 : dP0_dL = 0._dp
660 44 : dP0_dlnR = dlogg_dlnR*P0
661 44 : dP0_dlnM = dlogg_dlnM*P0
662 44 : dP0_dlnkap = -P0
663 :
664 44 : dPfactor_dL = Pextra/L
665 44 : dPfactor_dlnR = 0._dp
666 44 : dPfactor_dlnM = -Pextra
667 44 : dPfactor_dlnkap = Pextra
668 :
669 44 : dlnP_dL = (dP0_dL*Pfactor + P0*dPfactor_dL)/P
670 44 : dlnP_dlnR = (dP0_dlnR*Pfactor + P0*dPfactor_dlnR)/P
671 44 : dlnP_dlnM = (dP0_dlnM*Pfactor + P0*dPfactor_dlnM)/P
672 44 : dlnP_dlnkap = (dP0_dlnkap*Pfactor + P0*dPfactor_dlnkap)/P
673 :
674 44 : dlnTeff_dL = 0.25_dp/L
675 44 : dlnTeff_dlnR = -0.5_dp
676 :
677 44 : dlnT_dlnTeff = 1._dp
678 :
679 44 : dlnT_dL = dlnTeff_dL*dlnT_dlnTeff
680 44 : dlnT_dlnR = dlnTeff_dlnR*dlnT_dlnTeff
681 44 : dlnT_dlnM = 0._dp
682 44 : dlnT_dlnkap = 0._dp
683 :
684 : else
685 :
686 24 : dlnP_dL = 0._dp
687 24 : dlnP_dlnR = 0._dp
688 24 : dlnP_dlnM = 0._dp
689 24 : dlnP_dlnkap = 0._dp
690 :
691 24 : dlnT_dL = 0._dp
692 24 : dlnT_dlnR = 0._dp
693 24 : dlnT_dlnM = 0._dp
694 24 : dlnT_dlnkap = 0._dp
695 :
696 : end if
697 :
698 : return
699 :
700 68 : end subroutine eval_data
701 :
702 : end module atm_T_tau_uniform
|
