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 5 : subroutine set_up_network_information(rates)
31 : type(reaction_data), intent(inout) :: rates
32 : integer :: current_chapter, i
33 :
34 : ! set up bookmarks
35 5 : current_chapter = 0
36 5 : rates% nchapters_included = 0
37 60 : rates% chapters_present = 0
38 170 : rates% bookmarks = 0
39 242595 : do i = 1, rates% nreactions
40 242595 : new_chapter : if (rates% chapter(i) /= current_chapter) then
41 : ! close out the chapter we just left.
42 55 : if (current_chapter /= 0) rates% bookmarks(2,current_chapter) = i-1
43 : ! set up information on the new chapter
44 55 : current_chapter = rates% chapter(i)
45 55 : rates% nchapters_included = rates% nchapters_included + 1
46 55 : rates% chapters_present(rates% nchapters_included) = current_chapter
47 55 : rates% bookmarks(1,current_chapter) = i
48 : end if new_chapter
49 : end do
50 : ! mark the end of the last chapter
51 5 : rates% bookmarks(2,current_chapter) = rates% nreactions
52 5 : end subroutine set_up_network_information
53 :
54 :
55 5 : 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 5 : if (.not.associated(rates% weight)) then
63 : return
64 : end if
65 :
66 242595 : do i = 1, rates% nreactions
67 242590 : i1 = -1; i2 = -2; i3 = -3; i4 = -4
68 321865 : 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 79275 : i1 = rates% pspecies(1,i)
74 79275 : i2 = rates% pspecies(2,i)
75 : case (r_two_two)
76 75445 : i1 = rates% pspecies(1,i)
77 75445 : i2 = rates% pspecies(2,i)
78 : case (r_two_three)
79 115 : i1 = rates% pspecies(1,i)
80 115 : i2 = rates% pspecies(2,i)
81 : case (r_two_four)
82 20 : i1 = rates% pspecies(1,i)
83 20 : i2 = rates% pspecies(2,i)
84 : case (r_three_one)
85 35 : i1 = rates% pspecies(1,i)
86 35 : i2 = rates% pspecies(2,i)
87 35 : i3 = rates% pspecies(3,i)
88 : case (r_three_two)
89 5 : i1 = rates% pspecies(1,i)
90 5 : i2 = rates% pspecies(2,i)
91 5 : i3 = rates% pspecies(3,i)
92 : case (r_four_two)
93 10 : i1 = rates% pspecies(1,i)
94 10 : i2 = rates% pspecies(2,i)
95 10 : i3 = rates% pspecies(3,i)
96 242590 : i4 = rates% pspecies(4,i)
97 : case (r_one_four)
98 : end select
99 242595 : call set_weight(rates% weight(i))
100 : end do
101 :
102 242595 : do i = 1, rates% nreactions
103 242590 : i1 = -1; i2 = -2; i3 = -3; i4 = -4
104 269795 : select case (rates% chapter(i))
105 : case (r_one_one)
106 : case (r_one_two)
107 27205 : i1 = rates% pspecies(2,i)
108 27205 : i2 = rates% pspecies(3,i)
109 : case (r_one_three)
110 12900 : i1 = rates% pspecies(2,i)
111 12900 : i2 = rates% pspecies(3,i)
112 12900 : i3 = rates% pspecies(4,i)
113 : case (r_two_one)
114 : case (r_two_two)
115 75445 : i1 = rates% pspecies(3,i)
116 75445 : i2 = rates% pspecies(4,i)
117 : case (r_two_three)
118 115 : i1 = rates% pspecies(3,i)
119 115 : i2 = rates% pspecies(4,i)
120 115 : i3 = rates% pspecies(5,i)
121 : case (r_two_four)
122 20 : i1 = rates% pspecies(3,i)
123 20 : i2 = rates% pspecies(4,i)
124 20 : i3 = rates% pspecies(5,i)
125 20 : i4 = rates% pspecies(6,i)
126 : case (r_three_one)
127 : case (r_three_two)
128 5 : i1 = rates% pspecies(4,i)
129 5 : i2 = rates% pspecies(5,i)
130 : case (r_four_two)
131 10 : i1 = rates% pspecies(5,i)
132 10 : i2 = rates% pspecies(6,i)
133 : case (r_one_four)
134 10625 : i1 = rates% pspecies(2,i)
135 10625 : i2 = rates% pspecies(3,i)
136 10625 : i3 = rates% pspecies(4,i)
137 242590 : i4 = rates% pspecies(5,i)
138 : end select
139 :
140 242595 : call set_weight(rates% weight_reverse(i))
141 :
142 : end do
143 :
144 :
145 : contains
146 :
147 :
148 485180 : subroutine set_weight(w)
149 : ! nuclei are sorted, so if identical, then are adjacent in list
150 : real(dp), intent(out) :: w
151 485180 : if (i1 == i2 .and. i2 == i3 .and. i3 == i4) then
152 0 : w = 1d0/24d0
153 485180 : else if (i2 == i3 .and. (i1 == i2 .or. i3 == i4)) then
154 10665 : w = 1d0/6d0
155 474515 : else if (i1 == i2) then
156 13050 : if (i3 == i4) then
157 10 : w = 1d0/4d0
158 : else
159 13040 : w = 1d0/2d0
160 : end if
161 461465 : else if (i2 == i3 .or. i3 == i4) then
162 90 : w = 1d0/2d0
163 : else
164 461375 : w = 1d0
165 : end if
166 485180 : end subroutine set_weight
167 :
168 :
169 : end subroutine assign_weights
170 :
171 :
172 5 : 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 5 : real(dp) :: mp, mn
178 35 : real(dp), dimension(max_species_per_reaction) :: g ! statistical weights of nuclides
179 35 : 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 5 : real(dp) :: fac, massfac, sum1, sum2, tmp
183 :
184 :
185 : include 'formats'
186 :
187 : ! Get these consistently from the isotopes.data file
188 5 : mp=winvn%W(chem_get_iso_id('prot'))
189 5 : mn=winvn%W(chem_get_iso_id('neut'))
190 :
191 5 : fac = pow(1d9*kB/(2d0*pi*hbar*hbar*NA),1.5d0)/NA
192 5 : massfac = conv*NA/(c*c)
193 :
194 242595 : rates% weak_mask = 1d0
195 242595 : loop_over_rates: do i = 1,rates% nreactions
196 : ! weak rates don't have inverses, so set to innocuous values and mask out
197 242590 : if (rates% reaction_flag(i) == 'w' .or. rates% reaction_flag(i) == 'e') then
198 75440 : rates% weak_mask(i) = 0d0
199 226320 : rates% inverse_coefficients(:,i) = [-huge(1d0), 0d0]
200 75440 : rates% inverse_exp(i) = 0d0
201 1886000 : rates% inverse_part(:,i) = 1d0
202 : cycle loop_over_rates
203 : end if
204 167150 : Ni = Nin(rates% chapter(i))
205 167150 : No = Nout(rates% chapter(i))
206 167150 : Nt = Ni+No
207 744410 : ps(1:Nt) = rates% pspecies(1:Nt,i)
208 744410 : 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 744410 : mass(1:Nt) = winvn% W(ps(1:Nt))
211 :
212 : ! log(prefactor of reverse_ratio)
213 744410 : tmp = product(mass(1:Ni))/product(mass(Ni+1:Nt))
214 911560 : 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 489250 : sum1 = sum(winvn% binding_energy(ps(1:Ni)))
218 422310 : sum2 = sum(winvn% binding_energy(ps(Ni+1:Nt)))
219 167150 : 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 167150 : 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 167150 : 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 91705 : 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 75445 : 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 167150 : rates% inverse_coefficients(1,i) = log(rates% inverse_coefficients(1,i))
253 :
254 :
255 : rates% inverse_part(:,i) = product(winvn% pfcn(1:npart,ps(1:Ni)),dim=2)/ &
256 18032995 : & product(winvn% pfcn(1:npart,ps(Ni+1:Nt)),dim=2)
257 : end do loop_over_rates
258 5 : 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 232219 : 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 232219 : call get1_reaction_handle(num_in, num_out, iso_ids, chem_isos, reverse, reaction_flag, handle)
374 232219 : end subroutine reaction_handle
375 :
376 0 : subroutine reverse_reaction_handle(num_in, num_out, iso_ids, handle)
377 232219 : 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 242590 : 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 242590 : call get1_reaction_handle(num_in, num_out, pspecies, nuclides, reverse, reaction_flag, handle)
394 0 : end subroutine get_reaction_handle
395 :
396 167150 : 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 167150 : call get1_reaction_handle(num_in, num_out, pspecies, nuclides, reverse, reaction_flag, handle)
404 167150 : end subroutine get_reverse_reaction_handle
405 :
406 641959 : subroutine get1_reaction_handle( &
407 641959 : 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 1283918 : 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 641959 : num = num_in + num_out
422 2887928 : pspecies(1:num) = pspecies_in(1:num)
423 641959 : call sort(num_in, pspecies(1:num_in))
424 641959 : call sort(num_out, pspecies(num_in+1:num))
425 641959 : ec_flag = (reaction_flag == 'e')
426 641959 : wk_flag = (reaction_flag == 'w')
427 :
428 641959 : if (ec_flag) then ! special cases
429 12 : if (reverse) then
430 0 : handle = ''
431 0 : return
432 : end if
433 12 : if (num_in == 2 .and. num_out == 1) then
434 : if (nuclides% chem_id(pspecies(1)) == ih1 .and. &
435 5 : nuclides% chem_id(pspecies(2)) == ih1 .and. &
436 : nuclides% chem_id(pspecies(3)) == ih2) then
437 5 : handle = 'r_h1_h1_ec_h2'
438 5 : return
439 : end if
440 7 : else if (num_in == 1 .and. num_out == 1) then
441 7 : if (nuclides% chem_id(pspecies(1)) == ihe3 .and. &
442 : nuclides% chem_id(pspecies(2)) == ih3) then
443 5 : handle = 'r_he3_ec_h3'
444 5 : return
445 : end if
446 2 : if (nuclides% chem_id(pspecies(1)) == ibe7 .and. &
447 : nuclides% chem_id(pspecies(2)) == ili7) then
448 2 : handle = 'r_be7_wk_li7'
449 2 : return
450 : end if
451 : end if
452 641947 : else if (wk_flag) then
453 73105 : if (reverse) then
454 0 : handle = ''
455 0 : return
456 : end if
457 73105 : if (num_in == 2 .and. num_out == 1) then
458 : if (nuclides% chem_id(pspecies(1)) == ih1 .and. &
459 10 : nuclides% chem_id(pspecies(2)) == ih1 .and. &
460 : nuclides% chem_id(pspecies(3)) == ih2) then
461 5 : handle = 'r_h1_h1_wk_h2'
462 5 : return
463 : end if
464 73095 : else if (num_in == 2 .and. num_out == 1) then
465 : if (nuclides% chem_id(pspecies(1)) == ih1 .and. &
466 0 : nuclides% chem_id(pspecies(2)) == ihe3 .and. &
467 : 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 641942 : in1 = 0; in2 = 0; out1 = 0; out2 = 0
475 641942 : do_long_form = .true.
476 641942 : if (num_in == 1 .and. num_out == 1) then
477 36948 : call do_n_to_m(1,1)
478 36948 : do_long_form = one_one()
479 604994 : else if (num_in == 1 .and. num_out == 2) then
480 39465 : call do_n_to_m(1,2)
481 39465 : if (reverse) then
482 12260 : do_long_form = two_one()
483 : else
484 27205 : do_long_form = one_two()
485 : end if
486 565529 : else if (num_in == 2 .and. num_out == 1) then
487 249524 : call do_n_to_m(2,1)
488 249524 : if (reverse) then
489 79260 : do_long_form = one_two()
490 : else
491 170264 : do_long_form = two_one()
492 : end if
493 316005 : else if (num_in == 2 .and. num_out == 2) then
494 252034 : call do_n_to_m(2,2)
495 252034 : do_long_form = two_two()
496 : end if
497 :
498 1219930 : if (do_long_form) then
499 118991 : call long_form
500 522951 : else if (out2 /= 0) then
501 313513 : handle = trim(handle) // '_' // nuclides% name(out2)
502 : else
503 209438 : handle = trim(handle) // '_' // nuclides% name(out1)
504 : end if
505 :
506 :
507 : contains
508 :
509 1283918 : 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 2245969 : do i=1,n-1
515 962051 : cid = species(i)
516 962051 : if (cid <= 0) cycle
517 962051 : Zi = chem_isos% Z(cid)
518 962051 : Ni = chem_isos% N(cid)
519 962051 : isomer_i = chem_isos% isomeric_state(cid)
520 3313399 : do j=i+1,n
521 1067430 : cid = species(j)
522 1067430 : if (cid <= 0) cycle
523 1067430 : Zj = chem_isos% Z(cid)
524 1067430 : Nj = chem_isos% N(cid)
525 1067430 : isomer_j = chem_isos% isomeric_state(cid)
526 1067430 : if (Zj > Zi) cycle
527 166055 : if (Zj == Zi) then
528 166005 : if (Nj > Ni) cycle
529 165915 : if (Nj == Ni) then
530 165885 : if (isomer_j >= isomer_i) cycle
531 : end if
532 : end if
533 : ! exchange i and j
534 80 : species(j) = species(i)
535 80 : species(i) = cid
536 80 : Zi = Zj
537 80 : Ni = Nj
538 2029481 : isomer_i = isomer_j
539 : end do
540 : end do
541 641959 : end subroutine sort
542 :
543 118991 : subroutine long_form
544 : integer :: i, cid
545 : character (len=3) :: op
546 118991 : handle = 'r_'
547 118991 : if (wk_flag) then
548 37865 : op = 'wk_'
549 81126 : else if (ec_flag) then
550 0 : op = 'ec_'
551 : else
552 81126 : op = 'to_'
553 : end if
554 118991 : if (reverse) then
555 1255 : do i = num_in+1,num_in+num_out
556 875 : cid = pspecies(i)
557 1255 : handle = trim(handle) // trim(nuclides% name(cid)) // '_'
558 : end do
559 380 : handle = trim(handle) // op
560 1195 : do i = 1,num_in
561 815 : cid = pspecies(i)
562 815 : handle = trim(handle) // trim(nuclides% name(cid))
563 1195 : if (i < num_in) handle = trim(handle) // '_'
564 : end do
565 : else
566 338066 : do i = 1,num_in
567 219455 : cid = pspecies(i)
568 338066 : handle = trim(handle) // trim(nuclides% name(cid)) // '_'
569 : end do
570 118611 : handle = trim(handle) // op
571 390093 : do i = num_in+1,num_in+num_out
572 271482 : cid = pspecies(i)
573 271482 : handle = trim(handle) // trim(nuclides% name(cid))
574 390093 : if (i < num_in+num_out) handle = trim(handle) // '_'
575 : end do
576 : end if
577 118991 : end subroutine long_form
578 :
579 36948 : logical function one_one()
580 36948 : one_one = .true.
581 36948 : if (in1 == 0 .or. out1 == 0) return
582 36948 : one_one = .false.
583 36948 : if (nuclides% Z(out1) == nuclides% Z(in1) - 1 .and. &
584 : nuclides% N(out1) == nuclides% N(in1) + 1) then
585 10708 : handle = trim(handle) // '_wk'
586 26240 : else if (nuclides% Z(out1) == nuclides% Z(in1) + 1 .and. &
587 : nuclides% N(out1) == nuclides% N(in1) - 1) then
588 26240 : handle = trim(handle) // '_wk-minus'
589 : else
590 : one_one = .true.
591 : end if
592 : end function one_one
593 :
594 106465 : logical function one_two()
595 106465 : one_two = .true.
596 106465 : if (in1 == 0 .or. out1 == 0 .or. out2 == 0 .or. out1 == out2) return
597 106450 : one_two = .false.
598 : if (nuclides% Z(out1) == 0 .and. nuclides% N(out1) == 1 .and. &
599 106450 : nuclides% Z(out2) == nuclides% Z(in1) .and. &
600 : nuclides% N(out2) == nuclides% N(in1) - 1) then
601 36775 : handle = trim(handle) // '_gn'
602 : else if (nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &
603 69675 : nuclides% Z(out2) == nuclides% Z(in1) - 1 .and. &
604 : nuclides% N(out2) == nuclides% N(in1)) then
605 26265 : handle = trim(handle) // '_gp'
606 : else if (nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
607 43410 : nuclides% Z(out2) == nuclides% Z(in1) - 2 .and. &
608 : nuclides% N(out2) == nuclides% N(in1) - 2) then
609 28470 : handle = trim(handle) // '_ga'
610 : else if (nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &
611 14940 : nuclides% Z(out2) == nuclides% Z(in1) - 2 .and. &
612 : nuclides% N(out2) == nuclides% N(in1) + 1) then
613 540 : handle = trim(handle) // '_wk_h1'
614 : else if (nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
615 14400 : nuclides% Z(out2) == nuclides% Z(in1) - 3 .and. &
616 : nuclides% N(out2) == nuclides% N(in1) - 1) then
617 70 : handle = trim(handle) // '_wk_he4'
618 : else
619 : one_two = .true.
620 : end if
621 : end function one_two
622 :
623 182524 : logical function two_one()
624 : include 'formats'
625 182524 : two_one = .true.
626 182524 : if (in1 == 0 .or. in2 == 0 .or. out1 == 0 .or. in1 == in2) return
627 172500 : two_one = .false.
628 : if (nuclides% Z(in1) == 0 .and. nuclides% N(in1) == 1 .and. &
629 172500 : nuclides% Z(out1) == nuclides% Z(in2) .and. &
630 : nuclides% N(out1) == nuclides% N(in2) + 1) then
631 36937 : handle = trim(handle) // '_ng'
632 : else if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
633 135563 : nuclides% Z(out1) == nuclides% Z(in2) + 1 .and. &
634 : nuclides% N(out1) == nuclides% N(in2)) then
635 76972 : handle = trim(handle) // '_pg'
636 : else if (nuclides% Z(in1) == 2 .and. nuclides% N(in1) == 2 .and. &
637 58591 : nuclides% Z(out1) == nuclides% Z(in2) + 2 .and. &
638 : nuclides% N(out1) == nuclides% N(in2) + 2) then
639 58581 : handle = trim(handle) // '_ag'
640 : else
641 : two_one = .true.
642 : end if
643 : end function two_one
644 :
645 252034 : logical function two_two()
646 252034 : two_two = .true.
647 :
648 252034 : 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 : if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
653 : nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
654 : nuclides% Z(out2) == 2 .and. nuclides% N(out2) == 2 .and. &
655 252034 : nuclides% Z(in2) == 3 .and. nuclides% N(in2) == 4) then
656 20 : handle = trim(handle) // '_pa'
657 20 : two_two = .false.
658 20 : return
659 : end if
660 :
661 252014 : if(in1==in2 .or. out1==out2) return
662 :
663 231704 : 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 231704 : nuclides% Z(out2) == nuclides% Z(in2) + 1 .and. &
667 : nuclides% N(out2) == nuclides% N(in2) + 2) then
668 35632 : 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 196072 : nuclides% Z(out2) == nuclides% Z(in2) - 1 .and. &
672 : nuclides% N(out2) == nuclides% N(in2) - 2) then
673 85799 : 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 110273 : nuclides% Z(out2) == nuclides% Z(in2) + 2 .and. &
677 : nuclides% N(out2) == nuclides% N(in2) + 1) then
678 24995 : 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 85278 : nuclides% Z(out2) == nuclides% Z(in2) - 2 .and. &
682 : nuclides% N(out2) == nuclides% N(in2) - 1) then
683 25067 : 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 60211 : nuclides% Z(out2) == nuclides% Z(in2) + 1 .and. &
687 : nuclides% N(out2) == nuclides% N(in2) - 1) then
688 24940 : 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 35271 : nuclides% Z(out2) == nuclides% Z(in2) - 1 .and. &
692 : nuclides% N(out2) == nuclides% N(in2) + 1) then
693 24940 : handle = trim(handle) // '_np'
694 : else
695 : two_two = .true.
696 : end if
697 : end function two_two
698 :
699 577971 : subroutine do_n_to_m(n,m)
700 : integer, intent(in) :: n, m ! each is either 1 or 2
701 577971 : in1 = 0; in2 = 0; out1 = 0; out2 = 0
702 577971 : if (.not. reverse) then
703 411006 : in1 = pspecies(1)
704 411006 : if (n == 2) then
705 346853 : in2 = pspecies(2)
706 346853 : if (in2 == 0) then
707 0 : in1 = 0; return
708 : end if
709 346853 : call switch_if_necessary(in1,in2)
710 346853 : handle = 'r_' // nuclides% name(in2)
711 64153 : else if (n == 1) then
712 64153 : handle = 'r_' // nuclides% name(in1)
713 : else
714 0 : in1 = 0
715 0 : return
716 : end if
717 411006 : out1 = pspecies(n+1)
718 411006 : if (m == 2) then
719 203794 : out2 = pspecies(n+2)
720 203794 : call switch_if_necessary(out1,out2)
721 : end if
722 : else
723 166965 : in1 = pspecies(n+1)
724 166965 : if (m == 2) then
725 87705 : in2 = pspecies(n+2)
726 87705 : if (in2 == 0) then
727 0 : in1 = 0; return
728 : end if
729 87705 : call switch_if_necessary(in1,in2)
730 87705 : handle = 'r_' // nuclides% name(in2)
731 79260 : else if (m == 1) then
732 79260 : handle = 'r_' // nuclides% name(in1)
733 : else
734 0 : in1 = 0
735 0 : return
736 : end if
737 166965 : out1 = pspecies(1)
738 166965 : if (n == 2) then
739 154705 : out2 = pspecies(2)
740 154705 : call switch_if_necessary(out1,out2)
741 : end if
742 : end if
743 : end subroutine do_n_to_m
744 :
745 793057 : subroutine switch_if_necessary(iso1,iso2)
746 : integer, intent(inout) :: iso1, iso2
747 : integer :: j
748 793057 : if (nuclides% Z(iso2) == 1 .and. nuclides% N(iso2) == 0) then ! iso2 is ih1
749 30 : j = iso1; iso1 = iso2; iso2 = j; return
750 : end if
751 793027 : if (nuclides% Z(iso2) == 2 .and. nuclides% N(iso2) == 2) then ! iso2 is ihe4
752 10209 : if (nuclides% Z(iso1) == 1 .and. nuclides% N(iso1) == 0) return ! iso1 is ih1
753 170 : 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
|
