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 rates_initialize
21 : use const_def, only: dp, mesa_data_dir, ln2
22 : use math_lib
23 : use rates_def
24 :
25 : implicit none
26 :
27 : contains
28 :
29 10 : subroutine finish_rates_def_init
30 : use utils_lib, only: integer_dict_define, integer_dict_create_hash, integer_dict_size
31 : use reaclib_input, only: do_extract_rates
32 : use chem_def, only: nuclide_set, chem_isos, num_chem_isos, iso_name_length
33 : use chem_lib, only: generate_nuclide_set
34 : integer :: ierr
35 : character(len=iso_name_length), pointer :: names(:) =>null()
36 : type(nuclide_set), pointer :: set(:) =>null()
37 : integer, pointer :: chem_id(:) =>null()
38 : ! will be allocated by extract_nuclides_from_chem_isos
39 : logical :: use_weaklib
40 : include 'formats'
41 :
42 : ierr = 0
43 5 : call integer_dict_create_hash(reaction_names_dict, ierr)
44 5 : if (ierr /= 0) then
45 0 : write(*,*) 'FATAL ERROR: rates_def_init failed in integer_dict_create_hash'
46 0 : return
47 : end if
48 :
49 : ! set up reaclib info
50 5 : allocate(set(num_chem_isos))
51 :
52 5 : call generate_nuclide_set(chem_isos% name, set)
53 :
54 5 : use_weaklib = .true.
55 5 : call do_extract_rates(set, chem_isos, reaclib_rates, use_weaklib, ierr)
56 5 : deallocate(set)
57 5 : if (ierr /= 0) then
58 0 : write(*,*) 'FATAL ERROR: extract_reaclib_rates failed in rates_def_init'
59 0 : return
60 : end if
61 :
62 15 : end subroutine finish_rates_def_init
63 :
64 :
65 0 : subroutine do_add_reaction_for_handle(reaction_handle, ierr)
66 5 : use reaclib_support, only: do_parse_reaction_handle
67 : character (len=*), intent(in) :: reaction_handle ! to be added
68 : integer, intent(out) :: ierr
69 :
70 : integer :: ir, num_in, num_out
71 : integer :: particles_in, particles_out
72 : logical :: already_defined
73 : integer :: iso_ids(max_num_reaction_inputs+max_num_reaction_outputs)
74 : integer :: cin(max_num_reaction_inputs), cout(max_num_reaction_outputs)
75 : character (len=16) :: op ! e.g., 'pg', 'wk', 'to', or ...
76 :
77 : logical, parameter :: weak = .false.
78 : logical, parameter :: dbg = .false.
79 :
80 : include 'formats'
81 :
82 : ierr = 0
83 :
84 : if (dbg) write(*,*) 'do_add_reaction_for_handle ' // trim(reaction_handle)
85 :
86 : call do_parse_reaction_handle( &
87 0 : reaction_handle, particles_in, particles_out, iso_ids, op, ierr)
88 0 : if (ierr /= 0) then
89 : write(*,'(a)') 'add_reaction_for_handle failed in reaclib_parse_handle ' // &
90 0 : trim(reaction_handle)
91 0 : return
92 : end if
93 :
94 0 : call alloc_reaction_ir(reaction_handle, ir, already_defined, ierr)
95 0 : if (already_defined) return
96 0 : if (ierr /= 0) return
97 :
98 0 : reaction_inputs(:,ir) = 0
99 0 : reaction_outputs(:,ir) = 0
100 :
101 0 : cin(:) = 1
102 0 : cout(:) = 1
103 :
104 0 : call setup(reaction_inputs(:,ir), num_in, cin, 0, particles_in)
105 0 : call setup(reaction_outputs(:,ir), num_out, cout, particles_in, particles_out)
106 :
107 : call set_reaction_info( &
108 0 : ir, num_in, num_out, particles_in, particles_out, weak, reaction_handle, ierr)
109 0 : if (ierr /= 0) then
110 0 : write(*,*) 'failed in set_reaction_info for ' // trim(reaction_handle), ir
111 : end if
112 :
113 0 : if (dbg) write(*,*) 'done do_add_reaction_for_handle ' // trim(reaction_handle)
114 :
115 :
116 : contains
117 :
118 :
119 0 : subroutine setup(reaction_io, num, cnt, k, num_particles)
120 : integer :: reaction_io(:), num, cnt(:), k, num_particles
121 : integer :: j
122 0 : reaction_io(1) = 1
123 0 : reaction_io(2) = iso_ids(k+1)
124 0 : num = 1
125 0 : cnt(num) = 1
126 0 : do j=2,num_particles
127 0 : if (iso_ids(k+j) == iso_ids(k+j-1)) then
128 0 : cnt(num) = cnt(num) + 1
129 : else
130 0 : num = num + 1
131 0 : reaction_io(2*num) = iso_ids(k+j)
132 0 : cnt(num) = 1
133 : end if
134 0 : reaction_io(2*num-1) = cnt(num)
135 : end do
136 0 : end subroutine setup
137 :
138 :
139 : end subroutine do_add_reaction_for_handle
140 :
141 :
142 54 : subroutine do_add_reaction_from_reaclib(reaction_handle, reverse_handle, indx, ierr)
143 :
144 : character (len=*), intent(in) :: reaction_handle ! to be added
145 : character (len=*), intent(in) :: reverse_handle ! = '' if not a reverse
146 : integer, intent(in) :: indx ! index in reaclib rates
147 : integer, intent(out) :: ierr
148 :
149 : integer :: i, ir, chapter, num_in, num_out
150 : integer :: particles_in, particles_out
151 : logical :: weak, reverse, already_defined
152 : integer :: cin(max_num_reaction_inputs), cout(max_num_reaction_outputs)
153 : type (reaction_data), pointer :: r =>null()
154 :
155 : logical, parameter :: dbg = .false.
156 :
157 : include 'formats'
158 :
159 18 : ierr = 0
160 18 : i = indx
161 18 : reverse = (len_trim(reverse_handle) > 0)
162 18 : r => reaclib_rates
163 :
164 72 : cin(:) = 1
165 90 : cout(:) = 1
166 :
167 : if (dbg) write(*,'(a, 2x, i5)') 'do_add_reaction_from_reaclib ' // trim(reaction_handle), i
168 :
169 18 : if (reverse) then
170 2 : if (reverse_handle /= r% reaction_handle(i)) then
171 0 : write(*,'(a)') trim(reverse_handle) // ' ' // trim(r% reaction_handle(i))
172 0 : write(*,'(a,3x,i8)') 'bad reverse_handle for add_reaction_from_reaclib', indx
173 0 : ierr = -1
174 0 : return
175 : end if
176 : else
177 16 : if (reaction_handle /= r% reaction_handle(i)) then
178 0 : write(*,'(a)') trim(reaction_handle) // ' ' // trim(r% reaction_handle(i))
179 0 : write(*,'(a,3x,i8)') 'bad reaction_handle for add_reaction_from_reaclib', indx
180 0 : ierr = -1
181 0 : return
182 : end if
183 : end if
184 :
185 18 : chapter = r% chapter(i)
186 18 : weak = (adjustl(r% reaction_flag(i)) == 'w')
187 :
188 18 : call alloc_reaction_ir(reaction_handle, ir, already_defined, ierr)
189 18 : if (already_defined) return
190 18 : if (ierr /= 0) return
191 :
192 126 : reaction_inputs(:,ir) = 0
193 162 : reaction_outputs(:,ir) = 0
194 :
195 18 : if (.not. reverse) then
196 16 : particles_in = Nin(chapter)
197 16 : particles_out = Nout(chapter)
198 16 : call setup(reaction_inputs(:,ir), num_in, cin, 0, particles_in)
199 16 : call setup(reaction_outputs(:,ir), num_out, cout, particles_in, particles_out)
200 : else
201 2 : particles_out = Nin(chapter)
202 2 : particles_in = Nout(chapter)
203 2 : call setup(reaction_inputs(:,ir), num_in, cin, particles_out, particles_in)
204 2 : call setup(reaction_outputs(:,ir), num_out, cout, 0, particles_out)
205 : end if
206 :
207 : call set_reaction_info( &
208 18 : ir, num_in, num_out, particles_in, particles_out, weak, reaction_handle, ierr)
209 18 : if (ierr /= 0) then
210 0 : write(*,*) 'failed in set_reaction_info for ' // trim(reaction_handle), ir
211 : end if
212 :
213 : if (dbg) write(*,*) 'done do_add_reaction_from_reaclib'
214 :
215 :
216 : contains
217 :
218 :
219 36 : subroutine setup(reaction_test, num, cnt, k, num_particles)
220 : integer :: reaction_test(:), num, cnt(:), k, num_particles
221 : integer :: j
222 36 : reaction_test(1) = 1
223 36 : reaction_test(2) = r% pspecies(k+1,i)
224 36 : num = 1
225 36 : cnt(num) = 1
226 40 : do j=2,num_particles
227 4 : if (r% pspecies(k+j,i) == r% pspecies(k+j-1,i)) then
228 0 : cnt(num) = cnt(num) + 1
229 : else
230 4 : num = num + 1
231 4 : reaction_test(2*num) = r% pspecies(k+j,i)
232 4 : cnt(num) = 1
233 : end if
234 40 : reaction_test(2*num-1) = cnt(num)
235 : end do
236 36 : end subroutine setup
237 :
238 :
239 : end subroutine do_add_reaction_from_reaclib
240 :
241 :
242 18 : subroutine alloc_reaction_ir(reaction_handle, ir, already_defined, ierr)
243 : character (len=*), intent(in) :: reaction_handle
244 : integer, intent(out) :: ir
245 : logical, intent(out) :: already_defined
246 : integer, intent(out) :: ierr
247 : logical, parameter :: dbg = .false.
248 : include 'formats'
249 18 : ierr = 0
250 18 : ir = get_rates_reaction_id(reaction_handle)
251 18 : if (ir > 0) then
252 0 : already_defined = .true.
253 0 : return
254 : end if
255 18 : !$omp critical (lock_alloc_reaction_ir)
256 18 : ir = get_rates_reaction_id(reaction_handle)
257 18 : if (ir > 0) then
258 0 : already_defined = .true.
259 : else
260 18 : already_defined = .false.
261 : if (dbg) write(*,2) 'call increase_num_reactions', rates_reaction_id_max
262 18 : call increase_num_reactions(ierr)
263 18 : if (ierr == 0) then
264 18 : ir = rates_reaction_id_max
265 : if (dbg) write(*,2) 'after increase_num_reactions', rates_reaction_id_max
266 : end if
267 : end if
268 : !$omp end critical (lock_alloc_reaction_ir)
269 : end subroutine alloc_reaction_ir
270 :
271 :
272 18 : subroutine set_reaction_info( &
273 : ir, num_in, num_out, particles_in, particles_out, weak, reaction_handle, ierr)
274 : use chem_def
275 : use utils_lib, only: integer_dict_define
276 : integer, intent(in) :: ir, num_in, num_out, particles_in, particles_out
277 : logical, intent(in) :: weak
278 : character (len=*), intent(in) :: reaction_handle
279 : integer, intent(out) :: ierr
280 :
281 : integer :: j, iso_in, iso_out, weak_j, cin1
282 :
283 : logical, parameter :: dbg = .false.
284 :
285 : include 'formats'
286 :
287 18 : ierr = 0
288 :
289 18 : reaction_Name(ir) = reaction_handle
290 :
291 : if (dbg) write(*,*) 'call get_Qtotal'
292 18 : std_reaction_Qs(ir) = get_Qtotal(ir)
293 18 : std_reaction_neuQs(ir) = 0d0
294 18 : weak_lowT_rate(ir) = -1d99
295 :
296 18 : iso_in = reaction_inputs(num_in*2,ir)
297 18 : cin1 = reaction_inputs(num_in*2-1,ir)
298 18 : iso_out = reaction_outputs(num_out*2,ir)
299 :
300 18 : if (iso_in < 0 .or. iso_in > num_chem_isos) then
301 0 : write(*,3) 'bad iso_in', iso_in, num_in
302 0 : write(*,*) trim(reaction_handle)
303 0 : do j = 1, num_in
304 0 : write(*,3) 'reaction input', j, reaction_inputs(2*j,ir)
305 : end do
306 0 : call mesa_error(__FILE__,__LINE__,'set_reaction_info')
307 : end if
308 :
309 18 : if (iso_out < 0 .or. iso_out > num_chem_isos) then
310 0 : write(*,2) 'bad iso_out', iso_out
311 0 : write(*,2) 'num_out', num_out
312 0 : write(*,2) 'num_in', num_in
313 0 : write(*,*) trim(reaction_handle)
314 0 : do j = 1, num_out
315 0 : write(*,3) 'reaction output', j, reaction_outputs(2*j,ir)
316 : end do
317 0 : call mesa_error(__FILE__,__LINE__,'set_reaction_info')
318 : end if
319 :
320 18 : if (weak) then
321 14 : weak_reaction_info(1,ir) = iso_in
322 14 : weak_reaction_info(2,ir) = iso_out
323 14 : weak_j = do_get_weak_info_list_id(chem_isos% name(iso_in),chem_isos% name(iso_out))
324 14 : if (weak_j > 0) std_reaction_neuQs(ir) = weak_info_list_Qneu(weak_j)
325 14 : call set_weak_lowT_rate(ir, ierr)
326 14 : if (ierr /= 0) return
327 : else
328 12 : weak_reaction_info(1:2,ir) = 0
329 : end if
330 :
331 18 : reaction_ye_rho_exponents(1,ir) = 0 ! 1 for electron captures, 0 for rest.
332 :
333 18 : if (particles_in > 1) then
334 2 : reaction_screening_info(1,ir) = reaction_inputs(2,ir)
335 2 : if (cin1 > 1) then
336 0 : reaction_screening_info(2,ir) = reaction_screening_info(1,ir)
337 : else
338 2 : reaction_screening_info(2,ir) = reaction_inputs(4,ir)
339 : end if
340 2 : reaction_ye_rho_exponents(2,ir) = particles_in - 1
341 : else
342 48 : reaction_screening_info(1:2,ir) = 0
343 16 : reaction_ye_rho_exponents(2,ir) = 0
344 : end if
345 18 : reaction_screening_info(3,ir) = 0
346 18 : reaction_categories(ir) = iother
347 18 : if (is_pp_reaction(reaction_handle)) then
348 0 : reaction_categories(ir) = ipp
349 18 : else if (is_cno_reaction(reaction_handle)) then
350 12 : reaction_categories(ir) = icno
351 6 : else if (num_in == 1 .and. cin1 == 2) then
352 0 : if (iso_in == ic12) then
353 0 : reaction_categories(ir) = icc
354 0 : else if (iso_in == io16) then
355 0 : reaction_categories(ir) = ioo
356 : end if
357 6 : else if (num_in > 1) then
358 2 : select case(chem_isos% Z(iso_in))
359 : case (6)
360 0 : reaction_categories(ir) = i_burn_c
361 : case (7)
362 0 : reaction_categories(ir) = i_burn_n
363 : case (8)
364 0 : reaction_categories(ir) = i_burn_o
365 : case (10)
366 2 : reaction_categories(ir) = i_burn_ne
367 : case (11)
368 0 : reaction_categories(ir) = i_burn_na
369 : case (12)
370 0 : reaction_categories(ir) = i_burn_mg
371 : case (14)
372 0 : reaction_categories(ir) = i_burn_si
373 : case (16)
374 0 : reaction_categories(ir) = i_burn_s
375 : case (18)
376 0 : reaction_categories(ir) = i_burn_ar
377 : case (20)
378 0 : reaction_categories(ir) = i_burn_ca
379 : case (22)
380 0 : reaction_categories(ir) = i_burn_ti
381 : case (24)
382 0 : reaction_categories(ir) = i_burn_cr
383 : case (26)
384 2 : reaction_categories(ir) = i_burn_fe
385 : end select
386 4 : else if (particles_in == 1 .and. particles_out == 2 .and. .not. weak) then
387 2 : reaction_categories(ir) = iphoto
388 : end if
389 :
390 18 : reaction_Info(ir) = reaction_handle
391 :
392 : if (dbg) write(*,'(a,3x,i5)') 'call integer_dict_define ' // trim(reaction_handle), ir
393 :
394 18 : call integer_dict_define(reaction_names_dict, reaction_handle, ir, ierr)
395 18 : if (ierr /= 0) then
396 0 : write(*,*) 'FATAL ERROR: set_reaction_info failed in integer_dict_define'
397 0 : return
398 : end if
399 :
400 36 : end subroutine set_reaction_info
401 :
402 :
403 214 : subroutine set_weak_lowT_rate(ir, ierr)
404 18 : use chem_def
405 : use utils_lib, only: integer_dict_define
406 : use reaclib_eval, only: do_reaclib_indices_for_reaction
407 : use ratelib, only: reaclib_rate_and_dlnT
408 : use load_weak, only: do_get_weak_rate_id
409 :
410 : integer, intent(in) :: ir
411 : integer, intent(out) :: ierr
412 : integer :: i, lo, hi, weak_reaclib_id_i
413 214 : real(dp) :: half_life, lambda, dlambda_dlnT, rlambda, drlambda_dlnT
414 :
415 : include 'formats'
416 :
417 214 : weak_reaclib_id_i = 0
418 :
419 214 : if (weak_reaction_info(1,ir) == 0 .or. weak_reaction_info(2,ir) == 0) return
420 :
421 : call do_reaclib_indices_for_reaction( &
422 209 : reaction_Name(ir), reaclib_rates, lo, hi, ierr)
423 209 : if (ierr /= 0) then ! not in reaclib
424 195 : ierr = 0
425 : i = do_get_weak_info_list_id( &
426 : chem_isos% name(weak_reaction_info(1,ir)), &
427 195 : chem_isos% name(weak_reaction_info(2,ir)))
428 195 : if (i > 0) then
429 20 : half_life = weak_info_list_halflife(i)
430 20 : if (half_life > 0d0) then
431 20 : weak_lowT_rate(ir) = ln2/half_life
432 20 : return
433 : end if
434 : end if
435 175 : weak_lowT_rate(ir) = -1d-99
436 175 : return
437 : end if
438 :
439 : ! evaluate rates at T= 10^7 K (i.e. T9 = 0.01)
440 : call reaclib_rate_and_dlnT( &
441 : lo, hi, reaction_Name(ir), 0.01_dp, &
442 14 : lambda, dlambda_dlnT, rlambda, drlambda_dlnT, ierr)
443 14 : if (ierr /= 0) then
444 : write(*,2) 'set_reaction_info failed in reaclib_rate_and_dlnT ' // &
445 0 : trim(reaction_Name(ir)), ir
446 0 : call mesa_error(__FILE__,__LINE__)
447 : end if
448 14 : if (lo <= reaclib_rates% num_from_weaklib) then
449 8 : if (.not. reaclib_rates% also_in_reaclib(lo)) lambda = -1
450 : end if
451 14 : weak_lowT_rate(ir) = lambda
452 :
453 : weak_reaclib_id_i = do_get_weak_rate_id( &
454 : chem_isos% name(weak_reaction_info(1,ir)), &
455 14 : chem_isos% name(weak_reaction_info(2,ir)))
456 14 : if (weak_reaclib_id_i == 0) return
457 8 : weak_reaclib_id(weak_reaclib_id_i) = ir
458 :
459 214 : end subroutine set_weak_lowT_rate
460 :
461 :
462 18 : logical function is_pp_reaction(reaction_handle)
463 : character (len=*), intent(in) :: reaction_handle
464 : is_pp_reaction = &
465 : reaction_handle == 'r_h1_h1_wk_h2' .or. &
466 : reaction_handle == 'r_h1_h1_ec_h2' .or. &
467 : reaction_handle == 'r_h2_pg_he3' .or. &
468 : reaction_handle == 'r_he3_he3_to_h1_h1_he4' .or. &
469 : reaction_handle == 'r_he3_ag_be7' .or. &
470 : reaction_handle == 'r_be7_wk_li7' .or. &
471 : reaction_handle == 'r_h1_li7_to_he4_he4' .or. &
472 : reaction_handle == 'r_li7_pa_he4' .or. &
473 : reaction_handle == 'r_be7_pg_b8' .or. &
474 18 : reaction_handle == 'r_b8_wk_he4_he4'
475 214 : end function is_pp_reaction
476 :
477 :
478 18 : logical function is_cno_reaction(reaction_handle)
479 : character (len=*), intent(in) :: reaction_handle
480 : is_cno_reaction = &
481 : reaction_handle == 'r_c12_pg_n13' .or. &
482 : reaction_handle == 'r_c13_pg_n14' .or. &
483 : reaction_handle == 'r_n13_wk_c13' .or. &
484 : reaction_handle == 'r_n13_pg_o14' .or. &
485 : reaction_handle == 'r_n14_pg_o15' .or. &
486 : reaction_handle == 'r_n15_pg_o16' .or. &
487 : reaction_handle == 'r_n15_pa_c12' .or. &
488 : reaction_handle == 'r_o14_wk_n14' .or. &
489 : reaction_handle == 'r_o14_ap_f17' .or. &
490 : reaction_handle == 'r_o15_wk_n15' .or. &
491 : reaction_handle == 'r_o16_pg_f17' .or. &
492 : reaction_handle == 'r_o17_pg_f18' .or. &
493 : reaction_handle == 'r_o17_pa_n14' .or. &
494 : reaction_handle == 'r_o18_pg_f19' .or. &
495 : reaction_handle == 'r_o18_pa_n15' .or. &
496 : reaction_handle == 'r_f17_wk_o17' .or. &
497 : reaction_handle == 'r_f17_pg_ne18' .or. &
498 : reaction_handle == 'r_f18_pa_o15' .or. &
499 : reaction_handle == 'r_f18_wk_o18' .or. &
500 : reaction_handle == 'r_f19_pa_o16' .or. &
501 18 : reaction_handle == 'r_ne18_wk_f18'
502 18 : end function is_cno_reaction
503 :
504 :
505 3 : subroutine free_raw_rates_records
506 : type (rate_table_info), pointer :: ri =>null()
507 : integer :: i
508 3 : if (ASSOCIATED(raw_rates_records)) then
509 954 : do i = 1, rates_reaction_id_max
510 951 : ri => raw_rates_records(i)
511 951 : if (ASSOCIATED(ri% T8s)) deallocate(ri% T8s)
512 954 : if (ASSOCIATED(ri% f1)) deallocate(ri% f1)
513 : end do
514 3 : deallocate(raw_rates_records)
515 : end if
516 3 : end subroutine free_raw_rates_records
517 :
518 :
519 21 : subroutine init_raw_rates_records(ierr)
520 : use utils_lib
521 : use utils_def
522 : integer, intent(out) :: ierr
523 :
524 : type (rate_table_info), pointer :: ri =>null()
525 : integer :: i, iounit, n, t
526 : character (len=256) :: dir, rate_name, rate_fname, filename
527 : character (len=256) :: buffer
528 : logical :: okay
529 :
530 : logical, parameter :: dbg = .false.
531 :
532 : include 'formats'
533 :
534 : if (dbg) write(*,*) 'init_raw_rates_records'
535 :
536 21 : ierr = 0
537 :
538 : ! first try local rate_tables_dir
539 21 : dir = rates_table_dir
540 21 : filename = trim(dir) // '/rate_list.txt'
541 21 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
542 21 : if (ierr /= 0) then ! if don't find that file, look in rates_dir
543 3 : dir = trim(rates_dir) // '/rate_tables'
544 3 : filename = trim(dir) // '/rate_list.txt'
545 3 : ierr = 0
546 3 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
547 3 : if (ierr /= 0) then
548 0 : write(*,*) 'failed to open rates list file ' // trim(filename)
549 0 : return
550 : end if
551 : end if
552 21 : rates_table_dir = dir
553 :
554 21 : n = 0
555 21 : i = 0
556 :
557 : if (dbg) write(*,*) 'read rate list file ' // trim(filename)
558 :
559 69 : rate_loop: do
560 90 : t = token(iounit, n, i, buffer, rate_name)
561 90 : if (t == eof_token) exit rate_loop
562 69 : if (t /= name_token) then
563 0 : call error; return
564 : end if
565 : if (dbg) write(*,*) 'use rate table from file for ', trim(rate_name)
566 69 : t = token(iounit, n, i, buffer, rate_fname)
567 69 : if (t /= string_token) then
568 0 : call error; return
569 : end if
570 :
571 69 : call rate_from_file(rate_name, rate_fname, ierr)
572 90 : if (ierr /= 0) call error()
573 :
574 : end do rate_loop
575 :
576 21 : close(iounit)
577 :
578 : if (dbg) call mesa_error(__FILE__,__LINE__,'read rates')
579 : !call check
580 :
581 :
582 : contains
583 :
584 :
585 : subroutine check
586 : ! check that there are cases for all of the rates
587 : use ratelib, only: tfactors
588 : use raw_rates, only: set_raw_rate
589 : real(dp) :: logT, temp, raw_rate
590 : integer :: i, ierr
591 : type (T_Factors) :: tf
592 :
593 : logT = 8
594 : temp = exp10(logT)
595 : call tfactors(tf, logT, temp)
596 : ierr = 0
597 : okay = .true.
598 : do i = 1, rates_reaction_id_max
599 : call set_raw_rate(i, temp, tf, raw_rate, ierr)
600 : if (ierr /= 0) then
601 : write(*,'(a)') 'set_raw_rate failed for ' // reaction_Name(i)
602 : okay = .false.
603 : ierr = 0
604 : end if
605 : end do
606 : if (.not. okay) call mesa_error(__FILE__,__LINE__,'init_raw_rates_records')
607 : end subroutine check
608 :
609 :
610 : logical function failed(str)
611 : character (len=*), intent(in) :: str
612 : failed = .false.
613 : if (ierr == 0) return
614 : if (len_trim(str) > 0) then
615 : write(*,*) trim(str) // ' failed in reading ' // trim(filename)
616 : else ! non-fatal error, so just quietly stop reading
617 : ierr = 0
618 : end if
619 : failed = .true.
620 : return
621 : end function failed
622 :
623 :
624 0 : subroutine error
625 0 : ierr = -1
626 0 : close(iounit)
627 0 : end subroutine error
628 :
629 :
630 : end subroutine init_raw_rates_records
631 :
632 569 : subroutine rate_from_file(rate_name, filename, ierr)
633 : character(len=*) :: rate_name, filename
634 : integer, intent(out) :: ierr
635 :
636 : integer :: ir
637 : logical,parameter :: dbg=.false.
638 : type (rate_table_info), pointer :: ri =>null()
639 :
640 569 : ierr = 0
641 :
642 69 : if (len_trim(rate_name) == 0 .or. len_trim(filename)==0) return
643 :
644 69 : ir = lookup_rate_name(rate_name)
645 69 : if (ir <= 0) then
646 0 : call do_add_reaction_for_handle(rate_name, ierr)
647 0 : if (ierr == 0) ir = lookup_rate_name(rate_name)
648 0 : if (ierr /= 0 .or. ir <= 0) then
649 0 : write(*,*) 'invalid rate file ' // trim(rate_name)
650 0 : ierr = -1
651 0 : return
652 : end if
653 : end if
654 : if (dbg) write(*,*) 'rate_fname ', trim(filename)
655 69 : ri => raw_rates_records(ir)
656 69 : ri% use_rate_table = .true.
657 69 : ri% need_to_read = .true.
658 69 : ri% rate_fname = trim(filename)
659 : if (dbg) write(*,*) 'done'
660 :
661 :
662 569 : end subroutine rate_from_file
663 :
664 :
665 5 : subroutine read_reaction_parameters(reactionlist_filename, ierr)
666 : use utils_lib
667 : use chem_lib
668 : use chem_def, only: chem_isos, category_id, category_Name
669 : character (len=*), intent(in) :: reactionlist_filename
670 : integer, intent(out) :: ierr
671 :
672 : character (len=256) :: line, filename, rname, cname
673 : integer :: iounit, len, i, j, jj, k, cnt, ir, ic, ii, n, num_reactions
674 : character (len=maxlen_reaction_Info) :: info
675 : logical, parameter :: dbg = .false.
676 :
677 : include 'formats'
678 :
679 5 : ierr = 0
680 :
681 5 : call alloc_and_init_reaction_parameters(ierr)
682 5 : if (ierr /= 0) return
683 :
684 : ! first try the reaction_filename alone
685 5 : filename = trim(reactionlist_filename)
686 5 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
687 5 : if (ierr /= 0) then ! if don't find that file, look in rates_data
688 5 : filename = trim(mesa_data_dir) // '/rates_data/' // trim(reactionlist_filename)
689 5 : ierr = 0
690 5 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
691 5 : if (ierr /= 0) then
692 : write(*,*) &
693 0 : 'failed to open reaction_parameters file ' // trim(reactionlist_filename)
694 0 : return
695 : end if
696 : end if
697 :
698 5 : num_reactions = 0
699 1965 : do cnt = 1, rates_reaction_id_max*10 ! will stop when reach end of file
700 : if (dbg) write(*, *) 'cnt', cnt
701 :
702 1965 : read(iounit,'(a)',iostat=ierr) line
703 1965 : if (ierr /= 0) then
704 : if (dbg) write(*,*) 'reached end of file'
705 : exit
706 : end if
707 1960 : len = len_trim(line)
708 1960 : if (len == 0) then
709 : if (dbg) write(*,*) '(len == 0)'
710 : cycle
711 : end if
712 1770 : if (line(1:1) == '!') then
713 : if (dbg) write(*,*) '(line(1:1) == !)'
714 : cycle
715 : end if
716 :
717 1590 : i = 1; j = 35
718 1590 : rname = line(i:j)
719 :
720 : if (dbg) write(*,*) trim(rname)
721 :
722 1590 : if (line(i:i+1) == 'r ') then
723 0 : call increase_num_reactions(ierr)
724 0 : if (ierr /= 0) then
725 0 : write(*,*) 'FATAL ERROR: rates_def_init failed in increase_num_reactions'
726 0 : return
727 : end if
728 0 : ir = rates_reaction_id_max
729 : else
730 1590 : ir = get_rates_reaction_id(rname)
731 : end if
732 : if (dbg) write(*,*) 'name: ' // trim(rname), ir
733 1590 : if (ir == 0) then
734 50 : call increase_num_reactions(ierr)
735 50 : if (ierr /= 0) then
736 0 : write(*,*) 'FATAL ERROR: rates_def_init failed in increase_num_reactions'
737 0 : return
738 : end if
739 50 : ir = rates_reaction_id_max
740 : if (dbg) write(*,*) 'size(reaction_Name,dim=1), ir', size(reaction_Name,dim=1), ir
741 50 : reaction_Name(ir) = rname
742 50 : call integer_dict_define(reaction_names_dict, reaction_Name(ir), ir, ierr)
743 50 : if (ierr /= 0) then
744 0 : write(*,*) 'FATAL ERROR: rates_def_init failed in integer_dict_define'
745 0 : call mesa_error(__FILE__,__LINE__)
746 : end if
747 : end if
748 :
749 1590 : i = 36; j = 70
750 1590 : call read_inputs
751 1590 : if (ierr /= 0) return
752 :
753 1590 : i = 74; j = 108
754 1590 : call read_outputs
755 1590 : if (ierr /= 0) return
756 :
757 1590 : i = 110; j = 127
758 1590 : call read_Q
759 1590 : if (ierr /= 0) return
760 :
761 1590 : i = 128; j = 143
762 1590 : call read_Qneu
763 1590 : if (ierr /= 0) return
764 :
765 1590 : i = 144; j = 149
766 1590 : call read_ye_rho_exponent1
767 1590 : if (ierr /= 0) return
768 :
769 1590 : i = 150; j = 155
770 1590 : call read_ye_rho_exponent2
771 1590 : if (ierr /= 0) return
772 :
773 1590 : i = 156; j = 160
774 1590 : call read_screening_info(1)
775 1590 : if (ierr /= 0) return
776 :
777 1590 : i = 162; j = 166
778 1590 : call read_screening_info(2)
779 1590 : if (ierr /= 0) return
780 :
781 1590 : i = 168; j = 172
782 1590 : call read_screening_info(3)
783 1590 : if (ierr /= 0) return
784 :
785 1590 : i = 174; j = 178
786 1590 : call read_weak_info(1)
787 1590 : if (ierr /= 0) return
788 :
789 1590 : i = 180; j = 184
790 1590 : call read_weak_info(2)
791 1590 : if (ierr /= 0) return
792 :
793 1590 : if (std_reaction_neuQs(ir) > 0) then
794 200 : weak_reaction_info(1,ir) = reaction_inputs(2,ir)
795 200 : weak_reaction_info(2,ir) = reaction_outputs(2,ir)
796 : end if
797 :
798 : if (std_reaction_neuQs(ir) == 0 .and. &
799 1590 : weak_reaction_info(1,ir) > 0 .and. weak_reaction_info(2,ir) > 0) then
800 : j = do_get_weak_info_list_id( &
801 : chem_isos% name(weak_reaction_info(1,ir)), &
802 0 : chem_isos% name(weak_reaction_info(2,ir)))
803 0 : if (j > 0) std_reaction_neuQs(ir) = weak_info_list_Qneu(j)
804 0 : if (ierr /= 0) return
805 : end if
806 :
807 1590 : if (std_reaction_neuQs(ir) > 0) then
808 200 : call set_weak_lowT_rate(ir, ierr)
809 200 : if (ierr /= 0) return
810 : end if
811 :
812 1590 : i = 190; j = 203
813 1590 : call read_category_id
814 1590 : if (ierr /= 0) return
815 :
816 1590 : i = 208
817 1590 : call read_reaction_Info
818 1590 : if (ierr /= 0) return
819 :
820 : if (dbg) write(*,*)
821 :
822 1590 : num_reactions = cnt
823 :
824 : end do
825 :
826 : if (dbg) write(*,*) 'num_reactions', num_reactions
827 :
828 5 : ierr = 0
829 5 : close(iounit)
830 :
831 5 : num_reactions = rates_reaction_id_max
832 :
833 5 : call check_std_reaction_Qs
834 5 : call check_std_reaction_neuQs
835 5 : call check_reaction_categories
836 5 : call check_reaction_info
837 :
838 10 : if (dbg) call mesa_error(__FILE__,__LINE__,'read_reaction_parameters')
839 :
840 :
841 : contains
842 :
843 :
844 1590 : subroutine read_inputs
845 : if (dbg) write(*,*) ' inputs <' // line(i:j) // '>',i,j
846 :
847 1590 : jj = j
848 1590 : do k = 1, 2*max_num_reaction_inputs-1, 2
849 3565 : j = i; j = i+2
850 3565 : n = read_int()
851 3565 : if (n == 0) exit
852 : if (dbg) write(*,*) 'n <' // line(i:j) // '>', i, j
853 2035 : reaction_inputs(k,ir) = n
854 :
855 : ! fxt
856 2035 : i = j+1; j = i+7
857 :
858 : ! i = j+1; j = i+4
859 :
860 : if (dbg) write(*,*) 'iso <' // line(i:j) // '>', i, j
861 2035 : ii = read_iso()
862 2035 : if (ii <= 0) then
863 0 : write(*,'(a)') 'bad input iso in reaction_parameters file <' // line(i:j) // '>'
864 0 : write(*,'(a)') trim(line)
865 0 : call mesa_error(__FILE__,__LINE__,'read_reaction_parameters')
866 0 : ierr = -1
867 0 : return
868 : end if
869 2035 : i = j+2
870 2035 : reaction_inputs(k+1,ir) = ii
871 1590 : if (dbg) write(*,*) 'in', n, trim(chem_isos% name(ii)), i
872 : end do
873 5 : end subroutine read_inputs
874 :
875 :
876 1590 : subroutine read_outputs
877 : if (dbg) write(*,*) ' outputs: ' // line(i:j)
878 :
879 1590 : jj = j
880 1590 : do k = 1, 2*max_num_reaction_outputs-1, 2
881 3415 : j = i; j = i+2
882 3415 : n = read_int()
883 3415 : if (n == 0) exit
884 : if (dbg) write(*,*) 'n <' // line(i:j) // '>'
885 1825 : reaction_outputs(k,ir) = n
886 :
887 : ! fxt
888 1825 : i = j+1; j = i+7
889 :
890 : ! i = j+1; j = i+4
891 :
892 : if (dbg) write(*,*) 'iso <' // line(i:j) // '>'
893 1825 : ii = read_iso()
894 1825 : if (ii <= 0) then
895 0 : write(*,'(a)') 'bad output iso in reaction_parameters file <' // line(i:j) // '>'
896 0 : write(*,'(a)') trim(line)
897 0 : call mesa_error(__FILE__,__LINE__,'read_reaction_parameters')
898 0 : ierr = -1
899 0 : return
900 : end if
901 1825 : i = j+2
902 1825 : reaction_outputs(k+1,ir) = ii
903 1590 : if (dbg) write(*,*) 'out', n, trim(chem_isos% name(ii))
904 : end do
905 : end subroutine read_outputs
906 :
907 :
908 1590 : subroutine read_Q
909 1590 : if (missing_dbl()) then ! use standard Q from chem_lib
910 1540 : std_reaction_Qs(ir) = get_Qtotal(ir)
911 : else
912 50 : std_reaction_Qs(ir) = read_dbl()
913 : if (dbg) write(*,*) 'std_reaction_Qs(ir)', std_reaction_Qs(ir)
914 : end if
915 1590 : end subroutine read_Q
916 :
917 :
918 1590 : subroutine read_Qneu
919 : include 'formats'
920 1590 : std_reaction_neuQs(ir) = read_dbl()
921 : if (dbg) write(*,1) 'std_reaction_neuQs(ir)', std_reaction_neuQs(ir)
922 1590 : end subroutine read_Qneu
923 :
924 :
925 1590 : subroutine read_ye_rho_exponent1
926 1590 : reaction_ye_rho_exponents(1,ir) = read_int()
927 : if (dbg) write(*,*) 'ye', reaction_ye_rho_exponents(1,ir)
928 1590 : if (reaction_ye_rho_exponents(1,ir) > 2) then
929 0 : write(*,'(a)') 'ERROR: must revise rates_support for large ye exponent'
930 0 : write(*,'(a)') trim(line)
931 0 : call mesa_error(__FILE__,__LINE__,'read_reaction_parameters')
932 : end if
933 1590 : end subroutine read_ye_rho_exponent1
934 :
935 :
936 1590 : subroutine read_ye_rho_exponent2
937 1590 : ii = read_int()
938 1590 : if (ii > 0) ii = ii-1
939 1590 : reaction_ye_rho_exponents(2,ir) = ii
940 1590 : if (ii > 4) then
941 0 : write(*,'(a)') 'ERROR: must revise rates_support for large rho exponent'
942 0 : write(*,'(a)') trim(line)
943 0 : call mesa_error(__FILE__,__LINE__,'read_reaction_parameters')
944 : end if
945 : if (dbg) write(*,*) 'rho', reaction_ye_rho_exponents(2,ir)+1
946 1590 : end subroutine read_ye_rho_exponent2
947 :
948 :
949 4770 : subroutine read_screening_info(which)
950 : integer, intent(in) :: which
951 : logical :: empty
952 : integer :: jj
953 4770 : ii = read_iso()
954 4770 : reaction_screening_info(which,ir) = ii
955 :
956 : !Hack to get around a bug in ifort 17,18, which returns len_trim(line(i:j)) < 0
957 4770 : empty=.true.
958 28620 : do jj=i,j
959 28620 : if(len_trim(line(jj:jj)) /= 0) empty=.false.
960 : end do
961 :
962 4770 : if (ii <= 0 .and. .not. empty) then
963 0 : write(*,'(a)') 'bad iso name for screening in reaction_parameters file <' // line(i:j) // '>'
964 0 : write(*,'(a)') trim(line)
965 0 : call mesa_error(__FILE__,__LINE__,'read_screening_info')
966 0 : ierr = -1
967 0 : return
968 : end if
969 : if (dbg .and. ii > 0) &
970 : write(*,*) 'screen: ' // trim(chem_isos% name(ii)) // ' ' // line(i:j), ii
971 : end subroutine read_screening_info
972 :
973 :
974 3180 : subroutine read_weak_info(which)
975 : integer, intent(in) :: which
976 : logical :: empty
977 : integer :: jj
978 :
979 3180 : ii = read_iso()
980 3180 : weak_reaction_info(which,ir) = ii
981 :
982 : !Hack to get around a bug in ifort 7, 18 which returns len_trim(line(i:j)) < 0
983 3180 : empty=.true.
984 19080 : do jj=i,j
985 19080 : if(len_trim(line(jj:jj)) /= 0) empty=.false.
986 : end do
987 :
988 3180 : if (ii <= 0 .and. .not. empty) then
989 0 : write(*,'(a)') 'bad iso name for weak in reaction_parameters file <' // line(i:j) // '>'
990 0 : write(*,'(a)') trim(line)
991 0 : call mesa_error(__FILE__,__LINE__,'read_reaction_parameters len')
992 0 : ierr = -1
993 0 : return
994 : end if
995 3180 : if (ii > 0) then
996 0 : write(*,*) 'weak: ' // trim(chem_isos% name(ii)) // ' ' // line(i:j), ii
997 0 : write(*,'(a)') trim(line)
998 0 : write(*,'(a)') 'DO NOT USE reactions.list FOR WEAK ISOS; just give Qneu'
999 0 : call mesa_error(__FILE__,__LINE__,'read_reaction_parameters weak')
1000 0 : ierr = -1
1001 0 : return
1002 : end if
1003 : end subroutine read_weak_info
1004 :
1005 :
1006 1590 : subroutine read_category_id
1007 1590 : cname = line(i:j)
1008 1590 : ic = category_id(cname)
1009 1590 : if (ic == 0) then
1010 0 : write(*,*) 'bad category name in reaction_parameters file <' // line(i:j) // '>'
1011 0 : write(*,'(a)') trim(line)
1012 0 : call mesa_error(__FILE__,__LINE__,'read_category_id')
1013 0 : ierr = -1
1014 0 : return
1015 : end if
1016 1590 : reaction_categories(ir) = ic
1017 : if (dbg) write(*,*) 'category: ' // trim(cname), ic, trim(category_name(ic))
1018 : end subroutine read_category_id
1019 :
1020 :
1021 1590 : subroutine read_reaction_Info
1022 1590 : ii = min(maxlen_reaction_Info, len - i + 1)
1023 116070 : do j = 1, maxlen_reaction_Info
1024 116070 : if (j <= ii) then
1025 45655 : info(j:j) = line(i:i)
1026 45655 : i = i+1
1027 : else
1028 68825 : info(j:j) = ' '
1029 : end if
1030 : end do
1031 1590 : reaction_Info(ir) = info
1032 : if (dbg) write(*,*) 'info: ' // trim(reaction_Info(ir))
1033 1590 : end subroutine read_reaction_Info
1034 :
1035 :
1036 11810 : integer function read_iso()
1037 : use chem_def, only: iso_name_length
1038 : character (len=64) :: str
1039 11810 : str = line(i:j)
1040 11810 : read_iso = chem_get_iso_id(str)
1041 11810 : end function read_iso
1042 :
1043 :
1044 10160 : integer function read_int()
1045 : character (len=64) :: str
1046 : integer :: ierr
1047 10160 : str = line(i:j)
1048 10160 : read(str,fmt=*,iostat=ierr) read_int
1049 10160 : if (ierr /= 0) read_int = 0
1050 11810 : end function read_int
1051 :
1052 :
1053 1640 : real(dp) function read_dbl()
1054 : use math_lib, only: str_to_double
1055 : integer :: ierr
1056 1640 : ierr = 0
1057 1640 : if (len_trim(line(i:j)) > 0) then
1058 295 : call str_to_double(line(i:j),read_dbl,ierr)
1059 295 : if (ierr /= 0) read_dbl = 0d0
1060 : else
1061 1345 : read_dbl = 0d0
1062 : end if
1063 1640 : end function read_dbl
1064 :
1065 :
1066 1590 : logical function missing_dbl()
1067 : character (len=64) :: str
1068 1590 : str = line(i:j)
1069 1590 : missing_dbl = (len_trim(str) == 0)
1070 1640 : end function missing_dbl
1071 :
1072 :
1073 5 : subroutine check_std_reaction_Qs
1074 : integer :: i, cnt
1075 5 : cnt = 0
1076 1590 : do i=1,num_reactions
1077 1590 : if (std_reaction_Qs(i) < -1d50) then
1078 0 : write(*,*) 'missing reaction_Q for reaction ' // trim(reaction_Name(i)), i, cnt+1
1079 0 : write(*,*)
1080 0 : cnt = cnt+1
1081 : end if
1082 : end do
1083 5 : if (cnt > 0) call mesa_error(__FILE__,__LINE__,'check_std_reaction_Qs')
1084 5 : end subroutine check_std_reaction_Qs
1085 :
1086 :
1087 5 : subroutine check_std_reaction_neuQs
1088 : integer :: i, cnt
1089 5 : cnt = 0
1090 1590 : do i=1,num_reactions
1091 1590 : if (std_reaction_neuQs(i) < -1d50) then
1092 0 : write(*,*) 'missing std_reaction_neuQs for reaction ' // trim(reaction_Name(i))
1093 0 : write(*,*)
1094 0 : cnt = cnt+1
1095 : end if
1096 : end do
1097 5 : if (cnt > 0) call mesa_error(__FILE__,__LINE__,'check_std_reaction_neuQs')
1098 5 : end subroutine check_std_reaction_neuQs
1099 :
1100 :
1101 5 : subroutine check_reaction_categories
1102 : integer :: cnt, i
1103 5 : cnt = 0
1104 1590 : do i=1,num_reactions
1105 1590 : if (reaction_categories(i) < 0) then
1106 0 : write(*,*) 'missing reaction_category for reaction ' // trim(reaction_Name(i))
1107 0 : write(*,*)
1108 0 : cnt = cnt+1
1109 : end if
1110 : end do
1111 5 : if (cnt > 0) call mesa_error(__FILE__,__LINE__,'check_reaction_categories')
1112 5 : end subroutine check_reaction_categories
1113 :
1114 :
1115 5 : subroutine check_reaction_info
1116 : integer :: i, cnt
1117 5 : cnt = 0
1118 1590 : do i=1,num_reactions
1119 1590 : if (len_trim(reaction_Info(i)) == 0) then
1120 0 : write(*,*) 'missing info for reaction', i
1121 0 : if (i > 1) write(*,*) 'following ' // trim(reaction_Info(i-1))
1122 0 : write(*,*)
1123 0 : cnt = cnt+1
1124 : end if
1125 : end do
1126 5 : if (cnt > 0) call mesa_error(__FILE__,__LINE__,'check_reaction_info')
1127 5 : end subroutine check_reaction_info
1128 :
1129 :
1130 : end subroutine read_reaction_parameters
1131 :
1132 :
1133 68 : subroutine increase_num_reactions(ierr)
1134 : integer, intent(out) :: ierr
1135 :
1136 : integer :: old_max, new_max, i
1137 : type (rate_table_info), pointer :: old_raw_rates_records(:) =>null()
1138 : type (rate_table_info), pointer :: ri =>null()
1139 :
1140 : include 'formats'
1141 :
1142 68 : old_max = rates_reaction_id_max
1143 68 : rates_reaction_id_max = rates_reaction_id_max + 1
1144 :
1145 68 : if (rates_reaction_id_max > size(std_reaction_Qs,dim=1)) then
1146 :
1147 5 : new_max = rates_reaction_id_max*2 + 1000
1148 : !write(*,3) 'increase size', rates_reaction_id_max, new_max
1149 :
1150 5 : old_raw_rates_records => raw_rates_records
1151 5 : allocate(raw_rates_records(new_max))
1152 1540 : do i=1,old_max
1153 1540 : raw_rates_records(i) = old_raw_rates_records(i)
1154 : end do
1155 5 : deallocate(old_raw_rates_records)
1156 :
1157 5 : call grow_reactions_arrays(old_max, new_max, ierr)
1158 5 : if (ierr /= 0) return
1159 : end if
1160 :
1161 68 : ri => raw_rates_records(rates_reaction_id_max)
1162 68 : ri% nT8s = 0
1163 68 : ri% use_rate_table = .false.
1164 68 : ri% need_to_read = .false.
1165 68 : nullify(ri% T8s)
1166 68 : nullify(ri% f1)
1167 :
1168 : end subroutine increase_num_reactions
1169 :
1170 :
1171 5 : subroutine alloc_and_init_reaction_parameters(ierr)
1172 : integer, intent(out) :: ierr
1173 :
1174 : allocate( &
1175 : reaction_Info(rates_reaction_id_max), &
1176 : reaction_categories(rates_reaction_id_max), &
1177 : reaction_is_reverse(rates_reaction_id_max), &
1178 : reaction_reaclib_lo(rates_reaction_id_max), &
1179 : reaction_reaclib_hi(rates_reaction_id_max), &
1180 : reverse_reaction_id(rates_reaction_id_max), &
1181 : reaction_screening_info(3,rates_reaction_id_max), &
1182 : weak_reaction_info(2,rates_reaction_id_max), &
1183 : reaction_ye_rho_exponents(2,rates_reaction_id_max), &
1184 : reaction_inputs(2*max_num_reaction_inputs,rates_reaction_id_max), &
1185 : reaction_outputs(2*max_num_reaction_outputs,rates_reaction_id_max), &
1186 : std_reaction_Qs(rates_reaction_id_max), &
1187 : std_reaction_neuQs(rates_reaction_id_max), &
1188 : weak_lowT_rate(rates_reaction_id_max), &
1189 5 : stat=ierr)
1190 5 : if (ierr /= 0) return
1191 :
1192 1540 : reaction_Info(:) = ''
1193 1540 : reaction_categories(:) = -1
1194 1540 : reaction_is_reverse(:) = 0
1195 1540 : reaction_reaclib_lo(:) = 0
1196 1540 : reaction_reaclib_hi(:) = 0
1197 1540 : reverse_reaction_id(:) = 0
1198 6145 : reaction_screening_info(:,:) = 0
1199 4610 : weak_reaction_info(:,:) = 0
1200 4610 : reaction_ye_rho_exponents(:,:) = 0
1201 10750 : reaction_inputs(:,:) = 0
1202 13820 : reaction_outputs(:,:) = 0
1203 1540 : std_reaction_Qs(:) = -1d99
1204 1540 : std_reaction_neuQs(:) = -1d99
1205 1540 : weak_lowT_rate(:) = -1d99
1206 :
1207 : end subroutine alloc_and_init_reaction_parameters
1208 :
1209 :
1210 1558 : real(dp) function get_Qtotal(ir)
1211 : use chem_lib, only: reaction_Qtotal
1212 : use chem_def, only: chem_isos
1213 : integer, intent(in) :: ir
1214 :
1215 : integer :: num_in, num_out, reactants(100), k, n, i, ii, j
1216 : include 'formats'
1217 1558 : i = 0
1218 2010 : do k = 1, 2*max_num_reaction_inputs-1, 2
1219 3508 : n = reaction_inputs(k,ir)
1220 3508 : if (n == 0) exit
1221 2010 : ii = reaction_inputs(k+1,ir)
1222 5988 : do j = 1, n
1223 2420 : i = i+1
1224 4430 : reactants(i) = ii
1225 : end do
1226 : end do
1227 1558 : num_in = i
1228 3353 : do k = 1, 2*max_num_reaction_outputs-1, 2
1229 3353 : n = reaction_outputs(k,ir)
1230 3353 : if (n == 0) exit
1231 1795 : ii = reaction_outputs(k+1,ir)
1232 5258 : do j = 1, n
1233 1905 : i = i+1
1234 3700 : reactants(i) = ii
1235 : end do
1236 : end do
1237 1558 : num_out = i - num_in
1238 1558 : get_Qtotal = reaction_Qtotal(num_in,num_out,reactants,chem_isos)
1239 1558 : end function get_Qtotal
1240 :
1241 :
1242 20 : subroutine init_rates_info(reactionlist_filename, ierr)
1243 : character (len=*), intent(in) :: reactionlist_filename
1244 : integer, intent(out) :: ierr ! 0 means AOK.
1245 : include 'formats'
1246 :
1247 5 : ierr = 0
1248 5 : call init1_rates_info
1249 :
1250 5 : call start_rates_def_init(ierr)
1251 5 : if (ierr /= 0) then
1252 0 : write(*,*) 'start_rates_def_init failed'
1253 0 : return
1254 : end if
1255 :
1256 5 : call read_reaction_parameters(reactionlist_filename, ierr)
1257 5 : if (ierr /= 0) then
1258 0 : write(*,*) 'rates_init failed in read_reaction_parameters'
1259 0 : return
1260 : end if
1261 :
1262 5 : call finish_rates_def_init
1263 5 : call do_rates_init(ierr)
1264 5 : if (ierr /= 0) then
1265 0 : write(*,*) 'rates_init failed in do_rates_init'
1266 0 : return
1267 : end if
1268 :
1269 1558 : end subroutine init_rates_info
1270 :
1271 :
1272 5 : subroutine init1_rates_info
1273 : use rates_names, only: set_reaction_names
1274 : type (rate_table_info), pointer :: ri =>null()
1275 : integer :: i
1276 5 : rates_reaction_id_max = num_predefined_reactions
1277 : allocate( &
1278 : raw_rates_records(rates_reaction_id_max), &
1279 5 : reaction_Name(rates_reaction_id_max))
1280 1540 : do i = 1, rates_reaction_id_max
1281 1535 : ri => raw_rates_records(i)
1282 1535 : ri% nT8s = 0
1283 1535 : ri% use_rate_table = .false.
1284 1535 : ri% need_to_read = .false.
1285 1535 : nullify(ri% T8s)
1286 1540 : nullify(ri% f1)
1287 : end do
1288 5 : call set_reaction_names
1289 5 : end subroutine init1_rates_info
1290 :
1291 :
1292 69 : integer function lookup_rate_name(str) ! -1 if not found
1293 : use rates_def
1294 : character (len=*), intent(in) :: str
1295 : integer :: i
1296 69 : lookup_rate_name = -1
1297 21297 : do i = 1, rates_reaction_id_max
1298 21297 : if (trim(str) == trim(reaction_Name(i))) then
1299 69 : lookup_rate_name = i
1300 69 : return
1301 : end if
1302 : end do
1303 69 : end function lookup_rate_name
1304 :
1305 5 : subroutine grow_reactions_arrays(old_max, new_max, ierr)
1306 69 : use utils_lib
1307 : integer, intent(in) :: old_max, new_max
1308 : integer, intent(out) :: ierr
1309 :
1310 : character (len=maxlen_reaction_Name), pointer :: new_reaction_Name(:) =>null()
1311 : character (len=maxlen_reaction_Info), pointer :: new_reaction_Info(:) =>null()
1312 : integer :: i
1313 :
1314 : include 'formats'
1315 :
1316 5 : allocate(new_reaction_Name(new_max))
1317 1540 : do i=1,old_max
1318 1540 : new_reaction_Name(i) = reaction_Name(i)
1319 : end do
1320 5 : deallocate(reaction_Name)
1321 5 : reaction_Name => new_reaction_Name
1322 :
1323 5 : allocate(new_reaction_Info(new_max))
1324 1540 : do i=1,old_max
1325 1540 : new_reaction_Info(i) = reaction_Info(i)
1326 : end do
1327 5 : deallocate(reaction_Info)
1328 5 : reaction_Info => new_reaction_Info
1329 :
1330 5 : call realloc_integer(reaction_categories,new_max,ierr); if (ierr /= 0) return
1331 5 : call realloc_integer(reaction_is_reverse,new_max,ierr); if (ierr /= 0) return
1332 5 : call realloc_integer(reaction_reaclib_lo,new_max,ierr); if (ierr /= 0) return
1333 5 : call realloc_integer(reaction_reaclib_hi,new_max,ierr); if (ierr /= 0) return
1334 5 : call realloc_integer(reverse_reaction_id,new_max,ierr); if (ierr /= 0) return
1335 :
1336 5 : call realloc_double(std_reaction_Qs,new_max,ierr); if (ierr /= 0) return
1337 5 : call realloc_double(std_reaction_neuQs,new_max,ierr); if (ierr /= 0) return
1338 :
1339 5 : call realloc_double(weak_lowT_rate,new_max,ierr); if (ierr /= 0) return
1340 :
1341 : call realloc_integer2( &
1342 : reaction_screening_info,size( &
1343 5 : reaction_screening_info,dim=1),new_max,ierr); if (ierr /= 0) return
1344 : call realloc_integer2( &
1345 5 : weak_reaction_info,2,new_max,ierr); if (ierr /= 0) return
1346 : call realloc_integer2( &
1347 5 : reaction_ye_rho_exponents,2,new_max,ierr); if (ierr /= 0) return
1348 : call realloc_integer2( &
1349 5 : reaction_inputs,2*max_num_reaction_inputs,new_max,ierr); if (ierr /= 0) return
1350 : call realloc_integer2( &
1351 5 : reaction_outputs,2*max_num_reaction_outputs,new_max,ierr); if (ierr /= 0) return
1352 :
1353 6550 : reaction_Info(rates_reaction_id_max:new_max) = ''
1354 6550 : reaction_categories(rates_reaction_id_max:new_max) = -1
1355 :
1356 6550 : reaction_is_reverse(rates_reaction_id_max:new_max) = 0
1357 6550 : reaction_reaclib_lo(rates_reaction_id_max:new_max) = 0
1358 6550 : reaction_reaclib_hi(rates_reaction_id_max:new_max) = 0
1359 6550 : reverse_reaction_id(rates_reaction_id_max:new_max) = 0
1360 :
1361 26185 : reaction_screening_info(:,rates_reaction_id_max:new_max) = 0
1362 19640 : weak_reaction_info(:,rates_reaction_id_max:new_max) = 0
1363 19640 : reaction_ye_rho_exponents(:,rates_reaction_id_max:new_max) = 0
1364 45820 : reaction_inputs(:,rates_reaction_id_max:new_max) = 0
1365 58910 : reaction_outputs(:,rates_reaction_id_max:new_max) = 0
1366 6550 : std_reaction_Qs(rates_reaction_id_max:new_max) = -1d99
1367 6550 : std_reaction_neuQs(rates_reaction_id_max:new_max) = -1d99
1368 6550 : weak_lowT_rate(rates_reaction_id_max:new_max) = -1d99
1369 :
1370 40 : end subroutine grow_reactions_arrays
1371 :
1372 :
1373 3 : subroutine free_reaction_arrays()
1374 :
1375 3 : if (ASSOCIATED(reaction_Name)) deallocate(reaction_Name)
1376 3 : if (ASSOCIATED(reaction_Info)) deallocate(reaction_Info)
1377 :
1378 3 : if (ASSOCIATED(reaction_categories)) deallocate(reaction_categories)
1379 3 : if (ASSOCIATED(reaction_is_reverse)) deallocate(reaction_is_reverse)
1380 3 : if (ASSOCIATED(reaction_reaclib_lo)) deallocate(reaction_reaclib_lo)
1381 3 : if (ASSOCIATED(reaction_reaclib_hi)) deallocate(reaction_reaclib_hi)
1382 3 : if (ASSOCIATED(reverse_reaction_id)) deallocate(reverse_reaction_id)
1383 :
1384 3 : if (ASSOCIATED(std_reaction_Qs)) deallocate(std_reaction_Qs)
1385 3 : if (ASSOCIATED(std_reaction_neuQs)) deallocate(std_reaction_neuQs)
1386 :
1387 3 : if (ASSOCIATED(weak_lowT_rate)) deallocate(weak_lowT_rate)
1388 :
1389 3 : if (ASSOCIATED(reaction_screening_info)) deallocate(reaction_screening_info)
1390 3 : if (ASSOCIATED(weak_reaction_info)) deallocate(weak_reaction_info)
1391 3 : if (ASSOCIATED(reaction_ye_rho_exponents)) deallocate(reaction_ye_rho_exponents)
1392 3 : if (ASSOCIATED(reaction_inputs)) deallocate(reaction_inputs)
1393 3 : if (ASSOCIATED(reaction_outputs)) deallocate(reaction_outputs)
1394 :
1395 3 : nullify(reaction_Name)
1396 3 : nullify(reaction_Info)
1397 3 : nullify(reaction_is_reverse)
1398 3 : nullify(reaction_categories)
1399 3 : nullify(reaction_reaclib_lo)
1400 3 : nullify(reaction_reaclib_hi)
1401 3 : nullify(reverse_reaction_id)
1402 3 : nullify(std_reaction_Qs)
1403 3 : nullify(std_reaction_neuQs)
1404 3 : nullify(weak_lowT_rate)
1405 3 : nullify(reaction_screening_info)
1406 3 : nullify(reaction_ye_rho_exponents)
1407 3 : nullify(reaction_inputs)
1408 3 : nullify(reaction_outputs)
1409 :
1410 3 : end subroutine free_reaction_arrays
1411 :
1412 :
1413 5 : subroutine do_rates_init(ierr)
1414 : use ratelib, only: mazurek_init
1415 : integer, intent(out) :: ierr
1416 : ierr = 0
1417 : ! setup interpolation info for mazurek's 1973 fits for the ni56 electron capture rate
1418 5 : call mazurek_init(ierr)
1419 5 : if (ierr /= 0) return
1420 5 : end subroutine do_rates_init
1421 :
1422 : end module rates_initialize
|
