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 3 : 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 3 : call integer_dict_create_hash(reaction_names_dict, ierr)
44 3 : 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 3 : allocate(set(num_chem_isos))
51 :
52 3 : call generate_nuclide_set(chem_isos% name, set)
53 :
54 3 : use_weaklib = .true.
55 3 : call do_extract_rates(set, chem_isos, reaclib_rates, use_weaklib, ierr)
56 3 : deallocate(set)
57 3 : 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 6 : end subroutine finish_rates_def_init
63 :
64 :
65 0 : subroutine do_add_reaction_for_handle(reaction_handle, ierr)
66 : 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 : 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 27 : 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 9 : ierr = 0
160 9 : i = indx
161 9 : reverse = (len_trim(reverse_handle) > 0)
162 9 : r => reaclib_rates
163 :
164 36 : cin(:) = 1
165 45 : cout(:) = 1
166 :
167 : if (dbg) write(*,'(a, 2x, i5)') 'do_add_reaction_from_reaclib ' // trim(reaction_handle), i
168 :
169 9 : if (reverse) then
170 1 : 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 8 : 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 9 : chapter = r% chapter(i)
186 9 : weak = (adjustl(r% reaction_flag(i)) == 'w')
187 :
188 9 : call alloc_reaction_ir(reaction_handle, ir, already_defined, ierr)
189 9 : if (already_defined) return
190 9 : if (ierr /= 0) return
191 :
192 63 : reaction_inputs(:,ir) = 0
193 81 : reaction_outputs(:,ir) = 0
194 :
195 9 : if (.not. reverse) then
196 8 : particles_in = Nin(chapter)
197 8 : particles_out = Nout(chapter)
198 8 : call setup(reaction_inputs(:,ir), num_in, cin, 0, particles_in)
199 8 : call setup(reaction_outputs(:,ir), num_out, cout, particles_in, particles_out)
200 : else
201 1 : particles_out = Nin(chapter)
202 1 : particles_in = Nout(chapter)
203 1 : call setup(reaction_inputs(:,ir), num_in, cin, particles_out, particles_in)
204 1 : call setup(reaction_outputs(:,ir), num_out, cout, 0, particles_out)
205 : end if
206 :
207 : call set_reaction_info( &
208 9 : ir, num_in, num_out, particles_in, particles_out, weak, reaction_handle, ierr)
209 9 : 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 18 : subroutine setup(reaction_test, num, cnt, k, num_particles)
220 : integer :: reaction_test(:), num, cnt(:), k, num_particles
221 : integer :: j
222 18 : reaction_test(1) = 1
223 18 : reaction_test(2) = r% pspecies(k+1,i)
224 18 : num = 1
225 18 : cnt(num) = 1
226 20 : do j=2,num_particles
227 2 : if (r% pspecies(k+j,i) == r% pspecies(k+j-1,i)) then
228 0 : cnt(num) = cnt(num) + 1
229 : else
230 2 : num = num + 1
231 2 : reaction_test(2*num) = r% pspecies(k+j,i)
232 2 : cnt(num) = 1
233 : end if
234 20 : reaction_test(2*num-1) = cnt(num)
235 : end do
236 18 : end subroutine setup
237 :
238 :
239 : end subroutine do_add_reaction_from_reaclib
240 :
241 :
242 9 : 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 9 : ierr = 0
250 9 : ir = get_rates_reaction_id(reaction_handle)
251 9 : if (ir > 0) then
252 0 : already_defined = .true.
253 0 : return
254 : end if
255 9 : !$omp critical (lock_alloc_reaction_ir)
256 9 : ir = get_rates_reaction_id(reaction_handle)
257 9 : if (ir > 0) then
258 0 : already_defined = .true.
259 : else
260 9 : already_defined = .false.
261 : if (dbg) write(*,2) 'call increase_num_reactions', rates_reaction_id_max
262 9 : call increase_num_reactions(ierr)
263 9 : if (ierr == 0) then
264 9 : 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 9 : ierr = 0
288 :
289 9 : reaction_Name(ir) = reaction_handle
290 :
291 : if (dbg) write(*,*) 'call get_Qtotal'
292 9 : std_reaction_Qs(ir) = get_Qtotal(ir)
293 9 : std_reaction_neuQs(ir) = 0d0
294 9 : weak_lowT_rate(ir) = -1d99
295 :
296 9 : iso_in = reaction_inputs(num_in*2,ir)
297 9 : cin1 = reaction_inputs(num_in*2-1,ir)
298 9 : iso_out = reaction_outputs(num_out*2,ir)
299 :
300 9 : 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 9 : 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 9 : if (weak) then
321 7 : weak_reaction_info(1,ir) = iso_in
322 7 : weak_reaction_info(2,ir) = iso_out
323 7 : weak_j = do_get_weak_info_list_id(chem_isos% name(iso_in),chem_isos% name(iso_out))
324 7 : if (weak_j > 0) std_reaction_neuQs(ir) = weak_info_list_Qneu(weak_j)
325 7 : call set_weak_lowT_rate(ir, ierr)
326 7 : if (ierr /= 0) return
327 : else
328 6 : weak_reaction_info(1:2,ir) = 0
329 : end if
330 :
331 9 : reaction_ye_rho_exponents(1,ir) = 0 ! 1 for electron captures, 0 for rest.
332 :
333 9 : if (particles_in > 1) then
334 1 : reaction_screening_info(1,ir) = reaction_inputs(2,ir)
335 1 : if (cin1 > 1) then
336 0 : reaction_screening_info(2,ir) = reaction_screening_info(1,ir)
337 : else
338 1 : reaction_screening_info(2,ir) = reaction_inputs(4,ir)
339 : end if
340 1 : reaction_ye_rho_exponents(2,ir) = particles_in - 1
341 : else
342 24 : reaction_screening_info(1:2,ir) = 0
343 8 : reaction_ye_rho_exponents(2,ir) = 0
344 : end if
345 9 : reaction_screening_info(3,ir) = 0
346 9 : reaction_categories(ir) = iother
347 9 : if (is_pp_reaction(reaction_handle)) then
348 0 : reaction_categories(ir) = ipp
349 9 : else if (is_cno_reaction(reaction_handle)) then
350 6 : reaction_categories(ir) = icno
351 3 : 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 3 : else if (num_in > 1) then
358 1 : 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 1 : 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 1 : reaction_categories(ir) = i_burn_fe
385 : end select
386 2 : else if (particles_in == 1 .and. particles_out == 2 .and. .not. weak) then
387 1 : reaction_categories(ir) = iphoto
388 : end if
389 :
390 9 : reaction_Info(ir) = reaction_handle
391 :
392 : if (dbg) write(*,'(a,3x,i5)') 'call integer_dict_define ' // trim(reaction_handle), ir
393 :
394 9 : call integer_dict_define(reaction_names_dict, reaction_handle, ir, ierr)
395 9 : 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 9 : end subroutine set_reaction_info
401 :
402 :
403 251 : subroutine set_weak_lowT_rate(ir, ierr)
404 : 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 : real(dp) :: half_life, lambda, dlambda_dlnT, rlambda, drlambda_dlnT
414 :
415 : include 'formats'
416 :
417 127 : weak_reaclib_id_i = 0
418 :
419 247 : if (weak_reaction_info(1,ir) == 0 .or. weak_reaction_info(2,ir) == 0) return
420 :
421 : call do_reaclib_indices_for_reaction( &
422 124 : reaction_Name(ir), reaclib_rates, lo, hi, ierr)
423 124 : if (ierr /= 0) then ! not in reaclib
424 117 : ierr = 0
425 : i = do_get_weak_info_list_id( &
426 : chem_isos% name(weak_reaction_info(1,ir)), &
427 117 : chem_isos% name(weak_reaction_info(2,ir)))
428 117 : if (i > 0) then
429 12 : half_life = weak_info_list_halflife(i)
430 12 : if (half_life > 0d0) then
431 12 : weak_lowT_rate(ir) = ln2/half_life
432 12 : return
433 : end if
434 : end if
435 105 : weak_lowT_rate(ir) = -1d-99
436 105 : 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 7 : lambda, dlambda_dlnT, rlambda, drlambda_dlnT, ierr)
443 7 : 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 7 : if (lo <= reaclib_rates% num_from_weaklib) then
449 4 : if (.not. reaclib_rates% also_in_reaclib(lo)) lambda = -1
450 : end if
451 7 : 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 7 : chem_isos% name(weak_reaction_info(2,ir)))
456 7 : if (weak_reaclib_id_i == 0) return
457 4 : weak_reaclib_id(weak_reaclib_id_i) = ir
458 :
459 : end subroutine set_weak_lowT_rate
460 :
461 :
462 9 : 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 9 : reaction_handle == 'r_b8_wk_he4_he4'
475 9 : end function is_pp_reaction
476 :
477 :
478 9 : 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 9 : reaction_handle == 'r_ne18_wk_f18'
502 9 : end function is_cno_reaction
503 :
504 :
505 2 : subroutine free_raw_rates_records
506 : type (rate_table_info), pointer :: ri =>null()
507 : integer :: i
508 2 : if (ASSOCIATED(raw_rates_records)) then
509 636 : do i = 1, rates_reaction_id_max
510 634 : ri => raw_rates_records(i)
511 634 : if (ASSOCIATED(ri% T8s)) deallocate(ri% T8s)
512 636 : if (ASSOCIATED(ri% f1)) deallocate(ri% f1)
513 : end do
514 2 : deallocate(raw_rates_records)
515 : end if
516 2 : end subroutine free_raw_rates_records
517 :
518 :
519 11 : 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 11 : ierr = 0
537 :
538 : ! first try local rate_tables_dir
539 11 : dir = rates_table_dir
540 11 : filename = trim(dir) // '/rate_list.txt'
541 11 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
542 11 : if (ierr /= 0) then ! if don't find that file, look in rates_dir
543 2 : dir = trim(rates_dir) // '/rate_tables'
544 2 : filename = trim(dir) // '/rate_list.txt'
545 2 : ierr = 0
546 2 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
547 2 : 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 11 : rates_table_dir = dir
553 :
554 11 : n = 0
555 11 : i = 0
556 :
557 : if (dbg) write(*,*) 'read rate list file ' // trim(filename)
558 :
559 40 : rate_loop: do
560 51 : t = token(iounit, n, i, buffer, rate_name)
561 51 : if (t == eof_token) exit rate_loop
562 40 : 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 40 : t = token(iounit, n, i, buffer, rate_fname)
567 40 : if (t /= string_token) then
568 0 : call error; return
569 : end if
570 :
571 40 : call rate_from_file(rate_name, rate_fname, ierr)
572 51 : if (ierr /= 0) call error()
573 :
574 : end do rate_loop
575 :
576 11 : 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 540 : 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 540 : ierr = 0
641 :
642 40 : if (len_trim(rate_name) == 0 .or. len_trim(filename)==0) return
643 :
644 40 : ir = lookup_rate_name(rate_name)
645 40 : 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 40 : ri => raw_rates_records(ir)
656 40 : ri% use_rate_table = .true.
657 40 : ri% need_to_read = .true.
658 40 : ri% rate_fname = trim(filename)
659 : if (dbg) write(*,*) 'done'
660 :
661 :
662 540 : end subroutine rate_from_file
663 :
664 :
665 3 : 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 : ierr = 0
680 :
681 3 : call alloc_and_init_reaction_parameters(ierr)
682 3 : if (ierr /= 0) return
683 :
684 : ! first try the reaction_filename alone
685 3 : filename = trim(reactionlist_filename)
686 3 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
687 3 : if (ierr /= 0) then ! if don't find that file, look in rates_data
688 3 : filename = trim(mesa_data_dir) // '/rates_data/' // trim(reactionlist_filename)
689 3 : ierr = 0
690 3 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
691 3 : 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 3 : num_reactions = 0
699 1179 : do cnt = 1, rates_reaction_id_max*10 ! will stop when reach end of file
700 : if (dbg) write(*, *) 'cnt', cnt
701 :
702 1179 : read(iounit,'(a)',iostat=ierr) line
703 1179 : if (ierr /= 0) then
704 : if (dbg) write(*,*) 'reached end of file'
705 : exit
706 : end if
707 1176 : len = len_trim(line)
708 1176 : if (len == 0) then
709 : if (dbg) write(*,*) '(len == 0)'
710 : cycle
711 : end if
712 1062 : if (line(1:1) == '!') then
713 : if (dbg) write(*,*) '(line(1:1) == !)'
714 : cycle
715 : end if
716 :
717 954 : i = 1; j = 35
718 954 : rname = line(i:j)
719 :
720 : if (dbg) write(*,*) trim(rname)
721 :
722 954 : 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 954 : ir = get_rates_reaction_id(rname)
731 : end if
732 : if (dbg) write(*,*) 'name: ' // trim(rname), ir
733 954 : if (ir == 0) then
734 30 : call increase_num_reactions(ierr)
735 30 : if (ierr /= 0) then
736 0 : write(*,*) 'FATAL ERROR: rates_def_init failed in increase_num_reactions'
737 0 : return
738 : end if
739 30 : ir = rates_reaction_id_max
740 : if (dbg) write(*,*) 'size(reaction_Name,dim=1), ir', size(reaction_Name,dim=1), ir
741 30 : reaction_Name(ir) = rname
742 30 : call integer_dict_define(reaction_names_dict, reaction_Name(ir), ir, ierr)
743 30 : 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 954 : i = 36; j = 70
750 954 : call read_inputs
751 954 : if (ierr /= 0) return
752 :
753 954 : i = 74; j = 108
754 954 : call read_outputs
755 954 : if (ierr /= 0) return
756 :
757 954 : i = 110; j = 127
758 954 : call read_Q
759 954 : if (ierr /= 0) return
760 :
761 954 : i = 128; j = 143
762 954 : call read_Qneu
763 954 : if (ierr /= 0) return
764 :
765 954 : i = 144; j = 149
766 954 : call read_ye_rho_exponent1
767 954 : if (ierr /= 0) return
768 :
769 954 : i = 150; j = 155
770 954 : call read_ye_rho_exponent2
771 954 : if (ierr /= 0) return
772 :
773 954 : i = 156; j = 160
774 954 : call read_screening_info(1)
775 954 : if (ierr /= 0) return
776 :
777 954 : i = 162; j = 166
778 954 : call read_screening_info(2)
779 954 : if (ierr /= 0) return
780 :
781 954 : i = 168; j = 172
782 954 : call read_screening_info(3)
783 954 : if (ierr /= 0) return
784 :
785 954 : i = 174; j = 178
786 954 : call read_weak_info(1)
787 954 : if (ierr /= 0) return
788 :
789 954 : i = 180; j = 184
790 954 : call read_weak_info(2)
791 954 : if (ierr /= 0) return
792 :
793 954 : if (std_reaction_neuQs(ir) > 0) then
794 120 : weak_reaction_info(1,ir) = reaction_inputs(2,ir)
795 120 : weak_reaction_info(2,ir) = reaction_outputs(2,ir)
796 : end if
797 :
798 : if (std_reaction_neuQs(ir) == 0 .and. &
799 954 : 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 954 : if (std_reaction_neuQs(ir) > 0) then
808 120 : call set_weak_lowT_rate(ir, ierr)
809 120 : if (ierr /= 0) return
810 : end if
811 :
812 954 : i = 190; j = 203
813 954 : call read_category_id
814 954 : if (ierr /= 0) return
815 :
816 954 : i = 208
817 954 : call read_reaction_Info
818 954 : if (ierr /= 0) return
819 :
820 : if (dbg) write(*,*)
821 :
822 954 : num_reactions = cnt
823 :
824 : end do
825 :
826 : if (dbg) write(*,*) 'num_reactions', num_reactions
827 :
828 3 : ierr = 0
829 3 : close(iounit)
830 :
831 3 : num_reactions = rates_reaction_id_max
832 :
833 3 : call check_std_reaction_Qs
834 3 : call check_std_reaction_neuQs
835 3 : call check_reaction_categories
836 3 : call check_reaction_info
837 :
838 3 : if (dbg) call mesa_error(__FILE__,__LINE__,'read_reaction_parameters')
839 :
840 :
841 : contains
842 :
843 :
844 954 : subroutine read_inputs
845 : if (dbg) write(*,*) ' inputs <' // line(i:j) // '>',i,j
846 :
847 954 : jj = j
848 954 : do k = 1, 2*max_num_reaction_inputs-1, 2
849 2139 : j = i; j = i+2
850 2139 : n = read_int()
851 2139 : if (n == 0) exit
852 : if (dbg) write(*,*) 'n <' // line(i:j) // '>', i, j
853 1221 : reaction_inputs(k,ir) = n
854 :
855 : ! fxt
856 1221 : 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 1221 : ii = read_iso()
862 1221 : 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 : ierr = -1
867 : return
868 : end if
869 1221 : i = j+2
870 1221 : reaction_inputs(k+1,ir) = ii
871 954 : if (dbg) write(*,*) 'in', n, trim(chem_isos% name(ii)), i
872 : end do
873 : end subroutine read_inputs
874 :
875 :
876 954 : subroutine read_outputs
877 : if (dbg) write(*,*) ' outputs: ' // line(i:j)
878 :
879 954 : jj = j
880 954 : do k = 1, 2*max_num_reaction_outputs-1, 2
881 2049 : j = i; j = i+2
882 2049 : n = read_int()
883 2049 : if (n == 0) exit
884 : if (dbg) write(*,*) 'n <' // line(i:j) // '>'
885 1095 : reaction_outputs(k,ir) = n
886 :
887 : ! fxt
888 1095 : i = j+1; j = i+7
889 :
890 : ! i = j+1; j = i+4
891 :
892 : if (dbg) write(*,*) 'iso <' // line(i:j) // '>'
893 1095 : ii = read_iso()
894 1095 : 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 : ierr = -1
899 : return
900 : end if
901 1095 : i = j+2
902 1095 : reaction_outputs(k+1,ir) = ii
903 954 : if (dbg) write(*,*) 'out', n, trim(chem_isos% name(ii))
904 : end do
905 : end subroutine read_outputs
906 :
907 :
908 954 : subroutine read_Q
909 954 : if (missing_dbl()) then ! use standard Q from chem_lib
910 924 : std_reaction_Qs(ir) = get_Qtotal(ir)
911 : else
912 30 : std_reaction_Qs(ir) = read_dbl()
913 : if (dbg) write(*,*) 'std_reaction_Qs(ir)', std_reaction_Qs(ir)
914 : end if
915 954 : end subroutine read_Q
916 :
917 :
918 954 : subroutine read_Qneu
919 : include 'formats'
920 954 : std_reaction_neuQs(ir) = read_dbl()
921 : if (dbg) write(*,1) 'std_reaction_neuQs(ir)', std_reaction_neuQs(ir)
922 954 : end subroutine read_Qneu
923 :
924 :
925 954 : subroutine read_ye_rho_exponent1
926 954 : reaction_ye_rho_exponents(1,ir) = read_int()
927 : if (dbg) write(*,*) 'ye', reaction_ye_rho_exponents(1,ir)
928 954 : 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 954 : end subroutine read_ye_rho_exponent1
934 :
935 :
936 954 : subroutine read_ye_rho_exponent2
937 954 : ii = read_int()
938 954 : if (ii > 0) ii = ii-1
939 954 : reaction_ye_rho_exponents(2,ir) = ii
940 954 : 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 954 : end subroutine read_ye_rho_exponent2
947 :
948 :
949 2862 : subroutine read_screening_info(which)
950 : integer, intent(in) :: which
951 : logical :: empty
952 : integer :: jj
953 2862 : ii = read_iso()
954 2862 : 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 2862 : empty=.true.
958 17172 : do jj=i,j
959 17172 : if(len_trim(line(jj:jj)) /= 0) empty=.false.
960 : end do
961 :
962 2862 : 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 : ierr = -1
967 : 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 1908 : subroutine read_weak_info(which)
975 : integer, intent(in) :: which
976 : logical :: empty
977 : integer :: jj
978 :
979 1908 : ii = read_iso()
980 1908 : 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 1908 : empty=.true.
984 11448 : do jj=i,j
985 11448 : if(len_trim(line(jj:jj)) /= 0) empty=.false.
986 : end do
987 :
988 1908 : 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 : ierr = -1
993 : return
994 : end if
995 1908 : 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 : ierr = -1
1001 : return
1002 : end if
1003 : end subroutine read_weak_info
1004 :
1005 :
1006 954 : subroutine read_category_id
1007 954 : cname = line(i:j)
1008 954 : ic = category_id(cname)
1009 954 : 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 : ierr = -1
1014 : return
1015 : end if
1016 954 : 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 954 : subroutine read_reaction_Info
1022 954 : ii = min(maxlen_reaction_Info, len - i + 1)
1023 69642 : do j = 1, maxlen_reaction_Info
1024 69642 : if (j <= ii) then
1025 27393 : info(j:j) = line(i:i)
1026 27393 : i = i+1
1027 : else
1028 41295 : info(j:j) = ' '
1029 : end if
1030 : end do
1031 954 : reaction_Info(ir) = info
1032 : if (dbg) write(*,*) 'info: ' // trim(reaction_Info(ir))
1033 954 : end subroutine read_reaction_Info
1034 :
1035 :
1036 7086 : integer function read_iso()
1037 : use chem_def, only: iso_name_length
1038 : character (len=64) :: str
1039 7086 : str = line(i:j)
1040 7086 : read_iso = chem_get_iso_id(str)
1041 7086 : end function read_iso
1042 :
1043 :
1044 6096 : integer function read_int()
1045 : character (len=64) :: str
1046 : integer :: ierr
1047 6096 : str = line(i:j)
1048 6096 : read(str,fmt=*,iostat=ierr) read_int
1049 6096 : if (ierr /= 0) read_int = 0
1050 6096 : end function read_int
1051 :
1052 :
1053 984 : real(dp) function read_dbl()
1054 : use math_lib, only: str_to_double
1055 : integer :: ierr
1056 984 : ierr = 0
1057 984 : if (len_trim(line(i:j)) > 0) then
1058 177 : call str_to_double(line(i:j),read_dbl,ierr)
1059 177 : if (ierr /= 0) read_dbl = 0d0
1060 : else
1061 807 : read_dbl = 0d0
1062 : end if
1063 984 : end function read_dbl
1064 :
1065 :
1066 954 : logical function missing_dbl()
1067 : character (len=64) :: str
1068 954 : str = line(i:j)
1069 954 : missing_dbl = (len_trim(str) == 0)
1070 954 : end function missing_dbl
1071 :
1072 :
1073 3 : subroutine check_std_reaction_Qs
1074 : integer :: i, cnt
1075 3 : cnt = 0
1076 954 : do i=1,num_reactions
1077 954 : 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 3 : if (cnt > 0) call mesa_error(__FILE__,__LINE__,'check_std_reaction_Qs')
1084 3 : end subroutine check_std_reaction_Qs
1085 :
1086 :
1087 3 : subroutine check_std_reaction_neuQs
1088 : integer :: i, cnt
1089 3 : cnt = 0
1090 954 : do i=1,num_reactions
1091 954 : 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 3 : if (cnt > 0) call mesa_error(__FILE__,__LINE__,'check_std_reaction_neuQs')
1098 3 : end subroutine check_std_reaction_neuQs
1099 :
1100 :
1101 3 : subroutine check_reaction_categories
1102 : integer :: cnt, i
1103 3 : cnt = 0
1104 954 : do i=1,num_reactions
1105 954 : 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 3 : if (cnt > 0) call mesa_error(__FILE__,__LINE__,'check_reaction_categories')
1112 3 : end subroutine check_reaction_categories
1113 :
1114 :
1115 3 : subroutine check_reaction_info
1116 : integer :: i, cnt
1117 3 : cnt = 0
1118 954 : do i=1,num_reactions
1119 954 : 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 3 : if (cnt > 0) call mesa_error(__FILE__,__LINE__,'check_reaction_info')
1127 3 : end subroutine check_reaction_info
1128 :
1129 :
1130 : end subroutine read_reaction_parameters
1131 :
1132 :
1133 39 : 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 39 : old_max = rates_reaction_id_max
1143 39 : rates_reaction_id_max = rates_reaction_id_max + 1
1144 :
1145 39 : if (rates_reaction_id_max > size(std_reaction_Qs,dim=1)) then
1146 :
1147 3 : new_max = rates_reaction_id_max*2 + 1000
1148 : !write(*,3) 'increase size', rates_reaction_id_max, new_max
1149 :
1150 3 : old_raw_rates_records => raw_rates_records
1151 3 : allocate(raw_rates_records(new_max))
1152 924 : do i=1,old_max
1153 924 : raw_rates_records(i) = old_raw_rates_records(i)
1154 : end do
1155 3 : deallocate(old_raw_rates_records)
1156 :
1157 3 : call grow_reactions_arrays(old_max, new_max, ierr)
1158 3 : if (ierr /= 0) return
1159 : end if
1160 :
1161 39 : ri => raw_rates_records(rates_reaction_id_max)
1162 39 : ri% nT8s = 0
1163 39 : ri% use_rate_table = .false.
1164 39 : ri% need_to_read = .false.
1165 39 : nullify(ri% T8s)
1166 39 : nullify(ri% f1)
1167 :
1168 : end subroutine increase_num_reactions
1169 :
1170 :
1171 3 : 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 3 : stat=ierr)
1190 3 : if (ierr /= 0) return
1191 :
1192 924 : reaction_Info(:) = ''
1193 924 : reaction_categories(:) = -1
1194 924 : reaction_is_reverse(:) = 0
1195 924 : reaction_reaclib_lo(:) = 0
1196 924 : reaction_reaclib_hi(:) = 0
1197 924 : reverse_reaction_id(:) = 0
1198 3687 : reaction_screening_info(:,:) = 0
1199 2766 : weak_reaction_info(:,:) = 0
1200 2766 : reaction_ye_rho_exponents(:,:) = 0
1201 6450 : reaction_inputs(:,:) = 0
1202 8292 : reaction_outputs(:,:) = 0
1203 924 : std_reaction_Qs(:) = -1d99
1204 924 : std_reaction_neuQs(:) = -1d99
1205 924 : weak_lowT_rate(:) = -1d99
1206 :
1207 : end subroutine alloc_and_init_reaction_parameters
1208 :
1209 :
1210 933 : 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 933 : i = 0
1218 2137 : do k = 1, 2*max_num_reaction_inputs-1, 2
1219 2101 : n = reaction_inputs(k,ir)
1220 2101 : if (n == 0) exit
1221 1204 : ii = reaction_inputs(k+1,ir)
1222 3587 : do j = 1, n
1223 1450 : i = i+1
1224 2654 : reactants(i) = ii
1225 : end do
1226 : end do
1227 933 : num_in = i
1228 2008 : do k = 1, 2*max_num_reaction_outputs-1, 2
1229 2008 : n = reaction_outputs(k,ir)
1230 2008 : if (n == 0) exit
1231 1075 : ii = reaction_outputs(k+1,ir)
1232 3149 : do j = 1, n
1233 1141 : i = i+1
1234 2216 : reactants(i) = ii
1235 : end do
1236 : end do
1237 933 : num_out = i - num_in
1238 933 : get_Qtotal = reaction_Qtotal(num_in,num_out,reactants,chem_isos)
1239 933 : end function get_Qtotal
1240 :
1241 :
1242 12 : 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 3 : ierr = 0
1248 3 : call init1_rates_info
1249 :
1250 3 : call start_rates_def_init(ierr)
1251 3 : if (ierr /= 0) then
1252 0 : write(*,*) 'start_rates_def_init failed'
1253 0 : return
1254 : end if
1255 :
1256 3 : call read_reaction_parameters(reactionlist_filename, ierr)
1257 3 : if (ierr /= 0) then
1258 0 : write(*,*) 'rates_init failed in read_reaction_parameters'
1259 0 : return
1260 : end if
1261 :
1262 3 : call finish_rates_def_init
1263 3 : call do_rates_init(ierr)
1264 3 : if (ierr /= 0) then
1265 0 : write(*,*) 'rates_init failed in do_rates_init'
1266 0 : return
1267 : end if
1268 :
1269 : end subroutine init_rates_info
1270 :
1271 :
1272 3 : subroutine init1_rates_info
1273 : use rates_names, only: set_reaction_names
1274 : type (rate_table_info), pointer :: ri =>null()
1275 : integer :: i
1276 3 : rates_reaction_id_max = num_predefined_reactions
1277 : allocate( &
1278 : raw_rates_records(rates_reaction_id_max), &
1279 3 : reaction_Name(rates_reaction_id_max))
1280 924 : do i = 1, rates_reaction_id_max
1281 921 : ri => raw_rates_records(i)
1282 921 : ri% nT8s = 0
1283 921 : ri% use_rate_table = .false.
1284 921 : ri% need_to_read = .false.
1285 921 : nullify(ri% T8s)
1286 924 : nullify(ri% f1)
1287 : end do
1288 3 : call set_reaction_names
1289 3 : end subroutine init1_rates_info
1290 :
1291 :
1292 40 : integer function lookup_rate_name(str) ! -1 if not found
1293 : use rates_def
1294 : character (len=*), intent(in) :: str
1295 : integer :: i
1296 40 : lookup_rate_name = -1
1297 12359 : do i = 1, rates_reaction_id_max
1298 12359 : if (trim(str) == trim(reaction_Name(i))) then
1299 40 : lookup_rate_name = i
1300 40 : return
1301 : end if
1302 : end do
1303 : end function lookup_rate_name
1304 :
1305 3 : subroutine grow_reactions_arrays(old_max, new_max, ierr)
1306 : 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 3 : allocate(new_reaction_Name(new_max))
1317 924 : do i=1,old_max
1318 924 : new_reaction_Name(i) = reaction_Name(i)
1319 : end do
1320 3 : deallocate(reaction_Name)
1321 3 : reaction_Name => new_reaction_Name
1322 :
1323 3 : allocate(new_reaction_Info(new_max))
1324 924 : do i=1,old_max
1325 924 : new_reaction_Info(i) = reaction_Info(i)
1326 : end do
1327 3 : deallocate(reaction_Info)
1328 3 : reaction_Info => new_reaction_Info
1329 :
1330 3 : call realloc_integer(reaction_categories,new_max,ierr); if (ierr /= 0) return
1331 3 : call realloc_integer(reaction_is_reverse,new_max,ierr); if (ierr /= 0) return
1332 3 : call realloc_integer(reaction_reaclib_lo,new_max,ierr); if (ierr /= 0) return
1333 3 : call realloc_integer(reaction_reaclib_hi,new_max,ierr); if (ierr /= 0) return
1334 3 : call realloc_integer(reverse_reaction_id,new_max,ierr); if (ierr /= 0) return
1335 :
1336 3 : call realloc_double(std_reaction_Qs,new_max,ierr); if (ierr /= 0) return
1337 3 : call realloc_double(std_reaction_neuQs,new_max,ierr); if (ierr /= 0) return
1338 :
1339 3 : call realloc_double(weak_lowT_rate,new_max,ierr); if (ierr /= 0) return
1340 :
1341 : call realloc_integer2( &
1342 : reaction_screening_info,size( &
1343 3 : reaction_screening_info,dim=1),new_max,ierr); if (ierr /= 0) return
1344 : call realloc_integer2( &
1345 3 : weak_reaction_info,2,new_max,ierr); if (ierr /= 0) return
1346 : call realloc_integer2( &
1347 3 : reaction_ye_rho_exponents,2,new_max,ierr); if (ierr /= 0) return
1348 : call realloc_integer2( &
1349 3 : reaction_inputs,2*max_num_reaction_inputs,new_max,ierr); if (ierr /= 0) return
1350 : call realloc_integer2( &
1351 3 : reaction_outputs,2*max_num_reaction_outputs,new_max,ierr); if (ierr /= 0) return
1352 :
1353 3930 : reaction_Info(rates_reaction_id_max:new_max) = ''
1354 3930 : reaction_categories(rates_reaction_id_max:new_max) = -1
1355 :
1356 3930 : reaction_is_reverse(rates_reaction_id_max:new_max) = 0
1357 3930 : reaction_reaclib_lo(rates_reaction_id_max:new_max) = 0
1358 3930 : reaction_reaclib_hi(rates_reaction_id_max:new_max) = 0
1359 3930 : reverse_reaction_id(rates_reaction_id_max:new_max) = 0
1360 :
1361 15711 : reaction_screening_info(:,rates_reaction_id_max:new_max) = 0
1362 11784 : weak_reaction_info(:,rates_reaction_id_max:new_max) = 0
1363 11784 : reaction_ye_rho_exponents(:,rates_reaction_id_max:new_max) = 0
1364 27492 : reaction_inputs(:,rates_reaction_id_max:new_max) = 0
1365 35346 : reaction_outputs(:,rates_reaction_id_max:new_max) = 0
1366 3930 : std_reaction_Qs(rates_reaction_id_max:new_max) = -1d99
1367 3930 : std_reaction_neuQs(rates_reaction_id_max:new_max) = -1d99
1368 3930 : weak_lowT_rate(rates_reaction_id_max:new_max) = -1d99
1369 :
1370 24 : end subroutine grow_reactions_arrays
1371 :
1372 :
1373 2 : subroutine free_reaction_arrays()
1374 :
1375 2 : if (ASSOCIATED(reaction_Name)) deallocate(reaction_Name)
1376 2 : if (ASSOCIATED(reaction_Info)) deallocate(reaction_Info)
1377 :
1378 2 : if (ASSOCIATED(reaction_categories)) deallocate(reaction_categories)
1379 2 : if (ASSOCIATED(reaction_is_reverse)) deallocate(reaction_is_reverse)
1380 2 : if (ASSOCIATED(reaction_reaclib_lo)) deallocate(reaction_reaclib_lo)
1381 2 : if (ASSOCIATED(reaction_reaclib_hi)) deallocate(reaction_reaclib_hi)
1382 2 : if (ASSOCIATED(reverse_reaction_id)) deallocate(reverse_reaction_id)
1383 :
1384 2 : if (ASSOCIATED(std_reaction_Qs)) deallocate(std_reaction_Qs)
1385 2 : if (ASSOCIATED(std_reaction_neuQs)) deallocate(std_reaction_neuQs)
1386 :
1387 2 : if (ASSOCIATED(weak_lowT_rate)) deallocate(weak_lowT_rate)
1388 :
1389 2 : if (ASSOCIATED(reaction_screening_info)) deallocate(reaction_screening_info)
1390 2 : if (ASSOCIATED(weak_reaction_info)) deallocate(weak_reaction_info)
1391 2 : if (ASSOCIATED(reaction_ye_rho_exponents)) deallocate(reaction_ye_rho_exponents)
1392 2 : if (ASSOCIATED(reaction_inputs)) deallocate(reaction_inputs)
1393 2 : if (ASSOCIATED(reaction_outputs)) deallocate(reaction_outputs)
1394 :
1395 2 : nullify(reaction_Name)
1396 2 : nullify(reaction_Info)
1397 2 : nullify(reaction_is_reverse)
1398 2 : nullify(reaction_categories)
1399 2 : nullify(reaction_reaclib_lo)
1400 2 : nullify(reaction_reaclib_hi)
1401 2 : nullify(reverse_reaction_id)
1402 2 : nullify(std_reaction_Qs)
1403 2 : nullify(std_reaction_neuQs)
1404 2 : nullify(weak_lowT_rate)
1405 2 : nullify(reaction_screening_info)
1406 2 : nullify(reaction_ye_rho_exponents)
1407 2 : nullify(reaction_inputs)
1408 2 : nullify(reaction_outputs)
1409 :
1410 2 : end subroutine free_reaction_arrays
1411 :
1412 :
1413 3 : 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 3 : call mazurek_init(ierr)
1419 : if (ierr /= 0) return
1420 : end subroutine do_rates_init
1421 :
1422 : end module rates_initialize
|
