Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 Josiah Schwab & 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_derivs_support
21 : use net_def, only: Net_General_Info, Net_Info, net_test_partials, &
22 : net_test_partials_val, net_test_partials_dval_dx, net_test_partials_i, &
23 : net_test_partials_iother
24 : use const_def, only: dp, qp, avo, Qconv
25 : use chem_def
26 : use rates_def
27 :
28 : implicit none
29 :
30 :
31 : real(dp), parameter :: r_min = 1d-99
32 :
33 : logical, parameter :: checking_deriv_flags = .false.
34 :
35 :
36 : logical, parameter :: show_rate = .false.
37 : logical, parameter :: show_jac = .false.
38 : logical, parameter :: show_neuQs = .false.
39 :
40 : real(dp), parameter :: show_dydt_y = 7.0097206738032283D-02
41 : logical, parameter :: show_dydt = .false.
42 :
43 : logical, parameter :: show_d_dydt_dRho = .false.
44 : logical, parameter :: show_d_dydt_dT = .false.
45 :
46 : logical, parameter :: show_eps_nuc = .false.
47 : logical, parameter :: show_d_eps_nuc_dy = .false.
48 :
49 : logical, parameter :: checkQs = .false.
50 : real(dp), parameter :: checkQ_frac = 1d-4
51 :
52 :
53 : contains
54 :
55 :
56 0 : real(dp) function isoB(ci)
57 : use chem_lib, only: get_Q
58 : integer, intent(in) :: ci
59 :
60 0 : isoB = get_Q(chem_isos,ci)
61 :
62 0 : end function isoB
63 :
64 :
65 1440140 : subroutine do_in_out( &
66 1440140 : n, dydt, eps_nuc_MeV, i, r_in, &
67 : n_in, i_in, c_in, n_out, i_out, c_out, &
68 : idr, dr, &
69 : deriv_flgs, symbolic, just_dydt)
70 :
71 : type (Net_Info) :: n
72 : real(qp) :: dydt(:,:) ! (num_rvs, num_isos)
73 : real(qp), intent(out) :: eps_nuc_MeV(num_rvs)
74 : integer, intent(in) :: i ! the reaction number
75 : real(dp), intent(in) :: r_in ! coefficient of rate for the reaction
76 : integer, intent(in) :: n_in, n_out ! number of inputs and outputs
77 : integer, dimension(3), intent(in) :: i_in, i_out ! net isotope numbers for the reaction
78 : real(dp), dimension(3), intent(in) :: c_in, c_out ! isotope coefficients in reaction equation
79 : integer, dimension(3), intent(in) :: idr ! isotope number for dr
80 : real(dp), dimension(3), intent(in) :: dr ! coefficient for Jacobian entries d_dydt_dy(idr)
81 : logical, pointer :: deriv_flgs(:)
82 : logical, intent(in) :: symbolic, just_dydt
83 :
84 : ! the purpose of this wrapper is to automatically set the Q values associated with the reaction
85 :
86 : integer :: j
87 1440140 : j = n% g% reaction_id(i)
88 : call do_in_out_neu( &
89 : n, dydt, eps_nuc_MeV, i, r_in, &
90 : n_in, i_in, c_in, n_out, i_out, c_out, &
91 : idr, dr, &
92 : n% reaction_Qs(j), n% reaction_neuQs(j), 0d0, 0d0, &
93 1440140 : deriv_flgs, symbolic, just_dydt)
94 :
95 0 : end subroutine do_in_out
96 :
97 :
98 1440148 : subroutine do_in_out_neu( &
99 1440148 : n, dydt, eps_nuc_MeV, i, r_in, &
100 : n_in, i_in, c_in, n_out, i_out, c_out, &
101 : idr, dr, Q, Qneu, dQneu_dT, dQneu_dRho, &
102 : deriv_flgs, symbolic, just_dydt)
103 :
104 : ! this function handles reactions with 1-3 inputs going to 1-3 outputs
105 :
106 : type (Net_Info) :: n
107 : real(qp) :: dydt(:,:) ! (num_rvs, num_isos)
108 : real(qp), intent(out) :: eps_nuc_MeV(num_rvs)
109 : integer, intent(in) :: i ! the reaction number
110 : real(dp), intent(in) :: r_in ! coefficient of rate for the reaction
111 : integer, intent(in) :: n_in, n_out ! number of inputs and outputs
112 : integer, dimension(3), intent(in) :: i_in, i_out ! net isotope numbers for the reaction
113 : real(dp), dimension(3), intent(in) :: c_in, c_out ! isotope coefficients in reaction equation
114 : integer, dimension(3), intent(in) :: idr ! isotope number for dr
115 : real(dp), dimension(3), intent(in) :: dr ! coefficient for Jacobian entries d_dydt_dy(idr)
116 : real(dp), intent(in) :: Q, Qneu, dQneu_dT, dQneu_dRho
117 : logical, pointer :: deriv_flgs(:)
118 : logical, intent(in) :: symbolic, just_dydt
119 :
120 5760592 : real(dp) :: rvs(num_rvs), d, d1, d2, d3, lhs, rhs, r, checkQ
121 : type (Net_General_Info), pointer :: g
122 1440148 : integer, pointer :: chem_id(:)
123 : integer :: j, cid, icat, reaction_id
124 :
125 : logical :: condition ! for debugging output
126 :
127 : include 'formats'
128 :
129 : ! enforce minimum rate threshold
130 1440148 : r = r_in
131 1440148 : if (r < r_min .or. n% rate_screened(i) < r_min) r = 0
132 :
133 : ! identify reaction category
134 1440148 : icat = reaction_categories(n% g% reaction_id(i))
135 :
136 1440148 : g => n% g
137 1440148 : chem_id => g% chem_id
138 :
139 1440148 : d = n% rate_screened(i)
140 1440148 : d1 = dr(1) * d
141 1440148 : d2 = dr(2) * d
142 1440148 : d3 = dr(3) * d
143 1440148 : rvs(i_rate) = r * n% rate_screened(i)
144 1440148 : rvs(i_rate_dT) = r * n% rate_screened_dT(i)
145 1440148 : rvs(i_rate_dRho) = r * n% rate_screened_dRho(i)
146 :
147 1440148 : n% raw_rate(i) = n% rate_raw(i) * r * avo
148 1440148 : n% screened_rate(i) = n% rate_screened(i) * r * avo
149 1440148 : n% eps_nuc_rate(i) = n% rate_screened(i) * r * (Q - Qneu) * Qconv
150 1440148 : n% eps_neu_rate(i) = n% rate_screened(i) * r * Qneu * Qconv
151 :
152 : ! evaluate left hand side (inputs)
153 1440148 : lhs = 0
154 4320392 : do j = 1, n_in
155 2880244 : call check(i_in(j), 'input', j)
156 2880244 : cid = chem_id(i_in(j))
157 : call do_lhs_iso(n, dydt, i, c_in(j), i_in(j), rvs, idr(1), d1, idr(2), d2, idr(3), d3, &
158 2880244 : symbolic, just_dydt)
159 4320392 : lhs = lhs + c_in(j)*(chem_isos% Z(cid) + chem_isos% N(cid))
160 : end do
161 :
162 : ! evaluate right hand side (outputs)
163 1440148 : rhs = 0
164 3120360 : do j = 1, n_out
165 1680212 : call check(i_out(j), 'output', j)
166 1680212 : cid = chem_id(i_out(j))
167 : call do_rhs_iso(n, dydt, i, c_out(j), i_out(j), rvs, idr(1), d1, idr(2), d2, idr(3), d3, &
168 1680212 : symbolic, just_dydt)
169 3120360 : rhs = rhs + c_out(j)*(chem_isos% Z(cid) + chem_isos% N(cid))
170 : end do
171 :
172 : ! construct eps_nuc and its Rho & T derivatives
173 1440148 : eps_nuc_MeV(i_rate) = eps_nuc_MeV(i_rate) + rvs(i_rate)*(Q-Qneu)
174 1440148 : eps_nuc_MeV(i_rate_dT) = eps_nuc_MeV(i_rate_dT) + rvs(i_rate_dT)*(Q-Qneu) - rvs(i_rate)*dQneu_dT
175 1440148 : eps_nuc_MeV(i_rate_dRho) = eps_nuc_MeV(i_rate_dRho) + rvs(i_rate_dRho)*(Q-Qneu) - rvs(i_rate)*dQneu_dRho
176 :
177 :
178 : ! for debugging
179 1440148 : condition = abs(rvs(1)*Q) > 1d2
180 : if (show_eps_nuc .and. condition) &
181 : write(*,1) trim(reaction_Name(g% reaction_id(i))) // ' eps_nuc', rvs(1)*Q
182 :
183 :
184 : ! set categories and neutrinos
185 1440148 : n% eps_nuc_categories(icat) = n% eps_nuc_categories(icat) + rvs(i_rate)*(Q-Qneu)
186 1440148 : n% eps_neu_total = n% eps_neu_total + rvs(i_rate) * Qneu
187 :
188 : ! for debugging
189 : condition = n% reaction_neuQs(g% reaction_id(i))*rvs(i_rate) /= 0 .and. &
190 1440148 : abs(n% y(g% net_iso(ihe4)) - show_dydt_y) < 1d-20
191 : if (show_neuQs .and. condition) &
192 : write(*,1) trim(reaction_Name(g% reaction_id(i))) // ' neu', &
193 : n% reaction_neuQs(reaction_id)*rvs(:)
194 :
195 :
196 : ! set composition derivatives
197 : !write(*,*) trim(reaction_Name(g% reaction_id(i)))
198 : !write(*,*) idr(1), idr(2), idr(3), lbound(n% d_eps_nuc_dy),ubound(n% d_eps_nuc_dy),d1,d2,d3
199 1440148 : if(idr(1)>0) n% d_eps_nuc_dy(idr(1)) = n% d_eps_nuc_dy(idr(1)) + d1*(Q-Qneu)
200 1440148 : if(idr(2)>0) n% d_eps_nuc_dy(idr(2)) = n% d_eps_nuc_dy(idr(2)) + d2*(Q-Qneu)
201 1440148 : if(idr(3)>0) n% d_eps_nuc_dy(idr(3)) = n% d_eps_nuc_dy(idr(3)) + d3*(Q-Qneu)
202 :
203 : ! for debugging
204 : if(idr(1)>0) then
205 1440148 : condition = chem_id(idr(1)) == ic12
206 : if (show_d_eps_nuc_dy .and. d1 > 0 .and. condition) &
207 : write(*,1) trim(reaction_Name(g% reaction_id(i))) // ' d_epsnuc_dy', &
208 : d, dr(1), d1, Q, d1*Q, n% d_eps_nuc_dy(idr(1)), n% reaction_Qs(reaction_id), &
209 : n% reaction_neuQs(reaction_id)
210 : end if
211 :
212 : if(idr(2)>0) then
213 1440148 : condition = chem_id(idr(2)) == ic12
214 : if (show_d_eps_nuc_dy .and. d2 > 0 .and. condition) &
215 : write(*,1) trim(reaction_Name(g% reaction_id(i))) // ' d_epsnuc_dy', &
216 : d, dr(2), d2, Q, d2*Q, n% d_eps_nuc_dy(idr(2)), n% reaction_Qs(reaction_id), &
217 : n% reaction_neuQs(reaction_id)
218 : end if
219 :
220 : if(idr(3)>0) then
221 1440148 : condition = chem_id(idr(3)) == ic12
222 : if (show_d_eps_nuc_dy .and. d3 > 0 .and. condition) &
223 : write(*,1) trim(reaction_Name(g% reaction_id(i))) // ' d_epsnuc_dy', &
224 : d, dr(3), d3, Q, d3*Q, n% d_eps_nuc_dy(idr(3)), n% reaction_Qs(reaction_id), &
225 : n% reaction_neuQs(reaction_id)
226 : end if
227 :
228 1440148 : call check_balance(n, i, lhs, rhs)
229 : if (checking_deriv_flags) deriv_flgs(i) = .true.
230 :
231 1440148 : if (checkQs) then
232 : checkQ = 0
233 : do j = 1, n_out
234 : cid = chem_id(i_out(j))
235 : checkQ = checkQ + c_out(j)*isoB(cid)
236 : end do
237 : do j = 1, n_in
238 : cid = chem_id(i_in(j))
239 : checkQ = checkQ - c_in(j)*isoB(cid)
240 : end do
241 : if (abs(Q - checkQ) > checkQ_frac*abs(checkQ)) then
242 : write(*,1) 'do_in_out checkQ ' // trim(reaction_Name(g% reaction_id(i))), &
243 : Q, checkQ
244 : !stop
245 : end if
246 : end if
247 :
248 : contains
249 :
250 4560456 : subroutine check(ii, str, jj)
251 : integer, intent(in) :: ii, jj
252 : character (len=*), intent(in) :: str
253 4560456 : if (ii <= 0) then
254 : write(*,*) &
255 0 : 'do_in_out: bad iso num for ' // trim(str), jj, &
256 0 : ' in ' // trim(reaction_Name(g% reaction_id(i)))
257 0 : call mesa_error(__FILE__,__LINE__)
258 : end if
259 4560456 : end subroutine check
260 :
261 : end subroutine do_in_out_neu
262 :
263 :
264 2880244 : subroutine do_lhs_iso( &
265 2880244 : n, dydt, i, c, i1, rvs, i2, dr2, i3, dr3, i4, dr4, symbolic, just_dydt)
266 : type (Net_Info) :: n
267 : real(qp) :: dydt(:,:) ! (num_rvs, num_isos)
268 : integer, intent(in) :: i, i1, i2, i3, i4
269 : real(dp), intent(in) :: c, rvs(:), dr2, dr3, dr4
270 : logical, intent(in) :: symbolic, just_dydt
271 :
272 : ! i1, i2, i3, and 14 are isotope numbers
273 : ! dr2 = dr/dy(i2); dr3 = dr/dy(i3); dr4 = dr/dy(i4)
274 : ! -c * r = dydt(i1)
275 : ! -c * dr2 = d_dydt(i1)_dy(i2)
276 : ! -c * dr3 = d_dydt(i1)_dy(i3)
277 : ! -c * dr4 = d_dydt(i1)_dy(i4)
278 :
279 : integer :: j
280 : integer, pointer :: chem_id(:)
281 2880244 : chem_id => n% g% chem_id
282 :
283 : include 'formats'
284 :
285 2880244 : if (symbolic) then
286 0 : n% d_dydt_dy(i1, i2) = 1
287 0 : if (i3 <= 0) return
288 0 : n% d_dydt_dy(i1, i3) = 1
289 0 : if (i4 <= 0) return
290 0 : n% d_dydt_dy(i1, i4) = 1
291 0 : return
292 : end if
293 :
294 : ! update the dydt terms for i1
295 11520976 : do j=1,num_rvs
296 11520976 : dydt(j,i1) = dydt(j,i1) - c * rvs(j)
297 : end do
298 : if (chem_id(i1) == io16 .and. show_dydt .and. &
299 : abs(n% y(i1) - show_dydt_y) < 1d-20) &
300 : write(*,1) 'lhs ' // trim(reaction_Name(n% g% reaction_id(i))), &
301 : -c * rvs(i_rate), dydt(i_rate,i1), &
302 : n% rate_screened(i), &
303 : n% rate_raw(i), &
304 : n% y(i1)
305 : if (chem_id(i1) == ihe4 .and. show_d_dydt_dRho .and. &
306 : abs(n% y(i1) - show_dydt_y) < 1d-20) &
307 : write(*,1) trim(reaction_Name(n% g% reaction_id(i))) // ' d_dydt_dRho ' &
308 : // trim(chem_isos% name(chem_id(i1))), &
309 : -c * rvs(i_rate_dRho), n% rate_screened(i), &
310 : n% rate_screened_dRho(i), n% y(i1)
311 : if (chem_id(i1) == ihe4 .and. show_d_dydt_dT .and. &
312 : abs(n% y(i1) - show_dydt_y) < 1d-20) &
313 : write(*,1) trim(reaction_Name(n% g% reaction_id(i))) // ' d_dydt_dT ' &
314 : // trim(chem_isos% name(chem_id(i1))), &
315 : -c * rvs(i_rate_dT), n% rate_screened(i), &
316 : n% rate_screened_dT(i), n% y(i1)
317 :
318 2880244 : if (just_dydt) return
319 :
320 : ! update the Jacobian for d_dydt(i1)_dy(i2)
321 2880244 : n% d_dydt_dy(i1, i2) = n% d_dydt_dy(i1, i2) - c * dr2
322 :
323 : if (chem_id(i1) == ini56 .and. show_jac .and. c * dr2 /= 0) &
324 : write(*,1) trim(reaction_Name(n% g% reaction_id(i))) // ' d_dydt_dy l2 dr2 ' &
325 : // trim(chem_isos% name(chem_id(i1))) // ' ' // trim(chem_isos% name(chem_id(i2))), &
326 : - c * dr2, n% d_dydt_dy(i1, i2)
327 :
328 2880244 : if (i3 <= 0) return
329 :
330 : ! update the Jacobian for d_dydt(i1)_dy(i3)
331 2560188 : n% d_dydt_dy(i1, i3) = n% d_dydt_dy(i1, i3) - c * dr3
332 :
333 : if (chem_id(i1) == ini56 .and. show_jac .and. c * dr3 /= 0) &
334 : write(*,1) trim(reaction_Name(n% g% reaction_id(i))) // ' d_dydt_dy l3 dr3 ' &
335 : // trim(chem_isos% name(chem_id(i1))) // ' ' // trim(chem_isos% name(chem_id(i3))), &
336 : - c * dr3, n% d_dydt_dy(i1, i3)
337 :
338 2560188 : if (i4 <= 0) return
339 :
340 : !write(*,4) trim(reaction_Name(n% g% reaction_id(i))) // ' do_lhs_iso', i, i1, i4
341 : ! update the Jacobian for d_dydt(i1)_dy(i4)
342 960012 : n% d_dydt_dy(i1, i4) = n% d_dydt_dy(i1, i4) - c * dr4
343 :
344 : if (chem_id(i1) == ini56 .and. show_jac .and. c * dr4 /= 0) &
345 : write(*,1) trim(reaction_Name(n% g% reaction_id(i))) // ' d_dydt_dy l4 dr4 ' &
346 : // trim(chem_isos% name(chem_id(i1))) // ' ' // &
347 : trim(chem_isos% name(chem_id(i4))), &
348 : - c * dr4, n% d_dydt_dy(i1, i4)
349 :
350 : end subroutine do_lhs_iso
351 :
352 :
353 1680212 : subroutine do_rhs_iso( &
354 1680212 : n, dydt, i, c, i1, rvs, i2, dr2, i3, dr3, i4, dr4, symbolic, just_dydt)
355 : type (Net_Info) :: n
356 : real(qp) :: dydt(:,:) ! (num_rvs, num_isos)
357 : integer, intent(in) :: i, i1, i2, i3, i4
358 : real(dp), intent(in) :: c, rvs(:), dr2, dr3, dr4
359 : logical, intent(in) :: symbolic, just_dydt
360 :
361 : ! i1, i2, i3, and 14 are isotope numbers
362 : ! dr2 = dr/dy(i2); dr3 = dr/dy(i3); dr4 = dr/dy(i4)
363 : ! c * r = dydt(i1)
364 : ! c * dr2 = d_dydt(i1)_dy(i2)
365 : ! c * dr3 = d_dydt(i1)_dy(i3)
366 : ! c * dr4 = d_dydt(i1)_dy(i4)
367 :
368 : integer :: j
369 : integer, pointer :: chem_id(:)
370 1680212 : chem_id => n% g% chem_id
371 :
372 : include 'formats'
373 :
374 1680212 : if (symbolic) then
375 0 : n% d_dydt_dy(i1, i2) = 1
376 0 : if (i3 <= 0) return
377 0 : n% d_dydt_dy(i1, i3) = 1
378 0 : if (i4 <= 0) return
379 0 : n% d_dydt_dy(i1, i4) = 1
380 0 : return
381 : end if
382 :
383 : ! update the dydt terms for i1
384 6720848 : do j=1,num_rvs
385 6720848 : dydt(j,i1) = dydt(j,i1) + c * rvs(j)
386 : end do
387 : if (chem_id(i1) == io16 .and. show_dydt .and. &
388 : abs(n% y(i1) - show_dydt_y) < 1d-20) &
389 : write(*,1) 'rhs ' // trim(reaction_Name(n% g% reaction_id(i))), &
390 : c * rvs(i_rate), dydt(i_rate,i1), &
391 : n% rate_screened(i), n% rate_raw(i), n% y(i1)
392 : if (chem_id(i1) == ihe4 .and. show_d_dydt_dRho .and. &
393 : abs(n% y(i1) - show_dydt_y) < 1d-20) &
394 : write(*,1) trim(reaction_Name(n% g% reaction_id(i))) // ' d_dydt_dRho ' &
395 : // trim(chem_isos% name(chem_id(i1))), &
396 : c * rvs(i_rate_dRho), n% rate_screened(i), &
397 : n% rate_screened_dRho(i), n% y(i1)
398 : if (chem_id(i1) == ihe4 .and. show_d_dydt_dT .and. &
399 : abs(n% y(i1) - show_dydt_y) < 1d-20) &
400 : write(*,1) trim(reaction_Name(n% g% reaction_id(i))) // ' d_dydt_dT ' &
401 : // trim(chem_isos% name(chem_id(i1))), &
402 : c * rvs(i_rate_dT), n% rate_screened(i), &
403 : n% rate_screened_dT(i), n% y(i1)
404 :
405 1680212 : if (just_dydt) return
406 :
407 : ! update the Jacobian for d_dydt(i1)_dy(i2)
408 1680212 : n% d_dydt_dy(i1, i2) = n% d_dydt_dy(i1, i2) + c * dr2
409 :
410 : if (chem_id(i1) == ini56 .and. show_jac .and. c * dr2 /= 0) &
411 : write(*,1) trim(reaction_Name(n% g% reaction_id(i))) // ' d_dydt_dy r2 dr2 ' &
412 : // trim(chem_isos% name(chem_id(i1))) // ' ' // &
413 : trim(chem_isos% name(chem_id(i2))), &
414 : c * dr2, n% d_dydt_dy(i1, i2)
415 :
416 1680212 : if (i3 <= 0) return
417 :
418 : ! update the Jacobian for d_dydt(i1)_dy(i3)
419 1280126 : n% d_dydt_dy(i1, i3) = n% d_dydt_dy(i1, i3) + c * dr3
420 :
421 : if (chem_id(i1) == ini56 .and. show_jac .and. c * dr3 /= 0) &
422 : write(*,1) trim(reaction_Name(n% g% reaction_id(i))) // ' d_dydt_dy r3 dr3 ' &
423 : // trim(chem_isos% name(chem_id(i1))) // ' ' // trim(chem_isos% name(chem_id(i3))), &
424 : c * dr3, n% d_dydt_dy(i1, i3)
425 :
426 1280126 : if (i4 <= 0) return
427 :
428 : ! update the Jacobian for d_dydt(i1)_dy(i4)
429 400004 : n% d_dydt_dy(i1, i4) = n% d_dydt_dy(i1, i4) + c * dr4
430 :
431 : if (chem_id(i1) == ini56 .and. show_jac .and. c * dr4 /= 0) &
432 : write(*,1) trim(reaction_Name(n% g% reaction_id(i))) // ' d_dydt_dy r4 dr4 ' &
433 : // trim(chem_isos% name(chem_id(i1))) // ' ' // trim(chem_isos% name(chem_id(i4))), &
434 : c * dr4, n% d_dydt_dy(i1, i4)
435 :
436 : end subroutine do_rhs_iso
437 :
438 :
439 1440148 : subroutine check_balance(n, i, lhs, rhs) ! check conservation of nucleons
440 : type (Net_Info) :: n
441 : integer, intent(in) :: i
442 : real(dp), intent(in) :: lhs, rhs
443 1440148 : if (lhs == rhs) return
444 0 : if (abs(lhs-rhs) < 1d-6) return
445 0 : write(*,'(2a)') 'non-conservation of nucleons in reaction ', &
446 0 : reaction_Name(n% g% reaction_id(i))
447 0 : write(*,*) 'lhs aion sum', lhs
448 0 : write(*,*) 'rhs aion sum', rhs
449 0 : stop
450 : end subroutine check_balance
451 :
452 :
453 :
454 : end module net_derivs_support
|
