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 table_atm
21 : use atm_def
22 : use const_def, only: dp
23 : use math_lib
24 : use utils_lib, only: mesa_error
25 :
26 : implicit none
27 :
28 : logical, parameter :: dbg = .false.
29 :
30 : contains
31 :
32 : !reads in table_summary file from atm_data, initializes logZ, Teff_array,
33 : ! logg_array, and Teff_bound arrays; sets some flags
34 3 : subroutine table_atm_init(use_cache, ierr)
35 : logical, intent(in) :: use_cache
36 : integer, intent(out) :: ierr
37 :
38 : integer :: nZ, ng, nT, i, j, iounit
39 3 : integer, pointer :: ibound(:,:), tmp_version(:)
40 : character(len=256) :: filename
41 :
42 3 : if (table_atm_is_initialized) call table_atm_shutdown()
43 :
44 : ierr = 0
45 :
46 3 : call load_table_summary(ATM_TABLE_PHOTOSPHERE, 'table_summary.txt', ai_two_thirds, ierr)
47 3 : if (ierr /= 0) return
48 3 : call load_table_summary(ATM_TABLE_TAU_100, 'table100_summary.txt', ai_100, ierr)
49 3 : if (ierr /= 0) return
50 3 : call load_table_summary(ATM_TABLE_TAU_10, 'table10_summary.txt', ai_10, ierr)
51 3 : if (ierr /= 0) return
52 3 : call load_table_summary(ATM_TABLE_TAU_1, 'table1_summary.txt', ai_1, ierr)
53 3 : if (ierr /= 0) return
54 3 : call load_table_summary(ATM_TABLE_TAU_1M1, 'table1m1_summary.txt', ai_1m1, ierr)
55 3 : if (ierr /= 0) return
56 3 : call load_table_summary(ATM_TABLE_WD_TAU_25, 'table_wd_25_summary.txt', ai_wd_25, ierr)
57 3 : if (ierr /= 0) return
58 3 : call load_table_summary(ATM_TABLE_DB_WD_TAU_25, 'table_db_wd_25_summary.txt', ai_db_wd_25, ierr)
59 3 : if (ierr /= 0) return
60 :
61 3 : table_atm_is_initialized = .true.
62 :
63 :
64 : contains
65 :
66 :
67 21 : subroutine load_table_summary(id, fname, ai, ierr)
68 : use const_def, only: mesa_data_dir
69 : integer, intent(in) :: id
70 : character(len=*), intent(in) :: fname
71 : type(atm_info), intent(inout) :: ai
72 : integer, intent(out) :: ierr
73 :
74 441 : real(dp), target :: vec_ary(20)
75 21 : real(dp), pointer :: vec(:)
76 :
77 21 : vec => vec_ary
78 :
79 21 : filename = trim(mesa_data_dir)//'/atm_data/' // trim(fname)
80 :
81 : if (dbg) write(*,*) 'read ' // trim(filename)
82 :
83 21 : open(newunit=iounit,file=trim(filename),action='read',status='old',iostat=ierr)
84 21 : if (ierr/= 0) then
85 0 : write(*,*) 'table_atm_init: missing atm data'
86 0 : write(*,*) trim(filename)
87 0 : write(*,'(A)')
88 0 : write(*,'(A)')
89 0 : write(*,'(A)')
90 0 : write(*,'(A)')
91 0 : write(*,'(A)')
92 0 : write(*,*) 'FATAL ERROR: missing or bad atm data.'
93 0 : write(*,*) 'Please update by removing the directory mesa/data/atm_data,'
94 0 : write(*,*) 'and rerunning the mesa ./install script.'
95 0 : write(*,'(A)')
96 0 : call mesa_error(__FILE__,__LINE__)
97 : end if
98 :
99 : !read first line and (nZ, nT, ng)
100 21 : read(iounit,*) !first line is text, skip it
101 21 : read(iounit,*) nZ, nT, ng
102 :
103 21 : ai% nZ = nZ
104 21 : ai% nT = nT
105 21 : ai% ng = ng
106 21 : ai% id = id
107 :
108 : allocate( &
109 : ai% Teff_array(nT), ai% logg_array(ng), ai% Teff_bound(ng), &
110 : ai% logZ(nZ), ai% alphaFe(nZ), &
111 : ai% Pgas_interp1(4*ng*nT*nZ), ai% T_interp1(4*ng*nT*nZ), &
112 21 : ai% have_atm_table(nZ), ai% atm_mix(nZ), ai% table_atm_files(nZ))
113 :
114 21 : ai% Pgas_interp(1:4,1:ng,1:nT,1:nZ) => ai% Pgas_interp1(1:4*ng*nT*nZ)
115 21 : ai% T_interp(1:4,1:ng,1:nT,1:nZ) => ai% T_interp1(1:4*ng*nT*nZ)
116 :
117 21 : allocate(ibound(ng,nZ), tmp_version(nZ))
118 :
119 63 : ai% have_atm_table(:) = .false.
120 :
121 : !read filenames and headers
122 21 : read(iounit,*) !text
123 63 : do i=1,nZ
124 42 : read(iounit,'(a)') ai% table_atm_files(i)
125 42 : read(iounit,'(14x,i4)') tmp_version(i)
126 63 : read(iounit,1) ai% logZ(i), ai% alphaFe(i), ai% atm_mix(i), ibound(1:ng,i)
127 : end do
128 :
129 : !read Teff_array
130 21 : read(iounit,*) !text
131 21 : read(iounit,2) ai% Teff_array(:)
132 :
133 : !read logg_array
134 21 : read(iounit,*) !text
135 21 : read(iounit,3) ai% logg_array(:)
136 :
137 21 : close(iounit)
138 :
139 : 1 format(13x,f5.2,8x,f4.1,1x,a8,1x,15x,99i4)
140 : 2 format(13f7.0)
141 : 3 format(13f7.2)
142 :
143 : !determine table boundaries
144 393 : do i=1,ng ! -- for each logg, smallest Teff at which Pgas->0
145 372 : ai% Teff_bound(i) = ai% Teff_array(ibound(i,1))
146 666 : do j=2,nZ
147 645 : ai% Teff_bound(i) = min( ai% Teff_bound(i) , ai% Teff_array(ibound(i,j)) )
148 : end do
149 : end do
150 :
151 :
152 : if (dbg) write(*,*) 'ai% logg_array(:)', ai% logg_array(:)
153 :
154 :
155 21 : deallocate(ibound, tmp_version)
156 :
157 21 : end subroutine load_table_summary
158 :
159 :
160 : end subroutine table_atm_init
161 :
162 :
163 1 : subroutine table_atm_shutdown()
164 :
165 1 : if (.NOT. table_atm_is_initialized) return
166 :
167 1 : call free_table_summary(ai_two_thirds)
168 1 : call free_table_summary(ai_100)
169 1 : call free_table_summary(ai_10)
170 1 : call free_table_summary(ai_1)
171 1 : call free_table_summary(ai_1m1)
172 1 : call free_table_summary(ai_wd_25)
173 1 : call free_table_summary(ai_db_wd_25)
174 :
175 1 : table_atm_is_initialized = .false.
176 :
177 : contains
178 :
179 7 : subroutine free_table_summary(ai)
180 :
181 : type(atm_info), intent(inout) :: ai
182 :
183 0 : deallocate( &
184 0 : ai% Teff_array, ai% logg_array, ai% Teff_bound, &
185 0 : ai% logZ, ai% alphaFe, &
186 0 : ai% Pgas_interp1, ai% T_interp1, &
187 7 : ai% have_atm_table, ai% atm_mix, ai% table_atm_files)
188 :
189 7 : end subroutine free_table_summary
190 :
191 : end subroutine table_atm_shutdown
192 :
193 :
194 : !interpolate in Z, logg, & Teff: 4pt in Z; bicubic spline in Teff,logg
195 14 : subroutine get_table_values( &
196 : id, newZ, newlogg_in, newTeff_in, &
197 : newPgas, dPgas_dTeff, dPgas_dlogg, &
198 : newT, dT_dTeff, dT_dlogg, &
199 : ierr)
200 : use chem_def,only:GN93_Zsol
201 : use interp_1d_lib, only: interp_pm, interp_values
202 : use interp_1d_def
203 : use interp_2d_lib_db, only: interp_evbicub_db
204 : use utils_lib, only: is_bad
205 :
206 :
207 : integer, intent(in) :: id
208 : real(dp), intent(in) :: newZ, newlogg_in, newTeff_in
209 : real(dp), intent(out) :: newPgas, dPgas_dTeff, dPgas_dlogg
210 : real(dp), intent(out) :: newT, dT_dTeff, dT_dlogg
211 : integer, intent(out) :: ierr
212 : integer :: i, j, numZs, nZ, ng, nT, ict(6), Zindex, Zlo, Zhi
213 : integer, parameter :: num_res = 3 !number of results (Pgas, dPgas_dTeff, dPgas_dlogg)
214 42 : real(dp) :: newlogg, newTeff, newlogZ(1), result_Z(1)
215 238 : real(dp), target :: fZ1_ary(4*4)
216 : real(dp), pointer :: fZ1(:), fZ(:,:)
217 798 : real(dp), target :: work_ary(4*mp_work_size)
218 14 : real(dp), pointer :: result_2D(:,:), work(:)
219 : character(len=256) :: filename
220 : logical :: clip_Teff, clip_logg, gtv_dbg
221 :
222 : type (Atm_Info), pointer :: ai
223 :
224 : include 'formats'
225 :
226 14 : gtv_dbg = dbg
227 :
228 14 : fZ1 => fZ1_ary
229 14 : fZ(1:4,1:4) => fZ1(1:4*4)
230 :
231 14 : ierr = 0
232 14 : work => work_ary
233 :
234 : ! ensure data have been initialized
235 14 : if (.not.table_atm_is_initialized) then
236 0 : write(*,*) 'get_table_values: table_atm_init required to proceed'
237 0 : ierr=-1
238 0 : return
239 : end if
240 :
241 14 : newlogZ(1) = safe_log10(newZ/GN93_Zsol)
242 :
243 16 : select case (id)
244 : case (ATM_TABLE_PHOTOSPHERE)
245 2 : ai => ai_two_thirds
246 : case (ATM_TABLE_TAU_100)
247 2 : ai => ai_100
248 : case (ATM_TABLE_TAU_10)
249 2 : ai => ai_10
250 : case (ATM_TABLE_TAU_1)
251 2 : ai => ai_1
252 : case (ATM_TABLE_TAU_1M1)
253 2 : ai => ai_1m1
254 : case (ATM_TABLE_WD_TAU_25)
255 2 : ai => ai_wd_25
256 : case (ATM_TABLE_DB_WD_TAU_25)
257 2 : ai => ai_db_wd_25
258 : case default
259 0 : write(*,*) 'Invalid id in get_table_values:', id
260 0 : ierr = -1
261 14 : return
262 : end select
263 :
264 14 : nZ = ai% nZ
265 14 : nT = ai% nT
266 14 : ng = ai% ng
267 :
268 14 : clip_Teff = .false.
269 14 : if (newTeff_in < ai% Teff_array(1)) then
270 0 : newTeff = ai% Teff_array(1) ! clip to table in T
271 0 : clip_Teff = .true.
272 14 : else if (newTeff_in > ai% Teff_array(nT)) then
273 0 : newTeff = ai% Teff_array(nT) ! clip to table in T
274 0 : clip_Teff = .true.
275 : else
276 14 : newTeff = newTeff_in
277 : end if
278 :
279 14 : clip_logg = .false.
280 14 : if (newlogg_in < ai% logg_array(1)) then
281 0 : newlogg = ai% logg_array(1) ! clip to table in logg
282 0 : clip_logg = .true.
283 14 : else if (newlogg_in > ai% logg_array(ng)) then
284 0 : newlogg = ai% logg_array(ng) ! clip to table in logg
285 0 : clip_logg = .true.
286 : else
287 14 : newlogg = newlogg_in
288 : end if
289 :
290 14 : allocate(result_2D(6,nZ))
291 :
292 14 : if (gtv_dbg) write(*,*) 'loaded tables: ', ai% have_atm_table(:)
293 :
294 : if (.false.) then !check for Z within range of tables
295 : if (nz > 1 .and. (newlogZ(1) < ai% logZ(1) .or. newlogZ(1) > ai% logZ(nZ))) then
296 : write(*,*) 'get_table_values: Z outside range of atm tables'
297 : ierr=-1
298 : return
299 : end if
300 14 : else if (newlogZ(1) < ai% logZ(1)) then
301 0 : newlogZ(1) = ai% logZ(1)
302 14 : else if (newlogZ(1) > ai% logZ(nZ)) then
303 12 : newlogZ(1) = ai% logZ(nZ)
304 : end if
305 :
306 : !identify surrounding Z values and initialize interpolation
307 14 : if (newlogZ(1) <= ai% logZ(1)) then
308 12 : Zindex = 1
309 2 : else if (newlogZ(1) >= ai% logZ(nz)) then
310 0 : Zindex = nz
311 : else
312 2 : Zindex = nz
313 16 : do j=1,nz-1
314 16 : if (ai% logZ(j) <= newlogZ(1)) then
315 14 : Zindex = j
316 : end if
317 : end do
318 : end if
319 :
320 : !identify upper (Zhi) and lower (Zlo) bounds on Z
321 14 : Zlo = max( 1,Zindex-1)
322 14 : Zhi = min(nZ,Zindex+2)
323 14 : if (Zhi-Zlo /= 3) Zlo=max( 1,Zlo-1)
324 14 : if (Zhi-Zlo /= 3) Zhi=min(nZ,Zhi+1)
325 14 : if (gtv_dbg) write(*,*) 'Zlo, Zindex, Zhi = ', Zlo, Zindex, Zhi
326 14 : numZs = Zhi - Zlo + 1
327 :
328 : !do 2D Teff,logg interpolation on each of 4 Z tables (Zlo-Zhi)
329 : ! P(g,T) dP_dg, dP_dT, d2P_dg2, d2P_dT2, d2P_dg_dT
330 : ! 1 2 3 4 5 6
331 14 : ict(:) = [ 1, 1, 1, 0, 0, 0 ]
332 :
333 14 : if (gtv_dbg) write(*,*) 'do_interp for Pgas', id
334 14 : call do_interp(ai% Pgas_interp1, newPgas, dPgas_dlogg, dPgas_dTeff, ierr)
335 14 : if (ierr /= 0) then
336 0 : if (gtv_dbg) write(*,*) 'do_interp failed'
337 0 : return
338 : end if
339 14 : if (clip_logg) dPgas_dlogg = 0
340 14 : if (clip_Teff) dPgas_dTeff = 0
341 :
342 14 : if (id == ATM_TABLE_PHOTOSPHERE) then
343 2 : newT = newTeff
344 2 : if (clip_Teff) then
345 0 : dT_dTeff = 0
346 : else
347 2 : dT_dTeff = 1
348 : end if
349 2 : dT_dlogg = 0
350 2 : return
351 : end if
352 :
353 12 : if (gtv_dbg) write(*,*) 'do_interp for Temp', id
354 12 : call do_interp(ai% T_interp1, newT, dT_dlogg, dT_dTeff, ierr)
355 12 : if (ierr /= 0) then
356 0 : if (gtv_dbg) write(*,*) 'do_interp failed'
357 0 : return
358 : end if
359 12 : if (clip_logg) dT_dlogg = 0
360 12 : if (clip_Teff) dT_dTeff = 0
361 :
362 12 : if (dbg .or. is_bad(newPgas) .or. is_bad(newT)) then
363 0 : write(*,1) 'newPgas', newPgas
364 0 : write(*,1) 'dPgas_dTeff', dPgas_dTeff
365 0 : write(*,1) 'dPgas_dlogg', dPgas_dlogg
366 0 : write(*,'(A)')
367 0 : write(*,1) 'newT', newT
368 0 : write(*,1) 'dT_dTeff', dT_dTeff
369 0 : write(*,1) 'dT_dlogg', dT_dlogg
370 0 : write(*,'(A)')
371 0 : ierr = -1
372 0 : return
373 : !if (is_bad(newPgas) .or. is_bad(newT)) call mesa_error(__FILE__,__LINE__,'get_table_values')
374 : end if
375 :
376 : !if (dbg) write(*,*) 'loaded tables: ', ai% have_atm_table(:)
377 :
378 40 : deallocate(result_2D)
379 :
380 : contains
381 :
382 26 : subroutine do_interp(f1, newval, dval_dlogg, dval_dTeff, ierr)
383 : real(dp), dimension(:), pointer :: f1
384 : real(dp), intent(out) :: newval, dval_dlogg, dval_dTeff
385 : integer, intent(out) :: ierr
386 :
387 182 : real(dp) :: res(6)
388 : integer :: j
389 26 : real(dp), pointer :: f(:)
390 :
391 : include 'formats'
392 :
393 26 : ierr = 0
394 :
395 58 : do i = Zlo, Zhi
396 32 : if (.not. ai% have_atm_table(i)) then
397 20 : call load_atm_table(i,ierr) !<-load on demand
398 : end if
399 32 : if (ierr /= 0) then
400 0 : if (gtv_dbg) write(*,*) 'load_atm_table failed'
401 0 : return
402 : end if
403 :
404 32 : f(1:4*ng*nT) => f1(1+4*ng*nT*(i-1):4*ng*nT*i)
405 : call interp_evbicub_db(newlogg, newTeff, ai% logg_array, ng, ai% Teff_array, nT, &
406 32 : ai% iling, ai% ilinT, f, ng, ict, res, ierr)
407 224 : do j=1,6
408 224 : result_2D(j,i) = res(j)
409 : end do
410 90 : if (ierr /= 0) then
411 0 : if (gtv_dbg) write(*,*) 'interp_evbicub_db failed'
412 0 : if (newlogg < ai% logg_array(1)) then
413 0 : write(*,1) 'logg too small for atm table', newlogg, ai% logg_array(1)
414 : end if
415 0 : if (newlogg > ai% logg_array(ng)) then
416 0 : write(*,1) 'logg too large for atm table', newlogg, ai% logg_array(ng)
417 : end if
418 0 : if (newTeff < ai% Teff_array(1)) then
419 0 : write(*,1) 'Teff too small for atm table', newTeff, ai% Teff_array(1)
420 : end if
421 0 : if (newTeff > ai% Teff_array(ng)) then
422 0 : write(*,1) 'Teff too large for atm table', newTeff, ai% Teff_array(ng)
423 : end if
424 0 : return
425 : end if
426 : end do
427 :
428 : ! now we have val, dval_dTeff, and dval_dlogg in result_2D for each Z
429 :
430 26 : if (numZs == 1) then
431 :
432 24 : newval = result_2D(1,Zlo)
433 24 : dval_dlogg = result_2D(2,Zlo)
434 24 : dval_dTeff = result_2D(3,Zlo)
435 :
436 : else ! Z interpolation
437 :
438 18 : fZ(1,1:numZs) = result_2D(1,Zlo:Zhi)
439 2 : call interp_pm(ai% logZ(Zlo:Zhi),numZs,fZ1,mp_work_size,work,'atm',ierr)
440 2 : if (ierr /= 0) then
441 0 : if (gtv_dbg) write(*,*) 'interp_pm failed for Z interpolation'
442 0 : return
443 : end if
444 2 : call interp_values(ai% logZ(Zlo:Zhi),numZs,fZ1,1,newlogZ,result_Z(1:1),ierr)
445 2 : if (ierr /= 0) then
446 0 : if (gtv_dbg) write(*,*) 'interp_values failed for Z interpolation'
447 0 : return
448 : end if
449 2 : newval = result_Z(1)
450 :
451 18 : fZ(1,1:numZs) = result_2D(2,Zlo:Zhi)
452 2 : call interp_pm(ai% logZ(Zlo:Zhi),numZs,fZ1,mp_work_size,work,'atm',ierr)
453 2 : if (ierr /= 0) then
454 0 : if (gtv_dbg) write(*,*) 'interp_pm failed for Z interpolation of d_dlogg'
455 0 : return
456 : end if
457 2 : call interp_values(ai% logZ(Zlo:Zhi),numZs,fZ1,1,newlogZ,result_Z(1:1),ierr)
458 2 : if (ierr /= 0) then
459 0 : if (gtv_dbg) write(*,*) 'interp_values failed for Z interpolation of d_dlogg'
460 0 : return
461 : end if
462 2 : dval_dlogg = result_Z(1)
463 :
464 18 : fZ(1,1:numZs) = result_2D(3,Zlo:Zhi)
465 2 : call interp_pm(ai% logZ(Zlo:Zhi),numZs,fZ1,mp_work_size,work,'atm',ierr)
466 2 : if (ierr /= 0) then
467 0 : if (gtv_dbg) write(*,*) 'interp_pm failed for Z interpolation of d_dTeff'
468 0 : return
469 : end if
470 2 : call interp_values(ai% logZ(Zlo:Zhi),numZs,fZ1,1,newlogZ,result_Z(1:1),ierr)
471 2 : if (ierr /= 0) then
472 0 : if (gtv_dbg) write(*,*) 'interp_values failed for Z interpolation of d_dTeff'
473 0 : return
474 : end if
475 2 : dval_dTeff = result_Z(1)
476 :
477 : end if
478 :
479 40 : end subroutine do_interp
480 :
481 :
482 :
483 : !reads in the iZth atm pressure table and pre-processes it for 2D spline interpolation
484 : !optional: reads or writes binary version of the spline table to atm_data/cache
485 20 : subroutine load_atm_table(iZ,ierr)
486 : use utils_lib
487 : use interp_2D_lib_db, only: interp_mkbicub_db
488 : use const_def, only: mesa_data_dir
489 : integer, intent(in) :: iZ !index of Z table to be loaded
490 : integer, intent(out) :: ierr
491 40 : integer :: iounit, i, j, ibound_tmp(ng), ibcTmin, ibcTmax, ibcgmin, ibcgmax
492 : integer :: text_file_version, nvec
493 2996 : real(dp) :: Teff_tmp(nT), logg_tmp(ng)
494 63844 : real(dp) :: bcTmin(nT), bcTmax(ng), bcgmin(ng), bcgmax(ng), data_tmp(ng,nT)
495 2020 : real(dp), target :: vec_ary(100)
496 20 : real(dp), pointer :: f1(:), vec(:)
497 : character(len=2000) :: buf
498 :
499 20 : ierr = 0
500 20 : vec => vec_ary
501 :
502 : !read in and process the text files
503 20 : filename = trim(mesa_data_dir) // '/atm_data/' // trim(ai% table_atm_files(iZ))
504 20 : open(newunit=iounit,file=trim(filename),action='read',status='old',iostat=ierr)
505 20 : if (ierr/= 0) then
506 0 : write(*,*) 'load_atm_table: missing atm data:'
507 0 : write(*,*) trim(filename)
508 0 : write(*,'(A)')
509 0 : write(*,'(A)')
510 0 : write(*,'(A)')
511 0 : write(*,'(A)')
512 0 : write(*,'(A)')
513 0 : write(*,*) 'FATAL ERROR: missing or bad atm data.'
514 0 : write(*,*) 'Please update by removing the directory mesa/data/atm_data,'
515 0 : write(*,*) 'and rerunning the mesa ./install script.'
516 0 : write(*,'(A)')
517 0 : call mesa_error(__FILE__,__LINE__)
518 : end if
519 :
520 20 : read(iounit,'(14x,i4)',iostat=ierr) text_file_version
521 20 : if (failed(1)) return
522 20 : if (text_file_version /= table_atm_version) then
523 0 : write(*,*) 'load_atm_table: mismatch in table versions'
524 0 : write(*,*) 'expected version ', table_atm_version
525 0 : write(*,*) 'received version ', text_file_version
526 0 : write(*,'(A)')
527 0 : write(*,'(A)')
528 0 : write(*,'(A)')
529 0 : write(*,'(A)')
530 0 : write(*,'(A)')
531 0 : write(*,'(A)')
532 0 : write(*,*) 'FATAL ERROR: out-of-date version of atm data.'
533 0 : write(*,*) 'Please update by removing the directory mesa/data/atm_data,'
534 0 : write(*,*) 'and rerunning the mesa ./install script.'
535 0 : write(*,'(A)')
536 0 : write(*,'(A)')
537 0 : call mesa_error(__FILE__,__LINE__)
538 : end if
539 :
540 346 : ibound_tmp = -1
541 20 : read(iounit,1,iostat=ierr) ai% logZ(iZ), ai% alphaFe(iZ), ai% atm_mix(iZ), ibound_tmp(1:ng)
542 20 : if (ierr /= 0) then
543 0 : write(*,*) 'iZ', iZ
544 0 : write(*,*) 'ai% logZ(iZ)', ai% logZ(iZ)
545 0 : write(*,*) 'ai% alphaFe(iZ)', ai% alphaFe(iZ)
546 0 : write(*,*) 'ai% atm_mix(iZ)', ai% atm_mix(iZ)
547 0 : do j=1,ng
548 0 : write(*,*) j, ibound_tmp(j)
549 : end do
550 0 : write(*,'(A)')
551 : end if
552 20 : if (failed(2)) return
553 20 : read(iounit,2,iostat=ierr) logg_tmp(:)
554 20 : if (failed(3)) return
555 60216 : data_tmp(:,:) = -1
556 2670 : do j=1,nT
557 2650 : read(iounit,'(a)',iostat=ierr) buf
558 : !write(*,*) 'read ierr', ierr
559 2650 : if (ierr == 0) then
560 2650 : call str_to_vector(buf, vec, nvec, ierr)
561 : !write(*,*) 'str_to_vector ierr', ierr
562 2650 : if (nvec < ng + 1) ierr = -1
563 : !write(*,*) 'nvec ierr', ierr
564 : end if
565 2650 : if (ierr /= 0) then
566 0 : write(*,*) 'failed reading Pgas Teff', ai% Teff_array(j)
567 0 : write(*,'(a)') 'buf "' // trim(buf) // '"'
568 0 : write(*,*) 'len_trim(buf)', len_trim(buf)
569 0 : write(*,*) 'nvec', nvec
570 0 : write(*,*) 'ng', ng
571 0 : return
572 : !stop
573 : end if
574 2650 : if (failed(4)) return
575 2650 : Teff_tmp(j) = vec(1)
576 60216 : do i=1,ng
577 60196 : data_tmp(i,j) = vec(i+1)
578 : end do
579 : end do
580 60216 : ai% Pgas_interp(1,:,:,iZ) = data_tmp(:,:)
581 :
582 20 : if (ai% id /= ATM_TABLE_PHOTOSPHERE) then ! read T
583 12 : read(iounit,2,iostat=ierr) ! skip line
584 12 : if (failed(5)) return
585 2006 : do j=1,nT
586 1994 : read(iounit,'(a)',iostat=ierr) buf
587 1994 : if (ierr == 0) then
588 1994 : call str_to_vector(buf, vec, nvec, ierr)
589 1994 : if (nvec < ng + 1) ierr = -1
590 : end if
591 1994 : if (ierr /= 0) then
592 0 : write(*,*) 'failed reading T Teff', ai% Teff_array(j)
593 : end if
594 1994 : if (failed(6)) return
595 1994 : Teff_tmp(j) = vec(1)
596 51024 : do i=1,ng
597 51012 : data_tmp(i,j) = vec(i+1)
598 : end do
599 : end do
600 51024 : ai% T_interp(1,:,:,iZ) = data_tmp(:,:)
601 : end if
602 :
603 20 : close(iounit)
604 :
605 : 1 format(13x,f5.2,8x,f4.1,1x,a8,1x,15x,99i4)
606 : 2 format(15x,99(9x,f5.2,1x))
607 : 3 format(99e15.7)
608 :
609 : !verify that Teff and logg arrays are the same as globals
610 2670 : do i=1,nT
611 2670 : if (abs(Teff_tmp(i) - ai% Teff_array(i)) > 1d-10*Teff_tmp(i)) then
612 0 : ierr = -1
613 0 : if (gtv_dbg) then
614 0 : write(*,'(a30,i6,e24.12)') 'Teff_tmp(i)', i, Teff_tmp(i)
615 0 : write(*,'(a30,i6,e24.12)') 'ai% Teff_array(i)', i, ai% Teff_array(i)
616 : end if
617 0 : return
618 : end if
619 : end do
620 :
621 346 : do i=1,ng
622 : !write(*,*) 'logg', i, logg_tmp(i), ai% logg_array(i)
623 346 : if (abs(logg_tmp(i) - ai% logg_array(i)) > 1d-10*abs(logg_tmp(i)) ) then
624 0 : ierr=-1
625 0 : if (gtv_dbg) then
626 0 : write(*,'(a30,i6,e24.12)') 'logg_tmp(i)', i, logg_tmp(i)
627 0 : write(*,'(a30,i6,e24.12)') 'ai% logg_array(i)', i, ai% logg_array(i)
628 : end if
629 0 : return
630 : end if
631 : end do
632 : !stop
633 :
634 : !make sure boundaries set earlier are still valid
635 346 : do i=1,ng
636 346 : ai% Teff_bound(i) = min( ai% Teff_bound(i), Teff_tmp(ibound_tmp(i)) )
637 : end do
638 :
639 : ! use "not a knot" bc's
640 2670 : ibcTmin = 0; bcTmin(:) = 0d0
641 346 : ibcTmax = 0; bcTmax(:) = 0d0
642 346 : ibcgmin = 0; bcgmin(:) = 0d0
643 346 : ibcgmax = 0; bcgmax(:) = 0d0
644 :
645 : !generate Teff,logg 2D-interpolants
646 20 : f1(1:4*ng*nT) => ai% Pgas_interp1(1+4*ng*nT*(iZ-1):4*ng*nT*iZ)
647 : call interp_mkbicub_db(ai% logg_array, ng, ai% Teff_array, nT, &
648 : f1, ng, ibcgmin, bcgmin, ibcgmax, bcgmax, &
649 20 : ibcTmin, bcTmin, ibcTmax, bcTmax, ai% iling, ai% ilinT, ierr )
650 20 : if (ierr /= 0) then
651 0 : if (gtv_dbg) write(*,*) 'interp_mkbicub_db failed for Pgas_interp'
652 0 : return
653 : end if
654 :
655 20 : if (ai% id /= ATM_TABLE_PHOTOSPHERE) then
656 12 : f1(1:4*ng*nT) => ai% T_interp1(1+4*ng*nT*(iZ-1):4*ng*nT*iZ)
657 : call interp_mkbicub_db(ai% logg_array, ng, ai% Teff_array, nT, &
658 : f1, ng, ibcgmin, bcgmin, ibcgmax, bcgmax, &
659 12 : ibcTmin, bcTmin, ibcTmax, bcTmax, ai% iling, ai% ilinT, ierr )
660 12 : if (ierr /= 0) then
661 0 : if (gtv_dbg) write(*,*) 'interp_mkbicub_db failed for T_interp'
662 0 : return
663 : end if
664 : end if
665 :
666 : !this file has been loaded and processed
667 20 : ai% have_atm_table(iZ) = .true.
668 :
669 20 : end subroutine load_atm_table
670 :
671 :
672 4716 : logical function failed(i)
673 : integer, intent(in) :: i
674 4716 : failed = (ierr /= 0)
675 4716 : if (failed) then
676 : write(*,'(a,i6)') &
677 0 : 'failed in reading atm table ' // trim(filename), i
678 : !call mesa_error(__FILE__,__LINE__,'get_table_values')
679 : end if
680 4716 : end function failed
681 :
682 : end subroutine get_table_values
683 :
684 : end module table_atm
|
