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_CO_kap
21 :
22 : use kap_def
23 : use math_lib
24 : use const_def, only: dp, use_mesa_temp_cache
25 : use utils_lib, only: mesa_error, mv, switch_str
26 :
27 : implicit none
28 :
29 : private
30 : public :: min_version
31 : public :: setup_kap_co_tables
32 : public :: get_dx_lookup
33 : public :: load_one_CO
34 :
35 : integer, parameter :: min_version = 37
36 : logical, parameter :: CO_dbg = .false.
37 :
38 : contains
39 :
40 11 : subroutine Setup_Kap_CO_Tables(rq, co_z_tables, use_cache, load_on_demand, ierr)
41 : type (Kap_General_Info), pointer :: rq
42 : type (Kap_CO_Z_Table), dimension(:), pointer :: co_z_tables
43 : logical, intent(in) :: use_cache, load_on_demand
44 : integer, intent(out) :: ierr
45 : integer :: iz, ix, CO_option
46 : include 'formats'
47 11 : ierr = 0
48 11 : CO_option = rq% kap_CO_option
49 11 : if (.not. associated(co_z_tables)) then
50 81 : allocate(co_z_tables(num_kap_CO_Zs(CO_option)), STAT=ierr)
51 9 : if (ierr /= 0) return
52 81 : do iz = 1, num_kap_CO_Zs(CO_option)
53 72 : co_z_tables(iz)% Zbase = kap_CO_Zs(iz, CO_option)
54 72 : co_z_tables(iz)% log10_Zbase = safe_log10(kap_CO_Zs(iz, CO_option))
55 72 : co_z_tables(iz)% Zfrac_C = -1
56 72 : co_z_tables(iz)% Zfrac_N = -1
57 72 : co_z_tables(iz)% Zfrac_O = -1
58 72 : co_z_tables(iz)% Zfrac_Ne = -1
59 432 : allocate(co_z_tables(iz)% x_tables(num_kap_CO_Xs(CO_option)), STAT=ierr)
60 81 : if (ierr /= 0) return
61 : end do
62 : if (show_allocations) write(*,2) 'setup co_z_tables', &
63 : num_kap_CO_Zs(CO_option) + num_kap_CO_Zs(CO_option)*num_kap_CO_Xs(CO_option)
64 : end if
65 99 : do iz = 1, num_kap_CO_Zs(CO_option)
66 539 : do ix = 1, num_kap_CO_Xs(CO_option)
67 440 : call load_one_CO(rq,co_z_tables,iz,ix,load_on_demand,ierr)
68 528 : if (ierr /= 0) return
69 : end do
70 : end do
71 : end subroutine Setup_Kap_CO_Tables
72 :
73 :
74 888 : subroutine load_one_CO(rq, co_z_tables, iz, ix, read_later, ierr)
75 : use utils_lib
76 : type (Kap_General_Info), pointer :: rq
77 : type (Kap_CO_Z_Table), dimension(:), pointer :: co_z_tables
78 : integer, intent(in) :: iz, ix
79 : logical, intent(in) :: read_later
80 : integer, intent(out) :: ierr
81 :
82 : character (len=256) :: fname, filename, cache_filename, temp_cache_filename
83 :
84 : ierr = 0
85 :
86 : call Get_CO_Filenames(rq, &
87 : kap_CO_Zs(iz, rq% kap_CO_option), kap_CO_Xs(ix, rq% kap_CO_option), kap_dir, fname, filename, &
88 444 : cache_filename, temp_cache_filename, ierr)
89 444 : if (ierr /= 0) return
90 :
91 444 : if (rq% show_info) write(*,*) 'read filename <' // trim(filename) // '>'
92 :
93 : call Prepare_Kap_CO_X_Table(rq, &
94 : iz, co_z_tables, co_z_tables(iz)% x_tables, &
95 : ix, kap_CO_Xs(ix, rq% kap_CO_option), kap_CO_Zs(iz, rq% kap_CO_option), &
96 : read_later, fname, filename, &
97 444 : cache_filename, temp_cache_filename, ierr)
98 444 : if (ierr /= 0) return
99 :
100 444 : if (.not. read_later) co_z_tables(iz)% x_tables(ix)% not_loaded_yet = .false.
101 :
102 : end subroutine load_one_CO
103 :
104 :
105 888 : subroutine Prepare_Kap_CO_X_Table(rq, &
106 : iz, co_z_tables, x_tables, ix, X_in, Z_in, read_later, &
107 : fname, filename, cache_filename, temp_cache_filename, ierr)
108 : ! allocates the arrays and stores data
109 : type (Kap_General_Info), pointer :: rq
110 : type (Kap_CO_Z_Table), dimension(:), pointer :: co_z_tables
111 : type (Kap_CO_X_Table), dimension(:), pointer :: x_tables
112 : integer, intent(in) :: iz, ix
113 : real(dp), intent(in) :: X_in, Z_in ! expected contents of the Kap file
114 : logical, intent(in) :: read_later
115 : character (len=*), intent(in) :: fname, filename, cache_filename, temp_cache_filename
116 : integer, intent(out) :: ierr
117 :
118 : integer :: io_unit, cache_io_unit
119 :
120 : real(dp) :: X, Z
121 : integer :: ios, form, version, n, num_tables, num_logRs, num_logTs, status
122 444 : real(dp) :: xin, zz, Zfrac_C, Zfrac_N, Zfrac_O, Zfrac_Ne, &
123 888 : logR_min, logR_max, logT_min, logT_max, xErr, zErr
124 444 : type (Kap_CO_Table), dimension(:), pointer :: co_tables
125 : integer :: num_dXC_gt_dXO ! the number of tables with dXC > dXO
126 : integer :: CO_table_numbers(num_kap_CO_dXs,num_kap_CO_dXs)
127 : integer :: next_dXO_table(max_num_CO_tables)
128 : integer :: next_dXC_table(max_num_CO_tables)
129 : character (len=256) :: message
130 13764 : real(dp), target :: vec_ary(30)
131 : real(dp), pointer :: vec(:)
132 : integer :: nvec
133 : real(dp), parameter :: tiny = 1d-6
134 : include 'formats'
135 :
136 444 : ierr = 0
137 444 : Zfrac_C = 0.d0
138 444 : Zfrac_N = 0.d0
139 444 : Zfrac_O = 0.d0
140 444 : Zfrac_Ne = 0.d0
141 :
142 444 : vec => vec_ary
143 444 : X = X_in
144 444 : Z = Z_in
145 444 : form = 0
146 444 : nvec=-1
147 :
148 444 : num_dXC_gt_dXO = -1
149 32412 : CO_table_numbers = -1
150 31524 : next_dXO_table = -1
151 31524 : next_dXC_table = -1
152 :
153 444 : open(newunit=io_unit,file=trim(filename),action='read',status='old',iostat=ios)
154 444 : if (ios /= 0) then
155 0 : write(*,'(A)')
156 0 : write(*,'(A)')
157 0 : write(*,'(A)')
158 0 : write(*,'(A)')
159 0 : write(*,*) 'NOTICE: missing kap data ' // trim(filename)
160 0 : write(*,*)
161 0 : write(*,*) 'Please check the validity of the kap_prefix string for this file.'
162 0 : write(*,*)
163 0 : write(*,*) 'If it is okay, you may need to install new kap data.'
164 0 : write(*,*) 'To do that, remove the directory mesa/data/kap_data,'
165 0 : write(*,*) 'and rerun the mesa ./install script.'
166 0 : write(*,'(A)')
167 0 : call mesa_error(__FILE__,__LINE__)
168 : end if
169 :
170 444 : version = -1
171 444 : read(io_unit,*,iostat=ierr) ! skip the 1st line
172 444 : if (ierr == 0) then
173 444 : read(io_unit,*,iostat=ierr) ! skip the 2nd line
174 444 : if (ierr == 0) then
175 444 : read(io_unit,'(a)',iostat=ierr) message
176 444 : if (ierr == 0) call str_to_vector(message, vec, nvec, ierr)
177 444 : if (nvec < 15) ierr = -1
178 444 : if (ierr == 0) then
179 444 : form = int(vec(1))
180 444 : version = int(vec(2))
181 444 : num_tables = int(vec(3))
182 444 : xin = vec(4)
183 444 : zz = vec(5)
184 444 : Zfrac_C = vec(6)
185 444 : Zfrac_N = vec(7)
186 444 : Zfrac_O = vec(8)
187 444 : Zfrac_Ne = vec(9)
188 444 : num_logRs = int(vec(10))
189 444 : logR_min = vec(11)
190 444 : logR_max = vec(12)
191 444 : num_logTs = int(vec(13))
192 444 : logT_min = vec(14)
193 444 : logT_max = vec(15)
194 : end if
195 : end if
196 : end if
197 :
198 444 : if (ierr /= 0 .or. version < min_version) then
199 0 : write(*,'(A)')
200 0 : write(*,'(A)')
201 0 : write(*,'(A)')
202 0 : write(*,'(A)')
203 0 : write(*,'(A)')
204 0 : write(*,*) 'FATAL ERROR: out-of-date version of kap data ' // trim(filename)
205 0 : write(*,*) 'Please update by removing the directory mesa/data/kap_data,'
206 0 : write(*,*) 'and rerunning the mesa ./install script.'
207 0 : write(*,'(A)')
208 0 : call mesa_error(__FILE__,__LINE__)
209 : end if
210 :
211 444 : if (form /= kap_table_co_enhanced_form) then
212 0 : call mesa_error(__FILE__,__LINE__,'form /= kap_table_co_enhanced_form')
213 : end if
214 :
215 444 : if (max_num_CO_tables < num_tables) then
216 0 : ierr = -1
217 442 : return
218 : end if
219 :
220 444 : if (Zfrac_C < 0 .or. Zfrac_N < 0 .or. Zfrac_O < 0 .or. Zfrac_Ne < 0) then
221 0 : ierr = -1
222 0 : return
223 : end if
224 :
225 444 : co_z_tables(iz)% Zfrac_C = Zfrac_C
226 444 : co_z_tables(iz)% Zfrac_N = Zfrac_N
227 444 : co_z_tables(iz)% Zfrac_O = Zfrac_O
228 444 : co_z_tables(iz)% Zfrac_Ne = Zfrac_Ne
229 :
230 444 : call Setup_Kap_CO_X_Table(ierr)
231 444 : if (ierr /= 0) return
232 :
233 : if (CO_dbg) write(*,*) 'after Setup_Kap_CO_X_Table: associated co_tables', &
234 : associated(x_tables(ix)% co_tables), ix
235 :
236 444 : if (read_later) then
237 440 : close(io_unit)
238 440 : return
239 : end if
240 :
241 4 : if (kap_use_cache) then
242 : open(newunit=cache_io_unit,file=trim(cache_filename),action='read', &
243 4 : status='old',iostat=ios,form='unformatted')
244 4 : if (ios == 0) then ! try reading the cached data
245 2 : call Read_Kap_CO_X_Table(cache_io_unit, .true., ierr)
246 2 : close(cache_io_unit)
247 2 : if (ierr == 0) then
248 2 : close(io_unit)
249 2 : return
250 : end if
251 : ierr = 0
252 : else
253 : if (CO_dbg) write(*,'(a)') 'failed to open ' // trim(cache_filename)
254 : end if
255 : end if
256 :
257 : if (CO_dbg) write(*,*) 'before call: associated co_tables', &
258 : associated(x_tables(ix)% co_tables), ix
259 : if (show_allocations) write(*,2) 'loading ' // trim(filename)
260 2 : call Read_Kap_CO_X_Table(io_unit, .false., ierr)
261 2 : close(io_unit)
262 2 : if (ierr /= 0) then
263 0 : write(*,*) 'failed in Read_Kap_CO_X_Table ' // trim(filename)
264 0 : return
265 : end if
266 :
267 2 : if (.not. kap_use_cache) return
268 :
269 : open(newunit=cache_io_unit, file=trim(switch_str(temp_cache_filename, cache_filename, use_mesa_temp_cache)),&
270 2 : iostat=ios, action='write', form='unformatted')
271 2 : if (ios == 0) then
272 2 : write(*,'(a)') 'write ' // trim(cache_filename)
273 : call Write_Kap_CO_X_Table_Cache( &
274 : x_tables, ix, cache_io_unit, x_tables(ix)% num_logRs, &
275 2 : x_tables(ix)% num_logTs, version)
276 2 : close(cache_io_unit)
277 2 : if(use_mesa_temp_cache) call mv(temp_cache_filename, cache_filename,.true.)
278 2 : if (kap_read_after_write_cache) then
279 : open(newunit=cache_io_unit,file=trim(cache_filename),action='read', &
280 2 : status='old',iostat=ios,form='unformatted')
281 2 : if (ios == 0) then
282 2 : call Read_Kap_CO_X_Table(cache_io_unit, .true., ierr)
283 2 : close(cache_io_unit)
284 : else
285 0 : ierr = -1
286 : end if
287 : end if
288 : end if
289 :
290 :
291 : contains
292 :
293 :
294 444 : subroutine Setup_Kap_CO_X_Table(ierr)
295 : integer, intent(out) :: ierr
296 :
297 : integer :: i
298 : include 'formats'
299 :
300 444 : xErr = abs(xin - X); zErr = abs(zz - Z)
301 444 : if (num_tables <= 0 .or. xErr > tiny .or. zErr > tiny) then
302 0 : ierr = -1
303 0 : write(*,*) 'bug in file ' // trim(filename), num_tables, xErr, zErr
304 0 : write(*,*) 'num_tables <= 0', num_tables <= 0
305 0 : write(*,*) 'xErr > tiny', xErr > tiny
306 0 : write(*,*) 'zErr > tiny', zErr > tiny
307 0 : return
308 : end if
309 :
310 444 : x_tables(ix)% not_loaded_yet = .true.
311 :
312 444 : x_tables(ix)% X = X
313 444 : x_tables(ix)% Z = Z
314 :
315 444 : x_tables(ix)% logR_min = logR_min
316 444 : x_tables(ix)% logR_max = logR_max
317 444 : x_tables(ix)% num_logRs = num_logRs
318 444 : nullify(x_tables(ix)% logRs)
319 :
320 444 : x_tables(ix)% logT_min = logT_min
321 444 : x_tables(ix)% logT_max = logT_max
322 444 : x_tables(ix)% num_logTs = num_logTs
323 444 : nullify(x_tables(ix)% logTs)
324 :
325 444 : x_tables(ix)% num_CO_tables = num_tables
326 :
327 : if (show_allocations) write(*,2) 'co_tables', num_tables
328 22764 : allocate(co_tables(num_tables), STAT=status)
329 444 : if (status /= 0) then
330 0 : ierr = -1
331 : if (CO_dbg) write(*,*) 'InsufficientMemory for Prepare_Kap_CO_X_Table', iz, ix
332 0 : return
333 : end if
334 444 : x_tables(ix)% co_tables => co_tables
335 : if (CO_dbg) write(*,*) 'Setup_Kap_CO_X_Table: allocate co_tables', iz, ix, num_tables, status
336 : if (CO_dbg) write(*,*) 'Setup_Kap_CO_X_Table: associated(x_tables(ix)% co_tables)', &
337 : associated(x_tables(ix)% co_tables), associated(co_tables), ix
338 :
339 22764 : do i=1,num_tables
340 22764 : nullify(co_tables(i)% kap1) ! allocate when read the data
341 : end do
342 :
343 : end subroutine Setup_Kap_CO_X_Table
344 :
345 :
346 6 : subroutine Read_Kap_CO_X_Table(io_unit, reading_cache, ierr)
347 : integer, intent(in) :: io_unit ! use this for file access
348 : logical, intent(in) :: reading_cache
349 : integer, intent(out) :: ierr
350 :
351 : character (len=256) :: message
352 : character (len=1) :: char
353 : integer :: c_num_tables, c_num_logRs, c_num_logTs, c_version
354 6 : real(dp) :: c_xin, c_zz, c_logR_min, c_logR_max, c_logT_min, c_logT_max
355 :
356 : include 'formats'
357 6 : ierr = 0
358 :
359 6 : if (reading_cache) then
360 :
361 4 : ios = 0
362 4 : read(io_unit, iostat=ios) c_version
363 4 : if (ios /= 0 .or. c_version /= version) then
364 0 : ierr = 1
365 0 : if (ios /= 0) write(*,*) 'cache failed in read'
366 0 : if (c_version /= version) write(*,*) 'cache failed for c_version /= version'
367 0 : return
368 : end if
369 :
370 : read(io_unit, iostat=ios) &
371 4 : c_xin, c_zz, &
372 4 : c_logR_min, c_logR_max, c_num_logRs, x_tables(ix)% ili_logRs, &
373 8 : c_logT_min, c_logT_max, c_num_logTs, x_tables(ix)% ili_logTs
374 4 : if (ios /= 0 .or. c_num_logRs /= num_logRs .or. c_num_logTs /= num_logTs) then
375 0 : ierr = 1
376 0 : if (ios /= 0) write(*,*) 'cache failed in read'
377 0 : if (c_num_logRs /= num_logRs) write(*,*) 'cache failed for c_num_logRs /= num_logRs'
378 0 : if (c_num_logTs /= num_logTs) write(*,*) 'cache failed for c_num_logTs /= num_logTs'
379 0 : return
380 : end if
381 :
382 : read(io_unit, iostat=ios) &
383 4 : c_num_tables, num_dXC_gt_dXO, &
384 8 : CO_table_numbers, next_dXO_table, next_dXC_table
385 :
386 4 : if (ios /= 0) then
387 0 : ierr = 1
388 0 : if (ios /= 0) write(*,*) 'cache failed in read'
389 0 : return
390 : end if
391 :
392 : end if
393 :
394 6 : xErr = abs(xin - X); zErr = abs(zz - Z)
395 6 : if (num_tables <= 0 .or. xErr > tiny .or. zErr > tiny) then
396 0 : if (reading_cache) then
397 0 : if (num_tables <= 0) write(*,*) 'cache failed for num_tables <= 0'
398 0 : if (xErr > tiny) write(*,*) 'cache failed for xErr > tiny'
399 0 : if (zErr > tiny) write(*,*) 'cache failed for zErr > tiny'
400 0 : ierr = 1; return
401 : end if
402 0 : ierr = -1
403 0 : write(*,*) 'num_tables <= 0', num_tables <= 0
404 0 : write(*,*) 'xErr > tiny', xErr > tiny
405 0 : write(*,*) 'zErr > tiny', zErr > tiny
406 0 : return
407 : end if
408 :
409 : if (show_allocations) write(*,2) 'Read_Kap_CO_X_Table logRs logTs', &
410 : num_logRs + num_logTs
411 6 : allocate(x_tables(ix)% logRs(num_logRs), x_tables(ix)% logTs(num_logTs), STAT=status)
412 6 : if (status /= 0) then
413 0 : ierr = -1
414 0 : write(*,*) 'InsufficientMemory for Read_Kap_CO_X_Table'
415 0 : return
416 : end if
417 :
418 6 : if (.not. reading_cache) then
419 :
420 2 : read(io_unit,*,iostat=ierr) ! skip line
421 2 : if (ierr /= 0) return
422 2 : read(io_unit,'(i3)',iostat=ierr) num_dXC_gt_dXO
423 2 : if (ierr /= 0) return
424 2 : x_tables(ix)% num_dXC_gt_dXO = num_dXC_gt_dXO
425 : do ! skip to the start of the first table
426 134 : read(io_unit,'(a1)',iostat=ierr) char
427 134 : if (ierr /= 0) return
428 134 : if (char == '-') exit
429 : end do
430 2 : read(io_unit,*,iostat=ierr) ! skip line
431 2 : if (ierr /= 0) return
432 146 : x_tables(ix)% CO_table_numbers = -1
433 142 : x_tables(ix)% next_dXC_table = -1
434 142 : x_tables(ix)% next_dXO_table = -1
435 :
436 : else
437 :
438 : read(io_unit, iostat=ierr) &
439 152 : x_tables(ix)% logRs(1:num_logRs), &
440 560 : x_tables(ix)% logTs(1:num_logTs)
441 4 : if (ierr /= 0) return
442 :
443 4 : x_tables(ix)% num_dXC_gt_dXO = num_dXC_gt_dXO
444 292 : x_tables(ix)% CO_table_numbers = CO_table_numbers
445 284 : x_tables(ix)% next_dXO_table = next_dXO_table
446 284 : x_tables(ix)% next_dXC_table = next_dXC_table
447 :
448 : end if
449 :
450 354 : do n = 1, num_tables
451 : if (CO_dbg) write(*,*) 'call Read_Kap_CO_Table', ix, n, X, Z
452 : call Read_Kap_CO_Table(x_tables, ix, n, X, Z, &
453 348 : num_logRs, num_logTs, io_unit, reading_cache, ierr)
454 354 : if (ierr /= 0) return
455 : end do
456 :
457 6 : if (.not. reading_cache) then
458 118 : do n = 1, num_tables
459 116 : next_dXC_table(n) = find_next_dXC_table(n)
460 118 : next_dXO_table(n) = find_next_dXO_table(n)
461 : end do
462 142 : x_tables(ix)% next_dXO_table = next_dXO_table
463 142 : x_tables(ix)% next_dXC_table = next_dXC_table
464 : end if
465 :
466 : end subroutine Read_Kap_CO_X_Table
467 :
468 :
469 116 : integer function find_next_dXC_table(i)
470 : integer, intent(in) :: i
471 : integer :: j
472 116 : real(dp) :: dXC, dXO
473 116 : dXC = co_tables(i)% dXC
474 116 : dXO = co_tables(i)% dXO
475 116 : if (dXC > dXO) then
476 52 : find_next_dXC_table = i + 1; return
477 : end if
478 2042 : do j = 1, num_tables
479 2042 : if (co_tables(j)% dXC > dXC .and. co_tables(j)% dXO == dXO) then
480 50 : find_next_dXC_table = j; return
481 : end if
482 : end do
483 116 : find_next_dXC_table = -1
484 : end function find_next_dXC_table
485 :
486 :
487 116 : integer function find_next_dXO_table(i)
488 : integer, intent(in) :: i
489 : integer :: j
490 116 : real(dp) :: dXC, dXO
491 116 : dXC = co_tables(i)% dXC
492 116 : dXO = co_tables(i)% dXO
493 116 : if (dXC < dXO) then
494 52 : find_next_dXO_table = i + 1; return
495 : end if
496 2130 : do j = 1, num_tables
497 2130 : if (co_tables(j)% dXO > dXO .and. co_tables(j)% dXC == dXC) then
498 50 : find_next_dXO_table = j
499 50 : return
500 : end if
501 : end do
502 116 : find_next_dXO_table = -1
503 : end function find_next_dXO_table
504 :
505 :
506 : end subroutine Prepare_Kap_CO_X_Table
507 :
508 :
509 348 : subroutine Read_Kap_CO_Table( &
510 : x_tables, ix, n, X_in, Z_in, num_logRs, num_logTs, io_unit, reading_cache, ierr)
511 : type (Kap_CO_X_Table), dimension(:), pointer :: x_tables
512 : integer, intent(in) :: ix ! index in x_tables
513 : integer, intent(in) :: n ! index in co_tables
514 : real(dp), intent(in) :: X_in, Z_in
515 : integer, intent(in) :: io_unit ! use this for file access
516 : logical, intent(in) :: reading_cache
517 : integer, intent(in) :: num_logRs, num_logTs
518 : integer, intent(out) :: ierr ! return nonzero if had trouble
519 :
520 : type (Kap_CO_Table), pointer :: co_tables(:)
521 : integer :: table_num, i, j, ios, status, iXC, iXO
522 13920 : real(dp) :: X, Z, xin, zz, Y, dXC, dXO, err, logT
523 : real(dp), allocatable, target :: kap_table(:) ! data & spline coefficients
524 348 : real(dp), pointer :: kap(:,:,:)
525 26100 : real(dp) :: logKs(num_logRs), logRs(num_logRs)
526 : character (len=1000) :: message
527 17748 : real(dp), target :: vec_ary(50)
528 : real(dp), pointer :: vec(:)
529 : integer :: nvec
530 : real(dp), parameter :: tiny = 1d-6
531 :
532 : logical :: store_logRs, store_logTs
533 : include 'formats'
534 :
535 348 : ierr = 0
536 348 : vec => vec_ary
537 348 : X = X_in; Z = Z_in
538 348 : nvec=-1
539 :
540 : if (show_allocations) write(*,2) 'Read_Kap_CO_Table', &
541 : sz_per_Kap_point*num_logRs*num_logTs
542 348 : allocate(kap_table(sz_per_Kap_point*num_logRs*num_logTs))
543 :
544 348 : co_tables => x_tables(ix)% co_tables
545 : if (CO_dbg) write(*,*) 'Read_Kap_CO_Table associated(co_tables)', &
546 : associated(co_tables), ix, reading_cache
547 :
548 348 : if (.not. reading_cache) then
549 116 : read(io_unit,*,iostat=ierr)
550 116 : if (ierr /= 0) return
551 116 : read(io_unit,'(a)',iostat=ierr) message
552 116 : if (ierr == 0) call str_to_vector(message, vec, nvec, ierr)
553 116 : if (ierr /= 0 .or. nvec < 6) then
554 0 : write(*,*) 'Read_Kap_CO_Table str_to_vector'
555 0 : ierr = 1
556 0 : return
557 : end if
558 116 : table_num = int(vec(1))
559 116 : xin = vec(2)
560 116 : Y = vec(3)
561 116 : zz = vec(4)
562 116 : dXC = vec(5)
563 116 : dXO = vec(6)
564 116 : if (table_num /= n) then
565 0 : write(*,*) 'wrong num in opacity file for X', X, 'Z', Z, 'table', n
566 0 : ierr = -1
567 0 : return
568 : end if
569 : else ! reading_cache
570 232 : read(io_unit, iostat=ios) table_num, xin, Y, zz, dXC, dXO
571 232 : if (ios /= 0) then
572 0 : ierr = 1
573 0 : return
574 : end if
575 : end if
576 :
577 348 : if (abs(xin-X) > tiny .or. abs(zz-Z) > tiny) then
578 0 : if (reading_cache) then
579 0 : ierr = 1
580 0 : return
581 : end if
582 0 : write(*,*) 'error in opacity file for X', X, 'Z', Z, 'table', n
583 0 : ierr = -1
584 0 : return
585 : end if
586 348 : X = xin; Z = zz ! use the values from the file
587 :
588 348 : err = abs(1d0 - (X + Y + Z + dXC + dXO))
589 348 : if (err > tiny) then
590 0 : if (reading_cache) then
591 0 : ierr = 1; return
592 : end if
593 0 : write(*,*) 'abundance error in opacity table for X=', &
594 0 : X, 'Z=', Z, 'dXC=', dXC, 'dXO=', dXO
595 0 : ierr = -1
596 0 : return
597 : end if
598 :
599 348 : co_tables(n)% table_num = table_num
600 348 : co_tables(n)% X = X
601 348 : co_tables(n)% Z = Z
602 348 : co_tables(n)% dXC = dXC
603 348 : co_tables(n)% dXO = dXO
604 348 : co_tables(n)% dXC_lookup = get_dX_lookup(dXC, Z)
605 348 : co_tables(n)% dXO_lookup = get_dX_lookup(dXO, Z)
606 :
607 : if (show_allocations) write(*,2) 'co_tables', &
608 : sz_per_Kap_point*num_logRs*num_logTs
609 348 : allocate(co_tables(n)% kap1(sz_per_Kap_point*num_logRs*num_logTs), STAT=status)
610 348 : if (status /= 0) then
611 0 : ierr = -1
612 0 : return
613 : end if
614 :
615 348 : if (reading_cache) then
616 :
617 232 : read(io_unit, iostat=ios) kap_table
618 4738600 : do i=1,sz_per_Kap_point*num_logRs*num_logTs
619 4738600 : co_tables(n)% kap1(i) = kap_table(i)
620 : end do
621 232 : if (ios /= 0) then
622 0 : ierr = 4; return
623 : end if
624 :
625 : else
626 :
627 116 : iXC = find_in_dXs(dXC)
628 116 : if (iXC > 0) then
629 104 : iXO = find_in_dXs(dXO)
630 104 : if (iXO > 0) then
631 92 : x_tables(ix)% CO_table_numbers(iXC,iXO) = n
632 : end if
633 : end if
634 :
635 116 : read(io_unit,*,iostat=ierr)
636 116 : if (ierr /= 0) return
637 116 : read(io_unit,*,iostat=ierr)
638 116 : if (ierr /= 0) return
639 :
640 116 : read(io_unit,'(a)',iostat=ierr) message
641 116 : if (ierr == 0) call str_to_vector(message, vec, nvec, ierr)
642 116 : if (ierr /= 0) write(*,*) 'str_to_vector ierr', ierr
643 116 : if (nvec < num_logRs) ierr = -1
644 116 : if (ierr /= 0) then
645 0 : write(*,*) 'Read_Kap_CO_Table str_to_vector logRs: nvec, num_logRs', nvec, num_logRs
646 0 : write(*,'(a)') trim(message)
647 0 : stop
648 : return
649 : end if
650 4408 : do i=1,num_logRs
651 4408 : logRs(i) = vec(i)
652 : end do
653 :
654 116 : if (n == 1) then
655 76 : x_tables(ix)% logRs(:) = logRs(:)
656 : else ! check
657 4332 : do i=1, num_logRs
658 4332 : if (abs(logRs(i) - x_tables(ix)% logRs(i)) > 1d-7) then
659 0 : write(*,*) 'problem with inconsistent logRs'
660 0 : call mesa_error(__FILE__,__LINE__)
661 : end if
662 : end do
663 : end if
664 :
665 116 : read(io_unit,*,iostat=ierr)
666 116 : if (ierr /= 0) return
667 :
668 : kap(1:sz_per_Kap_point,1:num_logRs,1:num_logTs) => &
669 116 : co_tables(n)% kap1(1:sz_per_Kap_point*num_logRs*num_logTs)
670 :
671 16124 : do i = 1, num_logTs
672 :
673 16008 : read(io_unit,'(a)',iostat=ierr) message
674 16008 : if (ierr == 0) call str_to_vector(message, vec, nvec, ierr)
675 16008 : if (nvec < 1+num_logRs) ierr = -1
676 16008 : if (ierr /= 0) then
677 0 : write(*,*) 'Read_Kap_CO_Table str_to_vector logKs'
678 0 : return
679 : end if
680 16008 : logT = vec(1)
681 608304 : do j=1,num_logRs
682 592296 : logKs(j) = vec(j+1)
683 608304 : kap(1,j,i) = logKs(j)
684 : end do
685 :
686 16124 : if (n == 1) then
687 276 : x_tables(ix)% logTs(i) = logT
688 : else ! check
689 15732 : if (abs(logT - x_tables(ix)% logTs(i)) > 1d-7) then
690 0 : write(*,*) 'problem with inconsistent logTs'
691 0 : call mesa_error(__FILE__,__LINE__)
692 : end if
693 : end if
694 :
695 : end do
696 116 : read(io_unit,*,iostat=ierr)
697 116 : if (ierr /= 0) return
698 :
699 : call Make_CO_Interpolation_Data( &
700 : co_tables, n, x_tables(ix)% logRs, num_logRs, &
701 : x_tables(ix)% logTs, num_logTs, &
702 116 : x_tables(ix)% ili_logRs, x_tables(ix)% ili_logTs, ierr)
703 :
704 116 : if (ierr /= 0) then
705 0 : write(*,*) 'call Make_CO_Interpolation_Data ierr', ierr
706 0 : write(*,*) 'n', n
707 0 : write(*,*) 'num_logRs', num_logRs
708 0 : write(*,*) 'num_logTs', num_logTs
709 : end if
710 :
711 : end if
712 :
713 696 : end subroutine Read_Kap_CO_Table
714 :
715 :
716 116 : subroutine Make_CO_Interpolation_Data( &
717 : co_tables, n, logRs, num_logRs, logTs, num_logTs, ili_logRs, ili_logTs, ierr)
718 : use interp_2d_lib_db
719 : type (Kap_CO_Table), dimension(:), pointer :: co_tables
720 : integer, intent(in) :: n, num_logRs, num_logTs
721 : real(dp), intent(in), pointer :: logRs(:) ! =(num_logRs)
722 : real(dp), intent(in), pointer :: logTs(:) ! =(num_logTs)
723 : integer, intent(out) :: ili_logRs, ili_logTs, ierr
724 :
725 2369300 : real(dp), target :: table_ary(sz_per_kap_point*num_logRs*num_logTs)
726 116 : real(dp), pointer :: table(:,:,:), table1(:), kap(:,:,:)
727 : character (len=256) :: message
728 : integer :: ibcxmin ! bc flag for x=xmin
729 16124 : real(dp) :: bcxmin(num_logTs) ! bc data vs. y at x=xmin
730 : integer :: ibcxmax ! bc flag for x=xmax
731 16124 : real(dp) :: bcxmax(num_logTs) ! bc data vs. y at x=xmax
732 : integer :: ibcymin ! bc flag for y=ymin
733 4408 : real(dp) :: bcymin(num_logRs) ! bc data vs. x at y=ymin
734 : integer :: ibcymax ! bc flag for y=ymax
735 4408 : real(dp) :: bcymax(num_logRs) ! bc data vs. x at y=ymax
736 : integer :: ier ! =0 on exit if there is no error.
737 : integer :: i, j
738 : real(dp) :: tol
739 : real(dp), pointer :: x_out(:), y_out(:), f_out(:,:)
740 :
741 : ! just use "not a knot" bc's at edges of tables
742 16124 : ibcxmin = 0; bcxmin(1:num_logTs) = 0d0
743 16124 : ibcxmax = 0; bcxmax(1:num_logTs) = 0d0
744 4408 : ibcymin = 0; bcymin(1:num_logRs) = 0d0
745 4408 : ibcymax = 0; bcymax(1:num_logRs) = 0d0
746 :
747 116 : table1 => table_ary
748 : table(1:sz_per_kap_point,1:num_logRs,1:num_logTs) => &
749 116 : table_ary(1:sz_per_kap_point*num_logRs*num_logTs)
750 : kap(1:sz_per_kap_point,1:num_logRs,1:num_logTs) => &
751 116 : co_tables(n)% kap1(1:sz_per_kap_point*num_logRs*num_logTs)
752 :
753 16124 : do j=1,num_logTs
754 608420 : do i=1,num_logRs
755 608304 : table(1,i,j) = kap(1,i,j)
756 : end do
757 : end do
758 :
759 : call interp_mkbicub_db( &
760 : logRs, num_logRs, logTs, num_logTs, table1, num_logRs, &
761 : ibcxmin,bcxmin,ibcxmax,bcxmax, &
762 : ibcymin,bcymin,ibcymax,bcymax, &
763 116 : ili_logRs,ili_logTs,ier)
764 :
765 116 : if (ier /= 0) then
766 0 : write(*,*) 'interp_mkbicub_db error happened for Make_CO_Interpolation_Data for table', n
767 0 : ierr = -1
768 : return
769 : end if
770 :
771 116 : call Check_Interpolation_Data
772 :
773 2369300 : do i=1,sz_per_kap_point*num_logRs*num_logTs
774 2369300 : co_tables(n)% kap1(i) = table1(i)
775 : end do
776 :
777 232 : ierr = 0
778 :
779 :
780 : contains
781 :
782 116 : subroutine Check_Interpolation_Data
783 : use utils_lib,only:is_bad
784 : integer :: i, iR, jtemp
785 116 : real(dp) :: val
786 :
787 580 : do i = 1, sz_per_kap_point
788 17748 : do iR = 1, num_logRs
789 2386816 : do jtemp = 1, num_logTs
790 2369184 : val = table(i,iR,jtemp)
791 2386352 : if (is_bad(val)) then
792 : if (.true.) then
793 0 : write(*,*) 'bad value in xz', val, i, iR, jtemp
794 : write(*,'(99(a15,3x,f15.8,3x))') &
795 0 : 'logR', logRs(iR), 'logT', logTs(jtemp)
796 : end if
797 0 : table(i,iR,jtemp) = 0
798 : end if
799 : end do
800 : end do
801 : end do
802 :
803 116 : end subroutine Check_Interpolation_Data
804 :
805 :
806 : end subroutine Make_CO_Interpolation_Data
807 :
808 :
809 2 : subroutine Write_Kap_CO_X_Table_Cache(x_tables, ix, io_unit, num_logRs, num_logTs, version)
810 : type (Kap_CO_X_Table), dimension(:), pointer :: x_tables
811 : integer, intent(in) :: ix, io_unit, num_logRs, num_logTs, version
812 :
813 2 : type (Kap_CO_Table), dimension(:), pointer :: co_tables
814 :
815 : integer :: num_tables, n, i
816 2 : real(dp) :: X, Z, Y, dXC, dXO
817 : real(dp), allocatable, target :: kap_table(:) ! data & spline coefficients
818 : include 'formats'
819 :
820 : if (show_allocations) write(*,2) 'Write_Kap_CO_X_Table_Cache', &
821 : sz_per_Kap_point*num_logRs*num_logTs
822 2 : allocate(kap_table(sz_per_Kap_point*num_logRs*num_logTs))
823 :
824 2 : num_tables = x_tables(ix)% num_CO_tables
825 2 : X = x_tables(ix)% X
826 2 : Z = x_tables(ix)% Z
827 2 : co_tables => x_tables(ix)% co_tables
828 :
829 2 : write(io_unit) version
830 :
831 : write(io_unit) &
832 2 : X, Z, &
833 2 : x_tables(ix)% logR_min, &
834 2 : x_tables(ix)% logR_max, &
835 2 : num_logRs, x_tables(ix)% ili_logRs, &
836 2 : x_tables(ix)% logT_min, &
837 2 : x_tables(ix)% logT_max, &
838 4 : num_logTs, x_tables(ix)% ili_logTs
839 :
840 : write(io_unit) &
841 2 : num_tables, x_tables(ix)% num_dXC_gt_dXO, &
842 146 : x_tables(ix)% CO_table_numbers(:,:), &
843 142 : x_tables(ix)% next_dXO_table(:), &
844 144 : x_tables(ix)% next_dXC_table(:)
845 :
846 : write(io_unit) &
847 76 : x_tables(ix)% logRs(1:num_logRs), &
848 280 : x_tables(ix)% logTs(1:num_logTs)
849 :
850 118 : do n = 1, num_tables
851 116 : dXC = co_tables(n)% dXC
852 116 : dXO = co_tables(n)% dXO
853 116 : Y = 1 - (X + Z + dXC + dXO)
854 2369300 : do i=1,sz_per_Kap_point*num_logRs*num_logTs
855 2369300 : kap_table(i) = co_tables(n)% kap1(i)
856 : end do
857 116 : write(io_unit) co_tables(n)% table_num, X, Y, Z, dXC, dXO
858 118 : write(io_unit) kap_table
859 : end do
860 :
861 4 : end subroutine Write_Kap_CO_X_Table_Cache
862 :
863 :
864 356 : real(dp) function get_dX_lookup(dX, Z)
865 : real(dp), intent(in) :: dX, Z
866 348 : get_dX_lookup = log10(Z + 1d-3 + dX)
867 8 : end function get_dX_lookup
868 :
869 :
870 220 : integer function find_in_dXs(dX)
871 : real(dp), intent(in) :: dX
872 : integer :: i
873 : real(dp), parameter :: tiny = 1d-6
874 962 : do i = 1, num_kap_CO_dXs
875 962 : if (abs(dX-kap_CO_dXs(i)) < tiny) then
876 196 : find_in_dXs = i; return
877 : end if
878 : end do
879 220 : find_in_dXs = -1
880 : end function find_in_dXs
881 :
882 :
883 444 : subroutine Get_CO_Filenames(rq, &
884 : Z, X, data_dir, fname, filename, cache_filename, temp_cache_filename, ierr)
885 : type (Kap_General_Info), pointer :: rq
886 : real(dp), intent(in) :: Z, X
887 : character (*), intent(in) :: data_dir
888 : character (*), intent(out) :: fname, filename, cache_filename, temp_cache_filename
889 : integer, intent(out) :: ierr
890 : ierr=0
891 444 : call Create_fname(rq, Z, X, fname,ierr)
892 444 : filename = trim(data_dir) // '/' // trim(fname) // '.data'
893 444 : cache_filename = trim(kap_cache_dir) // '/' // trim(fname) // '.bin'
894 444 : temp_cache_filename = trim(kap_temp_cache_dir) // '/' // trim(fname) // '.bin'
895 444 : end subroutine Get_CO_Filenames
896 :
897 :
898 888 : subroutine Create_fname(rq, Z, X, fname,ierr)
899 : type (Kap_General_Info), pointer :: rq
900 : real(dp), intent(in) :: Z, X
901 : character (len=*),intent(out) :: fname
902 : integer, intent(out) :: ierr
903 : character (len=16) :: zstr, xstr
904 : ierr=0
905 444 : call get_output_string(Z, zstr,ierr)
906 444 : call get_output_string(X, xstr,ierr)
907 :
908 : write(fname,'(a)') trim(kap_CO_option_str(rq% kap_CO_option)) // '_z' // &
909 444 : trim(zstr) // '_x' // trim(xstr)
910 :
911 444 : end subroutine Create_fname
912 :
913 : end module load_CO_kap
|
