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 :
21 : ! last checked for consistency with Frank's ratlib.f on Sept 18, 2008
22 :
23 : module ratelib
24 : use rates_def
25 : use utils_lib
26 : use chem_def !, only: nuclide_data, chem_isos
27 : use chem_lib, only: chem_get_iso_id
28 : use const_def, only: dp, ln2, ln10, me, kerg, clight, pi2
29 : use math_lib
30 :
31 : implicit none
32 :
33 : real(dp), parameter :: lowT9_cutoff = 1d-3 ! all non-pp rates except decays go to 0 below this
34 : real(dp), parameter :: lowT9pp_cutoff = 1d-5 ! all pp rates except decays go to 0 below this
35 :
36 : real(dp) :: oneth, twoth, fourth, fiveth, elvnth, fivfour, onesix, &
37 : fivsix, sevsix, onefif, sixfif, onesev, twosev, foursev
38 : parameter (oneth = 1.0d0/3.0d0, &
39 : twoth = 2.0d0/3.0d0, &
40 : fourth = 4.0d0/3.0d0, &
41 : fiveth = 5.0d0/3.0d0, &
42 : elvnth = 11.0d0/3.0d0, &
43 : fivfour = 1.25d0, &
44 : onesix = 1.0d0/6.0d0, &
45 : fivsix = 5.0d0/6.0d0, &
46 : sevsix = 7.0d0/6.0d0, &
47 : onefif = 0.2d0, &
48 : sixfif = 1.2d0, &
49 : onesev = 1.0d0/7.0d0, &
50 : twosev = 2.0d0/7.0d0, &
51 : foursev = 4.0d0/7.0d0)
52 :
53 : integer, parameter :: nTs_rf18ap = 13
54 : logical :: have_f_rf18ap = .false.
55 : real(dp), target :: f_rf18ap_ary(4*nTs_rf18ap)
56 : real(dp), pointer :: f_rf18ap1(:), f_rf18ap(:,:)
57 :
58 : contains
59 : ! sources:
60 : ! cf88 Caughlin, G. R. & Fowler, W. A. 1988, Atom. Data and Nuc. Data Tables, 40, 283
61 : ! nacre C. Angulo et al., Nucl. Phys. A656 (1999)3-187
62 : ! wk82 wiescher and kettner, ap. j., 263, 891 (1982)
63 : ! c96 champagne 1996
64 :
65 :
66 : ! Hydrogen
67 :
68 : ! rpp, p(p,e+nu)h2
69 :
70 0 : subroutine rate_pp_fxt(tf, temp, fr, rr)
71 : type (T_Factors) :: tf
72 : real(dp), intent(in) :: temp
73 : real(dp), intent(out) :: fr, rr
74 0 : real(dp) :: term,aa,bb
75 :
76 0 : if (tf% t9 < lowT9pp_cutoff) then
77 0 : fr = 0; rr = 0; return
78 : end if
79 :
80 0 : if (tf% t9 <= 3d0) then
81 0 : aa = 4.01d-15 * tf% t9i23 * exp(-3.380d0*tf% t9i13)
82 0 : bb = 1.0d0 + 0.123d0*tf% t913 + 1.09d0*tf% t923 + 0.938d0*tf% t9
83 0 : term = aa * bb
84 : else
85 : term = 1.1581136d-15
86 : end if
87 0 : fr = term
88 0 : rr = 0.0d0
89 : end subroutine rate_pp_fxt
90 :
91 :
92 20038 : subroutine rate_pp_nacre(tf, temp, fr, rr)
93 : type (T_Factors) :: tf
94 : real(dp), intent(in) :: temp
95 : real(dp), intent(out) :: fr, rr
96 : real(dp) :: term
97 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
98 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
99 : ! + c0 T9i32 exp(-c1/T9)
100 : ! + d0 T9i32 exp(-d1/T9)
101 : ! + e0 T9^e1 exp(-e2/T9)
102 :
103 20038 : if (tf% t9 < lowT9pp_cutoff) then
104 0 : fr = 0; rr = 0; return
105 : end if
106 :
107 : call rnacre(tf, &
108 : 4.08d-15, 3.381d0, 0d0, & ! a0, a1, a2
109 : 3.82d0, 1.51d0, 0.144d0, -1.14d-02, 0d0, & ! b0, b1, b2, b3, b4
110 : 0d0, 0d0, & ! c0, c1
111 : 0d0, 0d0, & ! d0, d1
112 : 0d0, 0d0, 0d0, & ! e0, e1, e2
113 20038 : term)
114 20038 : fr = term
115 20038 : rr = 0.0d0
116 : end subroutine rate_pp_nacre
117 :
118 :
119 0 : subroutine rate_pp_jina(tf, temp, fr, rr) ! cf88
120 : type (T_Factors) :: tf
121 : real(dp), intent(in) :: temp
122 : real(dp), intent(out) :: fr, rr
123 : integer :: ierr
124 : include 'formats'
125 : ierr = 0
126 : ! p p d bet+w 1.44206d+00
127 0 : call reaclib_rate_for_handle('r_h1_h1_wk_h2', tf% T9, fr, rr, ierr)
128 0 : if (ierr /= 0) then
129 0 : write(*,'(a)') 'failed to get reaclib rate r_h1_h1_wk_h2'
130 0 : call rate_pp_fxt(tf, temp, fr, rr)
131 : end if
132 0 : end subroutine rate_pp_jina
133 :
134 :
135 : ! rpep, p(e-p, nu)h2
136 :
137 20038 : subroutine rate_pep_fxt(tf, temp, fr, rr)
138 : type (T_Factors) :: tf
139 : real(dp), intent(in) :: temp
140 : real(dp), intent(out) :: fr, rr
141 20038 : real(dp) :: term, aa, bb
142 :
143 20038 : if (tf% t9 < lowT9pp_cutoff) then
144 0 : fr = 0; rr = 0; return
145 : end if
146 :
147 20038 : if ((tf% T9) <= 3d0) then
148 16734 : aa = 1.36d-20 * (tf% T9i76) * exp(-3.380d0*(tf% T9i13))
149 16734 : bb = (1.0d0 - 0.729d0*(tf% T913) + 9.82d0*(tf% T923))
150 16734 : term = aa * bb
151 : else
152 : term = 7.3824387d-21
153 : end if
154 20038 : fr = term
155 20038 : rr = 0.0d0
156 : end subroutine rate_pep_fxt
157 :
158 :
159 0 : subroutine rate_pep_jina(tf, temp, fr, rr) ! cf88
160 : type (T_Factors) :: tf
161 : real(dp), intent(in) :: temp
162 : real(dp), intent(out) :: fr, rr
163 : integer :: ierr
164 : ierr = 0
165 : ! p p d ecw 1.44206d+00
166 0 : call reaclib_rate_for_handle('r_h1_h1_ec_h2', tf% T9, fr, rr, ierr)
167 0 : if (ierr /= 0) then
168 0 : write(*,'(a)') 'failed to get reaclib rate r_h1_h1_ec_h2'
169 0 : call rate_pep_fxt(tf, temp, fr, rr)
170 : end if
171 0 : end subroutine rate_pep_jina
172 :
173 :
174 : ! rdpg, h2(p,g)he3
175 :
176 :
177 36 : subroutine rate_dpg_fxt(tf, temp, fr, rr)
178 : type (T_Factors) :: tf
179 : real(dp), intent(in) :: temp
180 : real(dp), intent(out) :: fr, rr
181 36 : real(dp) :: term, rev, aa, bb
182 : include 'formats'
183 :
184 36 : if (tf% t9 < lowT9pp_cutoff) then
185 0 : fr = 0; rr = 0; return
186 : end if
187 :
188 36 : aa = 2.24d+03 * tf% t9i23 * exp(-3.720d0*tf% t9i13)
189 36 : bb = 1.0d0 + 0.112d0*tf% t913 + 3.38d0*tf% t923 + 2.65d0*tf% t9
190 36 : term = aa * bb
191 36 : fr = term
192 36 : rev = 1.63d+10 * tf% t932 * exp(-63.750d0*tf% t9i)
193 36 : rr = rev * term
194 : !if (temp > 3.1d6 .and. temp < 3.2d6) write(*,1) 'rates dpg', fr, temp
195 : end subroutine rate_dpg_fxt
196 :
197 :
198 20038 : subroutine rate_dpg_nacre(tf, temp, fr, rr)
199 : type (T_Factors) :: tf
200 : real(dp), intent(in) :: temp
201 : real(dp), intent(out) :: fr, rr
202 20038 : real(dp) :: term, rev
203 :
204 20038 : if (tf% t9 < lowT9pp_cutoff) then
205 0 : fr = 0; rr = 0; return
206 : end if
207 :
208 20038 : if (tf% T9 <= 0.11d0) then
209 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
210 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
211 : ! + c0 T9i32 exp(-c1/T9)
212 : ! + d0 T9i32 exp(-d1/T9)
213 : ! + e0 T9^e1 exp(-e2/T9)
214 : call rnacre(tf, &
215 : 1.81d3, 3.721d0, 0d0, & ! a0, a1, a2
216 : 14.3d0, -90.5d0, 395d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
217 : 0d0, 0d0, & ! c0, c1
218 : 0d0, 0d0, & ! d0, d1
219 : 0d0, 0d0, 0d0, & ! e0, e1, e2
220 10982 : term)
221 : else
222 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
223 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
224 : ! + c0 T9i32 exp(-c1/T9)
225 : ! + d0 T9i32 exp(-d1/T9)
226 : ! + e0 T9^e1 exp(-e2/T9)
227 : call rnacre(tf, &
228 : 2.58d3, 3.721d0, 0d0, & ! a0, a1, a2
229 : 3.96d0, 0.116d0, 0d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
230 : 0d0, 0d0, & ! c0, c1
231 : 0d0, 0d0, & ! d0, d1
232 : 0d0, 0d0, 0d0, & ! e0, e1, e2
233 9056 : term)
234 : end if
235 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
236 : 1.63d10, 63.749d0, & ! a0, a1
237 20038 : rev)
238 20038 : fr = term
239 20038 : rr = rev * term
240 : end subroutine rate_dpg_nacre
241 :
242 :
243 0 : subroutine rate_dpg_jina(tf, temp, fr, rr) ! cf88
244 : type (T_Factors) :: tf
245 : real(dp), intent(in) :: temp
246 : real(dp), intent(out) :: fr, rr
247 : ! p d he3 de04 5.49300d+00
248 0 : call jina_reaclib_2_1(ih1, ih2, ihe3, tf, fr, rr, 'rate_dpg_jina')
249 0 : end subroutine rate_dpg_jina
250 :
251 :
252 20074 : subroutine rate_png_fxt(tf, temp, fr, rr)
253 : type (T_Factors) :: tf
254 : real(dp), intent(in) :: temp
255 : real(dp), intent(out) :: fr, rr
256 20074 : real(dp) :: term, rev, aa
257 :
258 20074 : if (tf% t9 < lowT9pp_cutoff) then
259 0 : fr = 0; rr = 0; return
260 : end if
261 :
262 :
263 : ! p(n, g)d
264 : ! smith, kawano, malany 1992
265 :
266 : aa = 1.0d0 - 0.8504d0*(tf% T912) + 0.4895d0*(tf% T9) &
267 : - 0.09623d0*(tf% T932) + 8.471d-3*(tf% T92) &
268 20074 : - 2.80d-4*(tf% T952)
269 :
270 20074 : term = 4.742d4 * aa
271 :
272 : ! wagoner, schramm 1977
273 : ! aa = 1.0d0 - 0.86d0*(tf% T912) + 0.429d0*(tf% T9)
274 : ! daa = -0.5d0*0.86d0*(tf% T9i12) + 0.429
275 :
276 : ! term = 4.4d4 * aa
277 : ! dtermdt = 4.4d4 * daa
278 :
279 20074 : fr = term
280 20074 : rev = 4.71d+09 * (tf% T932) * exp(-25.82d0*(tf% T9i))
281 20074 : rr = rev * term
282 :
283 : end subroutine rate_png_fxt
284 :
285 :
286 10019 : subroutine rate_ddg_jina(tf, temp, fr, rr) ! cf88
287 : type (T_Factors) :: tf
288 : real(dp), intent(in) :: temp
289 : real(dp), intent(out) :: fr, rr
290 : ! d d he4 cf88n 2.38470d+01
291 10019 : call jina_reaclib_2_1(ih2, ih2, ihe4, tf, fr, rr, 'rate_ddg_jina')
292 10019 : end subroutine rate_ddg_jina
293 :
294 :
295 : ! Helium
296 :
297 : ! rhe3p, he3(p,e+nu)he4
298 :
299 :
300 0 : subroutine rate_hep_jina(tf, temp, fr, rr) ! cf88
301 : type (T_Factors) :: tf
302 : real(dp), intent(in) :: temp
303 : real(dp), intent(out) :: fr, rr
304 : integer :: ierr
305 : ierr = 0
306 : ! p he3 he4 bet+w 1.97960d+01
307 0 : call reaclib_rate_for_handle('r_h1_he3_wk_he4', tf% T9, fr, rr, ierr)
308 0 : if (ierr /= 0) then
309 0 : write(*,'(a)') 'failed to get reaclib rate r_h1_he3_wk_he4'
310 0 : call rate_hep_fxt(tf, temp, fr, rr)
311 : end if
312 0 : end subroutine rate_hep_jina
313 :
314 :
315 10019 : subroutine rate_hep_fxt(tf, temp, fr, rr)
316 : type (T_Factors) :: tf
317 : real(dp), intent(in) :: temp
318 : real(dp), intent(out) :: fr, rr
319 10019 : real(dp) :: term, aa
320 :
321 10019 : if (tf% t9 < lowT9pp_cutoff) then
322 0 : fr = 0; rr = 0; return
323 : end if
324 :
325 10019 : if ((tf% T9) <= 3d0) then
326 8367 : aa = 8.78d-13 * (tf% T9i23) * exp(-6.141d0*(tf% T9i13))
327 8367 : term = aa
328 : else
329 : term = 5.9733434d-15
330 : end if
331 10019 : fr = term
332 10019 : rr = 0.0d0
333 : end subroutine rate_hep_fxt
334 :
335 :
336 : ! rhe3d he3(d,p)he4 de04
337 :
338 10019 : subroutine rate_he3d_jina(tf, temp, fr, rr)
339 : type (T_Factors) :: tf
340 : real(dp), intent(in) :: temp
341 : real(dp), intent(out) :: fr, rr
342 : ! d he3 p he4 de04 1.83530d+01
343 10019 : call jina_reaclib_2_2(ih2, ihe3, ih1, ihe4, tf, fr, rr, 'rate_he3d_jina')
344 10019 : end subroutine rate_he3d_jina
345 :
346 :
347 : ! r33, he3(he3, 2p)he4
348 :
349 10037 : subroutine rate_he3he3_nacre(tf, temp, fr, rr)
350 : type (T_Factors) :: tf
351 : real(dp), intent(in) :: temp
352 : real(dp), intent(out) :: fr, rr
353 10037 : real(dp) :: term, rev
354 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
355 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
356 : ! + c0 T9i32 exp(-c1/T9)
357 : ! + d0 T9i32 exp(-d1/T9)
358 : ! + e0 T9^e1 exp(-e2/T9)
359 :
360 10037 : if (tf% t9 < lowT9pp_cutoff) then
361 0 : fr = 0; rr = 0; return
362 : end if
363 :
364 : call rnacre(tf, &
365 : 5.59d10, 12.277d0, 0d0, & ! a0, a1, a2
366 : -0.135d0, 2.54d-2, -1.29d-03, 0d0, 0d0, & ! b0, b1, b2, b3, b4
367 : 0d0, 0d0, & ! c0, c1
368 : 0d0, 0d0, & ! d0, d1
369 : 0d0, 0d0, 0d0, & ! e0, e1, e2
370 10037 : term)
371 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
372 : 3.392d-10, 149.23d0, & ! a0, a1
373 10037 : rev)
374 10037 : fr = term
375 10037 : rr = rev * term
376 : end subroutine rate_he3he3_nacre
377 :
378 0 : subroutine rate_he3he3_jina(tf, temp, fr, rr)
379 : type (T_Factors) :: tf
380 : real(dp), intent(in) :: temp
381 : real(dp), intent(out) :: fr, rr
382 : ! he3 he3 h1 h1 he4
383 0 : call jina_reaclib_2_3(ihe3, ihe3, ih1, ih1, ihe4, tf, fr, rr, 'rate_he3he3_jina')
384 0 : end subroutine rate_he3he3_jina
385 :
386 :
387 : ! r34, he4(he3,g)be7
388 :
389 :
390 0 : subroutine rate_he3he4_jina(tf, temp, fr, rr)
391 : type (T_Factors) :: tf
392 : real(dp), intent(in) :: temp
393 : real(dp), intent(out) :: fr, rr
394 : ! he4 he3 be7 de04 1.58700d+00
395 0 : call jina_reaclib_2_1(ihe4, ihe3, ibe7, tf, fr, rr, 'rate_he3he4_jina')
396 0 : end subroutine rate_he3he4_jina
397 :
398 :
399 30057 : subroutine rate_he3he4_nacre(tf, temp, fr, rr)
400 : type (T_Factors) :: tf
401 : real(dp), intent(in) :: temp
402 : real(dp), intent(out) :: fr, rr
403 30057 : real(dp) :: term, rev
404 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
405 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
406 : ! + c0 T9i32 exp(-c1/T9)
407 : ! + d0 T9i32 exp(-d1/T9)
408 : ! + e0 T9^e1 exp(-e2/T9)
409 :
410 30057 : if (tf% t9 < lowT9pp_cutoff) then
411 0 : fr = 0; rr = 0; return
412 : end if
413 :
414 : call rnacre(tf, &
415 : 5.46d6, 12.827d0, 0d0, & ! a0, a1, a2
416 : -0.307d0, 8.81d-2, -1.06d-2, 4.46d-4, 0d0, & ! b0, b1, b2, b3, b4
417 : 0d0, 0d0, & ! c0, c1
418 : 0d0, 0d0, & ! d0, d1
419 : 0d0, 0d0, 0d0, & ! e0, e1, e2
420 30057 : term)
421 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
422 : 1.113d10, 18.412d0, & ! a0, a1
423 30057 : rev)
424 30057 : fr = term
425 30057 : rr = rev * term
426 : end subroutine rate_he3he4_nacre
427 :
428 :
429 : ! r3a, triple alpha
430 :
431 :
432 20038 : subroutine rate_tripalf_jina(tf, temp, fr, rr)
433 : type (T_Factors) :: tf
434 : real(dp), intent(in) :: temp
435 : real(dp), intent(out) :: fr, rr
436 20038 : real(dp) :: fr1, rr1
437 : include 'formats'
438 :
439 : ! he4 he4 he4 c12 fy05r 7.27500d+00
440 20038 : call jina_reaclib_3_1(ihe4, ihe4, ihe4, ic12, tf, fr, rr, 'rate_tripalf_jina')
441 :
442 : return
443 : call rate_tripalf_reaclib(tf, temp, fr1, rr1)
444 :
445 : write(*,1) 'fr', fr
446 : write(*,1) 'fr1', fr1
447 : write(*,1) 'rr', rr
448 : write(*,1) 'rr1', rr1
449 : write(*,'(A)')
450 : call mesa_error(__FILE__,__LINE__,'rate_tripalf_jina')
451 : end subroutine rate_tripalf_jina
452 :
453 :
454 0 : subroutine rate_tripalf_reaclib(tf, temp, fr, rr)
455 : ! HE4(2A,G)C12 reaclib JINA - Fynbo et al. 2005 Nature 433, 136-139
456 : type (T_Factors) :: tf
457 : real(dp), intent(in) :: temp
458 : real(dp), intent(out) :: fr, rr
459 0 : real(dp) :: fr1, fr2, fr3, rev
460 0 : rr = 0
461 0 : if (temp < 1d7) then
462 0 : fr = 0; return
463 : end if
464 : call do_reaclib(tf, &
465 : -9.710520d-01, 0.000000d+00, -3.706000d+01, 2.934930d+01, &
466 : -1.155070d+02, -1.000000d+01, -1.333330d+00, &
467 0 : fr1)
468 : call do_reaclib(tf, &
469 : -2.435050d+01, -4.126560d+00, -1.349000d+01, 2.142590d+01, &
470 : -1.347690d+00, 8.798160d-02, -1.316530d+01, &
471 0 : fr2)
472 : call do_reaclib(tf, &
473 : -1.178840d+01, -1.024460d+00, -2.357000d+01, 2.048860d+01, &
474 : -1.298820d+01, -2.000000d+01, -2.166670d+00, &
475 0 : fr3)
476 0 : fr = fr1 + fr2 + fr3
477 : ! use the fxt reverse rate term
478 0 : rev = 2.00d+20*(tf% t93)*exp(-84.424d0*(tf% t9i))
479 0 : rr = fr * rev
480 : end subroutine rate_tripalf_reaclib
481 :
482 :
483 0 : subroutine rate_tripalf_nacre(tf, temp, fr, rr)
484 : type (T_Factors) :: tf
485 : real(dp), intent(in) :: temp
486 : real(dp), intent(out) :: fr, rr
487 0 : real(dp) :: r2abe, rbeac, bb, term, rev
488 : ! he4(a, g)be8
489 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
490 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
491 : ! + c0 T9i32 exp(-c1/T9)
492 : ! + d0 T9i32 exp(-d1/T9)
493 : ! + e0 T9^e1 exp(-e2/T9)
494 :
495 0 : if (tf% t9 < lowT9_cutoff) then
496 0 : fr = 0; rr = 0; return
497 : end if
498 :
499 : call rnacre(tf, &
500 : 2.43d9, 13.490d0, 1d0/0.15d0, & ! a0, a1, a2
501 : 74.5d0, 0d0, 0d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
502 : 6.09d5, 1.054d0, & ! c0, c1
503 : 0d0, 0d0, & ! d0, d1
504 : 0d0, 0d0, 0d0, & ! e0, e1, e2
505 0 : r2abe)
506 : ! be8(a, g)c12
507 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
508 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
509 : ! + c0 T9i32 exp(-c1/T9)
510 : ! + d0 T9i32 exp(-d1/T9)
511 : ! + e0 T9^e1 exp(-e2/T9)
512 : call rnacre(tf, &
513 : 2.76d7, 23.570d0, 1d0/0.4d0, & ! a0, a1, a2
514 : 5.47d0, 326d0, 0d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
515 : 130.7d0, 3.338d0, & ! c0, c1
516 : 2.51d4, 20.307d0, & ! d0, d1
517 : 0d0, 0d0, 0d0, & ! e0, e1, e2
518 0 : rbeac)
519 0 : if (tf% T9 <= 0.03d0) then
520 0 : bb = 3.07d-16*(1d0 - 29.1d0*(tf% T9) + 1308d0*(tf% T92))
521 0 : if (bb < 0) then
522 0 : bb = 0
523 : end if
524 : else
525 0 : bb = 3.44d-16*(1 + 0.0158d0*pow(tf% T9,-0.65d0))
526 : end if
527 0 : term = r2abe * rbeac * bb
528 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
529 : 2.003d20, 84.415d0, & ! a0, a1
530 0 : rev)
531 0 : fr = term
532 0 : rr = rev * term
533 : end subroutine rate_tripalf_nacre
534 :
535 :
536 0 : subroutine rate_tripalf_fxt(tf, temp, fr, rr)
537 : type (T_Factors) :: tf
538 : real(dp), intent(in) :: temp
539 : real(dp), intent(out) :: fr, rr
540 :
541 0 : real(dp) :: term, rev, r2abe, rbeac, &
542 0 : aa, bb, cc, dd, ee, ff, &
543 0 : xx, yy, zz, uu, vv, f1, rc28, &
544 : q1, q2
545 : parameter (rc28 = 0.1d0, &
546 : q1 = 1.0d0/0.009604d0, &
547 : q2 = 1.0d0/0.055225d0)
548 :
549 :
550 0 : if (tf% t9 < lowT9_cutoff) then
551 0 : fr = 0; rr = 0; return
552 : end if
553 :
554 : ! this is a(a, g)be8
555 0 : aa = 7.40d+05 * (tf% t9i32) * exp(-1.0663d0*(tf% t9i))
556 :
557 0 : bb = 4.164d+09 * (tf% t9i23) * exp(-13.49d0*(tf% t9i13) - (tf% t92)*q1)
558 :
559 : cc = 1.0d0 + 0.031d0*(tf% t913) + 8.009d0*(tf% t923) + 1.732d0*(tf% t9) &
560 0 : + 49.883d0*(tf% t943) + 27.426d0*(tf% t953)
561 :
562 0 : r2abe = aa + bb * cc
563 :
564 : ! this is be8(a, g)c12
565 0 : dd = 130.0d0 * (tf% t9i32) * exp(-3.3364d0*(tf% t9i))
566 :
567 0 : ee = 2.510d+07 * (tf% t9i23) * exp(-23.57d0*(tf% t9i13) - (tf% t92)*q2)
568 :
569 : ff = 1.0d0 + 0.018d0*(tf% t913) + 5.249d0*(tf% t923) + 0.650d0*(tf% t9) + &
570 0 : 19.176d0*(tf% t943) + 6.034d0*(tf% t953)
571 :
572 0 : rbeac = dd + ee * ff
573 :
574 : ! a factor
575 0 : xx = rc28 * 1.35d-07 * (tf% t9i32) * exp(-24.811d0*(tf% t9i))
576 :
577 :
578 : ! high temperature rate
579 0 : if ((tf% t9)>0.08d0) then
580 0 : term = 2.90d-16 * r2abe * rbeac + xx
581 :
582 :
583 : ! low temperature rate
584 : else
585 0 : uu = 0.8d0*exp(-pow(0.025d0*(tf% t9i),3.263d0))
586 0 : yy = 0.2d0 + uu ! fixes a typo in Frank's original
587 0 : vv = 4.0d0*exp(-pow((tf% t9)/0.025d0,9.227d0))
588 0 : zz = 1.0d0 + vv
589 0 : aa = 1.0d0/zz
590 0 : f1 = 0.01d0 + yy * aa ! fixes a typo in Frank's original
591 0 : term = 2.90d-16 * r2abe * rbeac * f1 + xx
592 : end if
593 :
594 0 : rev = 2.00d+20*(tf% t93)*exp(-84.424d0*(tf% t9i))
595 :
596 : ! term = 1.2d0 * term
597 : ! dtermdt = 1.2d0 * term
598 :
599 0 : fr = term
600 :
601 0 : rr = rev * term
602 :
603 : end subroutine rate_tripalf_fxt
604 :
605 :
606 20074 : subroutine rate_he3ng_fxt(tf, temp, fr, rr)
607 :
608 : type (T_Factors) :: tf
609 : real(dp), intent(in) :: temp
610 : real(dp), intent(out) :: fr, rr
611 :
612 20074 : real(dp) :: term, rev
613 :
614 :
615 20074 : if (tf% t9 < lowT9_cutoff) then
616 2796 : fr = 0; rr = 0; return
617 : end if
618 :
619 : ! he3(n, g)he4
620 17278 : term = 6.62d0 * (1.0d0 + 905.0d0*(tf% T9))
621 17278 : fr = term
622 17278 : rev = 2.61d+10 * (tf% T932) * exp(-238.81d0*(tf% T9i))
623 17278 : rr = rev * term
624 : end subroutine rate_he3ng_fxt
625 :
626 :
627 : ! Lithium
628 :
629 :
630 : ! rli7pa, li7(p,a)he4
631 :
632 10037 : subroutine rate_li7pa_nacre(tf, temp, fr, rr)
633 : type (T_Factors) :: tf
634 : real(dp), intent(in) :: temp
635 : real(dp), intent(out) :: fr, rr
636 10037 : real(dp) :: term, rev
637 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
638 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
639 : ! + c0 T9i32 exp(-c1/T9)
640 : ! + d0 T9i32 exp(-d1/T9)
641 : ! + e0 T9^e1 exp(-e2/T9)
642 :
643 10037 : if (tf% t9 < lowT9pp_cutoff) then
644 0 : fr = 0; rr = 0; return
645 : end if
646 :
647 : call rnacre(tf, &
648 : 7.20d8, 8.473d0, 1d0/6.5d0, & ! a0, a1, a2
649 : 1.05d0, -0.653d0, 0.185d0, -2.12d-2, 9.30d-4, & ! b0, b1, b2, b3, b4
650 : 0d0, 0d0, & ! c0, c1
651 : 0d0, 0d0, & ! d0, d1
652 : 9.85d6, 0.576d0, 10.415d0, & ! e0, e1, e2
653 10037 : term)
654 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
655 : 4.676d0, 201.30d0, & ! a0, a1
656 10037 : rev)
657 10037 : fr = term
658 10037 : rr = rev * term
659 : end subroutine rate_li7pa_nacre
660 :
661 :
662 0 : subroutine rate_li7pa_jina(tf, temp, fr, rr) ! jina reaclib
663 : type (T_Factors) :: tf
664 : real(dp), intent(in) :: temp
665 : real(dp), intent(out) :: fr, rr
666 0 : call jina_reaclib_2_2(ih1, ili7, ihe4, ihe4, tf, fr, rr, 'rate_li7pa_jina')
667 0 : end subroutine rate_li7pa_jina
668 :
669 :
670 : ! rli7pg, li7(p,g)be8 => 2 he4
671 :
672 :
673 : ! Beryllium
674 :
675 : ! rbe7ec, be7(e-, nu)li7
676 :
677 20038 : subroutine rate_be7em_fxt(tf, temp, fr, rr)
678 : type (T_Factors) :: tf
679 : real(dp), intent(in) :: temp
680 : real(dp), intent(out) :: fr, rr
681 20038 : real(dp) :: term, aa, bb
682 :
683 20038 : if (tf% t9 < lowT9pp_cutoff) then
684 0 : fr = 0; rr = 0; return
685 : end if
686 :
687 20038 : if (tf% T9 <= 3d0 .and. tf% T9 >= 1d-3) then
688 13938 : aa = 0.0027d0*(tf% T9i) * exp(2.515d-3*(tf% T9i))
689 13938 : bb = 1.0d0 - 0.537d0*(tf% T913) + 3.86d0*(tf% T923) + aa
690 13938 : term = 1.34d-10 * (tf% T9i12) * bb
691 : else
692 : term = 0.0d0
693 : end if
694 20038 : fr = term
695 20038 : rr = 0.0d0
696 : end subroutine rate_be7em_fxt
697 :
698 :
699 0 : subroutine rate_be7em_jina(tf, temp, fr, rr)
700 : type (T_Factors) :: tf
701 : real(dp), intent(in) :: temp
702 : real(dp), intent(out) :: fr, rr
703 : integer :: ierr
704 : ierr = 0
705 0 : call reaclib_rate_for_handle('r_be7_wk_li7', tf% T9, fr, rr, ierr)
706 0 : if (ierr /= 0) then
707 0 : write(*,'(a)') 'failed to get reaclib rate r_be7_wk_li7'
708 0 : call rate_b8ep(tf, temp, fr, rr)
709 : end if
710 0 : end subroutine rate_be7em_jina
711 :
712 :
713 : ! rbe7pg, be7(p,g)b8
714 :
715 30057 : subroutine rate_be7pg_nacre(tf, temp, fr, rr)
716 : type (T_Factors) :: tf
717 : real(dp), intent(in) :: temp
718 : real(dp), intent(out) :: fr, rr
719 30057 : real(dp) :: term, rev
720 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
721 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
722 : ! + c0 T9i32 exp(-c1/T9)
723 : ! + d0 T9i32 exp(-d1/T9)
724 : ! + e0 T9^e1 exp(-e2/T9)
725 :
726 30057 : if (tf% t9 < lowT9pp_cutoff) then
727 0 : fr = 0; rr = 0; return
728 : end if
729 :
730 : call rnacre(tf, &
731 : 2.61d5, 10.264d0, 0d0, & ! a0, a1, a2
732 : -5.11d-2, 4.68d-2, -6.60d-3, 3.12d-4, 0d0, & ! b0, b1, b2, b3, b4
733 : 2.05d3, 7.345d0, & ! c0, c1
734 : 0d0, 0d0, & ! d0, d1
735 : 0d0, 0d0, 0d0, & ! e0, e1, e2
736 30057 : term)
737 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
738 : 1.306d10, 1.594d0, & ! a0, a1
739 30057 : rev)
740 30057 : fr = term
741 30057 : rr = rev * term
742 : end subroutine rate_be7pg_nacre
743 :
744 :
745 0 : subroutine rate_be7pg_jina(tf, temp, fr, rr) ! jina reaclib cf88
746 : ! p be7 b8 cf88n 1.37000d-01
747 : type (T_Factors) :: tf
748 : real(dp), intent(in) :: temp
749 : real(dp), intent(out) :: fr, rr
750 0 : call jina_reaclib_2_1(ih1, ibe7, ib8, tf, fr, rr, 'rate_be7pg_jina')
751 0 : end subroutine rate_be7pg_jina
752 :
753 :
754 : ! rbe7dp be7(d,p)2he4
755 :
756 10019 : subroutine rate_be7dp_jina(tf, temp, fr, rr)
757 : ! d be7 p he4 he4 cf88n 1.67660d+01
758 : type (T_Factors) :: tf
759 : real(dp), intent(in) :: temp
760 : real(dp), intent(out) :: fr, rr
761 : !
762 : call jina_reaclib_2_3( &
763 10019 : ih2, ibe7, ih1, ihe4, ihe4, tf, fr, rr, 'rate_be7pg_jina')
764 10019 : end subroutine rate_be7dp_jina
765 :
766 :
767 : ! rbe7dp be7(he3,2p)2he4
768 :
769 :
770 10019 : subroutine rate_be7he3_jina(tf, temp, fr, rr)
771 : ! he3 be7 p p he4 he4 mafon 1.12721d+01
772 : type (T_Factors) :: tf
773 : real(dp), intent(in) :: temp
774 : real(dp), intent(out) :: fr, rr
775 : !
776 : call jina_reaclib_2_4( &
777 10019 : ihe3, ibe7, ih1, ih1, ihe4, ihe4, tf, fr, rr, 'rate_be7he3_jina')
778 10019 : end subroutine rate_be7he3_jina
779 :
780 :
781 : ! be9(p,d)be8 => 2a
782 :
783 :
784 : ! Boron
785 :
786 : ! rb8ep, b8(e+, nu)be8 => 2a
787 :
788 10019 : subroutine rate_b8ep(tf, temp, fr, rr)
789 : type (T_Factors) :: tf
790 : real(dp), intent(in) :: temp
791 : real(dp), intent(out) :: fr, rr
792 : real(dp), parameter :: halflife = 0.77d0 ! 770 ms
793 10019 : rr = 0.0d0
794 10019 : fr = ln2/halflife
795 : rr = 0
796 0 : end subroutine rate_b8ep
797 :
798 0 : subroutine rate_b8_wk_he4_he4_jina(tf, temp, fr, rr)
799 : type (T_Factors) :: tf
800 : real(dp), intent(in) :: temp
801 : real(dp), intent(out) :: fr, rr
802 : integer :: ierr
803 : ierr = 0
804 0 : call reaclib_rate_for_handle('r_b8_wk_he4_he4', tf% T9, fr, rr, ierr)
805 0 : if (ierr /= 0) then
806 0 : write(*,'(a)') 'failed to get reaclib rate r_b8_wk_he4_he4'
807 0 : call rate_b8ep(tf, temp, fr, rr)
808 : end if
809 0 : end subroutine rate_b8_wk_he4_he4_jina
810 :
811 :
812 : ! rb8gp, b8(g,p)be7
813 : ! see rbe7pg
814 :
815 :
816 : ! Carbon
817 :
818 : ! rc12pg, c12(p,g)n13
819 :
820 :
821 30075 : subroutine rate_c12pg_nacre(tf, temp, fr, rr)
822 : type (T_Factors) :: tf
823 : real(dp), intent(in) :: temp
824 : real(dp), intent(out) :: fr, rr
825 30075 : real(dp) :: term, rev
826 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
827 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
828 : ! + c0 T9i32 exp(-c1/T9)
829 : ! + d0 T9i32 exp(-d1/T9)
830 : ! + e0 T9^e1 exp(-e2/T9)
831 :
832 30075 : if (tf% t9 < lowT9_cutoff) then
833 4194 : fr = 0; rr = 0; return
834 : end if
835 :
836 : call rnacre(tf, &
837 : 2.00d7, 13.692d0, 1d0/0.46d0, & ! a0, a1, a2
838 : 9.89d0, -59.8d0, 266d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
839 : 1.00d5, 4.913d0, & ! c0, c1
840 : 4.24d5, 21.62d0, & ! d0, d1
841 : 0d0, 0d0, 0d0, & ! e0, e1, e2
842 25881 : term)
843 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
844 : 8.847d9, 22.553d0, & ! a0, a1
845 25881 : rev)
846 25881 : fr = term
847 25881 : rr = rev * term
848 : end subroutine rate_c12pg_nacre
849 :
850 :
851 0 : subroutine rate_c12pg_jina(tf, temp, fr, rr) ! jina reaclib nacre
852 : type (T_Factors) :: tf
853 : real(dp), intent(in) :: temp
854 : real(dp), intent(out) :: fr, rr
855 : ! p c12 n13 nacrn 1.94300d+00
856 0 : call jina_reaclib_2_1(ih1, ic12, in13, tf, fr, rr, 'rate_c12pg_jina')
857 0 : end subroutine rate_c12pg_jina
858 :
859 : ! rc12ap, c12(a,p)n15
860 :
861 30057 : subroutine rate_n15pa_jina(tf, temp, fr, rr)
862 : type (T_Factors) :: tf
863 : real(dp), intent(in) :: temp
864 : real(dp), intent(out) :: fr, rr
865 : ! p n15 he4 c12 nacrr 4.96600d+00
866 30057 : call jina_reaclib_2_2(ih1, in15, ihe4, ic12, tf, fr, rr, 'rate_n15pa_jina')
867 30057 : end subroutine rate_n15pa_jina
868 :
869 :
870 : ! rc12ag, c12(a,g)o16
871 0 : subroutine rate_c12ag_fxt(tf, temp, fr, rr)
872 : type (T_Factors) :: tf
873 : real(dp), intent(in) :: temp
874 : real(dp), intent(out) :: fr, rr
875 :
876 0 : real(dp) :: term, rev, aa, bb, cc, dd, ee, ff, gg, hh, f1, f2, zz, q1
877 : parameter (q1 = 1.0d0/12.222016d0)
878 :
879 0 : if (tf% t9 < lowT9_cutoff) then
880 0 : fr = 0; rr = 0; return
881 : end if
882 :
883 0 : aa = 1.0d0 + 0.0489d0*(tf% t9i23)
884 0 : bb = (tf% t92)*aa*aa
885 0 : cc = exp(-32.120d0*(tf% t9i13) - (tf% t92)*q1)
886 0 : dd = 1.0d0 + 0.2654d0*(tf% t9i23)
887 0 : ee = (tf% t92)*dd*dd
888 0 : ff = exp(-32.120d0*(tf% t9i13))
889 0 : gg = 1.25d3 * (tf% t9i32) * exp(-27.499d0*(tf% t9i))
890 0 : hh = 1.43d-2 * (tf% t95) * exp(-15.541d0*(tf% t9i))
891 0 : zz = 1.0d0/bb
892 0 : f1 = cc*zz
893 0 : zz = 1.0d0/ee
894 0 : f2 = ff*zz
895 0 : term = 1.04d8*f1 + 1.76d8*f2 + gg + hh
896 : ! 1.7 times cf88 value
897 0 : term = 1.7d0 * term
898 0 : rev = 5.13d10 * (tf% t932) * exp(-83.111d0*(tf% t9i))
899 0 : fr = term
900 0 : rr = rev * term
901 : end subroutine rate_c12ag_fxt
902 :
903 :
904 0 : subroutine rate_c12ag_nacre(tf, temp, fr, rr)
905 : type (T_Factors) :: tf
906 : real(dp), intent(in) :: temp
907 : real(dp), intent(out) :: fr, rr
908 0 : real(dp) :: term, rev, termE1, termE2, termRes, aa, bb, cc
909 : include 'formats'
910 :
911 0 : if (tf% t9 < lowT9_cutoff) then
912 0 : fr = 0; rr = 0; return
913 : end if
914 :
915 : ! note: uses T9i2 instead of T9i23, so special case it.
916 0 : aa = 6.66d7 * (tf% T9i2) * exp(-32.123d0*(tf% T9i13) - (tf% T92)/(4.6d0*4.6d0))
917 0 : bb = 1 + 2.54d0*(tf% T9) + 1.04d0*(tf% T92) - 0.226d0*(tf% T93)
918 0 : if (bb < 0) bb = 0
919 0 : cc = 1.39d3 * (tf% T9i32) * exp(-28.930d0*(tf% T9i))
920 0 : termE1 = aa * bb + cc
921 0 : aa = 6.56d7 * (tf% T9i2) * exp(-32.123d0*(tf% T9i13) - (tf% T92)/(1.3d0*1.3d0))
922 0 : bb = 1 + 9.23d0*(tf% T9) - 13.7d0*(tf% T92) + 7.4d0*(tf% T93)
923 0 : termE2 = aa * bb
924 0 : termRes = 19.2d0 * (tf% T92) * exp(-26.9d0*(tf% T9i))
925 0 : term = termE1 + termE2 + termRes
926 0 : rev = 5.132d10 * (tf% T932) * exp(-83.109d0*(tf% T9i))
927 0 : fr = term
928 0 : rr = rev * term
929 : end subroutine rate_c12ag_nacre
930 :
931 :
932 0 : subroutine rate_c12ag_kunz(tf, temp, fr, rr)
933 : ! kunz et al (2002)
934 : type (T_Factors) :: tf
935 : real(dp), intent(in) :: temp
936 : real(dp), intent(out) :: fr, rr
937 0 : real(dp) :: term, rev, aa, bb, cc, dd, ee
938 : real(dp), parameter :: &
939 : a0 = 1.21d8, a1 = 6.06d-2, a2 = 32.12d0, a3 = 1.7d0, a4 = 7.4d8, &
940 : a5 = 0.47d0, a6 = 32.12d0, a9tilda = 3.06d10, a11 = 38.534d0
941 : include 'formats'
942 :
943 0 : if (tf% t9 < lowT9_cutoff) then
944 0 : fr = 0; rr = 0; return
945 : end if
946 :
947 0 : aa = a0 * (tf% T9i2) * exp(-a2*(tf% T9i13) - (tf% T92)/(a3*a3))
948 0 : bb = 1 / pow(1 + a1*(tf% T9i23),2)
949 0 : cc = a4 * (tf% T9i2) * exp(-a6*(tf% T9i13))
950 0 : dd = 1 / pow(1 + a5*(tf% T9i23),2)
951 0 : ee = a9tilda * (tf% T9i13) * exp(-a11*(tf% T9i13))
952 0 : term = aa*bb + cc*dd + ee
953 0 : rev = 5.132d10 * (tf% T932) * exp(-83.109d0*(tf% T9i))
954 0 : fr = term
955 0 : rr = rev * term
956 : end subroutine rate_c12ag_kunz
957 :
958 :
959 10019 : subroutine rate_c12ag_jina(tf, temp, fr, rr)
960 : type (T_Factors) :: tf
961 : real(dp), intent(in) :: temp
962 : real(dp), intent(out) :: fr, rr
963 : ! he4 c12 o16 bu96n 7.16192d+00
964 10019 : call jina_reaclib_2_1(ihe4, ic12, io16, tf, fr, rr, 'rate_c12ag_jina')
965 10019 : end subroutine rate_c12ag_jina
966 :
967 :
968 : ! r1212p, c12(c12,p)na23
969 :
970 0 : subroutine rate_c12c12p_jina(tf, temp, fr, rr)
971 : type (T_Factors) :: tf
972 : real(dp), intent(in) :: temp
973 : real(dp), intent(out) :: fr, rr
974 : ! c12 c12 p na23 cf88r 2.24200d+00
975 0 : call jina_reaclib_2_2(ic12, ic12, ih1, ina23, tf, fr, rr, 'rate_c12c12p_jina')
976 0 : end subroutine rate_c12c12p_jina
977 :
978 : ! r1212a, c12(c12,a)ne20
979 :
980 0 : subroutine rate_c12c12_fxt(tf, temp, fr, rr)
981 : type (T_Factors) :: tf
982 : real(dp), intent(in) :: temp
983 : real(dp), intent(out) :: fr, rr
984 0 : real(dp) :: term, T9a, dt9a, T9a13, T9a56, aa, zz
985 0 : if (tf% t9 < lowT9_cutoff) then
986 0 : fr = 0; rr = 0; return
987 : end if
988 0 : aa = 1.0d0 + 0.0396d0*tf% t9
989 0 : zz = 1.0d0/aa
990 0 : t9a = tf% t9*zz
991 0 : dt9a = (1.0d0 - t9a*0.0396d0)*zz
992 0 : zz = dt9a/t9a
993 0 : t9a13 = pow(t9a,oneth)
994 0 : t9a56 = pow(t9a,fivsix)
995 0 : term = 4.27d+26 * t9a56 * tf% t9i32 * exp(-84.165d0/t9a13 - 2.12d-03*tf% t93)
996 0 : fr = term
997 0 : rr = 0.0d0
998 0 : return
999 : end subroutine rate_c12c12_fxt
1000 :
1001 :
1002 0 : subroutine rate_c12c12_fxt_basic(tf, temp, fr, rr)
1003 : type (T_Factors) :: tf
1004 : real(dp), intent(in) :: temp
1005 : real(dp), intent(out) :: fr, rr
1006 0 : real(dp) :: term,t9a,t9a13,t9a56,aa,zz
1007 :
1008 : include 'formats'
1009 :
1010 0 : if (tf% t9 < lowT9_cutoff) then
1011 0 : fr = 0; rr = 0; return
1012 : end if
1013 :
1014 0 : aa = 1.0d0 + 0.0396d0*tf% t9
1015 0 : zz = 1.0d0/aa
1016 0 : t9a = tf% t9*zz
1017 0 : t9a13 = pow(t9a,oneth)
1018 0 : t9a56 = pow(t9a,fivsix)
1019 0 : term = 4.27d+26 * t9a56 * tf% t9i32 * exp(-84.165d0/t9a13 - 2.12d-03*tf% t93)
1020 0 : fr = term
1021 0 : rr = 0.0d0
1022 :
1023 : end subroutine rate_c12c12_fxt_basic
1024 :
1025 :
1026 10019 : subroutine rate_c12c12_fxt_multi(tf, temp, fr, rr)
1027 : type (T_Factors) :: tf
1028 : real(dp), intent(in) :: temp
1029 : real(dp), intent(out) :: fr, rr
1030 10019 : real(dp) :: fr1, rr1, fr2, rr2
1031 10019 : call rate_c12c12npa(tf, temp, fr1, rr1, fr2, rr2, fr, rr)
1032 10019 : fr = fr + fr1 + fr2
1033 10019 : rr = rr + rr1 + rr2
1034 10019 : end subroutine rate_c12c12_fxt_multi
1035 :
1036 0 : subroutine rate_c12c12a_fxt(tf, temp, fr, rr)
1037 : type (T_Factors) :: tf
1038 : real(dp), intent(in) :: temp
1039 : real(dp), intent(out) :: fr, rr
1040 : real(dp) :: fr1, rr1, fr2, rr2
1041 0 : call rate_c12c12npa(tf, temp, fr1, rr1, fr2, rr2, fr, rr)
1042 0 : end subroutine rate_c12c12a_fxt
1043 :
1044 0 : subroutine rate_c12c12a_jina(tf, temp, fr, rr)
1045 : type (T_Factors) :: tf
1046 : real(dp), intent(in) :: temp
1047 : real(dp), intent(out) :: fr, rr
1048 : ! c12 c12 he4 ne20 cf88r 4.62100d+00
1049 0 : call jina_reaclib_2_2(ic12, ic12, ihe4, ine20, tf, fr, rr, 'rate_c12c12a_jina')
1050 0 : end subroutine rate_c12c12a_jina
1051 :
1052 10019 : subroutine rate_c12c12npa(tf, temp, &
1053 : fr1, rr1, & ! c12(c12,n)mg23
1054 : fr2, rr2, & ! c12(c12,p)na23
1055 : fr3, rr3) ! c12(c12,a)ne20
1056 :
1057 : type (T_Factors) :: tf
1058 : real(dp) :: temp, fr1, rr1, fr2, rr2, fr3, rr3
1059 :
1060 10019 : real(dp) :: term, rev, T9a, T9a13, &
1061 10019 : T9a56, aa, bb, cc, dd, &
1062 10019 : b24n, b24p, b24a
1063 :
1064 :
1065 10019 : if (tf% t9 < lowT9_cutoff) then
1066 1398 : fr1 = 0; rr1 = 0
1067 1398 : fr2 = 0; rr2 = 0
1068 1398 : fr3 = 0; rr3 = 0
1069 1398 : return
1070 : end if
1071 :
1072 8621 : aa = 1.0d0 + 0.0396d0*(tf% T9)
1073 :
1074 8621 : T9a = (tf% T9)/aa
1075 :
1076 8621 : T9a13 = pow(T9a,oneth)
1077 :
1078 8621 : T9a56 = pow(T9a,fivsix)
1079 :
1080 : aa = 4.27d+26 * T9a56 * (tf% T9i32) * &
1081 8621 : exp(-84.165d0/T9a13 - 2.12d-03*(tf% T93))
1082 :
1083 : ! neutron branching from dayras switkowski and woosley 1976
1084 8621 : if ((tf% T9) >= 1.5d0) then
1085 :
1086 2254 : bb = 0.055d0 * exp(0.976d0 - 0.789d0*(tf% T9))
1087 :
1088 2254 : b24n = 0.055d0 - bb
1089 :
1090 : else
1091 :
1092 6367 : bb = 1.0d0 + 0.0789d0*(tf% T9) + 7.74d0*(tf% T92)
1093 :
1094 6367 : cc = 0.766d0*(tf% T9i3)
1095 :
1096 6367 : dd = bb * cc
1097 :
1098 6367 : b24n = 0.859d0*exp(-dd)
1099 :
1100 : end if
1101 :
1102 : ! proton branching ratio
1103 8621 : if ((tf% T9)>3d0) then
1104 :
1105 1652 : b24p = oneth*(1.0d0 - b24n)
1106 :
1107 1652 : b24a = 2.0d0 * b24p
1108 :
1109 : else
1110 :
1111 6969 : b24p = 0.5d0*(1.0d0 - b24n)
1112 :
1113 6969 : b24a = b24p
1114 :
1115 : end if
1116 :
1117 : ! c12(c12, n)mg23
1118 8621 : term = aa * b24n
1119 8621 : fr1 = term
1120 :
1121 8621 : rev = 0.0d0
1122 8621 : if ((tf% T9) > 0.1d0) then
1123 4611 : rev = 3.93d0 * exp(30.16100515d0*(tf% T9i))
1124 : end if
1125 8621 : rr1 = rev * term
1126 :
1127 :
1128 : ! c12(c12, p)na23
1129 8621 : term = aa * b24p
1130 8621 : fr2 = term
1131 :
1132 8621 : rev = 0.0d0
1133 8621 : if ((tf% T9) > 0.1d0) then
1134 4611 : rev = 3.93d0 * exp(-25.98325915d0*(tf% T9i))
1135 : end if
1136 8621 : rr2 = rev * term
1137 :
1138 :
1139 : ! c12(c12, a)ne20
1140 8621 : term = aa * b24a
1141 8621 : fr3 = term
1142 :
1143 8621 : rev = 0.0d0
1144 8621 : if ((tf% T9) > 0.1d0) then
1145 4611 : rev = 2.42d0 * exp(-53.576110995d0*(tf% T9i))
1146 : end if
1147 8621 : rr3 = rev * term
1148 :
1149 : end subroutine rate_c12c12npa
1150 :
1151 :
1152 : ! r1216g, c12(o16,g)si28
1153 :
1154 0 : subroutine rate_c12o16_to_mg24_fxt(tf, temp, fr, rr)
1155 : ! this is a combined rate for reactions to mg24 and si28
1156 : type (T_Factors) :: tf
1157 : real(dp), intent(in) :: temp
1158 : real(dp), intent(out) :: fr, rr
1159 0 : call rate_c12o16_fxt(tf, temp, fr, rr)
1160 0 : fr = 0.5d0*fr
1161 0 : end subroutine rate_c12o16_to_mg24_fxt
1162 :
1163 :
1164 0 : subroutine rate_c12o16_to_si28_fxt(tf, temp, fr, rr)
1165 : ! this is a combined rate for reactions to mg24 and si28
1166 : type (T_Factors) :: tf
1167 : real(dp), intent(in) :: temp
1168 : real(dp), intent(out) :: fr, rr
1169 0 : call rate_c12o16_fxt(tf, temp, fr, rr)
1170 0 : fr = 0.5d0*fr
1171 0 : end subroutine rate_c12o16_to_si28_fxt
1172 :
1173 :
1174 10019 : subroutine rate_c12o16_fxt(tf, temp, fr, rr)
1175 : ! this is a combined rate for reactions to mg24 and si28
1176 : type (T_Factors) :: tf
1177 : real(dp), intent(in) :: temp
1178 : real(dp), intent(out) :: fr, rr
1179 10019 : real(dp) :: term, T9a, T9a13, T9a23, T9a56, aa, bb, cc
1180 :
1181 10019 : if (tf% t9 < lowT9_cutoff) then
1182 1398 : fr = 0; rr = 0; return
1183 : end if
1184 :
1185 : ! c12 + o16 reaction; see cf88 references 47-4
1186 8621 : if ((tf% T9)>=0.5d0) then
1187 3211 : aa = 1.0d0 + 0.055d0*(tf% T9)
1188 3211 : T9a = (tf% T9)/aa
1189 3211 : T9a13 = pow(T9a,oneth)
1190 3211 : T9a23 = T9a13*T9a13
1191 3211 : T9a56 = pow(T9a,fivsix)
1192 3211 : aa = exp(-0.18d0*T9a*T9a)
1193 3211 : bb = 1.06d-03*exp(2.562d0*T9a23)
1194 3211 : cc = aa + bb
1195 3211 : term = 1.72d+31 * T9a56 * (tf% T9i32) * exp(-106.594d0/T9a13)/cc
1196 : else
1197 : ! term = 2.6288035d-29
1198 : term = 0.0d0
1199 : end if
1200 8621 : fr = term
1201 8621 : rr = 0.0d0
1202 : end subroutine rate_c12o16_fxt
1203 :
1204 0 : subroutine rate_c12o16_jina(tf, temp, fr, rr)
1205 : type (T_Factors) :: tf
1206 : real(dp), intent(in) :: temp
1207 : real(dp), intent(out) :: fr, rr
1208 : real(dp) :: fr1, rr1
1209 0 : call rate_c12o16a_jina(tf, temp, fr, rr)
1210 0 : call rate_c12o16p_jina(tf, temp, fr1, rr1)
1211 0 : fr = fr + fr1
1212 0 : rr = rr + rr1
1213 0 : end subroutine rate_c12o16_jina
1214 :
1215 :
1216 0 : subroutine rate_c12o16p_fxt(tf, temp, fr, rr)
1217 : type (T_Factors) :: tf
1218 : real(dp), intent(in) :: temp
1219 : real(dp), intent(out) :: fr, rr
1220 0 : call rate_c12o16_fxt(tf, temp, fr, rr)
1221 0 : fr = 0.5d0*fr
1222 0 : end subroutine rate_c12o16p_fxt
1223 :
1224 :
1225 : ! r1216n, c12(o16,n)si27
1226 0 : subroutine rate_c12o16n(tf, temp, fr, rr)
1227 : type (T_Factors) :: tf
1228 : real(dp), intent(in) :: temp
1229 : real(dp), intent(out) :: fr, rr
1230 : real(dp) :: fr2, rr2, fr3, rr3
1231 0 : call rate_c12o16npa(tf, temp, fr, rr, fr2, rr2, fr3, rr3)
1232 0 : end subroutine rate_c12o16n
1233 :
1234 : ! r1216p, c12(o16,p)al27
1235 0 : subroutine rate_c12o16p(tf, temp, fr, rr)
1236 : type (T_Factors) :: tf
1237 : real(dp), intent(in) :: temp
1238 : real(dp), intent(out) :: fr, rr
1239 : real(dp) :: fr1, rr1, fr3, rr3
1240 0 : call rate_c12o16npa(tf, temp, fr1, rr1, fr, rr, fr3, rr3)
1241 0 : end subroutine rate_c12o16p
1242 :
1243 18 : subroutine rate_c12o16p_jina(tf, temp, fr, rr)
1244 : type (T_Factors) :: tf
1245 : real(dp), intent(in) :: temp
1246 : real(dp), intent(out) :: fr, rr
1247 : ! c12 o16 p al27 cf88r 5.17100d+00
1248 18 : call jina_reaclib_2_2(ic12, io16, ih1, ial27, tf, fr, rr, 'rate_c12o16p_jina')
1249 18 : end subroutine rate_c12o16p_jina
1250 :
1251 18 : subroutine rate_c12o16a_jina(tf, temp, fr, rr)
1252 : type (T_Factors) :: tf
1253 : real(dp), intent(in) :: temp
1254 : real(dp), intent(out) :: fr, rr
1255 : ! c12 o16 he4 mg24 cf88r 6.77100d+00
1256 18 : call jina_reaclib_2_2(ic12, io16, ihe4, img24, tf, fr, rr, 'rate_c12o16a_jina')
1257 18 : end subroutine rate_c12o16a_jina
1258 :
1259 :
1260 : ! r1216n, c12(o16,n)si27
1261 : ! r1216p, c12(o16,p)al27
1262 : ! r1216a, c12(o16,a)mg24
1263 :
1264 0 : subroutine rate_c12o16npa(tf, temp, &
1265 : fr1, rr1, & ! c12(o16,n)si27
1266 : fr2, rr2, & ! c12(o16,p)al27
1267 : fr3, rr3) ! c12(o16,a)mg24
1268 :
1269 : type (T_Factors) :: tf
1270 : real(dp) :: temp, fr1, rr1, fr2, rr2, fr3, rr3
1271 :
1272 0 : real(dp) :: term, rev, T9a, T9a13, &
1273 0 : T9a23, T9a56, aa, bb, cc, &
1274 0 : dd, b27n, b27p, b24a
1275 :
1276 :
1277 0 : if (tf% t9 < lowT9_cutoff) then
1278 0 : fr1 = 0; rr1 = 0
1279 0 : fr2 = 0; rr2 = 0
1280 0 : fr3 = 0; rr3 = 0
1281 0 : return
1282 : end if
1283 :
1284 0 : if ((tf% T9)>=0.5d0) then
1285 0 : aa = 1.0d0 + 0.055d0*(tf% T9)
1286 0 : T9a = (tf% T9)/aa
1287 0 : t9a13 = pow(t9a,oneth)
1288 0 : T9a23 = T9a13*T9a13
1289 0 : t9a56 = pow(t9a,fivsix)
1290 0 : aa = exp(-0.18d0*T9a*T9a)
1291 0 : bb = 1.06d-03*exp(2.562d0*T9a23)
1292 0 : cc = aa + bb
1293 0 : dd = 1.72d+31 * T9a56 * (tf% T9i32) * exp(-106.594d0/T9a13)/cc
1294 : else
1295 : ! dd = 2.6288035d-29
1296 : dd = 0.0d0
1297 : end if
1298 :
1299 : ! branching ratios from pwnsz data
1300 0 : b27n = 0.1d0
1301 0 : b27p = 0.5d0
1302 0 : b24a = 0.4d0
1303 :
1304 : ! c12(o16,n)si27
1305 0 : term = dd * b27n
1306 0 : fr1 = term
1307 :
1308 0 : rev = 0.0d0
1309 0 : if ((tf% T9) > 0.1d0) then
1310 0 : rev = 1.58d0 * exp(4.8972467d0*(tf% T9i))
1311 : end if
1312 0 : rr1 = rev * term
1313 :
1314 : ! c12(o16,p)al27
1315 0 : term = dd * b27p
1316 0 : fr2 = term
1317 :
1318 0 : rev = 0.0d0
1319 0 : if ((tf% T9) > 0.1d0) then
1320 0 : rev = 1.58d0 * exp(-59.9970745d0*(tf% T9i))
1321 : end if
1322 0 : rr2 = rev * term
1323 :
1324 : ! c12(o16,a)mg24
1325 0 : term = dd * b24a
1326 0 : fr3 = term
1327 0 : rev = 0.0d0
1328 0 : if ((tf% T9) > 0.1d0) then
1329 0 : rev = 2.83d0 * exp(-78.5648345d0*(tf% T9i))
1330 : end if
1331 0 : rr3 = rev * term
1332 :
1333 : end subroutine rate_c12o16npa
1334 :
1335 : ! rc13pg, c13(p,g)n14
1336 :
1337 20056 : subroutine rate_c13pg_nacre(tf, temp, fr, rr)
1338 : type (T_Factors) :: tf
1339 : real(dp), intent(in) :: temp
1340 : real(dp), intent(out) :: fr, rr
1341 20056 : real(dp) :: term, rev, bb, gs
1342 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
1343 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
1344 : ! + c0 T9i32 exp(-c1/T9)
1345 : ! + d0 T9i32 exp(-d1/T9)
1346 : ! + e0 T9^e1 exp(-e2/T9)
1347 :
1348 20056 : if (tf% t9 < lowT9_cutoff) then
1349 2796 : fr = 0; rr = 0; return
1350 : end if
1351 :
1352 : call rnacre(tf, &
1353 : 9.57d7, 13.72d0, 1d0, & ! a0, a1, a2
1354 : 3.56d0, 0d0, 0d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
1355 : 1.50d6, 5.930d0, & ! c0, c1
1356 : 0d0, 0d0, & ! d0, d1
1357 : 6.83d5, -0.864d0, 12.057d0, & ! e0, e1, e2
1358 17260 : gs)
1359 17260 : bb = 2.070d0 * exp(-37.938d0*(tf% T9i))
1360 17260 : if (bb > 1) then ! guard against rate going negative
1361 0 : bb = 1
1362 : end if
1363 17260 : term = gs * (1 - bb)
1364 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
1365 : 1.190d10, 87.619d0, & ! a0, a1
1366 17260 : rev)
1367 17260 : fr = term
1368 17260 : rr = rev * term
1369 : end subroutine rate_c13pg_nacre
1370 :
1371 :
1372 0 : subroutine rate_c13pg_jina(tf, temp, fr, rr)
1373 : type (T_Factors) :: tf
1374 : real(dp), intent(in) :: temp
1375 : real(dp), intent(out) :: fr, rr
1376 : ! p c13 n14 nacrr 7.55100d+00
1377 0 : call jina_reaclib_2_1(ih1, ic13, in14, tf, fr, rr, 'rate_c13pg_jina')
1378 0 : end subroutine rate_c13pg_jina
1379 :
1380 : ! rc13an, c13(a,n)o16
1381 0 : subroutine rate_c13an_jina(tf, temp, fr, rr)
1382 : type (T_Factors) :: tf
1383 : real(dp), intent(in) :: temp
1384 : real(dp), intent(out) :: fr, rr
1385 : ! he4 c13 n o16 nacrn 2.21600d+00
1386 0 : call jina_reaclib_2_2(ihe4, ic13, ineut, io16, tf, fr, rr, 'rate_c13an_jina')
1387 0 : end subroutine rate_c13an_jina
1388 :
1389 0 : subroutine rate_c13an_fxt(tf, temp, fr, rr)
1390 : type (T_Factors) :: tf
1391 : real(dp), intent(in) :: temp
1392 : real(dp), intent(out) :: fr, rr
1393 0 : real(dp) :: term, rev, aa, bb, cc, dd, ee, ff, gg, q1
1394 : parameter (q1 = 1.0d0/1.648656d0)
1395 :
1396 0 : if (tf% t9 < lowT9_cutoff) then
1397 0 : fr = 0; rr = 0; return
1398 : end if
1399 :
1400 0 : aa = 6.77d+15 * (tf% T9i23) * exp(-32.329d0*(tf% T9i13) - (tf% T92)*q1)
1401 0 : bb = 1.0d0 + 0.013d0*(tf% T913) + 2.04d0*(tf% T923) + 0.184d0*(tf% T9)
1402 0 : cc = aa * bb
1403 0 : dd = 3.82d+05 * (tf% T9i32) * exp(-9.373d0*(tf% T9i))
1404 0 : ee = 1.41d+06 * (tf% T9i32) * exp(-11.873d0*(tf% T9i))
1405 0 : ff = 2.0d+09 * (tf% T9i32) * exp(-20.409d0*(tf% T9i))
1406 0 : gg = 2.92d+09 * (tf% T9i32) * exp(-29.283d0*(tf% T9i))
1407 0 : term = cc + dd + ee + ff + gg
1408 0 : fr = term
1409 0 : rev = 5.79d+00 * exp(-25.711d0*(tf% T9i))
1410 0 : rr = rev * term
1411 : end subroutine rate_c13an_fxt
1412 :
1413 :
1414 : ! Nitrogen
1415 :
1416 :
1417 : ! rn13pg, n13(p,g)o14
1418 :
1419 20038 : subroutine rate_n13pg_nacre(tf, temp, fr, rr)
1420 : type (T_Factors) :: tf
1421 : real(dp), intent(in) :: temp
1422 : real(dp), intent(out) :: fr, rr
1423 20038 : real(dp) :: term, rev
1424 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
1425 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
1426 : ! + c0 T9i32 exp(-c1/T9)
1427 : ! + d0 T9i32 exp(-d1/T9)
1428 : ! + e0 T9^e1 exp(-e2/T9)
1429 :
1430 20038 : if (tf% t9 < lowT9_cutoff) then
1431 2796 : fr = 0; rr = 0; return
1432 : end if
1433 :
1434 : call rnacre(tf, &
1435 : 4.02d7, 15.205d0, 1d0/0.54d0, & ! a0, a1, a2
1436 : 3.81d0, 18.6d0, 32.3d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
1437 : 0d0, 0d0, & ! c0, c1
1438 : 0d0, 0d0, & ! d0, d1
1439 : 3.25d5, -1.35d0, 5.926d0, & ! e0, e1, e2
1440 17242 : term)
1441 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
1442 : 3.571d10, 53.705d0, & ! a0, a1
1443 17242 : rev)
1444 17242 : fr = term
1445 17242 : rr = rev * term
1446 : end subroutine rate_n13pg_nacre
1447 :
1448 0 : subroutine rate_n13pg_jina(tf, temp, fr, rr)
1449 : type (T_Factors) :: tf
1450 : real(dp), intent(in) :: temp
1451 : real(dp), intent(out) :: fr, rr
1452 : ! p n13 o14 lg06n 4.62797d+00
1453 0 : call jina_reaclib_2_1(ih1, in13, io14, tf, fr, rr, 'rate_n13pg_jina')
1454 0 : end subroutine rate_n13pg_jina
1455 :
1456 :
1457 : ! rn13ap n13(a,p)o16 cf88
1458 0 : subroutine rate_n13ap_jina(tf, temp, fr, rr)
1459 : type (T_Factors) :: tf
1460 : real(dp), intent(in) :: temp
1461 : real(dp), intent(out) :: fr, rr
1462 : ! he4 n13 p o16 cf88n 5.21800d+00
1463 0 : call jina_reaclib_2_2(ihe4, in13, ih1, io16, tf, fr, rr, 'rate_n13ap_jina')
1464 0 : end subroutine rate_n13ap_jina
1465 :
1466 : ! rn13gp, n13(g,p)c12
1467 : ! see c12pg
1468 :
1469 : ! n14(p,g)o15
1470 :
1471 0 : subroutine rate_n14pg_jina(tf, temp, fr, rr)
1472 : type (T_Factors) :: tf
1473 : real(dp), intent(in) :: temp
1474 : real(dp), intent(out) :: fr, rr
1475 : ! p n14 o15 im05n 7.29680d+00
1476 0 : call jina_reaclib_2_1(ih1, in14, io15, tf, fr, rr, 'rate_n14pg_jina')
1477 0 : end subroutine rate_n14pg_jina
1478 :
1479 :
1480 0 : subroutine rate_n14pg_fxt(tf, temp, fr, rr)
1481 : ! rn14pg ro15gp
1482 : ! n14(p, g)o15
1483 : type (T_Factors) :: tf
1484 : real(dp), intent(in) :: temp
1485 : real(dp), intent(out) :: fr, rr
1486 0 : real(dp) :: term, rev, aa, bb, cc, dd, ee, q1
1487 : parameter (q1 = 1.0d0/10.850436d0)
1488 :
1489 0 : if (tf% t9 < lowT9_cutoff) then
1490 0 : fr = 0; rr = 0; return
1491 : end if
1492 :
1493 0 : aa = 4.90d+07 * (tf% T9i23) * exp(-15.228d0*(tf% T9i13) - (tf% T92)*q1)
1494 : bb = 1.0d0 + 0.027d0*(tf% T913) - 0.778d0*(tf% T923) - 0.149d0*(tf% T9) &
1495 0 : + 0.261d0*(tf% T943) + 0.127d0*(tf% T953)
1496 0 : cc = aa * bb
1497 0 : dd = 2.37d+03 * (tf% T9i32) * exp(-3.011d0*(tf% T9i))
1498 0 : ee = 2.19d+04 * exp(-12.530d0*(tf% T9i))
1499 0 : term = cc + dd + ee
1500 0 : rev = 2.70d+10 * (tf% T932) * exp(-84.678d0*(tf% T9i))
1501 0 : fr = term
1502 0 : rr = rev * term
1503 : end subroutine rate_n14pg_fxt
1504 :
1505 :
1506 50113 : subroutine rate_n14pg_nacre(tf, temp, fr, rr)
1507 : type (T_Factors) :: tf
1508 : real(dp), intent(in) :: temp
1509 : real(dp), intent(out) :: fr, rr
1510 50113 : real(dp) :: term, rev
1511 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
1512 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
1513 : ! + c0 T9i32 exp(-c1/T9)
1514 : ! + d0 T9i32 exp(-d1/T9)
1515 : ! + e0 T9^e1 exp(-e2/T9)
1516 :
1517 50113 : if (tf% t9 < lowT9_cutoff) then
1518 6990 : fr = 0; rr = 0; return
1519 : end if
1520 :
1521 : call rnacre(tf, &
1522 : 4.83d7, 15.231d0, 1d0/0.8d0, & ! a0, a1, a2
1523 : -2.00d0, 3.41d0, -2.43d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
1524 : 2.36d3, 3.010d0, & ! c0, c1
1525 : 0d0, 0d0, & ! d0, d1
1526 : 6.72d3, 0.380d0, 9.530d0, & ! e0, e1, e2
1527 43123 : term)
1528 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
1529 : 2.699d10, 84.677d0, & ! a0, a1
1530 43123 : rev)
1531 43123 : fr = term
1532 43123 : rr = rev * term
1533 : end subroutine rate_n14pg_nacre
1534 :
1535 : ! rn14ap, n14(a,p)o17
1536 : ! see ro17pa
1537 :
1538 : ! rn14gp, n14(g,p)c13
1539 : ! see rc13pg
1540 :
1541 : ! rn14ag, n14(a,g)f18
1542 :
1543 : ! n14(a,g)f18
1544 20056 : subroutine rate_n14ag_jina(tf, temp, fr, rr)
1545 : type (T_Factors) :: tf
1546 : real(dp), intent(in) :: temp
1547 : real(dp), intent(out) :: fr, rr
1548 : ! he4 n14 f18 ga00r 4.41500d+00
1549 20056 : call jina_reaclib_2_1(ihe4, in14, if18, tf, fr, rr, 'rate_n14ag_jina')
1550 20056 : end subroutine rate_n14ag_jina
1551 :
1552 :
1553 : ! rn15pg, n15(p,g)o16
1554 :
1555 20056 : subroutine rate_n15pg_jina(tf, temp, fr, rr)
1556 : type (T_Factors) :: tf
1557 : real(dp), intent(in) :: temp
1558 : real(dp), intent(out) :: fr, rr
1559 20056 : call jina_reaclib_2_1(ih1, in15, io16, tf, fr, rr, 'rate_n15pg_jina')
1560 20056 : end subroutine rate_n15pg_jina
1561 :
1562 :
1563 : ! rn15pa, n15(p,a)c12
1564 : ! see rc12ap
1565 :
1566 : ! rn15ap, n15(a,p)o18
1567 : ! see ro18pa
1568 :
1569 :
1570 : ! Oxygen
1571 :
1572 :
1573 : ! ro14ap, o14(a,p)f17
1574 :
1575 0 : subroutine rate_o14ap_fxt(tf, temp, fr, rr)
1576 : type (T_Factors) :: tf
1577 : real(dp), intent(in) :: temp
1578 : real(dp), intent(out) :: fr, rr
1579 0 : real(dp) :: term, rev, aa, bb, cc, dd, ee, ff, q1
1580 : parameter (q1 = 1.0d0/0.514089d0)
1581 :
1582 0 : if (tf% t9 < lowT9_cutoff) then
1583 0 : fr = 0; rr = 0; return
1584 : end if
1585 :
1586 0 : aa = 1.68d+13 * (tf% T9i23) * exp(-39.388d0*(tf% T9i13)- (tf% T92)*q1)
1587 : bb = 1.0d0 + 0.011d0*(tf% T913) + 13.117d0*(tf% T923) + 0.971d0*(tf% T9) &
1588 0 : + 85.295d0*(tf% T943) + 16.061d0*(tf% T953)
1589 0 : cc = aa * bb
1590 0 : dd = 3.31d+04 * (tf% T9i32) * exp(-11.733d0*(tf% T9i))
1591 0 : ee = 1.79d+07 * (tf% T9i32) * exp(-22.609d0*(tf% T9i))
1592 0 : ff = 9.00d+03 * (tf% T9113) * exp(-12.517d0*(tf% T9i))
1593 0 : term = cc + dd + ee + ff
1594 0 : fr = term
1595 0 : rev = 4.93d-01*exp(-13.820d0*(tf% T9i))
1596 0 : rr = rev * term
1597 : end subroutine rate_o14ap_fxt
1598 :
1599 :
1600 0 : subroutine rate_o14ap_jina(tf, temp, fr, rr) ! Hahn 1996 PhRvC 54, 4, p1999-2013
1601 : type (T_Factors) :: tf
1602 : real(dp), intent(in) :: temp
1603 : real(dp), intent(out) :: fr, rr
1604 : ! he4 o14 p f17 Ha96r 1.19200d+00
1605 0 : call jina_reaclib_2_2(ihe4, io14, ih1, if17, tf, fr, rr, 'rate_o14ap_jina')
1606 0 : end subroutine rate_o14ap_jina
1607 :
1608 :
1609 : ! ro14ag, o14(a,g)ne18
1610 0 : subroutine rate_o14ag_jina(tf, temp, fr, rr)
1611 : type (T_Factors) :: tf
1612 : real(dp), intent(in) :: temp
1613 : real(dp), intent(out) :: fr, rr
1614 : ! he4 o14 ne18 wh87n 5.11400d+00
1615 0 : call jina_reaclib_2_1(ihe4, io14, ine18, tf, fr, rr, 'rate_o14ag_jina')
1616 0 : end subroutine rate_o14ag_jina
1617 :
1618 :
1619 : ! ro14gp, o14(g,p)n13
1620 : ! see rn13pg
1621 :
1622 :
1623 : ! ro15ap, o15(a,p)f18
1624 : ! see rf18pa
1625 :
1626 : ! ro15ag, o15(a,g)ne19
1627 :
1628 0 : subroutine rate_o15ag_fxt(tf, temp, fr, rr)
1629 : type (T_Factors) :: tf
1630 : real(dp), intent(in) :: temp
1631 : real(dp), intent(out) :: fr, rr
1632 0 : real(dp) :: term, rev, aa, bb, cc, dd, ee, ff, gg, hh, q1, q2, q3
1633 : parameter (q1 = 1.0d0/9.0d0, &
1634 : q2 = 1.0d0/3.751969d0, &
1635 : q3 = 1.0d0/64.0d0)
1636 :
1637 0 : if (tf% t9 < lowT9_cutoff) then
1638 0 : fr = 0; rr = 0; return
1639 : end if
1640 :
1641 0 : aa = 3.57d+11 * (tf% T9i23) * exp(-39.584d+0*(tf% T9i13) - (tf% T92)*q1)
1642 0 : bb = 1.0d0 + 0.011d0*(tf% T913) - 0.273d0*(tf% T923) - 0.020d0*(tf% T9)
1643 0 : cc = aa*bb
1644 0 : dd = 5.10d+10 * (tf% T9i23) * exp(-39.584d+0*(tf% T9i13) - (tf% T92)*q2)
1645 : ee = 1.0d0 + 0.011d0*(tf% T913) + 1.59d0*(tf% T923) + 0.117d0*(tf% T9) &
1646 0 : + 1.81d0*(tf% T943) + 0.338d0*(tf% T953)
1647 0 : ff = dd*ee
1648 0 : gg = 3.95d-1 * (tf% T9i32) * exp(-5.849d0*(tf% T9i))
1649 0 : hh = 1.90d+1 * pow((tf% T9),2.85d0) * exp(-7.356d0*(tf% T9i) - (tf% T92)*q3)
1650 0 : term = cc + ff + gg + hh
1651 0 : fr = term
1652 0 : rev = 5.54d+10 * (tf% T932) * exp(-40.957d0*(tf% T9i))
1653 0 : rr = rev * term
1654 : end subroutine rate_o15ag_fxt
1655 :
1656 :
1657 0 : subroutine rate_o15ag_jina(tf, temp, fr, rr)
1658 : type (T_Factors) :: tf
1659 : real(dp), intent(in) :: temp
1660 : real(dp), intent(out) :: fr, rr
1661 : ! he4 o15 ne19 Ha96n 3.52900d+00
1662 0 : call jina_reaclib_2_1(ihe4, io15, ine19, tf, fr, rr, 'rate_o15ag_jina')
1663 0 : end subroutine rate_o15ag_jina
1664 :
1665 :
1666 : ! ro15gp, o15(g,p)n14
1667 : ! see rn14pg
1668 :
1669 : ! ro16pg, o16(p,g)f17
1670 30075 : subroutine rate_o16pg_nacre(tf, temp, fr, rr)
1671 : type (T_Factors) :: tf
1672 : real(dp), intent(in) :: temp
1673 : real(dp), intent(out) :: fr, rr
1674 30075 : real(dp) :: term, rev, aa, bbm1, bb
1675 :
1676 30075 : if (tf% t9 < lowT9_cutoff) then
1677 4194 : fr = 0; rr = 0; return
1678 : end if
1679 :
1680 25881 : aa = 7.37d7 * pow((tf% T9),-0.82d0) * exp(-16.696d0*(tf% T9i13))
1681 25881 : bbm1 = 202d0 * exp(-70.348d0*(tf% T9i) - 0.161d0*(tf% T9))
1682 25881 : bb = 1 + bbm1
1683 25881 : term = aa * bb
1684 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
1685 : 3.037d9, 6.966d0, & ! a0, a1
1686 25881 : rev)
1687 25881 : fr = term
1688 25881 : rr = rev * term
1689 : end subroutine rate_o16pg_nacre
1690 :
1691 :
1692 0 : subroutine rate_o16pg_jina(tf, temp, fr, rr)
1693 : type (T_Factors) :: tf
1694 : real(dp), intent(in) :: temp
1695 : real(dp), intent(out) :: fr, rr
1696 : ! p o16 f17 nacrn 6.00000d-01
1697 0 : call jina_reaclib_2_1(ih1, io16, if17, tf, fr, rr, 'rate_o16pg_jina')
1698 0 : end subroutine rate_o16pg_jina
1699 :
1700 : ! ro16ap, o16(a,p)f19
1701 : ! see rf19pa
1702 :
1703 : ! ro16ag, o16(a,g)ne20
1704 :
1705 20038 : subroutine rate_o16ag_nacre(tf, temp, fr, rr)
1706 : type (T_Factors) :: tf
1707 : real(dp), intent(in) :: temp
1708 : real(dp), intent(out) :: fr, rr
1709 20038 : real(dp) :: term, rev
1710 :
1711 20038 : if (tf% t9 < lowT9_cutoff) then
1712 2796 : fr = 0; rr = 0; return
1713 : end if
1714 :
1715 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
1716 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
1717 : ! + c0 T9i32 exp(-c1/T9)
1718 : ! + d0 T9i32 exp(-d1/T9)
1719 : ! + e0 T9^e1 exp(-e2/T9)
1720 : call rnacre(tf, &
1721 : 2.68d10, 39.760d0, 1d0/1.6d0, & ! a0, a1, a2
1722 : 0d0, 0d0, 0d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
1723 : 51.1d0, 10.32d0, & ! c0, c1
1724 : 616.1d0, 12.200d0, & ! d0, d1
1725 : 0.41d0, 2.966d0, 11.900d0, & ! e0, e1, e2
1726 17242 : term)
1727 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
1728 : 5.653d10, 54.886d0, & ! a0, a1
1729 17242 : rev)
1730 17242 : fr = term
1731 17242 : rr = rev * term
1732 : end subroutine rate_o16ag_nacre
1733 :
1734 0 : subroutine rate_o16ag_jina(tf, temp, fr, rr) ! jina reaclib -- nacre
1735 : type (T_Factors) :: tf
1736 : real(dp), intent(in) :: temp
1737 : real(dp), intent(out) :: fr, rr
1738 : ! he4 o16 ne20 nacrr 4.73000d+00
1739 0 : call jina_reaclib_2_1(ihe4, io16, ine20, tf, fr, rr, 'rate_o16ag_jina')
1740 0 : end subroutine rate_o16ag_jina
1741 :
1742 :
1743 : ! ro16gp, o16(g,p)n15
1744 : ! see rn15pg
1745 :
1746 : ! ro16ga, o16(g,a)c12
1747 : ! see rc12ag
1748 :
1749 : ! r1616 cf88 fxt
1750 :
1751 18 : subroutine rate_o16o16_fxt(tf, temp, fr, rr)
1752 : type (T_Factors) :: tf
1753 : real(dp), intent(in) :: temp
1754 : real(dp), intent(out) :: fr, rr
1755 18 : real(dp) :: term
1756 :
1757 18 : if (tf% t9 < lowT9_cutoff) then
1758 0 : fr = 0; rr = 0; return
1759 : end if
1760 :
1761 : term = 7.10d36 * (tf% T9i23) * &
1762 : exp(-135.93d0*(tf% T9i13) - 0.629d0*(tf% T923) &
1763 18 : - 0.445d0*(tf% T943) + 0.0103d0*(tf% T9)*(tf% T9))
1764 18 : fr = term
1765 18 : rr = 0.0d0
1766 : end subroutine rate_o16o16_fxt
1767 :
1768 18 : subroutine rate_o16o16g_fxt(tf, temp, fr, rr)
1769 : type (T_Factors) :: tf
1770 : real(dp), intent(in) :: temp
1771 : real(dp), intent(out) :: fr, rr
1772 : !
1773 18 : call rate_o16o16_fxt(tf, temp, fr, rr)
1774 18 : fr = fr*0.10d0
1775 18 : end subroutine rate_o16o16g_fxt
1776 :
1777 : ! r1616n, o16(o16,n)s31
1778 :
1779 :
1780 : ! r1616p, o16(o16, p)p31
1781 10073 : subroutine rate_o16o16p_jina(tf, temp, fr, rr)
1782 : type (T_Factors) :: tf
1783 : real(dp), intent(in) :: temp
1784 : real(dp), intent(out) :: fr, rr
1785 : ! o16 o16 p p31 cf88r 7.67800d+00
1786 10073 : call jina_reaclib_2_2(io16, io16, ih1, ip31, tf, fr, rr, 'rate_o16o16p_jina')
1787 10073 : end subroutine rate_o16o16p_jina
1788 :
1789 :
1790 : ! r1616a, o16(o16, a)si28
1791 10037 : subroutine rate_o16o16a_jina(tf, temp, fr, rr)
1792 : type (T_Factors) :: tf
1793 : real(dp), intent(in) :: temp
1794 : real(dp), intent(out) :: fr, rr
1795 : ! o16 o16 he4 si28 cf88r 9.59300d+00
1796 10037 : call jina_reaclib_2_2(io16, io16, ihe4, isi28, tf, fr, rr, 'rate_o16o16a_jina')
1797 10037 : end subroutine rate_o16o16a_jina
1798 :
1799 0 : subroutine rate_o16o16a(tf, temp, fr, rr)
1800 : type (T_Factors) :: tf
1801 : real(dp), intent(in) :: temp
1802 : real(dp), intent(out) :: fr, rr
1803 : real(dp) :: fr1, rr1, fr2, rr2, fr4, rr4
1804 0 : call rate_o16o16npad(tf, temp, fr1, rr1, fr2, rr2, fr, rr, fr4, rr4)
1805 0 : end subroutine rate_o16o16a
1806 :
1807 0 : subroutine rate_o16o16a_fxt(tf, temp, fr, rr)
1808 : type (T_Factors) :: tf
1809 : real(dp), intent(in) :: temp
1810 : real(dp), intent(out) :: fr, rr
1811 : !
1812 0 : call rate_o16o16_fxt(tf, temp, fr, rr)
1813 0 : fr = fr*0.56d0
1814 0 : end subroutine rate_o16o16a_fxt
1815 :
1816 0 : subroutine rate_o16o16p_fxt(tf, temp, fr, rr)
1817 : type (T_Factors) :: tf
1818 : real(dp), intent(in) :: temp
1819 : real(dp), intent(out) :: fr, rr
1820 : !
1821 0 : call rate_o16o16_fxt(tf, temp, fr, rr)
1822 0 : fr = fr*0.34d0
1823 0 : end subroutine rate_o16o16p_fxt
1824 :
1825 20038 : subroutine rate_o16o16_jina(tf, temp, fr, rr)
1826 : type (T_Factors) :: tf
1827 : real(dp), intent(in) :: temp
1828 : real(dp), intent(out) :: fr, rr
1829 : real(dp) :: fr1, rr1, fr2, rr2
1830 10019 : call rate_o16o16a_jina(tf, temp, fr1, rr1)
1831 10019 : call rate_o16o16p_jina(tf, temp, fr2, rr2)
1832 10019 : fr = fr1 + fr2
1833 10019 : rr = rr1 + rr2
1834 10019 : end subroutine rate_o16o16_jina
1835 :
1836 0 : subroutine rate_o16o16npad(tf, temp, &
1837 : fr1, rr1, & ! o16(o16, n)s31 &
1838 : fr2, rr2, & ! o16(o16, p)p31 &
1839 : fr3, rr3, & ! o16(o16, a)si28 &
1840 : fr4, rr4) ! o16(o16, d)p30
1841 :
1842 : type (T_Factors) :: tf
1843 : real(dp) :: temp, fr1, rr1, fr2, rr2, fr3, rr3, fr4, rr4
1844 :
1845 0 : real(dp) :: term, rev, aa, daa, &
1846 0 : b32n, b32p, b32a, b32d, ezro, dlt, xxt, thrs
1847 :
1848 0 : if (tf% t9 < lowT9_cutoff) then
1849 0 : fr1 = 0; rr1 = 0
1850 0 : fr2 = 0; rr2 = 0
1851 0 : fr3 = 0; rr3 = 0
1852 0 : fr4 = 0; rr4 = 0
1853 0 : return
1854 : end if
1855 :
1856 :
1857 : aa = 7.10d36 * (tf% T9i23) * &
1858 : exp(-135.93d0*(tf% T9i13) - 0.629d0*(tf% T923) &
1859 0 : - 0.445d0*(tf% T943) + 0.0103d0*(tf% T9)*(tf% T9))
1860 :
1861 : daa = -twoth*aa*(tf% T9i) &
1862 : + aa * (oneth*135.93d0*(tf% T9i43) - twoth*0.629d0*(tf% T9i13) &
1863 0 : - fourth*0.445d0*(tf% T913) + 0.0206d0*(tf% T9))
1864 :
1865 :
1866 : ! branching ratios highly uncertain; guessed using fcz 1975
1867 : ! deuteron channel is endoergic. apply error function cut-off.
1868 0 : ezro = 3.9d0*(tf% T923)
1869 0 : dlt = 1.34d0*pow(tf% T9,fivsix)
1870 0 : xxt = 2.0d0*(2.406d0 - ezro)/dlt
1871 0 : call fowthrsh(xxt, thrs)
1872 0 : b32d = 0.05d0*thrs
1873 0 : b32n = 0.1d0
1874 0 : b32a = 0.25d0
1875 0 : b32p = 1.0d0 - b32d - b32a - b32n
1876 :
1877 : ! o16(o16, n)s31
1878 0 : term = aa * b32n
1879 0 : fr1 = term
1880 0 : rev = 0.0d0
1881 0 : if ((tf% T9) > 0.1d0) then
1882 0 : rev = 5.92d0 * exp(-16.8038228d0*(tf% T9i))
1883 : end if
1884 0 : rr1 = rev * term
1885 :
1886 : ! o16(o16, p)p31
1887 0 : term = aa * b32p
1888 0 : fr2 = term
1889 0 : rev = 0.0d0
1890 0 : if ((tf% T9) > 0.1d0) then
1891 0 : rev = 5.92d0*exp(-89.0788286d0*(tf% T9i))
1892 : end if
1893 0 : rr2 = rev * term
1894 :
1895 : ! o16(o16, a)si28
1896 0 : term = aa * b32a
1897 0 : fr3 = term
1898 0 : rev = 0.0d0
1899 0 : if ((tf% T9) > 0.1d0) then
1900 0 : rev = 3.46d0*exp(-111.3137212d0*(tf% T9i))
1901 : end if
1902 0 : rr3 = rev * term
1903 :
1904 : ! o16(o16, d)p30
1905 0 : term = aa * b32d
1906 0 : fr4 = term
1907 0 : rev = 0.0d0
1908 0 : if ((tf% T9) > 0.1d0) then
1909 0 : rev = 0.984d0*exp(27.9908982d0*(tf% T9i))
1910 : end if
1911 0 : rr4 = rev * term
1912 : end subroutine rate_o16o16npad
1913 :
1914 :
1915 0 : subroutine fowthrsh(x, thrs)
1916 :
1917 : ! fowler threshold fudge function.
1918 : ! err func rational (abramowitz p.299)7.1.25 and its derivative
1919 :
1920 : ! declare
1921 0 : real(dp) :: x, thrs, ag, z, z2, t, t2, t3, tt, er, aa
1922 0 : ag = sign(1.0d0, x)
1923 0 : z = abs(x)
1924 0 : z2 = z*z
1925 0 : aa = 1.0d0 + 0.47047d0*z
1926 0 : t = 1.0d0/aa
1927 0 : t2 = t*t
1928 0 : t3 = t2*t
1929 0 : tt = 0.3480242d0*t - 0.0958798d0*t2 + 0.7478556d0*t3
1930 0 : thrs = 0.5d0
1931 0 : if (z /= 0) then
1932 0 : aa = exp(-z2)
1933 0 : er = 1.0d0 - tt * aa
1934 0 : thrs = 0.5d0 * (1.0d0 - ag*er)
1935 : end if
1936 0 : end subroutine fowthrsh
1937 :
1938 :
1939 : ! o17(a,g)ne21
1940 :
1941 :
1942 : ! ro17pa, o17(p,a)n14
1943 :
1944 0 : subroutine rate_o17pa_jina(tf, temp, fr, rr)
1945 : type (T_Factors) :: tf
1946 : real(dp), intent(in) :: temp
1947 : real(dp), intent(out) :: fr, rr
1948 : ! p o17 he4 n14 ct07r 1.19164d+00
1949 0 : call jina_reaclib_2_2(ih1, io17, ihe4, in14, tf, fr, rr, 'rate_o17pa_jina')
1950 0 : end subroutine rate_o17pa_jina
1951 :
1952 :
1953 : ! ro17pg, o17(p,g)f18
1954 :
1955 18 : subroutine rate_o17pg_jina(tf, temp, fr, rr) ! jina reaclib Chafa et al. (2007)
1956 : type (T_Factors) :: tf
1957 : real(dp), intent(in) :: temp
1958 : real(dp), intent(out) :: fr, rr
1959 : ! p o17 f18 ct07n 5.60650d+00
1960 18 : call jina_reaclib_2_1(ih1, io17, if18, tf, fr, rr, 'rate_o17pg_jina')
1961 18 : end subroutine rate_o17pg_jina
1962 :
1963 : ! ro18pa, o18(p,a)n15
1964 :
1965 20038 : subroutine rate_o18pa_nacre(tf, temp, fr, rr)
1966 : type (T_Factors) :: tf
1967 : real(dp), intent(in) :: temp
1968 : real(dp), intent(out) :: fr, rr
1969 20038 : real(dp) :: term, rev, bb, gs
1970 :
1971 20038 : if (tf% t9 < lowT9_cutoff) then
1972 2796 : fr = 0; rr = 0; return
1973 : end if
1974 :
1975 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
1976 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
1977 : ! + c0 T9i32 exp(-c1/T9)
1978 : ! + d0 T9i32 exp(-d1/T9)
1979 : ! + e0 T9^e1 exp(-e2/T9)
1980 : call rnacre(tf, &
1981 : 5.58d11, 16.732d0, 1d0/0.51d0, & ! a0, a1, a2
1982 : 3.2d0, 21.8d0, 0d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
1983 : 9.91d-14, 0.232d0, & ! c0, c1
1984 : 2.58d4, 1.665d0, & ! d0, d1
1985 : 3.24d8, -0.378d0, 6.395d0, & ! e0, e1, e2
1986 17242 : gs)
1987 17242 : bb = 1.968d0 * exp(-25.673d0*(tf% T9i) - 0.083d0*(tf% T9))
1988 17242 : if (bb > 1) then ! guard against rate going negative
1989 0 : bb = 1
1990 : end if
1991 17242 : term = gs * (1 - bb)
1992 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
1993 : 1.660d-1, 46.192d0, & ! a0, a1
1994 17242 : rev)
1995 17242 : fr = term
1996 17242 : rr = rev * term
1997 : end subroutine rate_o18pa_nacre
1998 :
1999 :
2000 0 : subroutine rate_o18pa_jina(tf, temp, fr, rr) ! jina reaclib nacre
2001 : type (T_Factors) :: tf
2002 : real(dp), intent(in) :: temp
2003 : real(dp), intent(out) :: fr, rr
2004 : ! p o18 he4 n15 nacrn 3.98100d+00
2005 0 : call jina_reaclib_2_2(ih1, io18, ihe4, in15, tf, fr, rr, 'rate_o18pa_jina')
2006 0 : end subroutine rate_o18pa_jina
2007 :
2008 :
2009 : ! ro18pg, o18(p,g)f19
2010 :
2011 0 : subroutine rate_o18pg_jina(tf, temp, fr, rr)
2012 : type (T_Factors) :: tf
2013 : real(dp), intent(in) :: temp
2014 : real(dp), intent(out) :: fr, rr
2015 : ! p o18 f19 nacrr 7.99400d+00
2016 0 : call jina_reaclib_2_1(ih1, io18, if19, tf, fr, rr, 'rate_o18pg_jina')
2017 0 : end subroutine rate_o18pg_jina
2018 :
2019 :
2020 : ! ro18ag, o18(a,g)ne22
2021 :
2022 0 : subroutine rate_o18ag_jina(tf, temp, fr, rr)
2023 : type (T_Factors) :: tf
2024 : real(dp), intent(in) :: temp
2025 : real(dp), intent(out) :: fr, rr
2026 : ! he4 o18 ne22 dh03r 9.66900d+00
2027 0 : call jina_reaclib_2_1(ihe4, io18, ine22, tf, fr, rr, 'rate_o18ag_jina')
2028 0 : end subroutine rate_o18ag_jina
2029 :
2030 :
2031 : ! Fluorine
2032 :
2033 :
2034 : ! rf17pa, f17(p,a)o14
2035 : ! see ro14ap
2036 :
2037 :
2038 : ! rf17gp, f17(g,p)o16
2039 : ! see ro16pg
2040 :
2041 :
2042 : ! rf17ap f17(a,p)ne20
2043 0 : subroutine rate_f17ap_jina(tf, temp, fr, rr)
2044 : type (T_Factors) :: tf
2045 : real(dp), intent(in) :: temp
2046 : real(dp), intent(out) :: fr, rr
2047 : ! he4 f17 p ne20 nacr 4.13000d+00
2048 0 : call jina_reaclib_2_2(ihe4, if17, ih1, ine20, tf, fr, rr, 'rate_f17ap_jina')
2049 0 : end subroutine rate_f17ap_jina
2050 :
2051 : ! rf18pa, f18(p,a)o15
2052 :
2053 0 : subroutine rate_f18pa_wk82(tf, temp, fr, rr)
2054 : type (T_Factors) :: tf
2055 : real(dp), intent(in) :: temp
2056 : real(dp), intent(out) :: fr, rr
2057 0 : real(dp) :: term, rev, aa, bb, cc, dd, ee, ff
2058 :
2059 0 : if (tf% t9 < lowT9_cutoff) then
2060 0 : fr = 0; rr = 0; return
2061 : end if
2062 :
2063 0 : aa = 1.66d-10 * (tf% T9i32) * exp(-0.302d0*(tf% T9i))
2064 0 : bb = 1.56d+05 * (tf% T9i32) * exp(-3.84d0*(tf% T9i))
2065 0 : cc = 1.36d+06 * (tf% T9i32) * exp(-5.22d0*(tf% T9i))
2066 0 : dd = 8.1d-05 * (tf% T9i32) * exp(-1.05d0*(tf% T9i))
2067 0 : ee = 8.9d-04 * (tf% T9i32) * exp(-1.51d0*(tf% T9i))
2068 0 : ff = 3.0d+05 * (tf% T9i32) * exp(-4.29d0*(tf% T9i))
2069 0 : term = aa + bb + cc + dd + ee + ff
2070 0 : fr = term
2071 0 : rev = 4.93d-01 * exp(-33.433d0*(tf% T9i))
2072 0 : rr = rev * term
2073 : end subroutine rate_f18pa_wk82
2074 :
2075 :
2076 0 : subroutine rate_f18pa_jina(tf, temp, fr, rr)
2077 : type (T_Factors) :: tf
2078 : real(dp), intent(in) :: temp
2079 : real(dp), intent(out) :: fr, rr
2080 : ! p f18 he4 o15 sh03r 2.88215d+00
2081 0 : call jina_reaclib_2_2(ih1, if18, ihe4, io15, tf, fr, rr, 'rate_f18pa_jina')
2082 0 : end subroutine rate_f18pa_jina
2083 :
2084 :
2085 : ! rf18gp, f18(g,p)o17
2086 : ! see ro17pg
2087 :
2088 : ! rf19pg, f19(p,g)ne20
2089 :
2090 10055 : subroutine rate_f19pg_jina(tf, temp, fr, rr)
2091 : type (T_Factors) :: tf
2092 : real(dp), intent(in) :: temp
2093 : real(dp), intent(out) :: fr, rr
2094 : ! p f19 ne20 cf88r 1.28480d+01
2095 10055 : call jina_reaclib_2_1(ih1, if19, ine20, tf, fr, rr, 'rate_f19pg_jina')
2096 10055 : end subroutine rate_f19pg_jina
2097 :
2098 :
2099 : ! rf19pa, f19(p,a)o16
2100 :
2101 50095 : subroutine rate_f19pa_nacre(tf, temp, fr, rr)
2102 : type (T_Factors) :: tf
2103 : real(dp), intent(in) :: temp
2104 : real(dp), intent(out) :: fr, rr
2105 50095 : real(dp) :: term, rev, bb, dd, gs
2106 :
2107 50095 : if (tf% t9 < lowT9_cutoff) then
2108 6990 : fr = 0; rr = 0; return
2109 : end if
2110 :
2111 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
2112 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
2113 : ! + c0 T9i32 exp(-c1/T9)
2114 : ! + d0 T9i32 exp(-d1/T9)
2115 : ! + e0 T9^e1 exp(-e2/T9)
2116 : call rnacre(tf, &
2117 : 2.62d11, 18.116d0, 1d0/0.185d0, & ! a0, a1, a2
2118 : 6.26d-2, 0.285d0, 4.94d-3, 11.5d0, 7.40d4, & ! b0, b1, b2, b3, b4
2119 : 3.80d6, 3.752d0, & ! c0, c1
2120 : 0d0, 0d0, & ! d0, d1
2121 : 3.27d7, -0.193d0, 6.587d0, & ! e0, e1, e2
2122 43105 : gs)
2123 43105 : dd = 7.30d8 * pow(tf% T9,-0.201d0) * exp(-16.249d0*(tf% T9i))
2124 43105 : gs = gs + dd
2125 43105 : bb = 0.755d0 * exp(-1.755d0*(tf% T9i) - 0.174d0*(tf% T9))
2126 43105 : term = gs * (1 + bb)
2127 : call rnacre_rev(tf, & ! a0 T932 exp(-a1/T9)
2128 : 6.538d-1, 94.154d0, & ! a0, a1
2129 43105 : rev)
2130 43105 : fr = term
2131 43105 : rr = rev * term
2132 : end subroutine rate_f19pa_nacre
2133 :
2134 :
2135 0 : subroutine rate_f19pa_jina(tf, temp, fr, rr) ! jina reaclib
2136 : type (T_Factors) :: tf
2137 : real(dp), intent(in) :: temp
2138 : real(dp), intent(out) :: fr, rr
2139 0 : call jina_reaclib_2_2(ih1, if19, ihe4, io16, tf, fr, rr, 'rate_f19pa_jina')
2140 0 : end subroutine rate_f19pa_jina
2141 :
2142 :
2143 : ! rf19gp, f19(g,p)o18
2144 : ! see ro18pg
2145 :
2146 :
2147 : ! rf19ap, f19(a,p)ne22
2148 0 : subroutine rate_f19ap_cf88(tf, temp, fr, rr)
2149 : type (T_Factors) :: tf
2150 : real(dp), intent(in) :: temp
2151 : real(dp), intent(out) :: fr, rr
2152 0 : real(dp) :: term
2153 0 : if (tf% t9 < lowT9_cutoff) then
2154 0 : fr = 0; rr = 0; return
2155 : end if
2156 : term = 4.50d18*tf% T9i23*exp(-43.467d0*tf% T9i13-pow(tf% T9/0.637d0,2))+ &
2157 0 : 7.98d04*tf% T932*exp(-12.760d0*tf% T9i)
2158 0 : fr = term*6.36d00*exp(-19.439d0*tf% T9i)
2159 0 : rr = 0.0d0
2160 : end subroutine rate_f19ap_cf88
2161 :
2162 0 : subroutine rate_f19ap_jina(tf, temp, fr, rr) ! jina reaclib
2163 : type (T_Factors) :: tf
2164 : real(dp), intent(in) :: temp
2165 : real(dp), intent(out) :: fr, rr
2166 : ! he4 f19 p ne22 cf88r 1.67500d+00
2167 0 : call jina_reaclib_2_2(ihe4, if19, ih1, ine22, tf, fr, rr, 'rate_f19ap_jina')
2168 0 : end subroutine rate_f19ap_jina
2169 :
2170 :
2171 : ! Neon
2172 :
2173 :
2174 : ! rne18ap, ne18(a,p)na21
2175 :
2176 0 : subroutine rate_ne18ap_fxt(tf, temp, fr, rr)
2177 : type (T_Factors) :: tf
2178 : real(dp), intent(in) :: temp
2179 : real(dp), intent(out) :: fr, rr
2180 0 : real(dp) :: term, rev, aa, bb, cc, dd, ee, zz
2181 : real(dp) :: z1, a1, ztot, ared, r, c1, c2, c3, c4
2182 : parameter (z1 = 10.0d0, &
2183 : a1 = 18.0d0, &
2184 : ztot = 2.0d0 * z1, &
2185 : ared = 4.0d0*a1/(4.0d0 + a1), &
2186 : r = 5.1566081196876965d0, &
2187 : c1 = 4.9080044545315392d10, &
2188 : c2 = 4.9592784569936502d-2, &
2189 : c3 = 1.9288564401521285d1, &
2190 : c4 = 4.6477847042196437d1)
2191 :
2192 0 : if (tf% t9 < lowT9_cutoff) then
2193 0 : fr = 0; rr = 0; return
2194 : end if
2195 :
2196 : ! note:
2197 : ! r = 1.09 * a1**oneth + 2.3
2198 : ! c1 = 7.833e9 * 0.31 * ztot**fourth/(ared**fivsix)
2199 : ! c2 = 0.08617 * 0.1215 * sqrt(ared*r**3/ztot)
2200 : ! c3 = 2.0d0 * 0.52495 * sqrt(ared*r*ztot)
2201 : ! c4 = 4.2487 * (ztot**2*ared)**oneth
2202 : ! ne18ap(a, p)na21
2203 : ! was a call to aprate
2204 0 : aa = 1.0d0 + c2*(tf% T9)
2205 0 : zz = c2/aa
2206 0 : bb = pow(aa,fivsix)
2207 0 : cc = (tf% T923) * bb
2208 0 : dd = pow(aa,oneth)
2209 0 : ee = (tf% T9i13) * dd
2210 0 : term = c1*exp(c3 - c4*ee)/cc
2211 0 : fr = term
2212 0 : rev = 0.0d0
2213 0 : rr = 0.0d0
2214 : end subroutine rate_ne18ap_fxt
2215 :
2216 10037 : subroutine rate_ne18ap_jina(tf, temp, fr, rr) ! jina reaclib
2217 : type (T_Factors) :: tf
2218 : real(dp), intent(in) :: temp
2219 : real(dp), intent(out) :: fr, rr
2220 : ! he4 ne18 p na21 GW95r 2.62700d+00
2221 10037 : call jina_reaclib_2_2(ihe4, ine18, ih1, ina21, tf, fr, rr, 'rate_ne18ap_jina')
2222 10037 : end subroutine rate_ne18ap_jina
2223 :
2224 :
2225 : ! rne18ag ne18(a,g)mg22 rath
2226 18 : subroutine rate_ne18ag_jina(tf, temp, fr, rr)
2227 : type (T_Factors) :: tf
2228 : real(dp), intent(in) :: temp
2229 : real(dp), intent(out) :: fr, rr
2230 : ! he4 ne18 mg22 rath 8.14100d+00
2231 18 : call jina_reaclib_2_1(ihe4, ine18, img22, tf, fr, rr, 'rate_ne18ag_jina')
2232 18 : end subroutine rate_ne18ag_jina
2233 :
2234 :
2235 : ! rne18gp, ne18(g,p)f17
2236 : ! see rf17pg
2237 :
2238 :
2239 : ! rne19pg, ne19(p,g)na20
2240 :
2241 0 : subroutine rate_ne19pg_fxt(tf, temp, fr, rr)
2242 : type (T_Factors) :: tf
2243 : real(dp), intent(in) :: temp
2244 : real(dp), intent(out) :: fr, rr
2245 0 : real(dp) :: term, rev, aa, bb, cc, dd, ee, ff, gg, q1
2246 : parameter (q1 = 1.0d0/1.304164d0)
2247 :
2248 0 : if (tf% t9 < lowT9_cutoff) then
2249 0 : fr = 0; rr = 0; return
2250 : end if
2251 :
2252 0 : aa = 1.71d+6 * (tf% T9i23) * exp(-19.431d0*(tf% T9i13))
2253 : bb = 1.0d0 + 0.021d0*(tf% T913) + 0.130d0*(tf% T923) + 1.95d-2*(tf% T9) &
2254 0 : + 3.86d-2*(tf% T943) + 1.47d-02*(tf% T953)
2255 0 : cc = aa*bb
2256 0 : dd = 1.89d+5 * (tf% T9i23) * exp(-19.431d0*(tf% T9i13) - (tf% T92)*q1)
2257 : ee = 1.0d0 + 0.021d0*(tf% T913) + 2.13d0*(tf% T923) + 0.320d0*(tf% T9) &
2258 0 : + 2.80d0*(tf% T943) + 1.07d0*(tf% T953)
2259 0 : ff = dd*ee
2260 0 : gg = 8.45d+3 * (tf% T9i54) * exp(-7.64d0*(tf% T9i))
2261 0 : term = cc + ff + gg
2262 0 : fr = term
2263 0 : rev = 7.39d+09 * (tf% T932) * exp(-25.519d0*(tf% T9i))
2264 0 : rr = rev * term
2265 : end subroutine rate_ne19pg_fxt
2266 :
2267 :
2268 10037 : subroutine rate_ne19pg_jina(tf, temp, fr, rr)
2269 : type (T_Factors) :: tf
2270 : real(dp), intent(in) :: temp
2271 : real(dp), intent(out) :: fr, rr
2272 : ! p ne19 na20 cf88r 2.19900d+00
2273 10037 : call jina_reaclib_2_1(ih1, ine19, ina20, tf, fr, rr, 'rate_ne19pg_jina')
2274 10037 : end subroutine rate_ne19pg_jina
2275 :
2276 :
2277 : ! rne19ga, ne19(g,a)o15
2278 : ! see r016ag
2279 :
2280 : ! rne19gp, ne19(g,p)f18
2281 : ! see rf18pg
2282 :
2283 : ! rne20pg, ne20(p,g)na21
2284 :
2285 10037 : subroutine rate_ne20pg_nacre(tf, temp, fr, rr)
2286 : type (T_Factors) :: tf
2287 : real(dp), intent(in) :: temp
2288 : real(dp), intent(out) :: fr, rr
2289 10037 : real(dp) :: term, rev, aa, bb, gs
2290 :
2291 10037 : if (tf% t9 < lowT9_cutoff) then
2292 1398 : fr = 0; rr = 0; return
2293 : end if
2294 :
2295 8639 : aa = 2.35d7 * pow(tf% T9,-1.84d0) * exp(-19.451d0*(tf% T9i13)) * (1 + 10.80d0*(tf% T9))
2296 8639 : gs = aa
2297 8639 : aa = 18.0d0 * (tf% T9i32) * exp(-4.247d0*(tf% T9i))
2298 8639 : gs = gs + aa
2299 8639 : aa = 9.83d0 * (tf% T9i32) * exp(-4.619d0*(tf% T9i))
2300 8639 : gs = gs + aa
2301 8639 : aa = 6.76d4 * pow(tf% T9,-0.641d0) * exp(-11.922d0*(tf% T9i))
2302 8639 : gs = gs + aa
2303 8639 : bb = 7.929d0 * exp(-20.108d0*(tf% T9i) - 0.327d0*(tf% T9))
2304 8639 : if (bb > 1) then ! guard against rate going negative
2305 0 : bb = 1
2306 : end if
2307 8639 : term = gs * (1 - bb)
2308 8639 : rev = 4.637d9 * (tf% T932) * exp(-28.214d0*(tf% T9i))
2309 8639 : fr = term
2310 8639 : rr = rev * term
2311 : end subroutine rate_ne20pg_nacre
2312 :
2313 :
2314 0 : subroutine rate_ne20pg_jina(tf, temp, fr, rr)
2315 : type (T_Factors) :: tf
2316 : real(dp), intent(in) :: temp
2317 : real(dp), intent(out) :: fr, rr
2318 : ! p ne20 na21 nacrr 2.43100d+00
2319 0 : call jina_reaclib_2_1(ih1, ine20, ina21, tf, fr, rr, 'rate_ne20pg_jina')
2320 0 : end subroutine rate_ne20pg_jina
2321 :
2322 :
2323 : ! ne20(a,p)na23
2324 20038 : subroutine rate_ne20ap_jina(tf, temp, fr, rr)
2325 : type (T_Factors) :: tf
2326 : real(dp), intent(in) :: temp
2327 : real(dp), intent(out) :: fr, rr
2328 : include 'formats'
2329 : ! he4 ne20 p na23 ha04rv -2.37900d+00
2330 20038 : call jina_reaclib_2_2(ih1, ina23, ihe4, ine20, tf, rr, fr, 'rate_ne20ap_jina')
2331 20038 : end subroutine rate_ne20ap_jina
2332 :
2333 : ! rne20ag, ne20(a,g)mg24
2334 :
2335 0 : subroutine rate_ne20ag_jina(tf, temp, fr, rr)
2336 : type (T_Factors) :: tf
2337 : real(dp), intent(in) :: temp
2338 : real(dp), intent(out) :: fr, rr
2339 : ! he4 ne20 mg24 nacrr 9.31600d+00
2340 0 : call jina_reaclib_2_1(ihe4, ine20, img24, tf, fr, rr, 'rate_ne20ag_jina')
2341 0 : end subroutine rate_ne20ag_jina
2342 :
2343 :
2344 : ! rne20ga, ne20(g,a)o16
2345 : ! see ro16ag
2346 :
2347 : ! rne20gp, ne20(g,p)f19
2348 : ! see rf19pg
2349 :
2350 : ! rne22pg, ne22(p,g)na23
2351 :
2352 0 : subroutine rate_ne22pg_jina(tf, temp, fr, rr)
2353 : type (T_Factors) :: tf
2354 : real(dp), intent(in) :: temp
2355 : real(dp), intent(out) :: fr, rr
2356 : ! p ne22 na23 ha01r 8.79400d+00
2357 0 : call jina_reaclib_2_1(ih1, ine22, ina23, tf, fr, rr, 'rate_ne22pg_jina')
2358 0 : end subroutine rate_ne22pg_jina
2359 :
2360 : ! ne22(n,g)ne23
2361 :
2362 0 : subroutine rate_ne22ag_jina(tf, temp, fr, rr)
2363 : type (T_Factors) :: tf
2364 : real(dp), intent(in) :: temp
2365 : real(dp), intent(out) :: fr, rr
2366 : ! he4 ne22 mg26 nacr 1.06150d+01
2367 0 : call jina_reaclib_2_1(ihe4, ine22, img26, tf, fr, rr, 'rate_ne22ag_jina')
2368 0 : end subroutine rate_ne22ag_jina
2369 :
2370 10019 : subroutine rate_na23pa_jina(tf, temp, fr, rr)
2371 : type (T_Factors) :: tf
2372 : real(dp), intent(in) :: temp
2373 : real(dp), intent(out) :: fr, rr
2374 : ! p na23 he4 ne20 ha04n 2.37900d+00
2375 10019 : call jina_reaclib_2_2(ih1, ina23, ihe4, ine20, tf, fr, rr, 'rate_na23pa_jina')
2376 10019 : end subroutine rate_na23pa_jina
2377 :
2378 10055 : subroutine rate_na23pg_jina(tf, temp, fr, rr)
2379 : type (T_Factors) :: tf
2380 : real(dp), intent(in) :: temp
2381 : real(dp), intent(out) :: fr, rr
2382 : rr = 0
2383 : ! p na23 mg24 ha04r 1.16910d+01
2384 10055 : call jina_reaclib_2_1(ih1, ina23, img24, tf, fr, rr, 'rate_na23pg_jina')
2385 10055 : end subroutine rate_na23pg_jina
2386 :
2387 : ! Magnesium
2388 :
2389 : ! rmg24ag, mg24(a,g)si28
2390 :
2391 0 : subroutine rate_mg24ag_fxt(tf, temp, fr, rr)
2392 : type (T_Factors) :: tf
2393 : real(dp), intent(in) :: temp
2394 : real(dp), intent(out) :: fr, rr
2395 0 : real(dp) :: term, aa, bb, cc, dd, ee, ff, gg, hh, hhi, rev, rc121
2396 : parameter (rc121 = 0.1d0)
2397 :
2398 0 : if (tf% t9 < lowT9_cutoff) then
2399 0 : fr = 0; rr = 0; return
2400 : end if
2401 :
2402 0 : aa = 4.78d+01 * (tf% T9i32) * exp(-13.506d0*(tf% T9i))
2403 0 : bb = 2.38d+03 * (tf% T9i32) * exp(-15.218d0*(tf% T9i))
2404 0 : cc = 2.47d+02 * (tf% T932) * exp(-15.147d0*(tf% T9i))
2405 0 : dd = rc121 * 1.72d-09 * (tf% T9i32) * exp(-5.028d0*(tf% T9i))
2406 0 : ee = rc121* 1.25d-03 * (tf% T9i32) * exp(-7.929d0*(tf% T9i))
2407 0 : ff = rc121 * 2.43d+01 * (tf% T9i) * exp(-11.523d0*(tf% T9i))
2408 0 : gg = 5.0d0*exp(-15.882d0*(tf% T9i))
2409 0 : hh = 1.0d0 + gg
2410 0 : hhi = 1.0d0/hh
2411 0 : term = (aa + bb + cc + dd + ee + ff) * hhi
2412 0 : fr = term
2413 0 : rev = 6.27d+10 * (tf% T932) * exp(-115.862d0*(tf% T9i))
2414 0 : rr = rev * term
2415 : end subroutine rate_mg24ag_fxt
2416 :
2417 :
2418 0 : subroutine rate_mg24ag_jina(tf, temp, fr, rr)
2419 : type (T_Factors) :: tf
2420 : real(dp), intent(in) :: temp
2421 : real(dp), intent(out) :: fr, rr
2422 : ! he4 mg24 si28 cf88r 9.98400d+00
2423 0 : call jina_reaclib_2_1(ihe4, img24, isi28, tf, fr, rr, 'rate_mg24ag_jina')
2424 0 : end subroutine rate_mg24ag_jina
2425 :
2426 :
2427 : ! rmg24ga, mg24(g,a)ne20
2428 : ! see rne20ag
2429 :
2430 : ! rmg24ap, mg24(a,p)al27
2431 :
2432 0 : subroutine rate_mg24ap_fxt(tf, temp, fr, rr)
2433 : type (T_Factors) :: tf
2434 : real(dp), intent(in) :: temp
2435 : real(dp), intent(out) :: fr, rr
2436 0 : real(dp) :: term, aa, bb, cc, dd, ee, ff, gg, &
2437 0 : term1, term2, rev, rc148, q1
2438 : parameter (rc148 = 0.1d0, &
2439 : q1 = 1.0d0/0.024649d0)
2440 :
2441 0 : if (tf% t9 < lowT9_cutoff) then
2442 0 : fr = 0; rr = 0; return
2443 : end if
2444 :
2445 0 : aa = 1.10d+08 * (tf% T9i23) * exp(-23.261d0*(tf% T9i13) - (tf% T92)*q1)
2446 : bb = 1.0d0 + 0.018d0*(tf% T913) + 12.85d0*(tf% T923) + 1.61d0*(tf% T9) &
2447 0 : + 89.87d0*(tf% T943) + 28.66d0*(tf% T953)
2448 0 : term1 = aa * bb
2449 0 : aa = 129.0d0 * (tf% T9i32) * exp(-2.517d0*(tf% T9i))
2450 0 : bb = 5660.0d0 * (tf% T972) * exp(-3.421d0*(tf% T9i))
2451 0 : cc = rc148 * 3.89d-08 * (tf% T9i32) * exp(-0.853d0*(tf% T9i))
2452 0 : dd = rc148 * 8.18d-09 * (tf% T9i32) * exp(-1.001d0*(tf% T9i))
2453 0 : term2 = aa + bb + cc + dd
2454 0 : ee = oneth*exp(-9.792d0*(tf% T9i))
2455 0 : ff = twoth * exp(-11.773d0*(tf% T9i))
2456 0 : gg = 1.0d0 + ee + ff
2457 0 : term = (term1 + term2)/gg
2458 0 : rev = 1.81d0 * exp(-18.572d0*(tf% T9i))
2459 0 : fr = rev * term
2460 0 : rr = term
2461 : end subroutine rate_mg24ap_fxt
2462 :
2463 :
2464 54 : subroutine rate_mg24ap_jina(tf, temp, fr, rr)
2465 : type (T_Factors) :: tf
2466 : real(dp), intent(in) :: temp
2467 : real(dp), intent(out) :: fr, rr
2468 : rr = 0
2469 : ! he4 mg24 p al27 il01rv -1.60060d+00
2470 54 : call jina_reaclib_2_2(ih1, ial27, ihe4, img24, tf, rr, fr, 'rate_mg24ap_jina')
2471 54 : end subroutine rate_mg24ap_jina
2472 :
2473 : ! Aluminum
2474 :
2475 : ! ral27pg, al27(p,g)si28
2476 :
2477 0 : subroutine rate_al27pg_c96(tf, temp, fr, rr)
2478 : type (T_Factors) :: tf
2479 : real(dp), intent(in) :: temp
2480 : real(dp), intent(out) :: fr, rr
2481 0 : real(dp) :: term, rev, aa, bb, cc, dd, ee, ff, gg
2482 :
2483 0 : if (tf% t9 < lowT9_cutoff) then
2484 0 : fr = 0; rr = 0; return
2485 : end if
2486 :
2487 0 : aa = 1.32d+09 * (tf% T9i23) * exp(-23.26d0*(tf% T9i13))
2488 0 : bb = 3.22d-10 * (tf% T9i32) * exp(-0.836d0*(tf% T9i))*0.17d0
2489 0 : cc = 1.74d+00 * (tf% T9i32) * exp(-2.269d0*(tf% T9i))
2490 0 : dd = 9.92d+00 * (tf% T9i32) * exp(-2.492d0*(tf% T9i))
2491 0 : ee = 4.29d+01 * (tf% T9i32) * exp(-3.273d0*(tf% T9i))
2492 0 : ff = 1.34d+02 * (tf% T9i32) * exp(-3.654d0*(tf% T9i))
2493 0 : gg = 1.77d+04 * pow(tf% T9, 0.53d0) * exp(-4.588d0*(tf% T9i))
2494 0 : term = aa + bb + cc + dd + ee + ff + gg
2495 0 : fr = term
2496 0 : rev = 1.13d+11 * (tf% T932) * exp(-134.434d0*(tf% T9i))
2497 0 : rr = rev * term
2498 : end subroutine rate_al27pg_c96
2499 :
2500 54 : subroutine rate_al27pg_jina(tf, temp, fr, rr)
2501 : type (T_Factors) :: tf
2502 : real(dp), intent(in) :: temp
2503 : real(dp), intent(out) :: fr, rr
2504 : ! p al27 si28 il01r 1.15860d+01
2505 54 : call jina_reaclib_2_1(ih1, ial27, isi28, tf, fr, rr, 'rate_al27pg_jina')
2506 54 : end subroutine rate_al27pg_jina
2507 :
2508 : ! Silicon
2509 :
2510 : ! rsi28ag, si28(a,g)s32
2511 0 : subroutine rate_si28ag_jina(tf, temp, fr, rr)
2512 : type (T_Factors) :: tf
2513 : real(dp), intent(in) :: temp
2514 : real(dp), intent(out) :: fr, rr
2515 : ! he4 si28 s32 rath 6.94800d+00
2516 0 : call jina_reaclib_2_1(ihe4, isi28, is32, tf, fr, rr, 'rate_si28ag_jina')
2517 : !if (abs(temp - 3.0097470376051402D+09) < 1d2) then
2518 : ! include 'formats'
2519 : ! write(*,1) 'rate_si28ag_jina', fr, rr, temp
2520 : !end if
2521 0 : end subroutine rate_si28ag_jina
2522 :
2523 0 : subroutine rate_si28ag_fxt(tf, temp, fr, rr)
2524 : type (T_Factors) :: tf
2525 : real(dp), intent(in) :: temp
2526 : real(dp), intent(out) :: fr, rr
2527 0 : real(dp) :: term, aa, rev, z, z2, z3
2528 :
2529 : include 'formats'
2530 :
2531 0 : if (tf% t9 < lowT9_cutoff) then
2532 0 : fr = 0; rr = 0; return
2533 : end if
2534 :
2535 0 : z = min((tf% T9), 10.0d0)
2536 0 : z2 = z*z
2537 0 : z3 = z2*z
2538 0 : aa = 1.0d0 + 6.340d-2*z + 2.541d-3*z2 - 2.900d-4*z3
2539 0 : term = 4.82d+22 * (tf% T9i23) * exp(-61.015d0*(tf% T9i13) * aa)
2540 0 : fr = term
2541 0 : rev = 6.461d+10 * (tf% T932) * exp(-80.643d0*(tf% T9i))
2542 0 : rr = rev * term
2543 :
2544 : !if (abs(temp - 3.0097470376051402D+09) < 1d2) then
2545 : ! write(*,1) 'rate_si28ag_fxt', fr, rr, temp
2546 : !end if
2547 :
2548 : end subroutine rate_si28ag_fxt
2549 :
2550 :
2551 : ! rsi28ga, si28(g,a)mg24
2552 : ! see rmg24ag
2553 :
2554 : ! rsi28ap, si28(a,p)p31
2555 :
2556 0 : subroutine rate_si28ap_fxt(tf, temp, fr, rr)
2557 : type (T_Factors) :: tf
2558 : real(dp), intent(in) :: temp
2559 : real(dp), intent(out) :: fr, rr
2560 0 : real(dp) :: term, aa, rev, z, z2, z3
2561 :
2562 0 : if (tf% t9 < lowT9_cutoff) then
2563 0 : fr = 0; rr = 0; return
2564 : end if
2565 :
2566 0 : z = min((tf% T9), 10.0d0)
2567 0 : z2 = z*z
2568 0 : z3 = z2*z
2569 0 : aa = 1.0d0 + 2.798d-3*z + 2.763d-3*z2 - 2.341d-4*z3
2570 0 : term = 4.16d+13 * (tf% T9i23) * exp(-25.631d0*(tf% T9i13) * aa)
2571 0 : rev = 0.5825d0 * exp(-22.224d0*(tf% T9i))
2572 0 : fr = rev * term
2573 0 : rr = term
2574 : end subroutine rate_si28ap_fxt
2575 :
2576 :
2577 54 : subroutine rate_si28ap_jina(tf, temp, fr, rr)
2578 : type (T_Factors) :: tf
2579 : real(dp), intent(in) :: temp
2580 : real(dp), intent(out) :: fr, rr
2581 : ! he4 si28 p p31 il01rv -1.91710d+00
2582 54 : call jina_reaclib_2_2(ih1, ip31, ihe4, isi28, tf, rr, fr, 'rate_si28ap_jina')
2583 54 : end subroutine rate_si28ap_jina
2584 :
2585 :
2586 : ! rsi28gp, si28(g,p)al27
2587 : ! see ral27pg
2588 :
2589 : ! Phosphorus
2590 :
2591 : ! rp31pg, p31(p,g)s32
2592 :
2593 0 : subroutine rate_p31pg_fxt(tf, temp, fr, rr)
2594 : type (T_Factors) :: tf
2595 : real(dp), intent(in) :: temp
2596 : real(dp), intent(out) :: fr, rr
2597 0 : real(dp) :: term, aa, rev, z, z2, z3
2598 : include 'formats'
2599 :
2600 0 : if (tf% t9 < lowT9_cutoff) then
2601 0 : fr = 0; rr = 0; return
2602 : end if
2603 :
2604 0 : z = min((tf% T9), 10.0d0)
2605 0 : z2 = z*z
2606 0 : z3 = z2*z
2607 0 : aa = 1.0d0 + 1.928d-1*z - 1.540d-2*z2 + 6.444d-4*z3
2608 0 : term = 1.08d+16 * (tf% T9i23) * exp(-27.042d0*(tf% T9i13) * aa)
2609 0 : fr = term
2610 0 : rev = 3.764d+10 * (tf% T932) * exp(-102.865d0*(tf% T9i))
2611 0 : rr = rev * term
2612 : end subroutine rate_p31pg_fxt
2613 :
2614 :
2615 54 : subroutine rate_p31pg_jina(tf, temp, fr, rr)
2616 : type (T_Factors) :: tf
2617 : real(dp), intent(in) :: temp
2618 : real(dp), intent(out) :: fr, rr
2619 : ! p p31 s32 il01n 8.86400d+00
2620 54 : call jina_reaclib_2_1(ih1, ip31, is32, tf, fr, rr, 'rate_p31pg_jina')
2621 54 : end subroutine rate_p31pg_jina
2622 :
2623 :
2624 : ! rp31pa, p31(p,a)si28
2625 : ! see rsi28ap
2626 :
2627 :
2628 : ! Sulfur
2629 :
2630 :
2631 : ! rs32ag, s32(a,g)ar36
2632 0 : subroutine rate_s32ag_jina(tf, temp, fr, rr)
2633 : type (T_Factors) :: tf
2634 : real(dp), intent(in) :: temp
2635 : real(dp), intent(out) :: fr, rr
2636 : ! he4 s32 ar36 rath 6.63900d+00
2637 0 : call jina_reaclib_2_1(ihe4, is32, iar36, tf, fr, rr, 'rate_s32ag_jina')
2638 0 : end subroutine rate_s32ag_jina
2639 :
2640 : ! rs32ag, s32(a,g)ar36
2641 :
2642 0 : subroutine rate_s32ag_fxt(tf, temp, fr, rr)
2643 : type (T_Factors) :: tf
2644 : real(dp), intent(in) :: temp
2645 : real(dp), intent(out) :: fr, rr
2646 0 : real(dp) :: term, aa, rev, z, z2, z3
2647 :
2648 0 : if (tf% t9 < lowT9_cutoff) then
2649 0 : fr = 0; rr = 0; return
2650 : end if
2651 :
2652 0 : z = min((tf% T9), 10.0d0)
2653 0 : z2 = z*z
2654 0 : z3 = z2*z
2655 0 : aa = 1.0d0 + 4.913d-2*z + 4.637d-3*z2 - 4.067d-4*z3
2656 0 : term = 1.16d+24 * (tf% T9i23) * exp(-66.690d0*(tf% T9i13) * aa)
2657 0 : fr = term
2658 0 : rev = 6.616d+10 * (tf% T932) * exp(-77.080d0*(tf% T9i))
2659 0 : rr = rev * term
2660 : end subroutine rate_s32ag_fxt
2661 :
2662 :
2663 : ! rs32ga, s32(g,a)si28
2664 : ! see rsi28ag
2665 :
2666 : ! rs32ap, s32(a,p)cl35
2667 :
2668 0 : subroutine rate_s32ap_fxt(tf, temp, fr, rr)
2669 : type (T_Factors) :: tf
2670 : real(dp), intent(in) :: temp
2671 : real(dp), intent(out) :: fr, rr
2672 0 : real(dp) :: term, aa, rev, z, z2, z3
2673 :
2674 0 : if (tf% t9 < lowT9_cutoff) then
2675 0 : fr = 0; rr = 0; return
2676 : end if
2677 :
2678 0 : z = min((tf% T9), 10.0d0)
2679 0 : z2 = z*z
2680 0 : z3 = z2*z
2681 0 : aa = 1.0d0 + 1.041d-1*z - 1.368d-2*z2 + 6.969d-4*z3
2682 0 : term = 1.27d+16 * (tf% T9i23) * exp(-31.044d0*(tf% T9i13) * aa)
2683 0 : rev = 1.144d0 * exp(-21.643d0*(tf% T9i))
2684 0 : fr = rev * term
2685 0 : rr = term
2686 : end subroutine rate_s32ap_fxt
2687 :
2688 :
2689 54 : subroutine rate_s32ap_jina(tf, temp, fr, rr)
2690 : type (T_Factors) :: tf
2691 : real(dp), intent(in) :: temp
2692 : real(dp), intent(out) :: fr, rr
2693 : ! he4 s32 p cl35 il01rv -1.86700d+00
2694 54 : call jina_reaclib_2_2(ih1, icl35, ihe4, is32, tf, rr, fr, 'rate_s32ap_jina')
2695 54 : end subroutine rate_s32ap_jina
2696 :
2697 : ! rs32gp, s32(g,p)p31
2698 : ! see rp31pg
2699 :
2700 :
2701 : ! Chlorine
2702 :
2703 : ! rcl35pg, cl35(p,g)ar36
2704 :
2705 0 : subroutine rate_cl35pg_fxt(tf, temp, fr, rr)
2706 : type (T_Factors) :: tf
2707 : real(dp), intent(in) :: temp
2708 : real(dp), intent(out) :: fr, rr
2709 0 : real(dp) :: term, aa, rev
2710 :
2711 0 : if (tf% t9 < lowT9_cutoff) then
2712 0 : fr = 0; rr = 0; return
2713 : end if
2714 :
2715 0 : aa = 1.0d0 + 1.761d-1*(tf% T9) - 1.322d-2*(tf% T92) + 5.245d-4*(tf% T93)
2716 0 : term = 4.48d+16 * (tf% T9i23) * exp(-29.483d0*(tf% T9i13) * aa)
2717 0 : fr = term
2718 0 : rev = 7.568d+10*(tf% T932)*exp(-98.722d0*(tf% T9i))
2719 0 : rr = rev * term
2720 : end subroutine rate_cl35pg_fxt
2721 :
2722 :
2723 54 : subroutine rate_cl35pg_jina(tf, temp, fr, rr)
2724 : type (T_Factors) :: tf
2725 : real(dp), intent(in) :: temp
2726 : real(dp), intent(out) :: fr, rr
2727 : ! p cl35 ar36 il01r 8.50600d+00
2728 54 : call jina_reaclib_2_1(ih1, icl35, iar36, tf, fr, rr, 'rate_cl35pg_jina')
2729 54 : end subroutine rate_cl35pg_jina
2730 :
2731 : ! rcl35pa, cl35(p,a)s32
2732 : ! see rs32ap
2733 :
2734 : ! Argon
2735 :
2736 :
2737 : ! rar36ag, ar36(a,g)ca40
2738 0 : subroutine rate_ar36ag_jina(tf, temp, fr, rr)
2739 : type (T_Factors) :: tf
2740 : real(dp), intent(in) :: temp
2741 : real(dp), intent(out) :: fr, rr
2742 : ! he4 ar36 ca40 rath 7.04000d+00
2743 0 : call jina_reaclib_2_1(ihe4, iar36, ica40, tf, fr, rr, 'rate_ar36ag_jina')
2744 0 : end subroutine rate_ar36ag_jina
2745 :
2746 : ! rar36ag, ar36(a,g)ca40
2747 :
2748 0 : subroutine rate_ar36ag_fxt(tf, temp, fr, rr)
2749 : type (T_Factors) :: tf
2750 : real(dp), intent(in) :: temp
2751 : real(dp), intent(out) :: fr, rr
2752 0 : real(dp) :: term, aa, rev, z, z2, z3
2753 :
2754 0 : if (tf% t9 < lowT9_cutoff) then
2755 0 : fr = 0; rr = 0; return
2756 : end if
2757 :
2758 0 : z = min((tf% T9), 10.0d0)
2759 0 : z2 = z*z
2760 0 : z3 = z2*z
2761 0 : aa = 1.0d0 + 1.458d-1*z - 1.069d-2*z2 + 3.790d-4*z3
2762 0 : term = 2.81d+30 * (tf% T9i23) * exp(-78.271d0*(tf% T9i13) * aa)
2763 0 : fr = term
2764 0 : rev = 6.740d+10 * (tf% T932) * exp(-81.711d0*(tf% T9i))
2765 0 : rr = rev * term
2766 : end subroutine rate_ar36ag_fxt
2767 :
2768 : ! rar36ga, ar36(g,a)s32
2769 : ! see rs32ag
2770 :
2771 : ! rar36ap, ar36(a,p)k39
2772 :
2773 0 : subroutine rate_ar36ap_fxt(tf, temp, fr, rr)
2774 : type (T_Factors) :: tf
2775 : real(dp), intent(in) :: temp
2776 : real(dp), intent(out) :: fr, rr
2777 0 : real(dp) :: term, aa, rev, z, z2, z3
2778 :
2779 0 : if (tf% t9 < lowT9_cutoff) then
2780 0 : fr = 0; rr = 0; return
2781 : end if
2782 :
2783 0 : z = min((tf% T9), 10.0d0)
2784 0 : z2 = z*z
2785 0 : z3 = z2*z
2786 0 : aa = 1.0d0 + 4.826d-3*z - 5.534d-3*z2 + 4.021d-4*z3
2787 0 : term = 2.76d+13 * (tf% T9i23) * exp(-34.922d0*(tf% T9i13) * aa)
2788 0 : rev = 1.128d0*exp(-14.959d0*(tf% T9i))
2789 0 : fr = rev * term
2790 0 : rr = term
2791 : end subroutine rate_ar36ap_fxt
2792 :
2793 :
2794 54 : subroutine rate_ar36ap_jina(tf, temp, fr, rr)
2795 : type (T_Factors) :: tf
2796 : real(dp), intent(in) :: temp
2797 : real(dp), intent(out) :: fr, rr
2798 : ! he4 ar36 p k39 rath v -1.28800d+00
2799 54 : call jina_reaclib_2_2(ih1, ik39, ihe4, iar36, tf, rr, fr, 'rate_ar36ap_jina')
2800 54 : end subroutine rate_ar36ap_jina
2801 :
2802 : ! rar36gp, ar36(g,p)cl35
2803 : ! see rcl35pg
2804 :
2805 : ! Potassium
2806 :
2807 : ! rk39pg, k39(p,g)ca40
2808 :
2809 0 : subroutine rate_k39pg_fxt(tf, temp, fr, rr)
2810 : type (T_Factors) :: tf
2811 : real(dp), intent(in) :: temp
2812 : real(dp), intent(out) :: fr, rr
2813 0 : real(dp) :: term, aa, rev, z, z2, z3
2814 :
2815 0 : if (tf% t9 < lowT9_cutoff) then
2816 0 : fr = 0; rr = 0; return
2817 : end if
2818 :
2819 0 : z = min((tf% T9), 10.0d0)
2820 0 : z2 = z*z
2821 0 : z3 = z2*z
2822 0 : aa = 1.0d0 + 1.622d-1*z - 1.119d-2*z2 + 3.910d-4*z3
2823 0 : term = 4.09d+16 * (tf% T9i23) * exp(-31.727d0*(tf% T9i13) * aa)
2824 0 : fr = term
2825 0 : rev = 7.600d+10 * (tf% T932) * exp(-96.657d0*(tf% T9i))
2826 0 : rr = rev * term
2827 : end subroutine rate_k39pg_fxt
2828 :
2829 :
2830 54 : subroutine rate_k39pg_jina(tf, temp, fr, rr)
2831 : type (T_Factors) :: tf
2832 : real(dp), intent(in) :: temp
2833 : real(dp), intent(out) :: fr, rr
2834 : ! p k39 ca40 rath 8.32800d+00
2835 54 : call jina_reaclib_2_1(ih1, ik39, ica40, tf, fr, rr, 'rate_k39pg_jina')
2836 54 : end subroutine rate_k39pg_jina
2837 :
2838 : ! rk39pa, k39(p,a)ar36
2839 : ! see rar36ap
2840 :
2841 : ! Calcium
2842 :
2843 :
2844 : ! rca40ag, ca40(a,g)ti44
2845 0 : subroutine rate_ca40ag_jina(tf, temp, fr, rr)
2846 : type (T_Factors) :: tf
2847 : real(dp), intent(in) :: temp
2848 : real(dp), intent(out) :: fr, rr
2849 : include 'formats'
2850 : ! he4 ca40 ti44 rath 5.12700d+00
2851 0 : call jina_reaclib_2_1(ihe4, ica40, iti44, tf, fr, rr, 'rate_ca40ag_jina')
2852 0 : end subroutine rate_ca40ag_jina
2853 :
2854 : ! rca40ag, ca40(a,g)ti44
2855 :
2856 0 : subroutine rate_ca40ag_fxt(tf, temp, fr, rr)
2857 : type (T_Factors) :: tf
2858 : real(dp), intent(in) :: temp
2859 : real(dp), intent(out) :: fr, rr
2860 0 : real(dp) :: term, aa, rev, z, z2, z3
2861 :
2862 0 : if (tf% t9 < lowT9_cutoff) then
2863 0 : fr = 0; rr = 0; return
2864 : end if
2865 :
2866 0 : z = min((tf% T9), 10.0d0)
2867 0 : z2 = z*z
2868 0 : z3 = z2*z
2869 0 : aa = 1.0d0 + 1.650d-2*z + 5.973d-3*z2 - 3.889d-04*z3
2870 0 : term = 4.66d+24 * (tf% T9i23) * exp(-76.435d0*(tf% T9i13) * aa)
2871 0 : fr = term
2872 0 : rev = 6.843d+10 * (tf% T932) * exp(-59.510d0*(tf% T9i))
2873 0 : rr = rev * term
2874 : end subroutine rate_ca40ag_fxt
2875 :
2876 :
2877 : ! rca40ga, ca40(g,a)ar36
2878 : ! see rar36ag
2879 :
2880 : ! rca40ap, ca40(a,p)sc43(p,g)ti44
2881 :
2882 0 : subroutine rate_ca40ap_fxt(tf, temp, fr, rr)
2883 : type (T_Factors) :: tf
2884 : real(dp), intent(in) :: temp
2885 : real(dp), intent(out) :: fr, rr
2886 0 : real(dp) :: term, aa, rev, z, z2, z3
2887 :
2888 0 : if (tf% t9 < lowT9_cutoff) then
2889 0 : fr = 0; rr = 0; return
2890 : end if
2891 :
2892 0 : z = min((tf% T9), 10.0d0)
2893 0 : z2 = z*z
2894 0 : z3 = z2*z
2895 0 : aa = 1.0d0 - 1.206d-2*z + 7.753d-3*z2 - 5.071d-4*z3
2896 0 : term = 4.54d+14 * (tf% T9i23) * exp(-32.177d0*(tf% T9i13) * aa)
2897 0 : rev = 2.229d0 * exp(-40.966d0*(tf% T9i))
2898 0 : fr = rev * term
2899 0 : rr = term
2900 : end subroutine rate_ca40ap_fxt
2901 :
2902 :
2903 54 : subroutine rate_ca40ap_jina(tf, temp, fr, rr)
2904 : type (T_Factors) :: tf
2905 : real(dp), intent(in) :: temp
2906 : real(dp), intent(out) :: fr, rr
2907 : ! he4 ca40 p sc43 rath v -3.52300d+00
2908 54 : call jina_reaclib_2_2(ih1, isc43, ihe4, ica40, tf, rr, fr, 'rate_ca40ap_jina')
2909 54 : end subroutine rate_ca40ap_jina
2910 :
2911 :
2912 : ! rca40ap, ca40(a,p)sc43
2913 : ! see rsc43pa
2914 :
2915 : ! rca40gp, ca40(g,p)k39
2916 : ! see rk39pg
2917 :
2918 : ! Scandium
2919 :
2920 : ! rsc43pg, sc43(p,g)ti44
2921 :
2922 0 : subroutine rate_sc43pg_fxt(tf, temp, fr, rr)
2923 : type (T_Factors) :: tf
2924 : real(dp), intent(in) :: temp
2925 : real(dp), intent(out) :: fr, rr
2926 0 : real(dp) :: term, aa, rev, z, z2, z3
2927 :
2928 0 : if (tf% t9 < lowT9_cutoff) then
2929 0 : fr = 0; rr = 0; return
2930 : end if
2931 :
2932 0 : z = min((tf% T9), 10.0d0)
2933 0 : z2 = z*z
2934 0 : z3 = z2*z
2935 0 : aa = 1.0d0 + 1.023d-1*z - 2.242d-3*z2 - 5.463d-5*z3
2936 0 : term = 3.85d+16 * (tf% T9i23) * exp(-33.234d0*(tf% T9i13) * aa)
2937 0 : fr = term
2938 0 : rev = 1.525d+11 * (tf% T932) * exp(-100.475d0*(tf% T9i))
2939 0 : rr = rev * term
2940 : end subroutine rate_sc43pg_fxt
2941 :
2942 :
2943 54 : subroutine rate_sc43pg_jina(tf, temp, fr, rr)
2944 : type (T_Factors) :: tf
2945 : real(dp), intent(in) :: temp
2946 : real(dp), intent(out) :: fr, rr
2947 : ! p sc43 ti44 rath 8.65000d+00
2948 54 : call jina_reaclib_2_1(ih1, isc43, iti44, tf, fr, rr, 'rate_sc43pg_jina')
2949 54 : end subroutine rate_sc43pg_jina
2950 :
2951 :
2952 : ! rsc43pa, sc43(p,a)ca40
2953 : ! see rca40ap
2954 :
2955 :
2956 : ! Titanium
2957 :
2958 :
2959 : ! rti44ag, ti44(a,g)cr48
2960 0 : subroutine rate_ti44ag_jina(tf, temp, fr, rr)
2961 : type (T_Factors) :: tf
2962 : real(dp), intent(in) :: temp
2963 : real(dp), intent(out) :: fr, rr
2964 : ! he4 ti44 cr48 rath 7.69200d+00
2965 0 : call jina_reaclib_2_1(ihe4, iti44, icr48, tf, fr, rr, 'rate_ti44ag_jina')
2966 0 : end subroutine rate_ti44ag_jina
2967 :
2968 : ! rti44ag, ti44(a,g)cr48
2969 :
2970 0 : subroutine rate_ti44ag_fxt(tf, temp, fr, rr)
2971 : type (T_Factors) :: tf
2972 : real(dp), intent(in) :: temp
2973 : real(dp), intent(out) :: fr, rr
2974 0 : real(dp) :: term, aa, rev, z, z2, z3
2975 :
2976 0 : if (tf% t9 < lowT9_cutoff) then
2977 0 : fr = 0; rr = 0; return
2978 : end if
2979 :
2980 0 : z = min((tf% T9), 10.0d0)
2981 0 : z2 = z*z
2982 0 : z3 = z2*z
2983 0 : aa = 1.0d0 + 1.066d-1*z - 1.102d-2*z2 + 5.324d-4*z3
2984 0 : term = 1.37d+26 * (tf% T9i23) * exp(-81.227d0*(tf% T9i13) * aa)
2985 0 : fr = term
2986 0 : rev = 6.928d+10*(tf% T932)*exp(-89.289d0*(tf% T9i))
2987 0 : rr = rev * term
2988 : end subroutine rate_ti44ag_fxt
2989 :
2990 :
2991 : ! rti44ga, ti44(g,a)ca40
2992 : ! see rca40ag
2993 :
2994 : ! rti44ap, ti44(a,p)v47
2995 :
2996 0 : subroutine rate_ti44ap_fxt(tf, temp, fr, rr)
2997 : type (T_Factors) :: tf
2998 : real(dp), intent(in) :: temp
2999 : real(dp), intent(out) :: fr, rr
3000 0 : real(dp) :: term, aa, rev, z, z2, z3
3001 :
3002 0 : if (tf% t9 < lowT9_cutoff) then
3003 0 : fr = 0; rr = 0; return
3004 : end if
3005 :
3006 0 : z = min((tf% T9), 10.0d0)
3007 0 : z2 = z*z
3008 0 : z3 = z2*z
3009 0 : aa = 1.0d0 + 2.655d-2*z - 3.947d-3*z2 + 2.522d-4*z3
3010 0 : term = 6.54d+20 * (tf% T9i23) * exp(-66.678d0*(tf% T9i13) * aa)
3011 0 : rev = 1.104d0 * exp(-4.723d0*(tf% T9i))
3012 0 : fr = rev * term
3013 0 : rr = term
3014 : end subroutine rate_ti44ap_fxt
3015 :
3016 :
3017 54 : subroutine rate_ti44ap_jina(tf, temp, fr, rr)
3018 : type (T_Factors) :: tf
3019 : real(dp), intent(in) :: temp
3020 : real(dp), intent(out) :: fr, rr
3021 : ! he4 ti44 p v47 chw0r -4.10500d-01
3022 54 : call jina_reaclib_2_2(ihe4, iti44, ih1, iv47, tf, fr, rr, 'rate_ti44ap_jina')
3023 54 : end subroutine rate_ti44ap_jina
3024 :
3025 :
3026 : ! rti44gp, ti44(g,p)sc43
3027 : ! see rsc43pg
3028 :
3029 :
3030 : ! Vanadium
3031 :
3032 : ! rv47pg, v47(p,g)cr48
3033 :
3034 0 : subroutine rate_v47pg_fxt(tf, temp, fr, rr)
3035 : type (T_Factors) :: tf
3036 : real(dp), intent(in) :: temp
3037 : real(dp), intent(out) :: fr, rr
3038 0 : real(dp) :: term, aa, rev, z, z2, z3
3039 :
3040 0 : if (tf% t9 < lowT9_cutoff) then
3041 0 : fr = 0; rr = 0; return
3042 : end if
3043 :
3044 0 : z = min((tf% T9), 10.0d0)
3045 0 : z2 = z*z
3046 0 : z3 = z2*z
3047 0 : aa = 1.0d0 + 9.979d-2*z - 2.269d-3*z2 - 6.662d-5*z3
3048 0 : term = 2.05d+17 * (tf% T9i23) * exp(-35.568d0*(tf% T9i13) * aa)
3049 0 : fr = term
3050 0 : rev = 7.649d+10*(tf% T932)*exp(-93.999d0*(tf% T9i))
3051 0 : rr = rev * term
3052 : end subroutine rate_v47pg_fxt
3053 :
3054 :
3055 54 : subroutine rate_v47pg_jina(tf, temp, fr, rr)
3056 : type (T_Factors) :: tf
3057 : real(dp), intent(in) :: temp
3058 : real(dp), intent(out) :: fr, rr
3059 : ! p v47 cr48 nfisn 8.10607d+00
3060 54 : call jina_reaclib_2_1(ih1, iv47, icr48, tf, fr, rr, 'rate_v47pg_jina')
3061 54 : end subroutine rate_v47pg_jina
3062 :
3063 :
3064 : ! rv47pa, v47(p,a)ti44
3065 : ! see rti44ap
3066 :
3067 :
3068 : ! Chromium
3069 :
3070 :
3071 : ! rcr48ag, cr48(a,g)fe52
3072 0 : subroutine rate_cr48ag_jina(tf, temp, fr, rr)
3073 : type (T_Factors) :: tf
3074 : real(dp), intent(in) :: temp
3075 : real(dp), intent(out) :: fr, rr
3076 : ! he4 cr48 fe52 rath 7.93900d+00
3077 0 : call jina_reaclib_2_1(ihe4, icr48, ife52, tf, fr, rr, 'rate_cr48ag_jina')
3078 0 : end subroutine rate_cr48ag_jina
3079 :
3080 : ! rcr48ag, cr48(a,g)fe52
3081 :
3082 0 : subroutine rate_cr48ag_fxt(tf, temp, fr, rr)
3083 : type (T_Factors) :: tf
3084 : real(dp), intent(in) :: temp
3085 : real(dp), intent(out) :: fr, rr
3086 0 : real(dp) :: term, aa, rev, z, z2, z3
3087 :
3088 0 : if (tf% t9 < lowT9_cutoff) then
3089 0 : fr = 0; rr = 0; return
3090 : end if
3091 :
3092 0 : z = min((tf% T9), 10.0d0)
3093 0 : z2 = z*z
3094 0 : z3 = z2*z
3095 0 : aa = 1.0d0 + 6.325d-2*z - 5.671d-3*z2 + 2.848d-4*z3
3096 0 : term = 1.04d+23 * (tf% T9i23) * exp(-81.420d0*(tf% T9i13) * aa)
3097 0 : fr = term
3098 0 : rev = 7.001d+10 * (tf% T932) * exp(-92.177d0*(tf% T9i))
3099 0 : rr = rev * term
3100 : end subroutine rate_cr48ag_fxt
3101 :
3102 :
3103 : ! rcr48ga, cr48(g,a)ti44
3104 : ! see rti44ag
3105 :
3106 : ! rcr48ap, cr48(a,p)mn51
3107 :
3108 0 : subroutine rate_cr48ap_fxt(tf, temp, fr, rr)
3109 : type (T_Factors) :: tf
3110 : real(dp), intent(in) :: temp
3111 : real(dp), intent(out) :: fr, rr
3112 0 : real(dp) :: term, aa, rev, z, z2, z3
3113 :
3114 0 : if (tf% t9 < lowT9_cutoff) then
3115 0 : fr = 0; rr = 0; return
3116 : end if
3117 :
3118 0 : z = min((tf% T9), 10.0d0)
3119 0 : z2 = z*z
3120 0 : z3 = z2*z
3121 0 : aa = 1.0d0 + 1.384d-2*z + 1.081d-3*z2 - 5.933d-5*z3
3122 0 : term = 1.83d+26 * (tf% T9i23) * exp(-86.741d0*(tf% T9i13) * aa)
3123 0 : fr = term
3124 0 : rev = 0.6087d0*exp(-6.510d0*(tf% T9i))
3125 0 : rr = rev * term
3126 : end subroutine rate_cr48ap_fxt
3127 :
3128 :
3129 54 : subroutine rate_cr48ap_jina(tf, temp, fr, rr)
3130 : type (T_Factors) :: tf
3131 : real(dp), intent(in) :: temp
3132 : real(dp), intent(out) :: fr, rr
3133 : ! he4 cr48 p mn51 rath 5.58000d-01
3134 54 : call jina_reaclib_2_2(ihe4, icr48, ih1, imn51, tf, fr, rr, 'rate_cr48ap_jina')
3135 54 : end subroutine rate_cr48ap_jina
3136 :
3137 :
3138 : ! rcr48gp, cr48(g,p)v47
3139 : ! see rv47pg
3140 :
3141 :
3142 : ! Manganese
3143 :
3144 : ! rmn51pg, mn51(p,g)fe52
3145 :
3146 0 : subroutine rate_mn51pg_fxt(tf, temp, fr, rr)
3147 : type (T_Factors) :: tf
3148 : real(dp), intent(in) :: temp
3149 : real(dp), intent(out) :: fr, rr
3150 0 : real(dp) :: term, aa, rev, z, z2, z3
3151 :
3152 0 : if (tf% t9 < lowT9_cutoff) then
3153 0 : fr = 0; rr = 0; return
3154 : end if
3155 :
3156 0 : z = min((tf% T9), 10.0d0)
3157 0 : z2 = z*z
3158 0 : z3 = z2*z
3159 0 : aa = 1.0d0 + 8.922d-2*z - 1.256d-3*z2 - 9.453d-5*z3
3160 0 : term = 3.77d+17 * (tf% T9i23) * exp(-37.516d0*(tf% T9i13) * aa)
3161 0 : fr = term
3162 0 : rev = 1.150d+11*(tf% T932)*exp(-85.667d0*(tf% T9i))
3163 0 : rr = rev * term
3164 : end subroutine rate_mn51pg_fxt
3165 :
3166 :
3167 54 : subroutine rate_mn51pg_jina(tf, temp, fr, rr)
3168 : type (T_Factors) :: tf
3169 : real(dp), intent(in) :: temp
3170 : real(dp), intent(out) :: fr, rr
3171 : ! p mn51 fe52 rath 7.38100d+00
3172 54 : call jina_reaclib_2_1(ih1, imn51, ife52, tf, fr, rr, 'rate_mn51pg_jina')
3173 54 : end subroutine rate_mn51pg_jina
3174 :
3175 :
3176 : ! rmn51pa, mn51(p,a)cr48
3177 : ! see rcr48ap
3178 :
3179 :
3180 : ! rfe52ag, fe52(a,g)ni56
3181 :
3182 0 : subroutine rate_fe52ag_fxt(tf, temp, fr, rr)
3183 : type (T_Factors) :: tf
3184 : real(dp), intent(in) :: temp
3185 : real(dp), intent(out) :: fr, rr
3186 0 : real(dp) :: term, aa, rev, z, z2, z3
3187 :
3188 0 : if (tf% t9 < lowT9_cutoff) then
3189 0 : fr = 0; rr = 0; return
3190 : end if
3191 :
3192 0 : z = min((tf% T9), 10.0d0)
3193 0 : z2 = z*z
3194 0 : z3 = z2*z
3195 0 : aa = 1.0d0 + 7.846d-2*z - 7.430d-3*z2 + 3.723d-4*z3
3196 0 : term = 1.05d+27 * (tf% T9i23) * exp(-91.674d0*(tf% T9i13) * aa)
3197 0 : fr = term
3198 0 : rev = 7.064d+10*(tf% T932)*exp(-92.850d0*(tf% T9i))
3199 0 : rr = rev * term
3200 : end subroutine rate_fe52ag_fxt
3201 :
3202 :
3203 : ! rfe52ga, fe52(g,a)cr48
3204 : ! see rcr48ag
3205 :
3206 : ! rfe52ap, fe52(a,p)co55
3207 :
3208 0 : subroutine rate_fe52ap_fxt(tf, temp, fr, rr)
3209 : type (T_Factors) :: tf
3210 : real(dp), intent(in) :: temp
3211 : real(dp), intent(out) :: fr, rr
3212 0 : real(dp) :: term, aa, rev, z, z2, z3
3213 :
3214 0 : if (tf% t9 < lowT9_cutoff) then
3215 0 : fr = 0; rr = 0; return
3216 : end if
3217 :
3218 0 : z = min((tf% T9), 10.0d0)
3219 0 : z2 = z*z
3220 0 : z3 = z2*z
3221 0 : aa = 1.0d0 + 1.367d-2*z + 7.428d-4*z2 - 3.050d-5*z3
3222 0 : term = 1.30d+27 * (tf% T9i23) * exp(-91.674d0*(tf% T9i13) * aa)
3223 0 : fr = term
3224 0 : rev = 0.4597d0*exp(-9.470d0*(tf% T9i))
3225 0 : rr = rev * term
3226 : end subroutine rate_fe52ap_fxt
3227 :
3228 :
3229 : ! rfe52gp, fe52(g,p)mn51
3230 : ! see mg51pg
3231 :
3232 :
3233 36 : subroutine rate_fe52ng_jina(tf, temp, fr, rr)
3234 : type (T_Factors) :: tf
3235 : real(dp), intent(in) :: temp
3236 : real(dp), intent(out) :: fr, rr
3237 : ! n fe52 fe53 rath 1.06840d+01
3238 36 : call jina_reaclib_2_1(ineut, ife52, ife53, tf, fr, rr, 'rate_fe52ng_jina')
3239 36 : end subroutine rate_fe52ng_jina
3240 :
3241 :
3242 0 : subroutine rate_fe52ng_fxt(tf, temp, fr, rr)
3243 : type (T_Factors) :: tf
3244 : real(dp), intent(in) :: temp
3245 : real(dp), intent(out) :: fr, rr
3246 0 : real(dp) :: term, rev, tq2
3247 :
3248 0 : if (tf% t9 < lowT9_cutoff) then
3249 0 : fr = 0; rr = 0; return
3250 : end if
3251 :
3252 : ! fe52(n, g)fe53
3253 0 : tq2 = (tf% T9) - 0.348d0
3254 0 : term = 9.604d+05 * exp(-0.0626d0*tq2)
3255 0 : fr = term
3256 0 : rev = 2.43d+09 * (tf% T932) * exp(-123.951d0*(tf% T9i))
3257 0 : rr = rev * term
3258 : end subroutine rate_fe52ng_fxt
3259 :
3260 :
3261 0 : subroutine rate_fe53ng_fxt(tf, temp, fr, rr)
3262 : type (T_Factors) :: tf
3263 : real(dp), intent(in) :: temp
3264 : real(dp), intent(out) :: fr, rr
3265 0 : real(dp) :: term, rev, tq1, tq10, tq2
3266 :
3267 0 : if (tf% t9 < lowT9_cutoff) then
3268 0 : fr = 0; rr = 0; return
3269 : end if
3270 :
3271 : ! fe53(n, g)fe54
3272 0 : tq1 = (tf% T9)/0.348d0
3273 0 : tq10 = pow(tq1, 0.10d0)
3274 0 : tq2 = (tf% T9) - 0.348d0
3275 0 : term = 1.817d+06 * tq10 * exp(-0.06319d0*tq2)
3276 0 : fr = term
3277 0 : rev = 1.56d+11 * (tf% T932) * exp(-155.284d0*(tf% T9i))
3278 0 : rr = rev * term
3279 : end subroutine rate_fe53ng_fxt
3280 :
3281 :
3282 54 : subroutine rate_fe53ng_jina(tf, temp, fr, rr)
3283 : type (T_Factors) :: tf
3284 : real(dp), intent(in) :: temp
3285 : real(dp), intent(out) :: fr, rr
3286 : ! n fe53 fe54 rath 1.33780d+01
3287 54 : call jina_reaclib_2_1(ineut, ife53, ife54, tf, fr, rr, 'rate_fe53ng_jina')
3288 54 : end subroutine rate_fe53ng_jina
3289 :
3290 :
3291 0 : subroutine rate_fe54pg_fxt(tf, temp, fr, rr)
3292 : type (T_Factors) :: tf
3293 : real(dp), intent(in) :: temp
3294 : real(dp), intent(out) :: fr, rr
3295 0 : real(dp) :: term, rev, aa, z, z2, z3
3296 :
3297 0 : if (tf% t9 < lowT9_cutoff) then
3298 0 : fr = 0; rr = 0; return
3299 : end if
3300 :
3301 : ! fe54(p, g)co55
3302 0 : z = min((tf% T9), 10.0d0)
3303 0 : z2 = z*z
3304 0 : z3 = z2*z
3305 0 : aa = 1.0d0 + 9.593d-2*z - 3.445d-3*z2 + 8.594d-5*z3
3306 0 : term = 4.51d+17 * (tf% T9i23) * exp(-38.483d0*(tf% T9i13) * aa)
3307 0 : fr = term
3308 0 : rev = 2.400d+09 * (tf% T932) * exp(-58.605d0*(tf% T9i))
3309 0 : rr = rev * term
3310 : end subroutine rate_fe54pg_fxt
3311 :
3312 :
3313 54 : subroutine rate_fe54a_jina(tf, temp, fr, rr)
3314 : type (T_Factors) :: tf
3315 : real(dp), intent(in) :: temp
3316 : real(dp), intent(out) :: fr, rr
3317 : real(dp) :: fr1, fr2, fr3
3318 : include 'formats'
3319 18 : call rate_fe54ag_jina(tf, temp, fr1, rr)
3320 18 : call rate_fe54an_jina(tf, temp, fr2, rr)
3321 18 : call rate_fe54ap_jina(tf, temp, fr3, rr)
3322 18 : fr = fr1 + fr2 + fr3
3323 18 : end subroutine rate_fe54a_jina
3324 :
3325 72 : subroutine rate_fe54an_jina(tf, temp, fr, rr)
3326 : type (T_Factors) :: tf
3327 : real(dp), intent(in) :: temp
3328 : real(dp), intent(out) :: fr, rr
3329 72 : call jina_reaclib_2_2(ineut, ini57, ihe4, ife54, tf, rr, fr, 'rate_fe54an_jina')
3330 72 : end subroutine rate_fe54an_jina
3331 :
3332 :
3333 36 : subroutine rate_fe54ng_jina(tf, temp, fr, rr)
3334 : type (T_Factors) :: tf
3335 : real(dp), intent(in) :: temp
3336 : real(dp), intent(out) :: fr, rr
3337 36 : call jina_reaclib_2_1(ineut, ife54, ife55, tf, fr, rr, 'rate_fe54ng_jina')
3338 36 : end subroutine rate_fe54ng_jina
3339 :
3340 :
3341 36 : subroutine rate_fe55ng_jina(tf, temp, fr, rr)
3342 : type (T_Factors) :: tf
3343 : real(dp), intent(in) :: temp
3344 : real(dp), intent(out) :: fr, rr
3345 36 : call jina_reaclib_2_1(ineut, ife55, ife56, tf, fr, rr, 'rate_fe55ng_jina')
3346 36 : end subroutine rate_fe55ng_jina
3347 :
3348 :
3349 : ! Cobalt
3350 :
3351 : ! rco55pg, co55(p,g)ni56
3352 :
3353 0 : subroutine rate_co55pg_fxt(tf, temp, fr, rr)
3354 : type (T_Factors) :: tf
3355 : real(dp), intent(in) :: temp
3356 : real(dp), intent(out) :: fr, rr
3357 0 : real(dp) :: term, aa, rev, z, z2, z3
3358 :
3359 0 : if (tf% t9 < lowT9_cutoff) then
3360 0 : fr = 0; rr = 0; return
3361 : end if
3362 :
3363 0 : z = min((tf% T9), 10.0d0)
3364 0 : z2 = z*z
3365 0 : z3 = z2*z
3366 0 : aa = 1.0d0 + 9.894d-2*z - 3.131d-3*z2 - 2.160d-5*z3
3367 0 : term = 1.21d+18 * (tf% T9i23) * exp(-39.604d0*(tf% T9i13) * aa)
3368 0 : fr = term
3369 0 : rev = 1.537d+11*(tf% T932)*exp(-83.382d0*(tf% T9i))
3370 0 : rr = rev * term
3371 : end subroutine rate_co55pg_fxt
3372 :
3373 :
3374 : ! rco55pa, co55(p,a)fe52
3375 : ! see rfe52ap
3376 :
3377 : ! Nickel
3378 :
3379 : ! rni56ga, ni56(g,a)fe52
3380 : ! see rfe52ag
3381 :
3382 : ! rni56gp, ni56(g,p)co55
3383 : ! see rco55pg
3384 :
3385 0 : subroutine rate_v44pg_jina(tf, temp, fr, rr)
3386 : type (T_Factors) :: tf
3387 : real(dp), intent(in) :: temp
3388 : real(dp), intent(out) :: fr, rr
3389 : ! p v44 cr45 rath 3.10000d+00
3390 0 : call jina_reaclib_2_1(ih1, iv44, icr45, tf, fr, rr, 'rate_v44pg_jina')
3391 0 : end subroutine rate_v44pg_jina
3392 :
3393 :
3394 0 : subroutine rate_v45pg_jina(tf, temp, fr, rr)
3395 : type (T_Factors) :: tf
3396 : real(dp), intent(in) :: temp
3397 : real(dp), intent(out) :: fr, rr
3398 : ! p v45 cr46 rath 4.88600d+00
3399 0 : call jina_reaclib_2_1(ih1, iv45, icr46, tf, fr, rr, 'rate_v45pg_jina')
3400 0 : end subroutine rate_v45pg_jina
3401 :
3402 0 : subroutine rate_co53pg_jina(tf, temp, fr, rr)
3403 : type (T_Factors) :: tf
3404 : real(dp), intent(in) :: temp
3405 : real(dp), intent(out) :: fr, rr
3406 : ! p co53 ni54 rath 3.85600d+00
3407 0 : call jina_reaclib_2_1(ih1, ico53, ini54, tf, fr, rr, 'rate_co53pg_jina')
3408 0 : end subroutine rate_co53pg_jina
3409 :
3410 :
3411 0 : subroutine rate_co54pg_jina(tf, temp, fr, rr)
3412 : type (T_Factors) :: tf
3413 : real(dp), intent(in) :: temp
3414 : real(dp), intent(out) :: fr, rr
3415 : ! p co54 ni55 rath 4.61400d+00
3416 0 : call jina_reaclib_2_1(ih1, ico54, ini55, tf, fr, rr, 'rate_co54pg_jina')
3417 0 : end subroutine rate_co54pg_jina
3418 :
3419 0 : subroutine rate_ga62pg_jina(tf, temp, fr, rr)
3420 : type (T_Factors) :: tf
3421 : real(dp), intent(in) :: temp
3422 : real(dp), intent(out) :: fr, rr
3423 : ! p ga62 ge63 nfisn 2.23867d+00
3424 0 : call jina_reaclib_2_1(ih1, iga62, ige63, tf, fr, rr, 'rate_ga62pg_jina')
3425 0 : end subroutine rate_ga62pg_jina
3426 :
3427 :
3428 0 : subroutine rate_ga63pg_jina(tf, temp, fr, rr)
3429 : type (T_Factors) :: tf
3430 : real(dp), intent(in) :: temp
3431 : real(dp), intent(out) :: fr, rr
3432 : ! p ga63 ge64 rath 5.02500d+00
3433 0 : call jina_reaclib_2_1(ih1, iga63, ige64, tf, fr, rr, 'rate_ga63pg_jina')
3434 0 : end subroutine rate_ga63pg_jina
3435 :
3436 : ! ni56
3437 :
3438 0 : subroutine rate_fe52ag_jina(tf, temp, fr, rr)
3439 : type (T_Factors) :: tf
3440 : real(dp), intent(in) :: temp
3441 : real(dp), intent(out) :: fr, rr
3442 0 : call jina_reaclib_2_1(ihe4, ife52, ini56, tf, fr, rr, 'rate_fe52ag_jina')
3443 0 : end subroutine rate_fe52ag_jina
3444 :
3445 :
3446 126 : subroutine rate_fe52ap_jina(tf, temp, fr, rr)
3447 : type (T_Factors) :: tf
3448 : real(dp), intent(in) :: temp
3449 : real(dp), intent(out) :: fr, rr
3450 126 : call jina_reaclib_2_2(ihe4, ife52, ih1, ico55, tf, fr, rr, 'rate_fe52ap_jina')
3451 126 : end subroutine rate_fe52ap_jina
3452 :
3453 :
3454 72 : subroutine rate_co55pg_jina(tf, temp, fr, rr)
3455 : type (T_Factors) :: tf
3456 : real(dp), intent(in) :: temp
3457 : real(dp), intent(out) :: fr, rr
3458 72 : call jina_reaclib_2_1(ih1, ico55, ini56, tf, fr, rr, 'rate_co55pg_jina')
3459 72 : end subroutine rate_co55pg_jina
3460 :
3461 :
3462 36 : subroutine rate_fe54pg_jina(tf, temp, fr, rr)
3463 : type (T_Factors) :: tf
3464 : real(dp), intent(in) :: temp
3465 : real(dp), intent(out) :: fr, rr
3466 36 : call jina_reaclib_2_1(ih1, ife54, ico55, tf, fr, rr, 'rate_fe54pg_jina')
3467 36 : end subroutine rate_fe54pg_jina
3468 :
3469 : ! ni58
3470 :
3471 18 : subroutine rate_fe54ag_jina(tf, temp, fr, rr)
3472 : type (T_Factors) :: tf
3473 : real(dp), intent(in) :: temp
3474 : real(dp), intent(out) :: fr, rr
3475 18 : call jina_reaclib_2_1(ihe4, ife54, ini58, tf, fr, rr, 'rate_fe54ag_jina')
3476 18 : end subroutine rate_fe54ag_jina
3477 :
3478 :
3479 36 : subroutine rate_fe54ap_jina(tf, temp, fr, rr)
3480 : type (T_Factors) :: tf
3481 : real(dp), intent(in) :: temp
3482 : real(dp), intent(out) :: fr, rr
3483 36 : call jina_reaclib_2_2(ihe4, ife54, ih1, ico57, tf, fr, rr, 'rate_fe54ap_jina')
3484 36 : end subroutine rate_fe54ap_jina
3485 :
3486 0 : subroutine rate_fe56pg_jina(tf, temp, fr, rr)
3487 : type (T_Factors) :: tf
3488 : real(dp), intent(in) :: temp
3489 : real(dp), intent(out) :: fr, rr
3490 0 : call jina_reaclib_2_1(ih1, ife56, ico57, tf, fr, rr, 'rate_fe56pg_jina')
3491 0 : end subroutine rate_fe56pg_jina
3492 :
3493 :
3494 18 : subroutine rate_c12_c12_to_h1_na23_jina(tf, temp, fr, rr)
3495 : type (T_Factors) :: tf
3496 : real(dp), intent(in) :: temp
3497 : real(dp), intent(out) :: fr, rr
3498 : call jina_reaclib_2_2( &
3499 18 : ic12, ic12, ih1, ina23, tf, fr, rr, 'rate_c12_c12_to_h1_na23_jina')
3500 18 : end subroutine rate_c12_c12_to_h1_na23_jina
3501 :
3502 :
3503 36 : subroutine rate_he4_ne20_to_c12_c12_jina(tf, temp, fr, rr)
3504 : type (T_Factors) :: tf
3505 : real(dp), intent(in) :: temp
3506 : real(dp), intent(out) :: fr, rr
3507 : call jina_reaclib_2_2( &
3508 36 : ihe4, ine20, ic12, ic12, tf, fr, rr, 'rate_he4_ne20_to_c12_c12_jina')
3509 36 : end subroutine rate_he4_ne20_to_c12_c12_jina
3510 :
3511 :
3512 18 : subroutine rate_he4_mg24_to_c12_o16_jina(tf, temp, fr, rr)
3513 : type (T_Factors) :: tf
3514 : real(dp), intent(in) :: temp
3515 : real(dp), intent(out) :: fr, rr
3516 : call jina_reaclib_2_2( &
3517 18 : ihe4, img24, ic12, io16, tf, fr, rr, 'rate_he4_mg24_to_c12_o16_jina')
3518 18 : end subroutine rate_he4_mg24_to_c12_o16_jina
3519 :
3520 :
3521 60026 : subroutine tfactors(tf, logT_in, temp_in)
3522 : ! sets various popular temperature factors
3523 : ! this routine must be called before any of the rates are called
3524 :
3525 : type (T_Factors) :: tf
3526 : real(dp), intent(in) :: logT_in, temp_in
3527 :
3528 60026 : real(dp) :: logT, temp
3529 :
3530 60026 : logT = max(logT_in, 0d0)
3531 60026 : temp = max(temp_in, 1d0)
3532 :
3533 60026 : tf% lnT9 = (logT - 9)*ln10
3534 60026 : tf% T9 = temp * 1.0d-9
3535 60026 : tf% T92 = tf% T9 * tf% T9
3536 60026 : tf% T93 = tf% T9 * tf% T92
3537 60026 : tf% T94 = tf% T9 * tf% T93
3538 60026 : tf% T95 = tf% T9 * tf% T94
3539 60026 : tf% T96 = tf% T9 * tf% T95
3540 :
3541 60026 : tf% T912 = sqrt(tf% T9)
3542 60026 : tf% T932 = tf% T9 * tf% T912
3543 60026 : tf% T952 = tf% T9 * tf% T932
3544 60026 : tf% T972 = tf% T9 * tf% T952
3545 :
3546 60026 : tf% T913 = pow(tf% T9,oneth)
3547 60026 : tf% T923 = tf% T913 * tf% T913
3548 60026 : tf% T943 = tf% T9 * tf% T913
3549 60026 : tf% T953 = tf% T9 * tf% T923
3550 60026 : tf% T973 = tf% T953 * tf% T923
3551 60026 : tf% T9113 = tf% T973 * tf% T943
3552 :
3553 60026 : tf% T914 = pow(tf% T9, 0.25d0)
3554 60026 : tf% T934 = tf% T914 * tf% T914 * tf% T914
3555 60026 : tf% T954 = tf% T9 * tf% T914
3556 60026 : tf% T974 = tf% T9 * tf% T934
3557 :
3558 60026 : tf% T915 = pow(tf% T9,onefif)
3559 60026 : tf% T935 = tf% T915 * tf% T915 * tf% T915
3560 60026 : tf% T945 = tf% T915 * tf% T935
3561 60026 : tf% T965 = tf% T9 * tf% T915
3562 :
3563 60026 : tf% T916 = pow(tf% T9, onesix)
3564 60026 : tf% T976 = tf% T9 * tf% T916
3565 60026 : tf% T9i76 = 1.0d0 / tf% T976
3566 :
3567 60026 : tf% T917 = pow(tf% T9, onesev)
3568 60026 : tf% T927 = tf% T917 * tf% T917
3569 60026 : tf% T947 = tf% T927 * tf% T927
3570 :
3571 60026 : tf% T918 = sqrt(tf% T914)
3572 60026 : tf% T938 = tf% T918 * tf% T918 * tf% T918
3573 60026 : tf% T958 = tf% T938 * tf% T918 * tf% T918
3574 :
3575 60026 : tf% T9i = 1.0d0 / tf% T9
3576 60026 : tf% T9i2 = tf% T9i * tf% T9i
3577 60026 : tf% T9i3 = tf% T9i2 * tf% T9i
3578 :
3579 60026 : tf% T9i12 = 1.0d0 / tf% T912
3580 60026 : tf% T9i32 = tf% T9i * tf% T9i12
3581 60026 : tf% T9i52 = tf% T9i * tf% T9i32
3582 60026 : tf% T9i72 = tf% T9i * tf% T9i52
3583 :
3584 60026 : tf% T9i13 = 1.0d0 / tf% T913
3585 60026 : tf% T9i23 = tf% T9i13 * tf% T9i13
3586 60026 : tf% T9i43 = tf% T9i * tf% T9i13
3587 60026 : tf% T9i53 = tf% T9i * tf% T9i23
3588 :
3589 60026 : tf% T9i14 = 1.0d0 / tf% T914
3590 60026 : tf% T9i34 = tf% T9i14 * tf% T9i14 * tf% T9i14
3591 60026 : tf% T9i54 = tf% T9i * tf% T9i14
3592 :
3593 60026 : tf% T9i15 = 1.0d0 / tf% T915
3594 60026 : tf% T9i35 = tf% T9i15 * tf% T9i15 * tf% T9i15
3595 60026 : tf% T9i45 = tf% T9i15 * tf% T9i35
3596 60026 : tf% T9i65 = tf% T9i * tf% T9i15
3597 :
3598 60026 : tf% T9i17 = 1.0d0 / tf% T917
3599 60026 : tf% T9i27 = tf% T9i17 * tf% T9i17
3600 60026 : tf% T9i47 = tf% T9i27 * tf% T9i27
3601 :
3602 60026 : tf% T9i18 = 1.0d0 / tf% T918
3603 60026 : tf% T9i38 = tf% T9i18 * tf% T9i18 * tf% T9i18
3604 60026 : tf% T9i58 = tf% T9i38 * tf% T9i18 * tf% T9i18
3605 :
3606 60026 : end subroutine tfactors
3607 :
3608 :
3609 0 : subroutine show_nacre_terms( &
3610 : tf, a0, a1, a2, b0, b1, b2, b3, b4, c0, c1, d0, d1, e0, e1, e2)
3611 : type (T_Factors) :: tf
3612 : real(dp), intent(in) :: a0, a1, a2, b0, b1, b2, b3, b4, &
3613 : c0, c1, d0, d1, e0, e1, e2
3614 :
3615 : include 'formats'
3616 : real(dp) :: aa, bb, cc, dd, ee
3617 0 : aa = a0 * (tf% T9i23) * exp(-a1*(tf% T9i13) - (tf% T92)*(a2*a2))
3618 0 : bb = 1 + b0*(tf% T9) + b1*(tf% T92) + b2*(tf% T93) + b3*(tf% T94) + b4*(tf% T95)
3619 0 : cc = c0 * (tf% T9i32) * exp(-c1*(tf% T9i))
3620 0 : dd = d0 * (tf% T9i32) * exp(-d1*(tf% T9i))
3621 0 : ee = e0 * pow(tf% T9,e1) * exp(-e2*(tf% T9i))
3622 :
3623 0 : write(*,1) 'aa', aa
3624 0 : write(*,1) 'bb', bb
3625 0 : write(*,1) 'aa*bb', aa*bb
3626 0 : write(*,1) 'cc', cc
3627 0 : write(*,1) 'dd', dd
3628 0 : write(*,1) 'ee', ee
3629 :
3630 0 : end subroutine show_nacre_terms
3631 :
3632 :
3633 301359 : subroutine rnacre( &
3634 : ! a0 T9i23 exp(-a1 T9i13 - (T9*a2)^2)
3635 : ! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
3636 : ! + c0 T9i32 exp(-c1/T9)
3637 : ! + d0 T9i32 exp(-d1/T9)
3638 : ! + e0 T9^e1 exp(-e2/T9) &
3639 : tf, a0, a1, a2, b0, b1, b2, b3, b4, c0, c1, d0, d1, e0, e1, e2, term)
3640 : type (T_Factors) :: tf
3641 : real(dp), intent(in) :: a0, a1, a2, b0, b1, b2, b3, b4, &
3642 : c0, c1, d0, d1, e0, e1, e2
3643 : real(dp), intent(out) :: term
3644 301359 : real(dp) :: aa, bb, cc, dd, ee
3645 301359 : aa = a0 * (tf% T9i23) * exp(-a1*(tf% T9i13) - (tf% T92)*(a2*a2))
3646 301359 : bb = 1 + b0*(tf% T9) + b1*(tf% T92) + b2*(tf% T93) + b3*(tf% T94) + b4*(tf% T95)
3647 301359 : if (bb < 0) then
3648 13270 : bb = 0
3649 : end if
3650 301359 : cc = c0 * (tf% T9i32) * exp(-c1*(tf% T9i))
3651 301359 : dd = d0 * (tf% T9i32) * exp(-d1*(tf% T9i))
3652 301359 : ee = e0 * pow(tf% T9,e1) * exp(-e2*(tf% T9i))
3653 301359 : term = aa * bb + cc + dd + ee
3654 301359 : end subroutine rnacre
3655 :
3656 :
3657 307202 : subroutine rnacre_rev(tf, a0, a1, rev) ! a0 T932 exp(-a1/T9)
3658 : real(dp), intent(in) :: a0, a1
3659 : real(dp), intent(out) :: rev
3660 : type (T_Factors) :: tf
3661 307202 : rev = a0 * (tf% T932) * exp(-a1*(tf% T9i))
3662 307202 : end subroutine rnacre_rev
3663 :
3664 :
3665 0 : subroutine jina_reaclib_1_1(i1, o1, tf, fr, rr, str)
3666 : integer, intent(in) :: i1, o1
3667 : type (T_Factors) :: tf
3668 : real(dp), intent(out) :: fr, rr
3669 : character (len=*), intent(in) :: str
3670 : integer :: ierr
3671 0 : call try1_reaclib_1_1(i1, o1, tf, fr, rr, str, ierr)
3672 0 : if (ierr == 0) return
3673 : ierr = 0
3674 0 : call try1_reaclib_1_1(o1, i1, tf, rr, fr, str, ierr)
3675 0 : if (ierr /= 0) then
3676 0 : write(*,*) 'failed in jina_reaclib_1_1 ' // trim(str), tf% T9
3677 0 : call mesa_error(__FILE__,__LINE__)
3678 : end if
3679 : end subroutine jina_reaclib_1_1
3680 :
3681 :
3682 0 : subroutine jina_reaclib_1_2(i1, o1, o2, tf, fr, rr, str)
3683 : integer, intent(in) :: i1, o1, o2
3684 : type (T_Factors) :: tf
3685 : real(dp), intent(out) :: fr, rr
3686 : character (len=*), intent(in) :: str
3687 : integer :: ierr
3688 0 : call try1_reaclib_1_2(i1, o1, o2, tf, fr, rr, str, ierr)
3689 0 : if (ierr == 0) return
3690 : ierr = 0
3691 0 : call try1_reaclib_2_1(o1, o2, i1, tf, rr, fr, str, ierr)
3692 0 : if (ierr /= 0) then
3693 0 : write(*,*) 'failed in jina_reaclib_1_2 ' // trim(str), tf% T9
3694 0 : call mesa_error(__FILE__,__LINE__)
3695 : end if
3696 : end subroutine jina_reaclib_1_2
3697 :
3698 :
3699 0 : subroutine jina_reaclib_1_3(i1, o1, o2, o3, tf, fr, rr, str)
3700 : integer, intent(in) :: i1, o1, o2, o3
3701 : type (T_Factors) :: tf
3702 : real(dp), intent(out) :: fr, rr
3703 : character (len=*), intent(in) :: str
3704 : integer :: ierr
3705 0 : call try1_reaclib_1_3(i1, o1, o2, o3, tf, fr, rr, str, ierr)
3706 0 : if (ierr == 0) return
3707 : ierr = 0
3708 0 : call try1_reaclib_3_1(o1, o2, o3, i1, tf, rr, fr, str, ierr)
3709 0 : if (ierr /= 0) then
3710 0 : write(*,*) 'failed in jina_reaclib_1_3 ' // trim(str), tf% T9
3711 0 : call mesa_error(__FILE__,__LINE__)
3712 : end if
3713 : end subroutine jina_reaclib_1_3
3714 :
3715 :
3716 0 : subroutine jina_reaclib_1_4(i1, o1, o2, o3, o4, tf, fr, rr, str)
3717 : integer, intent(in) :: i1, o1, o2, o3, o4
3718 : type (T_Factors) :: tf
3719 : real(dp), intent(out) :: fr, rr
3720 : character (len=*), intent(in) :: str
3721 : integer :: ierr
3722 0 : call try1_reaclib_1_4(i1, o1, o2, o3, o4, tf, fr, rr, str, ierr)
3723 : !if (ierr == 0) return
3724 : !ierr = 0
3725 : !call try1_reaclib_4_1(o1, o2, o3, o4, i1, tf, rr, fr, str, ierr)
3726 0 : if (ierr /= 0) then
3727 0 : write(*,*) 'failed in jina_reaclib_1_4 ' // trim(str), tf% T9
3728 0 : call mesa_error(__FILE__,__LINE__)
3729 : end if
3730 0 : end subroutine jina_reaclib_1_4
3731 :
3732 :
3733 90999 : subroutine jina_reaclib_2_1(i1, i2, o1, tf, fr, rr, str)
3734 : integer, intent(in) :: i1, i2, o1
3735 : type (T_Factors) :: tf
3736 : real(dp), intent(out) :: fr, rr
3737 : character (len=*), intent(in) :: str
3738 : integer :: ierr
3739 90999 : call try1_reaclib_2_1(i1, i2, o1, tf, fr, rr, str, ierr)
3740 90999 : if (ierr == 0) return
3741 : ierr = 0
3742 0 : call try1_reaclib_1_2(o1, i1, i2, tf, rr, fr, str, ierr)
3743 0 : if (ierr /= 0) then
3744 0 : write(*,*) 'failed in jina_reaclib_2_1 ' // trim(str), tf% T9
3745 0 : call mesa_error(__FILE__,__LINE__)
3746 : end if
3747 : end subroutine jina_reaclib_2_1
3748 :
3749 :
3750 101144 : subroutine jina_reaclib_2_2(i1, i2, o1, o2, tf, fr, rr, str)
3751 : integer, intent(in) :: i1, i2, o1, o2
3752 : type (T_Factors) :: tf
3753 : real(dp), intent(out) :: fr, rr
3754 : character (len=*), intent(in) :: str
3755 : integer :: ierr
3756 101000 : call try1_reaclib_2_2(i1, i2, o1, o2, tf, fr, rr, str, ierr)
3757 101000 : if (ierr == 0) return
3758 : ierr = 0
3759 144 : call try1_reaclib_2_2(o1, o2, i1, i2, tf, rr, fr, str, ierr)
3760 144 : if (ierr /= 0) then
3761 0 : write(*,*) 'failed in jina_reaclib_2_2 ' // trim(str), tf% T9
3762 0 : call mesa_error(__FILE__,__LINE__)
3763 : end if
3764 : end subroutine jina_reaclib_2_2
3765 :
3766 :
3767 10019 : subroutine jina_reaclib_2_3(i1, i2, o1, o2, o3, tf, fr, rr, str)
3768 : integer, intent(in) :: i1, i2, o1, o2, o3
3769 : type (T_Factors) :: tf
3770 : real(dp), intent(out) :: fr, rr
3771 : character (len=*), intent(in) :: str
3772 : integer :: ierr
3773 10019 : call try1_reaclib_2_3(i1, i2, o1, o2, o3, tf, fr, rr, str, ierr)
3774 10019 : if (ierr == 0) return
3775 : ierr = 0
3776 0 : call try1_reaclib_3_2(o1, o2, o3, i1, i2, tf, rr, fr, str, ierr)
3777 0 : if (ierr /= 0) then
3778 0 : write(*,*) 'failed in jina_reaclib_3_2 ' // trim(str), tf% T9
3779 0 : call mesa_error(__FILE__,__LINE__)
3780 : end if
3781 : end subroutine jina_reaclib_2_3
3782 :
3783 :
3784 10019 : subroutine jina_reaclib_2_4(i1, i2, o1, o2, o3, o4, tf, fr, rr, str)
3785 : integer, intent(in) :: i1, i2, o1, o2, o3, o4
3786 : type (T_Factors) :: tf
3787 : real(dp), intent(out) :: fr, rr
3788 : character (len=*), intent(in) :: str
3789 : integer :: ierr
3790 10019 : call try1_reaclib_2_4(i1, i2, o1, o2, o3, o4, tf, fr, rr, str, ierr)
3791 10019 : if (ierr == 0) return
3792 : ierr = 0
3793 0 : call try1_reaclib_4_2(o1, o2, o3, o4, i1, i2, tf, rr, fr, str, ierr)
3794 0 : if (ierr /= 0) then
3795 0 : write(*,*) 'failed in jina_reaclib_2_4 ' // trim(str), tf% T9
3796 0 : call mesa_error(__FILE__,__LINE__)
3797 : end if
3798 : end subroutine jina_reaclib_2_4
3799 :
3800 :
3801 20038 : subroutine jina_reaclib_3_1(i1, i2, i3, o1, tf, fr, rr, str)
3802 : integer, intent(in) :: i1, i2, i3, o1
3803 : type (T_Factors) :: tf
3804 : real(dp), intent(out) :: fr, rr
3805 : character (len=*), intent(in) :: str
3806 : integer :: ierr
3807 20038 : call try1_reaclib_3_1(i1, i2, i3, o1, tf, fr, rr, str, ierr)
3808 20038 : if (ierr == 0) return
3809 : ierr = 0
3810 0 : call try1_reaclib_3_1(o1, i1, i2, i3, tf, rr, fr, str, ierr)
3811 0 : if (ierr /= 0) then
3812 0 : write(*,*) 'failed in jina_reaclib_3_1 ' // trim(str), tf% T9
3813 0 : call mesa_error(__FILE__,__LINE__)
3814 : end if
3815 : end subroutine jina_reaclib_3_1
3816 :
3817 :
3818 0 : subroutine jina_reaclib_3_2(i1, i2, i3, o1, o2, tf, fr, rr, str)
3819 : integer, intent(in) :: i1, i2, i3, o1, o2
3820 : type (T_Factors) :: tf
3821 : real(dp), intent(out) :: fr, rr
3822 : character (len=*), intent(in) :: str
3823 : integer :: ierr
3824 0 : call try1_reaclib_3_2(i1, i2, i3, o1, o2, tf, fr, rr, str, ierr)
3825 0 : if (ierr == 0) return
3826 : ierr = 0
3827 0 : call try1_reaclib_2_3(o1, o2, i1, i2, i3, tf, rr, fr, str, ierr)
3828 0 : if (ierr /= 0) then
3829 0 : write(*,*) 'failed in jina_reaclib_3_2 ' // trim(str), tf% T9
3830 0 : call mesa_error(__FILE__,__LINE__)
3831 : end if
3832 : end subroutine jina_reaclib_3_2
3833 :
3834 :
3835 0 : subroutine jina_reaclib_4_2(i1, i2, i3, i4, o1, o2, tf, fr, rr, str)
3836 : integer, intent(in) :: i1, i2, i3, i4, o1, o2
3837 : type (T_Factors) :: tf
3838 : real(dp), intent(out) :: fr, rr
3839 : character (len=*), intent(in) :: str
3840 : integer :: ierr
3841 0 : call try1_reaclib_4_2(i1, i2, i3, i4, o1, o2, tf, fr, rr, str, ierr)
3842 0 : if (ierr == 0) return
3843 : ierr = 0
3844 0 : call try1_reaclib_2_4(o1, o2, i1, i2, i3, i4, tf, rr, fr, str, ierr)
3845 0 : if (ierr /= 0) then
3846 0 : write(*,*) 'failed in jina_reaclib_4_2 ' // trim(str), tf% T9
3847 0 : call mesa_error(__FILE__,__LINE__)
3848 : end if
3849 : end subroutine jina_reaclib_4_2
3850 :
3851 :
3852 0 : subroutine try1_reaclib_1_1(i1, o1, tf, fr, rr, str, ierr)
3853 : integer, intent(in) :: i1, o1
3854 : type (T_Factors) :: tf
3855 : real(dp), intent(out) :: fr, rr
3856 : character (len=*), intent(in) :: str
3857 : integer, intent(out) :: ierr
3858 : integer :: num_in, num_out
3859 : integer, dimension(4) :: nuclides_in, nuclides_out
3860 : logical, parameter :: dbg = .false.
3861 : include 'formats'
3862 : ierr = 0
3863 0 : num_in = 1
3864 0 : nuclides_in(1) = i1
3865 0 : num_out = 1
3866 0 : nuclides_out(1) = o1
3867 : call reaclib_rate( &
3868 : str, num_in, nuclides_in, num_out, nuclides_out, tf% T9, &
3869 0 : fr, rr, ierr)
3870 0 : end subroutine try1_reaclib_1_1
3871 :
3872 :
3873 0 : subroutine try1_reaclib_1_2(i1, o1, o2, tf, fr, rr, str, ierr)
3874 : integer, intent(in) :: i1, o1, o2
3875 : type (T_Factors) :: tf
3876 : real(dp), intent(out) :: fr, rr
3877 : character (len=*), intent(in) :: str
3878 : integer, intent(out) :: ierr
3879 : integer :: num_in, num_out
3880 : integer, dimension(4) :: nuclides_in, nuclides_out
3881 : logical, parameter :: dbg = .false.
3882 : include 'formats'
3883 : ierr = 0
3884 0 : num_in = 1
3885 0 : nuclides_in(1) = i1
3886 0 : num_out = 2
3887 0 : nuclides_out(1) = o1
3888 0 : nuclides_out(2) = o2
3889 : call reaclib_rate( &
3890 : str, num_in, nuclides_in, num_out, nuclides_out, tf% T9, &
3891 0 : fr, rr, ierr)
3892 0 : end subroutine try1_reaclib_1_2
3893 :
3894 :
3895 0 : subroutine try1_reaclib_1_3(i1, o1, o2, o3, tf, fr, rr, str, ierr)
3896 : integer, intent(in) :: i1, o1, o2, o3
3897 : type (T_Factors) :: tf
3898 : real(dp), intent(out) :: fr, rr
3899 : character (len=*), intent(in) :: str
3900 : integer, intent(out) :: ierr
3901 : integer :: num_in, num_out
3902 : integer, dimension(4) :: nuclides_in, nuclides_out
3903 : logical, parameter :: dbg = .false.
3904 : include 'formats'
3905 : ierr = 0
3906 0 : num_in = 1
3907 0 : nuclides_in(1) = i1
3908 0 : num_out = 3
3909 0 : nuclides_out(1) = o1
3910 0 : nuclides_out(2) = o2
3911 0 : nuclides_out(3) = o3
3912 : call reaclib_rate( &
3913 : str, num_in, nuclides_in, num_out, nuclides_out, tf% T9, &
3914 0 : fr, rr, ierr)
3915 0 : end subroutine try1_reaclib_1_3
3916 :
3917 :
3918 0 : subroutine try1_reaclib_1_4(i1, o1, o2, o3, o4, tf, fr, rr, str, ierr)
3919 : integer, intent(in) :: i1, o1, o2, o3, o4
3920 : type (T_Factors) :: tf
3921 : real(dp), intent(out) :: fr, rr
3922 : character (len=*), intent(in) :: str
3923 : integer, intent(out) :: ierr
3924 : integer :: num_in, num_out
3925 : integer, dimension(4) :: nuclides_in, nuclides_out
3926 : logical, parameter :: dbg = .false.
3927 : include 'formats'
3928 : ierr = 0
3929 0 : num_in = 1
3930 0 : nuclides_in(1) = i1
3931 0 : num_out = 4
3932 0 : nuclides_out(1) = o1
3933 0 : nuclides_out(2) = o2
3934 0 : nuclides_out(3) = o3
3935 0 : nuclides_out(4) = o4
3936 : call reaclib_rate( &
3937 : str, num_in, nuclides_in, num_out, nuclides_out, tf% T9, &
3938 0 : fr, rr, ierr)
3939 0 : end subroutine try1_reaclib_1_4
3940 :
3941 :
3942 90999 : subroutine try1_reaclib_2_1(i1, i2, o1, tf, fr, rr, str, ierr)
3943 : integer, intent(in) :: i1, i2, o1
3944 : type (T_Factors) :: tf
3945 : real(dp), intent(out) :: fr, rr
3946 : character (len=*), intent(in) :: str
3947 : integer, intent(out) :: ierr
3948 : integer :: num_in, num_out
3949 : integer, dimension(4) :: nuclides_in, nuclides_out
3950 : logical, parameter :: dbg = .false.
3951 : include 'formats'
3952 : ierr = 0
3953 90999 : num_in = 2
3954 90999 : nuclides_in(1) = i1
3955 90999 : nuclides_in(2) = i2
3956 90999 : num_out = 1
3957 90999 : nuclides_out(1) = o1
3958 : call reaclib_rate( &
3959 : str, num_in, nuclides_in, num_out, nuclides_out, tf% T9, &
3960 90999 : fr, rr, ierr)
3961 90999 : end subroutine try1_reaclib_2_1
3962 :
3963 :
3964 101144 : subroutine try1_reaclib_2_2(i1, i2, o1, o2, tf, fr, rr, str, ierr)
3965 : integer, intent(in) :: i1, i2, o1, o2
3966 : type (T_Factors) :: tf
3967 : real(dp), intent(out) :: fr, rr
3968 : character (len=*), intent(in) :: str
3969 : integer, intent(out) :: ierr
3970 : integer :: num_in, num_out
3971 : integer, dimension(4) :: nuclides_in, nuclides_out
3972 : logical, parameter :: dbg = .false.
3973 : include 'formats'
3974 : ierr = 0
3975 101144 : num_in = 2
3976 101144 : nuclides_in(1) = i1
3977 101144 : nuclides_in(2) = i2
3978 101144 : num_out = 2
3979 101144 : nuclides_out(1) = o1
3980 101144 : nuclides_out(2) = o2
3981 : call reaclib_rate( &
3982 : str, num_in, nuclides_in, num_out, nuclides_out, tf% T9, &
3983 101144 : fr, rr, ierr)
3984 101144 : end subroutine try1_reaclib_2_2
3985 :
3986 :
3987 10019 : subroutine try1_reaclib_2_3(i1, i2, o1, o2, o3, tf, fr, rr, str, ierr)
3988 : integer, intent(in) :: i1, i2, o1, o2, o3
3989 : type (T_Factors) :: tf
3990 : real(dp), intent(out) :: fr, rr
3991 : character (len=*), intent(in) :: str
3992 : integer, intent(out) :: ierr
3993 : integer :: num_in, num_out
3994 : integer, dimension(4) :: nuclides_in, nuclides_out
3995 : logical, parameter :: dbg = .false.
3996 : include 'formats'
3997 : ierr = 0
3998 10019 : num_in = 2
3999 10019 : nuclides_in(1) = i1
4000 10019 : nuclides_in(2) = i2
4001 10019 : num_out = 3
4002 10019 : nuclides_out(1) = o1
4003 10019 : nuclides_out(2) = o2
4004 10019 : nuclides_out(3) = o3
4005 : call reaclib_rate( &
4006 : str, num_in, nuclides_in, num_out, nuclides_out, tf% T9, &
4007 10019 : fr, rr, ierr)
4008 10019 : end subroutine try1_reaclib_2_3
4009 :
4010 :
4011 10019 : subroutine try1_reaclib_2_4(i1, i2, o1, o2, o3, o4, tf, fr, rr, str, ierr)
4012 : integer, intent(in) :: i1, i2, o1, o2, o3, o4
4013 : type (T_Factors) :: tf
4014 : real(dp), intent(out) :: fr, rr
4015 : character (len=*), intent(in) :: str
4016 : integer, intent(out) :: ierr
4017 : integer :: num_in, num_out
4018 : integer, dimension(4) :: nuclides_in, nuclides_out
4019 : logical, parameter :: dbg = .false.
4020 : include 'formats'
4021 : ierr = 0
4022 10019 : num_in = 2
4023 10019 : nuclides_in(1) = i1
4024 10019 : nuclides_in(2) = i2
4025 10019 : num_out = 4
4026 10019 : nuclides_out(1) = o1
4027 10019 : nuclides_out(2) = o2
4028 10019 : nuclides_out(3) = o3
4029 10019 : nuclides_out(4) = o4
4030 : call reaclib_rate( &
4031 : str, num_in, nuclides_in, num_out, nuclides_out, tf% T9, &
4032 10019 : fr, rr, ierr)
4033 10019 : end subroutine try1_reaclib_2_4
4034 :
4035 :
4036 20038 : subroutine try1_reaclib_3_1(i1, i2, i3, o1, tf, fr, rr, str, ierr)
4037 : integer, intent(in) :: i1, i2, i3, o1
4038 : type (T_Factors) :: tf
4039 : real(dp), intent(out) :: fr, rr
4040 : character (len=*), intent(in) :: str
4041 : integer, intent(out) :: ierr
4042 : integer :: num_in, num_out
4043 : integer, dimension(4) :: nuclides_in, nuclides_out
4044 : logical, parameter :: dbg = .false.
4045 : include 'formats'
4046 : ierr = 0
4047 20038 : num_in = 3
4048 20038 : nuclides_in(1) = i1
4049 20038 : nuclides_in(2) = i2
4050 20038 : nuclides_in(3) = i3
4051 20038 : num_out = 1
4052 20038 : nuclides_out(1) = o1
4053 : call reaclib_rate( &
4054 : str, num_in, nuclides_in, num_out, nuclides_out, tf% T9, &
4055 20038 : fr, rr, ierr)
4056 20038 : end subroutine try1_reaclib_3_1
4057 :
4058 :
4059 0 : subroutine try1_reaclib_3_2(i1, i2, i3, o1, o2, tf, fr, rr, str, ierr)
4060 : integer, intent(in) :: i1, i2, i3, o1, o2
4061 : type (T_Factors) :: tf
4062 : real(dp), intent(out) :: fr, rr
4063 : character (len=*), intent(in) :: str
4064 : integer, intent(out) :: ierr
4065 : integer :: num_in, num_out
4066 : integer, dimension(4) :: nuclides_in, nuclides_out
4067 : logical, parameter :: dbg = .false.
4068 : include 'formats'
4069 : ierr = 0
4070 0 : num_in = 3
4071 0 : nuclides_in(1) = i1
4072 0 : nuclides_in(2) = i2
4073 0 : nuclides_in(3) = i3
4074 0 : num_out = 2
4075 0 : nuclides_out(1) = o1
4076 0 : nuclides_out(2) = o2
4077 : call reaclib_rate( &
4078 : str, num_in, nuclides_in, num_out, nuclides_out, tf% T9, &
4079 0 : fr, rr, ierr)
4080 0 : end subroutine try1_reaclib_3_2
4081 :
4082 :
4083 0 : subroutine try1_reaclib_4_2(i1, i2, i3, i4, o1, o2, tf, fr, rr, str, ierr)
4084 : integer, intent(in) :: i1, i2, i3, i4, o1, o2
4085 : type (T_Factors) :: tf
4086 : real(dp), intent(out) :: fr, rr
4087 : character (len=*), intent(in) :: str
4088 : integer, intent(out) :: ierr
4089 : integer :: num_in, num_out
4090 : integer, dimension(4) :: nuclides_in, nuclides_out
4091 : logical, parameter :: dbg = .false.
4092 : include 'formats'
4093 : ierr = 0
4094 0 : num_in = 4
4095 0 : nuclides_in(1) = i1
4096 0 : nuclides_in(2) = i2
4097 0 : nuclides_in(3) = i3
4098 0 : nuclides_in(4) = i4
4099 0 : num_out = 2
4100 0 : nuclides_out(1) = o1
4101 0 : nuclides_out(2) = o2
4102 : call reaclib_rate( &
4103 : str, num_in, nuclides_in, num_out, nuclides_out, tf% T9, &
4104 0 : fr, rr, ierr)
4105 0 : end subroutine try1_reaclib_4_2
4106 :
4107 :
4108 232219 : subroutine reaclib_rate( &
4109 232219 : str, num_in, nuclides_in, num_out, nuclides_out, T9, &
4110 : lambda, rlambda, ierr)
4111 : use reaclib_support, only: reaction_handle
4112 : character (len=*), intent(in) :: str
4113 : integer, intent(in) :: num_in, nuclides_in(:)
4114 : integer, intent(in) :: num_out, nuclides_out(:)
4115 : real(dp), intent(in) :: T9
4116 : real(dp), intent(out) :: lambda, rlambda
4117 : integer, intent(out) :: ierr
4118 : character (len=max_id_length) :: handle
4119 464438 : integer :: iso_ids(num_in+num_out)
4120 : include 'formats'
4121 232219 : ierr = 0
4122 716695 : iso_ids(1:num_in) = nuclides_in(1:num_in)
4123 615677 : iso_ids(1+num_in:num_out+num_in) = nuclides_out(1:num_out)
4124 232219 : call reaction_handle(num_in, num_out, iso_ids, '-', handle)
4125 232219 : call reaclib_rate_for_handle(handle, T9, lambda, rlambda, ierr)
4126 232219 : end subroutine reaclib_rate
4127 :
4128 :
4129 232219 : subroutine reaclib_rate_for_handle(handle, T9, lambda, rlambda, ierr)
4130 232219 : use reaclib_eval, only: do_reaclib_indices_for_reaction, do_reaclib_reaction_rates
4131 : character (len=*), intent(in) :: handle
4132 : real(dp), intent(in) :: T9
4133 : real(dp), intent(out) :: lambda, rlambda
4134 : integer, intent(out) :: ierr
4135 232219 : real(dp) :: dlambda_dlnT, drlambda_dlnT
4136 : integer :: lo, hi
4137 : logical, parameter :: forward_only = .false.
4138 : include 'formats'
4139 : ierr = 0
4140 232219 : call do_reaclib_indices_for_reaction(handle, reaclib_rates, lo, hi, ierr)
4141 232219 : if (ierr /= 0) then
4142 : return
4143 : write(*,'(a)') &
4144 : 'reaclib_rate_and_dlnT_for_handle: failed in reaclib_indices_for_reaction ' &
4145 : // trim(handle)
4146 : call mesa_error(__FILE__,__LINE__,'reaclib_rate_and_dlnT_for_handle')
4147 : end if
4148 232075 : if (T9 < reaclib_min_T9 .and. reaclib_rates% reaction_flag(lo) /= 'w' .and. &
4149 : reaclib_rates% reaction_flag(lo) /= 'e') then ! w or ec
4150 78610 : lambda = 0; rlambda = 0
4151 78610 : return
4152 : end if
4153 : call do_reaclib_reaction_rates( &
4154 : lo, hi, T9, reaclib_rates, chem_isos, forward_only, &
4155 : lambda, dlambda_dlnT, &
4156 : rlambda, drlambda_dlnT, &
4157 153465 : ierr)
4158 153465 : if (ierr /= 0) then
4159 : write(*,'(a)') &
4160 : 'reaclib_rate_for_handle: failed in reaclib_reaction_rates ' &
4161 0 : // trim(handle)
4162 0 : call mesa_error(__FILE__,__LINE__,'reaclib_rate_for_handle')
4163 0 : return
4164 : end if
4165 232219 : end subroutine reaclib_rate_for_handle
4166 :
4167 :
4168 0 : subroutine reaclib_rate_and_dlnT_for_handle( &
4169 : handle, T9, lambda, dlambda_dlnT, rlambda, drlambda_dlnT, ierr)
4170 232219 : use reaclib_eval, only: do_reaclib_indices_for_reaction, do_reaclib_reaction_rates
4171 : character (len=*), intent(in) :: handle
4172 : real(dp), intent(in) :: T9
4173 : real(dp), intent(out) :: lambda, dlambda_dlnT, rlambda, drlambda_dlnT
4174 : integer, intent(out) :: ierr
4175 : integer :: lo, hi
4176 : include 'formats'
4177 : ierr = 0
4178 0 : call do_reaclib_indices_for_reaction(handle, reaclib_rates, lo, hi, ierr)
4179 0 : if (ierr /= 0) then
4180 : return
4181 : write(*,'(a)') &
4182 : 'reaclib_rate_and_dlnT_for_handle: failed in reaclib_indices_for_reaction ' &
4183 : // trim(handle)
4184 : call mesa_error(__FILE__,__LINE__,'reaclib_rate_and_dlnT_for_handle')
4185 : end if
4186 : call reaclib_rate_and_dlnT( &
4187 0 : lo, hi, handle, T9, lambda, dlambda_dlnT, rlambda, drlambda_dlnT, ierr)
4188 0 : end subroutine reaclib_rate_and_dlnT_for_handle
4189 :
4190 :
4191 1021898 : subroutine reaclib_rate_and_dlnT( &
4192 : lo, hi, handle, T9, lambda, dlambda_dlnT, rlambda, drlambda_dlnT, ierr)
4193 0 : use reaclib_eval, only: do_reaclib_reaction_rates
4194 : integer, intent(in) :: lo,hi
4195 : character (len=*), intent(in) :: handle
4196 : real(dp), intent(in) :: T9
4197 : real(dp), intent(out) :: lambda, dlambda_dlnT, rlambda, drlambda_dlnT
4198 : integer, intent(out) :: ierr
4199 : logical, parameter :: forward_only = .false.
4200 : include 'formats'
4201 1021898 : ierr = 0
4202 1021898 : if (T9 < reaclib_min_T9 .and. reaclib_rates% reaction_flag(lo) /= 'w') then
4203 323206 : lambda = 0; dlambda_dlnT = 0
4204 323206 : rlambda = 0; drlambda_dlnT = 0
4205 323206 : return
4206 : end if
4207 : call do_reaclib_reaction_rates( &
4208 : lo, hi, T9, reaclib_rates, chem_isos, forward_only, &
4209 : lambda, dlambda_dlnT, rlambda, drlambda_dlnT, &
4210 698692 : ierr)
4211 698692 : if (ierr /= 0) then
4212 0 : write(*,*) 'failed in reaclib_reaction_rates ' // trim(handle)
4213 0 : return
4214 : end if
4215 : return
4216 : write(*,3) 'reaclib_rate_and_dlnT_for_handle ' // trim(handle), lo, hi
4217 : write(*,1) 'T9', T9
4218 : write(*,1) 'lambda', lambda
4219 : write(*,1) 'dlambda_dlnT', dlambda_dlnT
4220 : write(*,1) 'rlambda', rlambda
4221 : write(*,1) 'drlambda_dlnT', drlambda_dlnT
4222 : write(*,'(A)')
4223 1021898 : end subroutine reaclib_rate_and_dlnT
4224 :
4225 :
4226 0 : subroutine do_reaclib(tf, a1, a2, a3, a4, a5, a6, a7, term)
4227 : type (T_Factors) :: tf
4228 : real(dp), intent(in) :: a1, a2, a3, a4, a5, a6, a7
4229 : real(dp), intent(out) :: term
4230 : real(dp) :: exponent
4231 : include 'formats'
4232 : !rate = exp(a1 + a2/T9 + a3/T913 + a4*T913 + a5*T9 + a6*T953 + a7*ln(T9))
4233 0 : if (tf% T9 < reaclib_min_T9) then
4234 0 : term = 0; return
4235 : end if
4236 : exponent = a1 + a2*(tf% T9i) + a3*(tf% T9i13) + a4*(tf% T913) + &
4237 0 : a5*(tf% T9) + a6*(tf% T953) + a7*(tf% lnT9)
4238 0 : term = exp(exponent)
4239 0 : if (is_bad(term)) then
4240 0 : write(*,1) 'exponent', exponent
4241 0 : write(*,1) 'term', term
4242 0 : write(*,1) 'a1', a1
4243 0 : write(*,1) 'a2*(tf% T9i)', a2*(tf% T9i)
4244 0 : write(*,1) 'a3*(tf% T9i13)', a3*(tf% T9i13)
4245 0 : write(*,1) 'a4*(tf% T913)', a4*(tf% T913)
4246 0 : write(*,1) 'a5*(tf% T9)', a5*(tf% T9)
4247 0 : write(*,1) 'a6*(tf% T953)', a6*(tf% T953)
4248 0 : write(*,1) 'a7*(tf% lnT9)', a7*(tf% lnT9)
4249 0 : write(*,1) 'tf% lnT9/ln10', tf% lnT9/ln10
4250 0 : call mesa_error(__FILE__,__LINE__,'reaclib')
4251 : end if
4252 1021898 : end subroutine do_reaclib
4253 :
4254 :
4255 : ! returns the indx corresponding to Tpart just less than T9
4256 : ! T9 is the temperature in units of GK
4257 : ! returns a value of 0 or npart if value is less than the minimum or maximum of the partition function
4258 : ! temperature array Tpart
4259 0 : function get_partition_fcn_indx(T9) result(indx)
4260 : real(dp), intent(in) :: T9
4261 : integer :: indx
4262 : integer, parameter :: max_iterations = 8
4263 : integer :: low, high, mid, i
4264 0 : low = 1
4265 0 : high = npart
4266 0 : if (T9 < Tpart(low)) then
4267 0 : indx = low-1
4268 : return
4269 : end if
4270 : if (T9 > Tpart(high)) then
4271 0 : indx = high + 1
4272 : end if
4273 0 : do i = 1, max_iterations
4274 0 : if (high-low <= 1) then
4275 0 : indx = low
4276 0 : return
4277 : end if
4278 0 : mid = (high+low)/2
4279 0 : if (T9 < Tpart(mid)) then
4280 : high = mid
4281 : else
4282 0 : low = mid
4283 : end if
4284 : end do
4285 : ! should never get here
4286 0 : indx = low-1
4287 0 : end function get_partition_fcn_indx
4288 :
4289 :
4290 5 : subroutine mazurek_init(ierr)
4291 : use rates_def, only: tv,rv,rfdm,rfd0,rfd1,rfd2,tfdm,tfd0,tfd1,tfd2
4292 : integer, intent(out) :: ierr
4293 : integer :: k,j
4294 5 : rv(:) = [ 6D0, 7D0, 8D0, 9D0, 10D0, 11D0 ]
4295 5 : tv(:) = [ 2D0, 4D0, 6D0, 8D0, 10D0, 12D0, 14D0 ]
4296 5 : ierr = 0
4297 20 : do k=2,4
4298 15 : rfdm(k)=1.d0/((rv(k-1)-rv(k))*(rv(k-1)-rv(k+1))*(rv(k-1)-rv(k+2)))
4299 15 : rfd0(k)=1.d0/((rv(k)-rv(k-1))*(rv(k)-rv(k+1))*(rv(k)-rv(k+2)))
4300 15 : rfd1(k)=1.d0/((rv(k+1)-rv(k-1))*(rv(k+1)-rv(k))*(rv(k+1)-rv(k+2)))
4301 20 : rfd2(k)=1.d0/((rv(k+2)-rv(k-1))*(rv(k+2)-rv(k))*(rv(k+2)-rv(k+1)))
4302 : end do
4303 25 : do j=2,5
4304 20 : tfdm(j)=1.d0/((tv(j-1)-tv(j))*(tv(j-1)-tv(j+1))*(tv(j-1)-tv(j+2)))
4305 20 : tfd0(j)=1.d0/((tv(j)-tv(j-1))*(tv(j)-tv(j+1))*(tv(j)-tv(j+2)))
4306 20 : tfd1(j)=1.d0/((tv(j+1)-tv(j-1))*(tv(j+1)-tv(j))*(tv(j+1)-tv(j+2)))
4307 25 : tfd2(j)=1.d0/((tv(j+2)-tv(j-1))*(tv(j+2)-tv(j))*(tv(j+2)-tv(j+1)))
4308 : end do
4309 5 : end subroutine mazurek_init
4310 :
4311 0 : subroutine mazurek(btemp,bden,y56,ye,rn56ec,sn56ec)
4312 5 : use rates_def, only: tv,rv,rfdm,rfd0,rfd1,rfd2,tfdm,tfd0,tfd1,tfd2
4313 : real(dp), intent(in) :: btemp,bden,y56,ye
4314 : real(dp), intent(out) :: rn56ec,sn56ec
4315 :
4316 : ! this routine evaluates mazurek's 1973 fits for the ni56 electron
4317 : ! capture rate rn56ec and neutrino loss rate sn56ec
4318 :
4319 : ! input:
4320 : ! y56 = nickel56 molar abundance
4321 : ! ye = electron to baryon number, zbar/abar
4322 :
4323 : ! output:
4324 : ! rn56ec = ni56 electron capture rate
4325 : ! sn56ec = ni56 neutrino loss rate
4326 :
4327 : ! declare
4328 : integer :: jp,kp,jr,jd,ii,ik,ij
4329 0 : real(dp) :: rnt(2),rne(2,7),datn(2,6,7), &
4330 0 : t9,r,rfm,rf0,rf1,rf2,dfacm,dfac0,dfac1,dfac2, &
4331 0 : tfm,tf0,tf1,tf2,tfacm,tfac0,tfac1,tfac2
4332 :
4333 : ! initialize
4334 : data (((datn(ii,ik,ij),ik=1,6),ij=1,7),ii=1,1) / &
4335 : -3.98d0, -2.84d0, -1.41d0, 0.20d0, 1.89d0, 3.63d0, &
4336 : -3.45d0, -2.62d0, -1.32d0, 0.22d0, 1.89d0, 3.63d0, &
4337 : -2.68d0, -2.30d0, -1.19d0, 0.27d0, 1.91d0, 3.62d0, &
4338 : -2.04d0, -1.87d0, -1.01d0, 0.34d0, 1.94d0, 3.62d0, &
4339 : -1.50d0, -1.41d0, -0.80d0, 0.45d0, 1.99d0, 3.60d0, &
4340 : -1.00d0, -0.95d0, -0.54d0, 0.60d0, 2.06d0, 3.58d0, &
4341 : -0.52d0, -0.49d0, -0.21d0, 0.79d0, 2.15d0, 3.55d0 /
4342 : data (((datn(ii,ik,ij),ik=1,6),ij=1,7),ii=2,2) / &
4343 : -3.68d0, -2.45d0, -0.80d0, 1.12d0, 3.13d0, 5.19d0, &
4344 : -2.91d0, -2.05d0, -0.64d0, 1.16d0, 3.14d0, 5.18d0, &
4345 : -1.95d0, -1.57d0, -0.40d0, 1.24d0, 3.16d0, 5.18d0, &
4346 : -1.16d0, -0.99d0, -0.11d0, 1.37d0, 3.20d0, 5.18d0, &
4347 : -0.48d0, -0.40d0, 0.22d0, 1.54d0, 3.28d0, 5.16d0, &
4348 : 0.14d0, 0.19d0, 0.61d0, 1.78d0, 3.38d0, 5.14d0, &
4349 : 0.75d0, 0.78d0, 1.06d0, 2.07d0, 3.51d0, 5.11d0 /
4350 :
4351 : ! calculate ni56 electron capture and neutrino loss rates
4352 0 : rn56ec = 0.0d0
4353 0 : sn56ec = 0.0d0
4354 :
4355 0 : if (btemp*1d-9 < lowT9_cutoff) return
4356 :
4357 0 : if ( (btemp < 2.0d9) .or. (bden*ye < 1.0d6)) return
4358 0 : t9 = min(btemp, 1.4d10) * 1.0d-9
4359 0 : r = max(6.0d0,min(11.0d0,log10(bden*ye)))
4360 0 : jp = min(max(2,int(0.5d0*t9)), 5)
4361 0 : kp = min(max(2,int(r)-5), 4)
4362 0 : rfm = r - rv(kp-1)
4363 0 : rf0 = r - rv(kp)
4364 0 : rf1 = r - rv(kp+1)
4365 0 : rf2 = r - rv(kp+2)
4366 0 : dfacm = rf0*rf1*rf2*rfdm(kp)
4367 0 : dfac0 = rfm*rf1*rf2*rfd0(kp)
4368 0 : dfac1 = rfm*rf0*rf2*rfd1(kp)
4369 0 : dfac2 = rfm*rf0*rf1*rfd2(kp)
4370 0 : tfm = t9 - tv(jp-1)
4371 0 : tf0 = t9 - tv(jp)
4372 0 : tf1 = t9 - tv(jp+1)
4373 0 : tf2 = t9 - tv(jp+2)
4374 0 : tfacm = tf0*tf1*tf2*tfdm(jp)
4375 0 : tfac0 = tfm*tf1*tf2*tfd0(jp)
4376 0 : tfac1 = tfm*tf0*tf2*tfd1(jp)
4377 0 : tfac2 = tfm*tf0*tf1*tfd2(jp)
4378 :
4379 : ! evaluate the spline fits
4380 0 : do jr = 1,2
4381 0 : do jd = jp-1,jp+2
4382 : rne(jr,jd) = dfacm*datn(jr,kp-1,jd) + dfac0*datn(jr,kp,jd) &
4383 0 : + dfac1*datn(jr,kp+1,jd) + dfac2*datn(jr,kp+2,jd)
4384 : end do
4385 : rnt(jr) = tfacm*rne(jr,jp-1) + tfac0*rne(jr,jp) &
4386 0 : + tfac1*rne(jr,jp+1) + tfac2*rne(jr,jp+2)
4387 : end do
4388 :
4389 : ! set the output
4390 0 : rn56ec = exp10(rnt(1))
4391 0 : sn56ec = 6.022548d+23 * 8.18683d-7 * y56 * exp10(rnt(2))
4392 0 : return
4393 0 : end subroutine mazurek
4394 :
4395 :
4396 0 : subroutine n14_electron_capture_rate(T,Rho,UE,rate)
4397 : real(dp), intent(in) :: T ! temperature
4398 : real(dp), intent(in) :: Rho ! density
4399 : real(dp), intent(in) :: UE ! electron molecular weight
4400 : real(dp), intent(out) :: rate ! (s^-1)
4401 :
4402 0 : real(dp) :: Q, AMC2, AMULTIP, AL92, T8, X, XFER, EF, Y, AA, GUESS, ELCAP
4403 :
4404 : ! from Lars
4405 :
4406 :
4407 : ! Inputs are T in K, rho in gr/cm^3, and UE=electron mean mol. weight
4408 : !
4409 : ! Gives a reasonable estimate (i.e. within factor of 50% or so) of the
4410 : ! electron capture rate for electrons on 14N in a plasma assumed to be quite
4411 : ! degenerate.
4412 : !
4413 : ! x=KT/Q, y=E_FERMI/Q
4414 : !
4415 : ! ELCAP is the rate in 1/seconds
4416 : !
4417 : !
4418 : ! Let's start by putting in the Q value, electron rest mass and
4419 : ! temperature in units of keV.
4420 : !
4421 : !
4422 0 : Q = 667.479d0
4423 0 : AMC2 = 510.999d0
4424 0 : AMULTIP = 0.693d0/1.104d9
4425 0 : AL92 = LOG(9d0/2d0)
4426 0 : T8 = T/1d8
4427 0 : X = 8.617d0*T8/Q
4428 : !
4429 : ! For this value of the density, find the electron fermi momentum
4430 : ! assuming that the KT corrections to the electron EOS are not
4431 : ! important.
4432 : !
4433 0 : XFER = pow(RHO/(0.9739D6*UE),1d0/3d0)
4434 : !
4435 : ! The parameter we need that is used in the fitting formula is
4436 : ! the electron Fermi energy
4437 : !
4438 0 : EF = AMC2*SQRT(1.0D0 + XFER*XFER)
4439 0 : Y = EF/Q
4440 0 : IF(Y < (1.0D0 + AL92*X)) THEN
4441 0 : AA = (Y-1.0D0)/X
4442 0 : GUESS = 2.0D0*X*X*X*exp(AA)
4443 : ELSE
4444 0 : GUESS = pow3(Y-1.0D0+(3.0D0-AL92)*X)/3.0D0
4445 : end if
4446 : !
4447 : ! Now multiply by the prefactors .. .
4448 : !
4449 0 : ELCAP = GUESS*AMULTIP
4450 :
4451 :
4452 0 : rate = ELCAP
4453 :
4454 0 : end subroutine n14_electron_capture_rate
4455 :
4456 :
4457 9106 : subroutine ecapnuc(etakep,temp,rho,rpen,rnep,spen,snep)
4458 : real(dp), intent(in) :: etakep,temp,rho
4459 : real(dp), intent(out) :: rpen,rnep,spen,snep
4460 :
4461 : ! given the electron degeneracy parameter etakep (chemical potential
4462 : ! without the electron's rest mass divided by kt) and the temperature temp,
4463 : ! this routine calculates rates for
4464 : ! electron capture on protons rpen (captures/sec/proton),
4465 : ! positron capture on neutrons rnep (captures/sec/neutron),
4466 : ! and their associated neutrino energy loss rates
4467 : ! spen (ergs/sec/proton) and snep (ergs/sec/neutron)
4468 :
4469 : ! declare
4470 : integer :: iflag
4471 : real(dp), parameter :: c2me = me * clight * clight
4472 : real(dp), parameter :: cmk5 = (kerg / c2me) * (kerg / c2me) * (kerg / c2me) * (kerg / c2me) * (kerg / c2me)
4473 : real(dp), parameter :: cmk6 = cmk5 * (kerg / c2me)
4474 9106 : real(dp) :: t9,t5,qn,etaef,etael,zetan,eta,etael2, &
4475 9106 : etael3,etael4,f1l,f2l,f3l,f4l,f5l,f1g, &
4476 9106 : f2g,f3g,f4g,f5g,exmeta,eta2,eta3,eta4, &
4477 9106 : fac0,fac1,fac2,fac3,rie1,rie2,facv0,facv1, &
4478 9106 : facv2,facv3,facv4,rjv1,rjv2,spenc,snepc, &
4479 9106 : exeta,zetan2,f0,etael5, &
4480 : qn1,ft,qn2, &
4481 : qndeca,tmean, &
4482 : rho_low_cutoff, eta_low_cutoff
4483 : parameter (qn1 = -2.0716446d-06, &
4484 : ft = 1083.9269d0, &
4485 : qn2 = 2.0716446d-06, &
4486 : qndeca = 1.2533036d-06, &
4487 : tmean = 886.7d0, &
4488 : rho_low_cutoff = 1d-9, eta_low_cutoff = -50d0)
4489 :
4490 :
4491 : ! tmean and qndeca are the mean lifetime and decay energy of the neutron
4492 : ! c2me is the constant used to convert the neutrino energy
4493 : ! loss rate from mec2/s (as in the paper) to ergs/particle/sec.
4494 :
4495 : ! initialize
4496 9106 : rpen = 0.0d0
4497 9106 : rnep = 0.0d0
4498 9106 : spen = 0.0d0
4499 9106 : snep = 0.0d0
4500 9106 : t9 = temp * 1.0d-9
4501 :
4502 9106 : if (t9 < lowT9_cutoff .or. rho < rho_low_cutoff .or. etakep < eta_low_cutoff) return
4503 :
4504 9104 : iflag = 0
4505 9104 : qn = qn1
4506 :
4507 :
4508 : ! chemical potential including the electron rest mass
4509 9104 : etaef = etakep + c2me/kerg/temp
4510 :
4511 :
4512 : ! iflag=1 is for electrons, iflag=2 is for positrons
4513 18208 : 502 iflag = iflag + 1
4514 18208 : if (iflag==1) etael = qn2/kerg/temp
4515 18208 : if (iflag==2) etael = c2me/kerg/temp
4516 9104 : if (iflag==2) etaef = -etaef
4517 :
4518 18208 : t5 = temp*temp*temp*temp*temp
4519 18208 : zetan = qn/kerg/temp
4520 18208 : eta = etaef - etael
4521 :
4522 : ! protect from overflowing with large eta values
4523 18208 : if (eta <= 6.8d+02) then
4524 18208 : exeta = exp(eta)
4525 : else
4526 0 : exeta = 0.0d0
4527 : end if
4528 18208 : etael2 = etael*etael
4529 18208 : etael3 = etael2*etael
4530 18208 : etael4 = etael3*etael
4531 18208 : etael5 = etael4*etael
4532 18208 : zetan2 = zetan*zetan
4533 18208 : if (eta <= 6.8d+02) then
4534 18208 : f0 = log1p(exeta)
4535 : else
4536 : f0 = eta
4537 : end if
4538 :
4539 : ! if eta le. 0., the following fermi integrals apply
4540 18208 : f1l = exeta
4541 18208 : f2l = 2.0d0 * f1l
4542 18208 : f3l = 6.0d0 * f1l
4543 18208 : f4l = 24.0d0 * f1l
4544 18208 : f5l = 120.0d0 * f1l
4545 :
4546 : ! if eta gt. 0., the following fermi integrals apply:
4547 18208 : f1g = 0.0d0
4548 18208 : f2g = 0.0d0
4549 18208 : f3g = 0.0d0
4550 18208 : f4g = 0.0d0
4551 18208 : f5g = 0.0d0
4552 18208 : if (eta > 0.0d0) then
4553 1938 : exmeta = exp(-eta)
4554 1938 : eta2 = eta*eta
4555 1938 : eta3 = eta2*eta
4556 1938 : eta4 = eta3*eta
4557 1938 : f1g = 0.5d0*eta2 + 2.0d0 - exmeta
4558 1938 : f2g = eta3/3.0d0 + 4.0d0*eta + 2.0d0*exmeta
4559 1938 : f3g = 0.25d0*eta4 + 0.5d0*pi2*eta2 + 12.0d0 - 6.0d0*exmeta
4560 : f4g = 0.2d0*eta4*eta + 2.0d0*pi2/3.0d0*eta3 + 48.0d0*eta &
4561 1938 : + 24.0d0*exmeta
4562 : f5g = eta4*eta2/6.0d0 + 5.0d0/6.0d0*pi2*eta4 &
4563 1938 : + 7.0d0/6.0d0*pi2*eta2 + 240.0d0 -120.d0*exmeta
4564 : end if
4565 :
4566 : ! factors which are multiplied by the fermi integrals
4567 18208 : fac3 = 2.0d0*zetan + 4.0d0*etael
4568 18208 : fac2 = 6.0d0*etael2 + 6.0d0*etael*zetan + zetan2
4569 18208 : fac1 = 4.0d0*etael3 + 6.0d0*etael2*zetan + 2.0d0*etael*zetan2
4570 18208 : fac0 = etael4 + 2.0d0*zetan*etael3 + etael2*zetan2
4571 :
4572 : ! electron capture rates onto protons with no blocking
4573 18208 : rie1 = f4l + fac3*f3l + fac2*f2l + fac1*f1l + fac0*f0
4574 18208 : rie2 = f4g + fac3*f3g + fac2*f2g + fac1*f1g + fac0*f0
4575 :
4576 : ! neutrino emission rate for electron capture:
4577 18208 : facv4 = 5.0d0*etael + 3.0d0*zetan
4578 18208 : facv3 = 10.0d0*etael2 + 12.0d0*etael*zetan + 3.0d0*zetan2
4579 : facv2 = 10.0d0*etael3 + 18.0d0*etael2*zetan &
4580 18208 : + 9.0d0*etael*zetan2 + zetan2*zetan
4581 : facv1 = 5.0d0*etael4 + 12.0d0*etael3*zetan &
4582 18208 : + 9.0d0*etael2*zetan2 + 2.0d0*etael*zetan2*zetan
4583 : facv0 = etael5 + 3.0d0*etael4*zetan &
4584 18208 : + 3.0d0*etael3*zetan2 + etael2*zetan2*zetan
4585 : rjv1 = f5l + facv4*f4l + facv3*f3l &
4586 18208 : + facv2*f2l + facv1*f1l + facv0*f0
4587 : rjv2 = f5g + facv4*f4g + facv3*f3g &
4588 18208 : + facv2*f2g + facv1*f1g + facv0*f0
4589 :
4590 : ! for electrons capture onto protons
4591 18208 : if (iflag == 2) GOTO 503
4592 9104 : if (eta > 0d0) GOTO 505
4593 7166 : rpen = ln2*cmk5*t5*rie1/ft
4594 7166 : spen = ln2*cmk6*t5*temp*rjv1/ft
4595 7166 : spenc = ln2*cmk6*t5*temp*rjv1/ft*c2me
4596 9104 : GOTO 504
4597 1938 : 505 rpen = ln2*cmk5*t5*rie2/ft
4598 1938 : spen = ln2*cmk6*t5*temp*rjv2/ft
4599 1938 : spenc = ln2*cmk6*t5*temp*rjv2/ft*c2me
4600 : 504 continue
4601 : qn = qn2
4602 9104 : GOTO 502
4603 :
4604 : ! for positrons capture onto neutrons
4605 9104 : 503 if (eta>0d0) GOTO 507
4606 9104 : rnep = ln2*cmk5*t5*rie1/ft
4607 9104 : snep = ln2*cmk6*t5*temp*rjv1/ft
4608 9104 : snepc = ln2*cmk6*t5*temp*rjv1/ft*c2me
4609 : ! if (rho.lt.1.0d+06) snep=snep+qndeca*xn(9)/mn/tmean
4610 9104 : GOTO 506
4611 0 : 507 rnep = ln2*cmk5*t5*rie2/ft
4612 0 : snep = ln2*cmk6*t5*temp*rjv2/ft
4613 0 : snepc = ln2*cmk6*t5*temp*rjv2/ft*c2me
4614 : ! if (rho.lt.1.0d+06) snep=snep+qndeca*xn(9)/mn/tmean
4615 : 506 continue
4616 : return
4617 : end subroutine ecapnuc
4618 :
4619 : end module ratelib
|
