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 eosDT_load_tables
21 : use eos_def
22 : use utils_lib, only: is_bad, mesa_error, mv, switch_str
23 : use const_def, only: dp, use_mesa_temp_cache
24 : use math_lib
25 :
26 : implicit none
27 :
28 : ! the file EOS data
29 : integer, parameter :: jlogPgas = 1
30 : integer, parameter :: jlogE = 2
31 : integer, parameter :: jlogS = 3
32 : integer, parameter :: jchiRho = 4
33 : integer, parameter :: jchiT = 5
34 : integer, parameter :: jCp = 6
35 : integer, parameter :: jCv = 7
36 : integer, parameter :: jdE_dRho = 8
37 : integer, parameter :: jdS_dT = 9
38 : integer, parameter :: jdS_dRho = 10
39 : integer, parameter :: jmu = 11
40 : integer, parameter :: jlogfree_e = 12
41 : integer, parameter :: jgamma1 = 13
42 : integer, parameter :: jgamma3 = 14
43 : integer, parameter :: jgrad_ad = 15
44 : integer, parameter :: jeta = 16
45 : integer, parameter :: num_eos_file_vals = 16
46 :
47 : integer, parameter :: file_max_num_logQs = 1000
48 :
49 : contains
50 :
51 0 : subroutine request_user_to_reinstall
52 0 : write(*,'(A)')
53 0 : write(*,'(A)')
54 0 : write(*,'(A)')
55 0 : write(*,'(A)')
56 0 : write(*,'(A)')
57 0 : write(*,'(A)')
58 0 : write(*,*) 'NOTICE: you need to install a new version of the eos data.'
59 0 : write(*,*) 'Please update by removing the directory mesa/data/eosDT_data,'
60 0 : write(*,*) 'and rerunning the mesa ./install script.'
61 0 : write(*,'(A)')
62 0 : write(*,'(A)')
63 0 : call mesa_error(__FILE__,__LINE__)
64 0 : end subroutine request_user_to_reinstall
65 :
66 :
67 118 : subroutine check_for_error_in_eosDT_data(ierr, fname)
68 : integer, intent(in) :: ierr
69 : character (len=*) :: fname
70 118 : if (ierr == 0) return
71 0 : write(*,*) 'load eos tables ' // trim(fname)
72 0 : write(*,'(A)')
73 0 : write(*,'(A)')
74 0 : write(*,'(A)')
75 0 : write(*,'(A)')
76 0 : write(*,'(A)')
77 0 : write(*,'(a)') 'FATAL ERROR: missing or bad eos data.'
78 : write(*,'(a)') 'Please update by removing the directories ' &
79 : // 'mesa/data/eos*_data' &
80 0 : // ' and rerunning the mesa ./install script.'
81 0 : write(*,'(A)')
82 0 : call mesa_error(__FILE__,__LINE__)
83 : end subroutine check_for_error_in_eosDT_data
84 :
85 :
86 1685204 : subroutine load_single_eosDT_table_by_id( &
87 : rq, which_eosdt, ep, ix, iz, ierr)
88 : use utils_lib
89 : type (EoS_General_Info), pointer :: rq
90 : integer, intent(in) :: which_eosdt
91 : type (EosDT_XZ_Info), pointer :: ep
92 : integer,intent(in) :: iz, ix
93 : integer, intent(out) :: ierr
94 :
95 1685204 : if (which_eosdt == eosdt_max_FreeEOS) then
96 1684776 : ep => FreeEOS_XZ_data(ix,iz)
97 1684776 : if (FreeEOS_XZ_loaded(ix,iz)) return
98 428 : else if (which_eosdt == eosdt_OPAL_SCVH) then
99 428 : ep => eosDT_XZ_data(ix,iz)
100 428 : if (eosDT_XZ_loaded(ix,iz)) return
101 : else
102 0 : ierr = -1
103 0 : return
104 : end if
105 :
106 124 : !$OMP CRITICAL(eosDT_load)
107 124 : if (which_eosdt == eosdt_max_FreeEOS) then
108 94 : if (.not. FreeEOS_XZ_loaded(ix,iz)) call do_read
109 : else
110 30 : if (.not. eosDT_XZ_loaded(ix,iz)) call do_read
111 : end if
112 : !$OMP END CRITICAL(eosDT_load)
113 :
114 : contains
115 :
116 118 : subroutine do_read
117 118 : call read_one(ix,iz,ierr)
118 118 : if (ierr /= 0) return
119 118 : if (which_eosdt == eosdt_max_FreeEOS) then
120 88 : FreeEOS_XZ_loaded(ix,iz) = .true.
121 : else
122 30 : eosDT_XZ_loaded(ix,iz) = .true.
123 : end if
124 : end subroutine do_read
125 :
126 236 : subroutine read_one(ix,iz,ierr)
127 : use const_def, only: mesa_data_dir
128 : integer, intent(in) :: ix, iz
129 : integer, intent(out) :: ierr
130 : character (len=256) :: fname, cache_filename, temp_cache_filename
131 : integer :: iounit1, iounit2
132 118 : real(dp) :: X, Z
133 : type (DT_xz_Info), pointer :: xz
134 : include 'formats'
135 118 : iounit1 = alloc_iounit(ierr); if (ierr /= 0) return
136 118 : iounit2 = alloc_iounit(ierr); if (ierr /= 0) return
137 118 : if (which_eosdt == eosdt_max_FreeEOS) then
138 88 : xz => FreeEOS_xz_struct
139 : else
140 30 : xz => eosDT_xz_struct
141 : end if
142 : call Get_eosDT_Table_Filenames(rq, which_eosdt, xz, &
143 118 : ix, iz, mesa_data_dir, fname, cache_filename, temp_cache_filename)
144 : call Load1_eosDT_Table(rq, which_eosdt, ep, xz, &
145 118 : ix, iz, fname, cache_filename, temp_cache_filename, iounit1, iounit2, use_cache_for_eos, ierr)
146 118 : if (ierr /= 0) then
147 0 : write(*,*) 'Load1_eosDT_Table ierr', ierr, ix, iz, X, Z
148 : end if
149 118 : call free_iounit(iounit2)
150 118 : call free_iounit(iounit1)
151 : end subroutine read_one
152 :
153 : end subroutine load_single_eosDT_table_by_id
154 :
155 :
156 118 : subroutine Get_eosDT_Table_Filenames(rq, which_eosdt, xz, &
157 : ix, iz, data_dir, fname, cache_filename, temp_cache_filename)
158 : type (EoS_General_Info), pointer :: rq
159 : integer, intent(in) :: which_eosdt
160 : type (DT_xz_Info), pointer :: xz
161 : integer, intent(in) :: ix, iz
162 : character (*), intent(in) :: data_dir
163 : character (len=*), intent(out) :: fname, cache_filename, temp_cache_filename
164 : character (len=256) :: Zstr, Xstr, suffix, data_dir_name, data_prefix
165 : real(dp) :: X, Z
166 :
167 118 : Z = xz% Zs(iz)
168 118 : X = xz% Xs_for_Z(ix,iz)
169 :
170 118 : call setstr(Z,Zstr)
171 118 : call setstr(X,Xstr)
172 118 : suffix = ''
173 118 : if (which_eosdt == eosdt_max_FreeEOS) then
174 88 : data_dir_name = '/eosFreeEOS_data/'
175 88 : data_prefix = '-FreeEOS_'
176 88 : suffix = rq% suffix_for_FreeEOS_Z(iz)
177 : else
178 30 : data_dir_name = '/eosDT_data/'
179 30 : data_prefix = '-eosDT_'
180 : end if
181 :
182 : fname = trim(data_dir) // &
183 : trim(data_dir_name) // trim(rq% eosDT_file_prefix) // trim(data_prefix) // &
184 118 : trim(Zstr) // 'z' // trim(Xstr) // 'x' // trim(suffix) // '.data'
185 : cache_filename = trim(eosDT_cache_dir) // &
186 : '/' // trim(rq% eosDT_file_prefix) // trim(data_prefix) // &
187 118 : trim(Zstr) // 'z' // trim(Xstr) // 'x' // trim(suffix) // '.bin'
188 : temp_cache_filename = trim(eosDT_temp_cache_dir) // &
189 : '/' // trim(rq% eosDT_file_prefix) // trim(data_prefix) // &
190 236 : trim(Zstr) // 'z' // trim(Xstr) // 'x' // trim(suffix) // '.bin'
191 :
192 : contains
193 :
194 236 : subroutine setstr(v,str)
195 : real(dp), intent(in) :: v
196 : character (len=*) :: str
197 236 : if (v > 0.99999d0) then
198 6 : str = '100'
199 230 : else if (v > 0.09999d0) then
200 104 : write(str, '(i2)') floor(100d0 * v + 0.5d0)
201 : else
202 126 : write(str, '(a,i1)') '0', floor(100d0 * v + 0.5d0)
203 : end if
204 236 : end subroutine setstr
205 :
206 : end subroutine Get_eosDT_Table_Filenames
207 :
208 :
209 118 : subroutine Load1_eosDT_Table(rq, which_eosdt, ep, xz, &
210 : ix, iz, filename, cache_filename, temp_cache_filename, &
211 : io_unit, cache_io_unit, use_cache, info)
212 : type (EoS_General_Info), pointer :: rq
213 : integer, intent(in) :: which_eosdt
214 : type (EosDT_XZ_Info), pointer :: ep
215 : type (DT_xz_Info), pointer :: xz
216 : integer, intent(in) :: ix, iz
217 : character (*), intent(in) :: filename, cache_filename, temp_cache_filename
218 : integer, intent(in) :: io_unit, cache_io_unit
219 : logical, intent(in) :: use_cache
220 : integer, intent(out) :: info
221 :
222 236 : real(dp) :: X, Z, logQ, logT, X_in, Z_in
223 : integer :: j, i, k, iQ, ios, status
224 : character (len=1000) :: message
225 : real(dp), parameter :: tiny = 1d-10
226 118 : real(dp), pointer :: tbl(:,:,:,:) ! => ep% tbl1
227 118 : real(dp), pointer :: tbl2_1(:), tbl2(:,:,:)
228 6018 : real(dp), target :: vec_ary(50)
229 : real(dp), pointer :: vec(:)
230 : integer :: n
231 :
232 : include 'formats'
233 :
234 118 : info = 0
235 118 : vec => vec_ary
236 118 : Z = xz% Zs(iz)
237 118 : X = xz% Xs_for_Z(ix,iz)
238 :
239 118 : write(message,*) 'open ', trim(filename)
240 :
241 118 : open(UNIT=io_unit, FILE=trim(filename), ACTION='READ', STATUS='OLD', IOSTAT=ios)
242 118 : call check_for_error_in_eosDT_data(ios, filename)
243 :
244 118 : read(io_unit,*,iostat=info)
245 236 : if (info /= 0) return
246 :
247 118 : read(io_unit,'(a)',iostat=info) message
248 118 : if (info == 0) call str_to_vector(message, vec, n, info)
249 118 : if (info /= 0 .or. n < 11) then
250 0 : write(*,'(a)') 'failed while reading ' // trim(filename)
251 0 : close(io_unit)
252 0 : info = -1
253 0 : return
254 : end if
255 118 : ep% version = int(vec(1))
256 118 : X_in = vec(2)
257 118 : Z_in = vec(3)
258 118 : ep% num_logTs = int(vec(4))
259 118 : ep% logT_min = vec(5)
260 118 : ep% logT_max = vec(6)
261 118 : ep% del_logT = vec(7)
262 118 : ep% num_logQs = int(vec(8))
263 118 : ep% logQ_min = vec(9)
264 118 : ep% logQ_max = vec(10)
265 118 : ep% del_logQ = vec(11)
266 :
267 118 : read(io_unit,*,iostat=info)
268 118 : if (info /= 0) return
269 :
270 118 : if (abs(X-X_in) > tiny .or. abs(Z-Z_in) > tiny) then
271 0 : write(*,*) 'bad header info in ' // trim(filename)
272 0 : info = -1
273 0 : close(io_unit)
274 0 : if (abs(X-X_in) > tiny) then
275 0 : write(*,'(a50,l1)') 'abs(X-X_in) > tiny', abs(X-X_in) > tiny
276 : end if
277 0 : if (abs(Z-Z_in) > tiny) then
278 0 : write(*,'(a50,l1)') 'abs(Z-Z_in) > tiny', abs(Z-Z_in) > tiny
279 : end if
280 0 : write(*,'(A)')
281 0 : call request_user_to_reinstall
282 0 : return
283 : end if
284 :
285 : if (show_allocations) write(*,2) 'Load1_eosDT_Table ep% tbl1', &
286 : sz_per_eos_point*nv*ep% num_logQs*ep% num_logTs + ep% num_logQs + ep% num_logTs
287 : allocate(ep% tbl1(sz_per_eos_point*nv*ep% num_logQs*ep% num_logTs), &
288 : ep% logQs(ep% num_logQs), ep% logTs(ep% num_logTs), &
289 118 : STAT=info)
290 118 : if (info /= 0) then
291 0 : write(*,*) "Info: ",info
292 0 : call mesa_error(__FILE__,__LINE__, "Allocation in Load1_eosDT_Table failed, you're likely out of memory")
293 : end if
294 :
295 : tbl(1:sz_per_eos_point,1:nv,1:ep% num_logQs,1:ep% num_logTs) => &
296 118 : ep% tbl1(1:sz_per_eos_point*nv*ep% num_logQs*ep% num_logTs)
297 :
298 118 : ep% logQs(1) = ep% logQ_min
299 61494 : do i = 2, ep% num_logQs-1
300 61494 : ep% logQs(i) = ep% logQs(i-1) + ep% del_logQ
301 : end do
302 118 : ep% logQs(ep% num_logQs) = ep% logQ_max
303 :
304 118 : ep% logTs(1) = ep% logT_min
305 35990 : do i = 2, ep% num_logTs-1
306 35990 : ep% logTs(i) = ep% logTs(i-1) + ep% del_logT
307 : end do
308 118 : ep% logTs(ep% num_logTs) = ep% logT_max
309 :
310 118 : if (use_cache) then
311 0 : call Read_EoS_Cache(X, Z, ep, cache_filename, cache_io_unit, ios)
312 0 : if (ios == 0) then
313 0 : close(io_unit)
314 0 : return
315 : end if
316 : end if
317 :
318 118 : status = 0
319 118 : allocate(tbl2_1(num_eos_file_vals*ep% num_logQs*ep% num_logTs), STAT=status)
320 : if (status /= 0) then
321 0 : info = -1
322 0 : return
323 : end if
324 :
325 : tbl2(1:num_eos_file_vals,1:ep% num_logQs,1:ep% num_logTs) => &
326 118 : tbl2_1(1:num_eos_file_vals*ep% num_logQs*ep% num_logTs)
327 :
328 61612 : do iQ=1,ep% num_logQs
329 :
330 61612 : read(io_unit,*,iostat=info)
331 61612 : if (failed('skip line')) return
332 :
333 61612 : read(io_unit,'(a)',iostat=info) message
334 61612 : if (info == 0) call str_to_double(message, vec(1), info)
335 61612 : if (failed('read logQ')) return
336 61612 : logQ = vec(1)
337 :
338 61612 : read(io_unit,*,iostat=info)
339 61612 : if (failed('skip line')) return
340 :
341 61612 : read(io_unit,*,iostat=info)
342 61612 : if (failed('skip line')) return
343 :
344 18914884 : do i=1,ep% num_logTs
345 :
346 18853272 : read(io_unit,'(a)',iostat=info) message
347 18853272 : if (failed('read line')) then
348 0 : write(*,'(a)') trim(message)
349 0 : write(*,*) trim(filename)
350 0 : write(*,*) 'iQ, i', iQ, i
351 0 : write(*,*) 'logQ', logQ
352 0 : write(*,*) 'bad input line?'
353 0 : call mesa_error(__FILE__,__LINE__)
354 : end if
355 :
356 18853272 : call str_to_vector(message, vec, n, info)
357 18853272 : if (info /= 0 .or. n < 1+num_eos_file_vals) then
358 0 : write(*,'(a)') trim(message)
359 0 : write(*,*) trim(filename)
360 0 : write(*,*) 'iQ, i', iQ, i
361 0 : write(*,*) 'logQ', logQ
362 0 : write(*,*) 'bad input line?'
363 0 : call mesa_error(__FILE__,__LINE__)
364 : end if
365 18853272 : logT = vec(1)
366 339420508 : do j=1,num_eos_file_vals
367 320505624 : tbl2(j,iQ,i) = vec(1+j)
368 : end do
369 :
370 : end do
371 :
372 61612 : if(iQ == ep% num_logQs) exit
373 61494 : read(io_unit,*,iostat=info)
374 61494 : if (failed('skip line')) return
375 61494 : read(io_unit,*,iostat=info)
376 61612 : if (failed('skip line')) return
377 :
378 : end do
379 :
380 118 : close(io_unit)
381 :
382 118 : call Make_XEoS_Interpolation_Data(ep, tbl2_1, info)
383 118 : deallocate(tbl2_1)
384 118 : if (failed('Make_XEoS_Interpolation_Data')) return
385 :
386 118 : call Check_XEoS_Interpolation_Data(ep)
387 :
388 118 : if (.not. use_cache) return
389 :
390 : open(unit=cache_io_unit, &
391 : file=trim(switch_str(temp_cache_filename, cache_filename, use_mesa_temp_cache)), &
392 0 : iostat=ios,action='write', form='unformatted')
393 :
394 236 : if (ios == 0) then
395 0 : write(*,'(a)') 'write ' // trim(cache_filename)
396 : write(cache_io_unit) &
397 0 : X_in, Z_in, ep% num_logTs, ep% logT_min, ep% logT_max, ep% del_logT, &
398 0 : ep% num_logQs, ep% logQ_min, ep% logQ_max, ep% del_logQ, ep% version
399 : write(cache_io_unit) &
400 : ep% tbl1(&
401 0 : 1:sz_per_eos_point*nv*ep% num_logQs*ep% num_logTs)
402 0 : close(cache_io_unit)
403 0 : if(use_mesa_temp_cache) call mv(temp_cache_filename, cache_filename,.true.)
404 : end if
405 :
406 :
407 : contains
408 :
409 118 : subroutine Check_XEoS_Interpolation_Data(ep)
410 : use utils_lib,only:is_bad
411 : type (EosDT_XZ_Info), pointer :: ep
412 : ! for logT > 6.8 and logRho < -10, splines can get bogus higher order terms
413 : ! replace NaN's and Infinities with 0
414 : integer :: i, j, iQ, jtemp
415 590 : do i = 1, sz_per_eos_point
416 12862 : do j = 1, nv
417 6420392 : do iQ = 1, ep% num_logQs
418 1967160208 : do jtemp = 1, ep% num_logTs
419 1967147936 : if (is_bad(tbl(i,j,iQ,jtemp))) then
420 0 : tbl(i,j,iQ,jtemp) = 0
421 : end if
422 : end do
423 : end do
424 : end do
425 : end do
426 118 : end subroutine Check_XEoS_Interpolation_Data
427 :
428 19222826 : logical function failed(str)
429 : character (len=*), intent(in) :: str
430 19222826 : failed = (info /= 0)
431 19222826 : if (failed) write(*,*) 'Load1_eosDT_Table failed: ' // trim(str)
432 19222826 : end function failed
433 :
434 :
435 : end subroutine Load1_eosDT_Table
436 :
437 :
438 118 : subroutine Make_XEoS_Interpolation_Data(ep, tbl2_1, info)
439 : use interp_2d_lib_db
440 : use const_def, only: crad, ln10
441 :
442 : type (EosDT_XZ_Info), pointer :: ep
443 : real(dp), pointer :: tbl2_1(:) ! =(num_eos_file_vals, ep% num_logQs, ep% num_logTs)
444 : integer, intent(out) :: info
445 :
446 61730 : real(dp) :: logQs(ep% num_logQs) ! x vector, strict ascending
447 36226 : real(dp) :: logTs(ep% num_logTs) ! y vector, strict ascending
448 : real(dp) :: Ts(ep% num_logTs)
449 118 : real(dp), allocatable, target :: f1_ary(:) ! data & spline coefficients
450 118 : real(dp), pointer :: f1(:), f(:,:,:), ep_tbl(:,:,:,:), tbl2(:,:,:)
451 : integer :: ibcxmin ! bc flag for x=xmin
452 36226 : real(dp) :: bcxmin(ep% num_logTs) ! bc data vs. y at x=xmin
453 : integer :: ibcxmax ! bc flag for x=xmax
454 36226 : real(dp) :: bcxmax(ep% num_logTs) ! bc data vs. y at x=xmax
455 : integer :: ibcymin ! bc flag for y=ymin
456 61730 : real(dp) :: bcymin(ep% num_logQs) ! bc data vs. x at y=ymin
457 : integer :: ibcymax ! bc flag for y=ymax
458 61848 : real(dp) :: bcymax(ep% num_logQs) ! bc data vs. x at y=ymax
459 : integer :: ili_logQs ! =1: logRho grid is "nearly" equally spaced
460 : integer :: ili_logTs ! =1: logT grid is "nearly" equally spaced
461 : integer :: ier ! =0 on exit if there is no error.
462 : real(dp) :: logQ, Rho, logRho, T, P, Cv, chiRho, chiT, logT, logT0, logT1, logQ0, logQ1
463 : real(dp) :: gamma3, gamma1, grad_ad, Prad, E, S
464 : integer :: iQ, jtemp, ilogT, ilogQ
465 : real(dp) :: fval(num_eos_file_vals), df_dx(num_eos_file_vals), df_dy(num_eos_file_vals)
466 :
467 : real(dp) :: x, y, dlnT, energy, lnE, entropy, lnS, Pgas, lnPgas, dlogT, &
468 : dlnPgas_dlnd, dlnE_dlnd, dlnS_dlnd, dlnPgas_dlnT, dlnE_dlnT, dlnS_dlnT
469 :
470 : integer :: v, vlist(3), var, i, j, num_logQs, num_logTs, ii, jj
471 : character (len=256) :: message
472 :
473 : include 'formats'
474 :
475 118 : info = 0
476 :
477 : ! just use "not a knot" bc's at edges of tables
478 36226 : ibcxmin = 0; bcxmin(:) = 0
479 36226 : ibcxmax = 0; bcxmax(:) = 0
480 61730 : ibcymin = 0; bcymin(:) = 0
481 61730 : ibcymax = 0; bcymax(:) = 0
482 :
483 118 : num_logQs = ep% num_logQs
484 118 : num_logTs = ep% num_logTs
485 :
486 : ep_tbl(1:sz_per_eos_point,1:nv,1:num_logQs,1:num_logTs) => &
487 118 : ep% tbl1(1:sz_per_eos_point*nv*num_logQs*num_logTs)
488 :
489 : tbl2(1:num_eos_file_vals,1:num_logQs,1:num_logTs) => &
490 118 : tbl2_1(1:num_eos_file_vals*num_logQs*num_logTs)
491 :
492 118 : allocate(f1_ary(sz_per_eos_point * ep% num_logQs * ep% num_logTs))
493 :
494 118 : f1 => f1_ary
495 : f(1:sz_per_eos_point,1:num_logQs,1:num_logTs) => &
496 118 : f1_ary(1:sz_per_eos_point*num_logQs*num_logTs)
497 :
498 61730 : do iQ = 1, ep% num_logQs
499 61730 : logQs(iQ) = ep% logQ_min + (iQ-1) * ep% del_logQ
500 : end do
501 :
502 36226 : do jtemp = 1, ep% num_logTs
503 36226 : logTs(jtemp) = ep% logT_min + (jtemp-1) * ep% del_logT
504 : end do
505 :
506 : ! copy file eos variables to internal eos interpolation tables
507 36226 : do j=1,num_logTs
508 18889498 : do i=1,num_logQs
509 18853272 : ep_tbl(1,i_lnPgas,i,j) = tbl2(jlogPgas,i,j)*ln10
510 18853272 : ep_tbl(1,i_lnE,i,j) = tbl2(jlogE,i,j)*ln10
511 18853272 : ep_tbl(1,i_lnS,i,j) = tbl2(jlogS,i,j)*ln10
512 18853272 : ep_tbl(1,i_grad_ad,i,j) = tbl2(jgrad_ad,i,j)
513 18853272 : ep_tbl(1,i_chiRho,i,j) = tbl2(jchiRho,i,j)
514 18853272 : ep_tbl(1,i_chiT,i,j) = tbl2(jchiT,i,j)
515 18853272 : ep_tbl(1,i_Cp,i,j) = tbl2(jCp,i,j)
516 18853272 : ep_tbl(1,i_Cv,i,j) = tbl2(jCv,i,j)
517 18853272 : ep_tbl(1,i_dE_dRho,i,j) = tbl2(jdE_dRho,i,j)
518 18853272 : ep_tbl(1,i_dS_dT,i,j) = tbl2(jdS_dT,i,j)
519 18853272 : ep_tbl(1,i_dS_dRho,i,j) = tbl2(jdS_dRho,i,j)
520 18853272 : ep_tbl(1,i_mu,i,j) = tbl2(jmu,i,j)
521 18853272 : ep_tbl(1,i_lnfree_e,i,j) = max(-4d0,tbl2(jlogfree_e,i,j))*ln10
522 : ! to protect against non-monotonic interpolation caused by extreme values
523 18853272 : ep_tbl(1,i_gamma1,i,j) = tbl2(jgamma1,i,j)
524 18853272 : ep_tbl(1,i_gamma3,i,j) = tbl2(jgamma3,i,j)
525 18889380 : ep_tbl(1,i_eta,i,j) = tbl2(jeta,i,j)
526 : end do
527 : end do
528 :
529 : ! create tables for bicubic spline interpolation
530 3186 : do v = 1, nv
531 1604980 : do i=1,ep% num_logQs
532 491790052 : do j=1,ep% num_logTs
533 491786984 : f(1,i,j) = ep_tbl(1,v,i,j)
534 : end do
535 : end do
536 : call interp_mkbicub_db( &
537 : logQs,ep% num_logQs,logTs,ep% num_logTs,f1,ep% num_logQs, &
538 : ibcxmin,bcxmin,ibcxmax,bcxmax, &
539 : ibcymin,bcymin,ibcymax,bcymax, &
540 3068 : ili_logQs,ili_logTs,ier)
541 3068 : if (ier /= 0) then
542 0 : write(*,*) 'interp_mkbicub_db error happened for eos_value', v
543 0 : info = 3
544 0 : return
545 : end if
546 1605098 : do i=1,ep% num_logQs
547 491790052 : do j=1,ep% num_logTs
548 490185072 : ep_tbl(2,v,i,j) = f(2,i,j)
549 490185072 : ep_tbl(3,v,i,j) = f(3,i,j)
550 491786984 : ep_tbl(4,v,i,j) = f(4,i,j)
551 : end do
552 : end do
553 : end do
554 :
555 :
556 236 : end subroutine Make_XEoS_Interpolation_Data
557 :
558 :
559 0 : subroutine Read_EoS_Cache(X, Z, ep, cache_filename, io_unit, ios)
560 : real(dp), intent(in) :: X, Z
561 : type (EosDT_XZ_Info), pointer :: ep
562 : character (*), intent(in) :: cache_filename
563 : integer, intent(in) :: io_unit ! use this for file access
564 : integer, intent(out) :: ios
565 :
566 0 : real(dp) :: X_in, Z_in, logT_min_in, logT_max_in, del_logT_in, &
567 0 : logQ_min_in, logQ_max_in, del_logQ_in
568 : integer :: num_logQs_in, num_logTs_in, version_in
569 : real(dp), parameter :: tiny = 1d-10
570 :
571 : include 'formats'
572 :
573 0 : ios = 0
574 : open(unit=io_unit,file=trim(cache_filename),action='read', &
575 0 : status='old',iostat=ios,form='unformatted')
576 0 : if (ios /= 0) return
577 :
578 : read(io_unit, iostat=ios) &
579 0 : X_in, Z_in, num_logTs_in, logT_min_in, logT_max_in, del_logT_in, &
580 0 : num_logQs_in, logQ_min_in, logQ_max_in, del_logQ_in, version_in
581 0 : if (ios /= 0) return
582 :
583 0 : if (ep% version /= version_in) then
584 0 : ios = 1
585 0 : write(*,*) 'read cache failed for version_in'
586 : end if
587 0 : if (ep% num_logQs /= num_logQs_in) then
588 0 : ios = 1
589 0 : write(*,*) 'read cache failed for ep% num_logQs'
590 : end if
591 0 : if (ep% num_logTs /= num_logTs_in) then
592 0 : ios = 1
593 0 : write(*,*) 'read cache failed for ep% num_logTs'
594 : end if
595 0 : if (abs(X-X_in) > tiny) then
596 0 : ios = 1
597 0 : write(*,*) 'read cache failed for X_in'
598 : end if
599 0 : if (abs(Z-Z_in) > tiny) then
600 0 : ios = 1
601 0 : write(*,*) 'read cache failed for Z_in'
602 : end if
603 0 : if (abs(ep% logT_min-logT_min_in) > tiny) then
604 0 : ios = 1
605 0 : write(*,*) 'read cache failed for eos_logT_min'
606 : end if
607 0 : if (abs(ep% logT_max-logT_max_in) > tiny) then
608 0 : ios = 1
609 0 : write(*,*) 'read cache failed for eos_logT_max'
610 : end if
611 0 : if (abs(ep% del_logT-del_logT_in) > tiny) then
612 0 : ios = 1
613 0 : write(*,*) 'read cache failed for eos_del_logT'
614 : end if
615 0 : if (abs(ep% logQ_min-logQ_min_in) > tiny) then
616 0 : ios = 1
617 0 : write(*,*) 'read cache failed for eos_logQ_min'
618 : end if
619 0 : if (abs(ep% logQ_max-logQ_max_in) > tiny) then
620 0 : ios = 1
621 0 : write(*,*) 'read cache failed for eos_logQ_max'
622 : end if
623 0 : if (abs(ep% del_logQ-del_logQ_in) > tiny) then
624 0 : ios = 1
625 0 : write(*,*) 'read cache failed for eos_del_logQ'
626 : end if
627 :
628 0 : if (ios /= 0) then
629 0 : close(io_unit); return
630 : end if
631 :
632 : read(io_unit, iostat=ios) &
633 : ep% tbl1( &
634 0 : 1:sz_per_eos_point*nv*ep% num_logQs*ep% num_logTs)
635 0 : if (ios /= 0) then
636 0 : close(io_unit); return
637 : end if
638 :
639 0 : close(io_unit)
640 :
641 : end subroutine Read_EoS_Cache
642 :
643 : end module eosDT_load_tables
|
