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_varying
21 :
22 : use const_def, only: dp, ln10, cgas
23 : use math_lib
24 : use utils_lib, only: mesa_error
25 :
26 : implicit none
27 :
28 : private
29 : public :: eval_T_tau_varying
30 : public :: build_T_tau_varying
31 :
32 : contains
33 :
34 : ! Evaluate atmosphere data from T-tau relation with varying opacity
35 :
36 10 : subroutine eval_T_tau_varying( &
37 : tau_surf, L, R, M, cgrav, &
38 : T_tau_id, eos_proc, kap_proc, &
39 : errtol, max_steps, skip_partials, &
40 : Teff, &
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 :
47 : real(dp), intent(in) :: tau_surf
48 : real(dp), intent(in) :: L
49 : real(dp), intent(in) :: R
50 : real(dp), intent(in) :: M
51 : real(dp), intent(in) :: cgrav
52 : integer, intent(in) :: T_tau_id
53 : procedure(atm_eos_iface) :: eos_proc
54 : procedure(atm_kap_iface) :: kap_proc
55 : real(dp), intent(in) :: errtol
56 : integer, intent(in) :: max_steps
57 : logical, intent(in) :: skip_partials
58 : real(dp), intent(in) :: Teff
59 : real(dp), intent(out) :: lnT
60 : real(dp), intent(out) :: dlnT_dL
61 : real(dp), intent(out) :: dlnT_dlnR
62 : real(dp), intent(out) :: dlnT_dlnM
63 : real(dp), intent(out) :: dlnT_dlnkap
64 : real(dp), intent(out) :: lnP
65 : real(dp), intent(out) :: dlnP_dL
66 : real(dp), intent(out) :: dlnP_dlnR
67 : real(dp), intent(out) :: dlnP_dlnM
68 : real(dp), intent(out) :: dlnP_dlnkap
69 : integer, intent(out) :: ierr
70 :
71 : real(dp), parameter :: DLNTEFF = 1.E-4_dp
72 :
73 : real(dp) :: g
74 10 : real(dp) :: lnTeff
75 10 : real(dp) :: lnTeff_p
76 10 : real(dp) :: lnTeff_m
77 10 : real(dp) :: lnT_p
78 10 : real(dp) :: lnT_m
79 10 : real(dp) :: lnP_p
80 10 : real(dp) :: lnP_m
81 : integer :: ierr_p
82 : integer :: ierr_m
83 10 : real(dp) :: dlnTeff_dlnR
84 10 : real(dp) :: dlnTeff_dL
85 10 : real(dp) :: dlnT_dlnTeff
86 10 : real(dp) :: dlnP_dlnTeff
87 :
88 10 : ierr = 0
89 :
90 : ! Sanity checks
91 :
92 10 : if (L <= 0._dp .OR. R <= 0._dp .OR. M <= 0._dp) then
93 0 : ierr = -1
94 0 : return
95 : end if
96 :
97 : ! Evaluate the gravity
98 :
99 10 : g = cgrav*M/(R*R)
100 :
101 : ! Perform the necessary evaluations
102 :
103 10 : if (skip_partials) then
104 :
105 : ! No partial evaluation required
106 :
107 : call eval_data( &
108 : tau_surf, Teff, g, &
109 : T_tau_id, eos_proc, kap_proc, errtol, max_steps, &
110 10 : lnT, lnP, ierr)
111 :
112 10 : if (ierr /= 0) then
113 0 : write(*,*) 'atm: Call to eval_data failed in atm_t_tau_varying'
114 0 : return
115 : end if
116 :
117 10 : dlnT_dL = 0._dp
118 10 : dlnT_dlnR = 0._dp
119 10 : dlnT_dlnM = 0._dp
120 10 : dlnT_dlnkap = 0._dp
121 :
122 10 : dlnP_dL = 0._dp
123 10 : dlnP_dlnR = 0._dp
124 10 : dlnP_dlnM = 0._dp
125 10 : dlnP_dlnkap = 0._dp
126 :
127 : else
128 :
129 : ! Partials required, use finite differencing in Teff (we could
130 : ! in principle get dlnT_dlnTeff from the T-tau relation, but
131 : ! for consistency with dln_dlnTeff we use the same finite
132 : ! differencing)
133 :
134 0 : lnTeff = log(Teff)
135 :
136 0 : lnTeff_p = lnTeff + DLNTEFF
137 0 : lnTeff_m = lnTeff - DLNTEFF
138 :
139 : !$OMP SECTIONS
140 :
141 : !$OMP SECTION
142 :
143 : call eval_data( &
144 : tau_surf, exp(lnTeff), g, &
145 : T_tau_id, eos_proc, kap_proc, errtol, max_steps, &
146 0 : lnT, lnP, ierr)
147 :
148 : !$OMP SECTION
149 :
150 : call eval_data( &
151 : tau_surf, exp(lnTeff_p), g, &
152 : T_tau_id, eos_proc, kap_proc, errtol, max_steps, &
153 0 : lnT_p, lnP_p, ierr_p)
154 :
155 : !$OMP SECTION
156 :
157 : call eval_data( &
158 : tau_surf, exp(lnTeff_m), g, &
159 : T_tau_id, eos_proc, kap_proc, errtol, max_steps, &
160 0 : lnT_m, lnP_m, ierr_m)
161 :
162 : !$OMP END SECTIONS
163 :
164 0 : if (ierr_p /= 0) ierr = ierr_p
165 0 : if (ierr_m /= 0) ierr = ierr_m
166 0 : if (ierr /= 0) then
167 0 : write(*,*) 'Call to eval_data failed in atm_t_tau_varying_opacity'
168 0 : return
169 : end if
170 :
171 : ! Set up the partials
172 :
173 0 : dlnTeff_dlnR = -0.5_dp
174 0 : dlnTeff_dL = 0.25_dp/L
175 :
176 0 : dlnT_dlnTeff = (lnT_p - lnT_m) / (lnTeff_p - lnTeff_m)
177 0 : dlnT_dL = dlnT_dlnTeff*dlnTeff_dL
178 0 : dlnT_dlnR = dlnT_dlnTeff*dlnTeff_dlnR
179 0 : dlnT_dlnM = 0._dp
180 0 : dlnT_dlnkap = 0._dp
181 :
182 0 : dlnP_dlnTeff = (lnP_p - lnP_m) / (lnTeff_p - lnTeff_m)
183 0 : dlnP_dL = dlnP_dlnTeff*dlnTeff_dL
184 0 : dlnP_dlnR = dlnP_dlnTeff*dlnTeff_dlnR
185 0 : dlnP_dlnM = 0._dp
186 0 : dlnP_dlnkap = 0._dp
187 :
188 : end if
189 :
190 : return
191 :
192 : end subroutine eval_T_tau_varying
193 :
194 :
195 : ! Evaluate atmosphere data from T-tau relation with varying opacity
196 :
197 10 : subroutine eval_data( &
198 : tau_surf, Teff, g, &
199 : T_tau_id, eos_proc, kap_proc, errtol, max_steps, &
200 : lnT, lnP, ierr)
201 :
202 : use atm_def, only: atm_eos_iface, atm_kap_iface
203 :
204 : real(dp), intent(in) :: tau_surf
205 : real(dp), intent(in) :: Teff
206 : real(dp), intent(in) :: g
207 : integer, intent(in) :: T_tau_id
208 : procedure(atm_eos_iface) :: eos_proc
209 : procedure(atm_kap_iface) :: kap_proc
210 : real(dp), intent(in) :: errtol
211 : integer, intent(in) :: max_steps
212 : real(dp), intent(out) :: lnT
213 : real(dp), intent(out) :: lnP
214 : integer, intent(out) :: ierr
215 :
216 : real(dp), parameter :: TAU_OUTER_FACTOR = 1.E-5_dp
217 : integer, parameter :: MAX_TRIES = 3
218 : real(dp), parameter :: TRY_SCALE = 10._dp
219 :
220 : real(dp) :: tau_outer_curr
221 : real(dp) :: errtol_curr
222 : integer :: try
223 :
224 : ! Try doing the integration, relaxing tau_outer and errtol if failures occur
225 :
226 10 : tau_outer_curr = tau_surf*TAU_OUTER_FACTOR
227 10 : errtol_curr = errtol
228 :
229 10 : try_loop : do try = 1, MAX_TRIES
230 :
231 : call eval_data_try( &
232 : tau_surf, Teff, g, tau_outer_curr, &
233 : T_tau_id, eos_proc, kap_proc, errtol_curr, max_steps, &
234 10 : lnT, lnP, ierr)
235 10 : if (ierr == 0) exit try_loop
236 :
237 0 : tau_outer_curr = tau_outer_curr*TRY_SCALE
238 20 : errtol_curr = errtol*TRY_SCALE
239 :
240 : end do try_loop
241 :
242 10 : return
243 :
244 : end subroutine eval_data
245 :
246 :
247 10 : subroutine eval_data_try( &
248 : tau_surf, Teff, g, tau_outer, &
249 : T_tau_id, eos_proc, kap_proc, errtol, max_steps, &
250 : lnT, lnP, ierr)
251 :
252 : use eos_lib, only: radiation_pressure
253 : use atm_def, only: atm_eos_iface, atm_kap_iface
254 : use atm_T_tau_relations, only: eval_T_tau
255 : use num_lib, only: dopri5_work_sizes, dopri5
256 :
257 : real(dp), intent(in) :: tau_surf
258 : real(dp), intent(in) :: Teff
259 : real(dp), intent(in) :: g
260 : real(dp), intent(in) :: tau_outer
261 : integer, intent(in) :: T_tau_id
262 : procedure(atm_eos_iface) :: eos_proc
263 : procedure(atm_kap_iface) :: kap_proc
264 : real(dp), intent(in) :: errtol
265 : integer, intent(in) :: max_steps
266 : real(dp), intent(out) :: lnT
267 : real(dp), intent(out) :: lnP
268 : integer, intent(out) :: ierr
269 :
270 : integer, parameter :: NUM_VARS = 1
271 : integer, parameter :: NRDENS = 0
272 : integer, parameter :: LRPAR = 0
273 : integer, parameter :: LIPAR = 0
274 : integer, parameter :: IOUT = 0
275 : integer, parameter :: LOUT = 0
276 :
277 : real(dp), parameter :: RHO_OUTER = 1.E-10_dp
278 : real(dp), parameter :: DLNTAU_SURF = 1.E-3_dp
279 : real(dp), parameter :: DLNTAU_MAX = 0._dp
280 :
281 : integer :: liwork
282 : integer :: lwork
283 10 : real(dp), pointer :: work(:)
284 10 : integer, pointer :: iwork(:)
285 : real(dp), target :: rpar_ary(LRPAR)
286 : integer, target :: ipar_ary(LIPAR)
287 20 : real(dp), target :: y_ary(NUM_VARS)
288 10 : real(dp), pointer :: rpar(:)
289 10 : integer, pointer :: ipar(:)
290 10 : real(dp), pointer :: y(:)
291 10 : real(dp) :: lnTeff
292 10 : real(dp) :: T_outer
293 10 : real(dp) :: Pgas_outer
294 10 : real(dp) :: P_outer
295 10 : real(dp) :: lntau_outer
296 10 : real(dp) :: lntau_surf
297 10 : real(dp) :: dlntau
298 20 : real(dp) :: rtol(NUM_VARS)
299 20 : real(dp) :: atol(NUM_VARS)
300 : integer :: idid
301 :
302 10 : ierr = 0
303 :
304 : ! Allocate work arrays for the integrator
305 :
306 10 : call dopri5_work_sizes(NUM_VARS, NRDENS, liwork, lwork)
307 10 : allocate(work(lwork), iwork(liwork), stat=ierr)
308 10 : if (ierr /= 0) then
309 0 : write(*,*) 'allocate failed in eval_data_try'
310 0 : return
311 : end if
312 :
313 330 : work = 0._dp
314 220 : iwork = 0
315 :
316 : ! Set pointers (simply because dopri5 wants pointer args)
317 :
318 10 : rpar => rpar_ary
319 10 : ipar => ipar_ary
320 :
321 10 : y => y_ary
322 :
323 : ! Set values at the outer boundary (y(1) = lnP; estimate
324 : ! Pgas from low-density ideal gas law)
325 :
326 10 : lnTeff = log(Teff)
327 :
328 10 : call eval_T_tau(T_tau_id, tau_outer, Teff, lnT, ierr)
329 10 : if (ierr /= 0) then
330 0 : write(*,*) 'atm: Call to eval_T_tau failed in eval_data_try'
331 0 : return
332 : end if
333 :
334 10 : T_outer = exp(lnT)
335 10 : Pgas_outer = cgas*RHO_OUTER*T_outer
336 10 : P_outer = Pgas_outer + radiation_pressure(T_outer)
337 10 : lnP = log(P_outer)
338 :
339 10 : y(1) = lnP
340 :
341 : ! Integrate inward from tau_outer to tau_surf
342 :
343 10 : lntau_outer = log(tau_outer)
344 10 : lntau_surf = log(tau_surf)
345 :
346 10 : dlntau = DLNTAU_SURF
347 :
348 20 : rtol = errtol
349 20 : atol = errtol
350 :
351 : call dopri5( &
352 : NUM_VARS, eval_fcn, lntau_outer, y, lntau_surf, &
353 : dlntau, DLNTAU_MAX, max_steps, &
354 : rtol, atol, 1, &
355 : eval_solout, IOUT, &
356 : work, lwork, iwork, liwork, &
357 : LRPAR, rpar, LIPAR, ipar, &
358 10 : LOUT, idid)
359 10 : if (idid < 0) then
360 0 : write(*,*) 'Call to dopri5 failed in eval_data_try: idid=', idid
361 0 : ierr = -1
362 0 : return
363 : end if
364 :
365 : ! Store the final pressure and temperature
366 :
367 10 : lnP = y(1)
368 :
369 10 : call eval_T_tau(T_tau_id, tau_surf, Teff, lnT, ierr)
370 10 : if (ierr /= 0) then
371 0 : write(*,*) 'atm: Call to eval_T_tau failed in eval_data_try'
372 0 : return
373 : end if
374 :
375 : ! Deallocate arrays
376 :
377 10 : deallocate(work, iwork)
378 :
379 40 : return
380 :
381 : contains
382 :
383 6178 : subroutine eval_fcn(n, x, h, y, f, lr, rpar, li, ipar, ierr)
384 :
385 10 : use eos_def, only: num_eos_basic_results
386 : use atm_T_tau_relations, only: eval_T_tau
387 :
388 : integer, intent(in) :: n, lr, li
389 : real(dp), intent(in) :: x, h
390 : real(dp), intent(inout) :: y(:)
391 : real(dp), intent(inout) :: f(:)
392 : integer, intent(inout), pointer :: ipar(:)
393 : real(dp), intent(inout), pointer :: rpar(:)
394 : integer, intent(out) :: ierr
395 :
396 : real(dp) :: tau
397 6178 : real(dp) :: lnP
398 : real(dp) :: lnT
399 6178 : real(dp) :: lnRho
400 166806 : real(dp) :: res(num_eos_basic_results)
401 166806 : real(dp) :: dres_dlnRho(num_eos_basic_results)
402 166806 : real(dp) :: dres_dlnT(num_eos_basic_results)
403 6178 : real(dp) :: kap
404 6178 : real(dp) :: dlnkap_dlnRho
405 6178 : real(dp) :: dlnkap_dlnT
406 :
407 6178 : ierr = 0
408 :
409 : ! Calculate the temperature
410 :
411 6178 : tau = exp(x)
412 :
413 6178 : lnP = y(1)
414 :
415 6178 : call eval_T_tau(T_tau_id, tau, Teff, lnT, ierr)
416 6178 : if (ierr /= 0) then
417 0 : write(*,*) 'atm: Call to eval_T_tau failed in eval_fcn'
418 0 : return
419 : end if
420 :
421 : ! Calculate the density and eos results
422 :
423 : call eos_proc( &
424 : lnP, lnT, &
425 : lnRho, res, dres_dlnRho, dres_dlnT, &
426 6178 : ierr)
427 6178 : if (ierr /= 0) then
428 0 : write(*,*) 'atm: Call to eos_proc failed in eval_fcn'
429 0 : return
430 : end if
431 :
432 : ! Calculate the opacity
433 :
434 : call kap_proc( &
435 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
436 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
437 6178 : ierr)
438 6178 : if (ierr /= 0) then
439 0 : write(*,*) 'atm: Call to kap_proc failed in eval_fcn'
440 0 : return
441 : end if
442 :
443 : ! Set up the rhs for the hydrostatic eqm equation
444 : ! dlnP/dlntau = tau*g/P*kappa
445 :
446 6178 : f(1) = tau*g/(kap*exp(lnP))
447 :
448 6178 : return
449 :
450 6178 : end subroutine eval_fcn
451 :
452 :
453 0 : subroutine eval_solout( &
454 : nr, xold, x, n, y, rwork_y, iwork_y, interp_y, lrpar, rpar, lipar, ipar, irtrn)
455 :
456 : integer, intent(in) :: nr, n, lrpar, lipar
457 : real(dp), intent(in) :: xold, x
458 : real(dp), intent(inout) :: y(:)
459 : real(dp), intent(inout), target :: rwork_y(*)
460 : integer, intent(inout), target :: iwork_y(*)
461 : integer, intent(inout), pointer :: ipar(:)
462 : real(dp), intent(inout), pointer :: rpar(:)
463 : interface
464 : real(dp) function interp_y(i, s, rwork_y, iwork_y, ierr)
465 : use const_def, only: dp
466 : implicit none
467 : integer, intent(in) :: i
468 : real(dp), intent(in) :: s
469 : real(dp), intent(inout), target :: rwork_y(*)
470 : integer, intent(inout), target :: iwork_y(*)
471 : integer, intent(out) :: ierr
472 : end function interp_y
473 : end interface
474 : integer, intent(out) :: irtrn ! < 0 causes solver to return to calling program.
475 :
476 0 : irtrn = 0 ! for ifort
477 :
478 : ! Dummy routine that's never called
479 :
480 0 : call mesa_error(__FILE__,__LINE__,'Bogus call to eval_solout')
481 :
482 6178 : end subroutine eval_solout
483 :
484 : end subroutine eval_data_try
485 :
486 :
487 : ! Build atmosphere structure data from T-tau relation with varying
488 : ! opacity
489 :
490 0 : subroutine build_T_tau_varying( &
491 : tau_surf, L, R, Teff, M, cgrav, lnP_surf, tau_outer, &
492 : T_tau_id, eos_proc, kap_proc, errtol, dlogtau, &
493 : atm_structure_num_pts, atm_structure, &
494 : ierr)
495 :
496 : use atm_def
497 : use atm_utils, only: eval_Teff_g
498 : use num_lib, only: dopri5_work_sizes, dopri5
499 :
500 : real(dp), intent(in) :: tau_surf
501 : real(dp), intent(in) :: L
502 : real(dp), intent(in) :: R
503 : real(dp), intent(in) :: Teff
504 : real(dp), intent(in) :: M
505 : real(dp), intent(in) :: cgrav
506 : real(dp), intent(in) :: lnP_surf
507 : real(dp), intent(in) :: tau_outer
508 : integer, intent(in) :: T_tau_id
509 : procedure(atm_eos_iface) :: eos_proc
510 : procedure(atm_kap_iface) :: kap_proc
511 : real(dp), intent(in) :: errtol
512 : real(dp), intent(in) :: dlogtau
513 : integer, intent(out) :: atm_structure_num_pts
514 : real(dp), pointer :: atm_structure(:,:)
515 : integer, intent(out) :: ierr
516 :
517 : integer, parameter :: INIT_NUM_PTS = 100
518 : integer, parameter :: NUM_VARS = 2
519 : integer, parameter :: NRDENS = 0
520 : integer, parameter :: MAX_STEPS = 0
521 : integer, parameter :: LRPAR = 0
522 : integer, parameter :: LIPAR = 0
523 : integer, parameter :: IOUT = 1
524 : integer, parameter :: LOUT = 0
525 :
526 0 : real(dp) :: g
527 : integer :: liwork
528 : integer :: lwork
529 0 : real(dp), pointer :: work(:)
530 0 : integer, pointer :: iwork(:)
531 : real(dp), target :: rpar_ary(LRPAR)
532 : integer, target :: ipar_ary(LIPAR)
533 0 : real(dp), target :: y_ary(NUM_VARS)
534 0 : real(dp), pointer :: rpar(:)
535 0 : integer, pointer :: ipar(:)
536 0 : real(dp), pointer :: y(:)
537 0 : real(dp) :: lntau_surf
538 0 : real(dp) :: lntau_outer
539 0 : real(dp) :: rtol(NUM_VARS)
540 0 : real(dp) :: atol(NUM_VARS)
541 0 : real(dp) :: dlntau
542 0 : real(dp) :: dlntau_max
543 : integer :: idid
544 :
545 0 : ierr = 0
546 :
547 : ! Sanity checks
548 :
549 0 : if (L <= 0._dp .OR. R <= 0._dp .OR. M <= 0._dp) then
550 0 : ierr = -1
551 0 : return
552 : end if
553 :
554 0 : if (dlntau <= 0._dp) then
555 0 : write(*,*) 'Invalid dlntau in build_atm_uniform:', dlntau
556 0 : call mesa_error(__FILE__,__LINE__)
557 : end if
558 :
559 : ! Evaluate the gravity
560 :
561 0 : g = cgrav*M/(R*R)
562 :
563 : ! Allocate atm_structure at its initial size
564 :
565 0 : allocate(atm_structure(num_results_for_build_atm,INIT_NUM_PTS))
566 :
567 0 : atm_structure_num_pts = 0
568 :
569 0 : if (tau_outer > tau_surf) return
570 :
571 : ! Allocate work rrays for the integrator
572 :
573 0 : call dopri5_work_sizes(NUM_VARS, NRDENS, liwork, lwork)
574 0 : allocate(work(lwork), iwork(liwork), stat=ierr)
575 0 : if (ierr /= 0) then
576 0 : write(*,*) 'atm: allocate failed in build_T_tau_varying'
577 0 : deallocate(atm_structure)
578 0 : return
579 : end if
580 :
581 0 : work = 0._dp
582 0 : iwork = 0
583 :
584 : ! Set pointers (simply because dopri5 wants pointer args)
585 :
586 0 : rpar => rpar_ary
587 0 : ipar => ipar_ary
588 :
589 0 : y => y_ary
590 :
591 : ! Set starting values for the integrator (y(1) = delta_r; y(2) = lnP)
592 :
593 0 : y(1) = 0._dp
594 0 : y(2) = lnP_surf
595 :
596 : ! Integrate from the atmosphere base outward
597 :
598 0 : lntau_outer = log(tau_outer)
599 0 : lntau_surf = log(tau_surf)
600 :
601 0 : dlntau_max = -dlogtau*ln10
602 0 : dlntau = 0.5_dp*dlntau_max
603 :
604 0 : rtol = errtol
605 0 : atol = errtol
606 :
607 : call dopri5( &
608 : NUM_VARS, build_fcn, lntau_surf, y, lntau_outer, &
609 : dlntau, dlntau_max, MAX_STEPS, &
610 : rtol, atol, 1, &
611 : build_solout, IOUT, &
612 : work, lwork, iwork, liwork, &
613 : LRPAR, rpar, LIPAR, ipar, &
614 0 : LOUT, idid)
615 0 : if (idid < 0) then
616 0 : write(*,*) 'atm: Call to dopri5 failed in build_T_tau_varying: idid=', idid
617 0 : ierr = -1
618 : end if
619 :
620 : ! Reverse the data ordering (since by convention data in atm_structure
621 : ! should be ordered outward-in)
622 :
623 : atm_structure(:,:atm_structure_num_pts) = &
624 0 : atm_structure(:,atm_structure_num_pts:1:-1)
625 :
626 : ! Deallocate arrays
627 :
628 0 : deallocate(work, iwork)
629 :
630 : contains
631 :
632 0 : subroutine build_fcn(n, x, h, y, f, lr, rpar, li, ipar, ierr)
633 :
634 : integer, intent(in) :: n, lr, li
635 : real(dp), intent(in) :: x, h
636 : real(dp), intent(inout) :: y(:)
637 : real(dp), intent(inout) :: f(:)
638 : integer, intent(inout), pointer :: ipar(:)
639 : real(dp), intent(inout), pointer :: rpar(:)
640 : integer, intent(out) :: ierr
641 :
642 0 : real(dp) :: tau
643 0 : real(dp) :: P
644 0 : real(dp) :: rho
645 0 : real(dp) :: kap
646 :
647 0 : real(dp) :: atm_structure_sgl(num_results_for_build_atm)
648 :
649 : ierr = 0
650 :
651 : ! Evaluate structure data
652 :
653 0 : call build_data(x, y(1), y(2), atm_structure_sgl, ierr)
654 0 : if (ierr /= 0) then
655 : ! Non-zero ierr signals we must use a smaller stepsize;
656 : ! whereas in fact we want to terminate the integration.
657 : ! Needs fixing!
658 0 : return
659 : end if
660 :
661 : ! Set up the rhs for the optical depth and hydrostatic
662 : ! equilibrium equations
663 : ! dr/dlntau = -tau/(kappa*rho)
664 : !
665 :
666 0 : tau = exp(x)
667 0 : P = exp(y(2))
668 0 : rho = exp(atm_structure_sgl(atm_lnd))
669 0 : kap = atm_structure_sgl(atm_kap)
670 :
671 0 : f(1) = -tau/(kap*rho)
672 0 : f(2) = tau*g/(kap*P)
673 :
674 0 : end subroutine build_fcn
675 :
676 :
677 0 : subroutine build_solout( &
678 0 : nr, xold, x, n, y, rwork_y, iwork_y, interp_y, lrpar, rpar, lipar, ipar, irtrn)
679 :
680 : use utils_lib, only: realloc_double2
681 :
682 : integer, intent(in) :: nr, n, lrpar, lipar
683 : real(dp), intent(in) :: xold, x
684 : real(dp), intent(inout) :: y(:)
685 : real(dp), intent(inout), target :: rwork_y(*)
686 : integer, intent(inout), target :: iwork_y(*)
687 : integer, intent(inout), pointer :: ipar(:)
688 : real(dp), intent(inout), pointer :: rpar(:)
689 : interface
690 : real(dp) function interp_y(i, s, rwork_y, iwork_y, ierr)
691 : use const_def, only: dp
692 : implicit none
693 : integer, intent(in) :: i
694 : real(dp), intent(in) :: s
695 : real(dp), intent(inout), target :: rwork_y(*)
696 : integer, intent(inout), target :: iwork_y(*)
697 : integer, intent(out) :: ierr
698 : end function interp_y
699 : end interface
700 : integer, intent(out) :: irtrn
701 :
702 : integer :: ierr
703 : integer :: sz
704 0 : real(dp) :: atm_structure_sgl(num_results_for_build_atm)
705 :
706 : ierr = 0
707 :
708 : ! Evaluate structure data
709 :
710 0 : call build_data(x, y(1), y(2), atm_structure_sgl, ierr)
711 0 : if (ierr /= 0) then
712 0 : irtrn = -1
713 0 : return
714 : end if
715 :
716 : ! If necessary, expand arrays
717 :
718 0 : atm_structure_num_pts = atm_structure_num_pts + 1
719 0 : sz = size(atm_structure, dim=2)
720 :
721 0 : if (atm_structure_num_pts > sz) then
722 0 : sz = 2*sz + 100
723 : call realloc_double2( &
724 0 : atm_structure,num_results_for_build_atm,sz,ierr)
725 : end if
726 :
727 : ! Store data
728 :
729 0 : atm_structure(:,atm_structure_num_pts) = atm_structure_sgl
730 :
731 : return
732 :
733 : end subroutine build_solout
734 :
735 :
736 0 : subroutine build_data(lntau, delta_r, lnP, atm_structure_sgl, ierr)
737 :
738 : use eos_def
739 : use eos_lib, only: radiation_pressure
740 : use atm_T_tau_relations, only: eval_T_tau
741 : use atm_utils, only: eval_Paczynski_gradr
742 :
743 : real(dp), intent(in) :: lntau
744 : real(dp), intent(in) :: delta_r
745 : real(dp), intent(in) :: lnP
746 : real(dp), intent(out) :: atm_structure_sgl(:)
747 : integer, intent(out) :: ierr
748 :
749 0 : real(dp) :: tau
750 0 : real(dp) :: lnT
751 0 : real(dp) :: lnRho
752 0 : real(dp) :: res(num_eos_basic_results)
753 0 : real(dp) :: dres_dlnRho(num_eos_basic_results)
754 0 : real(dp) :: dres_dlnT(num_eos_basic_results)
755 0 : real(dp) :: kap
756 0 : real(dp) :: dlnkap_dlnRho
757 0 : real(dp) :: dlnkap_dlnT
758 0 : real(dp) :: gradr
759 :
760 0 : ierr = 0
761 :
762 : ! Calculate the temperature
763 :
764 0 : tau = exp(lntau)
765 :
766 0 : call eval_T_tau(T_tau_id, tau, Teff, lnT, ierr)
767 0 : if (ierr /= 0) then
768 0 : write(*,*) 'atm: Call to eval_T_tau failed in build_data'
769 0 : return
770 : end if
771 :
772 : ! Calculate the density and eos results
773 :
774 : call eos_proc( &
775 : lnP, lnT, &
776 : lnRho, res, dres_dlnRho, dres_dlnT, &
777 0 : ierr)
778 0 : if (ierr /= 0) then
779 0 : write(*,*) 'atm: Call to eos_proc failed in build_data'
780 0 : return
781 : end if
782 :
783 : ! Calculate the opacity
784 :
785 : call kap_proc( &
786 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
787 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
788 0 : ierr)
789 0 : if (ierr /= 0) then
790 0 : write(*,*) 'atm: Call to kap_proc failed in build_data'
791 0 : return
792 : end if
793 :
794 : ! Evaluate radiative temperature gradient
795 :
796 0 : gradr = eval_Paczynski_gradr(exp(lnT), exp(lnP), exp(lnRho), tau, kap, L, M, R, cgrav)
797 :
798 : ! Store data
799 :
800 0 : atm_structure_sgl(atm_xm) = 0._dp ! We assume negligible mass in the atmosphere
801 0 : atm_structure_sgl(atm_delta_r) = delta_r
802 0 : atm_structure_sgl(atm_lnP) = lnP
803 0 : atm_structure_sgl(atm_lnd) = lnRho
804 0 : atm_structure_sgl(atm_lnT) = lnT
805 0 : atm_structure_sgl(atm_gradT) = gradr ! by assumption, atm is radiative
806 0 : atm_structure_sgl(atm_kap) = kap
807 0 : atm_structure_sgl(atm_gamma1) = res(i_gamma1)
808 0 : atm_structure_sgl(atm_grada) = res(i_grad_ad)
809 0 : atm_structure_sgl(atm_chiT) = res(i_chiT)
810 0 : atm_structure_sgl(atm_chiRho) = res(i_chiRho)
811 0 : atm_structure_sgl(atm_cp) = res(i_Cp)
812 0 : atm_structure_sgl(atm_cv) = res(i_Cv)
813 0 : atm_structure_sgl(atm_tau) = tau
814 0 : atm_structure_sgl(atm_lnfree_e) = res(i_lnfree_e)
815 0 : atm_structure_sgl(atm_dlnkap_dlnT) = dlnkap_dlnT
816 0 : atm_structure_sgl(atm_dlnkap_dlnd) = dlnkap_dlnRho
817 0 : atm_structure_sgl(atm_lnPgas) = res(i_lnPgas)
818 0 : atm_structure_sgl(atm_gradr) = gradr
819 :
820 0 : end subroutine build_data
821 :
822 : end subroutine build_T_tau_varying
823 :
824 : end module atm_T_tau_varying
|
