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_lib
21 :
22 : use const_def, only: dp
23 : use utils_lib, only: mesa_error
24 :
25 : implicit none
26 :
27 :
28 : contains
29 :
30 :
31 : ! call this routine to initialize the rates module.
32 : ! only needs to be done once at start of run.
33 :
34 6 : subroutine rates_init( &
35 : reactionlist_filename, jina_reaclib_filename, &
36 : rates_table_dir_in, &
37 : use_suzuki_weak_rates, &
38 : use_special_weak_rates, &
39 : special_weak_states_file, &
40 : special_weak_transitions_file, &
41 : cache_dir, &
42 : ierr)
43 : use rates_def
44 : use reaclib_input, only: do_read_reaclib
45 : use load_weak, only: load_weak_data
46 : use load_ecapture, only: load_ecapture_data
47 : use rates_initialize, only: init_rates_info
48 :
49 : character (len=*), intent(in) :: reactionlist_filename, jina_reaclib_filename, rates_table_dir_in
50 : logical, intent(in) :: use_special_weak_rates, use_suzuki_weak_rates
51 : character (len=*), intent(in) :: special_weak_states_file, special_weak_transitions_file
52 : character (len=*), intent(in) :: cache_dir ! '' means use default
53 : integer, intent(out) :: ierr ! 0 means AOK.
54 :
55 : logical, parameter :: dbg = .false.
56 : include 'formats'
57 :
58 5 : ierr = 0
59 :
60 5 : rates_table_dir = rates_table_dir_in
61 :
62 5 : call set_rates_cache_dir(cache_dir, ierr)
63 5 : if (ierr /= 0) return
64 :
65 5 : use_suzuki_tables = use_suzuki_weak_rates
66 :
67 : if (dbg) write(*,*) 'call weaklib_init'
68 5 : call weaklib_init(ierr)
69 5 : if (ierr /= 0) return
70 5 : call load_weak_data(ierr)
71 5 : if (ierr /= 0) return
72 :
73 : ! map special weak rates controls to old names
74 5 : do_ecapture = use_special_weak_rates
75 5 : ecapture_states_file = special_weak_states_file
76 5 : ecapture_transitions_file = special_weak_transitions_file
77 :
78 : if (dbg) write(*,*) 'call ecapture_init'
79 5 : call ecapture_init(ierr)
80 5 : if (ierr /= 0) return
81 5 : if (do_ecapture) call load_ecapture_data(ierr)
82 5 : if (ierr /= 0) return
83 :
84 : if (dbg) write(*,*) 'call reaclib_init'
85 5 : call reaclib_init(jina_reaclib_filename)
86 : if (dbg) write(*,*) 'call do_read_reaclib'
87 5 : call do_read_reaclib(ierr)
88 5 : if (ierr /= 0) return
89 :
90 : if (dbg) write(*,*) 'call init_rates_info'
91 5 : call init_rates_info(reactionlist_filename, ierr)
92 5 : if (ierr /= 0) return
93 :
94 5 : have_finished_initialization = .true.
95 :
96 10 : end subroutine rates_init
97 :
98 :
99 5 : subroutine rates_warning_init( &
100 : warn_rates_for_high_temp_in, max_safe_logT_for_rates_in)
101 5 : use rates_def
102 : logical, intent(in) :: warn_rates_for_high_temp_in
103 : real(dp), intent(in) :: max_safe_logT_for_rates_in
104 : ! Setup warnings
105 5 : warn_rates_for_high_temp = warn_rates_for_high_temp_in
106 5 : max_safe_logT_for_rates = max_safe_logT_for_rates_in
107 5 : end subroutine rates_warning_init
108 :
109 :
110 21 : subroutine read_raw_rates_records(ierr)
111 5 : use rates_initialize, only: init_raw_rates_records
112 : integer, intent(out) :: ierr ! 0 means AOK.
113 : ierr = 0
114 21 : call init_raw_rates_records(ierr)
115 21 : if (ierr /= 0) then
116 0 : write(*,*) 'rates_init failed in init_raw_rates_records'
117 : return
118 : end if
119 21 : end subroutine read_raw_rates_records
120 :
121 :
122 0 : subroutine read_rate_from_file(rate_name, filename, ierr)
123 21 : use rates_initialize
124 : character(len=*), intent(in) :: rate_name, filename
125 : integer, intent(out) :: ierr
126 :
127 0 : call rate_from_file(rate_name, filename, ierr)
128 0 : if(ierr/=0) then
129 0 : write(*,*) 'read_rate_from_file for',trim(rate_name), trim(filename)
130 : return
131 : end if
132 :
133 0 : end subroutine read_rate_from_file
134 :
135 1 : subroutine read_rates_from_files(rate_names, filenames, ierr)
136 0 : use rates_initialize
137 : character(len=*), intent(in) :: rate_names(:), filenames(:)
138 : integer, intent(out) :: ierr
139 : integer :: i
140 :
141 1 : ierr = 0
142 502 : do i=1, ubound(rate_names,dim=1)
143 500 : call rate_from_file(rate_names(i), filenames(i), ierr)
144 501 : if(ierr/=0) then
145 0 : write(*,*) 'read_rates_from_files for num=',i,trim(rate_names(i)), trim(filenames(i))
146 : return
147 : end if
148 : end do
149 :
150 1 : end subroutine read_rates_from_files
151 :
152 :
153 3 : subroutine rates_shutdown
154 :
155 1 : use rates_def
156 : use rates_initialize, only: free_reaction_arrays, free_raw_rates_records
157 : use reaclib_input, only: reaclib
158 : use utils_lib
159 :
160 3 : call integer_dict_free(skip_warnings_dict)
161 3 : call integer_dict_free(reaction_names_dict)
162 :
163 3 : call free_weak_info()
164 3 : call free_ecapture_info()
165 3 : call free_reaclib_data(reaclib)
166 3 : call free_reaction_data(reaclib_rates)
167 3 : call free_reaction_arrays()
168 3 : call free_raw_rates_records()
169 :
170 3 : end subroutine rates_shutdown
171 :
172 :
173 18 : subroutine add_reaction_from_reaclib(reaction_handle, reverse_handle, indx, ierr)
174 3 : use rates_initialize, only: do_add_reaction_from_reaclib
175 : character (len=*), intent(in) :: reaction_handle ! to be added
176 : character (len=*), intent(in) :: reverse_handle ! = '' if not a reverse
177 : integer, intent(in) :: indx ! index in reaclib rates
178 : integer, intent(out) :: ierr
179 18 : call do_add_reaction_from_reaclib(reaction_handle, reverse_handle, indx, ierr)
180 18 : end subroutine add_reaction_from_reaclib
181 :
182 :
183 0 : subroutine add_reaction_for_handle(handle, ierr)
184 18 : use rates_initialize, only: do_add_reaction_for_handle
185 : character (len=*), intent(in) :: handle ! to be added
186 : integer, intent(out) :: ierr
187 0 : call do_add_reaction_for_handle(handle, ierr)
188 0 : end subroutine add_reaction_for_handle
189 :
190 :
191 19 : subroutine make_rate_tables( &
192 19 : num_reactions, cache_suffix, net_reaction_id, &
193 19 : rattab, rattab_f1, nT8s, ttab, logttab, ierr)
194 0 : use rates_support, only : do_make_rate_tables
195 : integer, intent(in) :: num_reactions, nT8s, net_reaction_id(:)
196 : character (len=*), intent(in) :: cache_suffix
197 : real(dp) :: rattab(:,:), ttab(:), logttab(:)
198 : real(dp), pointer :: rattab_f1(:)
199 : integer, intent(out) :: ierr
200 : call do_make_rate_tables( &
201 : num_reactions, cache_suffix, net_reaction_id, &
202 19 : rattab, rattab_f1, nT8s, ttab, logttab, ierr)
203 19 : end subroutine make_rate_tables
204 :
205 :
206 0 : subroutine show_reaction_rates_from_cache(cache_filename, ierr)
207 19 : use rates_support, only: do_show_reaction_from_cache
208 : character (len=*) :: cache_filename
209 : integer, intent(out) :: ierr
210 0 : call do_show_reaction_from_cache(cache_filename, ierr)
211 0 : end subroutine show_reaction_rates_from_cache
212 :
213 :
214 0 : subroutine extract_reaclib_rates(set,nuclides,rates,use_weaklib,ierr)
215 0 : use rates_def
216 : use chem_def, only: nuclide_set, nuclide_data
217 : use reaclib_input, only: do_extract_rates
218 : type(nuclide_set), dimension(:), intent(in) :: set
219 : type(nuclide_data), intent(in), target :: nuclides
220 : logical, intent(in) :: use_weaklib
221 : type(reaction_data), intent(out) :: rates
222 : integer, intent(out) :: ierr
223 0 : call do_extract_rates(set,nuclides,rates,use_weaklib,ierr)
224 0 : end subroutine extract_reaclib_rates
225 :
226 :
227 0 : subroutine output_reaclib_rates(unitno,rates,nuclides,format)
228 0 : use reaclib_print
229 : use rates_def
230 : integer, intent(in) :: unitno
231 : type(reaction_data),intent(in) :: rates
232 : type(nuclide_data), intent(in) :: nuclides
233 : integer, intent(in) :: format
234 : integer :: err
235 : err = 0
236 0 : select case (format)
237 : case(internal_format)
238 0 : call write_reaction_data(unitno,rates,err)
239 : case(pretty_print_format)
240 0 : call pretty_print_reactions(unitno,rates,nuclides,err)
241 : case(short_format)
242 0 : call print_short_format_reactions(unitno,rates,nuclides,err)
243 : end select
244 0 : end subroutine output_reaclib_rates
245 :
246 :
247 0 : subroutine reaclib_pretty_print_reaction(unitno, i, rates, nuclides, reverse, str, ierr)
248 0 : use reaclib_print, only: do_pretty_print_reaction
249 : use rates_def
250 : integer, intent(in) :: unitno, i
251 : type(reaction_data), intent(in) :: rates
252 : type(nuclide_data), intent(in) :: nuclides
253 : logical, intent(in) :: reverse
254 : character (len=100), intent(inout) :: str
255 : integer, intent(out) :: ierr
256 0 : call do_pretty_print_reaction(unitno, i, rates, nuclides, reverse, str, ierr)
257 0 : end subroutine reaclib_pretty_print_reaction
258 :
259 :
260 20 : subroutine eval_tfactors(tf, logT, temp)
261 0 : use rates_def, only : T_Factors
262 : use ratelib, only: tfactors
263 : type (T_Factors), pointer :: tf ! allocate this before calling
264 : real(dp), intent(in) :: logT, temp
265 20 : call tfactors(tf, logT, temp)
266 20 : end subroutine eval_tfactors
267 :
268 :
269 2 : subroutine get_raw_rate(ir, temp, tf, raw_rate, ierr)
270 20 : use rates_def, only : T_Factors
271 : use raw_rates
272 : integer, intent(in) :: ir
273 : real(dp), intent(in) :: temp
274 : type (T_Factors), pointer :: tf
275 : real(dp), intent(out) :: raw_rate
276 : integer, intent(out) :: ierr
277 2 : call set_raw_rate(ir, temp, tf, raw_rate, ierr)
278 2 : end subroutine get_raw_rate
279 :
280 :
281 18 : subroutine get_raw_rates(n, irs, temp, tf, rates, ierr)
282 2 : use rates_def, only : T_Factors
283 : use raw_rates, only: set_raw_rates
284 : integer, intent(in) :: n
285 : integer, intent(in) :: irs(:) ! (n) maps 1..n to reaction id
286 : real(dp), intent(in) :: temp
287 : type (T_Factors), pointer :: tf
288 : real(dp), intent(inout) :: rates(:)
289 : integer, intent(out) :: ierr
290 18 : call set_raw_rates(n, irs, temp, tf, rates, ierr)
291 18 : end subroutine get_raw_rates
292 :
293 2385 : integer function rates_reaction_id(rname)
294 18 : use rates_def, only: get_rates_reaction_id
295 : character (len=*), intent(in) :: rname ! reaction name such as 'rpp'
296 : ! returns id for the reaction if there is a matching entry in reaction_Name
297 : ! returns 0 otherwise.
298 2385 : rates_reaction_id = get_rates_reaction_id(rname)
299 2385 : end function rates_reaction_id
300 :
301 0 : integer function eval_num_reaction_inputs(ir)
302 2385 : use rates_def, only: get_num_reaction_inputs
303 : integer, intent(in) :: ir
304 0 : eval_num_reaction_inputs = get_num_reaction_inputs(ir)
305 0 : end function eval_num_reaction_inputs
306 :
307 0 : integer function eval_num_reaction_outputs(ir)
308 0 : use rates_def, only: get_num_reaction_outputs
309 : integer, intent(in) :: ir
310 0 : eval_num_reaction_outputs = get_num_reaction_outputs(ir)
311 0 : end function eval_num_reaction_outputs
312 :
313 :
314 0 : subroutine rates_eval_reaclib_21( &
315 0 : ir, temp, den, rate_raw, reverse_rate_raw, ierr)
316 0 : use rates_support, only: do_eval_reaclib_21
317 : integer, intent(in) :: ir ! reaction_id
318 : real(dp), intent(in) :: temp, den
319 : real(dp), intent(inout) :: rate_raw(:), reverse_rate_raw(:)
320 : integer, intent(out) :: ierr
321 : call do_eval_reaclib_21( &
322 0 : ir, temp, den, rate_raw, reverse_rate_raw, ierr)
323 0 : end subroutine rates_eval_reaclib_21
324 :
325 :
326 0 : subroutine rates_eval_reaclib_22( &
327 0 : ir, temp, den, rate_raw, reverse_rate_raw, ierr)
328 0 : use rates_support, only: do_eval_reaclib_22
329 : integer, intent(in) :: ir ! reaction_id
330 : real(dp), intent(in) :: temp, den
331 : real(dp), intent(inout) :: rate_raw(:), reverse_rate_raw(:)
332 : integer, intent(out) :: ierr
333 : call do_eval_reaclib_22( &
334 0 : ir, temp, den, rate_raw, reverse_rate_raw, ierr)
335 0 : end subroutine rates_eval_reaclib_22
336 :
337 :
338 0 : subroutine rates_two_to_one_coeffs_for_reverse_factor( &
339 : Q, iso_A, iso_B, iso_C, a, b, ierr)
340 0 : use chem_def, only: chem_isos
341 : use math_lib
342 : real(dp), intent(in) :: Q
343 : integer, intent(in) :: iso_A, iso_B, iso_C
344 : real(dp), intent(out) :: a, b
345 : integer, intent(out) :: ierr
346 0 : real(dp) :: W_A, W_B, W_C, g_A, g_B, g_C
347 0 : if (Q < 0) then
348 0 : ierr = -1
349 : return
350 : end if
351 0 : ierr = 0
352 0 : W_A = chem_isos% W(iso_A)
353 0 : W_B = chem_isos% W(iso_B)
354 0 : W_C = chem_isos% W(iso_C)
355 0 : g_A = 2d0*chem_isos% spin(iso_A) + 1d0
356 0 : g_B = 2d0*chem_isos% spin(iso_B) + 1d0
357 0 : g_C = 2d0*chem_isos% spin(iso_C) + 1d0
358 : ! Arnett, Supernovae and Nucleosynthesis, eqn 3.136
359 0 : a = 9.8678d9*(g_A*g_B/g_C)*pow(W_A*W_B/W_C,1.5d0)
360 0 : b = -11.605d0*Q
361 0 : end subroutine rates_two_to_one_coeffs_for_reverse_factor
362 :
363 :
364 : ! note: assumes ground state spins and requires Q > 0.
365 : ! i.e., A + B -> C exothermic
366 0 : subroutine rates_two_to_one_reverse_factor( &
367 : Q, T9, T932, iso_A, iso_B, iso_C, rev, d_rev_dT, ierr) ! A + B <-> C
368 0 : use math_lib
369 : real(dp), intent(in) :: Q, T9, T932
370 : integer, intent(in) :: iso_A, iso_B, iso_C
371 : real(dp), intent(out) :: rev, d_rev_dT
372 : integer, intent(out) :: ierr
373 : real(dp) :: a, b
374 : call rates_two_to_one_coeffs_for_reverse_factor( &
375 0 : Q, iso_A, iso_B, iso_C, a, b, ierr)
376 0 : if (ierr /= 0) return
377 0 : rev = a*T932*exp(b/T9)
378 0 : d_rev_dT = rev*(1.5d0*T9 - b)/(T9*T9*1d9)
379 0 : end subroutine rates_two_to_one_reverse_factor
380 :
381 :
382 0 : subroutine rates_two_to_two_coeffs_for_reverse_factor( &
383 : Q, iso_A, iso_B, iso_C, iso_D, a, b, ierr)
384 0 : use chem_def, only: chem_isos
385 : real(dp), intent(in) :: Q
386 : integer, intent(in) :: iso_A, iso_B, iso_C, iso_D
387 : real(dp), intent(out) :: a, b
388 : integer, intent(out) :: ierr
389 0 : real(dp) :: W_A, W_B, W_C, W_D, g_A, g_B, g_C, g_D, a1
390 0 : if (Q < 0) then
391 0 : ierr = -1
392 : return
393 : end if
394 0 : ierr = 0
395 0 : W_A = chem_isos% W(iso_A)
396 0 : W_B = chem_isos% W(iso_B)
397 0 : W_C = chem_isos% W(iso_C)
398 0 : W_D = chem_isos% W(iso_D)
399 0 : g_A = 2d0*chem_isos% spin(iso_A) + 1d0
400 0 : g_B = 2d0*chem_isos% spin(iso_B) + 1d0
401 0 : g_C = 2d0*chem_isos% spin(iso_C) + 1d0
402 0 : g_D = 2d0*chem_isos% spin(iso_D) + 0.5d0
403 : ! Arnett, Supernovae and Nucleosynthesis, eqn 3.137
404 0 : a1 = ((g_A*g_B)/(g_C*g_D))*((W_A*W_B)/(W_C*W_D))
405 0 : a = a1*sqrt(a1)
406 0 : b = -11.605d0*Q
407 0 : end subroutine rates_two_to_two_coeffs_for_reverse_factor
408 :
409 :
410 : ! note: assumes ground state spins and requires Q > 0.
411 : ! i.e., A + B -> C + D exothermic
412 0 : subroutine rates_two_to_two_reverse_factor( &
413 : Q, T9, iso_A, iso_B, iso_C, iso_D, rev, d_rev_dT, ierr) ! A + B <-> C + D
414 0 : use math_lib
415 : real(dp), intent(in) :: Q, T9
416 : integer, intent(in) :: iso_A, iso_B, iso_C, iso_D
417 : real(dp), intent(out) :: rev, d_rev_dT
418 : integer, intent(out) :: ierr
419 : real(dp) :: a, b
420 : call rates_two_to_two_coeffs_for_reverse_factor( &
421 0 : Q, iso_A, iso_B, iso_C, iso_D, a, b, ierr)
422 0 : if (ierr /= 0) return
423 0 : rev = a*exp(b/T9)
424 0 : d_rev_dT = -rev*b/(T9*T9*1d9)
425 0 : end subroutine rates_two_to_two_reverse_factor
426 :
427 :
428 4079 : logical function is_weak_reaction(ir) ! not just weaklib. any weak reaction.
429 0 : use rates_def, only: weak_reaction_info, std_reaction_neuQs
430 : integer, intent(in) :: ir ! reaction index
431 : is_weak_reaction = &
432 : (weak_reaction_info(1,ir) > 0 .and. weak_reaction_info(2,ir) > 0) .or. &
433 4079 : (std_reaction_neuQs(ir) > 0)
434 4079 : end function is_weak_reaction
435 :
436 :
437 : ! weaklib sources
438 : ! FFN: G.M. Fuller, W.A. Fowler, M.J. Newman, Ap. J. 293 (1985)
439 : ! OHMT: Oda, Hino, Muto, Takahara, and Sato. Atomic Data and Nuclear Data Tables, 1994.
440 : ! LMP: K. Langanke, G. Martínez-Pinedo / Nuclear Physics A 673 (2000) 481–508
441 :
442 175 : integer function get_weak_rate_id(lhs, rhs) ! returns 0 if reaction not found
443 4079 : use rates_def, only: do_get_weak_rate_id
444 : character (len=*), intent(in) :: lhs, rhs
445 175 : get_weak_rate_id = do_get_weak_rate_id(lhs, rhs)
446 175 : end function get_weak_rate_id
447 :
448 2 : integer function get_weak_info_list_id(lhs, rhs) ! returns 0 if reaction not found
449 : ! value can be used to index weak_info_list_halflife and weak_info_list_Qneu
450 175 : use rates_def, only: do_get_weak_info_list_id
451 : character (len=*), intent(in) :: lhs, rhs ! names as in weak_info.list file
452 2 : get_weak_info_list_id = do_get_weak_info_list_id(lhs, rhs)
453 2 : end function get_weak_info_list_id
454 :
455 :
456 : ! ecapture
457 :
458 0 : integer function get_ecapture_rate_id(lhs, rhs) ! returns 0 if reaction not found
459 2 : use rates_def
460 : use utils_lib
461 : character (len=*), intent(in) :: lhs, rhs
462 : ! names of the nuclides as given in ecapturereactions.tables (e.g. 'p', 'n', 'ca42', etc.)
463 : integer :: ierr, i
464 : character (len=2*iso_name_length+1) :: key
465 : character (len=iso_name_length) :: lhs_name, rhs_name
466 0 : ierr = 0
467 0 : get_ecapture_rate_id = 0
468 0 : lhs_name = adjustl(lhs)
469 0 : rhs_name = adjustl(rhs)
470 0 : call create_ecapture_dict_key(lhs_name, rhs_name, key)
471 0 : call integer_dict_lookup(ecapture_reactions_dict, key, i, ierr)
472 0 : if (ierr /= 0) then
473 : !write(*,*) 'failed in integer_dict_lookup for key ' // trim(key)
474 : return
475 : end if
476 0 : get_ecapture_rate_id = i
477 0 : end function get_ecapture_rate_id
478 :
479 0 : integer function get_ecapture_info_list_id(lhs, rhs) ! returns 0 if reaction not found
480 : ! value can be used to index ecapture_info_life_halflife and ecapture_info_list_Qneu
481 0 : use rates_def
482 : use utils_lib
483 : character (len=*), intent(in) :: lhs, rhs ! names as in ecapture_info.list file
484 : integer :: ierr, i
485 : character (len=2*iso_name_length+1) :: key
486 : character (len=iso_name_length) :: lhs_name, rhs_name
487 0 : ierr = 0
488 0 : get_ecapture_info_list_id = 0
489 0 : lhs_name = adjustl(lhs)
490 0 : rhs_name = adjustl(rhs)
491 0 : call create_ecapture_dict_key(lhs_name, rhs_name, key)
492 0 : call integer_dict_lookup(ecapture_reactions_dict, key, i, ierr)
493 0 : if (ierr /= 0) then
494 : !write(*,'(a)') 'get_ecapture_info_list_id failed for ' // trim(key)
495 : return
496 : end if
497 0 : get_ecapture_info_list_id = i
498 0 : end function get_ecapture_info_list_id
499 :
500 : ! reaclib
501 :
502 0 : subroutine reaclib_parse_handle(handle, num_in, num_out, iso_ids, op, ierr)
503 0 : use reaclib_support, only: do_parse_reaction_handle
504 : character (len=*), intent(in) :: handle
505 : integer, intent(out) :: num_in, num_out
506 : integer, intent(out) :: iso_ids(:) ! holds chem_ids for input and output species
507 : character (len=*), intent(out) :: op ! e.g., 'pg', 'wk', 'to', or ...
508 : integer, intent(out) :: ierr
509 0 : call do_parse_reaction_handle(handle, num_in, num_out, iso_ids, op, ierr)
510 0 : end subroutine reaclib_parse_handle
511 :
512 0 : subroutine reaclib_create_handle(num_in, num_out, iso_ids, handle)
513 0 : use reaclib_support, only: reaction_handle
514 : integer, intent(in) :: num_in, num_out
515 : integer, intent(in) :: iso_ids(:) ! holds chem_ids for input and output species
516 : character (len=*), intent(out) :: handle
517 : character (len=1), parameter :: reaction_flag = '-'
518 0 : call reaction_handle(num_in, num_out, iso_ids, reaction_flag, handle)
519 0 : end subroutine reaclib_create_handle
520 :
521 0 : subroutine reaclib_create_ec_handle(num_in, num_out, iso_ids, handle)
522 0 : use reaclib_support, only: reaction_handle
523 : integer, intent(in) :: num_in, num_out
524 : integer, intent(in) :: iso_ids(:) ! holds chem_ids for input and output species
525 : character (len=*), intent(out) :: handle
526 : character (len=1), parameter :: reaction_flag = 'e'
527 0 : call reaction_handle(num_in, num_out, iso_ids, reaction_flag, handle)
528 0 : end subroutine reaclib_create_ec_handle
529 :
530 0 : subroutine reaclib_create_wk_handle(num_in, num_out, iso_ids, handle)
531 0 : use reaclib_support, only: reaction_handle
532 : integer, intent(in) :: num_in, num_out
533 : integer, intent(in) :: iso_ids(:) ! holds chem_ids for input and output species
534 : character (len=*), intent(out) :: handle
535 : character (len=1), parameter :: reaction_flag = 'w'
536 0 : call reaction_handle(num_in, num_out, iso_ids, reaction_flag, handle)
537 0 : end subroutine reaclib_create_wk_handle
538 :
539 0 : subroutine reaclib_create_reverse_handle(num_in, num_out, iso_ids, handle)
540 0 : use reaclib_support, only: reverse_reaction_handle
541 : integer, intent(in) :: num_in, num_out
542 : integer, intent(in) :: iso_ids(:) ! holds chem_ids for input and output species
543 : character (len=*), intent(out) :: handle
544 0 : call reverse_reaction_handle(num_in, num_out, iso_ids, handle)
545 0 : end subroutine reaclib_create_reverse_handle
546 :
547 :
548 20 : integer function reaclib_lookup(handle, rates_dict) result(indx)
549 : ! returns first reaction index that matches handle.
550 : ! there may be several following that one having the same handle.
551 : ! returns 0 if handle doesn't match any of the reactions
552 0 : use rates_def
553 : use reaclib_eval, only: do_reaclib_lookup
554 : character(len=*), intent(in) :: handle ! as in rates% reaction_handle
555 : type (integer_dict), pointer :: rates_dict ! from create_reaclib_rates_dict
556 20 : indx = do_reaclib_lookup(handle, rates_dict)
557 20 : end function reaclib_lookup
558 :
559 0 : subroutine create_reaction_handle( &
560 0 : num_in, num_out, pspecies, nuclides, reverse, reaction_flag, handle)
561 20 : use reaclib_support, only: get1_reaction_handle
562 : use rates_def
563 : integer, intent(in) :: num_in, num_out
564 : integer, intent(in) :: pspecies(:)
565 : type(nuclide_data), intent(in) :: nuclides
566 : logical, intent(in) :: reverse
567 : character (len=*), intent(in) :: reaction_flag
568 : character (len=*), intent(out) :: handle
569 : call get1_reaction_handle( &
570 0 : num_in, num_out, pspecies, nuclides, reverse, reaction_flag, handle)
571 0 : end subroutine create_reaction_handle
572 :
573 :
574 0 : subroutine reaclib_indices_for_reaction(handle, rates, lo, hi, ierr)
575 0 : use reaclib_eval, only: do_reaclib_indices_for_reaction
576 : use rates_def
577 : character(len=*), intent(in) :: handle ! as in rates% reaction_handle
578 : type(reaction_data), intent(in) :: rates
579 : integer, intent(out) :: lo, hi
580 : integer, intent(out) :: ierr
581 0 : call do_reaclib_indices_for_reaction(handle, rates, lo, hi, ierr)
582 0 : end subroutine reaclib_indices_for_reaction
583 :
584 :
585 0 : subroutine reaclib_reaction_rates( &
586 : lo, hi, T9, rates, nuclides, forward_only, &
587 : lambda, dlambda_dlnT, &
588 : rlambda, drlambda_dlnT, &
589 : ierr)
590 0 : use rates_def
591 : use reaclib_eval, only: do_reaclib_reaction_rates
592 : integer, intent(in) :: lo, hi ! from reaclib_indices_for_reaction
593 : real(dp), intent(in) :: T9
594 : type(reaction_data), intent(in) :: rates
595 : type(nuclide_data), intent(in) :: nuclides
596 : logical, intent(in) :: forward_only
597 : real(dp), intent(out) :: lambda, dlambda_dlnT
598 : real(dp), intent(out) :: rlambda, drlambda_dlnT
599 : integer, intent(out) :: ierr
600 : call do_reaclib_reaction_rates( &
601 : lo, hi, T9, rates, nuclides, forward_only, &
602 : lambda, dlambda_dlnT, &
603 : rlambda, drlambda_dlnT, &
604 0 : ierr)
605 0 : end subroutine reaclib_reaction_rates
606 :
607 :
608 : ! screen
609 :
610 852 : subroutine screen_init_AZ_info( &
611 : a1, z1, a2, z2, &
612 : zs13, zhat, zhat2, lzav, aznut, zs13inv, &
613 : ierr)
614 0 : use screen5, only: screen5_init_AZ_info
615 : real(dp), intent(in) :: a1, z1, a2, z2
616 : real(dp), intent(out) :: zs13, zhat, zhat2, lzav, aznut, zs13inv
617 : integer, intent(out) :: ierr
618 852 : if (ierr /= 0) return
619 : call screen5_init_AZ_info( &
620 852 : zs13, zhat, zhat2, lzav, aznut, zs13inv, a1, z1, a2, z2, ierr)
621 852 : end subroutine screen_init_AZ_info
622 :
623 44 : integer function screening_option(which_screening_option, ierr)
624 852 : use rates_def
625 : use utils_lib, only: StrLowCase
626 : character (len=*), intent(in) :: which_screening_option
627 : integer, intent(out) :: ierr
628 : character (len=64) :: option
629 44 : ierr = 0
630 :
631 44 : option = StrLowCase(which_screening_option)
632 44 : if (associated(rates_other_screening)) then
633 : screening_option = other_screening
634 44 : else if (option == 'no_screening' .or. len_trim(option) == 0) then
635 : screening_option = no_screening
636 44 : else if (option == 'extended') then
637 : screening_option = extended_screening
638 44 : else if (option == 'salpeter') then
639 : screening_option = salpeter_screening
640 44 : else if (option == 'chugunov') then
641 : screening_option = chugunov_screening
642 : else
643 0 : ierr = -1
644 0 : screening_option = -1
645 : end if
646 44 : end function screening_option
647 :
648 8 : subroutine screening_option_str(which_screening_option, screening_option, ierr)
649 44 : use rates_def
650 : integer, intent(in) :: which_screening_option
651 : character (len=*), intent(out) :: screening_option
652 : integer, intent(out) :: ierr
653 8 : ierr = 0
654 :
655 8 : if (which_screening_option == other_screening) then
656 0 : screening_option = 'other_screening'
657 : else if (which_screening_option == no_screening) then
658 2 : screening_option = 'no_screening'
659 : else if (which_screening_option == extended_screening) then
660 2 : screening_option = 'extended'
661 : else if (which_screening_option == salpeter_screening) then
662 2 : screening_option = 'salpeter'
663 : else if (which_screening_option == chugunov_screening) then
664 2 : screening_option = 'chugunov'
665 : else
666 0 : ierr = -1
667 0 : screening_option = ''
668 : end if
669 16 : end subroutine screening_option_str
670 :
671 :
672 : ! note: if do_ecapture is true, then eval_weak_reaction_info calls this.
673 4 : subroutine eval_ecapture_reaction_info( &
674 4 : n, ids, cc, T9, YeRho, &
675 : eta, d_eta_dlnT, d_eta_dlnRho, &
676 4 : lambda, dlambda_dlnT, dlambda_dlnRho, &
677 4 : Q, dQ_dlnT, dQ_dlnRho, &
678 4 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
679 : ierr)
680 8 : use rates_def, only: Coulomb_Info
681 : use eval_ecapture, only: do_eval_ecapture_reaction_info
682 : integer, intent(in) :: n, ids(:)
683 : type(Coulomb_Info), intent(in) :: cc
684 : real(dp), intent(in) :: T9, YeRho, eta, d_eta_dlnT, d_eta_dlnRho
685 : ! lambda = combined rate (capture and decay)
686 : ! Q and Qneu are for combined rate of beta decay and electron capture.
687 : ! Q is total, so Q-Qneu is the actual thermal energy.
688 : ! note: lambdas include Ye Rho factors for electron captures.
689 : ! so treat the rates as if just beta decays
690 : real(dp), dimension(:), intent(inout) :: &
691 : lambda, dlambda_dlnT, dlambda_dlnRho, &
692 : Q, dQ_dlnT, dQ_dlnRho, &
693 : Qneu, dQneu_dlnT, dQneu_dlnRho
694 : integer, intent(out) :: ierr
695 : call do_eval_ecapture_reaction_info( &
696 : n, ids, cc, T9, YeRho, &
697 : eta, d_eta_dlnT, d_eta_dlnRho, &
698 : lambda, dlambda_dlnT, dlambda_dlnRho, &
699 : Q, dQ_dlnT, dQ_dlnRho, &
700 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
701 4 : ierr)
702 4 : end subroutine eval_ecapture_reaction_info
703 :
704 :
705 2 : subroutine eval_salpeter_screening(sc, z1, z2, scor, scordt, scordd, ierr)
706 : ! weak screening only. following Salpeter (1954),
707 : ! with equations (4-215) and (4-221) of Clayton (1968).
708 4 : use rates_def
709 : use math_lib
710 : type (Screen_Info) :: sc ! previously setup
711 : real(dp), intent(in) :: z1, z2
712 : real(dp), intent(out) :: scor ! screening factor
713 : real(dp), intent(out) :: scordt ! partial wrt temperature
714 : real(dp), intent(out) :: scordd ! partial wrt density
715 : integer, intent(out) :: ierr
716 2 : real(dp) :: zeta, lnf, rho, T, dlnf_dd, dlnf_dt
717 2 : ierr = 0
718 2 : rho = sc% den
719 2 : T = sc% temp
720 2 : zeta = (sc% z2bar + sc% zbar) / sc% abar
721 2 : lnf = 1.88d8*z1*z2*sqrt(rho*zeta/(T*T*T))
722 2 : dlnf_dd = lnf/(2*rho)
723 2 : dlnf_dt = -lnf*3/(2*T)
724 2 : scor = exp(lnf)
725 2 : scordd = scor*dlnf_dd
726 2 : scordt = scor*dlnf_dt
727 2 : end subroutine eval_salpeter_screening
728 :
729 98224 : subroutine eval_weak_reaction_info( &
730 98224 : n, ids, reaction_ids, cc, T9, YeRho, &
731 : eta, d_eta_dlnT, d_eta_dlnRho, &
732 98224 : lambda, dlambda_dlnT, dlambda_dlnRho, &
733 98224 : Q, dQ_dlnT, dQ_dlnRho, &
734 98224 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
735 : ierr)
736 : use rates_def, only: Coulomb_Info
737 : use eval_weak, only: do_eval_weak_reaction_info
738 2 : use rates_def, only : do_ecapture
739 : integer, intent(in) :: n, ids(:), reaction_ids(:)
740 : type(Coulomb_Info), pointer :: cc
741 : real(dp), intent(in) :: T9, YeRho, eta, d_eta_dlnT, d_eta_dlnRho
742 : ! lambda = combined rate (capture and decay)
743 : ! Q and Qneu are for combined rate of beta decay and electron capture.
744 : ! Q is total, so Q-Qneu is the actual thermal energy.
745 : ! note: lambdas include Ye Rho factors for electron captures.
746 : ! so treat the rates as if just beta decays
747 : real(dp), dimension(:),intent(inout) :: &
748 : lambda, dlambda_dlnT, dlambda_dlnRho, &
749 : Q, dQ_dlnT, dQ_dlnRho, &
750 : Qneu, dQneu_dlnT, dQneu_dlnRho
751 : integer, intent(out) :: ierr
752 : call do_eval_weak_reaction_info( &
753 : n, ids, reaction_ids, T9, YeRho, &
754 : eta, d_eta_dlnT, d_eta_dlnRho, &
755 : lambda, dlambda_dlnT, dlambda_dlnRho, &
756 : Q, dQ_dlnT, dQ_dlnRho, &
757 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
758 98224 : ierr)
759 98224 : if (ierr /= 0) return
760 98224 : if (.not. do_ecapture) return
761 : call eval_ecapture_reaction_info( &
762 : n, ids, cc, T9, YeRho, &
763 : eta, d_eta_dlnT, d_eta_dlnRho, &
764 : lambda, dlambda_dlnT, dlambda_dlnRho, &
765 : Q, dQ_dlnT, dQ_dlnRho, &
766 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
767 2 : ierr)
768 98224 : end subroutine eval_weak_reaction_info
769 :
770 89108 : subroutine eval_using_rate_tables( &
771 89108 : num_reactions, reaction_id, rattab, rattab_f1, nT8s, &
772 89108 : ye, logtemp, btemp, bden, raw_rate_factor, logttab, &
773 89108 : rate_raw, rate_raw_dT, rate_raw_dRho, ierr)
774 98224 : use rates_support, only : do_get_raw_rates
775 : integer, intent(in) :: num_reactions, reaction_id(:), nT8s
776 : real(dp), intent(in) :: &
777 : ye, logtemp, btemp, bden, raw_rate_factor(:), &
778 : rattab(:,:), logttab(:)
779 : real(dp), pointer :: rattab_f1(:)
780 : real(dp), intent(out), dimension(:) :: rate_raw, rate_raw_dT, rate_raw_dRho
781 : integer, intent(out) :: ierr
782 : call do_get_raw_rates(num_reactions, reaction_id, rattab, rattab_f1, nT8s, &
783 : ye, logtemp, btemp, bden, raw_rate_factor, logttab, &
784 89108 : rate_raw, rate_raw_dT, rate_raw_dRho, ierr)
785 89108 : end subroutine eval_using_rate_tables
786 :
787 : ! call this once before calling screen_pair for each reaction
788 : ! sets info that depends only on temp, den, and overall composition
789 89116 : subroutine screen_set_context( &
790 : sc, temp, den, logT, logRho, zbar, abar, z2bar, &
791 89116 : screening_mode, num_isos, y, iso_z158)
792 89108 : use rates_def
793 : use screen, only: do_screen_set_context
794 : type (Screen_Info) :: sc
795 : integer, intent(in) :: num_isos
796 : real(dp), intent(in) :: &
797 : temp, den, logT, logRho, zbar, abar, z2bar, y(:), iso_z158(:)
798 : ! y(:) = x(:)/chem_A(chem_id(:))
799 : ! iso_z(:) = chem_Z(chem_id(:))**1.58
800 : integer, intent(in) :: screening_mode
801 : call do_screen_set_context( &
802 : sc, temp, den, logT, logRho, zbar, abar, z2bar, &
803 89116 : screening_mode, num_isos, y, iso_z158)
804 :
805 89116 : end subroutine screen_set_context
806 :
807 :
808 : ! set jscr = 0 before 1st call.
809 : ! make calls in exactly the same order as for screen_init_AZ_info
810 2426148 : subroutine screen_pair( &
811 : sc, a1, z1, a2, z2, screening_mode, &
812 : zs13, zhat, zhat2, lzav, aznut, zs13inv, low_logT_lim, &
813 : scor, scordt, scordd, ierr)
814 89116 : use rates_def
815 : use screen5, only: fxt_screen5
816 : use screening_chugunov, only: eval_screen_chugunov
817 :
818 : type (Screen_Info) :: sc ! previously setup
819 : real(dp), intent(in) :: a1, z1, a2, z2
820 : integer, intent(in) :: screening_mode ! see screen_def.
821 : ! cached info
822 : real(dp), intent(in) :: zs13, zhat, zhat2, lzav, aznut, zs13inv
823 : real(dp), intent(in) :: low_logT_lim ! scor==0 for T < low_logT_lim
824 : ! outputs
825 : real(dp), intent(out) :: scor ! screening factor
826 : real(dp), intent(out) :: scordt ! partial wrt temperature
827 : real(dp), intent(out) :: scordd ! partial wrt density
828 : integer, intent(out) :: ierr
829 :
830 2426148 : if(sc% logT < low_logT_lim ) then
831 0 : scor = 0d0
832 0 : scordt=0d0
833 0 : scordd=0d0
834 : return
835 : end if
836 :
837 2426148 : if(screening_mode == other_screening) then
838 0 : call rates_other_screening(sc, z1, z2, a1, a2, scor, scordt, scordd, ierr)
839 346298 : else if (screening_mode == extended_screening) then
840 : call fxt_screen5( &
841 : sc, zs13, zhat, zhat2, lzav, aznut, zs13inv, &
842 346298 : a1, z1, a2, z2, scor, scordt, scordd, ierr)
843 2 : else if (screening_mode == salpeter_screening) then
844 2 : call eval_salpeter_screening(sc, z1, z2, scor, scordt, scordd, ierr)
845 2079846 : else if (screening_mode == chugunov_screening) then
846 2079846 : call eval_screen_chugunov(sc, z1, z2, a1, a2, scor, scordt, scordd, ierr)
847 : else if (screening_mode == no_screening) then
848 2 : scor = 1; scordt = 0; scordd = 0
849 : else
850 0 : ierr = -1
851 0 : write(*,*) 'screen_pair: unknown value for screening_mode', screening_mode
852 : end if
853 2426148 : end subroutine screen_pair
854 :
855 9106 : subroutine eval_ecapnuc_rate(etakep,temp,rho,rpen,rnep,spen,snep)
856 2426148 : use ratelib, only: ecapnuc
857 : real(dp), intent(in) :: etakep,temp,rho
858 : real(dp), intent(out) :: rpen,rnep,spen,snep
859 : ! given the electron degeneracy parameter etakep (chemical potential
860 : ! without the electron's rest mass divided by kt) and the temperature temp,
861 : ! this routine calculates rates for
862 : ! electron capture on protons rpen (captures/sec/proton),
863 : ! positron capture on neutrons rnep (captures/sec/neutron),
864 : ! and their associated neutrino energy loss rates
865 : ! spen (ergs/sec/proton) and snep (ergs/sec/neutron)
866 9106 : call ecapnuc(etakep,temp,rho,rpen,rnep,spen,snep)
867 9106 : end subroutine eval_ecapnuc_rate
868 :
869 0 : subroutine eval_mazurek_rate(btemp,bden,y56,ye,rn56ec,sn56ec)
870 9106 : use ratelib, only: mazurek
871 : real(dp), intent(in) :: btemp,bden,y56,ye
872 : real(dp), intent(out) :: rn56ec,sn56ec
873 0 : call mazurek(btemp,bden,y56,ye,rn56ec,sn56ec)
874 0 : end subroutine eval_mazurek_rate
875 :
876 2 : subroutine eval_FL_epsnuc_3alf(T, Rho, Y, UE, eps_nuc, deps_nuc_dT, deps_nuc_dRho)
877 : ! based on analytic expressions in Fushiki and Lamb, Apj, 317, 368-388, 1987.
878 :
879 : ! Note: if you plot the results of this, you'll see abrupt changes in rate at
880 : ! logRho about 9.74 and 10.25 -- these aren't bugs in the code.
881 : ! They are discussed in F&L, and show up as step functions in their expressions.
882 :
883 : ! They provide expressions for both pyconuclear regime and strong screening regime.
884 : ! The transition between the regimes happens at U = 1, where U is defined below.
885 : ! Unfortunately, at U = 1, their expressions for pycnonuclear rate and
886 : ! strong screening rate disagree!
887 : ! Bummer. For example, at logRho = 8.0, U = 1 for logT = 7.1955.
888 : ! For these values, and pure He,
889 : ! their strong screening expression is larger than their pycno expression
890 : ! by a factor of about 25.
891 :
892 : ! need to add transition region in U instead of having an abrupt change at U = 1
893 :
894 0 : use pycno, only: FL_epsnuc_3alf
895 : real(dp), intent(in) :: T ! temperature
896 : real(dp), intent(in) :: Rho ! density
897 : real(dp), intent(in) :: Y ! helium mass fraction
898 : real(dp), intent(in) :: UE ! electron molecular weight
899 : real(dp), intent(out) :: eps_nuc ! eps_nuc in ergs/g/sec
900 : real(dp), intent(out) :: deps_nuc_dT ! partial wrt temperature
901 : real(dp), intent(out) :: deps_nuc_dRho ! partial wrt density
902 2 : call FL_epsnuc_3alf(T, Rho, Y, UE, eps_nuc, deps_nuc_dT, deps_nuc_dRho)
903 2 : end subroutine eval_FL_epsnuc_3alf
904 :
905 0 : subroutine eval_n14_electron_capture_rate(T,Rho,UE,rate)
906 2 : use ratelib, only: n14_electron_capture_rate
907 : real(dp), intent(in) :: T ! temperature
908 : real(dp), intent(in) :: Rho ! density
909 : real(dp), intent(in) :: UE ! electron molecular weight
910 : real(dp), intent(out) :: rate ! (s^-1)
911 0 : call n14_electron_capture_rate(T,Rho,UE,rate)
912 0 : end subroutine eval_n14_electron_capture_rate
913 :
914 :
915 : ! call this once before calling eval_ecapture
916 : ! sets info that depends only on temp, den, and overall composition
917 89116 : subroutine coulomb_set_context( &
918 : cc, temp, den, logT, logRho, zbar, abar, z2bar)
919 0 : use rates_def, only: Coulomb_Info
920 : use coulomb, only: do_coulomb_set_context
921 : type (Coulomb_Info), pointer :: cc
922 : real(dp), intent(in) :: &
923 : temp, den, logT, logRho, zbar, abar, z2bar
924 : call do_coulomb_set_context( &
925 89116 : cc, temp, den, logT, logRho, zbar, abar, z2bar)
926 89116 : end subroutine coulomb_set_context
927 :
928 :
929 : ! translate from option string to integer
930 1 : integer function get_mui_value(option_string)
931 89116 : use eval_coulomb, only : None, DGC1973, I1993, PCR2009
932 :
933 : character(len=*), intent(in) :: option_string
934 :
935 1 : select case (trim(option_string))
936 : case('none')
937 : get_mui_value = None
938 : case('DGC1973')
939 : get_mui_value = DGC1973
940 : case('I1993')
941 : get_mui_value = I1993
942 : case('PCR2009')
943 0 : get_mui_value = PCR2009
944 : case DEFAULT
945 1 : call mesa_error(__FILE__,__LINE__,'Incorrect option for ion_coulomb_corrections')
946 : end select
947 :
948 : return
949 :
950 1 : end function get_mui_value
951 :
952 :
953 1 : integer function get_vs_value(option_string)
954 1 : use eval_coulomb, only : None, ThomasFermi, Itoh2002
955 :
956 : character(len=*), intent(in) :: option_string
957 :
958 1 : select case (trim(option_string))
959 : case('none')
960 : get_vs_value = None
961 : case('ThomasFermi')
962 : get_vs_value = ThomasFermi
963 : case('Itoh2002')
964 0 : get_vs_value = Itoh2002
965 : case DEFAULT
966 1 : call mesa_error(__FILE__,__LINE__,'Incorrect option for electron_coulomb_corrections')
967 : end select
968 :
969 : return
970 :
971 1 : end function get_vs_value
972 :
973 :
974 : end module rates_lib
|
