Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module reaclib_support
21 : use rates_def
22 : use math_lib
23 : use chem_lib
24 :
25 : implicit none
26 :
27 : contains
28 :
29 :
30 3 : subroutine set_up_network_information(rates)
31 : type(reaction_data), intent(inout) :: rates
32 : integer :: current_chapter, i
33 :
34 : ! set up bookmarks
35 3 : current_chapter = 0
36 3 : rates% nchapters_included = 0
37 36 : rates% chapters_present = 0
38 102 : rates% bookmarks = 0
39 145557 : do i = 1, rates% nreactions
40 145557 : new_chapter : if (rates% chapter(i) /= current_chapter) then
41 : ! close out the chapter we just left.
42 33 : if (current_chapter /= 0) rates% bookmarks(2,current_chapter) = i-1
43 : ! set up information on the new chapter
44 33 : current_chapter = rates% chapter(i)
45 33 : rates% nchapters_included = rates% nchapters_included + 1
46 33 : rates% chapters_present(rates% nchapters_included) = current_chapter
47 33 : rates% bookmarks(1,current_chapter) = i
48 : end if new_chapter
49 : end do
50 : ! mark the end of the last chapter
51 3 : rates% bookmarks(2,current_chapter) = rates% nreactions
52 3 : end subroutine set_up_network_information
53 :
54 :
55 3 : subroutine assign_weights(rates)
56 : type(reaction_data), intent(inout) :: rates
57 : integer :: i, i1, i2, i3, i4
58 :
59 : include 'formats'
60 :
61 : ! check for allocation
62 3 : if (.not.associated(rates% weight)) then
63 : return
64 : end if
65 :
66 145557 : do i = 1, rates% nreactions
67 145554 : i1 = -1; i2 = -2; i3 = -3; i4 = -4
68 193119 : select case (rates% chapter(i))
69 : case (r_one_one)
70 : case (r_one_two)
71 : case (r_one_three)
72 : case (r_two_one)
73 47565 : i1 = rates% pspecies(1,i)
74 47565 : i2 = rates% pspecies(2,i)
75 : case (r_two_two)
76 45267 : i1 = rates% pspecies(1,i)
77 45267 : i2 = rates% pspecies(2,i)
78 : case (r_two_three)
79 69 : i1 = rates% pspecies(1,i)
80 69 : i2 = rates% pspecies(2,i)
81 : case (r_two_four)
82 12 : i1 = rates% pspecies(1,i)
83 12 : i2 = rates% pspecies(2,i)
84 : case (r_three_one)
85 21 : i1 = rates% pspecies(1,i)
86 21 : i2 = rates% pspecies(2,i)
87 21 : i3 = rates% pspecies(3,i)
88 : case (r_three_two)
89 3 : i1 = rates% pspecies(1,i)
90 3 : i2 = rates% pspecies(2,i)
91 3 : i3 = rates% pspecies(3,i)
92 : case (r_four_two)
93 6 : i1 = rates% pspecies(1,i)
94 6 : i2 = rates% pspecies(2,i)
95 6 : i3 = rates% pspecies(3,i)
96 145554 : i4 = rates% pspecies(4,i)
97 : case (r_one_four)
98 : end select
99 145557 : call set_weight(rates% weight(i))
100 : end do
101 :
102 145557 : do i = 1, rates% nreactions
103 145554 : i1 = -1; i2 = -2; i3 = -3; i4 = -4
104 161877 : select case (rates% chapter(i))
105 : case (r_one_one)
106 : case (r_one_two)
107 16323 : i1 = rates% pspecies(2,i)
108 16323 : i2 = rates% pspecies(3,i)
109 : case (r_one_three)
110 7740 : i1 = rates% pspecies(2,i)
111 7740 : i2 = rates% pspecies(3,i)
112 7740 : i3 = rates% pspecies(4,i)
113 : case (r_two_one)
114 : case (r_two_two)
115 45267 : i1 = rates% pspecies(3,i)
116 45267 : i2 = rates% pspecies(4,i)
117 : case (r_two_three)
118 69 : i1 = rates% pspecies(3,i)
119 69 : i2 = rates% pspecies(4,i)
120 69 : i3 = rates% pspecies(5,i)
121 : case (r_two_four)
122 12 : i1 = rates% pspecies(3,i)
123 12 : i2 = rates% pspecies(4,i)
124 12 : i3 = rates% pspecies(5,i)
125 12 : i4 = rates% pspecies(6,i)
126 : case (r_three_one)
127 : case (r_three_two)
128 3 : i1 = rates% pspecies(4,i)
129 3 : i2 = rates% pspecies(5,i)
130 : case (r_four_two)
131 6 : i1 = rates% pspecies(5,i)
132 6 : i2 = rates% pspecies(6,i)
133 : case (r_one_four)
134 6375 : i1 = rates% pspecies(2,i)
135 6375 : i2 = rates% pspecies(3,i)
136 6375 : i3 = rates% pspecies(4,i)
137 145554 : i4 = rates% pspecies(5,i)
138 : end select
139 :
140 145557 : call set_weight(rates% weight_reverse(i))
141 :
142 : end do
143 :
144 :
145 : contains
146 :
147 :
148 291108 : subroutine set_weight(w)
149 : ! nuclei are sorted, so if identical, then are adjacent in list
150 : real(dp), intent(out) :: w
151 291108 : if (i1 == i2 .and. i2 == i3 .and. i3 == i4) then
152 0 : w = 1d0/24d0
153 291108 : else if (i2 == i3 .and. (i1 == i2 .or. i3 == i4)) then
154 6399 : w = 1d0/6d0
155 284709 : else if (i1 == i2) then
156 7830 : if (i3 == i4) then
157 6 : w = 1d0/4d0
158 : else
159 7824 : w = 1d0/2d0
160 : end if
161 276879 : else if (i2 == i3 .or. i3 == i4) then
162 54 : w = 1d0/2d0
163 : else
164 276825 : w = 1d0
165 : end if
166 291108 : end subroutine set_weight
167 :
168 :
169 : end subroutine assign_weights
170 :
171 :
172 3 : subroutine compute_rev_ratio(rates,winvn)
173 : use const_def, only : pi, kB=>boltzm, NA=>avo, hbar, &
174 : c=>clight, conv=>mev_to_ergs
175 : type(reaction_data), intent(inout) :: rates
176 : type(nuclide_data), intent(in) :: winvn
177 : real(dp) :: mp, mn
178 : real(dp), dimension(max_species_per_reaction) :: g ! statistical weights of nuclides
179 : real(dp), dimension(max_species_per_reaction) :: mass ! mass no's of nuclides
180 : integer, dimension(max_species_per_reaction) :: ps
181 : integer :: Ni,No,Nt,i
182 : real(dp) :: fac, massfac, sum1, sum2, tmp
183 :
184 :
185 : include 'formats'
186 :
187 : ! Get these consistently from the isotopes.data file
188 3 : mp=winvn%W(chem_get_iso_id('prot'))
189 3 : mn=winvn%W(chem_get_iso_id('neut'))
190 :
191 3 : fac = pow(1d9*kB/(2d0*pi*hbar*hbar*NA),1.5d0)/NA
192 3 : massfac = conv*NA/(c*c)
193 :
194 145557 : rates% weak_mask = 1d0
195 145557 : loop_over_rates: do i = 1,rates% nreactions
196 : ! weak rates don't have inverses, so set to innocuous values and mask out
197 145554 : if (rates% reaction_flag(i) == 'w' .or. rates% reaction_flag(i) == 'e') then
198 45264 : rates% weak_mask(i) = 0d0
199 135792 : rates% inverse_coefficients(:,i) = [-huge(1d0), 0d0]
200 45264 : rates% inverse_exp(i) = 0d0
201 1131600 : rates% inverse_part(:,i) = 1d0
202 : cycle loop_over_rates
203 : end if
204 100290 : Ni = Nin(rates% chapter(i))
205 100290 : No = Nout(rates% chapter(i))
206 100290 : Nt = Ni+No
207 446646 : ps(1:Nt) = rates% pspecies(1:Nt,i)
208 446646 : g(1:Nt) = 2d0*winvn% spin(ps(1:Nt)) + 1d0
209 : ! mass(1:Nt) = winvn% Z(ps(1:Nt))*mp + winvn% N(ps(1:Nt))*mn - winvn% binding_energy(ps(1:Nt))*massfac
210 446646 : mass(1:Nt) = winvn% W(ps(1:Nt))
211 :
212 : ! log(prefactor of reverse_ratio)
213 446646 : tmp = product(mass(1:Ni))/product(mass(Ni+1:Nt))
214 : ! The reactant/product mass ratio always enters with a fixed 3/2 power.
215 : ! This follows standard detailed-balance derivations such as
216 : ! Fowler et al. 1967 and Smith Clark et al. 2023 Appendix B.
217 446646 : rates% inverse_coefficients(1,i) = pow(tmp,1.5d0)*(product(g(1:Ni))/product(g(Ni+1:Nt)))
218 :
219 : ! -Q/(kB*10**9)
220 293550 : sum1 = sum(winvn% binding_energy(ps(1:Ni)))
221 253386 : sum2 = sum(winvn% binding_energy(ps(Ni+1:Nt)))
222 100290 : rates% inverse_coefficients(2,i) = (sum1-sum2)*conv/kB/1d9
223 :
224 : ! Reichert et al. 2023 Eq. 21 is the source for the fac^(Ni-No) term,
225 : ! but the fixed 3/2 mass exponent above follows Fowler et al. 1967 and
226 : ! Smith Clark et al. 2023 Appendix B. The current WinNet implementation
227 : ! also uses the fixed 3/2 mass exponent.
228 : ! delta_reactants/delta_products is handled in net.
229 : ! fac == (mu * kb * T / (2 * pi * hbar^2))^3/2 term without a T^3/2.
230 : ! The n-dependent phase-space factor enters as fac^(Ni-No), where
231 : ! n = Ni - No. The mass ratio above remains a separate fixed 3/2 term.
232 :
233 : ! The T makes its way back into our expression inside
234 : ! the subroutine compute_some_inverse_lambdas, in reaclib_eval.f90.
235 : ! It appears in log form as 1.5d0*rates% inverse_exp(i)*lnT9, where,
236 : ! rates% inverse_exp(i) = Ni-No, so,
237 : ! 1.5d0*rates% inverse_exp(i)*lnT9 == ln(T^(3n/2)), where n = Ni-No.
238 100290 : rates% inverse_exp(i) = Ni - No ! We use this term in a log() expression in reaclib_eval
239 :
240 : ! Ni-No should be 0 for non-photo-disintegration reverse rates
241 : ! and >=1 for photos's in the reverse channel
242 100290 : if (Ni-No /= 0) then
243 : ! whether endothermic or exothermic,
244 : ! Ni-No handles the sign of Q from rates% inverse_coefficients(2,i)
245 55023 : rates% inverse_coefficients(1,i) = rates% inverse_coefficients(1,i) * pow(fac, Ni-No)
246 : ! negative values of Q denote endothermic photodisintegrations
247 : ! We multiply by fac^|Ni-No| as we want to compute
248 : ! rate_photo/rate_forward
249 :
250 : ! positive values of Q denote exothermic photodisintegrations
251 : ! We DIVIDE by fac^|Ni-No| as we want to compute
252 : ! rate_reverse/rate_photo
253 : else ! Ni - No = 0
254 45267 : rates% inverse_exp(i) = 0 ! ensure this is 0.
255 : rates% inverse_coefficients(1,i) = rates% inverse_coefficients(1,i) ! no point calling pow(fac,0)
256 : end if
257 100290 : rates% inverse_coefficients(1,i) = log(rates% inverse_coefficients(1,i))
258 :
259 :
260 : rates% inverse_part(:,i) = product(winvn% pfcn(1:npart,ps(1:Ni)),dim=2)/ &
261 10819797 : & product(winvn% pfcn(1:npart,ps(Ni+1:Nt)),dim=2)
262 : end do loop_over_rates
263 3 : end subroutine compute_rev_ratio
264 :
265 0 : subroutine do_parse_reaction_handle(handle, num_in, num_out, iso_ids, op, ierr)
266 : use chem_def
267 : use chem_lib
268 : character (len=*), intent(in) :: handle
269 : integer, intent(out) :: num_in, num_out
270 : integer, intent(out) :: iso_ids(:) ! holds chem_ids for input and output species
271 : character (len=*), intent(out) :: op ! e.g., 'pg', 'wk', 'to', or ...
272 : integer, intent(out) :: ierr
273 :
274 : integer :: len, i, j, cnt, cid, extra_in, extra_out
275 : logical :: doing_inputs
276 :
277 0 : num_in = 0; num_out = 0; op = ''
278 0 : ierr = -1
279 0 : len = len_trim(handle)
280 0 : if (handle(1:2) /= 'r_') return
281 0 : i = 3
282 0 : cnt = 0
283 0 : doing_inputs = .true.
284 0 : do while (i <= len)
285 0 : call nxt ! set j to last char of token
286 0 : cid = chem_get_iso_id(handle(i:j))
287 0 : if (cid == nuclide_not_found) then
288 0 : if (doing_inputs) then
289 0 : op = handle(i:j)
290 0 : extra_in = -1
291 0 : extra_out = -1
292 0 : if (j == i+1) then ! check 2 character ops
293 0 : select case (op(1:1))
294 : case ('p')
295 0 : extra_in = ih1
296 : case ('a')
297 0 : extra_in = ihe4
298 : case ('n')
299 0 : extra_in = ineut
300 : case ('g')
301 0 : extra_in = 0
302 : case default
303 : end select
304 0 : if (extra_in >= 0) then
305 0 : if (extra_in > 0) then
306 0 : cnt = cnt+1
307 0 : if (cnt /= 2) then
308 : !write(*,*) 'failed to parse ' // &
309 : ! trim(handle) // ' -- problem with ' // handle(i:j)
310 : return
311 : end if
312 0 : if (chem_isos% Z(iso_ids(1)) >= chem_isos% Z(extra_in)) then
313 0 : iso_ids(2) = iso_ids(1)
314 0 : iso_ids(1) = extra_in
315 : else
316 0 : iso_ids(2) = extra_in
317 : end if
318 : end if
319 0 : select case (op(2:2))
320 : case ('p')
321 0 : extra_out = ih1
322 : case ('a')
323 0 : extra_out = ihe4
324 : case ('n')
325 0 : extra_out = ineut
326 : case ('g')
327 : extra_out = 0
328 : case default
329 : !write(*,*) 'failed to parse ' // &
330 : ! trim(handle) // ' -- problem with ' // handle(i:j)
331 0 : return
332 : end select
333 : end if
334 : end if
335 0 : num_in = cnt
336 0 : doing_inputs = .false.
337 0 : if (extra_out > 0) then
338 0 : cnt = cnt+1
339 0 : iso_ids(cnt) = extra_out
340 : end if
341 : else
342 : !write(*,*) 'failed to parse ' // &
343 : ! trim(handle) // ' -- problem with ' // handle(i:j)
344 : return
345 : end if
346 : else
347 0 : cnt = cnt+1
348 0 : iso_ids(cnt) = cid
349 : end if
350 0 : i = j+2
351 : end do
352 0 : num_out = cnt - num_in
353 0 : ierr = 0
354 :
355 : contains
356 :
357 0 : subroutine nxt
358 0 : j = i
359 : do
360 0 : if (j >= len) return
361 0 : j = j+1
362 0 : if (handle(j:j) == '_') then
363 0 : j = j-1; return
364 : end if
365 : end do
366 : end subroutine nxt
367 :
368 :
369 : end subroutine do_parse_reaction_handle
370 :
371 231121 : subroutine reaction_handle(num_in, num_out, iso_ids, reaction_flag, handle)
372 : use chem_def, only: chem_isos
373 : integer, intent(in) :: num_in, num_out
374 : integer, intent(in) :: iso_ids(:)
375 : character (len=*), intent(in) :: reaction_flag
376 : character (len=*), intent(out) :: handle
377 : logical, parameter :: reverse = .false.
378 231121 : call get1_reaction_handle(num_in, num_out, iso_ids, chem_isos, reverse, reaction_flag, handle)
379 231121 : end subroutine reaction_handle
380 :
381 0 : subroutine reverse_reaction_handle(num_in, num_out, iso_ids, handle)
382 : use chem_def, only: chem_isos
383 : integer, intent(in) :: num_in, num_out
384 : integer, intent(in) :: iso_ids(:)
385 : character (len=*), intent(out) :: handle
386 : logical, parameter :: reverse = .true.
387 : character (len=1), parameter :: reaction_flag = '-'
388 0 : call get1_reaction_handle(num_in, num_out, iso_ids, chem_isos, reverse, reaction_flag, handle)
389 0 : end subroutine reverse_reaction_handle
390 :
391 145554 : subroutine get_reaction_handle(num_in, num_out, pspecies, nuclides, reaction_flag, handle)
392 : integer, intent(in) :: num_in, num_out
393 : integer, intent(in) :: pspecies(:)
394 : type(nuclide_data), intent(in) :: nuclides
395 : character (len=*), intent(in) :: reaction_flag
396 : character (len=*), intent(out) :: handle
397 : logical, parameter :: reverse = .false.
398 145554 : call get1_reaction_handle(num_in, num_out, pspecies, nuclides, reverse, reaction_flag, handle)
399 145554 : end subroutine get_reaction_handle
400 :
401 100290 : subroutine get_reverse_reaction_handle(num_in, num_out, pspecies, nuclides, handle)
402 : integer, intent(in) :: num_in, num_out
403 : integer, intent(in) :: pspecies(:)
404 : type(nuclide_data), intent(in) :: nuclides
405 : character (len=*), intent(out) :: handle
406 : logical, parameter :: reverse = .true.
407 : character (len=1), parameter :: reaction_flag = '-'
408 100290 : call get1_reaction_handle(num_in, num_out, pspecies, nuclides, reverse, reaction_flag, handle)
409 100290 : end subroutine get_reverse_reaction_handle
410 :
411 476965 : subroutine get1_reaction_handle( &
412 476965 : num_in, num_out, pspecies_in, nuclides, reverse, reaction_flag, handle)
413 : use chem_def, only: ih1, ih2, ih3, ihe3, ihe4, ibe7, ili7, chem_isos
414 : integer, intent(in) :: num_in, num_out
415 : integer, intent(in) :: pspecies_in(:)
416 : type(nuclide_data), intent(in) :: nuclides
417 : logical, intent(in) :: reverse
418 : character (len=*), intent(in) :: reaction_flag
419 : character (len=*), intent(out) :: handle
420 :
421 476965 : integer :: in1, in2, out1, out2, num, pspecies(num_in + num_out)
422 : logical :: do_long_form, ec_flag, wk_flag
423 :
424 : include 'formats'
425 :
426 476965 : num = num_in + num_out
427 2167796 : pspecies(1:num) = pspecies_in(1:num)
428 476965 : call sort(num_in, pspecies(1:num_in))
429 476965 : call sort(num_out, pspecies(num_in+1:num))
430 476965 : ec_flag = (reaction_flag == 'e')
431 476965 : wk_flag = (reaction_flag == 'w')
432 :
433 476965 : if (ec_flag) then ! special cases
434 7 : if (reverse) then
435 0 : handle = ''
436 0 : return
437 : end if
438 7 : if (num_in == 2 .and. num_out == 1) then
439 : if (nuclides% chem_id(pspecies(1)) == ih1 .and. &
440 3 : nuclides% chem_id(pspecies(2)) == ih1 .and. &
441 : nuclides% chem_id(pspecies(3)) == ih2) then
442 3 : handle = 'r_h1_h1_ec_h2'
443 3 : return
444 : end if
445 4 : else if (num_in == 1 .and. num_out == 1) then
446 4 : if (nuclides% chem_id(pspecies(1)) == ihe3 .and. &
447 : nuclides% chem_id(pspecies(2)) == ih3) then
448 3 : handle = 'r_he3_ec_h3'
449 3 : return
450 : end if
451 1 : if (nuclides% chem_id(pspecies(1)) == ibe7 .and. &
452 : nuclides% chem_id(pspecies(2)) == ili7) then
453 1 : handle = 'r_be7_wk_li7'
454 1 : return
455 : end if
456 : end if
457 476958 : else if (wk_flag) then
458 43863 : if (reverse) then
459 0 : handle = ''
460 0 : return
461 : end if
462 43863 : if (num_in == 2 .and. num_out == 1) then
463 : if (nuclides% chem_id(pspecies(1)) == ih1 .and. &
464 6 : nuclides% chem_id(pspecies(2)) == ih1 .and. &
465 : nuclides% chem_id(pspecies(3)) == ih2) then
466 3 : handle = 'r_h1_h1_wk_h2'
467 3 : return
468 : end if
469 43857 : else if (num_in == 2 .and. num_out == 1) then
470 : if (nuclides% chem_id(pspecies(1)) == ih1 .and. &
471 0 : nuclides% chem_id(pspecies(2)) == ihe3 .and. &
472 : nuclides% chem_id(pspecies(3)) == ihe4) then
473 0 : handle = 'r_h1_he3_wk_he4'
474 0 : return
475 : end if
476 : end if
477 : end if
478 :
479 476955 : in1 = 0; in2 = 0; out1 = 0; out2 = 0
480 476955 : do_long_form = .true.
481 476955 : if (num_in == 1 .and. num_out == 1) then
482 22169 : call do_n_to_m(1,1)
483 22169 : do_long_form = one_one()
484 454786 : else if (num_in == 1 .and. num_out == 2) then
485 23679 : call do_n_to_m(1,2)
486 23679 : if (reverse) then
487 7356 : do_long_form = two_one()
488 : else
489 16323 : do_long_form = one_two()
490 : end if
491 431107 : else if (num_in == 2 .and. num_out == 1) then
492 185619 : call do_n_to_m(2,1)
493 185619 : if (reverse) then
494 47556 : do_long_form = one_two()
495 : else
496 138063 : do_long_form = two_one()
497 : end if
498 245488 : else if (num_in == 2 .and. num_out == 2) then
499 191111 : call do_n_to_m(2,2)
500 191111 : do_long_form = two_two()
501 : end if
502 :
503 422578 : if (do_long_form) then
504 103360 : call long_form
505 373595 : else if (out2 /= 0) then
506 216026 : handle = trim(handle) // '_' // nuclides% name(out2)
507 : else
508 157569 : handle = trim(handle) // '_' // nuclides% name(out1)
509 : end if
510 :
511 :
512 : contains
513 :
514 953930 : subroutine sort(n, species)
515 : integer :: n
516 : integer :: species(n)
517 : integer :: i, j, Zi, Ni, Zj, Nj, isomer_j, isomer_i, cid
518 : include 'formats'
519 1690831 : do i=1,n-1
520 736901 : cid = species(i)
521 736901 : if (cid <= 0) cycle
522 736901 : Zi = chem_isos% Z(cid)
523 736901 : Ni = chem_isos% N(cid)
524 736901 : isomer_i = chem_isos% isomeric_state(cid)
525 2514951 : do j=i+1,n
526 824120 : cid = species(j)
527 824120 : if (cid <= 0) cycle
528 824120 : Zj = chem_isos% Z(cid)
529 824120 : Nj = chem_isos% N(cid)
530 824120 : isomer_j = chem_isos% isomeric_state(cid)
531 824120 : if (Zj > Zi) cycle
532 147600 : if (Zj == Zi) then
533 147570 : if (Nj > Ni) cycle
534 147516 : if (Nj == Ni) then
535 147498 : if (isomer_j >= isomer_i) cycle
536 : end if
537 : end if
538 : ! exchange i and j
539 48 : species(j) = species(i)
540 48 : species(i) = cid
541 48 : Zi = Zj
542 48 : Ni = Nj
543 1561021 : isomer_i = isomer_j
544 : end do
545 : end do
546 953930 : end subroutine sort
547 :
548 103360 : subroutine long_form
549 : integer :: i, cid
550 : character (len=3) :: op
551 103360 : handle = 'r_'
552 103360 : if (wk_flag) then
553 22719 : op = 'wk_'
554 80641 : else if (ec_flag) then
555 0 : op = 'ec_'
556 : else
557 80641 : op = 'to_'
558 : end if
559 103360 : if (reverse) then
560 753 : do i = num_in+1,num_in+num_out
561 525 : cid = pspecies(i)
562 753 : handle = trim(handle) // trim(nuclides% name(cid)) // '_'
563 : end do
564 228 : handle = trim(handle) // op
565 717 : do i = 1,num_in
566 489 : cid = pspecies(i)
567 489 : handle = trim(handle) // trim(nuclides% name(cid))
568 717 : if (i < num_in) handle = trim(handle) // '_'
569 : end do
570 : else
571 306733 : do i = 1,num_in
572 203601 : cid = pspecies(i)
573 306733 : handle = trim(handle) // trim(nuclides% name(cid)) // '_'
574 : end do
575 103132 : handle = trim(handle) // op
576 329952 : do i = num_in+1,num_in+num_out
577 226820 : cid = pspecies(i)
578 226820 : handle = trim(handle) // trim(nuclides% name(cid))
579 329952 : if (i < num_in+num_out) handle = trim(handle) // '_'
580 : end do
581 : end if
582 103360 : end subroutine long_form
583 :
584 22169 : logical function one_one()
585 22169 : one_one = .true.
586 22169 : if (in1 == 0 .or. out1 == 0) return
587 22169 : one_one = .false.
588 22169 : if (nuclides% Z(out1) == nuclides% Z(in1) - 1 .and. &
589 : nuclides% N(out1) == nuclides% N(in1) + 1) then
590 6425 : handle = trim(handle) // '_wk'
591 15744 : else if (nuclides% Z(out1) == nuclides% Z(in1) + 1 .and. &
592 : nuclides% N(out1) == nuclides% N(in1) - 1) then
593 15744 : handle = trim(handle) // '_wk-minus'
594 : else
595 : one_one = .true.
596 : end if
597 : end function one_one
598 :
599 63879 : logical function one_two()
600 63879 : one_two = .true.
601 63879 : if (in1 == 0 .or. out1 == 0 .or. out2 == 0 .or. out1 == out2) return
602 63870 : one_two = .false.
603 : if (nuclides% Z(out1) == 0 .and. nuclides% N(out1) == 1 .and. &
604 63870 : nuclides% Z(out2) == nuclides% Z(in1) .and. &
605 : nuclides% N(out2) == nuclides% N(in1) - 1) then
606 22065 : handle = trim(handle) // '_gn'
607 : else if (nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &
608 41805 : nuclides% Z(out2) == nuclides% Z(in1) - 1 .and. &
609 : nuclides% N(out2) == nuclides% N(in1)) then
610 15759 : handle = trim(handle) // '_gp'
611 : else if (nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
612 26046 : nuclides% Z(out2) == nuclides% Z(in1) - 2 .and. &
613 : nuclides% N(out2) == nuclides% N(in1) - 2) then
614 17082 : handle = trim(handle) // '_ga'
615 : else if (nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &
616 8964 : nuclides% Z(out2) == nuclides% Z(in1) - 2 .and. &
617 : nuclides% N(out2) == nuclides% N(in1) + 1) then
618 324 : handle = trim(handle) // '_wk_h1'
619 : else if (nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
620 8640 : nuclides% Z(out2) == nuclides% Z(in1) - 3 .and. &
621 : nuclides% N(out2) == nuclides% N(in1) - 1) then
622 42 : handle = trim(handle) // '_wk_he4'
623 : else
624 : one_two = .true.
625 : end if
626 : end function one_two
627 :
628 145419 : logical function two_one()
629 : include 'formats'
630 145419 : two_one = .true.
631 145419 : if (in1 == 0 .or. in2 == 0 .or. out1 == 0 .or. in1 == in2) return
632 135406 : two_one = .false.
633 : if (nuclides% Z(in1) == 0 .and. nuclides% N(in1) == 1 .and. &
634 135406 : nuclides% Z(out1) == nuclides% Z(in2) .and. &
635 : nuclides% N(out1) == nuclides% N(in2) + 1) then
636 22146 : handle = trim(handle) // '_ng'
637 : else if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
638 113260 : nuclides% Z(out1) == nuclides% Z(in2) + 1 .and. &
639 : nuclides% N(out1) == nuclides% N(in2)) then
640 66115 : handle = trim(handle) // '_pg'
641 : else if (nuclides% Z(in1) == 2 .and. nuclides% N(in1) == 2 .and. &
642 47145 : nuclides% Z(out1) == nuclides% Z(in2) + 2 .and. &
643 : nuclides% N(out1) == nuclides% N(in2) + 2) then
644 47139 : handle = trim(handle) // '_ag'
645 : else
646 : two_one = .true.
647 : end if
648 : end function two_one
649 :
650 191111 : logical function two_two()
651 191111 : two_two = .true.
652 :
653 191111 : if (in1 == 0 .or. in2 == 0 .or. out1 == 0 .or. out2 == 0) return
654 :
655 : ! Special case r_li7_pa_he4, this must come first otherwise the out1==out2
656 : ! check will label this rate as a _to_ reaction instead of an _ap reaction
657 : if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
658 : nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
659 : nuclides% Z(out2) == 2 .and. nuclides% N(out2) == 2 .and. &
660 191111 : nuclides% Z(in2) == 3 .and. nuclides% N(in2) == 4) then
661 12 : handle = trim(handle) // '_pa'
662 12 : two_two = .false.
663 12 : return
664 : end if
665 :
666 191099 : if(in1==in2 .or. out1==out2) return
667 :
668 170932 : two_two = .false.
669 : if (nuclides% Z(in1) == 2 .and. nuclides% N(in1) == 2 .and. &
670 : nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &
671 170932 : nuclides% Z(out2) == nuclides% Z(in2) + 1 .and. &
672 : nuclides% N(out2) == nuclides% N(in2) + 2) then
673 25349 : handle = trim(handle) // '_ap'
674 : else if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
675 : nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
676 145583 : nuclides% Z(out2) == nuclides% Z(in2) - 1 .and. &
677 : nuclides% N(out2) == nuclides% N(in2) - 2) then
678 75435 : handle = trim(handle) // '_pa'
679 : else if (nuclides% Z(in1) == 2 .and. nuclides% N(in1) == 2 .and. &
680 : nuclides% Z(out1) == 0 .and. nuclides% N(out1) == 1 .and. &
681 70148 : nuclides% Z(out2) == nuclides% Z(in2) + 2 .and. &
682 : nuclides% N(out2) == nuclides% N(in2) + 1) then
683 14997 : handle = trim(handle) // '_an'
684 : else if (nuclides% Z(in1) == 0 .and. nuclides% N(in1) == 1 .and. &
685 : nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
686 55151 : nuclides% Z(out2) == nuclides% Z(in2) - 2 .and. &
687 : nuclides% N(out2) == nuclides% N(in2) - 1) then
688 15033 : handle = trim(handle) // '_na'
689 : else if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
690 : nuclides% Z(out1) == 0 .and. nuclides% N(out1) == 1 .and. &
691 40118 : nuclides% Z(out2) == nuclides% Z(in2) + 1 .and. &
692 : nuclides% N(out2) == nuclides% N(in2) - 1) then
693 14964 : handle = trim(handle) // '_pn'
694 : else if (nuclides% Z(in1) == 0 .and. nuclides% N(in1) == 1 .and. &
695 : nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &
696 25154 : nuclides% Z(out2) == nuclides% Z(in2) - 1 .and. &
697 : nuclides% N(out2) == nuclides% N(in2) + 1) then
698 14964 : handle = trim(handle) // '_np'
699 : else
700 : two_two = .true.
701 : end if
702 : end function two_two
703 :
704 422578 : subroutine do_n_to_m(n,m)
705 : integer, intent(in) :: n, m ! each is either 1 or 2
706 422578 : in1 = 0; in2 = 0; out1 = 0; out2 = 0
707 422578 : if (.not. reverse) then
708 322399 : in1 = pspecies(1)
709 322399 : if (n == 2) then
710 283907 : in2 = pspecies(2)
711 283907 : if (in2 == 0) then
712 0 : in1 = 0; return
713 : end if
714 283907 : call switch_if_necessary(in1,in2)
715 283907 : handle = 'r_' // nuclides% name(in2)
716 38492 : else if (n == 1) then
717 38492 : handle = 'r_' // nuclides% name(in1)
718 : else
719 0 : in1 = 0
720 0 : return
721 : end if
722 322399 : out1 = pspecies(n+1)
723 322399 : if (m == 2) then
724 162167 : out2 = pspecies(n+2)
725 162167 : call switch_if_necessary(out1,out2)
726 : end if
727 : else
728 100179 : in1 = pspecies(n+1)
729 100179 : if (m == 2) then
730 52623 : in2 = pspecies(n+2)
731 52623 : if (in2 == 0) then
732 0 : in1 = 0; return
733 : end if
734 52623 : call switch_if_necessary(in1,in2)
735 52623 : handle = 'r_' // nuclides% name(in2)
736 47556 : else if (m == 1) then
737 47556 : handle = 'r_' // nuclides% name(in1)
738 : else
739 0 : in1 = 0
740 0 : return
741 : end if
742 100179 : out1 = pspecies(1)
743 100179 : if (n == 2) then
744 92823 : out2 = pspecies(2)
745 92823 : call switch_if_necessary(out1,out2)
746 : end if
747 : end if
748 : end subroutine do_n_to_m
749 :
750 591520 : subroutine switch_if_necessary(iso1,iso2)
751 : integer, intent(inout) :: iso1, iso2
752 : integer :: j
753 591520 : if (nuclides% Z(iso2) == 1 .and. nuclides% N(iso2) == 0) then ! iso2 is ih1
754 18 : j = iso1; iso1 = iso2; iso2 = j; return
755 : end if
756 591502 : if (nuclides% Z(iso2) == 2 .and. nuclides% N(iso2) == 2) then ! iso2 is ihe4
757 10124 : if (nuclides% Z(iso1) == 1 .and. nuclides% N(iso1) == 0) return ! iso1 is ih1
758 102 : j = iso1; iso1 = iso2; iso2 = j; return
759 : end if
760 : end subroutine switch_if_necessary
761 :
762 :
763 : end subroutine get1_reaction_handle
764 :
765 :
766 : end module reaclib_support
|
