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 net_initialize
21 : use net_def
22 : use net_screen
23 : use chem_def
24 : use math_lib
25 :
26 : use net_approx21, only : num_reactions_func => num_reactions, &
27 : num_mesa_reactions_func => num_mesa_reactions
28 :
29 : implicit none
30 :
31 : integer, parameter :: max_num_special_case_reactants = 5
32 : integer, parameter :: max_num_special_case_reactions = 70
33 : integer :: num_special_case_reactions
34 : integer :: special_case_reactants( &
35 : max_num_special_case_reactants, max_num_special_case_reactions)
36 : character (len=maxlen_reaction_Name) :: &
37 : special_case_reactions(2,max_num_special_case_reactions)
38 :
39 :
40 : contains
41 :
42 9106 : subroutine set_ptrs_for_approx21(n)
43 : use utils_lib, only: fill_with_NaNs, fill_with_NaNs_2D
44 :
45 : type(net_info) :: n
46 :
47 : integer :: num_isos, num_reactions
48 :
49 9106 : num_reactions = num_reactions_func(n% g% add_co56_to_approx21)
50 9106 : if (n% g% add_co56_to_approx21) then
51 : num_isos = 22
52 : else
53 9102 : num_isos = 21
54 : end if
55 :
56 9106 : if(.not.allocated(n% dfdy)) allocate(n% dfdy(num_isos,num_isos))
57 9106 : if(.not.allocated(n% dratdumdy1)) allocate(n% dratdumdy1(num_reactions))
58 9106 : if(.not.allocated(n% dratdumdy2)) allocate(n% dratdumdy2(num_reactions))
59 9106 : if(.not.allocated(n% d_epsnuc_dy)) allocate(n% d_epsnuc_dy(num_isos))
60 9106 : if(.not.allocated(n% d_epsneu_dy)) allocate(n% d_epsneu_dy(num_isos))
61 9106 : if(.not.allocated(n% dydt1)) allocate(n% dydt1(num_isos))
62 9106 : if(.not.allocated(n% dfdt)) allocate(n% dfdT(num_isos))
63 9106 : if(.not.allocated(n% dfdRho)) allocate(n% dfdRho(num_isos))
64 :
65 :
66 9106 : if(n% g% fill_arrays_with_NaNs) then
67 0 : call fill_with_NaNs_2D(n% dfdy)
68 0 : call fill_with_NaNs(n% dratdumdy1)
69 0 : call fill_with_NaNs(n% dratdumdy2)
70 0 : call fill_with_NaNs(n% d_epsnuc_dy)
71 0 : call fill_with_NaNs(n% d_epsneu_dy)
72 0 : call fill_with_NaNs(n% dydt1)
73 0 : call fill_with_NaNs(n% dfdt)
74 0 : call fill_with_NaNs(n% dfdRho)
75 : end if
76 :
77 9106 : end subroutine set_ptrs_for_approx21
78 :
79 :
80 89129 : subroutine setup_net_info(n)
81 : use chem_def
82 : use utils_lib, only: fill_with_NaNs, fill_with_NaNs_2D
83 : type (Net_Info) :: n
84 :
85 : integer :: num_reactions, num_isos, num_wk_reactions
86 :
87 89129 : num_isos = n% g% num_isos
88 89129 : num_wk_reactions = n% g% num_wk_reactions
89 89129 : if (n% g% doing_approx21) then
90 9118 : num_reactions = num_reactions_func(n% g% add_co56_to_approx21)
91 : else
92 80011 : num_reactions = n% g% num_reactions
93 : end if
94 :
95 89129 : if(.not.allocated(n% eps_nuc_categories)) allocate(n% eps_nuc_categories(num_categories))
96 :
97 89129 : if(.not.allocated(n% rate_screened)) allocate(n% rate_screened(num_reactions))
98 89129 : if(.not.allocated(n% rate_screened_dT)) allocate(n% rate_screened_dT(num_reactions))
99 89129 : if(.not.allocated(n% rate_screened_dRho)) allocate(n% rate_screened_dRho(num_reactions))
100 :
101 89129 : if(.not.allocated(n% rate_raw)) allocate(n% rate_raw(num_reactions))
102 89129 : if(.not.allocated(n% rate_raw_dT)) allocate(n% rate_raw_dT(num_reactions))
103 89129 : if(.not.allocated(n% rate_raw_dRho)) allocate(n% rate_raw_dRho(num_reactions))
104 :
105 89129 : if(.not.allocated(n% rate_factors)) allocate(n% rate_factors(num_reactions))
106 :
107 89129 : if(.not.allocated(n% y)) allocate(n% y(num_isos))
108 89129 : if(.not.allocated(n% x)) allocate(n% x(num_isos))
109 :
110 89129 : if(.not.allocated(n% d_eps_nuc_dx)) allocate(n% d_eps_nuc_dx(num_isos))
111 89129 : if(.not.allocated(n% dxdt)) allocate(n% dxdt(num_isos))
112 89129 : if(.not.allocated(n% d_dxdt_dRho)) allocate(n% d_dxdt_dRho(num_isos))
113 89129 : if(.not.allocated(n% d_dxdt_dT)) allocate(n% d_dxdt_dT(num_isos))
114 89129 : if(.not.allocated(n% dydt)) allocate(n% dydt(num_rvs, num_isos))
115 170056 : if(.not.allocated(n% d_dxdt_dx)) allocate(n% d_dxdt_dx(num_isos, num_isos))
116 :
117 89129 : if(.not.allocated(n% d_eps_nuc_dy)) allocate(n% d_eps_nuc_dy(num_isos))
118 170056 : if(.not.allocated(n% d_dydt_dy)) allocate(n% d_dydt_dy(num_isos,num_isos))
119 :
120 89129 : if(.not.allocated(n% lambda)) allocate(n% lambda(num_wk_reactions))
121 89129 : if(.not.allocated(n% dlambda_dlnT)) allocate(n% dlambda_dlnT(num_wk_reactions))
122 89129 : if(.not.allocated(n% dlambda_dlnRho)) allocate(n% dlambda_dlnRho(num_wk_reactions))
123 :
124 89129 : if(.not.allocated(n% Q)) allocate(n% Q(num_wk_reactions))
125 89129 : if(.not.allocated(n% dQ_dlnT)) allocate(n% dQ_dlnT(num_wk_reactions))
126 89129 : if(.not.allocated(n% dQ_dlnRho)) allocate(n% dQ_dlnRho(num_wk_reactions))
127 :
128 89129 : if(.not.allocated(n% Qneu)) allocate(n% Qneu(num_wk_reactions))
129 89129 : if(.not.allocated(n% dQneu_dlnT)) allocate(n% dQneu_dlnT(num_wk_reactions))
130 89129 : if(.not.allocated(n% dQneu_dlnRho)) allocate(n% dQneu_dlnRho(num_wk_reactions))
131 :
132 89129 : if(.not.allocated(n% raw_rate)) allocate(n% raw_rate(num_reactions))
133 89129 : if(.not.allocated(n% screened_rate)) allocate(n% screened_rate(num_reactions))
134 89129 : if(.not.allocated(n% eps_nuc_rate)) allocate(n% eps_nuc_rate(num_reactions))
135 89129 : if(.not.allocated(n% eps_neu_rate)) allocate(n% eps_neu_rate(num_reactions))
136 :
137 89129 : if(n% g% fill_arrays_with_NaNs) then
138 0 : call fill_with_NaNs(n% eps_nuc_categories)
139 0 : call fill_with_NaNs(n% rate_screened)
140 0 : call fill_with_NaNs(n% rate_screened_dt)
141 0 : call fill_with_NaNs(n% rate_screened_drho)
142 0 : call fill_with_NaNs(n% rate_raw)
143 0 : call fill_with_NaNs(n% rate_raw_dt)
144 0 : call fill_with_NaNs(n% rate_raw_drho)
145 0 : call fill_with_NaNs(n% rate_factors)
146 0 : call fill_with_NaNs(n% y)
147 0 : call fill_with_NaNs(n% x)
148 0 : call fill_with_NaNs(n% d_eps_nuc_dx)
149 0 : call fill_with_NaNs(n% dxdt)
150 0 : call fill_with_NaNs(n% d_dxdt_dRho)
151 0 : call fill_with_NaNs(n% d_dxdt_dT)
152 0 : call fill_with_NaNs_2D(n% d_dxdt_dx)
153 0 : call fill_with_NaNs(n% d_eps_nuc_dy)
154 0 : call fill_with_NaNs_2D(n% d_dydt_dy)
155 0 : call fill_with_NaNs(n% lambda)
156 0 : call fill_with_NaNs(n% dlambda_dlnT)
157 0 : call fill_with_NaNs(n% dlambda_dlnRho)
158 0 : call fill_with_NaNs(n% Q)
159 0 : call fill_with_NaNs(n% dQ_dlnT)
160 0 : call fill_with_NaNs(n% dQ_dlnRho)
161 0 : call fill_with_NaNs(n% Qneu)
162 0 : call fill_with_NaNs(n% dQneu_dlnT)
163 0 : call fill_with_NaNs(n% dQneu_dlnRho)
164 0 : call fill_with_NaNs(n% raw_rate)
165 0 : call fill_with_NaNs(n% screened_rate)
166 0 : call fill_with_NaNs(n% eps_nuc_rate)
167 0 : call fill_with_NaNs(n% eps_neu_rate)
168 : end if
169 :
170 :
171 89129 : end subroutine setup_net_info
172 :
173 :
174 19 : subroutine alloc_net_general_info(handle, cache_suffix, ierr)
175 89129 : use rates_def, only: extended_screening
176 : use rates_lib, only: make_rate_tables
177 : use chem_def, only: chem_isos
178 :
179 : integer, intent(in) :: handle
180 : character (len=*), intent(in) :: cache_suffix
181 : integer, intent(out) :: ierr
182 :
183 19 : type (Net_Info) :: n
184 : integer :: ios, status, lwork, num_reactions, &
185 : num_isos, num_wk_reactions, i, iwork
186 : type (Net_General_Info), pointer :: g
187 :
188 : include 'formats'
189 :
190 : ierr = 0
191 :
192 19 : call get_net_ptr(handle, g, ierr)
193 19 : if (ierr /= 0) return
194 :
195 19 : n% g => g
196 :
197 19 : n% reaction_Qs => std_reaction_Qs
198 19 : n% reaction_neuQs => std_reaction_neuQs
199 :
200 19 : num_reactions = g% num_reactions
201 19 : num_isos = g% num_isos
202 19 : num_wk_reactions = g% num_wk_reactions
203 :
204 19 : call setup_net_info(n)
205 :
206 19 : g% cache_suffix = trim(cache_suffix)
207 :
208 : call make_rate_tables( &
209 : g% num_reactions, g% cache_suffix, g% reaction_id, &
210 19 : g% rate_table, g% rattab_f1, nrattab, g% ttab, g% logttab, ierr)
211 19 : if (ierr /= 0) then
212 0 : write(*,*) 'alloc_net_general_info failed in call on make_rate_tables'
213 0 : return
214 : end if
215 :
216 19 : call make_screening_tables(n, ierr)
217 19 : if (ierr /= 0) then
218 0 : write(*,*) 'alloc_net_general_info failed in make_screening_tables'
219 0 : return
220 : end if
221 :
222 :
223 38 : end subroutine alloc_net_general_info
224 :
225 :
226 19 : subroutine set_reaction_max_Z(g)
227 : type (Net_General_Info), pointer :: g
228 : integer :: i, ir, max_Z, max_Z_plus_N
229 : include 'formats'
230 1349 : do i=1, g% num_reactions
231 1311 : ir = g% reaction_id(i)
232 1311 : max_Z = 0
233 1311 : max_Z_plus_N = 0
234 1311 : call update_max_Z(reaction_inputs(2,ir))
235 1311 : call update_max_Z(reaction_inputs(4,ir))
236 1311 : call update_max_Z(reaction_inputs(6,ir))
237 1311 : call update_max_Z(reaction_outputs(2,ir))
238 1311 : call update_max_Z(reaction_outputs(4,ir))
239 1311 : call update_max_Z(reaction_outputs(6,ir))
240 1311 : g% reaction_max_Z(i) = max_Z
241 1330 : g% reaction_max_Z_plus_N_for_max_Z(i) = max_Z_plus_N
242 : !write(*,3) trim(reaction_name(ir)), max_Z, max_Z_plus_N
243 : end do
244 :
245 : contains
246 :
247 7866 : subroutine update_max_Z(iso)
248 : use chem_def, only: chem_isos
249 : integer, intent(in) :: iso
250 : integer :: Z, Z_plus_N
251 7866 : if (iso == 0) return
252 3773 : Z = chem_isos% Z(iso)
253 3773 : if (Z < max_Z) return
254 2036 : Z_plus_N = chem_isos% Z_plus_N(iso)
255 2036 : if (Z > max_Z) then
256 1820 : max_Z = Z
257 1820 : max_Z_plus_N = Z_plus_N
258 216 : else if (Z_plus_N > max_Z_plus_N) then
259 142 : max_Z_plus_N = Z_plus_N
260 : end if
261 7866 : end subroutine update_max_Z
262 :
263 : end subroutine set_reaction_max_Z
264 :
265 :
266 33 : recursive subroutine do_read_net_file(net_filename, handle, ierr)
267 : use utils_def
268 : use utils_lib
269 : use chem_lib, only: chem_get_iso_id
270 : use chem_def, only: chem_isos, ih1, ihe4, ineut
271 : character (len=*), intent(in) :: net_filename
272 : integer, intent(in) :: handle
273 : integer, intent(out) :: ierr
274 :
275 : integer :: iounit, n, i, j, k, t, id, h1, he4, neut
276 : character (len=256) :: buffer, string, filename
277 : logical, parameter :: dbg = .false.
278 : type (Net_General_Info), pointer :: g
279 :
280 33 : ierr = 0
281 33 : call get_net_ptr(handle, g, ierr)
282 33 : if (ierr /= 0) return
283 :
284 33 : if (len_trim(g% net_filename) == 0) g% net_filename = trim(net_filename)
285 :
286 : ! first look in local directory
287 33 : filename = trim(net_filename)
288 33 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
289 33 : if (ierr /= 0) then ! look in local nets directory
290 33 : filename = 'nets/' // trim(net_filename)
291 33 : ierr = 0
292 33 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
293 33 : if (ierr /= 0) then ! look in global nets directory
294 33 : filename = trim(net_dir) // '/nets/' // trim(net_filename)
295 33 : ierr = 0
296 33 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
297 33 : if (ierr /= 0) then
298 0 : write(*,*) 'failed to open net file ' // trim(filename)
299 0 : return
300 : end if
301 : end if
302 : end if
303 :
304 : if (dbg) then
305 : write(*,'(A)')
306 : write(*,*) 'read_net_file <' // trim(filename) // '>'
307 : write(*,*) 'net_filename <' // trim(net_filename) // '>'
308 : write(*,'(A)')
309 : end if
310 :
311 33 : n = 0
312 33 : i = 0
313 :
314 : do
315 99 : t = token(iounit, n, i, buffer, string)
316 66 : select case(t)
317 : case(name_token)
318 17 : select case(string)
319 :
320 : case ('add_isos')
321 17 : call do_isos(.true., ierr)
322 17 : if (ierr /= 0) return
323 :
324 : case ('add_iso')
325 0 : call do_isos(.true., ierr)
326 0 : if (ierr /= 0) return
327 :
328 : case ('remove_isos')
329 0 : call do_isos(.false., ierr)
330 0 : if (ierr /= 0) return
331 :
332 : case ('remove_iso')
333 0 : call do_isos(.false., ierr)
334 0 : if (ierr /= 0) return
335 :
336 : case ('add_reaction')
337 0 : call do_reactions(.true., ierr)
338 0 : if (ierr /= 0) return
339 :
340 : case ('add_reactions')
341 17 : call do_reactions(.true., ierr)
342 17 : if (ierr /= 0) return
343 :
344 : case ('remove_reaction')
345 2 : call do_reactions(.false., ierr)
346 2 : if (ierr /= 0) return
347 :
348 : case ('remove_reactions')
349 6 : call do_reactions(.false., ierr)
350 6 : if (ierr /= 0) return
351 :
352 : case ('add_iso_and_reactions')
353 0 : call do_basic_reactions_for_isos(ierr)
354 0 : if (ierr /= 0) return
355 :
356 : case ('add_isos_and_reactions')
357 0 : call do_basic_reactions_for_isos(ierr)
358 0 : if (ierr /= 0) return
359 :
360 : case ('approx21')
361 6 : call do_approx21(1,ierr)
362 6 : if (ierr /= 0) return
363 :
364 : case ('approx21_plus_co56')
365 4 : call do_approx21(2,ierr)
366 4 : if (ierr /= 0) return
367 :
368 : case ('include')
369 14 : t = token(iounit, n, i, buffer, string)
370 14 : if (t /= string_token) then
371 0 : call error; return
372 : end if
373 14 : call do_read_net_file(string, handle, ierr)
374 14 : if (ierr /= 0) return
375 :
376 : case default
377 80 : call error; return
378 :
379 : end select
380 : case(eof_token)
381 0 : exit
382 : case default
383 99 : call error; return
384 : end select
385 :
386 : end do
387 :
388 33 : close(iounit)
389 :
390 :
391 : ! network veracity checks go here
392 :
393 : ! fxt check that al26 and al26-1 or al26-2 are not both specified
394 33 : i = chem_get_iso_id('al26')
395 : ! if (dbg) write(6,*) i, trim(chem_isos% name(i)), g% net_iso(i)
396 33 : j = chem_get_iso_id('al26-1')
397 : ! if (dbg) write(6,*) j, trim(chem_isos% name(j)), g% net_iso(j)
398 33 : k = chem_get_iso_id('al26-2')
399 : ! if (dbg) write(6,*) k, trim(chem_isos% name(k)), g% net_iso(k)
400 :
401 33 : if ( (g% net_iso(i) == 1 .and. g% net_iso(j) == 1) .or. &
402 : (g% net_iso(i) == 1 .and. g% net_iso(k) == 1)) then
403 0 : string = 'cannot specify al26 and al26-1 or al26-2'
404 0 : call error ; return
405 : end if
406 :
407 : ! fxt check that both al26-1 and al26-2 is specified if one or the other is given
408 33 : if (g% net_iso(j) == 1 .and. g% net_iso(k) == 0) then
409 0 : string = 'must specify al26-2 if al26-1 is set'
410 0 : call error ; return
411 : end if
412 33 : if (g% net_iso(k) == 1 .and. g% net_iso(j) == 0) then
413 0 : string = 'must specify al26-1 if al26-2 is set'
414 0 : call error ; return
415 : end if
416 :
417 : ! done with network veracity checks
418 :
419 66 : if (dbg) then
420 : write(*,'(A)')
421 : write(*,*) 'done read_net_file ' // trim(filename)
422 : write(*,'(A)')
423 : end if
424 :
425 :
426 : contains
427 :
428 :
429 0 : subroutine error
430 : character (len=256) :: message
431 0 : ierr = -1
432 : write(*,'(a)') ' problem reading net file ' // trim(filename) &
433 : // ' error somewhere around here <' // trim(string) // &
434 0 : '>. please check for missing comma or other typo.'
435 0 : close(iounit)
436 33 : end subroutine error
437 :
438 :
439 10 : subroutine do_approx21(which_case, ierr) ! e.g. approx21(cr56)
440 : use chem_lib, only: chem_get_iso_id
441 : use chem_def, only: chem_isos
442 : integer, intent(in) :: which_case
443 : integer, intent(out) :: ierr
444 10 : ierr = 0
445 10 : t = token(iounit, n, i, buffer, string)
446 10 : if (t /= left_paren_token) then
447 0 : call error; return
448 : end if
449 10 : t = token(iounit, n, i, buffer, string)
450 10 : if (t /= name_token) then
451 0 : call error; return
452 : end if
453 10 : id = chem_get_iso_id(string)
454 10 : if (id <= 0) then
455 0 : call error; return
456 : end if
457 10 : t = token(iounit, n, i, buffer, string)
458 10 : if (t /= right_paren_token) then
459 0 : call error; return
460 : end if
461 10 : g% approx21_ye_iso = id
462 : g% fe56ec_n_neut = &
463 : (chem_isos% N(id) - chem_isos% N(ife56)) - &
464 10 : (chem_isos% Z(ife56) - chem_isos% Z(id))
465 10 : g% doing_approx21 = .true.
466 10 : g% add_co56_to_approx21 = (which_case == 2)
467 10 : end subroutine do_approx21
468 :
469 :
470 17 : subroutine do_isos(add_flag, ierr)
471 10 : use chem_lib, only: chem_get_element_id
472 : use chem_lib, only: lookup_ZN
473 : logical, intent(in) :: add_flag
474 : integer, intent(out) :: ierr
475 : logical :: have_next_token
476 : integer :: A, A1, A2, Z
477 17 : ierr = 0
478 17 : t = token(iounit, n, i, buffer, string)
479 17 : if (t /= left_paren_token) then
480 0 : call error; return
481 : end if
482 : have_next_token = .false.
483 : iso_loop: do
484 110 : if (.not. have_next_token) t = token(iounit, n, i, buffer, string)
485 110 : if (t /= name_token) then
486 0 : call error; return
487 : end if
488 110 : id = chem_get_iso_id(string)
489 110 : if (id > 0) then
490 110 : if (add_flag) then
491 : if (dbg) write(*,'(a30,i5,a30)') 'add_net_iso', id, trim(chem_isos% name(id))
492 110 : call add_net_iso(handle, id, ierr)
493 110 : if (ierr /= 0) then
494 0 : write(*,*) 'failed in add_net_iso for ' // trim(chem_isos% name(id))
495 0 : call error; return
496 : end if
497 : else
498 : if (dbg) write(*,'(a30,i5,a30)') 'remove_net_iso', &
499 : id, trim(chem_isos% name(id))
500 0 : call remove_net_iso(handle, id, ierr)
501 0 : if (ierr /= 0) then
502 0 : write(*,*) 'failed in add_net_iso for ' // trim(chem_isos% name(id))
503 0 : call error; return
504 : end if
505 : end if
506 : else
507 0 : Z = chem_get_element_id(string)
508 0 : if (Z < 0) then
509 0 : read(string,fmt=*,iostat=ierr) Z
510 0 : if (Z < 0) then
511 0 : write(*,*) 'unknown name ' // trim(string)
512 0 : call error; return
513 : end if
514 : end if
515 0 : t = token(iounit, n, i, buffer, string)
516 : !write(*,*) 'iso_loop 2nd token ' // trim(string)
517 0 : if (t /= name_token) then
518 : !write(*,*) 'e 2'
519 0 : call error; return
520 : end if
521 0 : read(string,fmt=*,iostat=ierr) A1
522 0 : if (ierr /= 0) then
523 : !write(*,*) 'string <' // trim(string) // '>'
524 0 : call error; return
525 : end if
526 0 : t = token(iounit, n, i, buffer, string)
527 : !write(*,*) 'iso_loop 3rd token ' // trim(string)
528 0 : if (t /= name_token) then
529 : !write(*,*) 'e 3'
530 0 : call error; return
531 : end if
532 0 : read(string,fmt=*,iostat=ierr) A2
533 0 : if (ierr /= 0) then
534 : !write(*,*) 'string <' // trim(string) // '>'
535 0 : call error; return
536 : end if
537 0 : do A = A1, A2
538 0 : id = lookup_ZN(Z, A-Z)
539 0 : if (id <= 0) then
540 0 : write(*,'(a30,5i5,a30)') 'bad iso for ' // trim(string), &
541 0 : Z, A-Z, A, A1, A2
542 0 : ierr = -1; call error; return
543 : stop
544 : end if
545 0 : if (g% net_iso(id) == 0) then
546 0 : if (add_flag) then
547 : if (dbg) write(*,'(a30,i5,a30)') 'add_net_iso', &
548 : id, trim(chem_isos% name(id))
549 0 : call add_net_iso(handle, id, ierr)
550 0 : if (ierr /= 0) then
551 0 : write(*,*) 'failed in add_net_iso for ' // trim(chem_isos% name(id))
552 0 : call error; return
553 : end if
554 : else
555 : if (dbg) write(*,'(a30,i5,a30)') 'remove_net_iso', &
556 : id, trim(chem_isos% name(id))
557 0 : call remove_net_iso(handle, id, ierr)
558 0 : if (ierr /= 0) then
559 0 : write(*,*) 'failed in add_net_iso for ' // trim(chem_isos% name(id))
560 0 : call error; return
561 : end if
562 : end if
563 :
564 : end if
565 : end do
566 : end if
567 110 : t = token(iounit, n, i, buffer, string)
568 110 : if (t == right_paren_token) exit iso_loop
569 110 : if (t /= comma_token) then
570 : have_next_token = .true.
571 : else
572 30 : have_next_token = .false.
573 : end if
574 : end do iso_loop
575 17 : end subroutine do_isos
576 :
577 :
578 25 : subroutine do_reactions(add_flag, ierr)
579 17 : use rates_lib, only: rates_reaction_id
580 : logical, intent(in) :: add_flag
581 : integer, intent(out) :: ierr
582 : integer :: cnt, ir
583 : character (len=16) :: str2
584 : logical :: have_next_token
585 25 : ierr = 0
586 25 : t = token(iounit, n, i, buffer, string)
587 25 : if (t /= left_paren_token) then
588 0 : call error; return
589 : end if
590 : cnt = 0
591 : have_next_token = .false.
592 : reaction_loop: do
593 477 : if (.not. have_next_token) t = token(iounit, n, i, buffer, string)
594 477 : if (t /= name_token) then
595 0 : call error; return
596 : end if
597 477 : id = rates_reaction_id(string)
598 477 : if (id == 0) then
599 18 : if (add_flag) then
600 18 : call add_this_reaction(string, ir, ierr)
601 18 : if (ierr /= 0) then
602 0 : ierr = 0
603 0 : write(*,*) 'add_reaction failed: unknown reaction name ' // trim(string)
604 0 : cnt = cnt+1
605 : end if
606 : else
607 0 : write(*,*) 'remove: unknown reaction name ' // trim(string)
608 0 : cnt = cnt+1
609 : end if
610 : else
611 : if (dbg) write(*,'(a30,i5,a30)') 'reaction ' // trim(string), &
612 : id, trim(reaction_Name(id))
613 459 : if (add_flag) then
614 415 : call add_net_reaction(handle, id, ierr)
615 : else
616 44 : call remove_net_reaction(handle, id, ierr)
617 : end if
618 459 : if (ierr /= 0) then
619 0 : call error; return
620 : end if
621 : end if
622 477 : t = token(iounit, n, i, buffer, string)
623 477 : if (t == right_paren_token) exit reaction_loop
624 477 : if (t /= comma_token) then
625 : have_next_token = .true.
626 : else
627 78 : have_next_token = .false.
628 : end if
629 : end do reaction_loop
630 25 : if (cnt > 0) ierr = -1
631 : end subroutine do_reactions
632 :
633 :
634 0 : subroutine do_basic_reactions_for_isos(ierr)
635 : use chem_lib, only: chem_get_element_id
636 : use chem_lib, only: lookup_ZN
637 : integer, intent(out) :: ierr
638 : integer :: A, A1, A2, Z
639 : logical :: have_next_token, dbg
640 0 : ierr = 0
641 0 : t = token(iounit, n, i, buffer, string)
642 0 : if (t /= left_paren_token) then
643 0 : call error; return
644 : end if
645 : have_next_token = .false.
646 : dbg = .false.
647 : iso_loop: do
648 0 : if (.not. have_next_token) t = token(iounit, n, i, buffer, string)
649 : !write(*,*) 'iso_loop 1st token ' // trim(string)
650 : !dbg = string == 'al26-1'
651 0 : if (t /= name_token) then
652 : !write(*,*) 'e 1'
653 0 : call error; return
654 : end if
655 0 : id = chem_get_iso_id(string)
656 : if (dbg) write(*,*) trim(string) // ' id', id
657 0 : if (id > 0) then
658 0 : Z = chem_isos% Z(id)
659 0 : A1 = chem_isos% N(id) + Z
660 0 : A2 = A1
661 0 : A = A1
662 0 : call do1_iso(id,Z,A,A1,A2,ierr)
663 0 : if (ierr /= 0) then
664 0 : call error; return
665 : end if
666 : else
667 0 : Z = chem_get_element_id(string)
668 0 : if (Z < 0) then
669 0 : read(string,fmt=*,iostat=ierr) Z
670 0 : if (Z < 0) then
671 0 : write(*,*) 'unknown name ' // trim(string)
672 0 : call error; return
673 : end if
674 : end if
675 0 : t = token(iounit, n, i, buffer, string)
676 : !write(*,*) 'iso_loop 2nd token ' // trim(string)
677 0 : if (t /= name_token) then
678 : !write(*,*) 'e 2'
679 0 : call error; return
680 : end if
681 0 : read(string,fmt=*,iostat=ierr) A1
682 0 : if (ierr /= 0) then
683 : !write(*,*) 'string <' // trim(string) // '>'
684 0 : call error; return
685 : end if
686 0 : t = token(iounit, n, i, buffer, string)
687 : !write(*,*) 'iso_loop 3rd token ' // trim(string)
688 0 : if (t /= name_token) then
689 : !write(*,*) 'e 3'
690 0 : call error; return
691 : end if
692 0 : read(string,fmt=*,iostat=ierr) A2
693 0 : if (ierr /= 0) then
694 : !write(*,*) 'string <' // trim(string) // '>'
695 0 : call error; return
696 : end if
697 0 : do A = A1, A2
698 0 : id = lookup_ZN(Z, A-Z)
699 0 : call do1_iso(id,Z,A,A1,A2,ierr)
700 0 : if (ierr /= 0) then
701 0 : call error; return
702 : end if
703 : end do
704 : end if
705 0 : t = token(iounit, n, i, buffer, string)
706 : !write(*,*) 'iso_loop final token ' // trim(string)
707 0 : if (t == right_paren_token) exit iso_loop
708 0 : if (t /= comma_token) then
709 : have_next_token = .true.
710 : else
711 0 : have_next_token = .false.
712 : end if
713 : end do iso_loop
714 :
715 0 : end subroutine do_basic_reactions_for_isos
716 :
717 :
718 0 : subroutine do1_iso(id,Z,A,A1,A2,ierr)
719 : integer, intent(in) :: id,Z,A,A1,A2
720 : integer, intent(out) :: ierr
721 0 : ierr = 0
722 0 : if (id <= 0) then
723 0 : write(*,'(a30,5i5,a30)') 'bad iso for ' // trim(string), Z, A-Z, A, A1, A2
724 0 : ierr = -1; call error; return
725 : stop
726 : end if
727 0 : if (g% net_iso(id) == 0) then
728 : if (dbg) write(*,'(a30,i5,a30)') 'add_net_iso', id, trim(chem_isos% name(id))
729 0 : call add_net_iso(handle, id, ierr)
730 0 : if (ierr /= 0) then
731 0 : write(*,*) 'failed in add_net_iso for ' // trim(chem_isos% name(id))
732 0 : call error; return
733 : end if
734 : end if
735 0 : call do_basic_reactions_for_iso(id, ierr)
736 0 : if (ierr /= 0) then
737 : write(*,*) 'failed in do_basic_reactions_for_iso for ' // &
738 0 : trim(chem_isos% name(id))
739 0 : call error; return
740 : end if
741 0 : end subroutine do1_iso
742 :
743 :
744 0 : subroutine do_basic_reactions_for_iso(id, ierr)
745 : ! pg X, X gp, ag X, X ag, ap X, X pa, ng X, X gn,
746 : ! X wk, X wk_minus, X wk_h1, X wk_he4
747 : integer, intent(in) :: id
748 : integer, intent(out) :: ierr
749 : integer :: Z, N, Z2, N2, other_id
750 : integer :: i, j, k, ir, r_ir
751 : logical :: matches
752 :
753 : include 'formats'
754 :
755 0 : ierr = 0
756 0 : special_loop: do j = 1, num_special_case_reactions
757 0 : matches = .false.
758 0 : do i = 1, max_num_special_case_reactants
759 0 : k = special_case_reactants(i,j)
760 0 : if (k == 0) exit
761 0 : if (k == id) then
762 : matches = .true.; exit
763 : end if
764 : end do
765 0 : if (.not. matches) cycle special_loop
766 : ! check that all are present
767 0 : do i = 1, max_num_special_case_reactants
768 0 : k = special_case_reactants(i,j)
769 0 : if (k /= 0) then
770 0 : if (g% net_iso(k) == 0) cycle special_loop
771 : end if
772 : end do
773 :
774 0 : if (len_trim(special_case_reactions(1,j)) > 0) then
775 0 : call add_this_reaction(special_case_reactions(1,j), ir, ierr)
776 0 : if (ierr /= 0) return
777 : else
778 0 : ir = 0
779 : end if
780 :
781 0 : if (len_trim(special_case_reactions(2,j)) > 0) then
782 0 : call add_this_reaction(special_case_reactions(2,j), r_ir, ierr)
783 0 : if (ierr /= 0) return
784 : else
785 0 : r_ir = 0
786 : end if
787 :
788 0 : if (ir > 0) then
789 0 : reverse_reaction_id(ir) = r_ir
790 : else if (r_ir > 0) then
791 : !write(*,*) 'failed to find reverse for special ' // trim(reaction_name(r_ir))
792 : end if
793 0 : if (r_ir > 0) then
794 0 : reverse_reaction_id(r_ir) = ir
795 : else if (ir > 0) then
796 : !write(*,*) 'failed to find reverse for special ' // trim(reaction_name(ir))
797 : end if
798 :
799 : end do special_loop
800 :
801 0 : Z = chem_isos% Z(id)
802 0 : N = chem_isos% N(id)
803 :
804 0 : if (g% net_iso(ih1) > 0) then
805 0 : other_id = get_iso_id(Z-1, N)
806 0 : call add_basic_reaction(other_id, id, 'pg', 'gp', .true., ir, ierr)
807 0 : if (ierr /= 0) return
808 0 : call add_basic_reaction(id, other_id, 'gp', 'pg', .true., r_ir, ierr)
809 0 : if (ierr /= 0) return
810 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
811 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
812 0 : other_id = get_iso_id(Z+1, N)
813 0 : call add_basic_reaction(id, other_id, 'pg', 'gp', .true., ir, ierr)
814 0 : if (ierr /= 0) return
815 0 : call add_basic_reaction(other_id, id, 'gp', 'pg', .true., r_ir, ierr)
816 0 : if (ierr /= 0) return
817 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
818 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
819 : end if
820 :
821 0 : if (g% net_iso(ihe4) > 0) then
822 0 : if (id == ib8) then
823 0 : call add_this_reaction('r_b8_wk_he4_he4', ir, ierr)
824 0 : if (ierr /= 0) return
825 : end if
826 0 : other_id = get_iso_id(Z-2, N-2)
827 0 : call add_basic_reaction(other_id, id, 'ag', 'ga', .true., ir, ierr)
828 0 : if (ierr /= 0) return
829 0 : call add_basic_reaction(id, other_id, 'ga', 'ag', .true., r_ir, ierr)
830 0 : if (ierr /= 0) return
831 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
832 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
833 0 : other_id = get_iso_id(Z+2, N+2)
834 0 : call add_basic_reaction(id, other_id, 'ag', 'ga', .true., ir, ierr)
835 0 : if (ierr /= 0) return
836 0 : call add_basic_reaction(other_id, id, 'ga', 'ag', .true., r_ir, ierr)
837 0 : if (ierr /= 0) return
838 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
839 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
840 : end if
841 :
842 0 : if (g% net_iso(ih1) > 0 .and. g% net_iso(ihe4) > 0) then
843 0 : if (id /= ihe4) then
844 0 : other_id = get_iso_id(Z-1, N-2)
845 0 : call add_basic_reaction(other_id, id, 'ap', 'pa', .true., ir, ierr)
846 0 : if (ierr /= 0) return
847 0 : call add_basic_reaction(id, other_id, 'pa', 'ap', .true., r_ir, ierr)
848 0 : if (ierr /= 0) return
849 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
850 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
851 0 : other_id = get_iso_id(Z+1, N+2)
852 0 : call add_basic_reaction(id, other_id, 'ap', 'pa', .true., ir, ierr)
853 0 : if (ierr /= 0) return
854 0 : call add_basic_reaction(other_id, id, 'pa', 'ap', .true., r_ir, ierr)
855 0 : if (ierr /= 0) return
856 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
857 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
858 : end if
859 : end if
860 :
861 0 : if (id /= ih1 .and. g% net_iso(ih1) > 0 .and. g% net_iso(ineut) > 0) then
862 0 : other_id = get_iso_id(Z-1, N+1)
863 0 : call add_basic_reaction(id, other_id, 'np', 'pn', .true., ir, ierr)
864 0 : if (ierr /= 0) return
865 0 : call add_basic_reaction(other_id, id, 'pn', 'np', .true., r_ir, ierr)
866 0 : if (ierr /= 0) return
867 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
868 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
869 0 : if (id /= ineut) then
870 0 : other_id = get_iso_id(Z+1, N-1)
871 0 : call add_basic_reaction(other_id, id, 'np', 'pn', .true., ir, ierr)
872 0 : if (ierr /= 0) return
873 0 : call add_basic_reaction(id, other_id, 'pn', 'np', .true., r_ir, ierr)
874 0 : if (ierr /= 0) return
875 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
876 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
877 : end if
878 : end if
879 :
880 0 : if (id /= ihe4 .and. g% net_iso(ihe4) > 0 .and. g% net_iso(ineut) > 0) then
881 0 : other_id = get_iso_id(Z-2, N-1)
882 0 : if (other_id /= ihe4) then
883 0 : call add_basic_reaction(id, other_id, 'na', 'an', .true., ir, ierr)
884 0 : if (ierr /= 0) return
885 0 : call add_basic_reaction(other_id, id, 'an', 'na', .true., r_ir, ierr)
886 0 : if (ierr /= 0) return
887 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
888 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
889 : end if
890 0 : other_id = get_iso_id(Z+2, N+1)
891 0 : if (other_id /= ihe4) then
892 0 : call add_basic_reaction(other_id, id, 'na', 'an', .true., ir, ierr)
893 0 : if (ierr /= 0) return
894 0 : call add_basic_reaction(id, other_id, 'an', 'na', .true., r_ir, ierr)
895 0 : if (ierr /= 0) return
896 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
897 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
898 : end if
899 : end if
900 :
901 0 : if (g% net_iso(ineut) > 0 .and. id /= ih2) then
902 0 : other_id = get_iso_id(Z, N-1)
903 0 : call add_basic_reaction(other_id, id, 'ng', 'gn', .true., ir, ierr)
904 0 : if (ierr /= 0) return
905 0 : call add_basic_reaction(id, other_id, 'gn', 'ng', .true., r_ir, ierr)
906 0 : if (ierr /= 0) return
907 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
908 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
909 : end if
910 :
911 0 : other_id = get_iso_id(Z-1, N+1)
912 0 : call add_basic_reaction(id, other_id, 'wk', '', .true., ir, ierr)
913 0 : if (ierr /= 0) return
914 0 : call add_basic_reaction(other_id, id, 'wk-minus', '', .false., r_ir, ierr)
915 0 : if (ierr /= 0) then
916 0 : write(*,3) 'failed in add_basic_reaction wk-minus', other_id, id
917 0 : return
918 : end if
919 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
920 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
921 :
922 0 : other_id = get_iso_id(Z+1, N-1)
923 0 : call add_basic_reaction(other_id, id, 'wk', '', .true., ir, ierr)
924 0 : if (ierr /= 0) return
925 0 : call add_basic_reaction(id, other_id, 'wk-minus', '', .false., r_ir, ierr)
926 0 : if (ierr /= 0) then
927 0 : write(*,3) 'failed in add_basic_reaction wk-minus', other_id, id
928 0 : return
929 : end if
930 0 : if (ir > 0) reverse_reaction_id(ir) = r_ir
931 0 : if (r_ir > 0) reverse_reaction_id(r_ir) = ir
932 :
933 0 : if (g% net_iso(ih1) > 0) then
934 0 : other_id = get_iso_id(Z-2, N+1)
935 0 : call add_basic_reaction(id, other_id, 'wk-h1', '', .false., ir, ierr)
936 0 : if (ierr /= 0) return
937 0 : other_id = get_iso_id(Z+2, N-1)
938 0 : call add_basic_reaction(other_id, id, 'wk-h1', '', .false., ir, ierr)
939 0 : if (ierr /= 0) return
940 : end if
941 :
942 0 : if (g% net_iso(ihe4) > 0) then
943 0 : other_id = get_iso_id(Z-3, N-1)
944 0 : call add_basic_reaction(id, other_id, 'wk-he4', '', .false., ir, ierr)
945 0 : if (ierr /= 0) return
946 0 : if (ir > 0) reverse_reaction_id(ir) = 0
947 0 : other_id = get_iso_id(Z+3, N+1)
948 0 : call add_basic_reaction(other_id, id, 'wk-he4', '', .false., ir, ierr)
949 0 : if (ierr /= 0) return
950 0 : if (ir > 0) reverse_reaction_id(ir) = 0
951 : end if
952 :
953 :
954 : ! fxt for al26 isomers
955 0 : if (trim(chem_isos% name(id)) == 'al26-1' .and. g% net_iso(id) > 0) then
956 0 : call add_this_reaction('r_al26-1_to_al26-2', ir, ierr)
957 : end if
958 0 : if (trim(chem_isos% name(id)) == 'al26-2' .and. g% net_iso(id) > 0) then
959 0 : call add_this_reaction('r_al26-2_to_al26-1', ir, ierr)
960 : end if
961 :
962 : end subroutine do_basic_reactions_for_iso
963 :
964 :
965 0 : integer function get_iso_id(Z, N)
966 : use chem_lib, only: lookup_ZN
967 : integer, intent(in) :: Z, N
968 0 : get_iso_id = lookup_ZN(Z, N)
969 0 : if (get_iso_id <= 0) return
970 0 : if (g% net_iso(get_iso_id) <= 0) get_iso_id = 0
971 0 : end function get_iso_id
972 :
973 :
974 18 : subroutine add_this_reaction(string, ir, ierr)
975 0 : use rates_lib, only: rates_reaction_id
976 : character (len=*), intent(in) :: string
977 : integer, intent(out) :: ir, ierr
978 : logical :: okay, dbg
979 18 : dbg = .false. ! (string == 'r_h1_li7_to_he4_he4')
980 18 : okay = add_this_reaclib_forward_reaction(string, ir, ierr)
981 18 : if (ierr /= 0) then
982 0 : write(*,*) 'ierr from add_this_reaclib_forward_reaction ' // trim(string)
983 : !stop
984 0 : return
985 : end if
986 18 : if (okay) then
987 : if (dbg) write(*,*) 'add_this_reaclib_forward_reaction ' // &
988 : trim(string), rates_reaction_id(string)
989 : return
990 : end if
991 2 : okay = add_this_reaclib_reverse_reaction(string, ir, ierr)
992 2 : if (ierr /= 0) then
993 0 : write(*,*) 'ierr from add_this_reaclib_reverse_reaction ' // trim(string)
994 : !stop
995 0 : return
996 : end if
997 2 : if (okay) then
998 : if (dbg) write(*,*) 'add_this_reaclib_reverse_reaction ' // trim(string)
999 : return
1000 : end if
1001 0 : okay = add_reaction_for_this_handle(string, ir, ierr)
1002 0 : if (ierr /= 0) then
1003 0 : write(*,*) 'ierr from add_reaction_for_this_handle ' // trim(string)
1004 : !stop
1005 0 : return
1006 : end if
1007 0 : if (.not. okay) then
1008 0 : write(*,*) 'WARNING: failed to add reaction ' // trim(string)
1009 0 : stop
1010 : end if
1011 : if (dbg) write(*,*) 'add_reaction_for_this_handle ' // trim(string)
1012 : if (dbg) call mesa_error(__FILE__,__LINE__,'add_this_reaction')
1013 : end subroutine add_this_reaction
1014 :
1015 :
1016 0 : logical function add_reaction_for_this_handle(string, ir, ierr)
1017 : use rates_lib, only: rates_reaction_id, add_reaction_for_handle
1018 : character (len=*), intent(in) :: string
1019 : integer, intent(out) :: ir, ierr
1020 : include 'formats'
1021 0 : ierr = 0
1022 0 : add_reaction_for_this_handle = .false.
1023 0 : ir = rates_reaction_id(string)
1024 0 : if (ir == 0) then ! check if reaction is defined in reaclib
1025 0 : call add_reaction_for_handle(string, ierr)
1026 0 : if (ierr /= 0) then
1027 0 : write(*,*) 'failed in add_reaction_for_handle for ' // trim(string)
1028 0 : return
1029 : end if
1030 0 : ir = rates_reaction_id(string)
1031 : end if
1032 0 : if (ir > 0) then
1033 0 : call add_net_reaction(handle, ir, ierr)
1034 0 : if (ierr /= 0) then
1035 0 : write(*,*) 'failed in add_reaction_for_this_handle ' // trim(string)
1036 0 : call error
1037 : end if
1038 : add_reaction_for_this_handle = .true.
1039 : end if
1040 : end function add_reaction_for_this_handle
1041 :
1042 :
1043 18 : logical function add_this_reaclib_forward_reaction(string, ir, ierr)
1044 : use rates_lib, only: rates_reaction_id, add_reaction_from_reaclib, reaclib_lookup
1045 : character (len=*), intent(in) :: string
1046 : integer, intent(out) :: ir, ierr
1047 : integer :: indx
1048 : include 'formats'
1049 18 : ierr = 0
1050 18 : add_this_reaclib_forward_reaction = .false.
1051 18 : ir = rates_reaction_id(string)
1052 18 : if (ir == 0) then ! check if reaction is defined in reaclib
1053 18 : indx = reaclib_lookup(string, reaclib_rates% reaction_dict)
1054 18 : if (indx > 0) then ! add a definition for it
1055 16 : call add_reaction_from_reaclib(string, '', indx, ierr)
1056 16 : if (ierr /= 0) then
1057 0 : write(*,*) 'failed in add_this_reaclib_forward_reaction for ' // trim(string)
1058 0 : return
1059 : end if
1060 16 : ir = rates_reaction_id(string)
1061 : end if
1062 : end if
1063 18 : if (ir > 0) then
1064 16 : call add_net_reaction(handle, ir, ierr)
1065 16 : if (ierr /= 0) then
1066 0 : write(*,*) 'failed in add_this_reaclib_forward_reaction ' // trim(string)
1067 0 : call error
1068 : end if
1069 : add_this_reaclib_forward_reaction = .true.
1070 : end if
1071 : end function add_this_reaclib_forward_reaction
1072 :
1073 :
1074 2 : logical function add_this_reaclib_reverse_reaction(string, ir, ierr)
1075 : use rates_lib, only: rates_reaction_id, add_reaction_from_reaclib, reaclib_lookup
1076 : character (len=*), intent(in) :: string
1077 : integer, intent(out) :: ir, ierr
1078 : integer :: indx
1079 : character (len=256) :: forward
1080 2 : ierr = 0
1081 2 : add_this_reaclib_reverse_reaction = .false.
1082 2 : ir = rates_reaction_id(string)
1083 2 : if (ir == 0) then ! check if reaction is defined in reaclib
1084 2 : indx = reaclib_lookup(string, reaclib_rates% reverse_dict)
1085 2 : if (indx > 0) then ! add a definition for it
1086 : call add_reaction_from_reaclib( &
1087 2 : string, reaclib_rates% reaction_handle(indx), indx, ierr)
1088 2 : if (ierr /= 0) return
1089 2 : ir = rates_reaction_id(string)
1090 : end if
1091 : end if
1092 2 : if (ir > 0) then
1093 2 : call add_net_reaction(handle, ir, ierr)
1094 2 : if (ierr /= 0) then
1095 0 : call error
1096 : end if
1097 2 : add_this_reaclib_reverse_reaction = .true.
1098 2 : return
1099 : end if
1100 : end function add_this_reaclib_reverse_reaction
1101 :
1102 :
1103 0 : subroutine add_basic_reaction(iso_in, iso_out, op, reverse_op, warn, ir, ierr)
1104 : use rates_lib, only: rates_reaction_id, add_reaction_from_reaclib, reaclib_lookup
1105 : use rates_def, only: reaction_inputs
1106 : integer, intent(in) :: iso_in, iso_out
1107 : character (len=*), intent(in) :: op, reverse_op
1108 : logical, intent(in) :: warn
1109 : integer, intent(out) :: ir, ierr
1110 : character (len=100) :: string, reverse
1111 : integer :: indx, i
1112 : logical, parameter :: dbg = .false.
1113 : include 'formats'
1114 0 : ierr = 0
1115 0 : ir = 0
1116 :
1117 0 : if (iso_in <= 0 .or. iso_out <= 0) return
1118 :
1119 : string = 'r_' // trim(chem_isos% name(iso_in)) // '_' // &
1120 0 : trim(op) // '_' // trim(chem_isos% name(iso_out))
1121 :
1122 0 : ir = rates_reaction_id(string)
1123 :
1124 0 : if (ir < 0 .or. ir > rates_reaction_id_max) then
1125 0 : write(*,*) 'failed in rates_reaction_id for ' // trim(string)
1126 0 : call mesa_error(__FILE__,__LINE__,'net_init')
1127 : end if
1128 :
1129 0 : if (ir == 0) then ! check if reaction or reverse is defined in reaclib
1130 0 : indx = reaclib_lookup(string, reaclib_rates% reaction_dict)
1131 0 : if (indx > 0) then ! add a definition for it
1132 0 : call add_reaction_from_reaclib(string, '', indx, ierr)
1133 0 : if (ierr /= 0) then
1134 0 : write(*,*) 'failed in add_reaction_from_reaclib for ' // trim(string)
1135 0 : return
1136 : end if
1137 0 : ir = rates_reaction_id(string)
1138 :
1139 0 : if (ir < 0 .or. ir > rates_reaction_id_max) then
1140 0 : write(*,*) 'failed in rates_reaction_id for ' // trim(string)
1141 0 : call mesa_error(__FILE__,__LINE__,'net_init')
1142 : end if
1143 :
1144 : else
1145 0 : ierr = 0
1146 0 : if (len_trim(reverse_op) > 0) then ! check for the reverse reaction in reaclib
1147 : reverse = 'r_' // trim(chem_isos% name(iso_out)) // '_' // &
1148 0 : trim(reverse_op) // '_' // trim(chem_isos% name(iso_in))
1149 0 : indx = reaclib_lookup(reverse, reaclib_rates% reaction_dict)
1150 0 : if (indx > 0) then ! add a definition for it
1151 0 : call add_reaction_from_reaclib(string, reverse, indx, ierr)
1152 0 : if (ierr /= 0) then
1153 : write(*,'(a)') 'failed in add_reaction_from_reaclib for ' &
1154 0 : // trim(string) // ' reverse ' // trim(reverse)
1155 0 : return
1156 : end if
1157 0 : ir = rates_reaction_id(string)
1158 : else
1159 0 : ierr = 0
1160 : end if
1161 : end if
1162 : end if
1163 : end if
1164 :
1165 0 : if (ir > 0) then
1166 0 : call add_net_reaction(handle, ir, ierr)
1167 0 : if (ierr /= 0) then
1168 0 : write(*,*) 'failed in add_net_reaction ' // trim(string)
1169 0 : call error
1170 : end if
1171 : if (dbg) write(*,*) 'add_basic_reaction ' // trim(string)
1172 0 : return
1173 : end if
1174 :
1175 0 : if (warn) then
1176 0 : call integer_dict_lookup(skip_warnings_dict, string, i, ierr)
1177 0 : if (ierr /= 0 .or. i <= 0) then
1178 0 : ierr = 0
1179 : end if
1180 : end if
1181 :
1182 0 : end subroutine add_basic_reaction
1183 :
1184 :
1185 :
1186 : end subroutine do_read_net_file
1187 :
1188 :
1189 19 : subroutine start_net_def(handle, ierr)
1190 : use rates_def, only: rates_reaction_id_max
1191 : integer, intent(in) :: handle
1192 : integer, intent(out) :: ierr
1193 :
1194 : type (Net_General_Info), pointer :: g
1195 :
1196 : ierr = 0
1197 19 : call get_net_ptr(handle, g, ierr)
1198 19 : if (ierr /= 0) return
1199 :
1200 19 : if (.not. associated(g% net_iso)) then
1201 19 : allocate(g% net_iso(num_chem_isos), stat=ierr)
1202 19 : if (ierr /= 0) return
1203 : end if
1204 149283 : g% net_iso(:) = 0
1205 :
1206 19 : if (.not. associated(g% net_reaction)) then
1207 19 : allocate(g% net_reaction(rates_reaction_id_max), stat=ierr)
1208 19 : if (ierr /= 0) return
1209 : end if
1210 6132 : g% net_reaction(:) = 0
1211 :
1212 19 : g% net_has_been_defined = .false.
1213 19 : g% net_filename = ''
1214 :
1215 :
1216 : !call show_scr3('start_net_def')
1217 :
1218 19 : end subroutine start_net_def
1219 :
1220 :
1221 0 : subroutine show_scr3(str)
1222 19 : use rates_def
1223 : use chem_def
1224 : character (len=*), intent(in) :: str
1225 : integer :: ir
1226 : include 'formats'
1227 0 : write(*,*) trim(str)
1228 0 : do ir = 1, rates_reaction_id_max
1229 0 : if (reaction_screening_info(3,ir) <= 0) cycle
1230 : write(*,2) 'scr 3 ' // trim(reaction_Name(ir)) &
1231 : // ' ' // trim(chem_isos% name(reaction_screening_info(1,ir))) &
1232 : // ' ' // trim(chem_isos% name(reaction_screening_info(2,ir))) &
1233 0 : // ' ' // trim(chem_isos% name(reaction_screening_info(3,ir)))
1234 : end do
1235 0 : write(*,'(A)')
1236 0 : end subroutine show_scr3
1237 :
1238 :
1239 19 : subroutine finish_net_def(handle, ierr)
1240 0 : use const_def, only: mev2gr
1241 : use rates_def, only: reaction_names_dict
1242 : use utils_lib, only: integer_dict_create_hash, realloc_integer
1243 : use chem_def, only: chem_isos
1244 : use chem_lib, only: get_mass_excess
1245 : integer, intent(in) :: handle
1246 : integer, intent(out) :: ierr
1247 :
1248 : type (Net_General_Info), pointer :: g
1249 : integer :: i
1250 :
1251 : logical, parameter :: dbg = .false.
1252 :
1253 : include 'formats'
1254 :
1255 : ierr = 0
1256 :
1257 : if (dbg) write(*,*) 'finish_net_def'
1258 :
1259 : ! may have defined some new reactions as part of loading the net
1260 : ! so recreate the dict hash
1261 19 : call integer_dict_create_hash(reaction_names_dict, ierr)
1262 19 : if (ierr /= 0) then
1263 0 : write(*,*) 'FATAL ERROR: finish_net_def failed in integer_dict_create_hash'
1264 0 : return
1265 : end if
1266 :
1267 19 : call get_net_ptr(handle, g, ierr)
1268 19 : if (ierr /= 0) return
1269 :
1270 19 : if (.not. associated(g% net_iso)) then
1271 0 : ierr = -1
1272 0 : return
1273 : end if
1274 :
1275 19 : if (.not. associated(g% net_reaction)) then
1276 0 : ierr = -1
1277 0 : return
1278 : end if
1279 :
1280 19 : if (size(g% net_reaction, dim=1) > rates_reaction_id_max) then
1281 : if (dbg) write(*,*) 'call realloc_integer'
1282 2 : call realloc_integer(g% net_reaction, rates_reaction_id_max, ierr)
1283 2 : if (ierr /= 0) then
1284 0 : write(*,*) 'finish_net_def failed in realloc_integer'
1285 0 : return
1286 : end if
1287 : end if
1288 :
1289 19 : if (g% doing_approx21) then
1290 : if (dbg) write(*,*) 'call mark_approx21'
1291 10 : call do_mark_approx21(handle, ierr)
1292 10 : if (ierr /= 0) return
1293 : end if
1294 :
1295 : if (dbg) write(*,*) 'call setup_iso_info'
1296 19 : call setup_iso_info(g, ierr)
1297 19 : if (ierr /= 0) return
1298 :
1299 : if (dbg) write(*,*) 'call setup_reaction_info'
1300 19 : call setup_reaction_info(g, ierr)
1301 19 : if (ierr /= 0) return
1302 :
1303 : allocate( &
1304 : g% reaction_max_Z(g% num_reactions), &
1305 : g% reaction_max_Z_plus_N_for_max_Z(g% num_reactions), &
1306 19 : stat=ierr)
1307 19 : if (ierr /= 0) return
1308 :
1309 : allocate( &
1310 : g% z158(g% num_isos), &
1311 : g% z52(g% num_isos), &
1312 : g% mion(g% num_isos), &
1313 19 : stat=ierr)
1314 19 : if (ierr /= 0) return
1315 : ! Precompute some powers of Z for screening and coulomb
1316 343 : do i=1, g% num_isos
1317 324 : g% z158(i) = pow(real(chem_isos% Z(g% chem_id(i)),kind=dp),1.58d0)
1318 343 : g% z52(i) = pow(real(chem_isos% Z(g% chem_id(i)),kind=dp),2.50d0) ! 5.d0/2.d0)
1319 : end do
1320 :
1321 343 : do i=1, g% num_isos
1322 343 : g% mion(i) = get_mass_excess(chem_isos, g% chem_id(i)) * mev2gr
1323 : end do
1324 :
1325 19 : call set_reaction_max_Z(g)
1326 :
1327 19 : g% net_has_been_defined = .true.
1328 :
1329 19 : if (g% doing_approx21) then
1330 : if (dbg) write(*,*) 'call set_approx21'
1331 10 : call do_set_approx21(handle, ierr)
1332 10 : if (ierr /= 0) return
1333 : end if
1334 :
1335 95 : if (dbg) write(*,*) 'done finish_net_def'
1336 :
1337 : contains
1338 :
1339 20 : subroutine do_mark_approx21(handle, ierr)
1340 19 : use net_approx21, only: mark_approx21
1341 : integer, intent(in) :: handle
1342 : integer, intent(out) :: ierr
1343 10 : call mark_approx21(handle, ierr)
1344 10 : end subroutine do_mark_approx21
1345 :
1346 :
1347 10 : subroutine do_set_approx21(handle, ierr)
1348 10 : use net_approx21, only: set_approx21
1349 : integer, intent(in) :: handle
1350 : integer, intent(out) :: ierr
1351 10 : call set_approx21(handle, ierr)
1352 10 : end subroutine do_set_approx21
1353 :
1354 :
1355 : end subroutine finish_net_def
1356 :
1357 :
1358 19 : subroutine check_for_hardwired_pairs ! especially for approx21
1359 19 : call check_pair('r_he4_he4_he4_to_c12', 'r_c12_to_he4_he4_he4')
1360 19 : call check_pair('rhe4_breakup', 'rhe4_rebuild')
1361 : ! << should treat this as he4 + n + p -> 3 n + 3 p ???
1362 19 : call check_pair('rprot_to_neut', 'rneut_to_prot')
1363 19 : call check_pair('r_c12_ag_o16', 'r_o16_ga_c12')
1364 19 : call check_pair('rc12ap_to_o16', 'ro16gp_to_c12')
1365 19 : call check_pair('r_o16_ag_ne20', 'r_ne20_ga_o16')
1366 19 : call check_pair('ro16ap_to_ne20', 'rne20gp_to_o16')
1367 19 : call check_pair('r_ne20_ag_mg24', 'r_mg24_ga_ne20')
1368 19 : call check_pair('rne20ap_to_mg24', 'rmg24gp_to_ne20')
1369 19 : call check_pair('r_mg24_ag_si28', 'r_si28_ga_mg24')
1370 19 : call check_pair('rmg24ap_to_si28', 'rsi28gp_to_mg24')
1371 19 : call check_pair('r_si28_ag_s32', 'r_s32_ga_si28')
1372 19 : call check_pair('rsi28ap_to_s32', 'rs32gp_to_si28')
1373 19 : call check_pair('r_s32_ag_ar36', 'r_ar36_ga_s32')
1374 19 : call check_pair('rs32ap_to_ar36', 'rar36gp_to_s32')
1375 19 : call check_pair('r_ar36_ag_ca40', 'r_ca40_ga_ar36')
1376 19 : call check_pair('rar36ap_to_ca40', 'rca40gp_to_ar36')
1377 19 : call check_pair('r_ca40_ag_ti44', 'r_ti44_ga_ca40')
1378 19 : call check_pair('rca40ap_to_ti44', 'rti44gp_to_ca40')
1379 19 : call check_pair('r_ti44_ag_cr48', 'r_cr48_ga_ti44')
1380 19 : call check_pair('rti44ap_to_cr48', 'rcr48gp_to_ti44')
1381 19 : call check_pair('r_cr48_ag_fe52', 'r_fe52_ga_cr48')
1382 19 : call check_pair('rcr48ap_to_fe52', 'rfe52gp_to_cr48')
1383 19 : call check_pair('rfe52aprot_to_fe54', 'rfe54prot_to_fe52')
1384 19 : call check_pair('rfe52neut_to_fe54', 'rfe54g_to_fe52')
1385 19 : call check_pair('rfe54ng_to_fe56', 'rfe56gn_to_fe54')
1386 19 : call check_pair('r_fe52_ag_ni56', 'r_ni56_ga_fe52')
1387 19 : call check_pair('rfe52aprot_to_ni56', 'rni56gprot_to_fe52')
1388 19 : call check_pair('rfe54prot_to_ni56', 'rni56gprot_to_fe54')
1389 19 : call check_pair('rfe54prot_to_ni56', 'rni56gprot_to_fe54')
1390 19 : end subroutine check_for_hardwired_pairs
1391 :
1392 :
1393 570 : subroutine check_pair(s1,s2)
1394 : use rates_def, only: get_rates_reaction_id, reverse_reaction_id
1395 : character (len=*), intent(in) :: s1, s2
1396 : integer :: ir, r_ir
1397 570 : ir = get_rates_reaction_id(s1)
1398 570 : if (ir <= 0) then
1399 0 : write(*,*) 'missing ' // trim(reaction_name(ir))
1400 0 : return
1401 : end if
1402 570 : r_ir = get_rates_reaction_id(s2)
1403 570 : if (r_ir <= 0) then
1404 0 : write(*,*) 'missing ' // trim(reaction_name(r_ir))
1405 0 : return
1406 : end if
1407 570 : reverse_reaction_id(ir) = r_ir
1408 570 : reverse_reaction_id(r_ir) = ir
1409 570 : end subroutine check_pair
1410 :
1411 :
1412 :
1413 19 : subroutine set_reaction_kinds(g)
1414 570 : use chem_def, only: chem_isos
1415 : use rates_lib, only: &
1416 : reaclib_create_handle, reaclib_parse_handle, is_weak_reaction
1417 : type (Net_General_Info), pointer :: g
1418 :
1419 : integer :: iso_in1, iso_in2, iso_out1, iso_out2, &
1420 : num_in1, num_in2, num_out1, num_out2, in_a, in_b, out_c, out_d, &
1421 : num_basic_ng_kind, num_basic_pn_kind, num_basic_pg_kind, &
1422 : num_basic_ap_kind, num_basic_an_kind, num_basic_ag_kind, &
1423 : num_general_one_one_kind, &
1424 : num_general_two_one_kind, num_general_two_two_kind
1425 : integer :: i, ir, num_in, num_out, iso_ids(10), r_ir, r_i, ierr
1426 : character (len=100) :: reverse_handle, op
1427 : character (len=4) :: rstr
1428 19 : integer, pointer :: kind(:), reverse_id(:)
1429 :
1430 : include 'formats'
1431 :
1432 19 : kind => g% reaction_reaclib_kind
1433 19 : reverse_id => g% reverse_id_for_kind_ne_other
1434 :
1435 19 : num_basic_ng_kind = 0
1436 19 : num_basic_pn_kind = 0
1437 19 : num_basic_pg_kind = 0
1438 19 : num_basic_ap_kind = 0
1439 19 : num_basic_an_kind = 0
1440 19 : num_basic_ag_kind = 0
1441 19 : num_general_one_one_kind = 0
1442 19 : num_general_two_one_kind = 0
1443 19 : num_general_two_two_kind = 0
1444 :
1445 1330 : do i=1, g% num_reactions
1446 :
1447 1311 : kind(i) = other_kind
1448 1311 : reverse_id(i) = 0
1449 1311 : ir = g% reaction_id(i)
1450 :
1451 1311 : if (is_weak_reaction(ir)) cycle ! don't do weak reactions
1452 1186 : if (reaction_name(ir)(1:2) /= 'r_') cycle ! don't mess with special ones.
1453 :
1454 1001 : r_ir = reverse_reaction_id(ir)
1455 1001 : if (r_ir <= 0) then
1456 : cycle ! no reverse reaction
1457 : end if
1458 368 : reverse_id(i) = r_ir
1459 368 : r_i = g% net_reaction(r_ir)
1460 368 : if (r_i <= 0) then
1461 : cycle ! reverse reaction not in this net
1462 : end if
1463 :
1464 332 : if (reaction_inputs(6,ir) /= 0) then
1465 : cycle ! more than 2 input species
1466 : end if
1467 332 : if (reaction_outputs(6,ir) /= 0) then
1468 : cycle ! more than 2 output species
1469 : end if
1470 :
1471 332 : num_in1 = reaction_inputs(1,ir)
1472 332 : iso_in1 = reaction_inputs(2,ir)
1473 :
1474 332 : num_in2 = reaction_inputs(3,ir)
1475 332 : iso_in2 = reaction_inputs(4,ir)
1476 :
1477 332 : num_out1 = reaction_outputs(1,ir)
1478 332 : iso_out1 = reaction_outputs(2,ir)
1479 :
1480 332 : num_out2 = reaction_outputs(3,ir)
1481 332 : iso_out2 = reaction_outputs(4,ir)
1482 :
1483 332 : if (iso_in1 == 0 .or. iso_out1 == 0) then
1484 : cycle ! non-standard reaction
1485 : end if
1486 332 : if (num_in1 > 3 .or. num_out1 > 3) then
1487 : cycle ! non-standard reaction
1488 : end if
1489 :
1490 332 : if (iso_in2 == 0) then ! only 1 species on lhs
1491 186 : if (iso_out2 > 0) cycle ! 1 to many is treated as reverse of many to 1
1492 : ! 1 species to 1 species reaction
1493 20 : if (num_in1 == 1 .and. num_out1 == 1) cycle
1494 20 : kind(i) = general_one_one_kind
1495 20 : if (r_ir < ir) then
1496 10 : kind(i) = other_kind
1497 : else
1498 1311 : num_general_one_one_kind = num_general_one_one_kind + 1
1499 : end if
1500 : cycle
1501 : end if
1502 :
1503 146 : if (is_weak_reaction(ir)) cycle
1504 :
1505 146 : if (iso_out2 == 0) then ! 2 to 1
1506 110 : if (num_in1 > 1 .or. num_in2 > 1 .or. num_out1 > 1) then
1507 0 : kind(i) = general_two_one_kind
1508 : if (r_ir <= 0) then
1509 : kind(i) = other_kind
1510 : else
1511 0 : num_general_two_one_kind = num_general_two_one_kind + 1
1512 : end if
1513 : cycle
1514 : end if
1515 : else ! 2 to 2
1516 36 : if (num_in1 > 1 .or. num_in2 > 1 .or. num_out1 > 1 .or. num_out2 > 1) then
1517 0 : kind(i) = general_two_two_kind
1518 : if (r_ir <= 0) then
1519 : kind(i) = other_kind
1520 : else
1521 0 : num_general_two_two_kind = num_general_two_two_kind + 1
1522 : end if
1523 : cycle
1524 : end if
1525 : end if
1526 :
1527 : ! if get here, have a + b on left with different species a and b
1528 : ! and only have 1 of each species on left and right
1529 :
1530 146 : if (chem_isos% Z(iso_in1) <= chem_isos% Z(iso_in2)) then
1531 : in_a = iso_in1
1532 146 : in_b = iso_in2
1533 : else
1534 146 : in_a = iso_in2
1535 146 : in_b = iso_in1
1536 : end if
1537 : ! a + b => c; Z(a) <= Z(b)
1538 :
1539 146 : if (iso_out2 == 0) then ! 2-to-1
1540 :
1541 110 : if (in_a == ineut) then
1542 0 : kind(i) = ng_kind
1543 0 : num_basic_ng_kind = num_basic_ng_kind + 1
1544 110 : else if (in_a == ihe4) then
1545 110 : kind(i) = ag_kind
1546 110 : num_basic_ag_kind = num_basic_ag_kind + 1
1547 0 : else if (in_a == ih1) then
1548 0 : kind(i) = pg_kind
1549 0 : num_basic_pg_kind = num_basic_pg_kind + 1
1550 : else
1551 0 : kind(i) = general_two_one_kind
1552 : if (r_ir <= 0) then
1553 : kind(i) = other_kind
1554 : cycle
1555 : else
1556 0 : num_general_two_one_kind = num_general_two_one_kind + 1
1557 0 : reverse_id(i) = r_ir
1558 : !write(*,*) 'general 2-1 ' // trim(reaction_name(ir))
1559 : end if
1560 : end if
1561 :
1562 : else ! 2-to-2. only do ap,an,pn. skip pa,na,np.
1563 :
1564 36 : if (chem_isos% Z(iso_out1) <= chem_isos% Z(iso_out2)) then
1565 : out_c = iso_out1
1566 36 : out_d = iso_out2
1567 : else
1568 36 : out_c = iso_out2
1569 36 : out_d = iso_out1
1570 : end if
1571 : ! a + b => c + d; Z(a) <= Z(b) and Z(c) <= Z(d)
1572 36 : if (in_a == ihe4) then
1573 28 : if (out_c == ih1) then
1574 28 : kind(i) = ap_kind
1575 28 : num_basic_ap_kind = num_basic_ap_kind + 1
1576 0 : else if (out_c == ineut) then
1577 0 : kind(i) = an_kind
1578 0 : num_basic_an_kind = num_basic_an_kind + 1
1579 : else
1580 0 : kind(i) = general_two_two_kind
1581 : end if
1582 8 : else if (in_a == ih1) then
1583 8 : if (out_c == ineut) then
1584 0 : kind(i) = pn_kind
1585 0 : num_basic_pn_kind = num_basic_pn_kind + 1
1586 8 : else if (out_c == ihe4) then
1587 8 : kind(i) = other_kind
1588 8 : cycle
1589 : else
1590 0 : kind(i) = general_two_two_kind
1591 : end if
1592 0 : else if (in_a == ineut) then
1593 0 : if (out_c == ih1 .or. out_c == ihe4) then
1594 0 : kind(i) = other_kind
1595 0 : cycle
1596 : else
1597 0 : kind(i) = general_two_two_kind
1598 : end if
1599 : else
1600 0 : kind(i) = general_two_two_kind
1601 : end if
1602 :
1603 : end if
1604 :
1605 157 : if (kind(i) == general_two_two_kind) then
1606 0 : if (r_ir > ir) then
1607 0 : kind(i) = other_kind
1608 0 : cycle
1609 : else
1610 1311 : num_general_two_two_kind = num_general_two_two_kind + 1
1611 : end if
1612 : end if
1613 :
1614 : end do
1615 :
1616 : return
1617 :
1618 : if (.true.) then
1619 : do i=1, g% num_reactions
1620 : ir = g% reaction_id(i)
1621 : write(*,3) trim(reaction_name(ir)), i, kind(i)
1622 : !r_ir = reverse_reaction_id(ir)
1623 : !if (kind(i) == ag_kind) write(*,'(a60,i5)') trim(reaction_name(ir)) // &
1624 : ! ' ' // trim(reaction_name(r_ir)), kind(i)
1625 : !if (kind(i) /= other_kind) write(*,'(a60,i5)') trim(reaction_name(ir)) // &
1626 : ! ' ' // trim(reaction_name(r_ir)), kind(i)
1627 : end do
1628 : end if
1629 :
1630 : stop
1631 :
1632 : !return
1633 :
1634 : write(*,2) 'num_basic_ag_kind', num_basic_ag_kind
1635 : write(*,2) 'num_basic_pg_kind', num_basic_pg_kind
1636 : write(*,2) 'num_basic_ng_kind', num_basic_ng_kind
1637 : write(*,2) 'num_basic_pn_kind', num_basic_pn_kind
1638 : write(*,2) 'num_basic_ap_kind', num_basic_ap_kind
1639 : write(*,2) 'num_basic_an_kind', num_basic_an_kind
1640 : write(*,2) 'num_general_one_one_kind', num_general_one_one_kind
1641 : write(*,2) 'num_general_two_one_kind', num_general_two_one_kind
1642 : write(*,2) 'num_general_two_two_kind', num_general_two_two_kind
1643 : write(*,'(A)')
1644 :
1645 : !stop
1646 :
1647 38 : end subroutine set_reaction_kinds
1648 :
1649 :
1650 110 : subroutine add_net_iso(handle, iso_id, ierr)
1651 19 : use chem_def, only: chem_isos
1652 : integer, intent(in) :: handle
1653 : integer, intent(in) :: iso_id
1654 : integer, intent(out) :: ierr
1655 :
1656 : type (Net_General_Info), pointer :: g
1657 :
1658 : ierr = 0
1659 110 : call get_net_ptr(handle, g, ierr)
1660 110 : if (ierr /= 0) return
1661 :
1662 110 : if (.not. associated(g% net_iso)) then
1663 0 : ierr = -1
1664 0 : return
1665 : end if
1666 :
1667 110 : if (g% net_has_been_defined) then
1668 0 : ierr = -1
1669 0 : return
1670 : end if
1671 :
1672 110 : if (iso_id < 0 .or. iso_id > num_chem_isos) then
1673 0 : ierr = -1
1674 0 : return
1675 : end if
1676 :
1677 110 : g% net_iso(iso_id) = 1 ! mark as added
1678 :
1679 110 : end subroutine add_net_iso
1680 :
1681 :
1682 0 : subroutine add_net_isos(handle, num_isos, iso_ids, ierr)
1683 : integer, intent(in) :: handle
1684 : integer, intent(in) :: num_isos, iso_ids(num_isos)
1685 : integer, intent(out) :: ierr
1686 : integer :: i
1687 0 : ierr = 0
1688 0 : do i = 1, num_isos
1689 0 : call add_net_iso(handle, iso_ids(i), ierr)
1690 0 : if (ierr /= 0) return
1691 : end do
1692 110 : end subroutine add_net_isos
1693 :
1694 :
1695 0 : subroutine remove_net_iso(handle, iso_id, ierr)
1696 : integer, intent(in) :: handle
1697 : integer, intent(in) :: iso_id
1698 : integer, intent(out) :: ierr
1699 :
1700 : type (Net_General_Info), pointer :: g
1701 :
1702 : ierr = 0
1703 0 : call get_net_ptr(handle, g, ierr)
1704 0 : if (ierr /= 0) return
1705 :
1706 0 : if (.not. associated(g% net_iso)) then
1707 0 : ierr = -1
1708 0 : return
1709 : end if
1710 :
1711 0 : if (g% net_has_been_defined) then
1712 0 : ierr = -1
1713 0 : return
1714 : end if
1715 :
1716 0 : if (iso_id < 0 .or. iso_id > num_chem_isos) then
1717 0 : ierr = -1
1718 0 : return
1719 : end if
1720 :
1721 0 : g% net_iso(iso_id) = 0 ! mark as removed
1722 :
1723 0 : call remove_reactions_for_iso(g, iso_id, ierr)
1724 :
1725 : end subroutine remove_net_iso
1726 :
1727 :
1728 0 : subroutine remove_net_isos(handle, num_isos, iso_ids, ierr)
1729 : integer, intent(in) :: handle
1730 : integer, intent(in) :: num_isos, iso_ids(:)
1731 : integer, intent(out) :: ierr
1732 : integer :: i
1733 0 : ierr = 0
1734 0 : do i = 1, num_isos
1735 0 : call remove_net_iso(handle, iso_ids(i), ierr)
1736 0 : if (ierr /= 0) return
1737 : end do
1738 : end subroutine remove_net_isos
1739 :
1740 :
1741 433 : subroutine add_net_reaction(handle, reaction_id, ierr)
1742 : use utils_lib, only: realloc_integer
1743 : integer, intent(in) :: handle
1744 : integer, intent(in) :: reaction_id
1745 : integer, intent(out) :: ierr
1746 :
1747 : type (Net_General_Info), pointer :: g
1748 : integer :: old_sz, new_sz
1749 :
1750 : ierr = 0
1751 433 : call get_net_ptr(handle, g, ierr)
1752 433 : if (ierr /= 0) return
1753 :
1754 433 : if (.not. associated(g% net_reaction)) then
1755 0 : ierr = -1
1756 0 : write(*,*) 'must call net_start_def before calling net_add_reaction'
1757 0 : return
1758 : end if
1759 :
1760 433 : if (g% net_has_been_defined) then
1761 0 : ierr = -1
1762 0 : write(*,*) 'must call net_start_def before calling net_add_reaction'
1763 0 : return
1764 : end if
1765 :
1766 433 : if (reaction_id < 0 .or. reaction_id > rates_reaction_id_max) then
1767 0 : ierr = -1
1768 0 : write(*,*) 'invalid reaction_id for net_add_reaction'
1769 0 : return
1770 : end if
1771 :
1772 433 : old_sz = size(g% net_reaction, dim=1)
1773 433 : if (reaction_id > old_sz) then
1774 2 : new_sz = (reaction_id*11)/10 + 100
1775 2 : call realloc_integer(g% net_reaction,new_sz,ierr)
1776 2 : if (ierr /= 0) then
1777 0 : write(*,*) 'add_net_reaction failed in realloc_integer'
1778 0 : return
1779 : end if
1780 266 : g% net_reaction(old_sz+1:new_sz) = 0
1781 : end if
1782 :
1783 433 : g% net_reaction(reaction_id) = 1 ! mark as added
1784 :
1785 : end subroutine add_net_reaction
1786 :
1787 :
1788 0 : subroutine add_net_reactions(handle, num_reactions, reaction_ids, ierr)
1789 : integer, intent(in) :: handle
1790 : integer, intent(in) :: num_reactions, reaction_ids(:)
1791 : integer, intent(out) :: ierr
1792 : integer :: i
1793 0 : ierr = 0
1794 0 : do i = 1, num_reactions
1795 0 : call add_net_reaction(handle, reaction_ids(i), ierr)
1796 0 : if (ierr /= 0) return
1797 : end do
1798 : end subroutine add_net_reactions
1799 :
1800 :
1801 44 : subroutine remove_net_reaction(handle, reaction_id, ierr)
1802 : integer, intent(in) :: handle
1803 : integer, intent(in) :: reaction_id
1804 : integer, intent(out) :: ierr
1805 :
1806 : type (Net_General_Info), pointer :: g
1807 :
1808 : ierr = 0
1809 44 : call get_net_ptr(handle, g, ierr)
1810 44 : if (ierr /= 0) return
1811 :
1812 44 : if (.not. associated(g% net_reaction)) then
1813 0 : ierr = -1
1814 0 : return
1815 : end if
1816 :
1817 44 : if (g% net_has_been_defined) then
1818 0 : ierr = -1
1819 0 : return
1820 : end if
1821 :
1822 44 : if (reaction_id <= 0 .or. reaction_id > rates_reaction_id_max) then
1823 0 : ierr = -1
1824 0 : return
1825 : end if
1826 :
1827 44 : g% net_reaction(reaction_id) = 0 ! mark as removed
1828 :
1829 : end subroutine remove_net_reaction
1830 :
1831 :
1832 0 : subroutine remove_net_reactions(handle, num_reactions, reaction_ids, ierr)
1833 : integer, intent(in) :: handle
1834 : integer, intent(in) :: num_reactions, reaction_ids(:)
1835 : integer, intent(out) :: ierr
1836 : integer :: i
1837 0 : ierr = 0
1838 0 : do i = 1, num_reactions
1839 0 : call remove_net_reaction(handle, reaction_ids(i), ierr)
1840 0 : if (ierr /= 0) return
1841 : end do
1842 : end subroutine remove_net_reactions
1843 :
1844 :
1845 19 : subroutine setup_iso_info(g, ierr)
1846 : type (Net_General_Info), pointer :: g
1847 : integer, intent(out) :: ierr
1848 :
1849 : integer :: i, iso_num, num_isos
1850 19 : integer, pointer :: itab(:), chem_id(:)
1851 :
1852 19 : ierr = 0
1853 19 : itab => g% net_iso
1854 :
1855 149283 : num_isos = sum(itab(:))
1856 19 : if (num_isos <= 0) then
1857 0 : ierr = -1
1858 0 : return
1859 : end if
1860 :
1861 19 : g% num_isos = num_isos
1862 19 : allocate(g% chem_id(num_isos), stat=ierr)
1863 19 : if (ierr /= 0) return
1864 19 : chem_id => g% chem_id
1865 :
1866 19 : iso_num = 0
1867 149283 : do i = 1, num_chem_isos
1868 149264 : if (itab(i) == 0) cycle
1869 324 : iso_num = iso_num + 1
1870 324 : chem_id(iso_num) = i
1871 149283 : itab(i) = iso_num
1872 : end do
1873 :
1874 19 : end subroutine setup_iso_info
1875 :
1876 :
1877 0 : subroutine remove_reactions_for_iso(g, target_iso, ierr)
1878 : type (Net_General_Info), pointer :: g
1879 : integer, intent(in) :: target_iso
1880 : integer, intent(out) :: ierr
1881 :
1882 0 : integer, pointer, dimension(:) :: itab, rtab
1883 : integer :: i, j, ir, &
1884 : num_reaction_inputs, num_reaction_outputs
1885 :
1886 0 : ierr = 0
1887 0 : rtab => g% net_reaction
1888 0 : itab => g% net_iso
1889 :
1890 : !if (target_iso /= 0) &
1891 : ! write(*,*) 'remove_reactions_for_iso ' // trim(chem_isos% name(target_iso))
1892 :
1893 0 : do ir = 1, size(rtab,dim=1)
1894 0 : if (rtab(ir) == 0) cycle ! not used
1895 0 : num_reaction_inputs = get_num_reaction_inputs(ir)
1896 0 : do i = 1, num_reaction_inputs
1897 0 : j = reaction_inputs(2*i,ir)
1898 0 : if (j == target_iso) then
1899 0 : rtab(ir) = 0 ! mark as removed
1900 : !write(*,*) 'remove reaction ' // trim(reaction_Name(ir))
1901 0 : exit
1902 : end if
1903 : end do
1904 0 : num_reaction_outputs = get_num_reaction_outputs(ir)
1905 0 : do i = 1, num_reaction_outputs
1906 0 : j = reaction_outputs(2*i,ir)
1907 0 : if (j == target_iso) then
1908 0 : rtab(ir) = 0 ! mark as removed
1909 : !write(*,*) 'remove reaction ' // trim(reaction_Name(ir))
1910 0 : exit
1911 : end if
1912 : end do
1913 : end do
1914 :
1915 0 : end subroutine remove_reactions_for_iso
1916 :
1917 :
1918 19 : subroutine setup_reaction_info(g, ierr) ! assumes have already called setup_iso_info
1919 : use chem_def, only: chem_isos
1920 : use rates_lib, only: get_weak_rate_id, is_weak_reaction
1921 : use num_lib, only: qsort_string_index
1922 : type (Net_General_Info), pointer :: g
1923 : integer, intent(out) :: ierr
1924 :
1925 : integer :: i, j, k, kind, reaction_num, num_reactions, &
1926 : r1, r2, r3, &
1927 : num_wk_reactions, ir, cid_lhs, cid_rhs, r_ir, r_i, &
1928 : num_reaction_inputs, num_reaction_outputs, cid
1929 : logical :: has_neut, has_prot, found_one
1930 : integer, pointer, dimension(:) :: &
1931 19 : rtab, index, ids, reaction_kind, &
1932 19 : reaction_reaclib_kind, reaction_id, reverse_id_for_kind_ne_other
1933 :
1934 : include 'formats'
1935 :
1936 19 : ierr = 0
1937 19 : rtab => g% net_reaction
1938 :
1939 6150 : num_reactions = sum(rtab(:))
1940 19 : g% num_reactions = num_reactions
1941 : allocate( &
1942 : g% reaction_kind(num_reactions), &
1943 : g% reaction_reaclib_kind(num_reactions), &
1944 : g% reverse_id_for_kind_ne_other(num_reactions), &
1945 : g% reaction_id(num_reactions), &
1946 : g% zs13(num_reactions), &
1947 : g% zhat(num_reactions), &
1948 : g% zhat2(num_reactions), &
1949 : g% lzav(num_reactions), &
1950 : g% aznut(num_reactions), &
1951 : g% zs13inv(num_reactions), &
1952 : g% rate_table(num_reactions, nrattab), &
1953 : g% ttab(nrattab), &
1954 : g% logttab(nrattab), &
1955 : g% rattab_f1(4*nrattab*num_reactions), &
1956 38 : stat=ierr)
1957 19 : if (ierr /= 0) return
1958 19 : reaction_reaclib_kind => g% reaction_reaclib_kind
1959 19 : reaction_id => g% reaction_id
1960 19 : reverse_id_for_kind_ne_other => g% reverse_id_for_kind_ne_other
1961 :
1962 19 : reaction_num = 0
1963 19 : num_wk_reactions = 0
1964 6150 : do ir = 1, rates_reaction_id_max
1965 6131 : if (rtab(ir) == 0) cycle ! reaction not in this net
1966 1311 : reaction_num = reaction_num + 1
1967 1311 : reaction_id(reaction_num) = ir
1968 1311 : if (is_weak_reaction(ir)) &
1969 144 : num_wk_reactions = num_wk_reactions + 1
1970 : end do
1971 :
1972 19 : if (reaction_num /= num_reactions) then
1973 0 : write(*,*) 'reaction_num /= num_reactions', reaction_num, num_reactions
1974 0 : call mesa_error(__FILE__,__LINE__,'setup_reaction_info')
1975 : end if
1976 :
1977 19 : call check_for_hardwired_pairs
1978 :
1979 : ! need to order reactions to ensure bit-for-bit results
1980 : ! so sort reaction_id by reaction_Name
1981 19 : index(1:num_reactions) => g% reaction_kind(1:num_reactions)
1982 19 : ids(1:num_reactions) => g% reverse_id_for_kind_ne_other(1:num_reactions)
1983 1330 : do i=1, num_reactions
1984 1330 : ids(i) = reaction_id(i)
1985 : end do
1986 19 : call qsort_string_index(index,num_reactions,ids,reaction_Name)
1987 : ! use reverse_id_for_kind_ne_other as temp for reordering reaction_id
1988 1330 : do i=1, num_reactions
1989 1330 : reaction_id(i) = ids(index(i))
1990 : !write(*,*) trim(reaction_Name(reaction_id(i)))
1991 : end do
1992 :
1993 19 : call set_reaction_kinds(g)
1994 :
1995 1330 : do i = 1, num_reactions
1996 1311 : ir = reaction_id(i)
1997 1311 : if (ir < 1 .or. ir > rates_reaction_id_max) then
1998 0 : write(*,*) '(ir < 1 .or. ir > rates_reaction_id_max)', &
1999 0 : ir, i, rates_reaction_id_max
2000 0 : call mesa_error(__FILE__,__LINE__,'setup_reaction_info')
2001 : end if
2002 1330 : rtab(ir) = i
2003 : end do
2004 :
2005 19 : i = 1
2006 190 : kind_loop: do kind = 1, max_kind ! reorder by kind of reaction; other_kind goes last
2007 171 : do while (reaction_reaclib_kind(i) == kind)
2008 0 : i = i+1
2009 171 : if (i > num_reactions) exit kind_loop
2010 : end do
2011 19 : do ! move all of the instances of current kind
2012 319 : found_one = .false.
2013 17195 : do j=i+1,num_reactions ! locate the next instance of current kind of reaction
2014 17024 : if (reaction_reaclib_kind(j) /= kind) cycle
2015 :
2016 : ! exchange
2017 :
2018 148 : r1 = reaction_reaclib_kind(j)
2019 148 : r2 = reverse_id_for_kind_ne_other(j)
2020 148 : r3 = reaction_id(j)
2021 :
2022 148 : reaction_reaclib_kind(j) = reaction_reaclib_kind(i)
2023 148 : reverse_id_for_kind_ne_other(j) = reverse_id_for_kind_ne_other(i)
2024 148 : reaction_id(j) = reaction_id(i)
2025 :
2026 148 : reaction_reaclib_kind(i) = r1
2027 148 : reverse_id_for_kind_ne_other(i) = r2
2028 148 : reaction_id(i) = r3
2029 :
2030 148 : rtab(reaction_id(i)) = i
2031 148 : rtab(reaction_id(j)) = j
2032 :
2033 : found_one = .true.
2034 17047 : exit
2035 : end do
2036 :
2037 : if (.not. found_one) exit ! look for another kind
2038 148 : i = i+1 ! next destination
2039 148 : if (i > num_reactions) exit
2040 : end do
2041 : end do kind_loop
2042 :
2043 19 : i = 0
2044 : !g% have_all_reverses = .true.
2045 : do ! reorder so that forward + reverse pairs are adjacent.
2046 1163 : i = i+1
2047 1163 : if (i >= num_reactions) exit
2048 1144 : if (reaction_reaclib_kind(i) == other_kind) cycle
2049 148 : r_ir = reverse_id_for_kind_ne_other(i)
2050 148 : if (r_ir <= 0) then
2051 0 : write(*,*) trim(reaction_name(reaction_id(i)))
2052 0 : write(*,2) 'reaction kind', reaction_reaclib_kind(i)
2053 0 : call mesa_error(__FILE__,__LINE__,'setup_reaction_info: missing reverse id')
2054 : end if
2055 148 : r_i = rtab(r_ir)
2056 148 : if (r_i == 0) then ! reverse reaction not in net
2057 : !g% have_all_reverses = .false.
2058 : cycle
2059 : end if
2060 148 : if (reaction_id(r_i) /= r_ir) call mesa_error(__FILE__,__LINE__,'setup_reaction_info: bad reverse')
2061 148 : ir = reaction_id(i)
2062 148 : if (r_i <= i) then
2063 0 : write(*,2) 'r_i', r_i
2064 0 : write(*,2) 'i', i
2065 0 : write(*,2) 'reverse_id_for_kind_ne_other(r_i)', reverse_id_for_kind_ne_other(r_i)
2066 0 : write(*,2) trim(reaction_Name(ir)), i
2067 0 : write(*,2) trim(reaction_Name(r_ir)), r_i
2068 0 : write(*,*) 'r_i <= i'
2069 0 : stop
2070 : end if
2071 :
2072 148 : i = i+1
2073 :
2074 7278 : do k=r_i-1,i,-1
2075 7130 : reaction_reaclib_kind(k+1) = reaction_reaclib_kind(k)
2076 7130 : reverse_id_for_kind_ne_other(k+1) = reverse_id_for_kind_ne_other(k)
2077 7130 : reaction_id(k+1) = reaction_id(k)
2078 7278 : rtab(reaction_id(k+1)) = k+1
2079 : end do
2080 :
2081 148 : reaction_reaclib_kind(i) = other_kind
2082 148 : reverse_id_for_kind_ne_other(i) = ir
2083 148 : reaction_id(i) = r_ir
2084 1144 : rtab(r_ir) = i
2085 :
2086 : end do
2087 :
2088 19 : g% num_wk_reactions = num_wk_reactions
2089 : allocate( &
2090 : g% weak_reaction_num(num_wk_reactions), &
2091 : g% reaction_id_for_weak_reactions(num_wk_reactions), &
2092 : g% weaklib_ids(num_wk_reactions), &
2093 : g% weak_reaction_index(num_reactions), &
2094 19 : stat=ierr)
2095 19 : if (ierr /= 0) return
2096 :
2097 19 : j = 0
2098 1330 : do i = 1, reaction_num
2099 1311 : ir = reaction_id(i)
2100 1311 : if (.not. is_weak_reaction(ir)) then
2101 1186 : g% weak_reaction_index(i) = 0
2102 1186 : num_reaction_inputs = get_num_reaction_inputs(ir)
2103 1186 : num_reaction_outputs = get_num_reaction_outputs(ir)
2104 1186 : has_neut = .false.
2105 1186 : has_prot = .false.
2106 2946 : do k=1, num_reaction_inputs
2107 1760 : cid = reaction_inputs(2*k,ir)
2108 1760 : if (cid == ineut) has_neut = .true.
2109 2946 : if (cid == iprot .or. cid == ih1) has_prot = .true.
2110 : end do
2111 2850 : do k=1, num_reaction_outputs
2112 1664 : cid = reaction_outputs(2*k,ir)
2113 1664 : if (cid == ineut) has_neut = .true.
2114 2850 : if (cid == iprot .or. cid == ih1) has_prot = .true.
2115 : end do
2116 1186 : if (has_neut .and. .not. has_prot) then
2117 100 : g% reaction_kind(i) = neut_kind
2118 1086 : else if (has_prot) then
2119 607 : g% reaction_kind(i) = prot_kind
2120 : else
2121 479 : g% reaction_kind(i) = other_strong_kind
2122 : end if
2123 : cycle
2124 : end if
2125 125 : g% reaction_kind(i) = weak_kind
2126 125 : j = j+1
2127 125 : g% weak_reaction_index(i) = j
2128 125 : g% weak_reaction_num(j) = i
2129 125 : g% reaction_id_for_weak_reactions(j) = ir
2130 125 : cid_lhs = weak_reaction_info(1,ir)
2131 125 : cid_rhs = weak_reaction_info(2,ir)
2132 : g% weaklib_ids(j) = &
2133 144 : get_weak_rate_id(chem_isos% name(cid_lhs), chem_isos% name(cid_rhs))
2134 : end do
2135 19 : if (j /= num_wk_reactions) then
2136 0 : write(*,3) 'problem with num_wk_reactions in setup_reaction_info', j, num_wk_reactions
2137 0 : call mesa_error(__FILE__,__LINE__,'setup_reaction_info')
2138 : end if
2139 :
2140 38 : end subroutine setup_reaction_info
2141 :
2142 :
2143 : ! Fowler, Caughlan, Zimmerman, Annual Review Astro. Astrophys., 1975.12:69-112. eqn (1).
2144 0 : real(dp) function neutrino_Q(i1, i2)
2145 : integer, intent(in) :: i1, i2 ! i1 decays to i2
2146 0 : real(dp) :: sum, sum2
2147 0 : sum = chem_isos% binding_energy(i2) - chem_isos% binding_energy(i1) - 0.782d0 - 1.022d0
2148 0 : sum = 1.0d0 + sum/0.511d0
2149 0 : sum2 = sum*sum
2150 : neutrino_Q = 0.5d0 * sum * 0.511d0 * (1.0d0 - 1.0d0/sum2) &
2151 0 : * (1.0d0 - 1.0d0/(4.0d0*sum) - 1.0d0/(9.0d0*sum2))
2152 19 : end function neutrino_Q
2153 :
2154 :
2155 19 : subroutine init_special_case_reaction_info
2156 : integer :: i
2157 : include 'formats'
2158 :
2159 19 : i = 0
2160 :
2161 19 : call set(ih1, ih2, 0, 0, 0, 'r_h1_h1_wk_h2', 'r_h1_h1_ec_h2')
2162 :
2163 19 : call set(ineut, ih1, ih2, 0, 0, 'r_neut_h1_h1_to_h1_h2', 'r_h1_h2_to_neut_h1_h1')
2164 :
2165 19 : call set(ih1, ih2, ih3, 0, 0, 'r_h2_h2_to_h1_h3', 'r_h1_h3_to_h2_h2')
2166 :
2167 19 : call set(ih3, ihe3, 0, 0, 0, 'r_he3_ec_h3', '')
2168 :
2169 19 : call set(ineut, ih2, ihe3, 0, 0, 'r_h2_h2_to_neut_he3', 'r_neut_he3_to_h2_h2')
2170 :
2171 19 : call set(ih1, ihe3, ihe4, 0, 0, 'r_h1_he3_wk_he4', '')
2172 :
2173 19 : call set(ih2, ihe4, 0, 0, 0, 'r_h2_h2_to_he4', 'r_he4_to_h2_h2')
2174 :
2175 19 : call set(ineut, ih2, ih3, ihe4, 0, 'r_h2_h3_to_neut_he4', 'r_neut_he4_to_h2_h3')
2176 :
2177 19 : call set(ih1, ih2, ihe3, ihe4, 0, 'r_h2_he3_to_h1_he4', 'r_h1_he4_to_h2_he3')
2178 :
2179 19 : call set(ih2, ih3, ihe3, ihe4, 0, 'r_h3_he3_to_h2_he4', 'r_h2_he4_to_h3_he3')
2180 :
2181 19 : call set(ineut, ih3, ihe4, 0, 0, 'r_h3_h3_to_neut_neut_he4', 'r_neut_neut_he4_to_h3_h3')
2182 :
2183 19 : call set(ih1, ihe3, ihe4, 0, 0, 'r_he3_he3_to_h1_h1_he4', 'r_h1_h1_he4_to_he3_he3')
2184 :
2185 19 : call set(ineut, ih1, ihe4, ili6, 0, 'r_li6_to_neut_h1_he4', 'r_neut_h1_he4_to_li6')
2186 :
2187 : call set(ineut, ih2, ihe4, ili7, 0, 'r_h2_li7_to_neut_he4_he4', &
2188 19 : 'r_neut_he4_he4_to_h2_li7')
2189 :
2190 : call set(ineut, ih3, ihe4, ili7, 0, &
2191 19 : 'r_h3_li7_to_neut_neut_he4_he4', 'r_neut_neut_he4_he4_to_h3_li7')
2192 :
2193 19 : call set(ih1, ih2, ili6, ili7, 0, 'r_h1_li7_to_h2_li6', 'r_h2_li6_to_h1_li7')
2194 :
2195 19 : call set(ibe7, ili7, 0, 0, 0, 'r_be7_wk_li7', '')
2196 :
2197 : ! might also need reverse for r_n13_wk_c13, r_o15_wk_n15
2198 :
2199 19 : call set(ih1, ih2, ihe4, ibe7, 0, 'r_h2_be7_to_h1_he4_he4', 'r_h1_he4_he4_to_h2_be7')
2200 :
2201 19 : call set(ineut, ihe4, ibe7, 0, 0, 'r_neut_be7_to_he4_he4', 'r_he4_he4_to_neut_be7')
2202 :
2203 : call set(ih1, ihe3, ihe4, ibe7, 0, 'r_he3_be7_to_h1_h1_he4_he4', &
2204 19 : 'r_h1_h1_he4_he4_to_he3_be7')
2205 :
2206 19 : call set(ineut, ih2, ili6, ibe7, 0, 'r_neut_be7_to_h2_li6', 'r_h2_li6_to_neut_be7')
2207 :
2208 19 : call set(ineut, ih3, ili7, ibe9, 0, 'r_neut_be9_to_h3_li7', 'r_h3_li7_to_neut_be9')
2209 :
2210 19 : call set(ineut, ihe4, ibe9, 0, 0, 'r_be9_to_neut_he4_he4', 'r_neut_he4_he4_to_be9')
2211 :
2212 19 : call set(ih1, ih2, ihe4, ibe9, 0, 'r_h1_be9_to_h2_he4_he4', 'r_h2_he4_he4_to_h1_be9')
2213 :
2214 19 : call set(ineut, ih1, ihe4, ib8, 0, 'r_neut_b8_to_h1_he4_he4', 'r_h1_he4_he4_to_neut_b8')
2215 :
2216 19 : call set(ih1, ihe4, ib11, 0, 0, 'r_h1_b11_to_he4_he4_he4', 'r_he4_he4_he4_to_h1_b11')
2217 :
2218 19 : call set(ineut, ih3, ibe9, ib11, 0, 'r_neut_b11_to_h3_be9', 'r_h3_be9_to_neut_b11')
2219 :
2220 19 : call set(ih1, ihe4, ic9, 0, 0, 'r_c9_wk_h1_he4_he4', '')
2221 :
2222 : call set(ineut, ihe4, ic11, 0, 0, 'r_neut_c11_to_he4_he4_he4', &
2223 19 : 'r_he4_he4_he4_to_neut_c11')
2224 :
2225 19 : call set(ihe4, ic12, 0, 0, 0, 'r_he4_he4_he4_to_c12', 'r_c12_to_he4_he4_he4')
2226 :
2227 19 : call set(ineut, ih2, ic13, in14, 0, 'r_neut_n14_to_h2_c13', 'r_h2_c13_to_neut_n14')
2228 :
2229 19 : call set(ineut, ih2, ic14, in15, 0, 'r_neut_n15_to_h2_c14', 'r_h2_c14_to_neut_n15')
2230 :
2231 19 : call set(ih1, ihe4, io13, io15, 0, 'r_he4_o13_to_h1_h1_o15', 'r_h1_h1_o15_to_he4_o13')
2232 :
2233 19 : call set(ihe4, ic12, ine20, 0, 0, 'r_c12_c12_to_he4_ne20', 'r_he4_ne20_to_c12_c12')
2234 :
2235 19 : call set(ih1, ic12, ina23, 0, 0, 'r_c12_c12_to_h1_na23', 'r_h1_na23_to_c12_c12')
2236 :
2237 19 : call set(ineut, ic12, img23, 0, 0, 'r_neut_mg23_to_c12_c12', 'r_c12_c12_to_neut_mg23')
2238 :
2239 19 : call set(ihe4, ic12, io16, img24, 0, 'r_c12_o16_to_he4_mg24', 'r_he4_mg24_to_c12_o16')
2240 :
2241 19 : call set(ih1, ic12, io16, ial27, 0, 'r_c12_o16_to_h1_al27', 'r_h1_al27_to_c12_o16')
2242 :
2243 19 : call set(ineut, ic12, io16, isi27, 0, 'r_neut_si27_to_c12_o16', 'r_c12_o16_to_neut_si27')
2244 :
2245 19 : call set(ihe4, io16, isi28, 0, 0, 'r_o16_o16_to_he4_si28', 'r_he4_si28_to_o16_o16')
2246 :
2247 19 : call set(ihe4, ic12, ine20, isi28, 0, 'r_c12_ne20_to_he4_si28', 'r_he4_si28_to_c12_ne20')
2248 :
2249 19 : call set(ih1, io16, ip31, 0, 0, 'r_o16_o16_to_h1_p31', 'r_h1_p31_to_o16_o16')
2250 :
2251 : ! call set(ih2, io16, ip30, 0, 0, 'r_o16_o16_to_h2_p30', 'r_h2_p30_to_o16_o16')
2252 :
2253 19 : call set(ih1, ic12, ine20, ip31, 0, 'r_c12_ne20_to_h1_p31', 'r_h1_p31_to_c12_ne20')
2254 :
2255 19 : call set(ih1, ial25, is27, 0, 0, 'r_s27_wk_h1_h1_al25', '')
2256 :
2257 19 : call set(ineut, io16, is31, 0, 0, 'r_o16_o16_to_neut_s31', 'r_neut_s31_to_o16_o16')
2258 :
2259 19 : call set(ineut, ic12, ine20, is31, 0, 'r_c12_ne20_to_neut_s31', 'r_neut_s31_to_c12_ne20')
2260 :
2261 19 : call set(ih1, ip29, iar31, 0, 0, 'r_ar31_wk_h1_h1_p29', '')
2262 :
2263 : call set(ih1, ihe4, ica36, ica38, 0, 'r_he4_ca36_to_h1_h1_ca38', &
2264 19 : 'r_h1_h1_ca38_to_he4_ca36')
2265 :
2266 19 : call set(ineut, ih1, ih3, ihe3, ihe4, 'r_h3_he3_to_neut_h1_he4', '')
2267 :
2268 19 : call set(ineut, ih1, ihe3, ihe4, ili7, 'r_he3_li7_to_neut_h1_he4_he4', '')
2269 :
2270 19 : call set(ineut, ih1, ih3, ihe4, ibe7, 'r_h3_be7_to_neut_h1_he4_he4', '')
2271 :
2272 19 : num_special_case_reactions = i
2273 :
2274 : !write(*,2) 'num_special_case_reactions', num_special_case_reactions
2275 :
2276 :
2277 : contains
2278 :
2279 969 : subroutine set(i1, i2, i3, i4, i5, s1_in, s2_in)
2280 : integer, intent(in) :: i1, i2, i3, i4, i5
2281 : character (len=*), intent(in) :: s1_in, s2_in
2282 : character (len=maxlen_reaction_Name) :: s1, s2
2283 969 : i = i+1
2284 969 : special_case_reactants(1,i) = i1
2285 969 : special_case_reactants(2,i) = i2
2286 969 : special_case_reactants(3,i) = i3
2287 969 : special_case_reactants(4,i) = i4
2288 969 : special_case_reactants(5,i) = i5
2289 969 : s1 = s1_in
2290 969 : s2 = s2_in
2291 969 : special_case_reactions(1,i) = s1
2292 969 : special_case_reactions(2,i) = s2
2293 969 : end subroutine set
2294 :
2295 : end subroutine init_special_case_reaction_info
2296 :
2297 :
2298 :
2299 : end module net_initialize
2300 :
|
