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 load_kap
21 :
22 : use kap_def
23 : use load_CO_kap, only: Setup_Kap_CO_Tables, min_version
24 : use kapcn, only: kapcn_init
25 : use kap_aesopus, only: kap_aesopus_init
26 : use const_def, only: dp, use_mesa_temp_cache
27 : use math_lib
28 : use utils_lib, only: mesa_error, mv, switch_str
29 :
30 : implicit none
31 :
32 : private
33 : public :: load_one
34 : public :: setup_kap_tables
35 :
36 : logical, parameter :: dbg = .false.
37 : logical, parameter :: dbg_cache = .false.
38 :
39 : contains
40 :
41 11 : subroutine setup_kap_tables(rq, use_cache, load_on_demand, ierr)
42 : use const_def, only: mesa_data_dir
43 : use condint, only: init_potekhin
44 : type (Kap_General_Info), pointer :: rq
45 : logical, intent(in) :: use_cache, load_on_demand
46 : integer, intent(out) :: ierr
47 :
48 : integer :: id, iz, ix, num_Xs
49 11 : real(dp) :: X, Z
50 :
51 : include 'formats'
52 :
53 11 : ierr = 0
54 :
55 11 : kap_dir = trim(mesa_data_dir) // '/kap_data'
56 :
57 11 : kap_use_cache = use_cache
58 :
59 11 : call Setup_Kap_CO_Tables(rq, kap_CO_z_tables(rq% kap_CO_option)% ar, use_cache, load_on_demand, ierr)
60 11 : if (ierr /= 0) return
61 :
62 11 : call setup_lowT(kap_lowT_z_tables(rq% kap_lowT_option)% ar)
63 11 : call setup(kap_z_tables(rq% kap_option)% ar)
64 :
65 : rq% logT_Compton_blend_hi = & ! apply whichever is minimum for Type1 and Type2 options
66 : min(kap_z_tables(rq% kap_option)% ar(1)% x_tables(1)% logT_max - 0.01d0, &
67 11 : kap_CO_z_tables(rq% kap_CO_option)% ar(1)% x_tables(1)% logT_max - 0.01d0)
68 : !rq% kap_z_tables(1)% x_tables(1)% logT_max - 0.01d0
69 11 : rq% logR_Compton_blend_lo = kap_z_tables(rq% kap_option)% ar(1)% x_tables(1)% logR_min + 0.01d0
70 : !rq% kap_z_tables(1)% x_tables(1)% logR_min + 0.01d0
71 :
72 11 : call init_potekhin(ierr)
73 :
74 : contains
75 :
76 11 : subroutine setup(z_tables)
77 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
78 : integer :: iz
79 : include 'formats'
80 11 : if (.not. associated(z_tables)) then
81 98 : allocate(z_tables(num_kap_Zs(rq% kap_option)), STAT=ierr)
82 7 : if (ierr /= 0) return
83 98 : do iz = 1, num_kap_Zs(rq% kap_option)
84 91 : z_tables(iz)% lowT_flag = .false.
85 91 : z_tables(iz)% Z = kap_Zs(iz, rq% kap_option)
86 1001 : allocate(z_tables(iz)% x_tables(num_kap_Xs(rq% kap_option)), STAT=ierr)
87 98 : if (ierr /= 0) return
88 : end do
89 : if (show_allocations) write(*,2) 'setup z_tables', &
90 : num_kap_Zs(rq% kap_option) + num_kap_Zs(rq% kap_option)*num_kap_Xs(rq% kap_option)
91 : end if
92 154 : do iz = 1, num_kap_Zs(rq% kap_option)
93 143 : Z = kap_Zs(iz,rq% kap_option)
94 143 : num_Xs = num_kap_Xs_for_this_Z(iz,rq% kap_option)
95 143 : z_tables(iz)% num_Xs = num_Xs
96 1540 : do ix = 1, num_Xs
97 1386 : X = kap_Xs(ix,rq% kap_option)
98 1386 : if (X + Z > 1) X = 1-Z
99 : call load_one(rq, z_tables, iz, ix, X, Z, &
100 1386 : load_on_demand, ierr)
101 1529 : if (ierr /= 0) return
102 : end do
103 : end do
104 :
105 11 : end subroutine setup
106 :
107 11 : subroutine setup_lowT(z_tables)
108 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
109 : integer :: iz
110 : include 'formats'
111 11 : select case (rq% kap_lowT_option)
112 : case (kap_lowT_kapCN)
113 0 : if (rq% show_info) write(*,*) 'Loading lowT kapCN tables'
114 0 : call kapCN_init(rq% show_info, ierr)
115 0 : if (ierr /= 0) then
116 0 : write(*,*) 'Failed to load kapCN opacities'
117 0 : call mesa_error(__FILE__, __LINE__)
118 : end if
119 : case (kap_lowT_AESOPUS)
120 2 : if (rq% show_info) write(*,*) 'Loading lowT AESOPUS tables'
121 2 : call kap_aesopus_init(rq, ierr)
122 2 : if (ierr /= 0) then
123 0 : write(*,*) 'Failed to load AESOPUS opacities'
124 0 : call mesa_error(__FILE__, __LINE__)
125 : end if
126 : case default
127 9 : if (rq% kap_lowT_option > kap_lowT_options_max .or. &
128 : rq% kap_lowT_option <= 0) then
129 0 : write(*,*) 'invalid lowT opacity option'
130 0 : call mesa_error(__FILE__, __LINE__)
131 : end if
132 9 : if (rq% show_info) then
133 0 : if (rq% kap_lowT_option == kap_lowT_test) then
134 0 : write(*,*) 'Loading lowT test tables'
135 0 : else if (rq% kap_lowT_option == kap_lowT_Freedman11) then
136 0 : write(*,*) 'Loading lowT Freedman tables'
137 : else
138 0 : write(*,*) 'Loading lowT Fergusson tables'
139 : end if
140 : end if
141 9 : if (.not. associated(z_tables)) then
142 86 : allocate(z_tables(num_kap_lowT_Zs(rq% kap_lowT_option)), STAT=ierr)
143 7 : if (ierr /= 0) return
144 86 : do iz = 1, num_kap_lowT_Zs(rq% kap_lowT_option)
145 79 : z_tables(iz)% lowT_flag = .true.
146 79 : z_tables(iz)% Z = kap_lowT_Zs(iz, rq% kap_lowT_option)
147 743 : allocate(z_tables(iz)% x_tables(num_kap_lowT_Xs(rq% kap_lowT_option)), STAT=ierr)
148 86 : if (ierr /= 0) return
149 : end do
150 : if (show_allocations) write(*,2) 'setup_lowT z_tables', &
151 : num_kap_lowT_Zs(rq% kap_lowT_option) + num_kap_lowT_Zs(rq% kap_lowT_option)*num_kap_lowT_Xs(rq% kap_lowT_option)
152 : end if
153 127 : do iz = 1, num_kap_lowT_Zs(rq% kap_lowT_option)
154 105 : Z = kap_lowT_Zs(iz, rq% kap_lowT_option)
155 105 : num_Xs = num_kap_lowT_Xs_for_this_Z(iz, rq% kap_lowT_option)
156 105 : z_tables(iz)% num_Xs = num_Xs
157 1010 : do ix = 1, num_Xs
158 896 : X = kap_lowT_Xs(ix, rq% kap_lowT_option)
159 896 : if (X + Z > 1) X = 1-Z
160 : call load_one(rq, z_tables, iz, ix, X, Z, &
161 896 : load_on_demand, ierr)
162 1001 : if (ierr /= 0) then
163 0 : if (rq% kap_lowT_option == kap_lowT_test) then
164 0 : write(*,*) 'Failed to load lowT test tables'
165 0 : else if (rq% kap_lowT_option == kap_lowT_Freedman11) then
166 0 : write(*,*) 'Failed to load lowT Freedman tables'
167 : else
168 0 : write(*,*) 'Failed to load lowT Fergusson tables'
169 : end if
170 0 : call mesa_error(__FILE__, __LINE__)
171 : end if
172 : end do
173 : end do
174 : end select
175 :
176 : end subroutine setup_lowT
177 :
178 :
179 : end subroutine setup_kap_tables
180 :
181 :
182 4600 : subroutine load_one(rq, &
183 : z_tables, iz, ix, X, Z, read_later, ierr)
184 : type (Kap_General_Info), pointer :: rq
185 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
186 : integer, intent(in) :: iz, ix
187 : real(dp), intent(in) :: X, Z
188 : logical, intent(in) :: read_later
189 : integer, intent(out) :: ierr
190 :
191 : character (len=256) :: fname, filename, cache_filename, temp_cache_filename
192 :
193 : ierr = 0
194 :
195 : call Get_Filenames(rq, z_tables, X, Z, &
196 2300 : kap_dir, fname, filename, cache_filename, temp_cache_filename, ierr)
197 2300 : if (ierr /= 0) return
198 :
199 2300 : if (rq% show_info) write(*,*) 'read filename <' // trim(filename) // '>'
200 :
201 : call Prepare_Kap_X_Table(rq, &
202 : z_tables, iz, z_tables(iz)% x_tables, ix, X, Z, &
203 : read_later, &
204 : fname, filename, cache_filename, temp_cache_filename,&
205 2300 : ierr)
206 2300 : if (ierr /= 0) return
207 :
208 2300 : if (.not. read_later) z_tables(iz)% x_tables(ix)% not_loaded_yet = .false.
209 :
210 : end subroutine load_one
211 :
212 :
213 4600 : subroutine Prepare_Kap_X_Table(rq, &
214 : z_tables, iz, x_tables, ix, X_in, Z_in, &
215 : read_later, &
216 : fname, filename, cache_filename, temp_cache_filename, &
217 : ierr)
218 : type (Kap_General_Info), pointer :: rq
219 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
220 : ! allocates the arrays and stores data
221 : type (Kap_X_Table), dimension(:), pointer :: x_tables
222 : integer, intent(in) :: iz, ix
223 : real(dp), intent(in) :: X_in, Z_in
224 : logical, intent(in) :: read_later
225 : character (*), intent(in) :: fname, filename, cache_filename, temp_cache_filename
226 : integer, intent(out) :: ierr
227 :
228 : integer :: io_unit, cache_io_unit
229 :
230 : real(dp) :: X, Z
231 : integer :: ios, form, version, n, num_logRs, num_logTs, status, nvec, i
232 2300 : real(dp) :: xin, zz, logR_min, logR_max, logT_min, logT_max, xErr, zErr
233 : character (len=1000) :: message
234 117300 : real(dp), target :: vec_ary(50)
235 : real(dp), pointer :: vec(:)
236 : real(dp), parameter :: tiny = 1d-6
237 : include 'formats'
238 :
239 2300 : ierr = 0
240 2300 : vec => vec_ary
241 2300 : X = X_in
242 2300 : Z = Z_in
243 2300 : nvec=-1
244 :
245 2300 : open(newunit=io_unit,file=trim(filename),action='read',status='old',iostat=ios)
246 2300 : if (ios /= 0) then
247 : !write(*,*) 'load kap tables ' // trim(filename)
248 0 : write(*,'(A)')
249 0 : write(*,'(A)')
250 0 : write(*,'(A)')
251 0 : write(*,'(A)')
252 0 : write(*,'(A)')
253 0 : write(*,*) 'NOTICE: missing kap data ' // trim(filename)
254 0 : write(*,*)
255 0 : write(*,*) 'Please check the validity of the kap_prefix string for this file.'
256 0 : write(*,*)
257 0 : write(*,*) 'If it is okay, you may need to install new kap data.'
258 0 : write(*,*) 'To do that, remove the directory mesa/data/kap_data,'
259 0 : write(*,*) 'and rerun the mesa ./install script.'
260 0 : write(*,'(A)')
261 0 : call mesa_error(__FILE__,__LINE__)
262 : end if
263 :
264 2300 : version = -1
265 :
266 2300 : read(io_unit,*,iostat=ierr) ! skip
267 2300 : if (ierr == 0) then
268 2300 : read(io_unit,*,iostat=ierr) ! skip
269 2300 : if (ierr == 0) then
270 2300 : read(io_unit,'(a)',iostat=ierr) message
271 2300 : if (ierr == 0) call str_to_vector(message, vec, nvec, ierr)
272 2300 : if (nvec < 10) ierr = -1
273 2300 : if (ierr == 0) then
274 : i = 0
275 2300 : i = i+1; form = int(vec(i))
276 2300 : i = i+1; version = int(vec(i))
277 2300 : i = i+1; xin = vec(i)
278 2300 : i = i+1; zz = vec(i)
279 2300 : i = i+1; num_logRs = int(vec(i))
280 2300 : i = i+1; logR_min = vec(i)
281 2300 : i = i+1; logR_max = vec(i)
282 2300 : i = i+1; num_logTs = int(vec(i))
283 2300 : i = i+1; logT_min = vec(i)
284 2300 : i = i+1; logT_max = vec(i)
285 : end if
286 : end if
287 : end if
288 :
289 :
290 :
291 2300 : if (ierr /= 0 .or. version < min_version) then
292 0 : write(*,*) 'load kap tables ' // trim(filename)
293 0 : write(*,'(A)')
294 0 : write(*,'(A)')
295 0 : write(*,'(A)')
296 0 : write(*,'(A)')
297 0 : write(*,'(A)')
298 0 : write(*,*) 'NOTICE: you need to install a new version of the kap data.'
299 0 : write(*,*) 'Please remove the directory mesa/data/kap_data,'
300 0 : write(*,*) 'and rerun the mesa ./install script.'
301 0 : write(*,'(A)')
302 0 : if (ierr /= 0) write(*,*) 'ierr', ierr
303 0 : if (version < min_version) &
304 0 : write(*,*) 'version < min_version', version, min_version
305 0 : write(*,*) 'form', form
306 0 : write(*,*) 'xin', xin
307 0 : write(*,*) 'zz', zz
308 0 : write(*,*) 'num_logRs', num_logRs
309 0 : write(*,*) 'logR_min', logR_min
310 0 : write(*,*) 'logR_max', logR_max
311 0 : write(*,*) 'num_logTs', num_logTs
312 0 : write(*,*) 'logT_min', logT_min
313 0 : write(*,*) 'logT_max', logT_max
314 0 : call mesa_error(__FILE__,__LINE__)
315 : end if
316 :
317 2300 : if (form /= kap_table_fixed_metal_form) then
318 0 : call mesa_error(__FILE__,__LINE__,'form /= kap_table_fixed_metal_form')
319 : end if
320 :
321 2300 : call Setup_Kap_X_Table(ierr)
322 4592 : if (ierr /= 0) return
323 :
324 2300 : if (read_later) then
325 2282 : close(io_unit)
326 2282 : return
327 : end if
328 :
329 18 : if (kap_use_cache) then
330 18 : ios = 0
331 : if (dbg_cache) then
332 : open(newunit=cache_io_unit,file=trim(cache_filename),action='read', &
333 : status='old',iostat=ios)
334 : else
335 : open(newunit=cache_io_unit,file=trim(cache_filename),action='read', &
336 18 : status='old',iostat=ios,form='unformatted')
337 : end if
338 18 : if (ios == 0) then ! try reading the cached data
339 : !write(*,'(a)') 'loading ' // trim(cache_filename)
340 10 : call Read_Kap_X_Table(cache_io_unit, .true., ierr)
341 10 : close(cache_io_unit)
342 10 : if (ierr == 0) then
343 10 : close(io_unit)
344 10 : return
345 : end if
346 : ierr = 0
347 : else
348 : !write(*,*) 'failed to open ' // trim(cache_filename)
349 : end if
350 : end if
351 :
352 : if (show_allocations) write(*,2) 'loading ' // trim(filename)
353 8 : call Read_Kap_X_Table(io_unit, .false., ierr)
354 8 : close(io_unit)
355 8 : if (ierr /= 0) then
356 0 : write(*,*) 'failed in Read_Kap_X_Table ' // trim(filename)
357 0 : call mesa_error(__FILE__,__LINE__)
358 0 : return
359 : end if
360 :
361 8 : if (.not. kap_use_cache) return
362 :
363 : if (dbg_cache) then
364 : open(newunit=cache_io_unit, file=trim(switch_str(temp_cache_filename, cache_filename, use_mesa_temp_cache)),&
365 : iostat=ios, action='write')
366 : else
367 : open(newunit=cache_io_unit, file=trim(switch_str(temp_cache_filename, cache_filename, use_mesa_temp_cache)),&
368 8 : iostat=ios, action='write', form='unformatted')
369 : end if
370 :
371 8 : if (ios == 0) then
372 8 : write(*,'(a)') 'write ' // trim(cache_filename)
373 : call Write_Kap_X_Table_Cache( &
374 8 : z_tables(iz)% x_tables, ix, cache_io_unit, version)
375 8 : close(cache_io_unit)
376 : ! Atomic move temp cache file to cache folder
377 8 : if(use_mesa_temp_cache) call mv(temp_cache_filename,cache_filename,.true.)
378 8 : if (kap_read_after_write_cache) then
379 : open(newunit=cache_io_unit,file=trim(cache_filename),action='read', &
380 8 : status='old',iostat=ios,form='unformatted')
381 8 : if (ios == 0) then
382 8 : call Read_Kap_X_Table(cache_io_unit, .true., ierr)
383 8 : close(cache_io_unit)
384 : end if
385 : end if
386 : end if
387 :
388 :
389 : contains
390 :
391 :
392 2300 : subroutine Setup_Kap_X_Table(ierr)
393 : integer, intent(out) :: ierr
394 :
395 2300 : xErr = abs(xin - X); zErr = abs(zz - Z)
396 2300 : if (xErr > tiny .or. zErr > tiny) then
397 0 : ierr = -1
398 0 : write(*,*) 'bug in file ' // trim(filename), xErr, zErr
399 0 : write(*,*) 'xErr > tiny', xErr > tiny
400 0 : write(*,*) 'zErr > tiny', zErr > tiny
401 0 : write(*,*) 'xin', xin
402 0 : write(*,*) 'X', X
403 0 : write(*,*) 'zz', zz
404 0 : write(*,*) 'Z', Z
405 0 : return
406 : end if
407 :
408 2300 : x_tables(ix)% not_loaded_yet = .true.
409 2300 : x_tables(ix)% X = X
410 2300 : x_tables(ix)% Z = Z
411 2300 : x_tables(ix)% logR_min = logR_min
412 2300 : x_tables(ix)% logR_max = logR_max
413 :
414 2300 : x_tables(ix)% num_logRs = num_logRs
415 2300 : nullify(x_tables(ix)% logRs)
416 :
417 2300 : x_tables(ix)% logT_min = logT_min
418 2300 : x_tables(ix)% logT_max = logT_max
419 2300 : x_tables(ix)% num_logTs = num_logTs
420 2300 : nullify(x_tables(ix)% logTs)
421 :
422 2300 : nullify(x_tables(ix)% kap1) ! allocate when read the data
423 :
424 : end subroutine Setup_Kap_X_Table
425 :
426 :
427 26 : subroutine Read_Kap_X_Table(io_unit, reading_cache, ierr)
428 : integer, intent(in) :: io_unit ! use this for file access
429 : logical, intent(in) :: reading_cache
430 : integer, intent(out) :: ierr
431 :
432 : character (len=1000) :: message
433 : character (len=1) :: char
434 : integer :: i, j, c_version, c_num_logRs, c_num_logTs
435 26 : real(dp) :: c_xin, c_zz, c_logR_min, c_logR_max, &
436 26 : c_logT_min, c_logT_max
437 1362 : real(dp) :: kap_logKs(num_logRs), logT
438 26 : real(dp), pointer :: kap(:,:,:), kap1(:)
439 2626 : real(dp), target :: vec_ary(100)
440 : real(dp), pointer :: vec(:)
441 : integer :: nvec
442 :
443 : include 'formats'
444 :
445 26 : vec => vec_ary
446 26 : nvec=-1
447 :
448 26 : if (reading_cache) then
449 :
450 18 : ios = 0
451 : if (dbg_cache) then
452 : write(*,*) 'io_unit', io_unit
453 : read(io_unit, *, iostat=ios) c_version, c_num_logRs, c_num_logTs
454 : else
455 18 : read(io_unit, iostat=ios) c_version, c_num_logRs, c_num_logTs
456 : end if
457 18 : if (ios /= 0 .or. c_version /= version) then
458 0 : ierr = 1
459 0 : if (ios /= 0) then
460 0 : write(*,*) 'cache failed in read 1'
461 0 : else if (c_version /= version) then
462 0 : write(*,*) 'cache failed for c_version /= version'
463 0 : write(*,*) 'c_version', c_version
464 0 : write(*,*) 'version', version
465 0 : else if (c_num_logRs /= num_logRs .or. c_num_logTs /= num_logTs) then
466 0 : if (c_num_logRs /= num_logRs) write(*,*) 'cache failed for c_num_logRs /= num_logRs'
467 0 : if (c_num_logTs /= num_logTs) write(*,*) 'cache failed for c_num_logTs /= num_logTs'
468 : end if
469 0 : return
470 : end if
471 :
472 : read(io_unit, iostat=ios) &
473 18 : c_xin, c_zz, c_logR_min, c_logR_max, &
474 36 : c_logT_min, c_logT_max
475 18 : if (ios /= 0) then
476 0 : ierr = 1
477 0 : if (ios /= 0) write(*,*) 'cache failed in read 2'
478 : end if
479 :
480 : end if
481 :
482 26 : xErr = abs(xin - X); zErr = abs(zz - Z)
483 26 : if (xErr > tiny .or. zErr > tiny) then
484 0 : if (reading_cache) then
485 0 : if (xErr > tiny) write(*,*) 'cache failed for xErr > tiny'
486 0 : if (zErr > tiny) write(*,*) 'cache failed for zErr > tiny'
487 0 : ierr = 1; return
488 : end if
489 0 : ierr = -1
490 0 : return
491 : end if
492 :
493 : if (show_allocations) write(*,2) 'x_tables ' // trim(filename), &
494 : num_logRs + num_logTs + sz_per_Kap_point*num_logRs*num_logTs
495 : allocate(x_tables(ix)% logRs(num_logRs), x_tables(ix)% logTs(num_logTs), &
496 26 : x_tables(ix)% kap1(sz_per_Kap_point*num_logRs*num_logTs), STAT=status)
497 26 : if (status /= 0) then
498 0 : ierr = -1
499 0 : return
500 : end if
501 :
502 26 : kap1 => x_tables(ix)% kap1
503 : kap(1:sz_per_Kap_point,1:num_logRs,1:num_logTs) => &
504 26 : kap1(1:sz_per_Kap_point*num_logRs*num_logTs)
505 :
506 26 : if (.not. reading_cache) then
507 :
508 8 : read(io_unit,*,iostat=ierr) ! skip
509 8 : if (ierr /= 0) return
510 8 : read(io_unit,*,iostat=ierr) ! skip
511 8 : if (ierr /= 0) return
512 :
513 8 : read(io_unit,'(a)',iostat=ierr) message
514 8 : if (ierr == 0) call str_to_vector(message, vec, nvec, ierr)
515 8 : if (nvec < num_logRs) ierr = -1
516 8 : if (ierr /= 0) return
517 420 : do j=1,num_logRs
518 420 : x_tables(ix)% logRs(j) = vec(j)
519 : !write(*,*) 'read logR', j, vec(j)
520 : end do
521 : !write(*,*) "input line: <" // trim(message) // ">"
522 :
523 8 : read(io_unit,*,iostat=ierr) ! skip
524 8 : if (ierr /= 0) return
525 :
526 1020 : do i = 1, num_logTs
527 1012 : read(io_unit,'(a)',iostat=ierr) message
528 1012 : if (ierr /= 0) then
529 0 : write(*,*) 'failed in read for logT i', i
530 0 : return
531 : !stop
532 : end if
533 1012 : call str_to_vector(message, vec, nvec, ierr)
534 1012 : if (ierr /= 0) then
535 0 : write(*,*) 'failed in str_to_vector'
536 0 : write(*,*) 'nvec', nvec
537 0 : write(*,'(a)') 'message "' // trim(message) // '"'
538 0 : return
539 : !stop
540 : end if
541 :
542 1012 : if (nvec < 1+num_logRs) ierr = -1
543 1012 : if (ierr /= 0) return
544 1012 : x_tables(ix)% logTs(i) = vec(1)
545 52596 : do j=1,num_logRs
546 51584 : kap_logKs(j) = vec(j+1)
547 52596 : kap(1,j,i) = kap_logKs(j)
548 : end do
549 :
550 1020 : if (.false.) then
551 : write(*,2) 'logT', i, x_tables(ix)% logTs(i)
552 : do j=1,num_logRs
553 : write(*,3) 'kap', j, i, kap(1,j,i)
554 : end do
555 : write(*,'(a)') 'message "' // trim(message) // '"'
556 : end if
557 :
558 : end do
559 :
560 : call Make_Interpolation_Data( &
561 : kap1, x_tables(ix)% logRs, num_logRs, &
562 : x_tables(ix)% logTs, num_logTs, &
563 8 : x_tables(ix)% ili_logRs, x_tables(ix)% ili_logTs, ierr)
564 :
565 8 : if (ierr /= 0) write(*,*) 'Read_Kap_X_Table failed in Make_Interpolation_Data'
566 :
567 : else ! reading_cache
568 :
569 : read(io_unit, iostat=ierr) &
570 18 : x_tables(ix)% ili_logRs, x_tables(ix)% ili_logTs
571 18 : if (ierr /= 0) return
572 :
573 : read(io_unit, iostat=ierr) &
574 916 : x_tables(ix)% logRs(1:num_logRs), &
575 2336 : x_tables(ix)% logTs(1:num_logTs)
576 18 : if (ierr /= 0) return
577 :
578 453538 : read(io_unit, iostat=ierr) kap1
579 18 : if (ierr /= 0) return
580 :
581 : end if
582 :
583 26 : end subroutine Read_Kap_X_Table
584 :
585 :
586 : end subroutine Prepare_Kap_X_Table
587 :
588 :
589 8 : subroutine Make_Interpolation_Data( &
590 : kap1, logRs, num_logRs, logTs, num_logTs, ili_logRs, ili_logTs, ierr)
591 : use interp_2d_lib_db
592 : real(dp), pointer :: kap1(:)
593 : integer, intent(in) :: num_logRs, num_logTs
594 : real(dp), intent(in), pointer :: logRs(:) ! (num_logRs)
595 : real(dp), intent(in), pointer :: logTs(:) ! (num_logTs)
596 : integer, intent(out) :: ili_logRs, ili_logTs, ierr
597 :
598 : character (len=256) :: message
599 : integer :: ibcxmin ! bc flag for x=xmin
600 1020 : real(dp) :: bcxmin(num_logTs) ! bc data vs. y at x=xmin
601 : integer :: ibcxmax ! bc flag for x=xmax
602 1020 : real(dp) :: bcxmax(num_logTs) ! bc data vs. y at x=xmax
603 : integer :: ibcymin ! bc flag for y=ymin
604 420 : real(dp) :: bcymin(num_logRs) ! bc data vs. x at y=ymin
605 : integer :: ibcymax ! bc flag for y=ymax
606 428 : real(dp) :: bcymax(num_logRs) ! bc data vs. x at y=ymax
607 : integer :: ier ! =0 on exit if there is no error.
608 : integer :: i, j
609 : real(dp) :: tol
610 : real(dp), pointer :: kap(:,:,:)
611 : real(dp), pointer :: x_out(:), y_out(:), f_out(:,:)
612 :
613 : kap(1:sz_per_kap_point,1:num_logRs,1:num_logTs) => &
614 8 : kap1(1:sz_per_kap_point*num_logRs*num_logTs)
615 :
616 : ! just use "not a knot" bc's at edges of tables
617 1020 : ibcxmin = 0; bcxmin(1:num_logTs) = 0d0
618 1020 : ibcxmax = 0; bcxmax(1:num_logTs) = 0d0
619 420 : ibcymin = 0; bcymin(1:num_logRs) = 0d0
620 420 : ibcymax = 0; bcymax(1:num_logRs) = 0d0
621 :
622 : call interp_mkbicub_db( &
623 : logRs, num_logRs, logTs, num_logTs, kap1, num_logRs, &
624 : ibcxmin,bcxmin,ibcxmax,bcxmax, &
625 : ibcymin,bcymin,ibcymax,bcymax, &
626 8 : ili_logRs,ili_logTs,ier)
627 :
628 8 : if (ier /= 0) then
629 0 : write(*,*) 'interp_mkbicub_db error happened for Make_Interpolation_Data for table'
630 0 : write(*,*) 'num_logRs', num_logRs
631 0 : do i=1,num_logRs
632 0 : write(*,*) 'logR', i, logRs(i)
633 : end do
634 0 : write(*,*) 'num_logTs', num_logTs
635 0 : do i=1,num_logTs
636 0 : write(*,*) 'logT', i, logTs(i)
637 : end do
638 0 : write(*,'(A)')
639 0 : call mesa_error(__FILE__,__LINE__,'kap interp error')
640 0 : ierr = -1
641 : return
642 : end if
643 :
644 8 : call Check_Interpolation_Data
645 :
646 8 : ierr = 0
647 :
648 : contains
649 :
650 8 : subroutine Check_Interpolation_Data
651 : use utils_lib,only:is_bad
652 : integer :: i, iR, jtemp
653 8 : real(dp) :: val
654 :
655 40 : do i = 1, sz_per_kap_point
656 1688 : do iR = 1, num_logRs
657 208016 : do jtemp = 1, num_logTs
658 206336 : val = kap(i,iR,jtemp)
659 207984 : if (is_bad(val)) then
660 : if (.true.) then
661 0 : write(*,*) 'bad value in xz', val, i, iR, jtemp
662 : write(*,'(99(a15,3x,f15.8,3x))') &
663 0 : 'logR', logRs(iR), 'logT', logTs(jtemp)
664 : end if
665 0 : kap(i,iR,jtemp) = 0
666 : end if
667 : end do
668 : end do
669 : end do
670 :
671 8 : end subroutine Check_Interpolation_Data
672 :
673 :
674 : end subroutine Make_Interpolation_Data
675 :
676 :
677 8 : subroutine Write_Kap_X_Table_Cache(x_tables, ix, io_unit, version)
678 : type (Kap_X_Table), dimension(:), pointer :: x_tables
679 : integer, intent(in) :: ix, io_unit, version
680 :
681 :
682 : if (dbg_cache) then
683 : write(*,*) 'write cache plain text', io_unit
684 : write(io_unit,*) version, x_tables(ix)% num_logTs, x_tables(ix)% num_logRs
685 : write(io_unit,'(999(1pe26.16))') &
686 : x_tables(ix)% X, &
687 : x_tables(ix)% Z, &
688 : x_tables(ix)% logR_min, &
689 : x_tables(ix)% logR_max, &
690 : x_tables(ix)% logT_min, &
691 : x_tables(ix)% logT_max
692 : !write(io_unit) &
693 : ! x_tables(ix)% logRs(:), &
694 : ! x_tables(ix)% logTs(:)
695 : !write(io_unit) x_tables(ix)% kap(:,:,:)
696 : else
697 8 : write(io_unit) version, x_tables(ix)% num_logTs, x_tables(ix)% num_logRs
698 : write(io_unit) &
699 8 : x_tables(ix)% X, &
700 8 : x_tables(ix)% Z, &
701 8 : x_tables(ix)% logR_min, &
702 8 : x_tables(ix)% logR_max, &
703 8 : x_tables(ix)% logT_min, &
704 16 : x_tables(ix)% logT_max
705 8 : write(io_unit) x_tables(ix)% ili_logRs, x_tables(ix)% ili_logTs
706 : write(io_unit) &
707 420 : x_tables(ix)% logRs(:), &
708 1028 : x_tables(ix)% logTs(:)
709 206344 : write(io_unit) x_tables(ix)% kap1(:)
710 : end if
711 :
712 8 : end subroutine Write_Kap_X_Table_Cache
713 :
714 :
715 2300 : subroutine Get_Filenames(rq, &
716 : z_tables, X, Z, data_dir, &
717 : fname, filename, cache_filename, temp_cache_filename, ierr)
718 : type (Kap_General_Info), pointer :: rq
719 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
720 : real(dp), intent(in) :: X, Z
721 : character (*), intent(in) :: data_dir
722 : character (*), intent(out) :: fname, filename, cache_filename, temp_cache_filename
723 : integer, intent(out) :: ierr
724 : character (len=256) :: cache_fname
725 : ierr=0
726 2300 : call create_fname(rq, z_tables, X, Z, fname, cache_fname,ierr)
727 2300 : filename = trim(data_dir) // '/' // fname
728 2300 : cache_filename = trim(kap_cache_dir) // '/' // cache_fname
729 2300 : temp_cache_filename=trim(kap_temp_cache_dir) // '/' // cache_fname
730 2300 : end subroutine Get_Filenames
731 :
732 :
733 : ! this must match the preprocessor naming scheme
734 6864 : subroutine create_fname(rq, z_tables, X, Z, fname, cache_fname, ierr)
735 : type (Kap_General_Info), pointer :: rq
736 : type (Kap_Z_Table), dimension(:), pointer :: z_tables
737 : real(dp), intent(in) :: X, Z
738 : character (len=*), intent(out) :: fname, cache_fname
739 : integer, intent(out) :: ierr
740 : character (len=256) :: zstr, xstr, prefix
741 :
742 : ierr=0
743 :
744 904 : if (z_tables(1)% lowT_flag .and. rq% kap_lowT_option == kap_lowT_Freedman11) then
745 18 : call get_output_string(Z,zstr,ierr)
746 18 : fname = trim(kap_lowT_option_str(rq% kap_lowT_option)) // '_z' // trim(zstr) // '.data'
747 18 : cache_fname = trim(kap_lowT_option_str(rq% kap_lowT_option))// '_z' // trim(zstr) // '.bin'
748 18 : return
749 : end if
750 :
751 2282 : call get_output_string(Z,zstr,ierr)
752 2282 : call get_output_string(X,xstr,ierr)
753 2282 : if (z_tables(1)% lowT_flag) then
754 886 : prefix = kap_lowT_option_str(rq% kap_lowT_option)
755 : else
756 1396 : prefix = kap_option_str(rq% kap_option)
757 : end if
758 :
759 : fname = &
760 2282 : trim(prefix) // '_z' // trim(zstr) // '_x' // trim(xstr) // '.data'
761 : cache_fname = &
762 2282 : trim(prefix) // '_z' // trim(zstr) // '_x' // trim(xstr) // '.bin'
763 2300 : end subroutine create_fname
764 :
765 : end module load_kap
|
