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 92 : 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 : real(dp) :: lnRho
78 : real(dp) :: res(num_eos_basic_results)
79 : real(dp) :: dres_dlnRho(num_eos_basic_results)
80 : real(dp) :: dres_dlnT(num_eos_basic_results)
81 : real(dp) :: dlnkap_dlnRho
82 : real(dp) :: dlnkap_dlnT
83 : real(dp) :: kap_prev
84 : real(dp) :: err
85 : real(dp) :: chiRho
86 : real(dp) :: chiT
87 : real(dp) :: dlnkap_dlnP_T
88 : real(dp) :: dlnkap_dlnT_P
89 : real(dp) :: dlnP_dL_
90 : real(dp) :: dlnT_dL_
91 : real(dp) :: dlnP_dlnR_
92 : real(dp) :: dlnT_dlnR_
93 : real(dp) :: dlnP_dlnM_
94 : real(dp) :: dlnT_dlnM_
95 :
96 : include 'formats'
97 :
98 : ierr = 0
99 :
100 : ! Sanity checks
101 :
102 46 : 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 46 : g = cgrav*M/(R*R)
111 :
112 : ! Evaluate atmosphere data at optical depth tau_surf,
113 : ! using kap_guess as the opacity
114 :
115 46 : 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 46 : ierr)
123 46 : 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 46 : 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 46 : 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 56 : 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 11 : ierr)
148 11 : 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 11 : kap_prev = kap
157 :
158 : call kap_proc( &
159 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
160 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
161 11 : ierr)
162 11 : 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 11 : err = abs(kap_prev - kap)/(errtol + errtol*kap)
170 :
171 11 : if (err < 1._dp) exit iterate_loop
172 :
173 10 : 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 10 : ierr)
183 56 : 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 46 : 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 46 : 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 : 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 : 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 : 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 : 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 : real(dp) :: lntau_surf
293 : real(dp) :: lntau_outer
294 : real(dp) :: rtol(NUM_VARS)
295 : real(dp) :: atol(NUM_VARS)
296 : real(dp) :: dlntau
297 : 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 : 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 : real(dp) :: atm_structure_sgl(num_results_for_build_atm)
394 : real(dp) :: tau
395 : 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 : real(dp) :: atm_structure_sgl(num_results_for_build_atm)
448 :
449 : ierr = 0
450 0 : irtrn = 0
451 :
452 : ! Evaluate structure data
453 :
454 0 : call build_data(x, y(1), atm_structure_sgl, ierr)
455 0 : if (ierr /= 0) then
456 0 : irtrn = -1
457 0 : return
458 : end if
459 :
460 : ! If necessary, expand arrays
461 :
462 0 : atm_structure_num_pts = atm_structure_num_pts + 1
463 0 : sz = size(atm_structure, dim=2)
464 :
465 0 : if (atm_structure_num_pts > sz) then
466 0 : sz = 2*sz + 100
467 : call realloc_double2( &
468 0 : atm_structure,num_results_for_build_atm,sz,ierr)
469 : end if
470 :
471 : ! Store data
472 :
473 0 : atm_structure(:,atm_structure_num_pts) = atm_structure_sgl
474 :
475 : return
476 :
477 : end subroutine build_solout
478 :
479 :
480 0 : subroutine build_data(lntau, delta_r, atm_structure_sgl, ierr)
481 :
482 : use atm_def
483 : use atm_utils, only: eval_Paczynski_gradr
484 : use eos_def
485 :
486 : real(dp), intent(in) :: lntau
487 : real(dp), intent(in) :: delta_r
488 : real(dp), intent(out) :: atm_structure_sgl(:)
489 : integer, intent(out) :: ierr
490 :
491 : real(dp) :: tau
492 : real(dp) :: lnT
493 : real(dp) :: dlnT_dL
494 : real(dp) :: dlnT_dlnR
495 : real(dp) :: dlnT_dlnM
496 : real(dp) :: dlnT_dlnkap
497 : real(dp) :: lnP
498 : real(dp) :: dlnP_dL
499 : real(dp) :: dlnP_dlnR
500 : real(dp) :: dlnP_dlnM
501 : real(dp) :: dlnP_dlnkap
502 : real(dp) :: lnRho
503 : real(dp) :: res(num_eos_basic_results)
504 : real(dp) :: dres_dlnRho(num_eos_basic_results)
505 : real(dp) :: dres_dlnT(num_eos_basic_results)
506 : real(dp) :: gradr
507 :
508 0 : ierr = 0
509 :
510 : ! Evaluate temperature and pressure at optical depth tau
511 :
512 0 : tau = exp(lntau)
513 :
514 : call eval_data( &
515 : tau, Teff, g, L, M, cgrav, &
516 : kap, Pextra_factor, T_tau_id, .TRUE., &
517 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
518 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
519 0 : ierr)
520 0 : if (ierr /= 0) then
521 0 : write(*,*) 'atm: Call to eval_data failed in build_data'
522 0 : return
523 : end if
524 :
525 : ! Calculate the density & eos results
526 :
527 : call eos_proc( &
528 : lnP, lnT, &
529 : lnRho, res, dres_dlnRho, dres_dlnT, &
530 0 : ierr)
531 0 : if (ierr /= 0) then
532 0 : write(*,*) 'atm: Call to eos_proc failed in build_data'
533 0 : return
534 : end if
535 :
536 : ! Evaluate radiative temperature gradient
537 :
538 0 : gradr = eval_Paczynski_gradr(exp(lnT), exp(lnP), exp(lnRho), tau, kap, L, M, R, cgrav)
539 :
540 : ! Store data
541 :
542 0 : atm_structure_sgl(atm_xm) = 0._dp ! We assume negligible mass in the atmosphere
543 0 : atm_structure_sgl(atm_delta_r) = delta_r
544 0 : atm_structure_sgl(atm_lnP) = lnP
545 0 : atm_structure_sgl(atm_lnd) = lnRho
546 0 : atm_structure_sgl(atm_lnT) = lnT
547 0 : atm_structure_sgl(atm_gradT) = gradr ! by assumption, atm is radiative
548 0 : atm_structure_sgl(atm_kap) = kap
549 0 : atm_structure_sgl(atm_gamma1) = res(i_gamma1)
550 0 : atm_structure_sgl(atm_grada) = res(i_grad_ad)
551 0 : atm_structure_sgl(atm_chiT) = res(i_chiT)
552 0 : atm_structure_sgl(atm_chiRho) = res(i_chiRho)
553 0 : atm_structure_sgl(atm_cp) = res(i_Cp)
554 0 : atm_structure_sgl(atm_cv) = res(i_Cv)
555 0 : atm_structure_sgl(atm_tau) = tau
556 0 : atm_structure_sgl(atm_lnfree_e) = res(i_lnfree_e)
557 0 : atm_structure_sgl(atm_dlnkap_dlnT) = 0._dp
558 0 : atm_structure_sgl(atm_dlnkap_dlnd) = 0._dp
559 0 : atm_structure_sgl(atm_lnPgas) = res(i_lnPgas)
560 0 : atm_structure_sgl(atm_gradr) = gradr
561 :
562 : end subroutine build_data
563 :
564 : end subroutine build_T_tau_uniform
565 :
566 :
567 : ! Evaluate atmosphere data
568 :
569 112 : subroutine eval_data( &
570 : tau, Teff, g, L, M, cgrav, &
571 : kap, Pextra_factor, T_tau_id, skip_partials, &
572 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
573 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
574 : ierr)
575 :
576 : use atm_t_tau_relations, only: eval_T_tau
577 :
578 : real(dp), intent(in) :: tau
579 : real(dp), intent(in) :: Teff
580 : real(dp), intent(in) :: g
581 : real(dp), intent(in) :: L
582 : real(dp), intent(in) :: M
583 : real(dp), intent(in) :: cgrav
584 : real(dp), intent(in) :: kap
585 : real(dp), intent(in) :: Pextra_factor
586 : integer, intent(in) :: T_tau_id
587 : logical, intent(in) :: skip_partials
588 : real(dp), intent(out) :: lnT
589 : real(dp), intent(out) :: dlnT_dL
590 : real(dp), intent(out) :: dlnT_dlnR
591 : real(dp), intent(out) :: dlnT_dlnM
592 : real(dp), intent(out) :: dlnT_dlnkap
593 : real(dp), intent(out) :: lnP
594 : real(dp), intent(out) :: dlnP_dL
595 : real(dp), intent(out) :: dlnP_dlnR
596 : real(dp), intent(out) :: dlnP_dlnM
597 : real(dp), intent(out) :: dlnP_dlnkap
598 : integer, intent(out) :: ierr
599 :
600 : real(dp) :: P0
601 : real(dp) :: Pextra
602 : real(dp) :: Pfactor
603 : real(dp) :: P
604 : real(dp) :: dlogg_dlnR
605 : real(dp) :: dlogg_dlnM
606 : real(dp) :: dP0_dlnR
607 : real(dp) :: dP0_dlnkap
608 : real(dp) :: dP0_dL
609 : real(dp) :: dP0_dlnM
610 : real(dp) :: dPfactor_dlnR
611 : real(dp) :: dPfactor_dlnkap
612 : real(dp) :: dPfactor_dL
613 : real(dp) :: dPfactor_dlnM
614 : real(dp) :: dlnTeff_dL
615 : real(dp) :: dlnTeff_dlnR
616 : real(dp) :: dlnT_dlnTeff
617 :
618 : include 'formats'
619 :
620 56 : ierr = 0
621 :
622 : ! The analytic P(tau) relation is
623 : ! P = (tau*g/kap)*[1 + (kap/tau)*(L/M)/(6*pi*c*G)]
624 : ! The factor in square brackets comes from including nonzero Prad at tau=0
625 : ! see, e.g., Cox & Giuli, Section 20.1
626 :
627 56 : P0 = tau*g/kap
628 :
629 56 : Pextra = Pextra_factor*(kap/tau)*(L/M)/(6._dp*pi*clight*cgrav)
630 :
631 56 : Pfactor = 1._dp + Pextra
632 56 : P = P0*Pfactor
633 56 : lnP = log(P)
634 56 : if (is_bad(lnP)) then
635 0 : ierr = -1
636 0 : write(*,1) 'bad logP in atm_t_tau_uniform eval_data', lnP/ln10
637 0 : write(*,1) 'tau', tau
638 0 : write(*,1) 'g', g
639 0 : write(*,1) 'kap', kap
640 0 : write(*,1) 'P0', P0
641 0 : write(*,1) 'Pextra', Pextra
642 0 : write(*,1) 'P', P
643 : !call mesa_error(__FILE__,__LINE__,'atm_t_tau_uniform')
644 0 : return
645 : end if
646 :
647 56 : call eval_T_tau(T_tau_id, tau, Teff, lnT, ierr)
648 56 : if (ierr /= 0) then
649 0 : write(*,*) 'atm: Call to eval_T_tau failed in eval_data'
650 0 : return
651 : end if
652 :
653 : ! Set up partials
654 :
655 56 : if (.NOT. skip_partials) then
656 :
657 44 : dlogg_dlnR = -2._dp
658 44 : dlogg_dlnM = 1._dp
659 :
660 44 : dP0_dL = 0._dp
661 44 : dP0_dlnR = dlogg_dlnR*P0
662 44 : dP0_dlnM = dlogg_dlnM*P0
663 44 : dP0_dlnkap = -P0
664 :
665 44 : dPfactor_dL = Pextra/L
666 44 : dPfactor_dlnR = 0._dp
667 44 : dPfactor_dlnM = -Pextra
668 44 : dPfactor_dlnkap = Pextra
669 :
670 44 : dlnP_dL = (dP0_dL*Pfactor + P0*dPfactor_dL)/P
671 44 : dlnP_dlnR = (dP0_dlnR*Pfactor + P0*dPfactor_dlnR)/P
672 44 : dlnP_dlnM = (dP0_dlnM*Pfactor + P0*dPfactor_dlnM)/P
673 44 : dlnP_dlnkap = (dP0_dlnkap*Pfactor + P0*dPfactor_dlnkap)/P
674 :
675 44 : dlnTeff_dL = 0.25_dp/L
676 44 : dlnTeff_dlnR = -0.5_dp
677 :
678 44 : dlnT_dlnTeff = 1._dp
679 :
680 44 : dlnT_dL = dlnTeff_dL*dlnT_dlnTeff
681 44 : dlnT_dlnR = dlnTeff_dlnR*dlnT_dlnTeff
682 44 : dlnT_dlnM = 0._dp
683 44 : dlnT_dlnkap = 0._dp
684 :
685 : else
686 :
687 12 : dlnP_dL = 0._dp
688 12 : dlnP_dlnR = 0._dp
689 12 : dlnP_dlnM = 0._dp
690 12 : dlnP_dlnkap = 0._dp
691 :
692 12 : dlnT_dL = 0._dp
693 12 : dlnT_dlnR = 0._dp
694 12 : dlnT_dlnM = 0._dp
695 12 : dlnT_dlnkap = 0._dp
696 :
697 : end if
698 :
699 : return
700 :
701 : end subroutine eval_data
702 :
703 : end module atm_T_tau_uniform
|
