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_table
21 :
22 : use const_def, only: dp, ln10
23 : use math_lib
24 : use utils_lib, only: mesa_error
25 :
26 : implicit none
27 :
28 : private
29 : public :: eval_table
30 : public :: get_table_alfa_beta
31 : public :: get_table_base
32 :
33 : contains
34 :
35 14 : subroutine eval_table( &
36 : L, R, M, cgrav, id, Z, skip_partials, &
37 : Teff, &
38 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
39 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
40 : ierr)
41 :
42 : use atm_def
43 : use eos_lib, only: radiation_pressure
44 : use table_atm, only: get_table_values
45 : use utils_lib, only: is_bad
46 :
47 : real(dp), intent(in) :: L
48 : real(dp), intent(in) :: R
49 : real(dp), intent(in) :: M
50 : real(dp), intent(in) :: cgrav
51 : integer, intent(in) :: id
52 : real(dp), intent(in) :: Z
53 : logical, intent(in) :: skip_partials
54 : real(dp), intent(in) :: Teff
55 : real(dp), intent(out) :: lnT
56 : real(dp), intent(out) :: dlnT_dL
57 : real(dp), intent(out) :: dlnT_dlnR
58 : real(dp), intent(out) :: dlnT_dlnM
59 : real(dp), intent(out) :: dlnT_dlnkap
60 : real(dp), intent(out) :: lnP
61 : real(dp), intent(out) :: dlnP_dL
62 : real(dp), intent(out) :: dlnP_dlnR
63 : real(dp), intent(out) :: dlnP_dlnM
64 : real(dp), intent(out) :: dlnP_dlnkap
65 : integer, intent(out) :: ierr
66 :
67 : logical, parameter :: dbg = .false.
68 :
69 : real(dp) :: g
70 14 : real(dp) :: logg
71 14 : real(dp) :: Pgas
72 14 : real(dp) :: dPgas_dTeff
73 14 : real(dp) :: dPgas_dlogg
74 14 : real(dp) :: T
75 14 : real(dp) :: dT_dTeff
76 14 : real(dp) :: dT_dlogg
77 14 : real(dp) :: Prad
78 14 : real(dp) :: P
79 14 : real(dp) :: dlnTeff_dL
80 14 : real(dp) :: dlnTeff_dlnR
81 14 : real(dp) :: dTeff_dlnR
82 14 : real(dp) :: dlogg_dlnR
83 14 : real(dp) :: dlogg_dlnM
84 14 : real(dp) :: dPgas_dlnR
85 14 : real(dp) :: dPgas_dlnM
86 14 : real(dp) :: dPgas_dlnTeff
87 14 : real(dp) :: dPrad_dlnTeff
88 14 : real(dp) :: dPrad_dlnR
89 14 : real(dp) :: dPrad_dL
90 14 : real(dp) :: dT_dlnR
91 14 : real(dp) :: dT_dlnM
92 14 : real(dp) :: dT_dlnTeff
93 :
94 : include 'formats'
95 :
96 14 : ierr = 0
97 :
98 : ! Sanity checks
99 :
100 14 : if (L <= 0._dp .OR. R <= 0._dp .OR. M <= 0._dp) then
101 0 : ierr = -1
102 0 : return
103 : end if
104 :
105 : ! Evaluate the gravity
106 :
107 14 : g = cgrav*M/(R*R)
108 14 : logg = log10(g)
109 :
110 : ! Perform the table lookup
111 :
112 : if (dbg) write(*,*) 'call get_table_values', id
113 : call get_table_values( &
114 : id, Z, logg, Teff, &
115 : Pgas, dPgas_dTeff, dPgas_dlogg, &
116 : T, dT_dTeff, dT_dlogg, &
117 14 : ierr)
118 14 : if (ierr /= 0) then
119 : if (dbg) write(*,*) 'get_table_values(_at_fixed_Z) ierr', ierr
120 : return
121 : end if
122 :
123 : ! Set up lnT and lnP
124 :
125 14 : lnT = log(T)
126 :
127 14 : Prad = radiation_pressure(T)
128 14 : P = Pgas + Prad
129 14 : lnP = log(P)
130 :
131 : ! Set up partials
132 :
133 14 : if (.NOT. skip_partials) then
134 :
135 0 : dlnTeff_dlnR = -0.5_dp
136 0 : dlnTeff_dL = 0.25_dp/L
137 0 : dTeff_dlnR = Teff*dlnTeff_dlnR
138 :
139 0 : dlogg_dlnR = -2._dp/ln10
140 0 : dlogg_dlnM = 1._dp/ln10
141 :
142 0 : dPgas_dlnR = dPgas_dlogg*dlogg_dlnR + dPgas_dTeff*dTeff_dlnR
143 0 : dPgas_dlnM = dPgas_dlogg*dlogg_dlnM
144 0 : dPgas_dlnTeff = dPgas_dTeff*Teff
145 :
146 0 : dPrad_dlnTeff = 4._dp*Prad*dT_dTeff*Teff/T
147 0 : dPrad_dlnR = dPrad_dlnTeff*dlnTeff_dlnR
148 0 : dPrad_dL = dPrad_dlnTeff*dlnTeff_dL
149 :
150 0 : dlnP_dL = (dPgas_dlnTeff + dPrad_dlnTeff)*dlnTeff_dL / P
151 0 : dlnP_dlnR = (dPgas_dlnR + dPrad_dlnTeff*dlnTeff_dlnR) / P
152 0 : dlnP_dlnM = dPgas_dlnM/P
153 0 : dlnP_dlnkap = 0._dp
154 :
155 0 : dT_dlnTeff = dT_dTeff*Teff
156 0 : dT_dlnR = dT_dlogg*dlogg_dlnR + dT_dlnTeff*dlnTeff_dlnR
157 0 : dT_dlnM = dT_dlogg*dlogg_dlnM
158 :
159 0 : dlnT_dL = dT_dlnTeff*dlnTeff_dL / T
160 0 : dlnT_dlnR = dT_dlnR / T
161 0 : dlnT_dlnM = dT_dlnM / T
162 0 : dlnT_dlnkap = 0._dp
163 :
164 : else
165 :
166 14 : dlnP_dL = 0._dp
167 14 : dlnP_dlnR = 0._dp
168 14 : dlnP_dlnM = 0._dp
169 14 : dlnP_dlnkap = 0._dp
170 :
171 14 : dlnT_dL = 0._dp
172 14 : dlnT_dlnR = 0._dp
173 14 : dlnT_dlnM = 0._dp
174 14 : dlnT_dlnkap = 0._dp
175 :
176 : end if
177 :
178 14 : if (dbg .or. is_bad(lnP) .or. is_bad(lnT)) then
179 0 : write(*,*) 'eval_table'
180 0 : write(*,1) 'Teff', Teff
181 0 : write(*,1) 'T', T
182 0 : write(*,1) 'dT_dTeff', dT_dTeff
183 0 : write(*,1) 'dT_dlogg', dT_dlogg
184 0 : write(*,'(A)')
185 0 : ierr = -1
186 0 : return
187 : !if (is_bad(lnP) .or. is_bad(lnT)) call mesa_error(__FILE__,__LINE__,'eval_table')
188 : end if
189 :
190 : return
191 :
192 14 : end subroutine eval_table
193 :
194 :
195 0 : subroutine get_table_alfa_beta( &
196 : L, Teff, R, M, cgrav, id, alfa, beta, ierr)
197 :
198 14 : use atm_def
199 : use table_atm, only: &
200 : ai_two_thirds, ai_100, ai_10, ai_1, ai_1m1, ai_wd_25, ai_db_wd_25
201 :
202 : real(dp), intent(in) :: L
203 : real(dp), intent(in) :: Teff
204 : real(dp), intent(in) :: R
205 : real(dp), intent(in) :: M
206 : real(dp), intent(in) :: cgrav
207 : integer, intent(in) :: id
208 : real(dp), intent(out) :: alfa
209 : real(dp), intent(out) :: beta
210 : integer, intent(out) :: ierr
211 :
212 : integer, parameter :: PURE_GREY = 1
213 : integer, parameter :: PURE_TABLE = 2
214 : integer, parameter :: BLEND_IN_X = 3
215 : integer, parameter :: BLEND_IN_Y = 4
216 : integer, parameter :: BLEND_CORNER_OUT = 5
217 :
218 : logical, parameter :: dbg = .false.
219 :
220 0 : real(dp) :: g
221 0 : real(dp) :: logTeff
222 0 : real(dp) :: logg
223 : type(Atm_Info), pointer :: ai
224 0 : real(dp) :: logg_max
225 0 : real(dp) :: logg_min
226 0 : real(dp) :: logTeff_max
227 0 : real(dp) :: logTeff_min
228 0 : real(dp) :: Teff_max
229 0 : real(dp) :: Teff_min
230 0 : real(dp) :: logg_min_margin
231 0 : real(dp) :: logg_max_margin
232 0 : real(dp) :: logTeff_min_margin
233 0 : real(dp) :: logTeff_max_margin
234 0 : real(dp) :: logg1
235 0 : real(dp) :: logg2
236 0 : real(dp) :: logg3
237 0 : real(dp) :: logg4
238 0 : real(dp) :: logTeff1
239 0 : real(dp) :: logTeff2
240 0 : real(dp) :: logTeff3
241 0 : real(dp) :: logTeff4
242 0 : real(dp) :: c_dx
243 0 : real(dp) :: c_dy
244 : integer :: iregion
245 : integer :: logg_index
246 : integer :: ng
247 : integer :: j
248 :
249 : include 'formats'
250 :
251 0 : ierr = 0
252 :
253 : ! Sanity checks
254 :
255 0 : if (L <= 0._dp .OR. R <= 0._dp .OR. M <= 0._dp) then
256 0 : ierr = -1
257 0 : return
258 : end if
259 :
260 : ! Evaluate the gravity
261 :
262 0 : g = cgrav*M/(R*R)
263 0 : logTeff = log10(Teff)
264 0 : logg = log10(g)
265 :
266 : ! Set up the table pointer
267 :
268 0 : select case (id)
269 : case (ATM_TABLE_PHOTOSPHERE)
270 : ai => ai_two_thirds
271 : case (ATM_TABLE_TAU_100)
272 : ai => ai_100
273 : case (ATM_TABLE_TAU_10)
274 : ai => ai_10
275 : case (ATM_TABLE_TAU_1)
276 : ai => ai_1
277 : case (ATM_TABLE_TAU_1M1)
278 : ai => ai_1m1
279 : case (ATM_TABLE_WD_TAU_25)
280 : ai => ai_wd_25
281 : case (ATM_TABLE_DB_WD_TAU_25)
282 0 : ai => ai_db_wd_25
283 : case default
284 0 : write(*,*) 'Invalid id in get_table_alfa_beta:', id
285 0 : call mesa_error(__FILE__,__LINE__)
286 : end select
287 :
288 0 : ng = ai% ng
289 0 : logg_max = ai% logg_array(ng)
290 0 : logg_min = ai% logg_array(1)
291 :
292 : ! First, locate current point logg array for use with Teff_bound
293 :
294 0 : if (logg <= logg_min) then
295 : logg_index = 1
296 0 : else if (logg >= logg_max) then
297 : logg_index = ng
298 : else
299 : logg_index = ng
300 0 : do j = 1, ng-1
301 0 : if (ai% logg_array(j) <= logg) then
302 0 : logg_index = j
303 : end if
304 : end do
305 : end if
306 :
307 0 : Teff_max = ai% Teff_bound(logg_index)
308 0 : Teff_min = ai% Teff_array(1)
309 0 : logTeff_max = log10(Teff_max)
310 0 : logTeff_min = log10(Teff_min)
311 :
312 : ! Set up margins for blending
313 :
314 0 : logg_max_margin = 0.01d0
315 0 : logg_min_margin = 0.5d0
316 0 : logTeff_max_margin = 0.01d0
317 0 : logTeff_min_margin = 0.01d0
318 :
319 : ! if (id == ATM_TABLE_WD_TAU_25) then
320 : ! logg_max_margin = 0.5d0
321 : ! logg_min_margin = 0.5d0
322 : ! logTeff_max_margin = 0.2d0
323 : ! logTeff_min_margin = 1d0
324 : ! logg1 = logg_max + logg_max_margin
325 : ! logg2 = logg_max
326 : ! logg3 = logg_min
327 : ! logg4 = logg_min - logg_min_margin
328 : ! logTeff1 = logTeff_max + logTeff_max_margin
329 : ! logTeff2 = logTeff_max
330 : ! logTeff3 = logTeff_min
331 : ! logTeff4 = logTeff_min - logTeff_min_margin
332 : ! else
333 0 : logg1 = logg_max
334 0 : logg2 = logg_max - logg_max_margin
335 0 : logg3 = logg_min + logg_min_margin
336 0 : logg4 = logg_min
337 0 : logTeff1 = logTeff_max
338 0 : logTeff2 = logTeff_max - logTeff_max_margin
339 0 : logTeff3 = logTeff_min
340 0 : logTeff4 = logTeff_min
341 : ! end if
342 :
343 : ! Decide on what sort of region we're in
344 :
345 0 : if (logg < logg4 .or. logg > logg1 .or. logTeff > logTeff1) then
346 0 : iregion = pure_grey
347 0 : else if (logTeff > logTeff2) then
348 0 : c_dy = (logTeff - logTeff2) / (logTeff1 - logTeff2)
349 0 : if (logg > logg2) then
350 0 : c_dx = (logg - logg2) / (logg1 - logg2)
351 0 : iregion = blend_corner_out
352 0 : else if (logg > logg3) then
353 0 : iregion = blend_in_y
354 : else ! logg > logg4
355 0 : c_dx = (logg - logg3) / (logg4 - logg3)
356 0 : iregion = blend_corner_out
357 : end if
358 0 : else if (logTeff >= logTeff3) then
359 0 : if (logg > logg2) then
360 0 : c_dx = (logg - logg2) / (logg1 - logg2)
361 0 : iregion = blend_in_x
362 0 : else if (logg > logg3) then
363 0 : iregion = pure_table
364 : else ! logg > logg4
365 0 : c_dx = (logg - logg3) / (logg4 - logg3)
366 0 : iregion = blend_in_x
367 : end if
368 0 : else if (logTeff > logTeff4) then
369 0 : c_dy = (logTeff - logTeff3) / (logTeff4 - logTeff3)
370 0 : if (logg > logg2) then
371 0 : c_dx = (logg - logg2) / (logg1 - logg2)
372 0 : iregion = blend_corner_out
373 0 : else if (logg > logg3) then
374 0 : iregion = blend_in_y
375 : else ! logg > logg4
376 0 : c_dx = (logg - logg3) / (logg4 - logg3)
377 0 : iregion = blend_corner_out
378 : end if
379 : else ! logTeff <= logTeff4
380 0 : iregion = pure_grey
381 : end if
382 :
383 : ! Set up alfa and beta
384 :
385 0 : select case (iregion)
386 : case (pure_grey)
387 0 : alfa = 1._dp
388 0 : beta = 0._dp
389 : case (pure_table)
390 0 : alfa = 0._dp
391 0 : beta = 1._dp
392 : case (blend_in_y)
393 0 : alfa = c_dy
394 0 : beta = 1._dp - alfa
395 : case (blend_in_x)
396 0 : alfa = c_dx
397 0 : beta = 1._dp - alfa
398 : case (blend_corner_out)
399 0 : alfa = min(1d0, sqrt(c_dx*c_dx + c_dy*c_dy))
400 0 : beta = 1 - alfa
401 : case default
402 0 : write(*,*) 'Invalid iregion in get_table_alfa_beta:', iregion
403 0 : call mesa_error(__FILE__,__LINE__)
404 : end select
405 :
406 : return
407 :
408 0 : end subroutine get_table_alfa_beta
409 :
410 :
411 0 : subroutine get_table_base (id, tau_base, ierr)
412 :
413 0 : use atm_def
414 :
415 : integer, intent(in) :: id
416 : real(dp), intent(out) :: tau_base
417 : integer, intent(out) :: ierr
418 :
419 0 : ierr = 0
420 :
421 : ! Get the base optical depth for the atmosphere
422 :
423 0 : select case (id)
424 : case (ATM_TABLE_PHOTOSPHERE)
425 0 : tau_base = 2._dp/3._dp
426 : case (ATM_TABLE_TAU_100)
427 0 : tau_base = 100._dp
428 : case (ATM_TABLE_TAU_10)
429 0 : tau_base = 10._dp
430 : case (ATM_TABLE_TAU_1)
431 0 : tau_base = 1._dp
432 : case (ATM_TABLE_TAU_1M1)
433 0 : tau_base = 0.1_dp
434 : case (ATM_TABLE_WD_TAU_25)
435 0 : tau_base = 25.1188_dp
436 : case (ATM_TABLE_DB_WD_TAU_25)
437 0 : tau_base = 25._dp
438 : case default
439 0 : write(*,*) 'Invalid id in get_table_base:', id
440 0 : call mesa_error(__FILE__,__LINE__)
441 : end select
442 :
443 0 : return
444 :
445 : end subroutine get_table_base
446 :
447 : end module atm_table
|
