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 145554 : 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 145560 : 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 145554 : 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 151929 : 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 546936 : tmp = product(mass(1:Ni))/product(mass(Ni+1:Nt))
214 546936 : rates% inverse_coefficients(1,i) = pow(tmp,1.5d0*(Ni-No))*(product(g(1:Ni))/product(g(Ni+1:Nt)))
215 :
216 : ! -Q/(kB*10**9)
217 293550 : sum1 = sum(winvn% binding_energy(ps(1:Ni)))
218 253386 : sum2 = sum(winvn% binding_energy(ps(Ni+1:Nt)))
219 100290 : rates% inverse_coefficients(2,i) = (sum1-sum2)*conv/kB/1d9
220 :
221 : ! see equation 21 in Reichert et al. 2023 (https://doi.org/10.3847/1538-4365/acf033)
222 : ! delta_reactants/delta_products is handled in net.
223 : ! fac == (mu * kb * T / (2 * pi * hbar^2))^3/2 term with out a T^3/2.
224 : ! fac shows up as fac^(Ni-No) in rates% inverse_coefficients(1,i)
225 : ! so rates% inverse_coefficients(1,i) contains terms for
226 : ! fac^(n) == fac^(Ni-No), where n = Ni - No.
227 :
228 : ! The T makes its way back into our expression inside
229 : ! the subroutine compute_some_inverse_lambdas, in reaclib_eval.f90.
230 : ! It appears in log form as 1.5d0*rates% inverse_exp(i)*lnT9, where,
231 : ! rates% inverse_exp(i) = Ni-No, so,
232 : ! 1.5d0*rates% inverse_exp(i)*lnT9 == ln(T^(3n/2)), where n = Ni-No.
233 100290 : rates% inverse_exp(i) = Ni - No ! We use this term in a log() expression in reaclib_eval
234 :
235 : ! Ni-No should be 0 for non-photo-disintegration reverse rates
236 : ! and >=1 for photos's in the reverse channel
237 100290 : if (Ni-No /= 0) then
238 : ! whether endothermic or exothermic,
239 : ! Ni-No handles the sign of Q from rates% inverse_coefficients(2,i)
240 55023 : rates% inverse_coefficients(1,i) = rates% inverse_coefficients(1,i) * pow(fac, Ni-No)
241 : ! negative values of Q denote endothermic photodisintegrations
242 : ! We multiply by fac^|Ni-No| as we want to compute
243 : ! rate_photo/rate_forward
244 :
245 : ! positive values of Q denote exothermic photodisintegrations
246 : ! We DIVIDE by fac^|Ni-No| as we want to compute
247 : ! rate_reverse/rate_photo
248 : else ! Ni - No = 0
249 45267 : rates% inverse_exp(i) = 0 ! ensure this is 0.
250 : rates% inverse_coefficients(1,i) = rates% inverse_coefficients(1,i) ! no point calling pow(fac,0)
251 : end if
252 100290 : rates% inverse_coefficients(1,i) = log(rates% inverse_coefficients(1,i))
253 :
254 :
255 100290 : rates% inverse_part(:,i) = product(winvn% pfcn(1:npart,ps(1:Ni)),dim=2)/ &
256 10920087 : & product(winvn% pfcn(1:npart,ps(Ni+1:Nt)),dim=2)
257 : end do loop_over_rates
258 3 : end subroutine compute_rev_ratio
259 :
260 0 : subroutine do_parse_reaction_handle(handle, num_in, num_out, iso_ids, op, ierr)
261 : use chem_def
262 : use chem_lib
263 : character (len=*), intent(in) :: handle
264 : integer, intent(out) :: num_in, num_out
265 : integer, intent(out) :: iso_ids(:) ! holds chem_ids for input and output species
266 : character (len=*), intent(out) :: op ! e.g., 'pg', 'wk', 'to', or ...
267 : integer, intent(out) :: ierr
268 :
269 : integer :: len, i, j, cnt, cid, extra_in, extra_out
270 : logical :: doing_inputs
271 :
272 0 : num_in = 0; num_out = 0; op = ''
273 0 : ierr = -1
274 0 : len = len_trim(handle)
275 0 : if (handle(1:2) /= 'r_') return
276 0 : i = 3
277 0 : cnt = 0
278 0 : doing_inputs = .true.
279 0 : do while (i <= len)
280 0 : call nxt ! set j to last char of token
281 0 : cid = chem_get_iso_id(handle(i:j))
282 0 : if (cid == nuclide_not_found) then
283 0 : if (doing_inputs) then
284 0 : op = handle(i:j)
285 0 : extra_in = -1
286 0 : extra_out = -1
287 0 : if (j == i+1) then ! check 2 character ops
288 0 : select case (op(1:1))
289 : case ('p')
290 0 : extra_in = ih1
291 : case ('a')
292 0 : extra_in = ihe4
293 : case ('n')
294 0 : extra_in = ineut
295 : case ('g')
296 0 : extra_in = 0
297 : case default
298 : end select
299 0 : if (extra_in >= 0) then
300 0 : if (extra_in > 0) then
301 0 : cnt = cnt+1
302 0 : if (cnt /= 2) then
303 : !write(*,*) 'failed to parse ' // &
304 : ! trim(handle) // ' -- problem with ' // handle(i:j)
305 : return
306 : end if
307 0 : if (chem_isos% Z(iso_ids(1)) >= chem_isos% Z(extra_in)) then
308 0 : iso_ids(2) = iso_ids(1)
309 0 : iso_ids(1) = extra_in
310 : else
311 0 : iso_ids(2) = extra_in
312 : end if
313 : end if
314 0 : select case (op(2:2))
315 : case ('p')
316 0 : extra_out = ih1
317 : case ('a')
318 0 : extra_out = ihe4
319 : case ('n')
320 0 : extra_out = ineut
321 : case ('g')
322 : extra_out = 0
323 : case default
324 : !write(*,*) 'failed to parse ' // &
325 : ! trim(handle) // ' -- problem with ' // handle(i:j)
326 0 : return
327 : end select
328 : end if
329 : end if
330 0 : num_in = cnt
331 0 : doing_inputs = .false.
332 0 : if (extra_out > 0) then
333 0 : cnt = cnt+1
334 0 : iso_ids(cnt) = extra_out
335 : end if
336 : else
337 : !write(*,*) 'failed to parse ' // &
338 : ! trim(handle) // ' -- problem with ' // handle(i:j)
339 : return
340 : end if
341 : else
342 0 : cnt = cnt+1
343 0 : iso_ids(cnt) = cid
344 : end if
345 0 : i = j+2
346 : end do
347 0 : num_out = cnt - num_in
348 0 : ierr = 0
349 :
350 : contains
351 :
352 0 : subroutine nxt
353 0 : j = i
354 : do
355 0 : if (j >= len) return
356 0 : j = j+1
357 0 : if (handle(j:j) == '_') then
358 0 : j = j-1; return
359 : end if
360 : end do
361 0 : end subroutine nxt
362 :
363 :
364 : end subroutine do_parse_reaction_handle
365 :
366 231121 : subroutine reaction_handle(num_in, num_out, iso_ids, reaction_flag, handle)
367 : use chem_def, only: chem_isos
368 : integer, intent(in) :: num_in, num_out
369 : integer, intent(in) :: iso_ids(:)
370 : character (len=*), intent(in) :: reaction_flag
371 : character (len=*), intent(out) :: handle
372 : logical, parameter :: reverse = .false.
373 231121 : call get1_reaction_handle(num_in, num_out, iso_ids, chem_isos, reverse, reaction_flag, handle)
374 231121 : end subroutine reaction_handle
375 :
376 0 : subroutine reverse_reaction_handle(num_in, num_out, iso_ids, handle)
377 231121 : use chem_def, only: chem_isos
378 : integer, intent(in) :: num_in, num_out
379 : integer, intent(in) :: iso_ids(:)
380 : character (len=*), intent(out) :: handle
381 : logical, parameter :: reverse = .true.
382 : character (len=1), parameter :: reaction_flag = '-'
383 0 : call get1_reaction_handle(num_in, num_out, iso_ids, chem_isos, reverse, reaction_flag, handle)
384 0 : end subroutine reverse_reaction_handle
385 :
386 145554 : subroutine get_reaction_handle(num_in, num_out, pspecies, nuclides, reaction_flag, handle)
387 : integer, intent(in) :: num_in, num_out
388 : integer, intent(in) :: pspecies(:)
389 : type(nuclide_data), intent(in) :: nuclides
390 : character (len=*), intent(in) :: reaction_flag
391 : character (len=*), intent(out) :: handle
392 : logical, parameter :: reverse = .false.
393 145554 : call get1_reaction_handle(num_in, num_out, pspecies, nuclides, reverse, reaction_flag, handle)
394 0 : end subroutine get_reaction_handle
395 :
396 100290 : subroutine get_reverse_reaction_handle(num_in, num_out, pspecies, nuclides, handle)
397 : integer, intent(in) :: num_in, num_out
398 : integer, intent(in) :: pspecies(:)
399 : type(nuclide_data), intent(in) :: nuclides
400 : character (len=*), intent(out) :: handle
401 : logical, parameter :: reverse = .true.
402 : character (len=1), parameter :: reaction_flag = '-'
403 100290 : call get1_reaction_handle(num_in, num_out, pspecies, nuclides, reverse, reaction_flag, handle)
404 100290 : end subroutine get_reverse_reaction_handle
405 :
406 476965 : subroutine get1_reaction_handle( &
407 476965 : num_in, num_out, pspecies_in, nuclides, reverse, reaction_flag, handle)
408 : use chem_def, only: ih1, ih2, ih3, ihe3, ihe4, ibe7, ili7, chem_isos
409 : integer, intent(in) :: num_in, num_out
410 : integer, intent(in) :: pspecies_in(:)
411 : type(nuclide_data), intent(in) :: nuclides
412 : logical, intent(in) :: reverse
413 : character (len=*), intent(in) :: reaction_flag
414 : character (len=*), intent(out) :: handle
415 :
416 953930 : integer :: in1, in2, out1, out2, num, pspecies(num_in + num_out)
417 : logical :: do_long_form, ec_flag, wk_flag
418 :
419 : include 'formats'
420 :
421 476965 : num = num_in + num_out
422 2167796 : pspecies(1:num) = pspecies_in(1:num)
423 476965 : call sort(num_in, pspecies(1:num_in))
424 476965 : call sort(num_out, pspecies(num_in+1:num))
425 476965 : ec_flag = (reaction_flag == 'e')
426 476965 : wk_flag = (reaction_flag == 'w')
427 :
428 476965 : if (ec_flag) then ! special cases
429 7 : if (reverse) then
430 0 : handle = ''
431 0 : return
432 : end if
433 7 : if (num_in == 2 .and. num_out == 1) then
434 6 : if (nuclides% chem_id(pspecies(1)) == ih1 .and. &
435 6 : nuclides% chem_id(pspecies(2)) == ih1 .and. &
436 3 : nuclides% chem_id(pspecies(3)) == ih2) then
437 3 : handle = 'r_h1_h1_ec_h2'
438 3 : return
439 : end if
440 4 : else if (num_in == 1 .and. num_out == 1) then
441 8 : if (nuclides% chem_id(pspecies(1)) == ihe3 .and. &
442 4 : nuclides% chem_id(pspecies(2)) == ih3) then
443 3 : handle = 'r_he3_ec_h3'
444 3 : return
445 : end if
446 1 : if (nuclides% chem_id(pspecies(1)) == ibe7 .and. &
447 : nuclides% chem_id(pspecies(2)) == ili7) then
448 1 : handle = 'r_be7_wk_li7'
449 1 : return
450 : end if
451 : end if
452 476958 : else if (wk_flag) then
453 43863 : if (reverse) then
454 0 : handle = ''
455 0 : return
456 : end if
457 43863 : if (num_in == 2 .and. num_out == 1) then
458 12 : if (nuclides% chem_id(pspecies(1)) == ih1 .and. &
459 12 : nuclides% chem_id(pspecies(2)) == ih1 .and. &
460 6 : nuclides% chem_id(pspecies(3)) == ih2) then
461 3 : handle = 'r_h1_h1_wk_h2'
462 3 : return
463 : end if
464 43857 : else if (num_in == 2 .and. num_out == 1) then
465 0 : if (nuclides% chem_id(pspecies(1)) == ih1 .and. &
466 0 : nuclides% chem_id(pspecies(2)) == ihe3 .and. &
467 0 : nuclides% chem_id(pspecies(3)) == ihe4) then
468 0 : handle = 'r_h1_he3_wk_he4'
469 0 : return
470 : end if
471 : end if
472 : end if
473 :
474 476955 : in1 = 0; in2 = 0; out1 = 0; out2 = 0
475 476955 : do_long_form = .true.
476 476955 : if (num_in == 1 .and. num_out == 1) then
477 22169 : call do_n_to_m(1,1)
478 22169 : do_long_form = one_one()
479 454786 : else if (num_in == 1 .and. num_out == 2) then
480 23679 : call do_n_to_m(1,2)
481 23679 : if (reverse) then
482 7356 : do_long_form = two_one()
483 : else
484 16323 : do_long_form = one_two()
485 : end if
486 431107 : else if (num_in == 2 .and. num_out == 1) then
487 185619 : call do_n_to_m(2,1)
488 185619 : if (reverse) then
489 47556 : do_long_form = one_two()
490 : else
491 138063 : do_long_form = two_one()
492 : end if
493 245488 : else if (num_in == 2 .and. num_out == 2) then
494 191111 : call do_n_to_m(2,2)
495 191111 : do_long_form = two_two()
496 : end if
497 :
498 899543 : if (do_long_form) then
499 103360 : call long_form
500 373595 : else if (out2 /= 0) then
501 216026 : handle = trim(handle) // '_' // nuclides% name(out2)
502 : else
503 157569 : handle = trim(handle) // '_' // nuclides% name(out1)
504 : end if
505 :
506 :
507 : contains
508 :
509 953930 : subroutine sort(n, species)
510 : integer :: n
511 : integer :: species(n)
512 : integer :: i, j, Zi, Ni, Zj, Nj, isomer_j, isomer_i, cid
513 : include 'formats'
514 1690831 : do i=1,n-1
515 736901 : cid = species(i)
516 736901 : if (cid <= 0) cycle
517 736901 : Zi = chem_isos% Z(cid)
518 736901 : Ni = chem_isos% N(cid)
519 736901 : isomer_i = chem_isos% isomeric_state(cid)
520 2514951 : do j=i+1,n
521 824120 : cid = species(j)
522 824120 : if (cid <= 0) cycle
523 824120 : Zj = chem_isos% Z(cid)
524 824120 : Nj = chem_isos% N(cid)
525 824120 : isomer_j = chem_isos% isomeric_state(cid)
526 824120 : if (Zj > Zi) cycle
527 147600 : if (Zj == Zi) then
528 147570 : if (Nj > Ni) cycle
529 147516 : if (Nj == Ni) then
530 147498 : if (isomer_j >= isomer_i) cycle
531 : end if
532 : end if
533 : ! exchange i and j
534 48 : species(j) = species(i)
535 48 : species(i) = cid
536 48 : Zi = Zj
537 48 : Ni = Nj
538 1561021 : isomer_i = isomer_j
539 : end do
540 : end do
541 476965 : end subroutine sort
542 :
543 103360 : subroutine long_form
544 : integer :: i, cid
545 : character (len=3) :: op
546 103360 : handle = 'r_'
547 103360 : if (wk_flag) then
548 22719 : op = 'wk_'
549 80641 : else if (ec_flag) then
550 0 : op = 'ec_'
551 : else
552 80641 : op = 'to_'
553 : end if
554 103360 : if (reverse) then
555 753 : do i = num_in+1,num_in+num_out
556 525 : cid = pspecies(i)
557 753 : handle = trim(handle) // trim(nuclides% name(cid)) // '_'
558 : end do
559 228 : handle = trim(handle) // op
560 717 : do i = 1,num_in
561 489 : cid = pspecies(i)
562 489 : handle = trim(handle) // trim(nuclides% name(cid))
563 717 : if (i < num_in) handle = trim(handle) // '_'
564 : end do
565 : else
566 306733 : do i = 1,num_in
567 203601 : cid = pspecies(i)
568 306733 : handle = trim(handle) // trim(nuclides% name(cid)) // '_'
569 : end do
570 103132 : handle = trim(handle) // op
571 329952 : do i = num_in+1,num_in+num_out
572 226820 : cid = pspecies(i)
573 226820 : handle = trim(handle) // trim(nuclides% name(cid))
574 329952 : if (i < num_in+num_out) handle = trim(handle) // '_'
575 : end do
576 : end if
577 103360 : end subroutine long_form
578 :
579 22169 : logical function one_one()
580 22169 : one_one = .true.
581 22169 : if (in1 == 0 .or. out1 == 0) return
582 22169 : one_one = .false.
583 44338 : if (nuclides% Z(out1) == nuclides% Z(in1) - 1 .and. &
584 22169 : nuclides% N(out1) == nuclides% N(in1) + 1) then
585 6425 : handle = trim(handle) // '_wk'
586 15744 : else if (nuclides% Z(out1) == nuclides% Z(in1) + 1 .and. &
587 : nuclides% N(out1) == nuclides% N(in1) - 1) then
588 15744 : handle = trim(handle) // '_wk-minus'
589 : else
590 : one_one = .true.
591 : end if
592 : end function one_one
593 :
594 63879 : logical function one_two()
595 63879 : one_two = .true.
596 63879 : if (in1 == 0 .or. out1 == 0 .or. out2 == 0 .or. out1 == out2) return
597 63870 : one_two = .false.
598 63870 : if (nuclides% Z(out1) == 0 .and. nuclides% N(out1) == 1 .and. &
599 191610 : nuclides% Z(out2) == nuclides% Z(in1) .and. &
600 63870 : nuclides% N(out2) == nuclides% N(in1) - 1) then
601 22065 : handle = trim(handle) // '_gn'
602 : else if (nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &
603 41805 : nuclides% Z(out2) == nuclides% Z(in1) - 1 .and. &
604 : nuclides% N(out2) == nuclides% N(in1)) then
605 15759 : handle = trim(handle) // '_gp'
606 : else if (nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
607 26046 : nuclides% Z(out2) == nuclides% Z(in1) - 2 .and. &
608 : nuclides% N(out2) == nuclides% N(in1) - 2) then
609 17082 : handle = trim(handle) // '_ga'
610 : else if (nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &
611 8964 : nuclides% Z(out2) == nuclides% Z(in1) - 2 .and. &
612 : nuclides% N(out2) == nuclides% N(in1) + 1) then
613 324 : handle = trim(handle) // '_wk_h1'
614 : else if (nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
615 8640 : nuclides% Z(out2) == nuclides% Z(in1) - 3 .and. &
616 : nuclides% N(out2) == nuclides% N(in1) - 1) then
617 42 : handle = trim(handle) // '_wk_he4'
618 : else
619 : one_two = .true.
620 : end if
621 : end function one_two
622 :
623 145419 : logical function two_one()
624 : include 'formats'
625 145419 : two_one = .true.
626 145419 : if (in1 == 0 .or. in2 == 0 .or. out1 == 0 .or. in1 == in2) return
627 135406 : two_one = .false.
628 135406 : if (nuclides% Z(in1) == 0 .and. nuclides% N(in1) == 1 .and. &
629 406218 : nuclides% Z(out1) == nuclides% Z(in2) .and. &
630 135406 : nuclides% N(out1) == nuclides% N(in2) + 1) then
631 22146 : handle = trim(handle) // '_ng'
632 : else if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
633 113260 : nuclides% Z(out1) == nuclides% Z(in2) + 1 .and. &
634 : nuclides% N(out1) == nuclides% N(in2)) then
635 66115 : handle = trim(handle) // '_pg'
636 : else if (nuclides% Z(in1) == 2 .and. nuclides% N(in1) == 2 .and. &
637 47145 : nuclides% Z(out1) == nuclides% Z(in2) + 2 .and. &
638 : nuclides% N(out1) == nuclides% N(in2) + 2) then
639 47139 : handle = trim(handle) // '_ag'
640 : else
641 : two_one = .true.
642 : end if
643 : end function two_one
644 :
645 191111 : logical function two_two()
646 191111 : two_two = .true.
647 :
648 191111 : if (in1 == 0 .or. in2 == 0 .or. out1 == 0 .or. out2 == 0) return
649 :
650 : ! Special case r_li7_pa_he4, this must come first otherwise the out1==out2
651 : ! check will label this rate as a _to_ reaction instead of an _ap reaction
652 191111 : if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
653 382222 : nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
654 382222 : nuclides% Z(out2) == 2 .and. nuclides% N(out2) == 2 .and. &
655 764444 : nuclides% Z(in2) == 3 .and. nuclides% N(in2) == 4) then
656 12 : handle = trim(handle) // '_pa'
657 12 : two_two = .false.
658 12 : return
659 : end if
660 :
661 191099 : if(in1==in2 .or. out1==out2) return
662 :
663 170932 : two_two = .false.
664 : if (nuclides% Z(in1) == 2 .and. nuclides% N(in1) == 2 .and. &
665 : nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &
666 170932 : nuclides% Z(out2) == nuclides% Z(in2) + 1 .and. &
667 : nuclides% N(out2) == nuclides% N(in2) + 2) then
668 25349 : handle = trim(handle) // '_ap'
669 : else if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
670 : nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
671 145583 : nuclides% Z(out2) == nuclides% Z(in2) - 1 .and. &
672 : nuclides% N(out2) == nuclides% N(in2) - 2) then
673 75435 : handle = trim(handle) // '_pa'
674 : else if (nuclides% Z(in1) == 2 .and. nuclides% N(in1) == 2 .and. &
675 : nuclides% Z(out1) == 0 .and. nuclides% N(out1) == 1 .and. &
676 70148 : nuclides% Z(out2) == nuclides% Z(in2) + 2 .and. &
677 : nuclides% N(out2) == nuclides% N(in2) + 1) then
678 14997 : handle = trim(handle) // '_an'
679 : else if (nuclides% Z(in1) == 0 .and. nuclides% N(in1) == 1 .and. &
680 : nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
681 55151 : nuclides% Z(out2) == nuclides% Z(in2) - 2 .and. &
682 : nuclides% N(out2) == nuclides% N(in2) - 1) then
683 15033 : handle = trim(handle) // '_na'
684 : else if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
685 : nuclides% Z(out1) == 0 .and. nuclides% N(out1) == 1 .and. &
686 40118 : nuclides% Z(out2) == nuclides% Z(in2) + 1 .and. &
687 : nuclides% N(out2) == nuclides% N(in2) - 1) then
688 14964 : handle = trim(handle) // '_pn'
689 : else if (nuclides% Z(in1) == 0 .and. nuclides% N(in1) == 1 .and. &
690 : nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &
691 25154 : nuclides% Z(out2) == nuclides% Z(in2) - 1 .and. &
692 : nuclides% N(out2) == nuclides% N(in2) + 1) then
693 14964 : handle = trim(handle) // '_np'
694 : else
695 : two_two = .true.
696 : end if
697 : end function two_two
698 :
699 422578 : subroutine do_n_to_m(n,m)
700 : integer, intent(in) :: n, m ! each is either 1 or 2
701 422578 : in1 = 0; in2 = 0; out1 = 0; out2 = 0
702 422578 : if (.not. reverse) then
703 322399 : in1 = pspecies(1)
704 322399 : if (n == 2) then
705 283907 : in2 = pspecies(2)
706 283907 : if (in2 == 0) then
707 0 : in1 = 0; return
708 : end if
709 283907 : call switch_if_necessary(in1,in2)
710 283907 : handle = 'r_' // nuclides% name(in2)
711 38492 : else if (n == 1) then
712 38492 : handle = 'r_' // nuclides% name(in1)
713 : else
714 0 : in1 = 0
715 0 : return
716 : end if
717 322399 : out1 = pspecies(n+1)
718 322399 : if (m == 2) then
719 162167 : out2 = pspecies(n+2)
720 162167 : call switch_if_necessary(out1,out2)
721 : end if
722 : else
723 100179 : in1 = pspecies(n+1)
724 100179 : if (m == 2) then
725 52623 : in2 = pspecies(n+2)
726 52623 : if (in2 == 0) then
727 0 : in1 = 0; return
728 : end if
729 52623 : call switch_if_necessary(in1,in2)
730 52623 : handle = 'r_' // nuclides% name(in2)
731 47556 : else if (m == 1) then
732 47556 : handle = 'r_' // nuclides% name(in1)
733 : else
734 0 : in1 = 0
735 0 : return
736 : end if
737 100179 : out1 = pspecies(1)
738 100179 : if (n == 2) then
739 92823 : out2 = pspecies(2)
740 92823 : call switch_if_necessary(out1,out2)
741 : end if
742 : end if
743 : end subroutine do_n_to_m
744 :
745 591520 : subroutine switch_if_necessary(iso1,iso2)
746 : integer, intent(inout) :: iso1, iso2
747 : integer :: j
748 591520 : if (nuclides% Z(iso2) == 1 .and. nuclides% N(iso2) == 0) then ! iso2 is ih1
749 18 : j = iso1; iso1 = iso2; iso2 = j; return
750 : end if
751 591502 : if (nuclides% Z(iso2) == 2 .and. nuclides% N(iso2) == 2) then ! iso2 is ihe4
752 10124 : if (nuclides% Z(iso1) == 1 .and. nuclides% N(iso1) == 0) return ! iso1 is ih1
753 102 : j = iso1; iso1 = iso2; iso2 = j; return
754 : end if
755 : end subroutine switch_if_necessary
756 :
757 :
758 : end subroutine get1_reaction_handle
759 :
760 :
761 : end module reaclib_support
|
