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_weak
21 :
22 : use rates_def
23 : use const_def, only: dp
24 : use utils_lib, only: mesa_error
25 : use weaklib_tables, only: weaklib_rate_table
26 : use suzuki_tables, only: private_load_suzuki_tables
27 :
28 : implicit none
29 :
30 : private :: private_load_weak_tables
31 :
32 : contains
33 :
34 :
35 15 : subroutine load_weak_data(ierr)
36 : integer, intent(out) :: ierr
37 : ierr = 0
38 5 : call private_load_weak_tables(ierr)
39 5 : if (ierr /= 0) return
40 :
41 5 : call load_user_weak_tables(ierr)
42 5 : if (ierr /= 0) return
43 5 : if (use_suzuki_tables) then
44 2 : call private_load_suzuki_tables(ierr)
45 2 : if (ierr /= 0) return
46 : end if
47 :
48 5 : call load_weak_info_list(ierr)
49 : end subroutine load_weak_data
50 :
51 :
52 5 : subroutine load_weak_info_list(ierr)
53 : use utils_lib
54 : use math_lib, only: str_to_vector
55 : integer, intent(out) :: ierr
56 :
57 : integer :: iounit, i, nvec
58 : character (len=256) :: filename, string
59 : character(len=iso_name_length) :: lhs, rhs
60 : character(len=2*iso_name_length+1) :: key
61 15 : real(dp), target :: vec_ary(2)
62 : real(dp), pointer :: vec(:)
63 : integer, parameter :: max_num_weak_info = 1000
64 :
65 : logical, parameter :: dbg = .false.
66 :
67 : include 'formats'
68 :
69 5 : ierr = 0
70 5 : vec => vec_ary
71 :
72 5 : filename = trim(weak_data_dir) // '/weak_info.list'
73 5 : ierr = 0
74 5 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
75 5 : if (ierr /= 0) then
76 0 : write(*,*) 'failed to open ' // trim(filename)
77 0 : return
78 : end if
79 :
80 : if (dbg) then
81 : write(*,'(A)')
82 : write(*,*) 'weak info filename <' // trim(filename) // '>'
83 : write(*,'(A)')
84 : end if
85 :
86 : do ! skip to line starting with 'from '
87 80 : read(iounit,'(a)',iostat=ierr) string
88 80 : if (failed('read weak info comments')) return
89 85 : if (len_trim(string) > 4) then
90 60 : if (string(1:5) == 'from ') exit
91 : end if
92 : end do
93 :
94 5 : nullify(weak_info_list_dict)
95 5 : allocate(weak_info_list_halflife(max_num_weak_info))
96 5 : allocate(weak_info_list_Qneu(max_num_weak_info))
97 5 : num_weak_info_list_reactions = 0
98 1190 : do i = 1, max_num_weak_info ! keep reading until end of file
99 1190 : read(iounit,fmt='(a5,a5,a)',iostat=ierr) lhs, rhs, string
100 1190 : if (ierr == 0) then
101 1185 : call str_to_vector(string, vec, nvec, ierr)
102 1185 : if (nvec < 2) ierr = -1
103 : end if
104 1190 : if (ierr /= 0) then
105 5 : ierr = 0; exit
106 : end if
107 1185 : weak_info_list_halflife(i) = vec(1)
108 1185 : weak_info_list_Qneu(i) = vec(2)
109 1185 : call create_weak_dict_key(lhs, rhs, key)
110 : !write(*,'(a)') 'weak info list key ' // trim(key)
111 1185 : call integer_dict_define(weak_info_list_dict, key, i, ierr)
112 1185 : if (failed('integer_dict_define')) return
113 2370 : num_weak_info_list_reactions = i
114 : end do
115 :
116 5 : close(iounit)
117 :
118 5 : if (num_weak_info_list_reactions == 0) then
119 0 : ierr = -1
120 0 : write(*,*) 'failed trying to read weak_info.list -- no reactions?'
121 0 : return
122 : end if
123 :
124 5 : if (num_weak_info_list_reactions == max_num_weak_info) then
125 0 : ierr = -1
126 0 : write(*,*) 'failed trying to read weak_info.list -- too many reactions?'
127 0 : return
128 : end if
129 :
130 5 : call integer_dict_create_hash(weak_info_list_dict, ierr)
131 5 : if (ierr /= 0) return
132 :
133 5 : call realloc_double(weak_info_list_halflife, num_weak_info_list_reactions, ierr)
134 5 : if (ierr /= 0) return
135 :
136 5 : call realloc_double(weak_info_list_Qneu, num_weak_info_list_reactions, ierr)
137 15 : if (ierr /= 0) return
138 :
139 :
140 : contains
141 :
142 1265 : logical function failed(str)
143 : character (len=*) :: str
144 1265 : failed = (ierr /= 0)
145 1265 : if (failed) then
146 0 : write(*,*) 'failed: ' // trim(str)
147 : end if
148 5 : end function failed
149 :
150 :
151 : end subroutine load_weak_info_list
152 :
153 :
154 5 : subroutine private_load_weak_tables(ierr)
155 : use utils_lib
156 : use chem_lib, only: chem_get_iso_id
157 : use chem_def, only: iso_name_length
158 : integer, intent(out) :: ierr
159 :
160 : integer :: iounit, i, ios, id
161 : character (len=256) :: filename, cache_filename, string
162 : character(len=iso_name_length) :: lhs1, rhs1, lhs2, rhs2, weak_lhs, weak_rhs
163 : character(len=2*iso_name_length+1) :: key
164 :
165 : integer, parameter :: weak_num_T9 = 12, weak_num_lYeRho = 11
166 : real(dp), parameter :: weak_reaction_T9s(weak_num_T9) = &
167 : [ 0.01d0, 0.1d0, 0.2d0, 0.4d0, 0.7d0, 1.0d0, 1.5d0, 2.0d0, 3.0d0, 5.0d0, 10.0d0, 30.0d0 ]
168 : real(dp), parameter :: weak_reaction_lYeRhos(weak_num_lYeRho) = &
169 : [ 1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 6.0d0, 7.0d0, 8.0d0, 9.0d0, 10.0d0, 11.0d0 ]
170 :
171 : integer, parameter :: i_ldecay = 1, i_lcapture = 2, i_lneutrino = 3
172 :
173 : logical, parameter :: dbg = .false.
174 :
175 : include 'formats'
176 :
177 5 : ierr = 0
178 :
179 5 : ios = -1
180 5 : if (rates_use_cache) then
181 5 : cache_filename = trim(rates_cache_dir) // '/weakreactions.bin'
182 5 : ios = 0
183 : open(newunit=iounit,file=trim(cache_filename),action='read', &
184 5 : status='old',iostat=ios,form='unformatted')
185 5 : if (ios == 0) then ! opened it okay
186 : call read_weak_cache(iounit,ios)
187 4 : close(iounit)
188 : end if
189 : end if
190 :
191 5 : if (ios /= 0) then ! need to read data file
192 :
193 1 : filename = trim(weak_data_dir) // '/weakreactions.tables'
194 1 : ierr = 0
195 1 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
196 1 : if (ierr /= 0) then
197 0 : write(*,*) 'failed to open ' // trim(filename)
198 0 : return
199 : end if
200 :
201 : if (dbg) then
202 : write(*,'(A)')
203 : write(*,*) 'weaklib filename <' // trim(filename) // '>'
204 : write(*,'(A)')
205 : end if
206 :
207 : do ! skip to after line starting with '='
208 22 : read(iounit,'(a)',iostat=ierr) string
209 22 : if (failed('read header')) return
210 23 : if (len_trim(string) > 0) then
211 17 : if (string(1:1) == '=') exit
212 : end if
213 : end do
214 :
215 1 : if (.not. skip_line()) return
216 :
217 1 : read(iounit,*,iostat=ierr) num_weak_reactions
218 1 : if (failed('read num_weak_reactions')) return
219 :
220 : if (dbg) write(*,2) 'num_weak_reactions', num_weak_reactions
221 :
222 1 : call alloc
223 1 : if (failed('allocate')) return
224 :
225 463 : do i = 1, num_weak_reactions
226 462 : if (.not. skip_line()) return
227 462 : if (mod(i,2)==1) then ! first of pair
228 231 : if (.not. skip_line()) return
229 231 : if (.not. skip_line()) return
230 231 : read(iounit,fmt='(2a5)',iostat=ierr) lhs1, rhs1
231 231 : if (failed('read lhs1, rhs1')) return
232 231 : if (lhs1 == 'al-6') lhs1 = 'al26-1'
233 231 : if (rhs1 == 'al-6') rhs1 = 'al26-1'
234 231 : if (lhs1 == 'al*6') lhs1 = 'al26-2'
235 231 : if (rhs1 == 'al*6') rhs1 = 'al26-2'
236 231 : read(iounit,fmt='(2a5)',iostat=ierr) lhs2, rhs2
237 231 : if (failed('read lhs2, rhs2')) return
238 231 : if (lhs2 == 'al-6') lhs2 = 'al26-1'
239 231 : if (rhs2 == 'al-6') rhs2 = 'al26-1'
240 231 : if (lhs2 == 'al*6') lhs2 = 'al26-2'
241 231 : if (rhs2 == 'al*6') rhs2 = 'al26-2'
242 231 : if (.not. skip_line()) return
243 231 : if (.not. skip_line()) return
244 231 : weak_lhs = lhs1
245 231 : weak_rhs = rhs1
246 : else
247 231 : weak_lhs = lhs2
248 231 : weak_rhs = rhs2
249 : end if
250 462 : call adjust_name(weak_lhs)
251 462 : call adjust_name(weak_rhs)
252 462 : id = chem_get_iso_id(weak_lhs)
253 462 : if (id <= 0) then
254 0 : write(*,*) 'weaklib FATAL ERROR: unknown nuclide ' // weak_lhs
255 0 : call mesa_error(__FILE__,__LINE__)
256 : end if
257 462 : weak_lhs_nuclide_id(i) = id
258 462 : id = chem_get_iso_id(weak_rhs)
259 462 : if (id <= 0) then
260 0 : write(*,*) 'weaklib FATAL ERROR: unknown nuclide ' // weak_rhs
261 0 : call mesa_error(__FILE__,__LINE__)
262 : end if
263 462 : weak_reaclib_id(i) = 0
264 462 : weak_rhs_nuclide_id(i) = id
265 462 : weak_lhs_nuclide_name(i) = weak_lhs
266 462 : weak_rhs_nuclide_name(i) = weak_rhs
267 462 : if (.not. skip_line()) return
268 462 : call read_table(i,i_ldecay)
269 462 : if (failed('read ldecay')) return
270 462 : if (.not. skip_line()) return
271 462 : call read_table(i,i_lcapture)
272 462 : if (failed('read lcapture')) return
273 462 : if (.not. skip_line()) return
274 462 : call read_table(i,i_lneutrino)
275 463 : if (failed('read lneutrino')) return
276 : end do
277 :
278 1 : close(iounit)
279 :
280 1 : if (rates_use_cache) then
281 : open(newunit=iounit, file=trim(cache_filename), iostat=ios, &
282 1 : action='write', form='unformatted')
283 1 : if (ios == 0) then
284 1 : call write_weak_cache(iounit)
285 1 : close(iounit)
286 : end if
287 : end if
288 :
289 : end if
290 :
291 5 : nullify(weak_reactions_dict)
292 2315 : do i = 1, num_weak_reactions
293 2310 : call create_weak_dict_key(weak_lhs_nuclide_name(i), weak_rhs_nuclide_name(i), key)
294 2310 : call integer_dict_define(weak_reactions_dict, key, i, ierr)
295 2315 : if (failed('integer_dict_define')) return
296 : end do
297 :
298 5 : call integer_dict_create_hash(weak_reactions_dict, ierr)
299 5 : if (failed('integer_dict_create_hash')) return
300 :
301 2315 : do i = 1, num_weak_reactions
302 : associate(t => weak_reactions_tables(i) % t)
303 2310 : if (ierr == 0) call t% setup(ierr)
304 : end associate
305 2315 : if (failed('setup')) return
306 : end do
307 :
308 5 : if (dbg) write(*,*) 'finished load_weak_tables'
309 :
310 :
311 : contains
312 :
313 :
314 4 : subroutine read_weak_cache(iounit,ios)
315 : integer, intent(in) :: iounit
316 : integer, intent(out) :: ios
317 : integer :: n, i
318 :
319 : include 'formats'
320 :
321 4 : read(iounit,iostat=ios) num_weak_reactions
322 4 : if (ios /= 0) return
323 :
324 : if (dbg) write(*,2) 'num_weak_reactions', num_weak_reactions
325 :
326 4 : call alloc
327 4 : if (failed('allocate')) return
328 :
329 4 : n = num_weak_reactions
330 :
331 : read(iounit,iostat=ios) &
332 4 : weak_lhs_nuclide_id(1:n), &
333 4 : weak_rhs_nuclide_id(1:n), &
334 4 : weak_reaclib_id(1:n), &
335 4 : weak_lhs_nuclide_name(1:n), &
336 8 : weak_rhs_nuclide_name(1:n)
337 :
338 1852 : do i = 1, n
339 : read(iounit, iostat=ios) &
340 800188 : weak_reactions_tables(i) % t % data(1,1:weak_num_T9,1:weak_num_lYeRho,1:3)
341 : end do
342 :
343 5 : end subroutine read_weak_cache
344 :
345 :
346 1 : subroutine write_weak_cache(iounit)
347 : integer, intent(in) :: iounit
348 : integer :: n, i
349 :
350 : include 'formats'
351 :
352 1 : write(iounit) num_weak_reactions
353 :
354 1 : n = num_weak_reactions
355 :
356 : write(iounit) &
357 1 : weak_lhs_nuclide_id(1:n), &
358 1 : weak_rhs_nuclide_id(1:n), &
359 1 : weak_reaclib_id(1:n), &
360 1 : weak_lhs_nuclide_name(1:n), &
361 2 : weak_rhs_nuclide_name(1:n)
362 :
363 463 : do i = 1, n
364 : write(iounit, iostat=ios) &
365 200047 : weak_reactions_tables(i) % t % data(1,1:weak_num_T9,1:weak_num_lYeRho,1:3)
366 : end do
367 :
368 1 : end subroutine write_weak_cache
369 :
370 :
371 5 : subroutine alloc
372 :
373 : integer :: i
374 5 : type(weaklib_rate_table) :: table
375 :
376 : allocate( &
377 : weak_reaclib_id(num_weak_reactions), &
378 : weak_lhs_nuclide_name(num_weak_reactions), &
379 : weak_rhs_nuclide_name(num_weak_reactions), &
380 : weak_lhs_nuclide_id(num_weak_reactions), &
381 : weak_rhs_nuclide_id(num_weak_reactions), &
382 : weak_reactions_tables(num_weak_reactions), &
383 2315 : stat=ierr)
384 :
385 2315 : do i = 1, num_weak_reactions
386 2310 : table = weaklib_rate_table(weak_reaction_T9s, weak_reaction_lYeRhos)
387 2315 : allocate(weak_reactions_tables(i)% t, source=table)
388 : end do
389 :
390 5 : end subroutine alloc
391 :
392 :
393 924 : subroutine adjust_name(nm)
394 : character(len=iso_name_length) :: nm
395 924 : nm = adjustl(nm)
396 924 : if (nm == 'p') then
397 2 : nm = 'h1'
398 922 : else if (nm == 'n') then
399 2 : nm = 'neut'
400 : end if
401 924 : end subroutine adjust_name
402 :
403 :
404 1386 : subroutine read_table(i,ii)
405 : use math_lib, only: str_to_vector
406 : integer, intent(in) :: i, ii
407 : integer :: k, j, skip, nvec
408 : !real :: buffer(weak_num_T9)
409 : character (len=256) :: buf
410 70686 : real(dp), target :: vec_ary(50)
411 : real(dp), pointer :: vec(:)
412 : logical, parameter :: dbg = .false.
413 1386 : vec => vec_ary
414 1386 : skip = -1
415 16632 : do j = 1, weak_num_lYeRho
416 : !read(iounit,fmt=*,iostat=ierr) skip, buffer
417 15246 : read(iounit,fmt='(a)',iostat=ierr) buf
418 15246 : if (ierr == 0) then
419 15246 : call str_to_vector(buf, vec, nvec, ierr)
420 15246 : skip = int(vec(1))
421 15246 : if (nvec < weak_num_T9+1) ierr = -1
422 : end if
423 15246 : if (ierr /= 0 .or. j /= skip) then
424 : if (dbg) then
425 : write(*,*) 'error in reading table', j, skip
426 : write(*,*) 'these are the NEXT lines after the error'
427 : do k=1,20
428 : read(iounit,fmt='(a)') string
429 : write(*,'(a)') trim(string)
430 : end do
431 : write(*,'(A)')
432 : call mesa_error(__FILE__,__LINE__,'read_table')
433 : end if
434 : return
435 : end if
436 199584 : do k=1,weak_num_T9
437 198198 : weak_reactions_tables(i) % t % data(1,k,j,ii) = vec(k+1)
438 : end do
439 : !if (dbg) write(*,'(a,2i6,99f9.3)') 'read_table', j, skip, buffer
440 : end do
441 1386 : end subroutine read_table
442 :
443 :
444 9274 : logical function failed(str)
445 : character (len=*) :: str
446 9274 : failed = (ierr /= 0)
447 9274 : if (failed) then
448 0 : write(*,*) 'failed: ' // trim(str)
449 : end if
450 1386 : end function failed
451 :
452 :
453 2773 : logical function skip_line()
454 : logical, parameter :: dbg = .false.
455 : if (dbg) then
456 : read(iounit,fmt='(a)') string
457 : write(*,'(a)') 'skip line ' // trim(string)
458 : else
459 2773 : read(iounit,'(a)',iostat=ierr)
460 : end if
461 2773 : skip_line = .not. (failed('skip line'))
462 2773 : end function skip_line
463 :
464 :
465 : end subroutine private_load_weak_tables
466 :
467 5 : subroutine load_user_weak_tables(ierr)
468 : use utils_def
469 : use utils_lib
470 : use chem_lib, only: chem_get_iso_id
471 : use chem_def, only: iso_name_length
472 : use weak_support, only: parse_weak_rate_name
473 :
474 : integer, intent(out) :: ierr
475 :
476 : character (len=256) :: filename
477 : character(len=iso_name_length) :: lhs, rhs
478 : character(len=2*iso_name_length+1) :: key
479 :
480 : integer :: i, iounit, n, t, ir, id
481 : character (len=256) :: dir, rate_name, rate_fname, buffer
482 :
483 : logical, parameter :: dbg = .false.
484 :
485 : include 'formats'
486 :
487 5 : ierr = 0
488 :
489 : if (dbg) write(*,*) 'load_user_weak_tables'
490 :
491 : ! first try local rate_tables_dir
492 5 : dir = rates_table_dir
493 5 : filename = trim(dir) // '/weak_rate_list.txt'
494 5 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
495 5 : if (ierr /= 0) then ! if don't find that file, look in rates_dir
496 3 : dir = trim(rates_dir) // '/rate_tables'
497 3 : filename = trim(dir) // '/weak_rate_list.txt'
498 3 : ierr = 0
499 3 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
500 3 : if (ierr /= 0) then
501 0 : write(*,*) 'failed to open weak rates list file ' // trim(filename)
502 0 : return
503 : end if
504 : end if
505 :
506 5 : n = 0
507 5 : i = 0
508 :
509 : if (dbg) write(*,*) 'read rate list file ' // trim(filename)
510 :
511 9 : rate_loop: do
512 8 : t = token(iounit, n, i, buffer, rate_name)
513 8 : if (t == eof_token) exit rate_loop
514 3 : if (t /= name_token) then
515 0 : call error; return
516 : end if
517 : if (dbg) write(*,*) 'use rate table from file for ', trim(rate_name)
518 :
519 : ! first, parse the weak rate id
520 3 : call parse_weak_rate_name(rate_name, lhs, rhs, ierr)
521 : if (dbg) write(*,*) 'parse_weak_rate_name gives ', trim(lhs), ' ', trim(rhs), ierr
522 :
523 : ! check if we already have a rate with this name
524 3 : call create_weak_dict_key(lhs, rhs, key)
525 : if (dbg) write(*,*) 'key is ', trim(key), ierr
526 :
527 : !write(*,'(a)') 'weak info list key ' // trim(key)
528 3 : call integer_dict_lookup(weak_reactions_dict, key, ir, ierr)
529 : if (dbg) write(*,*) ir, ierr
530 :
531 :
532 3 : if (ierr /= 0 .or. ir <= 0) then
533 3 : call extend
534 3 : ir = num_weak_reactions
535 :
536 3 : id = chem_get_iso_id(lhs)
537 3 : if (id <= 0) then
538 0 : write(*,*) 'weaklib FATAL ERROR: unknown nuclide ' // lhs
539 0 : call mesa_error(__FILE__,__LINE__)
540 : end if
541 3 : weak_lhs_nuclide_id(ir) = id
542 :
543 3 : id = chem_get_iso_id(rhs)
544 3 : if (id <= 0) then
545 0 : write(*,*) 'weaklib FATAL ERROR: unknown nuclide ' // rhs
546 0 : call mesa_error(__FILE__,__LINE__)
547 : end if
548 3 : weak_reaclib_id(ir) = 0
549 3 : weak_rhs_nuclide_id(ir) = id
550 3 : weak_lhs_nuclide_name(ir) = lhs
551 3 : weak_rhs_nuclide_name(ir) = rhs
552 :
553 : end if
554 :
555 3 : t = token(iounit, n, i, buffer, rate_fname)
556 3 : if (t /= string_token) then
557 0 : call error; return
558 : end if
559 : if (dbg) write(*,*) 'rate_fname ', trim(rate_fname)
560 :
561 3 : call read_hd5_file
562 3 : if (ierr < 0) then
563 0 : write(*,*) 'failed to read hdf5 file in load_user_weak_tables'
564 0 : call mesa_error(__FILE__,__LINE__)
565 : end if
566 :
567 3 : nullify(weak_reactions_dict)
568 1392 : do i = 1, num_weak_reactions
569 1389 : call create_weak_dict_key(weak_lhs_nuclide_name(i), weak_rhs_nuclide_name(i), key)
570 1389 : call integer_dict_define(weak_reactions_dict, key, i, ierr)
571 1392 : if (failed('integer_dict_define')) return
572 : end do
573 :
574 3 : call integer_dict_create_hash(weak_reactions_dict, ierr)
575 8 : if (failed('integer_dict_create_hash')) return
576 :
577 : end do rate_loop
578 :
579 5 : close(iounit)
580 :
581 5 : if (dbg) write(*,*) 'finished load_weak_tables'
582 :
583 : contains
584 :
585 3 : subroutine read_hd5_file
586 :
587 5 : use forum_m, only: hdf5io_t, OPEN_FILE_RO
588 :
589 : character (len=256) :: filename
590 3 : type(hdf5io_t) :: hi
591 3 : real(dp), allocatable, dimension(:) :: T9s, lYeRhos
592 : integer :: num_T9, num_lYeRho
593 3 : type(weaklib_rate_table) :: table
594 :
595 : logical, parameter :: dbg = .false.
596 :
597 3 : filename = trim(dir) // '/' // trim(rate_fname)
598 3 : write(*,*) 'reading user weak rate file ', trim(filename)
599 :
600 : ! open file (read-only)
601 :
602 3 : hi = hdf5io_t(filename, OPEN_FILE_RO)
603 :
604 : ! read axis data
605 :
606 3 : call hi% alloc_read_dset('T9s', T9s)
607 3 : num_T9 = SIZE(T9s)
608 :
609 3 : call hi% alloc_read_dset('lYeRhos', lYeRhos)
610 3 : num_lYeRho = SIZE(lYeRhos)
611 :
612 : ! create the table
613 :
614 3 : table = weaklib_rate_table(T9s, lYeRhos)
615 :
616 : ! read data into it
617 :
618 3 : call hi% read_dset('ldecay', table% data(1, 1:num_T9, 1:num_lYeRho, table% i_ldecay))
619 3 : call hi% read_dset('lcapture', table% data(1, 1:num_T9, 1:num_lYeRho, table% i_lcapture))
620 3 : call hi% read_dset('lneutrino', table% data(1, 1:num_T9, 1:num_lYeRho, table% i_lneutrino))
621 :
622 : ! store the table
623 :
624 3 : allocate(weak_reactions_tables(ir)% t, source=table)
625 : associate(t => weak_reactions_tables(ir) % t)
626 3 : if (ierr == 0) call t% setup(ierr)
627 6 : if (failed('setup')) return
628 : end associate
629 :
630 : ! close file
631 :
632 3 : call hi% final()
633 :
634 6 : end subroutine read_hd5_file
635 :
636 :
637 3 : subroutine extend
638 : integer :: n
639 :
640 3 : type(table_c), dimension(:), allocatable :: tmp_weak_reactions_tables
641 :
642 : integer, allocatable, dimension(:) :: &
643 3 : tmp_weak_lhs_nuclide_id, tmp_weak_rhs_nuclide_id, tmp_weak_reaclib_id
644 : character(len=iso_name_length), dimension(:), allocatable :: &
645 3 : tmp_weak_lhs_nuclide_name, tmp_weak_rhs_nuclide_name
646 :
647 3 : n = num_weak_reactions + 1
648 :
649 : allocate( &
650 : tmp_weak_reaclib_id(n), &
651 : tmp_weak_lhs_nuclide_name(n), &
652 : tmp_weak_rhs_nuclide_name(n), &
653 : tmp_weak_lhs_nuclide_id(n), &
654 : tmp_weak_rhs_nuclide_id(n), &
655 3 : stat=ierr)
656 :
657 1389 : tmp_weak_reaclib_id(1:num_weak_reactions) = weak_reaclib_id
658 1389 : tmp_weak_lhs_nuclide_name(1:num_weak_reactions) = weak_lhs_nuclide_name
659 1389 : tmp_weak_rhs_nuclide_name(1:num_weak_reactions) = weak_rhs_nuclide_name
660 1389 : tmp_weak_lhs_nuclide_id(1:num_weak_reactions) = weak_lhs_nuclide_id
661 1389 : tmp_weak_rhs_nuclide_id(1:num_weak_reactions) = weak_rhs_nuclide_id
662 :
663 : allocate( &
664 : weak_reaclib_id(n), &
665 : weak_lhs_nuclide_name(n), &
666 : weak_rhs_nuclide_name(n), &
667 : weak_lhs_nuclide_id(n), &
668 : weak_rhs_nuclide_id(n), &
669 3 : stat=ierr)
670 :
671 1389 : weak_reaclib_id(1:num_weak_reactions) = tmp_weak_reaclib_id(1:num_weak_reactions)
672 1389 : weak_lhs_nuclide_name(1:num_weak_reactions) = tmp_weak_lhs_nuclide_name(1:num_weak_reactions)
673 1389 : weak_rhs_nuclide_name(1:num_weak_reactions) = tmp_weak_rhs_nuclide_name(1:num_weak_reactions)
674 1389 : weak_lhs_nuclide_id(1:num_weak_reactions) = tmp_weak_lhs_nuclide_id(1:num_weak_reactions)
675 1389 : weak_rhs_nuclide_id(1:num_weak_reactions) = tmp_weak_rhs_nuclide_id(1:num_weak_reactions)
676 :
677 : deallocate( &
678 : tmp_weak_reaclib_id, &
679 : tmp_weak_lhs_nuclide_name, &
680 : tmp_weak_rhs_nuclide_name, &
681 : tmp_weak_lhs_nuclide_id, &
682 : tmp_weak_rhs_nuclide_id, &
683 3 : stat=ierr)
684 :
685 1392 : allocate(tmp_weak_reactions_tables(n))
686 1389 : tmp_weak_reactions_tables(1:num_weak_reactions) = weak_reactions_tables
687 3 : call move_alloc(tmp_weak_reactions_tables, weak_reactions_tables)
688 :
689 3 : num_weak_reactions = n
690 :
691 3 : end subroutine extend
692 :
693 0 : subroutine error
694 0 : ierr = -1
695 0 : close(iounit)
696 0 : end subroutine error
697 :
698 1395 : logical function failed(str)
699 : character (len=*) :: str
700 1395 : failed = (ierr /= 0)
701 1395 : if (failed) then
702 0 : write(*,*) 'failed: ' // trim(str)
703 : end if
704 1395 : end function failed
705 :
706 :
707 : end subroutine load_user_weak_tables
708 :
709 : end module load_weak
710 :
|
