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_def
21 : use utils_def, only: integer_dict
22 : use const_def, only: dp
23 : use chem_def, only: iso_name_length, nuclide_data, npart
24 : use auto_diff
25 :
26 : implicit none
27 :
28 : ! weaklib
29 :
30 : character (len=256) :: weak_data_dir
31 :
32 : ! ecapture
33 :
34 : character (len=1000) :: ecapture_states_file
35 : character (len=1000) :: ecapture_transitions_file
36 :
37 : ! reaclib
38 :
39 :
40 : integer, parameter :: max_nreaclib=100000
41 : integer, parameter :: max_species_per_reaction=6
42 : integer, parameter :: ncoefficients=7
43 : integer, parameter :: nchapters=11
44 : integer, parameter :: ninverse_coeff = 2
45 : integer, parameter :: max_terms_per_rate = 20
46 : integer, parameter :: max_id_length = 36
47 :
48 :
49 : ! i/o parameters
50 : integer, parameter :: internal_format = 0
51 : integer, parameter :: pretty_print_format = 1
52 : integer, parameter :: short_format = 2
53 :
54 :
55 : ! flags for reaction types -- values are chapter numbers
56 : integer, parameter :: &
57 : r_one_one = 1, &
58 : r_one_two = 2, &
59 : r_one_three = 3, &
60 : r_two_one = 4, &
61 : r_two_two = 5, &
62 : r_two_three = 6, &
63 : r_two_four = 7, &
64 : r_three_one = 8, &
65 : r_three_two = 9, &
66 : r_four_two = 10, &
67 : r_one_four = 11
68 :
69 :
70 : integer, dimension(nchapters), parameter :: Nin = [1, 1, 1, 2, 2, 2, 2, 3, 3, 4, 1]
71 : integer, dimension(nchapters), parameter :: Nout = [1, 2, 3, 1, 2, 3, 4, 1, 2, 2, 4]
72 :
73 :
74 : type reaclib_data
75 : integer, dimension(:), allocatable :: chapter
76 : character(len=iso_name_length), dimension(:,:), allocatable :: species
77 : character(len=iso_name_length), dimension(:), allocatable :: label
78 : character, dimension(:), allocatable :: reaction_flag
79 : character, dimension(:), allocatable :: reverse_flag
80 : real(dp), dimension(:), allocatable :: Qvalue
81 : real(dp), dimension(:,:), allocatable :: coefficients
82 : end type reaclib_data
83 :
84 :
85 : type reaction_data
86 :
87 : integer :: nreactions
88 : integer :: nchapters_included
89 : integer, dimension(nchapters) :: chapters_present
90 : type(nuclide_data), pointer :: nuclides=>NULL()
91 :
92 : integer :: num_from_weaklib
93 : integer, dimension(:), pointer :: weaklib_ids=>NULL() ! (num_from_weaklib)
94 : logical, dimension(:), pointer :: also_in_reaclib=>NULL() ! (num_from_weaklib)
95 :
96 : integer, dimension(2, nchapters) :: bookmarks
97 : character(len=max_id_length), dimension(:), pointer :: reaction_handle=>NULL(), reverse_handle=>NULL()
98 : integer, dimension(:), pointer :: category=>NULL()
99 : integer, dimension(:), pointer :: chapter=>NULL()
100 : integer, dimension(:, :), pointer :: pspecies=>NULL()
101 : character, dimension(:), pointer :: reaction_flag=>NULL()
102 : real(dp), dimension(:), pointer :: weight=>NULL()
103 : real(dp), dimension(:), pointer :: weight_reverse=>NULL()
104 : real(dp), dimension(:, :), pointer :: coefficients=>NULL()
105 : real(dp), dimension(:), pointer :: weak_mask=>NULL()
106 : real(dp), dimension(:, :), pointer :: inverse_coefficients=>NULL()
107 : integer, dimension(:), pointer :: inverse_exp=>NULL()
108 : real(dp), dimension(:, :), pointer :: inverse_part=>NULL()
109 : real(dp), dimension(:), pointer :: Q=>NULL()
110 : real(dp), dimension(:), pointer :: Qneu=>NULL()
111 : type (integer_dict), pointer :: reaction_dict=>NULL(), reverse_dict=>NULL()
112 :
113 : end type reaction_data
114 :
115 :
116 : character (len=256) :: reaclib_dir, reaclib_filename
117 :
118 :
119 : ! reactions information
120 :
121 : integer, parameter :: maxlen_reaction_Name = 32
122 : character (len=maxlen_reaction_Name), pointer :: reaction_Name(:)=>NULL() ! (rates_reaction_id_max)
123 :
124 : integer, parameter :: maxlen_reaction_Info = 72
125 : character (len=maxlen_reaction_Info), pointer :: reaction_Info(:)=>NULL() ! (rates_reaction_id_max)
126 :
127 : real(dp), pointer :: std_reaction_Qs(:)=>NULL() ! (rates_reaction_id_max)
128 : ! set at initialization; read-only afterwards.
129 : ! avg energy including neutrinos
130 :
131 : real(dp), pointer :: std_reaction_neuQs(:)=>NULL() ! (rates_reaction_id_max)
132 : ! set at initialization; read-only afterwards.
133 : ! avg neutrino loss
134 :
135 : real(dp), pointer :: weak_lowT_rate(:)=>NULL() ! (rates_reaction_id_max)
136 : ! these are from reaclib or weak_info.list
137 : ! set at initialization; read-only afterwards.
138 :
139 : integer, pointer :: reaction_screening_info(:,:)=>NULL() !(3,rates_reaction_id_max)
140 : ! reaction_screen_info(1:2,i) = [chem_id1, chem_id2] for screening. 0's if no screening.
141 :
142 : integer, pointer :: weak_reaction_info(:,:)=>NULL() ! (2,rates_reaction_id_max)
143 : ! weak_reaction_info(1:2,i) = [chem_id_in, chem_id_out]. 0's if not a weak reaction.
144 :
145 : integer, pointer :: reaction_ye_rho_exponents(:,:)=>NULL() ! (2,rates_reaction_id_max)
146 : ! multiply T dependent rate by Ye^a(i) * Rho^b(i)
147 : ! reaction_ye_rho_coeffs(1,i) is a(i)
148 : ! reaction_ye_rho_coeffs(2,i) is b(i)
149 : ! (0,0) for photodisintegrations and decays
150 : ! (0,1) for standard 2 body reactions
151 : ! (0,2) for 3 body reactions such as triple alpha
152 : ! (1,1) for 2 body electron captures
153 : ! (1,2) for 3 body electron captures (e.g., pep)
154 :
155 : integer, parameter :: max_num_reaction_inputs = 3
156 : integer, pointer :: reaction_inputs(:,:)=>NULL() ! (2*max_num_reaction_inputs,rates_reaction_id_max)
157 : ! up to max_num_reaction_inputs pairs of coefficients and chem id's, terminated by 0's.
158 : ! e.g., o16(p,g)f17 would be (/ 1, io16, 1, ih1, 0 /)
159 : ! triple alpha would be (/ 3, ihe4, 0 /)
160 : ! he3(he4, g)be7(e-,nu)li7(p,a)he4 would be (/ 1, ihe3, 1, ihe4, i, ih1, 0 /)
161 :
162 : integer, parameter :: max_num_reaction_outputs = 4
163 : integer, pointer :: reaction_outputs(:,:)=>NULL() ! (2*max_num_reaction_outputs,rates_reaction_id_max)
164 : ! up to max_num_reaction_outputs pairs of coefficients and chem id's, terminated by 0's.
165 : ! e.g., o16(p,g)f17 would be (/ 1, if17, 0 /)
166 : ! c12(a, p)n15 would be (/ 1, in15, 1, ih1, 0 /)
167 :
168 : ! weak_info_list
169 :
170 : integer :: num_weak_info_list_reactions
171 : real(dp), pointer :: weak_info_list_halflife(:)=>NULL(), weak_info_list_Qneu(:)=>NULL()
172 : type (integer_dict), pointer :: weak_info_list_dict=>NULL()
173 :
174 : ! weak
175 :
176 : real(dp) :: &
177 : T9_weaklib_full_off = 0.01d0, & ! use pure reaclib for T <= this
178 : T9_weaklib_full_on = 0.02d0 ! use pure weaklib for T >= this
179 : ! blend for intermediate temperatures
180 :
181 : ! for high Z elements, switch to reaclib at temp where no longer fully ionized
182 : ! as rough approximation for this, we switch at Fe to higher values of T9
183 : integer :: weaklib_blend_hi_Z = 26
184 : ! if input element has Z >= weaklib_blend_hi_Z, then use the following T9 limits
185 : real(dp) :: &
186 : T9_weaklib_full_off_hi_Z = 0.063d0, & ! use pure reaclib for T <= this
187 : T9_weaklib_full_on_hi_Z = 0.073d0 ! use pure weaklib for T >= this
188 :
189 : integer :: num_weak_reactions
190 :
191 : type, abstract :: weak_rate_table
192 : integer :: num_T9
193 : integer :: num_lYeRho
194 :
195 : real(dp), allocatable :: T9s(:)
196 : real(dp), allocatable :: lYeRhos(:)
197 :
198 : real(dp), allocatable :: data(:,:,:,:) ! (4, num_T9, num_lYeRho, 3)
199 : contains
200 :
201 : procedure(setup_weak_table), deferred :: setup
202 : procedure(interpolate_weak_table), deferred :: interpolate
203 :
204 : end type weak_rate_table
205 :
206 : abstract interface
207 : subroutine setup_weak_table(table, ierr)
208 : import weak_rate_table
209 : implicit none
210 : class(weak_rate_table), intent(inout) :: table
211 : integer, intent(out) :: ierr
212 : end subroutine setup_weak_table
213 : end interface
214 :
215 : abstract interface
216 : subroutine interpolate_weak_table(table, T9, lYeRho, &
217 : lambda, dlambda_dlnT, dlambda_dlnRho, &
218 : Qneu, dQneu_dlnT, dQneu_dlnRho, ierr)
219 : use const_def, only : dp
220 : import weak_rate_table
221 : implicit none
222 : class(weak_rate_table), intent(inout) :: table
223 : real(dp), intent(in) :: T9, lYeRho
224 : real(dp), intent(out) :: lambda, dlambda_dlnT, dlambda_dlnRho
225 : real(dp), intent(out) :: Qneu, dQneu_dlnT, dQneu_dlnRho
226 : integer, intent(out) :: ierr
227 : end subroutine interpolate_weak_table
228 : end interface
229 :
230 : type :: table_c
231 : class(weak_rate_table), pointer :: t=>NULL()
232 : end type table_c
233 :
234 : type(table_c), dimension(:), allocatable :: weak_reactions_tables
235 :
236 : integer, pointer, dimension(:) :: & ! (num_weak_reactions)
237 : weak_lhs_nuclide_id=>NULL(), weak_rhs_nuclide_id=>NULL(), weak_reaclib_id=>NULL()
238 : character(len=iso_name_length), dimension(:), pointer :: &
239 : weak_lhs_nuclide_name=>NULL(), weak_rhs_nuclide_name=>NULL() ! (num_weak_reactions)
240 : type (integer_dict), pointer :: weak_reactions_dict=>NULL()
241 :
242 : logical :: weak_bicubic = .false.
243 : ! true means do bicubic splines for interpolation
244 : ! false means just do bilinear
245 : ! bilinear is safe; bicubic can overshoot near jumps
246 :
247 : ! Suzuki et al. (2016)
248 : logical :: use_suzuki_tables = .false.
249 :
250 : ! ecapture
251 :
252 : logical :: do_ecapture = .false.
253 :
254 : type (integer_dict), pointer :: ecapture_states_number_dict=>NULL()
255 : type (integer_dict), pointer :: ecapture_states_offset_dict=>NULL()
256 :
257 : type (integer_dict), pointer :: ecapture_transitions_number_dict=>NULL()
258 : type (integer_dict), pointer :: ecapture_transitions_offset_dict=>NULL()
259 :
260 : integer, parameter :: max_num_ecapture_nuclei = 25
261 : integer, parameter :: max_num_ecapture_states = 25
262 :
263 : integer, parameter :: max_num_ecapture_reactions = 25
264 : integer, parameter :: max_num_ecapture_transitions = 25
265 :
266 : integer, parameter :: num_transitions_data = 2
267 : integer, parameter :: i_Si = 1, i_Sf = 2
268 : integer :: num_ecapture_transitions
269 : integer, pointer :: ecapture_transitions_data(:,:)=>NULL()
270 : real(dp), pointer :: ecapture_logft_data(:)=>NULL()
271 :
272 : integer, parameter :: num_states_data = 2
273 : integer, parameter :: i_E = 1, i_J = 2
274 : integer :: num_ecapture_states
275 : real(dp), pointer :: ecapture_states_data(:,:)=>NULL()
276 :
277 : integer :: num_ecapture_nuclei, num_ecapture_reactions
278 : integer, pointer, dimension(:) :: ecapture_nuclide_id=>NULL(), &
279 : ecapture_lhs_nuclide_id=>NULL(), ecapture_rhs_nuclide_id=>NULL() ! (num_ecapture_reactions)
280 : character(len=iso_name_length), dimension(:), pointer :: ecapture_nuclide_name=>NULL(), &
281 : ecapture_lhs_nuclide_name=>NULL(), ecapture_rhs_nuclide_name=>NULL() ! (num_ecapture_reactions)
282 : type (integer_dict), pointer :: ecapture_reactions_dict=>NULL()
283 :
284 :
285 : integer, pointer :: reaction_categories(:)=>NULL() ! (rates_reaction_id_max) set by net using reactions.list info
286 :
287 :
288 : integer, pointer,dimension(:) :: &
289 : reaction_is_reverse, reaction_reaclib_lo, reaction_reaclib_hi, reverse_reaction_id
290 : ! caches for get_reaclib_rate_and_dlnT (in raw_rates.f)
291 : ! for all of these, 0 means "cache entry not yet set -- don't have the information"
292 :
293 :
294 : ! for tabular evaluation of the raw reaction rates
295 : real(dp) :: rattab_thi != 10.301029995664d0 ! log10(highest temp = 2e10)
296 : real(dp) :: rattab_tlo != 5.30102999566398d0 ! log10(lowest temp = 2e5)
297 : real(dp) :: rattab_temp_hi != 10**rattab_thi
298 : real(dp) :: rattab_temp_lo != 10**rattab_tlo
299 :
300 : integer :: rattab_points_per_decade = 2000
301 : integer :: nrattab ! number of reaction rate table temperatures
302 : ! nrattab = <points per decade>*(rattab_thi - rattab_tlo) + 1
303 :
304 : real(dp) :: rattab_tstp != (rattab_thi-rattab_tlo)/(nrattab-1)! step size
305 :
306 : ! reactions for hardwired nets and reactions with multiple choices for rates
307 : integer, parameter :: ir1212 = 1
308 : integer, parameter :: ir1216 = ir1212+1
309 : integer, parameter :: ir1216_to_mg24 = ir1216+1
310 : integer, parameter :: ir1216_to_si28 = ir1216_to_mg24+1
311 : integer, parameter :: ir1616 = ir1216_to_si28+1
312 : integer, parameter :: ir1616a = ir1616+1
313 : integer, parameter :: ir1616g = ir1616a+1
314 : integer, parameter :: ir1616p_aux = ir1616g+1
315 : integer, parameter :: ir1616ppa = ir1616p_aux+1
316 : integer, parameter :: ir1616ppg = ir1616ppa+1
317 : integer, parameter :: ir_he3_ag_be7 = ir1616ppg+1
318 : integer, parameter :: ir34_pp2 = ir_he3_ag_be7+1
319 : integer, parameter :: ir34_pp3 = ir34_pp2+1
320 : integer, parameter :: ir_al27_pa_mg24 = ir34_pp3+1
321 : integer, parameter :: ir_ar36_ag_ca40 = ir_al27_pa_mg24+1
322 : integer, parameter :: ir_ar36_ga_s32 = ir_ar36_ag_ca40+1
323 : integer, parameter :: ir_b8_gp_be7 = ir_ar36_ga_s32+1
324 : integer, parameter :: ir_be7_pg_b8 = ir_b8_gp_be7+1
325 : integer, parameter :: ir_c12_ag_o16 = ir_be7_pg_b8+1
326 : integer, parameter :: ir_c12_ap_n15 = ir_c12_ag_o16+1
327 : integer, parameter :: ir_c12_pg_n13 = ir_c12_ap_n15+1
328 : integer, parameter :: ir_c12_to_he4_he4_he4 = ir_c12_pg_n13+1
329 : integer, parameter :: ir_c13_an_o16 = ir_c12_to_he4_he4_he4+1
330 : integer, parameter :: ir_c13_pg_n14 = ir_c13_an_o16+1
331 : integer, parameter :: ir_ca40_ag_ti44 = ir_c13_pg_n14+1
332 : integer, parameter :: ir_ca40_ga_ar36 = ir_ca40_ag_ti44+1
333 : integer, parameter :: ir_cr48_ag_fe52 = ir_ca40_ga_ar36+1
334 : integer, parameter :: ir_cr48_ga_ti44 = ir_cr48_ag_fe52+1
335 : integer, parameter :: ir_f17_ap_ne20 = ir_cr48_ga_ti44+1
336 : integer, parameter :: ir_f17_gp_o16 = ir_f17_ap_ne20+1
337 : integer, parameter :: ir_f17_pa_o14 = ir_f17_gp_o16+1
338 : integer, parameter :: ir_f18_gp_o17 = ir_f17_pa_o14+1
339 : integer, parameter :: ir_f18_pa_o15 = ir_f18_gp_o17+1
340 :
341 : integer, parameter :: ir_f17_pg_ne18 = ir_f18_pa_o15+1
342 : integer, parameter :: ir_f18_pg_ne19 = ir_f17_pg_ne18+1
343 :
344 : integer, parameter :: ir_f19_ap_ne22 = ir_f18_pg_ne19+1
345 : integer, parameter :: ir_f19_gp_o18 = ir_f19_ap_ne22+1
346 : integer, parameter :: ir_f19_pa_o16 = ir_f19_gp_o18+1
347 : integer, parameter :: ir_f19_pg_ne20 = ir_f19_pa_o16+1
348 : integer, parameter :: ir_fe52_ag_ni56 = ir_f19_pg_ne20+1
349 : integer, parameter :: ir_fe52_ga_cr48 = ir_fe52_ag_ni56+1
350 : integer, parameter :: ir_h2_be7_to_h1_he4_he4 = ir_fe52_ga_cr48+1
351 : integer, parameter :: ir_h2_h2_to_he4 = ir_h2_be7_to_h1_he4_he4+1
352 : integer, parameter :: ir_h2_he3_to_h1_he4 = ir_h2_h2_to_he4+1
353 :
354 : integer, parameter :: ir_he3_be7_to_h1_h1_he4_he4 = ir_h2_he3_to_h1_he4+1
355 : integer, parameter :: ir_he3_he3_to_h1_h1_he4 = ir_he3_be7_to_h1_h1_he4_he4+1
356 : integer, parameter :: ir_h1_h1_he4_to_he3_he3 = ir_he3_he3_to_h1_h1_he4+1
357 : integer, parameter :: ir_he4_he4_he4_to_c12 = ir_h1_h1_he4_to_he3_he3+1
358 : integer, parameter :: ir_li7_pa_he4 = ir_he4_he4_he4_to_c12+1
359 : integer, parameter :: ir_he4_ap_li7 = ir_li7_pa_he4+1
360 : integer, parameter :: ir_mg24_ag_si28 = ir_he4_ap_li7+1
361 : integer, parameter :: ir_mg24_ap_al27 = ir_mg24_ag_si28+1
362 : integer, parameter :: ir_mg24_ga_ne20 = ir_mg24_ap_al27+1
363 : integer, parameter :: ir_n13_ap_o16 = ir_mg24_ga_ne20+1
364 : integer, parameter :: ir_n13_gp_c12 = ir_n13_ap_o16+1
365 : integer, parameter :: ir_n13_pg_o14 = ir_n13_gp_c12+1
366 : integer, parameter :: ir_n14_ag_f18 = ir_n13_pg_o14+1
367 : integer, parameter :: ir_n14_ap_o17 = ir_n14_ag_f18+1
368 : integer, parameter :: ir_n14_gp_c13 = ir_n14_ap_o17+1
369 : integer, parameter :: ir_n14_pg_o15 = ir_n14_gp_c13+1
370 : integer, parameter :: ir_n15_ag_f19 = ir_n14_pg_o15+1
371 : integer, parameter :: ir_n15_ap_o18 = ir_n15_ag_f19+1
372 : integer, parameter :: ir_n15_pa_c12 = ir_n15_ap_o18+1
373 : integer, parameter :: ir_n15_pg_o16 = ir_n15_pa_c12+1
374 : integer, parameter :: ir_na23_pa_ne20 = ir_n15_pg_o16+1
375 : integer, parameter :: ir_ne18_gp_f17 = ir_na23_pa_ne20+1
376 : integer, parameter :: ir_ne19_ga_o15 = ir_ne18_gp_f17+1
377 : integer, parameter :: ir_ne19_gp_f18 = ir_ne19_ga_o15+1
378 : integer, parameter :: ir_ne20_ag_mg24 = ir_ne19_gp_f18+1
379 : integer, parameter :: ir_ne20_ap_na23 = ir_ne20_ag_mg24+1
380 : integer, parameter :: ir_ne20_ga_o16 = ir_ne20_ap_na23+1
381 : integer, parameter :: ir_ne20_gp_f19 = ir_ne20_ga_o16+1
382 : integer, parameter :: ir_ne22_ag_mg26 = ir_ne20_gp_f19+1
383 : integer, parameter :: ir_ne22_pg_na23 = ir_ne22_ag_mg26+1
384 : integer, parameter :: ir_ni56_ga_fe52 = ir_ne22_pg_na23+1
385 : integer, parameter :: ir_o14_ag_ne18 = ir_ni56_ga_fe52+1
386 : integer, parameter :: ir_o14_ap_f17 = ir_o14_ag_ne18+1
387 : integer, parameter :: ir_o14_gp_n13 = ir_o14_ap_f17+1
388 : integer, parameter :: ir_o15_ag_ne19 = ir_o14_gp_n13+1
389 : integer, parameter :: ir_o15_ap_f18 = ir_o15_ag_ne19+1
390 : integer, parameter :: ir_o15_gp_n14 = ir_o15_ap_f18+1
391 : integer, parameter :: ir_o16_ag_ne20 = ir_o15_gp_n14+1
392 : integer, parameter :: ir_o16_ap_f19 = ir_o16_ag_ne20+1
393 : integer, parameter :: ir_o16_ga_c12 = ir_o16_ap_f19+1
394 : integer, parameter :: ir_o16_gp_n15 = ir_o16_ga_c12+1
395 : integer, parameter :: ir_o16_pg_f17 = ir_o16_gp_n15+1
396 : integer, parameter :: ir_o17_pa_n14 = ir_o16_pg_f17+1
397 : integer, parameter :: ir_o17_pg_f18 = ir_o17_pa_n14+1
398 : integer, parameter :: ir_o18_ag_ne22 = ir_o17_pg_f18+1
399 : integer, parameter :: ir_o18_pa_n15 = ir_o18_ag_ne22+1
400 : integer, parameter :: ir_o18_pg_f19 = ir_o18_pa_n15+1
401 : integer, parameter :: ir_s32_ag_ar36 = ir_o18_pg_f19+1
402 : integer, parameter :: ir_s32_ga_si28 = ir_s32_ag_ar36+1
403 : integer, parameter :: ir_si28_ag_s32 = ir_s32_ga_si28+1
404 : integer, parameter :: ir_si28_ga_mg24 = ir_si28_ag_s32+1
405 : integer, parameter :: ir_ti44_ag_cr48 = ir_si28_ga_mg24+1
406 : integer, parameter :: ir_ti44_ga_ca40 = ir_ti44_ag_cr48+1
407 : integer, parameter :: iral27pa_aux = ir_ti44_ga_ca40+1
408 : integer, parameter :: iral27pg_aux = iral27pa_aux+1
409 : integer, parameter :: irar36ap_aux = iral27pg_aux+1
410 : integer, parameter :: irar36ap_to_ca40 = irar36ap_aux+1
411 : integer, parameter :: irar36gp_aux = irar36ap_to_ca40+1
412 : integer, parameter :: irar36gp_to_s32 = irar36gp_aux+1
413 : integer, parameter :: irbe7ec_li7_aux = irar36gp_to_s32+1
414 : integer, parameter :: irbe7pg_b8_aux = irbe7ec_li7_aux+1
415 : integer, parameter :: irc12_to_c13 = irbe7pg_b8_aux+1
416 : integer, parameter :: irc12_to_n14 = irc12_to_c13+1
417 : integer, parameter :: irc12ap_aux = irc12_to_n14+1
418 : integer, parameter :: irc12ap_to_o16 = irc12ap_aux+1
419 : integer, parameter :: irca40ap_aux = irc12ap_to_o16+1
420 : integer, parameter :: irca40ap_to_ti44 = irca40ap_aux+1
421 : integer, parameter :: irca40gp_aux = irca40ap_to_ti44+1
422 : integer, parameter :: irca40gp_to_ar36 = irca40gp_aux+1
423 : integer, parameter :: ircl35pa_aux = irca40gp_to_ar36+1
424 : integer, parameter :: ircl35pg_aux = ircl35pa_aux+1
425 : integer, parameter :: irco55gprot_aux = ircl35pg_aux+1
426 : integer, parameter :: irco55pg_aux = irco55gprot_aux+1
427 : integer, parameter :: irco55protg_aux = irco55pg_aux+1
428 : integer, parameter :: ircr48ap_aux = irco55protg_aux+1
429 : integer, parameter :: ircr48ap_to_fe52 = ircr48ap_aux+1
430 : integer, parameter :: ircr48gp_aux = ircr48ap_to_fe52+1
431 : integer, parameter :: ircr48gp_to_ti44 = ircr48gp_aux+1
432 : integer, parameter :: irf19pg_aux = ircr48gp_to_ti44+1
433 : integer, parameter :: irfe52ap_aux = irf19pg_aux+1
434 : integer, parameter :: irfe52ap_to_ni56 = irfe52ap_aux+1
435 : integer, parameter :: irfe52aprot_aux = irfe52ap_to_ni56+1
436 : integer, parameter :: irfe52aprot_to_fe54 = irfe52aprot_aux+1
437 : integer, parameter :: irfe52aprot_to_ni56 = irfe52aprot_to_fe54+1
438 : integer, parameter :: irfe52gp_aux = irfe52aprot_to_ni56+1
439 : integer, parameter :: irfe52gp_to_cr48 = irfe52gp_aux+1
440 : integer, parameter :: irfe52neut_to_fe54 = irfe52gp_to_cr48+1
441 : integer, parameter :: irfe52ng_aux = irfe52neut_to_fe54+1
442 : integer, parameter :: irfe53gn_aux = irfe52ng_aux+1
443 : integer, parameter :: irfe53ng_aux = irfe53gn_aux+1
444 : integer, parameter :: irfe54a_to_ni56 = irfe53ng_aux+1
445 : integer, parameter :: irfe54an_aux = irfe54a_to_ni56+1
446 : integer, parameter :: irfe54an_to_ni56 = irfe54an_aux+1
447 : integer, parameter :: irfe54aprot_to_fe56 = irfe54an_to_ni56+1
448 : integer, parameter :: irfe54g_to_fe52 = irfe54aprot_to_fe56+1
449 : integer, parameter :: irfe54ng_aux = irfe54g_to_fe52+1
450 : integer, parameter :: irfe54ng_to_fe56 = irfe54ng_aux+1
451 : integer, parameter :: irfe54prot_to_fe52 = irfe54ng_to_fe56+1
452 : integer, parameter :: irfe54prot_to_ni56 = irfe54prot_to_fe52+1
453 : integer, parameter :: irfe54protg_aux = irfe54prot_to_ni56+1
454 : integer, parameter :: irfe55gn_aux = irfe54protg_aux+1
455 : integer, parameter :: irfe55ng_aux = irfe55gn_aux+1
456 :
457 : integer, parameter :: irfe56ec_fake_to_mn56 = irfe55ng_aux+1
458 : integer, parameter :: irfe56ec_fake_to_mn57 = irfe56ec_fake_to_mn56+1
459 : integer, parameter :: irfe56ec_fake_to_cr56 = irfe56ec_fake_to_mn57+1
460 : integer, parameter :: irfe56ec_fake_to_cr57 = irfe56ec_fake_to_cr56+1
461 : integer, parameter :: irfe56ec_fake_to_cr58 = irfe56ec_fake_to_cr57+1
462 : integer, parameter :: irfe56ec_fake_to_cr59 = irfe56ec_fake_to_cr58+1
463 : integer, parameter :: irfe56ec_fake_to_cr60 = irfe56ec_fake_to_cr59+1
464 : integer, parameter :: irfe56ec_fake_to_cr61 = irfe56ec_fake_to_cr60+1
465 : integer, parameter :: irfe56ec_fake_to_cr62 = irfe56ec_fake_to_cr61+1
466 : integer, parameter :: irfe56ec_fake_to_cr63 = irfe56ec_fake_to_cr62+1
467 : integer, parameter :: irfe56ec_fake_to_cr64 = irfe56ec_fake_to_cr63+1
468 : integer, parameter :: irfe56ec_fake_to_cr65 = irfe56ec_fake_to_cr64+1
469 : integer, parameter :: irfe56ec_fake_to_cr66 = irfe56ec_fake_to_cr65+1
470 :
471 : integer, parameter :: irfe56ee_to_ni56 = irfe56ec_fake_to_cr66+1
472 : integer, parameter :: irfe56gn_aux = irfe56ee_to_ni56+1
473 : integer, parameter :: irfe56gn_to_fe54 = irfe56gn_aux+1
474 : integer, parameter :: irfe56prot_to_fe54 = irfe56gn_to_fe54+1
475 : integer, parameter :: irh2_protg_aux = irfe56prot_to_fe54+1
476 : integer, parameter :: irh2g_neut_aux = irh2_protg_aux+1
477 : integer, parameter :: irhe3_neutg_aux = irh2g_neut_aux+1
478 : integer, parameter :: irhe3gprot_aux = irhe3_neutg_aux+1
479 : integer, parameter :: irhe4_breakup = irhe3gprot_aux+1
480 : integer, parameter :: irhe4_rebuild = irhe4_breakup+1
481 : integer, parameter :: irhe4g_neut_aux = irhe4_rebuild+1
482 : integer, parameter :: irk39pa_aux = irhe4g_neut_aux+1
483 : integer, parameter :: irk39pg_aux = irk39pa_aux+1
484 : integer, parameter :: irmg24ap_aux = irk39pg_aux+1
485 : integer, parameter :: irmg24ap_to_si28 = irmg24ap_aux+1
486 : integer, parameter :: irmg24gp_aux = irmg24ap_to_si28+1
487 : integer, parameter :: irmg24gp_to_ne20 = irmg24gp_aux+1
488 : integer, parameter :: irmn51pg_aux = irmg24gp_to_ne20+1
489 : integer, parameter :: irn14_to_c12 = irmn51pg_aux+1
490 : integer, parameter :: irn14_to_n15 = irn14_to_c12+1
491 : integer, parameter :: irn14_to_o16 = irn14_to_n15+1
492 : integer, parameter :: irn14ag_lite = irn14_to_o16+1
493 : integer, parameter :: irn14gc12 = irn14ag_lite+1
494 : integer, parameter :: irn14pg_aux = irn14gc12+1
495 : integer, parameter :: irn15pa_aux = irn14pg_aux+1
496 : integer, parameter :: irn15pg_aux = irn15pa_aux+1
497 : integer, parameter :: irna23pa_aux = irn15pg_aux+1
498 : integer, parameter :: irna23pg_aux = irna23pa_aux+1
499 : integer, parameter :: irne18ag_to_mg24 = irna23pg_aux+1
500 : integer, parameter :: irne18ap_to_mg22 = irne18ag_to_mg24+1
501 : integer, parameter :: irne18ap_to_mg24 = irne18ap_to_mg22+1
502 : integer, parameter :: irne19pg_to_mg22 = irne18ap_to_mg24+1
503 : integer, parameter :: irne19pg_to_mg24 = irne19pg_to_mg22+1
504 : integer, parameter :: irne20ap_aux = irne19pg_to_mg24+1
505 : integer, parameter :: irne20ap_to_mg24 = irne20ap_aux+1
506 : integer, parameter :: irne20gp_aux = irne20ap_to_mg24+1
507 : integer, parameter :: irne20gp_to_o16 = irne20gp_aux+1
508 : integer, parameter :: irne20pg_to_mg22 = irne20gp_to_o16+1
509 : integer, parameter :: irne20pg_to_mg24 = irne20pg_to_mg22+1
510 : integer, parameter :: irneut_to_prot = irne20pg_to_mg24+1
511 : integer, parameter :: irni56ec_to_fe54 = irneut_to_prot+1
512 : integer, parameter :: irni56ec_to_fe56 = irni56ec_to_fe54+1
513 : integer, parameter :: irni56ec_to_co56 = irni56ec_to_fe56+1
514 : integer, parameter :: irco56ec_to_fe56 = irni56ec_to_co56+1
515 : integer, parameter :: irni56gp_aux = irco56ec_to_fe56+1
516 : integer, parameter :: irni56gp_to_fe52 = irni56gp_aux+1
517 : integer, parameter :: irni56gprot_aux = irni56gp_to_fe52+1
518 : integer, parameter :: irni56gprot_to_fe52 = irni56gprot_aux+1
519 : integer, parameter :: irni56gprot_to_fe54 = irni56gprot_to_fe52+1
520 : integer, parameter :: irni56ng_to_fe54 = irni56gprot_to_fe54+1
521 : integer, parameter :: irni57na_aux = irni56ng_to_fe54+1
522 : integer, parameter :: iro16_to_n14 = irni57na_aux+1
523 : integer, parameter :: iro16_to_o17 = iro16_to_n14+1
524 : integer, parameter :: iro16ap_aux = iro16_to_o17+1
525 : integer, parameter :: iro16ap_to_ne20 = iro16ap_aux+1
526 : integer, parameter :: iro16gp_aux = iro16ap_to_ne20+1
527 : integer, parameter :: iro16gp_to_c12 = iro16gp_aux+1
528 : integer, parameter :: iro17_to_o18 = iro16gp_to_c12+1
529 : integer, parameter :: irp31pa_aux = iro17_to_o18+1
530 : integer, parameter :: irp31pg_aux = irp31pa_aux+1
531 : integer, parameter :: irpep_to_he3 = irp31pg_aux+1
532 : integer, parameter :: irpp_to_he3 = irpep_to_he3+1
533 : integer, parameter :: irprot_neutg_aux = irpp_to_he3+1
534 : integer, parameter :: irprot_to_neut = irprot_neutg_aux+1
535 : integer, parameter :: irs32ap_aux = irprot_to_neut+1
536 : integer, parameter :: irs32ap_to_ar36 = irs32ap_aux+1
537 : integer, parameter :: irs32gp_aux = irs32ap_to_ar36+1
538 : integer, parameter :: irs32gp_to_si28 = irs32gp_aux+1
539 : integer, parameter :: irsc43pa_aux = irs32gp_to_si28+1
540 : integer, parameter :: irsc43pg_aux = irsc43pa_aux+1
541 : integer, parameter :: irsi28ap_aux = irsc43pg_aux+1
542 : integer, parameter :: irsi28ap_to_s32 = irsi28ap_aux+1
543 : integer, parameter :: irsi28gp_aux = irsi28ap_to_s32+1
544 : integer, parameter :: irsi28gp_to_mg24 = irsi28gp_aux+1
545 : integer, parameter :: irti44ap_aux = irsi28gp_to_mg24+1
546 : integer, parameter :: irti44ap_to_cr48 = irti44ap_aux+1
547 : integer, parameter :: irti44gp_aux = irti44ap_to_cr48+1
548 : integer, parameter :: irti44gp_to_ca40 = irti44gp_aux+1
549 : integer, parameter :: irv47pa_aux = irti44gp_to_ca40+1
550 : integer, parameter :: irv47pg_aux = irv47pa_aux+1
551 : integer, parameter :: ir_h1_h1_wk_h2 = irv47pg_aux+1
552 : integer, parameter :: ir_h1_h1_ec_h2 = ir_h1_h1_wk_h2+1
553 : integer, parameter :: irn14ag_to_ne22 = ir_h1_h1_ec_h2+1
554 : integer, parameter :: irf19pa_aux = irn14ag_to_ne22+1
555 : integer, parameter :: ir_b8_wk_he4_he4 = irf19pa_aux+1
556 : integer, parameter :: irmn51pa_aux = ir_b8_wk_he4_he4+1
557 : integer, parameter :: irfe54gn_aux = irmn51pa_aux+1
558 : integer, parameter :: irco55pa_aux = irfe54gn_aux+1
559 : integer, parameter :: irco55prota_aux = irco55pa_aux+1
560 : integer, parameter :: irn14ag_to_o18 = irco55prota_aux+1
561 : integer, parameter :: ir_h1_he3_wk_he4 = irn14ag_to_o18+1
562 : integer, parameter :: ir_be7_wk_li7 = ir_h1_he3_wk_he4+1
563 :
564 : ! rates added for approx21
565 : integer, parameter :: ir_al27_pg_si28 = ir_be7_wk_li7+1
566 : integer, parameter :: ir_si28_gp_al27 = ir_al27_pg_si28+1
567 : integer, parameter :: ir_si28_ap_p31 = ir_si28_gp_al27+1
568 : integer, parameter :: ir_p31_pa_si28 = ir_si28_ap_p31+1
569 : integer, parameter :: ir_p31_pg_s32 = ir_p31_pa_si28+1
570 : integer, parameter :: ir_s32_gp_p31 = ir_p31_pg_s32+1
571 : integer, parameter :: ir_s32_ap_cl35 = ir_s32_gp_p31+1
572 : integer, parameter :: ir_cl35_pa_s32 = ir_s32_ap_cl35+1
573 : integer, parameter :: ir_cl35_pg_ar36 = ir_cl35_pa_s32+1
574 : integer, parameter :: ir_ar36_gp_cl35 = ir_cl35_pg_ar36+1
575 : integer, parameter :: ir_ar36_ap_k39 = ir_ar36_gp_cl35+1
576 : integer, parameter :: ir_k39_pa_ar36 = ir_ar36_ap_k39+1
577 : integer, parameter :: ir_k39_pg_ca40 = ir_k39_pa_ar36+1
578 : integer, parameter :: ir_ca40_gp_k39 = ir_k39_pg_ca40+1
579 : integer, parameter :: ir_ca40_ap_sc43 = ir_ca40_gp_k39+1
580 : integer, parameter :: ir_sc43_pa_ca40 = ir_ca40_ap_sc43+1
581 : integer, parameter :: ir_sc43_pg_ti44 = ir_sc43_pa_ca40+1
582 : integer, parameter :: ir_ti44_gp_sc43 = ir_sc43_pg_ti44+1
583 : integer, parameter :: ir_ti44_ap_v47 = ir_ti44_gp_sc43+1
584 : integer, parameter :: ir_v47_pa_ti44 = ir_ti44_ap_v47+1
585 : integer, parameter :: ir_v47_pg_cr48 = ir_v47_pa_ti44+1
586 : integer, parameter :: ir_cr48_gp_v47 = ir_v47_pg_cr48+1
587 : integer, parameter :: ir_cr48_ap_mn51 = ir_cr48_gp_v47+1
588 : integer, parameter :: ir_mn51_pa_cr48 = ir_cr48_ap_mn51+1
589 : integer, parameter :: ir_mn51_pg_fe52 = ir_mn51_pa_cr48+1
590 : integer, parameter :: ir_fe52_gp_mn51 = ir_mn51_pg_fe52+1
591 : integer, parameter :: ir_fe52_ap_co55 = ir_fe52_gp_mn51+1
592 : integer, parameter :: ir_co55_pa_fe52 = ir_fe52_ap_co55+1
593 : integer, parameter :: ir_co55_pg_ni56 = ir_co55_pa_fe52+1
594 : integer, parameter :: ir_ni56_gp_co55 = ir_co55_pg_ni56+1
595 : integer, parameter :: ir_fe52_ng_fe53 = ir_ni56_gp_co55+1
596 : integer, parameter :: ir_fe53_gn_fe52 = ir_fe52_ng_fe53+1
597 : integer, parameter :: ir_fe53_ng_fe54 = ir_fe53_gn_fe52+1
598 : integer, parameter :: ir_fe54_gn_fe53 = ir_fe53_ng_fe54+1
599 : integer, parameter :: ir_fe54_pg_co55 = ir_fe54_gn_fe53+1
600 : integer, parameter :: ir_co55_gp_fe54 = ir_fe54_pg_co55+1
601 : integer, parameter :: ir_he3_ng_he4 = ir_co55_gp_fe54+1
602 : integer, parameter :: ir_he4_gn_he3 = ir_he3_ng_he4+1
603 : integer, parameter :: ir_h1_ng_h2 = ir_he4_gn_he3+1
604 : integer, parameter :: ir_h2_gn_h1 = ir_h1_ng_h2+1
605 : integer, parameter :: ir_h2_pg_he3 = ir_h2_gn_h1+1
606 : integer, parameter :: ir_he3_gp_h2 = ir_h2_pg_he3+1
607 : integer, parameter :: ir_fe54_ng_fe55 = ir_he3_gp_h2+1
608 : integer, parameter :: ir_fe55_gn_fe54 = ir_fe54_ng_fe55+1
609 : integer, parameter :: ir_fe55_ng_fe56 = ir_fe55_gn_fe54+1
610 : integer, parameter :: ir_fe56_gn_fe55 = ir_fe55_ng_fe56+1
611 : integer, parameter :: ir_fe54_ap_co57 = ir_fe56_gn_fe55+1
612 : integer, parameter :: ir_co57_pa_fe54 = ir_fe54_ap_co57+1
613 : integer, parameter :: ir_fe56_pg_co57 = ir_co57_pa_fe54+1
614 : integer, parameter :: ir_co57_gp_fe56 = ir_fe56_pg_co57+1
615 :
616 : integer, parameter :: ir_c12_c12_to_h1_na23 = ir_co57_gp_fe56+1
617 : integer, parameter :: ir_he4_ne20_to_c12_c12 = ir_c12_c12_to_h1_na23+1
618 : integer, parameter :: ir_c12_c12_to_he4_ne20 = ir_he4_ne20_to_c12_c12+1
619 : integer, parameter :: ir_he4_mg24_to_c12_o16 = ir_c12_c12_to_he4_ne20+1
620 :
621 :
622 : ! fxt for al26 isomers
623 : integer, parameter :: ir_al26_1_to_al26_2 = ir_he4_mg24_to_c12_o16 + 1
624 : integer, parameter :: ir_al26_2_to_al26_1 = ir_al26_1_to_al26_2 + 1
625 :
626 :
627 : integer, parameter :: num_predefined_reactions = ir_al26_2_to_al26_1
628 :
629 : integer :: rates_reaction_id_max
630 :
631 :
632 : ! for mazurek's ni56 electron capture rate interpolation
633 : real(dp) :: tv(7),rv(6),rfdm(4),rfd0(4),rfd1(4),rfd2(4),tfdm(5),tfd0(5),tfd1(5),tfd2(5)
634 :
635 :
636 : type T_Factors
637 : real(dp) :: lnT9
638 : real(dp) :: T9
639 : real(dp) :: T92
640 : real(dp) :: T93
641 : real(dp) :: T94
642 : real(dp) :: T95
643 : real(dp) :: T96
644 : real(dp) :: T912
645 : real(dp) :: T932
646 : real(dp) :: T952
647 : real(dp) :: T972
648 : real(dp) :: T913
649 : real(dp) :: T923
650 : real(dp) :: T943
651 : real(dp) :: T953
652 : real(dp) :: T973
653 : real(dp) :: T9113
654 : real(dp) :: T914
655 : real(dp) :: T934
656 : real(dp) :: T954
657 : real(dp) :: T974
658 : real(dp) :: T915
659 : real(dp) :: T935
660 : real(dp) :: T945
661 : real(dp) :: T965
662 : real(dp) :: T917
663 : real(dp) :: T927
664 : real(dp) :: T947
665 : real(dp) :: T918
666 : real(dp) :: T938
667 : real(dp) :: T958
668 : real(dp) :: T9i
669 : real(dp) :: T9i2
670 : real(dp) :: T9i3
671 : real(dp) :: T9i12
672 : real(dp) :: T9i32
673 : real(dp) :: T9i52
674 : real(dp) :: T9i72
675 : real(dp) :: T9i13
676 : real(dp) :: T9i23
677 : real(dp) :: T9i43
678 : real(dp) :: T9i53
679 : real(dp) :: T9i14
680 : real(dp) :: T9i34
681 : real(dp) :: T9i54
682 : real(dp) :: T9i15
683 : real(dp) :: T9i35
684 : real(dp) :: T9i45
685 : real(dp) :: T9i65
686 : real(dp) :: T9i17
687 : real(dp) :: T9i27
688 : real(dp) :: T9i47
689 : real(dp) :: T9i18
690 : real(dp) :: T9i38
691 : real(dp) :: T9i58
692 : real(dp) :: T916
693 : real(dp) :: T976
694 : real(dp) :: T9i76
695 : end type T_Factors
696 :
697 :
698 : ! rate results components
699 :
700 : integer, parameter :: i_rate = 1
701 : integer, parameter :: i_rate_dT = 2
702 : integer, parameter :: i_rate_dRho = 3
703 : integer, parameter :: num_rvs = 3
704 :
705 :
706 : ! screening
707 :
708 : integer, parameter :: no_screening = 0
709 : ! 1 was graboske screening so leave undefined so its an error if people keep trying to use it
710 : integer, parameter :: extended_screening = 2
711 : ! based on code from Frank Timmes
712 : ! extends the Graboske method using results from Alastuey and Jancovici (1978),
713 : ! along with plasma parameters from Itoh et al (1979) for strong screening.
714 : integer, parameter :: salpeter_screening = 3
715 : ! weak screening only. following Salpeter (1954),
716 : ! with equations (4-215) and (4-221) of Clayton (1968).
717 : integer, parameter :: chugunov_screening = 4
718 : ! based on code from Sam Jones
719 : ! Implements screening from Chugunov et al (2007)
720 : integer, parameter :: other_screening = 5 ! User defined screening
721 :
722 : type Screen_Info
723 : real(dp) :: temp
724 : real(dp) :: den
725 : real(dp) :: logT
726 : real(dp) :: logRho
727 : real(dp) :: zbar
728 : real(dp) :: abar
729 : real(dp) :: z2bar
730 : real(dp) :: zbar13
731 : real(dp) :: z1pt58bar
732 : real(dp) :: ytot
733 : real(dp) :: rr
734 : real(dp) :: tempi
735 : real(dp) :: dtempi
736 : real(dp) :: deni
737 : real(dp) :: pp
738 : real(dp) :: dppdt
739 : real(dp) :: dppdd
740 : real(dp) :: qlam0z
741 : real(dp) :: qlam0zdt
742 : real(dp) :: qlam0zdd
743 : real(dp) :: taufac
744 : real(dp) :: taufacdt
745 : real(dp) :: xni
746 : real(dp) :: dxnidd
747 : real(dp) :: aa
748 : real(dp) :: daadt
749 : real(dp) :: daadd
750 : real(dp) :: ntot, a_e
751 : integer :: num_calls, num_cache_hits
752 : end type Screen_Info
753 :
754 :
755 : interface
756 : subroutine other_screening_interface(sc, z1, z2, a1, a2, screen, dscreendt, dscreendd, ierr)
757 : import dp, screen_info
758 : implicit none
759 :
760 : type (Screen_Info) :: sc ! See rates_def
761 : real(dp),intent(in) :: z1, z2 !< charge numbers of reactants
762 : real(dp),intent(in) :: a1, a2 !< mass numbers of reactants
763 : real(dp),intent(out) :: screen !< on return, screening factor for this reaction
764 : real(dp),intent(out) :: dscreendt !< on return, temperature derivative of the screening factor
765 : real(dp),intent(out) :: dscreendd !< on return, density derivative of the screening factor
766 : integer, intent(out) :: ierr
767 :
768 : end subroutine other_screening_interface
769 :
770 : subroutine other_rate_get_interface(ir, temp, tf, raw_rate, ierr)
771 : import dp, t_factors
772 : implicit none
773 :
774 : integer :: ir ! Rate id
775 : real(dp),intent(in) :: temp !< Temperature
776 : type (T_Factors) :: tf !< Various temperature factors
777 : real(dp),intent(inout) :: raw_rate !< Unscreened reaction_rate, note this will have the default mesa rate on entry
778 : integer, intent(out) :: ierr
779 :
780 : end subroutine other_rate_get_interface
781 :
782 :
783 : end interface
784 :
785 :
786 : real(dp) :: reaclib_min_T9 ! for T9 < this, return 0 for reaclib strong rates
787 :
788 : type (integer_dict), pointer :: reaction_names_dict
789 :
790 : logical :: have_finished_initialization = .false.
791 : logical :: rates_use_cache = .true.
792 :
793 : procedure (other_screening_interface), pointer :: &
794 : rates_other_screening => null()
795 :
796 : procedure (other_rate_get_interface), pointer :: &
797 : rates_other_rate_get => null()
798 :
799 :
800 : ! choices for various rates
801 : ! NOTE: if change these, must edit raw_rates to match.
802 :
803 : ! NACRE = Angulo et al. 1999 Nucl. Phys. A, 656, 3
804 : ! This is for reactions that care about their values below 10^7K
805 : ! JR = jina reaclib -- (Sakharuk et al. 2006)
806 : ! This is for everything else
807 : ! CF88 = Frank Timmes' version of
808 : ! Caughlin, G. R. & Fowler, W. A. 1988, Atom. Data and Nuc. Data Tables, 40, 283
809 : ! This fills in some of the gaps for rates not in REACLIB or NACRE
810 :
811 : ! Here be dragons, higher temperatures
812 : ! will generate possibly un-physical values as the partition table cuts off at 1d10
813 : ! and the polynomial fits to the rates cuts off at 1d10. We truncate the rates whether
814 : ! the warn flag is on or off, to stop the truncation set a higher max_safe_logT_for_rates
815 :
816 : ! Warn if rates exceed the max usable temperature
817 : logical :: warn_rates_for_high_temp = .true.
818 : real(dp) :: max_safe_logT_for_rates = 10d0
819 :
820 : ! Maximum sensible value for a reaction rate
821 : ! This is to try and catch rates that go bad
822 : real(dp),parameter :: max_safe_rate_for_any_temp = 1d40
823 :
824 : ! info for rates being evaluated using tables (rate_list.txt)
825 : type rate_table_info
826 : logical :: use_rate_table
827 : logical :: need_to_read
828 : character (len=132) :: rate_fname
829 : integer :: nT8s
830 : real(dp), pointer :: T8s(:) ! (nT8s)
831 : real(dp), pointer :: f1(:) ! =(4,nT8s)
832 : end type rate_table_info
833 :
834 : type (rate_table_info), pointer :: raw_rates_records(:)
835 : character (len=1000) :: rates_table_dir
836 :
837 : type (integer_dict), pointer :: skip_warnings_dict
838 :
839 : type (reaction_data), target :: reaclib_rates
840 :
841 : character (len=1000) :: rates_dir, rates_cache_dir, rates_temp_cache_dir
842 :
843 :
844 : ! coulomb corrections for weak reactions
845 : integer :: which_vs_coulomb = 0
846 : integer :: which_mui_coulomb = 0
847 :
848 : type Coulomb_Info
849 : real(dp) :: temp
850 : real(dp) :: den
851 : real(dp) :: logT
852 : real(dp) :: logRho
853 : real(dp) :: zbar
854 : real(dp) :: abar
855 : real(dp) :: z2bar
856 : real(dp) :: ye
857 : type(auto_diff_real_2var_order1) :: rs
858 : type(auto_diff_real_2var_order1) :: gamma_e
859 : end type Coulomb_Info
860 :
861 :
862 : logical :: star_debugging_rates_flag
863 : real(dp) :: rates_test_partials_val, rates_test_partials_dval_dx
864 : real(dp) :: rates_test_partials_logT_lo, rates_test_partials_logT_hi
865 : real(dp) :: rates_test_partials_logRho_lo, rates_test_partials_logRho_hi
866 :
867 :
868 : contains
869 :
870 :
871 5 : subroutine set_rates_cache_dir(rates_cache_dir_in, ierr)
872 : use const_def, only: mesa_data_dir, mesa_caches_dir, mesa_temp_caches_dir, use_mesa_temp_cache
873 : use utils_lib, only : mkdir, switch_str
874 : character (len=*), intent(in) :: rates_cache_dir_in
875 : integer, intent(out) :: ierr
876 5 : ierr = 0
877 5 : rates_dir = trim(mesa_data_dir) // '/rates_data'
878 5 : if (len_trim(rates_cache_dir_in) > 0) then
879 0 : rates_cache_dir = rates_cache_dir_in
880 5 : else if (len_trim(mesa_caches_dir) > 0) then
881 0 : rates_cache_dir = trim(mesa_caches_dir) // '/rates_cache'
882 : else
883 5 : rates_cache_dir = trim(rates_dir) // '/cache'
884 : end if
885 5 : if (rates_use_cache) call mkdir(rates_cache_dir)
886 :
887 5 : rates_temp_cache_dir=trim(mesa_temp_caches_dir)//'/rates_cache'
888 5 : if (use_mesa_temp_cache) call mkdir(rates_temp_cache_dir)
889 :
890 5 : end subroutine set_rates_cache_dir
891 :
892 :
893 5 : subroutine start_rates_def_init(ierr)
894 : use utils_lib, only: integer_dict_define
895 : use math_lib
896 : integer, intent(out) :: ierr
897 :
898 : integer :: i
899 : ierr = 0
900 5 : star_debugging_rates_flag = .false.
901 5 : call create_skip_warnings_dict(ierr)
902 5 : if (ierr /= 0) return
903 5 : nullify(reaction_names_dict)
904 1540 : do i=1,rates_reaction_id_max
905 1535 : call integer_dict_define(reaction_names_dict, reaction_Name(i), i, ierr)
906 1540 : if (ierr /= 0) then
907 0 : write(*,*) 'FATAL ERROR: rates_def_init failed in integer_dict_define'
908 0 : return
909 : end if
910 : end do
911 5 : call do_start_rates_def_init(ierr)
912 :
913 5 : end subroutine start_rates_def_init
914 :
915 :
916 10 : subroutine create_skip_warnings_dict(ierr)
917 5 : use utils_lib
918 : use utils_def
919 : integer, intent(out) :: ierr
920 :
921 : integer :: iounit, n, i, t
922 : character (len=256) :: buffer, string, filename, list_filename
923 :
924 5 : ierr = 0
925 :
926 5 : nullify(skip_warnings_dict)
927 :
928 5 : list_filename = 'skip_warnings.list'
929 : ! first try the local directory
930 5 : filename = trim(list_filename)
931 5 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
932 5 : if (ierr /= 0) then ! if don't find that file, look in rates_data
933 5 : filename = trim(rates_dir) // '/' // trim(list_filename)
934 5 : ierr = 0
935 5 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
936 5 : if (ierr /= 0) then
937 0 : write(*,*) 'failed to open file ' // trim(list_filename)
938 0 : return
939 : end if
940 : end if
941 :
942 5 : n = 0
943 5 : i = 0
944 :
945 : reaction_loop: do
946 575 : t = token(iounit, n, i, buffer, string)
947 575 : if (t == eof_token) exit reaction_loop
948 570 : if (t /= name_token) then
949 0 : call error; return
950 : end if
951 570 : call integer_dict_define(skip_warnings_dict, string, 1, ierr)
952 575 : if (ierr /= 0) then
953 0 : write(*,*) 'FATAL ERROR: create_skip_warnings_dict failed in integer_dict_define'
954 0 : return
955 : end if
956 : end do reaction_loop
957 :
958 5 : close(iounit)
959 :
960 5 : call integer_dict_create_hash(skip_warnings_dict, ierr)
961 5 : if (ierr /= 0) then
962 0 : write(*,*) 'FATAL ERROR: create_skip_warnings_dict failed'
963 0 : return
964 : end if
965 :
966 : contains
967 :
968 0 : subroutine error
969 0 : ierr = -1
970 0 : close(iounit)
971 0 : end subroutine error
972 :
973 : end subroutine create_skip_warnings_dict
974 :
975 :
976 0 : integer function reaclib_index(handle) result(indx)
977 : use utils_lib, only: integer_dict_lookup
978 : character(len=*), intent(in) :: handle ! as in rates% reaction_handle
979 : integer :: ierr
980 : ierr = 0
981 0 : call integer_dict_lookup(reaclib_rates% reaction_dict, handle, indx, ierr)
982 0 : if (ierr /= 0) indx = 0
983 0 : end function reaclib_index
984 :
985 :
986 300 : integer function reaclib_reverse(handle) result(indx)
987 : use utils_lib, only: integer_dict_lookup
988 : character(len=*), intent(in) :: handle ! as in rates% reaction_handle
989 : integer :: ierr
990 : ierr = 0
991 300 : call integer_dict_lookup(reaclib_rates% reverse_dict, handle, indx, ierr)
992 300 : if (ierr /= 0) indx = 0
993 300 : end function reaclib_reverse
994 :
995 :
996 5 : subroutine weaklib_init(ierr)
997 : use const_def, only: mesa_data_dir
998 : integer, intent(out) :: ierr
999 5 : ierr = 0
1000 5 : weak_data_dir = trim(mesa_data_dir) // '/rates_data'
1001 5 : nullify(weak_reactions_dict)
1002 5 : nullify(weak_info_list_dict)
1003 5 : end subroutine weaklib_init
1004 :
1005 :
1006 3 : subroutine free_weak_info
1007 :
1008 : use utils_lib, only: integer_dict_free
1009 :
1010 : integer :: i
1011 :
1012 3 : if (ALLOCATED(weak_reactions_tables)) then
1013 1392 : do i = 1, num_weak_reactions
1014 2781 : if (ASSOCIATED(weak_reactions_tables(i)%t)) deallocate(weak_reactions_tables(i)%t)
1015 : end do
1016 3 : deallocate(weak_reactions_tables)
1017 : end if
1018 :
1019 3 : if (ASSOCIATED(weak_lhs_nuclide_id)) deallocate(weak_lhs_nuclide_id)
1020 3 : if (ASSOCIATED(weak_rhs_nuclide_id)) deallocate(weak_rhs_nuclide_id)
1021 3 : if (ASSOCIATED(weak_lhs_nuclide_name)) deallocate(weak_lhs_nuclide_name)
1022 3 : if (ASSOCIATED(weak_rhs_nuclide_name)) deallocate(weak_rhs_nuclide_name)
1023 3 : if (ASSOCIATED(weak_reaclib_id)) deallocate(weak_reaclib_id)
1024 :
1025 3 : if (ASSOCIATED(weak_info_list_halflife)) deallocate(weak_info_list_halflife)
1026 3 : if (ASSOCIATED(weak_info_list_Qneu)) deallocate(weak_info_list_Qneu)
1027 :
1028 3 : if (ASSOCIATED(weak_reactions_dict)) call integer_dict_free(weak_reactions_dict)
1029 3 : if (ASSOCIATED(weak_info_list_dict)) call integer_dict_free(weak_info_list_dict)
1030 :
1031 3 : end subroutine free_weak_info
1032 :
1033 :
1034 5 : subroutine ecapture_init(ierr)
1035 :
1036 : integer, intent(out) :: ierr
1037 :
1038 5 : ierr = 0
1039 :
1040 5 : nullify(ecapture_reactions_dict)
1041 5 : nullify(ecapture_transitions_number_dict)
1042 5 : nullify(ecapture_transitions_offset_dict)
1043 5 : nullify(ecapture_states_number_dict)
1044 5 : nullify(ecapture_states_offset_dict)
1045 :
1046 5 : end subroutine ecapture_init
1047 :
1048 :
1049 3 : subroutine free_ecapture_info
1050 :
1051 : use utils_lib, only: integer_dict_free
1052 :
1053 3 : if (ASSOCIATED(ecapture_transitions_data)) deallocate(ecapture_transitions_data)
1054 3 : if (ASSOCIATED(ecapture_states_data)) deallocate(ecapture_states_data)
1055 3 : if (ASSOCIATED(ecapture_logft_data)) deallocate(ecapture_logft_data)
1056 3 : if (ASSOCIATED(ecapture_nuclide_id)) deallocate(ecapture_nuclide_id)
1057 3 : if (ASSOCIATED(ecapture_lhs_nuclide_id)) deallocate(ecapture_lhs_nuclide_id)
1058 3 : if (ASSOCIATED(ecapture_rhs_nuclide_id)) deallocate(ecapture_rhs_nuclide_id)
1059 3 : if (ASSOCIATED(ecapture_nuclide_name)) deallocate(ecapture_nuclide_name)
1060 3 : if (ASSOCIATED(ecapture_lhs_nuclide_name)) deallocate(ecapture_lhs_nuclide_name)
1061 3 : if (ASSOCIATED(ecapture_rhs_nuclide_name)) deallocate(ecapture_rhs_nuclide_name)
1062 :
1063 3 : if (ASSOCIATED(ecapture_reactions_dict)) call integer_dict_free(ecapture_reactions_dict)
1064 3 : if (ASSOCIATED(ecapture_transitions_number_dict)) call integer_dict_free(ecapture_transitions_number_dict)
1065 3 : if (ASSOCIATED(ecapture_transitions_offset_dict)) call integer_dict_free(ecapture_transitions_offset_dict)
1066 3 : if (ASSOCIATED(ecapture_states_number_dict)) call integer_dict_free(ecapture_states_number_dict)
1067 3 : if (ASSOCIATED(ecapture_states_offset_dict)) call integer_dict_free(ecapture_states_offset_dict)
1068 :
1069 3 : end subroutine free_ecapture_info
1070 :
1071 :
1072 5 : subroutine reaclib_init(jina_reaclib_filename)
1073 : use const_def, only: mesa_data_dir
1074 : character (len=*), intent(in) :: jina_reaclib_filename
1075 5 : reaclib_dir = trim(mesa_data_dir) // '/rates_data'
1076 : !reaclib_filename = 'jina_reaclib_results05301331'
1077 5 : reaclib_filename = jina_reaclib_filename
1078 5 : if (len_trim(reaclib_filename) == 0) &
1079 5 : reaclib_filename = 'jina_reaclib_results_20171020_default'
1080 5 : end subroutine reaclib_init
1081 :
1082 :
1083 5 : subroutine allocate_reaclib_data(r, n, ierr)
1084 : type(reaclib_data), intent(inout) :: r
1085 : integer, intent(in) :: n
1086 : integer, intent(out) :: ierr
1087 5 : ierr = 0
1088 : allocate( &
1089 : r% chapter(n), r% species(max_species_per_reaction,1:n), &
1090 : r% label(n), r% reaction_flag(n), r% reverse_flag(n), &
1091 10 : r% Qvalue(n), r% coefficients(ncoefficients,1:n), stat=ierr)
1092 5 : end subroutine allocate_reaclib_data
1093 :
1094 :
1095 3 : subroutine free_reaclib_data(reaclib)
1096 : type(reaclib_data), intent(inout) :: reaclib
1097 3 : if (allocated(reaclib% chapter)) &
1098 : deallocate( &
1099 0 : reaclib% chapter, reaclib% species, reaclib% label, reaclib% reaction_flag, &
1100 3 : reaclib% reverse_flag, reaclib% Qvalue, reaclib% coefficients)
1101 3 : end subroutine free_reaclib_data
1102 :
1103 :
1104 10 : subroutine allocate_reaction_data(r, n, nweak, ierr)
1105 : type(reaction_data), intent(out) :: r
1106 : integer, intent(in) :: n ! number of rates
1107 : integer, intent(in) :: nweak ! number of weaklib rates
1108 : integer, intent(out) :: ierr
1109 : allocate( &
1110 : r% reaction_handle(n), r% reverse_handle(n), r% category(n), r% chapter(n), &
1111 : r% weaklib_ids(nweak), r% also_in_reaclib(nweak), &
1112 : r% pspecies(max_species_per_reaction, n), r% reaction_flag(n), &
1113 : r% weight(n), r% weight_reverse(n), r% coefficients(ncoefficients, n), &
1114 : r% weak_mask(n), r% inverse_coefficients(ninverse_coeff, n), &
1115 20 : r% inverse_exp(n), r% inverse_part(npart, n), r% Q(n), r% Qneu(n), stat=ierr)
1116 : nullify(r% reaction_dict)
1117 : nullify(r% reverse_dict)
1118 :
1119 10 : end subroutine allocate_reaction_data
1120 :
1121 :
1122 8 : subroutine free_reaction_data(r)
1123 : use utils_lib, only: integer_dict_free
1124 : type(reaction_data), intent(inout) :: r
1125 8 : if (associated(r% chapter)) then
1126 0 : deallocate( &
1127 0 : r% reaction_handle, r% reverse_handle, r% category, r% chapter, &
1128 0 : r% weaklib_ids, r% also_in_reaclib, r% pspecies, &
1129 0 : r% reaction_flag, r% weight, r% weight_reverse, r% coefficients, r% weak_mask, &
1130 8 : r% inverse_coefficients, r% inverse_exp, r% inverse_part, r% Q, r% Qneu)
1131 : end if
1132 8 : if (associated(r% reaction_dict)) call integer_dict_free(r% reaction_dict)
1133 8 : if (associated(r% reverse_dict)) call integer_dict_free(r% reverse_dict)
1134 8 : end subroutine free_reaction_data
1135 :
1136 :
1137 5 : subroutine do_start_rates_def_init(ierr)
1138 : use math_lib
1139 : integer, intent(out) :: ierr
1140 5 : ierr = 0
1141 5 : call set_rattab_range(5.30102999566398d0, 10.301029995664d0)
1142 :
1143 5 : reaclib_min_T9 = 1d-2
1144 : ! need <= 2d-3 for pre-ms li7 burning
1145 : ! pre-ms deuterium burning needs much lower (4d-4)
1146 : ! but that seems to cause problems during advanced burning.
1147 :
1148 5 : end subroutine do_start_rates_def_init
1149 :
1150 :
1151 5 : subroutine set_rattab_range(tlo, thi)
1152 5 : use math_lib
1153 : real(dp), intent(in) :: tlo, thi
1154 5 : if (abs(thi - tlo) < 1d-6) then
1155 0 : rattab_tlo = tlo
1156 0 : rattab_temp_lo = exp10(rattab_tlo)
1157 0 : rattab_thi = rattab_tlo
1158 0 : rattab_temp_hi = rattab_temp_lo
1159 0 : nrattab = 1
1160 0 : rattab_tstp = 0
1161 : return
1162 : end if
1163 5 : rattab_thi = thi
1164 5 : rattab_tlo = tlo
1165 5 : rattab_temp_hi = exp10(rattab_thi)
1166 5 : rattab_temp_lo = exp10(rattab_tlo)
1167 5 : nrattab = rattab_points_per_decade*(rattab_thi - rattab_tlo) + 1
1168 5 : if (nrattab <= 1) then
1169 0 : rattab_thi = rattab_tlo
1170 0 : nrattab = 1
1171 0 : rattab_tstp = 0
1172 : else
1173 5 : rattab_tstp = (rattab_thi-rattab_tlo)/(nrattab-1)
1174 : end if
1175 5 : end subroutine set_rattab_range
1176 :
1177 :
1178 5276 : integer function get_rates_reaction_id(reaction_name) result(value)
1179 5 : use utils_lib, only: integer_dict_lookup
1180 : character (len=*), intent(in) :: reaction_name
1181 : integer :: ierr
1182 : ierr = 0
1183 5276 : call integer_dict_lookup(reaction_names_dict, reaction_name, value, ierr)
1184 5276 : if (ierr /= 0) value = 0
1185 5276 : end function get_rates_reaction_id
1186 :
1187 :
1188 934 : subroutine create_ecapture_dict_key(ecapture_lhs, ecapture_rhs, key)
1189 : character(len=iso_name_length), intent(in) :: ecapture_lhs, ecapture_rhs
1190 : character(len=2*iso_name_length+1), intent(out) :: key
1191 934 : key = trim(ecapture_lhs) // ' ' // trim(ecapture_rhs)
1192 934 : end subroutine create_ecapture_dict_key
1193 :
1194 :
1195 224576 : subroutine create_weak_dict_key(weak_lhs, weak_rhs, key)
1196 : character(len=iso_name_length), intent(in) :: weak_lhs, weak_rhs
1197 : character(len=2*iso_name_length+1), intent(out) :: key
1198 224576 : key = trim(weak_lhs) // ' ' // trim(weak_rhs)
1199 224576 : end subroutine create_weak_dict_key
1200 :
1201 :
1202 290612 : integer function do_get_weak_rate_id(lhs, rhs) ! returns 0 if reaction not found
1203 : use utils_lib
1204 : character (len=*), intent(in) :: lhs, rhs
1205 : integer :: ierr, i
1206 : character(len=2*iso_name_length+1) :: key
1207 : character (len=iso_name_length) :: lhs_name, rhs_name
1208 145306 : ierr = 0
1209 145306 : do_get_weak_rate_id = 0
1210 145306 : lhs_name = adjustl(lhs)
1211 145306 : rhs_name = adjustl(rhs)
1212 145306 : call create_weak_dict_key(lhs_name, rhs_name, key)
1213 145306 : call integer_dict_lookup(weak_reactions_dict, key, i, ierr)
1214 145306 : if (ierr /= 0) then
1215 : !write(*,*) 'failed in integer_dict_lookup for key ' // trim(key)
1216 : return
1217 : end if
1218 3540 : do_get_weak_rate_id = i
1219 145306 : end function do_get_weak_rate_id
1220 :
1221 :
1222 146626 : integer function do_get_weak_info_list_id(lhs, rhs) ! returns 0 if reaction not found
1223 : ! value can be used to index weak_info_list_halflife and weak_info_list_Qneu
1224 : use utils_lib
1225 : character (len=*), intent(in) :: lhs, rhs ! names as in weak_info.list file
1226 : integer :: ierr, i
1227 : character(len=2*iso_name_length+1) :: key
1228 : character (len=iso_name_length) :: lhs_name, rhs_name
1229 73313 : ierr = 0
1230 73313 : do_get_weak_info_list_id = 0
1231 73313 : lhs_name = adjustl(lhs)
1232 73313 : rhs_name = adjustl(rhs)
1233 73313 : call create_weak_dict_key(lhs_name, rhs_name, key)
1234 73313 : call integer_dict_lookup(weak_info_list_dict, key, i, ierr)
1235 73313 : if (ierr /= 0) then
1236 : !write(*,'(a)') 'get_weak_info_list_id failed for ' // trim(key)
1237 : return
1238 : end if
1239 148 : do_get_weak_info_list_id = i
1240 73313 : end function do_get_weak_info_list_id
1241 :
1242 :
1243 1361336 : integer function get_num_reaction_inputs(ir)
1244 : integer, intent(in) :: ir
1245 : integer :: j
1246 : include 'formats'
1247 : if (max_num_reaction_inputs == 3) then
1248 1361336 : if (reaction_inputs(5,ir) /= 0) then
1249 : get_num_reaction_inputs = 3
1250 1041330 : else if (reaction_inputs(3,ir) /= 0) then
1251 : get_num_reaction_inputs = 2
1252 320562 : else if (reaction_inputs(1,ir) /= 0) then
1253 : get_num_reaction_inputs = 1
1254 : else
1255 108 : get_num_reaction_inputs = 0
1256 : end if
1257 : return
1258 : end if
1259 : get_num_reaction_inputs = max_num_reaction_inputs
1260 : do j = 1, 2*max_num_reaction_inputs-1, 2
1261 : if (reaction_inputs(j,ir) == 0) then
1262 : get_num_reaction_inputs = (j-1)/2
1263 : exit
1264 : end if
1265 : end do
1266 : end function get_num_reaction_inputs
1267 :
1268 :
1269 1361336 : integer function get_num_reaction_outputs(ir)
1270 : integer, intent(in) :: ir
1271 : integer :: j
1272 : if (max_num_reaction_outputs == 3) then
1273 : if (reaction_outputs(5,ir) /= 0) then
1274 : get_num_reaction_outputs = 3
1275 : else if (reaction_outputs(3,ir) /= 0) then
1276 : get_num_reaction_outputs = 2
1277 : else if (reaction_outputs(1,ir) /= 0) then
1278 : get_num_reaction_outputs = 1
1279 : else
1280 : get_num_reaction_outputs = 0
1281 : end if
1282 : return
1283 : end if
1284 1361336 : get_num_reaction_outputs = max_num_reaction_outputs
1285 2963214 : do j = 1, 2*max_num_reaction_outputs-1, 2
1286 2963214 : if (reaction_outputs(j,ir) == 0) then
1287 1361336 : get_num_reaction_outputs = (j-1)/2
1288 1361336 : exit
1289 : end if
1290 : end do
1291 : end function get_num_reaction_outputs
1292 :
1293 :
1294 0 : end module rates_def
1295 :
|
