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
21 : use net_def
22 : use const_def, only: dp
23 : use chem_def
24 : use net_derivs_support
25 : use rates_def
26 :
27 : implicit none
28 :
29 : real(dp), parameter :: tiny_rate = 1d-50
30 :
31 : contains
32 :
33 80002 : subroutine get_derivs( &
34 80002 : n, dydt, eps_nuc_MeV, eta, ye, logtemp, temp, den, abar, zbar, &
35 80002 : num_reactions, rate_factors, &
36 : symbolic, just_dydt, ierr)
37 : type (Net_Info) :: n
38 : real(qp), intent(inout) :: dydt(:,:)
39 : real(qp), intent(out) :: eps_nuc_MeV(:)
40 : integer, intent(in) :: num_reactions
41 : real(dp), intent(in) ::eta, ye, logtemp, temp, den, abar, zbar, &
42 : rate_factors(:)
43 : logical, intent(in) :: symbolic, just_dydt
44 : integer, intent(out) :: ierr
45 :
46 : logical :: all_okay, show_derivs_dydt
47 : integer, pointer, dimension(:) :: &
48 80002 : reaction_kind, reverse_id, itab, rtab, reaction_id
49 : integer :: ir, r_ir, r_i, i, j, op_err, kind, jmax, icat_f, icat_r
50 160004 : logical, target :: deriv_flgs_data(num_reactions)
51 80002 : logical, pointer :: deriv_flgs(:)
52 : type (Net_General_Info), pointer :: g
53 80002 : real(dp) :: T9, T932, eps_nuc_cancel_factor, eps_factor, &
54 80002 : old_eps_nuc_categories_val
55 :
56 : include 'formats'
57 :
58 80002 : ierr = 0
59 :
60 80002 : T9 = temp*1d-9
61 80002 : T932 = T9*sqrt(T9)
62 :
63 80002 : g => n% g
64 :
65 : if (.true. .or. logtemp <= g% logT_lo_eps_nuc_cancel) then
66 80002 : eps_nuc_cancel_factor = 1d0
67 : else if (logtemp >= g% logT_hi_eps_nuc_cancel) then
68 : eps_nuc_cancel_factor = 0d0
69 : else
70 : eps_nuc_cancel_factor = &
71 : (g% logT_hi_eps_nuc_cancel - logtemp)/&
72 : (g% logT_hi_eps_nuc_cancel - g% logT_lo_eps_nuc_cancel)
73 : end if
74 :
75 : !write(*,1) 'eps_nuc_cancel_factor', eps_nuc_cancel_factor, logtemp
76 :
77 80002 : itab => g% net_iso
78 80002 : rtab => g% net_reaction
79 80002 : reaction_kind => g% reaction_reaclib_kind
80 80002 : reverse_id => g% reverse_id_for_kind_ne_other
81 80002 : reaction_id => g% reaction_id
82 :
83 80002 : show_derivs_dydt = .false.
84 : if (show_dydt) then
85 : i = itab(io16)
86 : if (i > 0) &
87 : show_derivs_dydt = (abs(n% y(i) - show_dydt_y) < 1d-14)
88 : end if
89 :
90 :
91 80002 : deriv_flgs => deriv_flgs_data
92 : if (checking_deriv_flags) deriv_flgs(:) = .false.
93 :
94 2640218 : dydt = 0
95 5841170 : if (.not. just_dydt) n% d_dydt_dy = 0
96 :
97 : ierr = 0
98 320008 : eps_nuc_MeV = 0d0
99 :
100 : if (just_dydt) then
101 : jmax = 1 ! = i_rate
102 : else
103 : jmax = num_rvs
104 : end if
105 :
106 : ! Update special rates that depend on the composition
107 2560162 : do i=1,num_reactions
108 : call update_special_rates(n, dydt, eps_nuc_MeV, i, eta, ye, temp, den, abar, zbar, &
109 : num_reactions, rate_factors, rtab, itab, &
110 : deriv_flgs, symbolic, just_dydt, &
111 2560162 : ierr)
112 : end do
113 :
114 80002 : old_eps_nuc_categories_val = 0
115 :
116 80002 : i = 1
117 2560162 : do while (i <= num_reactions)
118 :
119 2480160 : if (ierr /= 0) exit
120 :
121 2480160 : ir = reaction_id(i)
122 2480160 : icat_f = reaction_categories(ir)
123 :
124 2480160 : if (i < num_reactions) then
125 2400158 : r_i = i+1 ! reactions are ordered so reverse immediately follows forward
126 2400158 : r_ir = reaction_id(r_i)
127 2400158 : icat_r = reaction_categories(r_ir)
128 : else
129 80002 : r_i = 0
130 80002 : r_ir = 0
131 80002 : icat_r = 0
132 : end if
133 :
134 :
135 2480160 : kind = reaction_kind(i)
136 2480160 : if ( &
137 : !kind == ng_kind .or. &
138 : !kind == pn_kind .or. &
139 : !kind == pg_kind .or. &
140 : !kind == ap_kind .or. &
141 : !kind == an_kind .or. &
142 : !kind == ag_kind .or. &
143 : !kind == general_one_one_kind .or. &
144 : !kind == general_two_one_kind .or. &
145 : !kind == general_two_two_kind .or. &
146 : kind == -1 &
147 : ) kind = other_kind
148 :
149 : !kind = other_kind ! TESTING
150 :
151 : !if (reaction_name(ir) == 'rfe52aprot_to_ni56') then
152 : ! write(*,2) 'rfe52aprot_to_ni56 kind', kind
153 : ! stop
154 : !end if
155 :
156 0 : eps_factor = 1d0
157 : !write(*,*) trim(reaction_name(ir)),kind
158 2480160 : select case(kind)
159 : case (other_kind)
160 : call get1_derivs( &
161 : n, dydt, eps_nuc_MeV, i, eta, ye, temp, den, abar, zbar, &
162 : num_reactions, rate_factors, rtab, itab, &
163 : deriv_flgs, symbolic, just_dydt, &
164 2480160 : ierr)
165 2480160 : i = i+1
166 2480160 : cycle
167 : case (ng_kind)
168 : eps_factor = eps_nuc_cancel_factor
169 0 : call get_basic_2_to_1_derivs(i,ierr)
170 : case (pn_kind)
171 : eps_factor = eps_nuc_cancel_factor
172 0 : call get_basic_2_to_2_derivs(i,ierr)
173 : case (pg_kind)
174 : eps_factor = eps_nuc_cancel_factor
175 0 : call get_basic_2_to_1_derivs(i,ierr)
176 : case (ap_kind)
177 0 : call get_basic_2_to_2_derivs(i,ierr)
178 : case (an_kind)
179 0 : call get_basic_2_to_2_derivs(i,ierr)
180 : case (ag_kind)
181 0 : call get_basic_2_to_1_derivs(i,ierr)
182 : case (general_one_one_kind)
183 0 : call get_general_1_to_1_derivs(i,ierr)
184 : case (general_two_one_kind)
185 0 : call get_general_2_to_1_derivs(i,ierr)
186 : case (general_two_two_kind)
187 0 : call get_general_2_to_2_derivs(i,ierr)
188 : case default
189 2480160 : call mesa_error(__FILE__,__LINE__,'confusion in net wrt reaction kind')
190 : end select
191 0 : i = i+2
192 :
193 : ! icat 16 is burn_fe
194 80002 : if (.false. .and. icat_f == 16 .and. n% logT >= 9.4864903d0 .and. n% logT <= 9.48649039d0) then
195 : write(*,1) trim(category_name(icat_f)) // ' ' // trim(reaction_name(ir)), &
196 : Qconv*(n% eps_nuc_categories(icat_f) - old_eps_nuc_categories_val), &
197 : Qconv*n% eps_nuc_categories(icat_f)
198 : old_eps_nuc_categories_val = n% eps_nuc_categories(icat_f)
199 : end if
200 :
201 :
202 : end do
203 :
204 160004 : if(associated(net_other_net_derivs)) then
205 : call net_other_net_derivs(n, dydt, eps_nuc_MeV, eta, ye, logtemp, temp, den, abar, zbar, &
206 : num_reactions, rate_factors, &
207 0 : symbolic, just_dydt, ierr)
208 : end if
209 :
210 :
211 : contains
212 :
213 0 : subroutine get_general_1_to_1_derivs(i,ierr) ! e.g., 2 c12 -> mg24, 3 he4 -> c12
214 : integer, intent(in) :: i
215 : integer, intent(out) :: ierr
216 :
217 0 : real(qp) :: b, b_f, b_r, rate
218 0 : real(dp) :: d, d_f, d_r, d1, d2, Q, ys_f, ys_r, &
219 0 : d_ysf_dy1, d_ysr_dy2
220 : integer :: c1, c2, i1, i2, o1, o3
221 :
222 : include 'formats'
223 :
224 0 : ierr = 0
225 :
226 : ! forward reaction is c1 i1 -> c2 i2
227 0 : c1 = reaction_inputs(1,ir)
228 0 : i1 = itab(reaction_inputs(2,ir))
229 0 : c2 = reaction_outputs(1,ir)
230 0 : i2 = itab(reaction_outputs(2,ir))
231 :
232 0 : if (symbolic) then
233 0 : n% d_dydt_dy(i1,i1) = 1
234 0 : n% d_dydt_dy(i1,i2) = 1
235 0 : n% d_dydt_dy(i2,i1) = 1
236 0 : n% d_dydt_dy(i2,i2) = 1
237 0 : return
238 : end if
239 :
240 0 : select case(c1)
241 : case (1)
242 0 : ys_f = n% y(i1)
243 0 : d_ysf_dy1 = 1d0
244 : case (2)
245 0 : ys_f = n% y(i1)*n% y(i1)/2d0
246 0 : d_ysf_dy1 = n% y(i1)
247 : case (3)
248 0 : ys_f = n% y(i1)*n% y(i1)*n% y(i1)/6d0
249 0 : d_ysf_dy1 = n% y(i1)*n% y(i1)/2d0
250 : case default
251 0 : write(*,2) 'c1 bad for ' // trim(reaction_name(ir)), c1
252 0 : call mesa_error(__FILE__,__LINE__,'get_general_1_to_1_derivs')
253 : end select
254 0 : d_f = ys_f
255 :
256 0 : select case(c2)
257 : case (1)
258 0 : ys_r = n% y(i2)
259 0 : d_ysr_dy2 = 1d0
260 : case (2)
261 0 : ys_r = n% y(i2)*n% y(i2)/2d0
262 0 : d_ysr_dy2 = n% y(i2)
263 : case (3)
264 0 : ys_r = n% y(i2)*n% y(i2)*n% y(i2)/6d0
265 0 : d_ysr_dy2 = n% y(i2)*n% y(i2)/2d0
266 : case default
267 0 : write(*,2) 'c2 bad for ' // trim(reaction_name(ir)), c2
268 0 : call mesa_error(__FILE__,__LINE__,'get_general_1_to_1_derivs')
269 : end select
270 0 : d_r = ys_r
271 :
272 0 : rate = n% rate_screened(i)
273 0 : b_f = d_f*rate
274 0 : rate = n% rate_screened(r_i)
275 0 : b_r = d_r*rate
276 0 : b = b_f - b_r
277 0 : dydt(i_rate,i1) = dydt(i_rate,i1) - c1*b
278 0 : dydt(i_rate,i2) = dydt(i_rate,i2) + c2*b
279 :
280 0 : n% raw_rate(i) = d_f * n% rate_raw(i) * avo
281 0 : n% raw_rate(r_i) = d_r * n% rate_raw(r_i) * avo
282 :
283 0 : n% screened_rate(i) = d_f * n% rate_screened(i) * avo
284 0 : n% screened_rate(r_i) = d_r * n% rate_screened(r_i) * avo
285 :
286 0 : if (just_dydt) return
287 :
288 0 : Q = n% reaction_Qs(ir)*eps_factor
289 0 : b_f = Q*b_f
290 0 : b_r = -Q*b_r
291 0 : b = b_f + b_r
292 0 : eps_nuc_MeV(i_rate) = eps_nuc_MeV(i_rate) + b
293 0 : n% eps_nuc_categories(icat_f) = n% eps_nuc_categories(icat_f) + b_f
294 0 : n% eps_nuc_categories(icat_r) = n% eps_nuc_categories(icat_r) + b_r
295 : if (show_eps_nuc .and. abs(b) > 1d2) &
296 : write(*,1) trim(reaction_Name(ir)) // ' eps_nuc', b, b_f, b_r
297 :
298 0 : n% eps_nuc_rate(i) = b_f * Qconv
299 0 : n% eps_nuc_rate(r_i) = b_r * Qconv
300 0 : n% eps_neu_rate(i) = 0d0
301 0 : n% eps_neu_rate(r_i) = 0d0
302 :
303 :
304 0 : rate = n% rate_screened_dT(i)
305 0 : b_f = d_f*rate
306 0 : rate = n% rate_screened_dT(r_i)
307 0 : b_r = d_r*rate
308 0 : b = b_f - b_r
309 0 : dydt(i_rate_dT,i1) = dydt(i_rate_dT,i1) - c1*b
310 0 : dydt(i_rate_dT,i2) = dydt(i_rate_dT,i2) + c2*b
311 0 : b_f = Q*b_f
312 0 : b_r = -Q*b_r
313 0 : b = b_f + b_r
314 0 : eps_nuc_MeV(i_rate_dT) = eps_nuc_MeV(i_rate_dT) + b
315 :
316 0 : rate = n% rate_screened_dRho(i)
317 0 : b_f = d_f*rate
318 0 : rate = n% rate_screened_dRho(r_i)
319 0 : b_r = d_r*rate
320 0 : b = b_f - b_r
321 0 : dydt(i_rate_dRho,i1) = dydt(i_rate_dRho,i1) - c1*b
322 0 : dydt(i_rate_dRho,i2) = dydt(i_rate_dRho,i2) + c2*b
323 0 : b_f = Q*b_f
324 0 : b_r = -Q*b_r
325 0 : b = b_f + b_r
326 0 : eps_nuc_MeV(i_rate_dRho) = eps_nuc_MeV(i_rate_dRho) + b
327 :
328 : if (checking_deriv_flags) then
329 : deriv_flgs(i) = .true.
330 : deriv_flgs(r_i) = .true.
331 : end if
332 :
333 0 : d1 = d_ysf_dy1*n% rate_screened(i) ! d(rate_f)/d(y1)
334 0 : d2 = d_ysr_dy2*n% rate_screened(r_i) ! d(rate_r)/d(y2)
335 :
336 0 : n% d_eps_nuc_dy(i1) = n% d_eps_nuc_dy(i1) + Q*d1
337 0 : n% d_eps_nuc_dy(i2) = n% d_eps_nuc_dy(i2) - Q*d2
338 :
339 0 : n% d_dydt_dy(i1,i1) = n% d_dydt_dy(i1,i1) - c1*d1
340 0 : n% d_dydt_dy(i2,i1) = n% d_dydt_dy(i2,i1) + c2*d1
341 :
342 0 : n% d_dydt_dy(i1,i2) = n% d_dydt_dy(i1,i2) + c1*d2
343 0 : n% d_dydt_dy(i2,i2) = n% d_dydt_dy(i2,i2) - c2*d2
344 :
345 : end subroutine get_general_1_to_1_derivs
346 :
347 0 : subroutine get_general_2_to_1_derivs(i,ierr) ! e.g., r_he4_si28_to_o16_o16
348 : integer, intent(in) :: i
349 : integer, intent(out) :: ierr
350 :
351 0 : real(qp) :: b, b_f, b_r, rate
352 0 : real(dp) :: d, d_f, d_r, d1, d2, d3, Q, ys_f, ys_r, &
353 0 : d_ysf_dy1, d_ysf_dy2, d_ysr_dy3, y1, y2, y3
354 : integer :: c1, c2, c3, i1, i2, i3, o1, o3
355 :
356 : include 'formats'
357 :
358 0 : ierr = 0
359 :
360 : ! forward reaction is c1 i1 + c2 i2 -> c3 i3
361 0 : c1 = reaction_inputs(1,ir)
362 0 : i1 = itab(reaction_inputs(2,ir))
363 0 : c2 = reaction_inputs(3,ir)
364 0 : i2 = itab(reaction_inputs(4,ir))
365 0 : c3 = reaction_outputs(1,ir)
366 0 : i3 = itab(reaction_outputs(2,ir))
367 :
368 0 : if (symbolic) then
369 0 : n% d_dydt_dy(i1,i1) = 1
370 0 : n% d_dydt_dy(i1,i2) = 1
371 0 : n% d_dydt_dy(i1,i3) = 1
372 0 : n% d_dydt_dy(i2,i1) = 1
373 0 : n% d_dydt_dy(i2,i2) = 1
374 0 : n% d_dydt_dy(i2,i3) = 1
375 0 : n% d_dydt_dy(i3,i1) = 1
376 0 : n% d_dydt_dy(i3,i2) = 1
377 0 : n% d_dydt_dy(i3,i3) = 1
378 0 : return
379 : end if
380 :
381 0 : y1 = n% y(i1)
382 0 : y2 = n% y(i2)
383 0 : y3 = n% y(i3)
384 :
385 : select case(c1)
386 : case (1)
387 : ys_f = y1
388 0 : d_ysf_dy1 = 1d0
389 : case (2)
390 0 : ys_f = y1*y1/2d0
391 0 : d_ysf_dy1 = y1
392 : case (3)
393 0 : ys_f = y1*y1*y1/6d0
394 0 : d_ysf_dy1 = y1*y1/2d0
395 : case default
396 0 : write(*,2) 'c1 too big for ' // trim(reaction_name(ir))
397 0 : call mesa_error(__FILE__,__LINE__,'get_general_2_to_1_derivs')
398 : end select
399 : ! at this point, ys_f and d_ysf_dy1 only have y1 terms
400 : ! now combine with y2 info
401 0 : select case(c2)
402 : case (1)
403 0 : d_ysf_dy2 = ys_f
404 0 : ys_f = ys_f*y2
405 0 : d_ysf_dy1 = d_ysf_dy1*y2
406 : case (2)
407 0 : d_ysf_dy2 = ys_f*y2
408 0 : ys_f = ys_f*(y2*y2/2d0)
409 0 : d_ysf_dy1 = d_ysf_dy1*(y2*y2/2d0)
410 : case (3)
411 0 : d_ysf_dy2 = ys_f*(y2*y2/2d0)
412 0 : ys_f = ys_f*(y2*y2*y2/6d0)
413 0 : d_ysf_dy1 = d_ysf_dy1*(y2*y2*y2/6d0)
414 : case default
415 0 : write(*,2) 'c1 too big for ' // trim(reaction_name(ir))
416 0 : call mesa_error(__FILE__,__LINE__,'get_general_2_to_1_derivs')
417 : end select
418 :
419 : select case(c3)
420 : case (1)
421 : ys_r = y3
422 0 : d_ysr_dy3 = 1d0
423 : case (2)
424 0 : ys_r = y3*y3/2d0
425 0 : d_ysr_dy3 = y3
426 : case (3)
427 0 : ys_r = y3*y3*y3/6d0
428 0 : d_ysr_dy3 = y3*y3/2d0
429 : case default
430 0 : write(*,2) 'c3 too big for ' // trim(reaction_name(ir))
431 0 : call mesa_error(__FILE__,__LINE__,'get_general_2_to_1_derivs')
432 : end select
433 :
434 0 : d_f = ys_f
435 0 : d_r = ys_r
436 :
437 0 : rate = n% rate_screened(i)
438 0 : b_f = d_f*rate
439 0 : rate = n% rate_screened(r_i)
440 0 : b_r = d_r*rate
441 0 : b = b_f - b_r
442 0 : dydt(i_rate,i1) = dydt(i_rate,i1) - c1*b
443 0 : dydt(i_rate,i2) = dydt(i_rate,i2) - c2*b
444 0 : dydt(i_rate,i3) = dydt(i_rate,i3) + c3*b
445 :
446 0 : n% raw_rate(i) = d_f * n% rate_raw(i) * avo
447 0 : n% raw_rate(r_i) = d_r * n% rate_raw(r_i) * avo
448 :
449 0 : n% screened_rate(i) = d_f * n% rate_screened(i) * avo
450 0 : n% screened_rate(r_i) = d_r * n% rate_screened(r_i) * avo
451 :
452 0 : if (just_dydt) return
453 :
454 0 : Q = n% reaction_Qs(ir)*eps_factor
455 0 : b_f = Q*b_f
456 0 : b_r = -Q*b_r
457 0 : b = b_f + b_r
458 0 : eps_nuc_MeV(i_rate) = eps_nuc_MeV(i_rate) + b
459 0 : n% eps_nuc_categories(icat_f) = n% eps_nuc_categories(icat_f) + b_f
460 0 : n% eps_nuc_categories(icat_r) = n% eps_nuc_categories(icat_r) + b_r
461 : if (show_eps_nuc .and. abs(b) > 1d2) &
462 : write(*,1) trim(reaction_Name(ir)) // ' eps_nuc', b, b_f, b_r
463 :
464 0 : n% eps_nuc_rate(i) = b_f * Qconv
465 0 : n% eps_nuc_rate(r_i) = b_r * Qconv
466 0 : n% eps_neu_rate(i) = 0d0
467 0 : n% eps_neu_rate(r_i) = 0d0
468 :
469 0 : rate = n% rate_screened_dT(i)
470 0 : b_f = d_f*rate
471 0 : rate = n% rate_screened_dT(r_i)
472 0 : b_r = d_r*rate
473 0 : b = b_f - b_r
474 0 : dydt(i_rate_dT,i1) = dydt(i_rate_dT,i1) - c1*b
475 0 : dydt(i_rate_dT,i2) = dydt(i_rate_dT,i2) - c2*b
476 0 : dydt(i_rate_dT,i3) = dydt(i_rate_dT,i3) + c3*b
477 0 : b_f = Q*b_f
478 0 : b_r = -Q*b_r
479 0 : b = b_f + b_r
480 0 : eps_nuc_MeV(i_rate_dT) = eps_nuc_MeV(i_rate_dT) + b
481 :
482 0 : rate = n% rate_screened_dRho(i)
483 0 : b_f = d_f*rate
484 0 : rate = n% rate_screened_dRho(r_i)
485 0 : b_r = d_r*rate
486 0 : b = b_f - b_r
487 0 : dydt(i_rate_dRho,i1) = dydt(i_rate_dRho,i1) - c1*b
488 0 : dydt(i_rate_dRho,i2) = dydt(i_rate_dRho,i2) - c2*b
489 0 : dydt(i_rate_dRho,i3) = dydt(i_rate_dRho,i3) + c3*b
490 0 : b_f = Q*b_f
491 0 : b_r = -Q*b_r
492 0 : b = b_f + b_r
493 0 : eps_nuc_MeV(i_rate_dRho) = eps_nuc_MeV(i_rate_dRho) + b
494 :
495 : if (checking_deriv_flags) then
496 : deriv_flgs(i) = .true.
497 : deriv_flgs(r_i) = .true.
498 : end if
499 :
500 0 : d1 = d_ysf_dy1*n% rate_screened(i) ! d(rate_f)/d(y1)
501 0 : d2 = d_ysf_dy2*n% rate_screened(i) ! d(rate_f)/d(y2)
502 0 : d3 = d_ysr_dy3*n% rate_screened(r_i) ! d(rate_r)/d(y3)
503 :
504 0 : n% d_eps_nuc_dy(i1) = n% d_eps_nuc_dy(i1) + Q*d1
505 0 : n% d_eps_nuc_dy(i2) = n% d_eps_nuc_dy(i2) + Q*d2
506 0 : n% d_eps_nuc_dy(i3) = n% d_eps_nuc_dy(i3) - Q*d3
507 :
508 : ! dydt(1,i1) = dydt(1,i1) - c1*(ys_f*n% rate_screened(i) - ys_r*n% rate_screened(r_i))
509 0 : n% d_dydt_dy(i1,i1) = n% d_dydt_dy(i1,i1) - c1*d1
510 0 : n% d_dydt_dy(i1,i2) = n% d_dydt_dy(i1,i2) - c1*d2
511 0 : n% d_dydt_dy(i1,i3) = n% d_dydt_dy(i1,i3) + c1*d3
512 :
513 : ! dydt(1,i2) = dydt(1,i2) - c2*(ys_f*n% rate_screened(i) - ys_r*n% rate_screened(r_i))
514 0 : n% d_dydt_dy(i2,i1) = n% d_dydt_dy(i2,i1) - c2*d1
515 0 : n% d_dydt_dy(i2,i2) = n% d_dydt_dy(i2,i2) - c2*d2
516 0 : n% d_dydt_dy(i2,i3) = n% d_dydt_dy(i2,i3) + c2*d3
517 :
518 : ! dydt(1,i3) = dydt(1,i3) + c3*(ys_f*n% rate_screened(i) - ys_r*n% rate_screened(r_i))
519 0 : n% d_dydt_dy(i3,i1) = n% d_dydt_dy(i3,i1) + c3*d1
520 0 : n% d_dydt_dy(i3,i2) = n% d_dydt_dy(i3,i2) + c3*d2
521 0 : n% d_dydt_dy(i3,i3) = n% d_dydt_dy(i3,i3) - c3*d3
522 :
523 : if (.false. .and. reaction_name(ir) == 'r_he4_si28_to_o16_o16') then ! .and. &
524 : !y1 > 1d-20 .and. y2 > 1d-20 .and. y3 > 1d-20) then
525 : write(*,'(A)')
526 : write(*,2) trim(reaction_name(ir))
527 : write(*,2) trim(reaction_name(r_ir))
528 : write(*,'(A)')
529 : write(*,2) 'c1', c1
530 : write(*,2) 'c2', c2
531 : write(*,2) 'c3', c3
532 : write(*,'(A)')
533 : write(*,1) 'd1', d1
534 : write(*,1) 'd2', d2
535 : write(*,'(A)')
536 : write(*,1) 'dr1', d_ysf_dy1
537 : write(*,1) 'dr2', d_ysf_dy2
538 : write(*,'(A)')
539 : write(*,1) 'd_ysf_dy1', d_ysf_dy1
540 : write(*,1) 'd_ysf_dy2', d_ysf_dy2
541 : write(*,'(A)')
542 : write(*,1) 'y1', y1
543 : write(*,1) 'y2', y2
544 : write(*,1) 'y3', y3
545 : write(*,'(A)')
546 : write(*,1) 'd_dydt_dy(i1,i1)', - c1*d1
547 : write(*,1) 'd_dydt_dy(i1,i2)', - c1*d2
548 : write(*,1) 'd_dydt_dy(i1,i3)', c1*d3
549 : write(*,'(A)')
550 : write(*,1) 'd_dydt_dy(i2,i1)', - c2*d1
551 : write(*,1) 'd_dydt_dy(i2,i2)', - c2*d2
552 : write(*,1) 'd_dydt_dy(i2,i3)', c2*d3
553 : write(*,'(A)')
554 : write(*,1) 'd_dydt_dy(i3,i1)', c3*d1
555 : write(*,1) 'd_dydt_dy(i3,i2)', c3*d2
556 : write(*,1) 'd_dydt_dy(i3,i3)', -c3*d3
557 : write(*,'(A)')
558 : call mesa_error(__FILE__,__LINE__,'get_general_2_to_1_derivs')
559 : end if
560 :
561 : end subroutine get_general_2_to_1_derivs
562 :
563 0 : subroutine get_general_2_to_2_derivs(i,ierr)
564 : integer, intent(in) :: i
565 : integer, intent(out) :: ierr
566 :
567 0 : real(qp) :: b, b_f, b_r, rate
568 0 : real(dp) :: d, d_f, d_r, d1, d2, d3, d4, Q, ys_f, ys_r, &
569 0 : d_ysf_dy1, d_ysf_dy2, d_ysr_dy3, d_ysr_dy4, y1, y2, y3, y4
570 : integer :: c1, c2, c3, c4, i1, i2, i3, i4
571 :
572 : include 'formats'
573 :
574 0 : ierr = 0
575 :
576 : ! forward reaction is c1 i1 + c2 i2 -> c3 i3 + c4 i4
577 0 : c1 = reaction_inputs(1,ir)
578 0 : i1 = itab(reaction_inputs(2,ir))
579 0 : c2 = reaction_inputs(3,ir)
580 0 : i2 = itab(reaction_inputs(4,ir))
581 0 : c3 = reaction_outputs(1,ir)
582 0 : i3 = itab(reaction_outputs(2,ir))
583 0 : c4 = reaction_outputs(3,ir)
584 0 : i4 = itab(reaction_outputs(4,ir))
585 :
586 0 : if (symbolic) then
587 0 : n% d_dydt_dy(i1,i1) = 1
588 0 : n% d_dydt_dy(i1,i2) = 1
589 0 : n% d_dydt_dy(i1,i3) = 1
590 0 : n% d_dydt_dy(i1,i4) = 1
591 :
592 0 : n% d_dydt_dy(i2,i1) = 1
593 0 : n% d_dydt_dy(i2,i2) = 1
594 0 : n% d_dydt_dy(i2,i3) = 1
595 0 : n% d_dydt_dy(i2,i4) = 1
596 :
597 0 : n% d_dydt_dy(i3,i1) = 1
598 0 : n% d_dydt_dy(i3,i2) = 1
599 0 : n% d_dydt_dy(i3,i3) = 1
600 0 : n% d_dydt_dy(i3,i4) = 1
601 :
602 0 : n% d_dydt_dy(i4,i1) = 1
603 0 : n% d_dydt_dy(i4,i2) = 1
604 0 : n% d_dydt_dy(i4,i3) = 1
605 0 : n% d_dydt_dy(i4,i4) = 1
606 0 : return
607 : end if
608 :
609 0 : y1 = n% y(i1)
610 0 : y2 = n% y(i2)
611 0 : y3 = n% y(i3)
612 0 : y4 = n% y(i4)
613 :
614 : select case(c1)
615 : case (1)
616 : ys_f = y1
617 0 : d_ysf_dy1 = 1d0
618 : case (2)
619 0 : ys_f = y1*y1/2d0
620 0 : d_ysf_dy1 = y1
621 : case (3)
622 0 : ys_f = y1*y1*y1/6d0
623 0 : d_ysf_dy1 = y1*y1/2d0
624 : case default
625 0 : write(*,2) 'c1 too big for ' // trim(reaction_name(ir)), c1
626 0 : call mesa_error(__FILE__,__LINE__,'get_general_2_to_2_derivs')
627 : end select
628 : ! at this point, ys_f and d_ysf_dy1 only have y1 terms
629 : ! now combine with y2 info
630 0 : select case(c2)
631 : case (1)
632 0 : d_ysf_dy2 = ys_f
633 0 : ys_f = ys_f*y2
634 0 : d_ysf_dy1 = d_ysf_dy1*y2
635 : case (2)
636 0 : d_ysf_dy2 = ys_f*y2
637 0 : ys_f = ys_f*(y2*y2/2d0)
638 0 : d_ysf_dy1 = d_ysf_dy1*(y2*y2/2d0)
639 : case (3)
640 0 : d_ysf_dy2 = ys_f*(y2*y2/2d0)
641 0 : ys_f = ys_f*(y2*y2*y2/6d0)
642 0 : d_ysf_dy1 = d_ysf_dy1*(y2*y2*y2/6d0)
643 : case default
644 0 : write(*,2) 'c2 too big for ' // trim(reaction_name(ir)), c2
645 0 : call mesa_error(__FILE__,__LINE__,'get_general_2_to_2_derivs')
646 : end select
647 :
648 : select case(c3)
649 : case (1)
650 : ys_r = y3
651 0 : d_ysr_dy3 = 1d0
652 : case (2)
653 0 : ys_r = y3*y3/2d0
654 0 : d_ysr_dy3 = y3
655 : case (3)
656 0 : ys_r = y3*y3*y3/6d0
657 0 : d_ysr_dy3 = y3*y3/2d0
658 : case default
659 0 : write(*,2) 'c3 too big for ' // trim(reaction_name(ir)), c3
660 0 : call mesa_error(__FILE__,__LINE__,'get_general_2_to_2_derivs')
661 : end select
662 : ! at this point, ys_r and d_ysf_dy3 only have y3 terms
663 : ! now combine with y4 info
664 0 : select case(c4)
665 : case (1)
666 0 : d_ysr_dy4 = ys_r
667 0 : ys_r = ys_r*y4
668 0 : d_ysr_dy3 = d_ysr_dy3*y4
669 : case (2)
670 0 : d_ysr_dy4 = ys_r*y4
671 0 : ys_r = ys_r*(y4*y4/2d0)
672 0 : d_ysr_dy3 = d_ysr_dy3*(y4*y4/2d0)
673 : case (3)
674 0 : d_ysr_dy4 = ys_r*(y4*y4/2d0)
675 0 : ys_r = ys_r*(y4*y4*y4/6d0)
676 0 : d_ysr_dy3 = d_ysr_dy3*(y4*y4*y4/6d0)
677 : case default
678 0 : write(*,2) 'c4 too big for ' // trim(reaction_name(ir))
679 0 : call mesa_error(__FILE__,__LINE__,'get_general_2_to_2_derivs')
680 : end select
681 :
682 0 : d_f = ys_f
683 0 : d_r = ys_r
684 :
685 0 : rate = n% rate_screened(i)
686 0 : b_f = d_f*rate
687 0 : rate = n% rate_screened(r_i)
688 0 : b_r = d_r*rate
689 0 : b = b_f - b_r
690 0 : dydt(i_rate,i1) = dydt(i_rate,i1) - c1*b
691 0 : dydt(i_rate,i2) = dydt(i_rate,i2) - c2*b
692 0 : dydt(i_rate,i3) = dydt(i_rate,i3) + c3*b
693 0 : dydt(i_rate,i4) = dydt(i_rate,i4) + c4*b
694 :
695 0 : n% raw_rate(i) = d_f * n% rate_raw(i) * avo
696 0 : n% raw_rate(r_i) = d_r * n% rate_raw(r_i) * avo
697 :
698 0 : n% screened_rate(i) = d_f * n% rate_screened(i) * avo
699 0 : n% screened_rate(r_i) = d_r * n% rate_screened(r_i) * avo
700 :
701 0 : if (just_dydt) return
702 :
703 0 : Q = n% reaction_Qs(ir)*eps_factor
704 0 : b_f = Q*b_f
705 0 : b_r = -Q*b_r
706 0 : b = b_f + b_r
707 0 : eps_nuc_MeV(i_rate) = eps_nuc_MeV(i_rate) + b
708 0 : n% eps_nuc_categories(icat_f) = n% eps_nuc_categories(icat_f) + b_f
709 0 : n% eps_nuc_categories(icat_r) = n% eps_nuc_categories(icat_r) + b_r
710 : if (show_eps_nuc .and. abs(b) > 1d2) &
711 : write(*,1) trim(reaction_Name(ir)) // ' eps_nuc', b, b_f, b_r
712 :
713 0 : n% eps_nuc_rate(i) = b_f * Qconv
714 0 : n% eps_nuc_rate(r_i) = b_r * Qconv
715 0 : n% eps_neu_rate(i) = 0d0
716 0 : n% eps_neu_rate(r_i) = 0d0
717 :
718 0 : rate = n% rate_screened_dT(i)
719 0 : b_f = d_f*rate
720 0 : rate = n% rate_screened_dT(r_i)
721 0 : b_r = d_r*rate
722 0 : b = b_f - b_r
723 0 : dydt(i_rate_dT,i1) = dydt(i_rate_dT,i1) - c1*b
724 0 : dydt(i_rate_dT,i2) = dydt(i_rate_dT,i2) - c2*b
725 0 : dydt(i_rate_dT,i3) = dydt(i_rate_dT,i3) + c3*b
726 0 : dydt(i_rate_dT,i4) = dydt(i_rate_dT,i4) + c4*b
727 0 : b_f = Q*b_f
728 0 : b_r = -Q*b_r
729 0 : b = b_f + b_r
730 0 : eps_nuc_MeV(i_rate_dT) = eps_nuc_MeV(i_rate_dT) + b
731 :
732 0 : rate = n% rate_screened_dRho(i)
733 0 : b_f = d_f*rate
734 0 : rate = n% rate_screened_dRho(r_i)
735 0 : b_r = d_r*rate
736 0 : b = b_f - b_r
737 0 : dydt(i_rate_dRho,i1) = dydt(i_rate_dRho,i1) - c1*b
738 0 : dydt(i_rate_dRho,i2) = dydt(i_rate_dRho,i2) - c2*b
739 0 : dydt(i_rate_dRho,i3) = dydt(i_rate_dRho,i3) + c3*b
740 0 : dydt(i_rate_dRho,i4) = dydt(i_rate_dRho,i4) + c4*b
741 0 : b_f = Q*b_f
742 0 : b_r = -Q*b_r
743 0 : b = b_f + b_r
744 0 : eps_nuc_MeV(i_rate_dRho) = eps_nuc_MeV(i_rate_dRho) + b
745 :
746 : if (checking_deriv_flags) then
747 : deriv_flgs(i) = .true.
748 : deriv_flgs(r_i) = .true.
749 : end if
750 :
751 0 : d1 = d_ysf_dy1*n% rate_screened(i) ! d(rate_f)/d(y1)
752 0 : d2 = d_ysf_dy2*n% rate_screened(i) ! d(rate_f)/d(y2)
753 0 : d3 = d_ysr_dy3*n% rate_screened(r_i) ! d(rate_r)/d(y3)
754 0 : d4 = d_ysr_dy4*n% rate_screened(r_i) ! d(rate_r)/d(y4)
755 :
756 0 : n% d_eps_nuc_dy(i1) = n% d_eps_nuc_dy(i1) + Q*d1
757 0 : n% d_eps_nuc_dy(i2) = n% d_eps_nuc_dy(i2) + Q*d2
758 0 : n% d_eps_nuc_dy(i3) = n% d_eps_nuc_dy(i3) - Q*d3
759 0 : n% d_eps_nuc_dy(i4) = n% d_eps_nuc_dy(i4) - Q*d4
760 :
761 0 : n% d_dydt_dy(i1,i1) = n% d_dydt_dy(i1,i1) - c1*d1
762 0 : n% d_dydt_dy(i1,i2) = n% d_dydt_dy(i1,i2) - c1*d2
763 0 : n% d_dydt_dy(i1,i3) = n% d_dydt_dy(i1,i3) + c1*d3
764 0 : n% d_dydt_dy(i1,i4) = n% d_dydt_dy(i1,i4) + c1*d4
765 :
766 0 : n% d_dydt_dy(i2,i1) = n% d_dydt_dy(i2,i1) - c2*d1
767 0 : n% d_dydt_dy(i2,i2) = n% d_dydt_dy(i2,i2) - c2*d2
768 0 : n% d_dydt_dy(i2,i3) = n% d_dydt_dy(i2,i3) + c2*d3
769 0 : n% d_dydt_dy(i2,i4) = n% d_dydt_dy(i2,i4) + c2*d4
770 :
771 0 : n% d_dydt_dy(i3,i1) = n% d_dydt_dy(i3,i1) + c3*d1
772 0 : n% d_dydt_dy(i3,i2) = n% d_dydt_dy(i3,i2) + c3*d2
773 0 : n% d_dydt_dy(i3,i3) = n% d_dydt_dy(i3,i3) - c3*d3
774 0 : n% d_dydt_dy(i3,i4) = n% d_dydt_dy(i3,i4) - c3*d4
775 :
776 0 : n% d_dydt_dy(i4,i1) = n% d_dydt_dy(i4,i1) + c4*d1
777 0 : n% d_dydt_dy(i4,i2) = n% d_dydt_dy(i4,i2) + c4*d2
778 0 : n% d_dydt_dy(i4,i3) = n% d_dydt_dy(i4,i3) - c4*d3
779 0 : n% d_dydt_dy(i4,i4) = n% d_dydt_dy(i4,i4) - c4*d4
780 :
781 : end subroutine get_general_2_to_2_derivs
782 :
783 0 : subroutine get_basic_2_to_2_derivs(i,ierr)
784 : integer, intent(in) :: i
785 : integer, intent(out) :: ierr
786 :
787 0 : real(qp) :: b, b_f, b_r, rate
788 0 : real(dp) :: d_f, d_r, e_f, e_r, d1, d2, d3, d4, &
789 0 : ys_f, ys_r, Q, y1, y2, y3, y4
790 : integer :: i1, i2, i3, i4, o1, o2, o3, o4
791 :
792 : include 'formats'
793 :
794 0 : ierr = 0
795 :
796 : ! forward reaction is i1 + i2 -> i3 + i4
797 0 : i1 = itab(reaction_inputs(2,ir))
798 0 : i2 = itab(reaction_inputs(4,ir))
799 0 : i3 = itab(reaction_outputs(2,ir))
800 0 : i4 = itab(reaction_outputs(4,ir))
801 :
802 0 : if (symbolic) then
803 0 : n% d_dydt_dy(i1,i1) = 1
804 0 : n% d_dydt_dy(i1,i2) = 1
805 0 : n% d_dydt_dy(i1,i3) = 1
806 0 : n% d_dydt_dy(i1,i4) = 1
807 0 : n% d_dydt_dy(i2,i1) = 1
808 0 : n% d_dydt_dy(i2,i2) = 1
809 0 : n% d_dydt_dy(i2,i3) = 1
810 0 : n% d_dydt_dy(i2,i4) = 1
811 0 : n% d_dydt_dy(i3,i1) = 1
812 0 : n% d_dydt_dy(i3,i2) = 1
813 0 : n% d_dydt_dy(i3,i3) = 1
814 0 : n% d_dydt_dy(i3,i4) = 1
815 0 : n% d_dydt_dy(i4,i1) = 1
816 0 : n% d_dydt_dy(i4,i2) = 1
817 0 : n% d_dydt_dy(i4,i3) = 1
818 0 : n% d_dydt_dy(i4,i4) = 1
819 0 : return
820 : end if
821 :
822 0 : y1 = n% y(i1)
823 0 : y2 = n% y(i2)
824 0 : y3 = n% y(i3)
825 0 : y4 = n% y(i4)
826 :
827 0 : ys_f = y1*y2
828 0 : d_f = ys_f
829 :
830 0 : ys_r = y3*y4
831 0 : d_r = ys_r
832 :
833 0 : rate = n% rate_screened(i)
834 0 : b_f = d_f*rate
835 0 : rate = n% rate_screened(r_i)
836 0 : b_r = d_r*rate
837 0 : b = b_f - b_r
838 0 : dydt(i_rate,i1) = dydt(i_rate,i1) - b
839 0 : dydt(i_rate,i2) = dydt(i_rate,i2) - b
840 0 : dydt(i_rate,i3) = dydt(i_rate,i3) + b
841 0 : dydt(i_rate,i4) = dydt(i_rate,i4) + b
842 :
843 0 : n% raw_rate(i) = d_f * n% rate_raw(i) * avo
844 0 : n% raw_rate(r_i) = d_r * n% rate_raw(r_i) * avo
845 :
846 0 : n% screened_rate(i) = d_f * n% rate_screened(i) * avo
847 0 : n% screened_rate(r_i) = d_r * n% rate_screened(r_i) * avo
848 :
849 0 : if (just_dydt) return
850 :
851 0 : Q = n% reaction_Qs(ir)*eps_factor
852 0 : eps_nuc_MeV(i_rate) = eps_nuc_MeV(i_rate) + b*Q
853 0 : n% eps_nuc_categories(icat_f) = n% eps_nuc_categories(icat_f) + b_f*Q
854 0 : n% eps_nuc_categories(icat_r) = n% eps_nuc_categories(icat_r) - b_r*Q
855 : if (show_eps_nuc .and. abs(b) > 1d2) &
856 : write(*,1) trim(reaction_Name(ir)) // ' eps_nuc', b, b_f, b_r
857 :
858 0 : n% eps_nuc_rate(i) = b_f * Q * Qconv
859 0 : n% eps_nuc_rate(r_i) = -b_r * Q * Qconv
860 0 : n% eps_neu_rate(i) = 0d0
861 0 : n% eps_neu_rate(r_i) = 0d0
862 :
863 :
864 0 : rate = n% rate_screened_dT(i)
865 0 : b_f = d_f*rate
866 0 : rate = n% rate_screened_dT(r_i)
867 0 : b_r = d_r*rate
868 0 : b = b_f - b_r
869 0 : dydt(i_rate_dT,i1) = dydt(i_rate_dT,i1) - b
870 0 : dydt(i_rate_dT,i2) = dydt(i_rate_dT,i2) - b
871 0 : dydt(i_rate_dT,i3) = dydt(i_rate_dT,i3) + b
872 0 : dydt(i_rate_dT,i4) = dydt(i_rate_dT,i4) + b
873 0 : eps_nuc_MeV(i_rate_dT) = eps_nuc_MeV(i_rate_dT) + b*Q
874 :
875 0 : rate = n% rate_screened_dRho(i)
876 0 : b_f = d_f*rate
877 0 : rate = n% rate_screened_dRho(r_i)
878 0 : b_r = d_r*rate
879 0 : b = b_f - b_r
880 0 : dydt(i_rate_dRho,i1) = dydt(i_rate_dRho,i1) - b
881 0 : dydt(i_rate_dRho,i2) = dydt(i_rate_dRho,i2) - b
882 0 : dydt(i_rate_dRho,i3) = dydt(i_rate_dRho,i3) + b
883 0 : dydt(i_rate_dRho,i4) = dydt(i_rate_dRho,i4) + b
884 0 : eps_nuc_MeV(i_rate_dRho) = eps_nuc_MeV(i_rate_dRho) + b*Q
885 :
886 : if (checking_deriv_flags) then
887 : deriv_flgs(i) = .true.
888 : deriv_flgs(r_i) = .true.
889 : end if
890 :
891 0 : e_f = n% rate_screened(i)
892 0 : e_r = n% rate_screened(r_i)
893 :
894 0 : d1 = y2*e_f ! d(rate_f)/d(y1)
895 0 : d2 = y1*e_f ! d(rate_f)/d(y2)
896 0 : d3 = y4*e_r ! d(rate_r)/d(y3)
897 0 : d4 = y3*e_r ! d(rate_r)/d(y4)
898 :
899 0 : n% d_eps_nuc_dy(i1) = n% d_eps_nuc_dy(i1) + Q*d1
900 0 : n% d_eps_nuc_dy(i2) = n% d_eps_nuc_dy(i2) + Q*d2
901 0 : n% d_eps_nuc_dy(i3) = n% d_eps_nuc_dy(i3) - Q*d3
902 0 : n% d_eps_nuc_dy(i4) = n% d_eps_nuc_dy(i4) - Q*d4
903 :
904 0 : n% d_dydt_dy(i1,i1) = n% d_dydt_dy(i1,i1) - d1
905 0 : n% d_dydt_dy(i2,i1) = n% d_dydt_dy(i2,i1) - d1
906 0 : n% d_dydt_dy(i3,i1) = n% d_dydt_dy(i3,i1) + d1
907 0 : n% d_dydt_dy(i4,i1) = n% d_dydt_dy(i4,i1) + d1
908 :
909 0 : n% d_dydt_dy(i1,i2) = n% d_dydt_dy(i1,i2) - d2
910 0 : n% d_dydt_dy(i2,i2) = n% d_dydt_dy(i2,i2) - d2
911 0 : n% d_dydt_dy(i3,i2) = n% d_dydt_dy(i3,i2) + d2
912 0 : n% d_dydt_dy(i4,i2) = n% d_dydt_dy(i4,i2) + d2
913 :
914 0 : n% d_dydt_dy(i1,i3) = n% d_dydt_dy(i1,i3) + d3
915 0 : n% d_dydt_dy(i2,i3) = n% d_dydt_dy(i2,i3) + d3
916 0 : n% d_dydt_dy(i3,i3) = n% d_dydt_dy(i3,i3) - d3
917 0 : n% d_dydt_dy(i4,i3) = n% d_dydt_dy(i4,i3) - d3
918 :
919 0 : n% d_dydt_dy(i1,i4) = n% d_dydt_dy(i1,i4) + d4
920 0 : n% d_dydt_dy(i2,i4) = n% d_dydt_dy(i2,i4) + d4
921 0 : n% d_dydt_dy(i3,i4) = n% d_dydt_dy(i3,i4) - d4
922 0 : n% d_dydt_dy(i4,i4) = n% d_dydt_dy(i4,i4) - d4
923 :
924 : end subroutine get_basic_2_to_2_derivs
925 :
926 0 : subroutine get_basic_2_to_1_derivs(i,ierr)
927 : integer, intent(in) :: i
928 : integer, intent(out) :: ierr
929 :
930 0 : real(qp) :: b, b_f, b_r, rate
931 0 : real(dp) :: d_f, d_r, e_f, e_r, d1, d2, d3, &
932 0 : ys_f, ys_r, Q, y1, y2, y3
933 : integer :: i1, i2, i3, o1, o2, o3, k
934 :
935 : include 'formats'
936 :
937 0 : ierr = 0
938 :
939 : ! forward reaction is i1 + i2 -> i3
940 0 : i1 = itab(reaction_inputs(2,ir))
941 0 : i2 = itab(reaction_inputs(4,ir))
942 0 : i3 = itab(reaction_outputs(2,ir))
943 :
944 : ! if (reaction_inputs(1,ir) /= 1 .or. &
945 : ! reaction_inputs(3,ir) /= 1 .or. &
946 : ! reaction_outputs(1,ir) /= 1 .or. &
947 : ! reaction_inputs(6,ir) /= 0 .or. &
948 : ! reaction_outputs(4,ir) /= 0) then
949 : ! write(*,*) 'bad reaction for _ag_ ' // trim(reaction_name(ir))
950 : ! stop
951 : ! end if
952 :
953 0 : if (symbolic) then
954 0 : n% d_dydt_dy(i1,i1) = 1
955 0 : n% d_dydt_dy(i1,i2) = 1
956 0 : n% d_dydt_dy(i1,i3) = 1
957 0 : n% d_dydt_dy(i2,i1) = 1
958 0 : n% d_dydt_dy(i2,i2) = 1
959 0 : n% d_dydt_dy(i2,i3) = 1
960 0 : n% d_dydt_dy(i3,i1) = 1
961 0 : n% d_dydt_dy(i3,i2) = 1
962 0 : n% d_dydt_dy(i3,i3) = 1
963 0 : return
964 : end if
965 :
966 0 : y1 = n% y(i1)
967 0 : y2 = n% y(i2)
968 0 : y3 = n% y(i3)
969 :
970 0 : ys_f = y1*y2
971 0 : d_f = ys_f
972 :
973 0 : ys_r = y3
974 0 : d_r = ys_r
975 :
976 0 : rate = n% rate_screened(i)
977 0 : b_f = d_f*rate
978 0 : rate = n% rate_screened(r_i)
979 0 : b_r = d_r*rate
980 0 : b = b_f - b_r
981 0 : dydt(i_rate,i1) = dydt(i_rate,i1) - b
982 0 : dydt(i_rate,i2) = dydt(i_rate,i2) - b
983 0 : dydt(i_rate,i3) = dydt(i_rate,i3) + b
984 :
985 0 : n% raw_rate(i) = d_f * n% rate_raw(i) * avo
986 0 : n% raw_rate(r_i) = d_r * n% rate_raw(r_i) * avo
987 :
988 0 : n% screened_rate(i) = d_f * n% rate_screened(i) * avo
989 0 : n% screened_rate(r_i) = d_r * n% rate_screened(r_i) * avo
990 :
991 0 : if (just_dydt) return
992 :
993 0 : Q = n% reaction_Qs(ir)*eps_factor
994 0 : eps_nuc_MeV(i_rate) = eps_nuc_MeV(i_rate) + b*Q
995 0 : n% eps_nuc_categories(icat_f) = n% eps_nuc_categories(icat_f) + b_f*Q
996 0 : n% eps_nuc_categories(icat_r) = n% eps_nuc_categories(icat_r) - b_r*Q
997 : if (show_eps_nuc .and. abs(b) > 1d2) &
998 : write(*,1) trim(reaction_Name(ir)) // ' eps_nuc', b, b_f, b_r
999 :
1000 0 : n% eps_nuc_rate(i) = b_f * Q * Qconv
1001 0 : n% eps_nuc_rate(r_i) = -b_r * Q * Qconv
1002 0 : n% eps_neu_rate(i) = 0d0
1003 0 : n% eps_neu_rate(r_i) = 0d0
1004 :
1005 0 : rate = n% rate_screened_dT(i)
1006 0 : b_f = d_f*rate
1007 0 : rate = n% rate_screened_dT(r_i)
1008 0 : b_r = d_r*rate
1009 0 : b = b_f - b_r
1010 0 : dydt(i_rate_dT,i1) = dydt(i_rate_dT,i1) - b
1011 0 : dydt(i_rate_dT,i2) = dydt(i_rate_dT,i2) - b
1012 0 : dydt(i_rate_dT,i3) = dydt(i_rate_dT,i3) + b
1013 0 : eps_nuc_MeV(i_rate_dT) = eps_nuc_MeV(i_rate_dT) + b*Q
1014 :
1015 0 : rate = n% rate_screened_dRho(i)
1016 0 : b_f = d_f*rate
1017 0 : rate = n% rate_screened_dRho(r_i)
1018 0 : b_r = d_r*rate
1019 0 : b = b_f - b_r
1020 0 : dydt(i_rate_dRho,i1) = dydt(i_rate_dRho,i1) - b
1021 0 : dydt(i_rate_dRho,i2) = dydt(i_rate_dRho,i2) - b
1022 0 : dydt(i_rate_dRho,i3) = dydt(i_rate_dRho,i3) + b
1023 0 : eps_nuc_MeV(i_rate_dRho) = eps_nuc_MeV(i_rate_dRho) + b*Q
1024 :
1025 : if (checking_deriv_flags) then
1026 : deriv_flgs(i) = .true.
1027 : deriv_flgs(r_i) = .true.
1028 : end if
1029 :
1030 0 : e_f = n% rate_screened(i)
1031 0 : e_r = n% rate_screened(r_i)
1032 :
1033 0 : d1 = y2*e_f ! d(rate_f)/d(y1)
1034 0 : d2 = y1*e_f ! d(rate_f)/d(y2)
1035 0 : d3 = e_r ! d(rate_r)/d(y3)
1036 :
1037 0 : n% d_eps_nuc_dy(i1) = n% d_eps_nuc_dy(i1) + Q*d1
1038 0 : n% d_eps_nuc_dy(i2) = n% d_eps_nuc_dy(i2) + Q*d2
1039 0 : n% d_eps_nuc_dy(i3) = n% d_eps_nuc_dy(i3) - Q*d3
1040 :
1041 0 : n% d_dydt_dy(i1,i1) = n% d_dydt_dy(i1,i1) - d1
1042 0 : n% d_dydt_dy(i2,i1) = n% d_dydt_dy(i2,i1) - d1
1043 0 : n% d_dydt_dy(i3,i1) = n% d_dydt_dy(i3,i1) + d1
1044 :
1045 0 : n% d_dydt_dy(i1,i2) = n% d_dydt_dy(i1,i2) - d2
1046 0 : n% d_dydt_dy(i2,i2) = n% d_dydt_dy(i2,i2) - d2
1047 0 : n% d_dydt_dy(i3,i2) = n% d_dydt_dy(i3,i2) + d2
1048 :
1049 0 : n% d_dydt_dy(i1,i3) = n% d_dydt_dy(i1,i3) + d3
1050 0 : n% d_dydt_dy(i2,i3) = n% d_dydt_dy(i2,i3) + d3
1051 0 : n% d_dydt_dy(i3,i3) = n% d_dydt_dy(i3,i3) - d3
1052 :
1053 : end subroutine get_basic_2_to_1_derivs
1054 :
1055 : subroutine Check
1056 : integer :: nrates
1057 : nrates = n% g% num_reactions
1058 :
1059 : do ir = 1, nrates
1060 : if (.not. deriv_flgs(ir)) then
1061 : all_okay = .false.
1062 : write(*,'(a,i4,2x,a)') 'missing derivs for ', ir, &
1063 : trim(reaction_Name(g% reaction_id(ir)))
1064 : end if
1065 : end do
1066 :
1067 : end subroutine Check
1068 :
1069 :
1070 : end subroutine get_derivs
1071 :
1072 :
1073 2480160 : subroutine get1_derivs( &
1074 2480160 : n, dydt, eps_nuc_MeV, i, eta, ye, temp, den, abar, zbar, &
1075 : num_reactions, rate_factors, rtab, itab, &
1076 : deriv_flgs, symbolic, just_dydt, ierr)
1077 : use rates_lib, only: eval_n14_electron_capture_rate
1078 : type (Net_Info) :: n
1079 : integer, intent(in) :: i, num_reactions
1080 : real(qp), intent(inout) :: dydt(:,:)
1081 : real(qp), intent(out) :: eps_nuc_MeV(num_rvs)
1082 : real(dp), intent(in) :: eta, ye, temp, den, abar, zbar, rate_factors(:)
1083 : integer, pointer, intent(in) :: rtab(:), itab(:)
1084 : logical, pointer :: deriv_flgs(:)
1085 : logical, intent(in) :: symbolic, just_dydt
1086 : integer, intent(out) :: ierr
1087 :
1088 : integer :: ir, j, prot, neut, h1, he4, c14, n14, ne20, ne22, &
1089 : mg21, mg22, mg23, mg24, al23, al24, si24, si25, si26, &
1090 : s28, s29, s30, cl31, ar32, ar33, ar34, k35, ca36, ca37, ca38, &
1091 : ti41, ti42, v43, cr44, cr45, cr46, cr56, &
1092 : fe48, fe49, fe50, fe51, fe52, fe54, fe56, fe58, fe60, fe62, fe64, &
1093 : ni52, ni53, ni54, ni55, ni56, ni58, ni60, ni62, &
1094 : zn57, zn58, zn59, zn60, ge62, ge63, ge64, &
1095 : se68, kr72, sr76, mo84, sn104
1096 : integer :: weak_id, num_reaction_inputs, in1, in2, in3, in4, in5
1097 : integer :: cin1, cin2, cin3, cin4, cin5
1098 9920640 : real(dp) :: din1, din2, din3, din4, din5
1099 : integer :: num_reaction_outputs, out1, out2, out3, out4, out5
1100 : integer :: cout1, cout2, cout3, cout4, cout5
1101 9920640 : real(dp) :: dout1, dout2, dout3, dout4, dout5
1102 : type (Net_General_Info), pointer :: g
1103 2480160 : integer, pointer :: reaction_id(:)
1104 : integer :: i1, i2, i3, idr1, idr2, idr3, o1, o2, o3
1105 2480160 : real(dp) :: r, dr1, dr2, dr3, rate, d_rate_dlnT, d_rate_dlnRho, Q, Qneu, rn14ec
1106 :
1107 : integer, dimension(3) :: i_in, i_out, idr
1108 24801600 : real(dp), dimension(3) :: c_in, c_out, dr
1109 :
1110 : logical :: done, has_prot, has_neut, has_h1, switch_to_prot
1111 : integer :: max_Z, Z_plus_N_for_max_Z
1112 : integer, parameter :: min_Z_for_switch_to_prot = 12
1113 :
1114 : logical, parameter :: dbg = .false.
1115 :
1116 : include 'formats'
1117 :
1118 2480160 : g => n% g
1119 2480160 : reaction_id => g% reaction_id
1120 :
1121 2480160 : ierr = 0
1122 2480160 : prot = itab(iprot)
1123 2480160 : neut = itab(ineut)
1124 2480160 : h1 = itab(ih1)
1125 2480160 : he4 = itab(ihe4)
1126 :
1127 2480160 : ir = reaction_id(i)
1128 :
1129 2480160 : if (reaction_outputs(1,ir) == 0) then
1130 1040012 : n% raw_rate(i) = 0d0
1131 1040012 : n% screened_rate(i) = 0d0
1132 1040012 : n% eps_nuc_rate(i) = 0d0
1133 1040012 : n% eps_neu_rate(i) = 0d0
1134 1120010 : return ! skip aux reactions
1135 : end if
1136 :
1137 :
1138 : if (dbg) &
1139 : write(*,'(/,a,2i6)') ' reaction name <' // trim(reaction_Name(ir)) // '>', i, ir
1140 :
1141 1440148 : max_Z = g% reaction_max_Z(i)
1142 1440148 : Z_plus_N_for_max_Z = g% reaction_max_Z_plus_N_for_max_Z(i)
1143 :
1144 1440148 : din1 = 0d0
1145 1440148 : din2 = 0d0
1146 1440148 : din3 = 0d0
1147 1440148 : din4 = 0d0
1148 1440148 : din5 = 0d0
1149 1440148 : dout1 = 0d0
1150 1440148 : dout2 = 0d0
1151 1440148 : dout3 = 0d0
1152 1440148 : dout4 = 0d0
1153 1440148 : dout5 = 0d0
1154 :
1155 : ! These rates are setup in update_special_rates
1156 80002 : select case(ir)
1157 :
1158 : case(ir_he4_he4_he4_to_c12) ! triple alpha
1159 80002 : if (g% use_3a_fl87) then
1160 : return
1161 : end if
1162 :
1163 : case(irn14ag_lite) ! n14 + 1.5 alpha => ne20
1164 1440148 : return
1165 :
1166 : end select
1167 :
1168 1360150 : num_reaction_inputs = get_num_reaction_inputs(ir)
1169 1360150 : num_reaction_outputs = get_num_reaction_outputs(ir)
1170 :
1171 : if (dbg) write(*,*) 'num_reaction_inputs', num_reaction_inputs
1172 : if (dbg) write(*,*) 'num_reaction_outputs', num_reaction_outputs
1173 : if (dbg) write(*,*)
1174 :
1175 1360150 : switch_to_prot = .false.
1176 1360150 : cout1 = 0; out1 = 0; o1 = 0
1177 1360150 : cout2 = 0; out2 = 0; o2 = 0
1178 1360150 : cout3 = 0; out3 = 0; o3 = 0
1179 1360150 : cin1 = 0; in1 = 0; i1 = 0
1180 1360150 : cin2 = 0; in2 = 0; i2 = 0
1181 1360150 : cin3 = 0; in3 = 0; i3 = 0
1182 1360150 : idr1 = 0; idr2 = 0; idr3 = 0
1183 1360150 : dr1 = 0; dr2 = 0; dr3 = 0
1184 :
1185 1360150 : if (num_reaction_outputs >= 1) then
1186 1360150 : cout1 = reaction_outputs(1,ir); dout1 = cout1
1187 1360150 : out1 = reaction_outputs(2,ir)
1188 1360150 : o1 = itab(out1)
1189 1360150 : if (o1 == 0) then
1190 0 : write(*,*) trim(reaction_Name(ir))
1191 0 : call mesa_error(__FILE__,__LINE__,'get1_derivs: itab(out1) = 0')
1192 : end if
1193 : end if
1194 :
1195 1360150 : if (num_reaction_outputs >= 2) then
1196 240064 : cout2 = reaction_outputs(3,ir); dout2 = cout2
1197 240064 : out2 = reaction_outputs(4,ir)
1198 240064 : o2 = itab(out2)
1199 240064 : if (o2 == 0) then
1200 0 : write(*,*) trim(reaction_Name(ir))
1201 0 : call mesa_error(__FILE__,__LINE__,'get1_derivs: itab(out2) = 0')
1202 : end if
1203 : end if
1204 :
1205 1360150 : if (num_reaction_outputs >= 3) then
1206 0 : cout3 = reaction_outputs(5,ir); dout3 = cout3
1207 0 : out3 = reaction_outputs(6,ir)
1208 0 : o3 = itab(out3)
1209 0 : if (o3 == 0) then
1210 0 : write(*,*) trim(reaction_Name(ir))
1211 0 : call mesa_error(__FILE__,__LINE__,'get1_derivs: itab(out3) = 0')
1212 : end if
1213 : end if
1214 :
1215 1360150 : if (num_reaction_outputs >= 4) then
1216 0 : write(*,*) trim(reaction_Name(ir))
1217 0 : call mesa_error(__FILE__,__LINE__,'get1_derivs: num_reaction_outputs >= 4')
1218 : end if
1219 :
1220 1360150 : if (num_reaction_inputs == 1) then
1221 320056 : cin1 = reaction_inputs(1,ir); din1 = cin1
1222 320056 : in1 = reaction_inputs(2,ir)
1223 320056 : i1 = itab(in1)
1224 1040094 : else if (num_reaction_inputs == 2 .or. num_reaction_inputs == 3) then
1225 1040094 : cin1 = reaction_inputs(1,ir); din1 = cin1
1226 1040094 : in1 = reaction_inputs(2,ir)
1227 1040094 : i1 = itab(in1)
1228 1040094 : cin2 = reaction_inputs(3,ir); din2 = cin2
1229 1040094 : in2 = reaction_inputs(4,ir)
1230 1040094 : i2 = itab(in2)
1231 1040094 : if (num_reaction_inputs == 3) then
1232 320004 : cin3 = reaction_inputs(5,ir); din3 = cin3
1233 320004 : in3 = reaction_inputs(6,ir)
1234 320004 : i3 = itab(in3)
1235 : end if
1236 : end if
1237 :
1238 1360150 : switch_to_prot = (prot /= 0) .and. (max_Z >= min_Z_for_switch_to_prot)
1239 1360150 : if (switch_to_prot) then
1240 0 : if (i1 == h1) then
1241 0 : in1 = iprot
1242 0 : i1 = prot
1243 0 : else if (i2 == h1) then
1244 0 : in2 = iprot
1245 0 : i2 = prot
1246 0 : else if (i3 == h1) then
1247 0 : in3 = iprot
1248 0 : i3 = prot
1249 : end if
1250 0 : if (o1 == h1) then
1251 1360150 : out1 = iprot
1252 1360150 : o1 = prot
1253 0 : else if (o2 == h1) then
1254 1360150 : out2 = iprot
1255 1360150 : o2 = prot
1256 0 : else if (o3 == h1) then
1257 0 : out3 = iprot
1258 0 : o3 = prot
1259 : end if
1260 : end if
1261 :
1262 1360150 : if (num_reaction_inputs == 1) then
1263 :
1264 320056 : if (i1 == 0) then
1265 0 : write(*,*) trim(reaction_Name(ir))
1266 0 : write(*,2) 'num_reaction_inputs', num_reaction_inputs
1267 0 : call mesa_error(__FILE__,__LINE__,'get1_derivs: itab(in1) = 0')
1268 : end if
1269 :
1270 320056 : if (cin1 == 1) then
1271 46 : r = n% y(i1)
1272 46 : idr1 = i1
1273 46 : dr1 = 1
1274 320010 : else if (cin1 == 3 .and. in1 /= ih1) then ! 3 he4
1275 : !write(*,'(/,a)') '1/6*r reaction name <' // trim(reaction_Name(ir)) // '>'
1276 80002 : r = (1d0/6d0)*n% y(i1)*n% y(i1)*n% y(i1)
1277 80002 : idr1 = i1
1278 80002 : dr1 = 0.5d0*n% y(i1)*n% y(i1)
1279 : else ! 2 body
1280 : !write(*,'(/,a)') '1/2*r reaction name <' // trim(reaction_Name(ir)) // '>'
1281 : !write(*,'(i3,3x,99e20.10)') i, n% rate_raw(i), n% rate_screened(i)
1282 240008 : r= 0.5d0*n% y(i1)*n% y(i1)
1283 240008 : idr1 = i1
1284 240008 : dr1 = n% y(i1)
1285 : !stop
1286 : end if
1287 :
1288 1040094 : else if (num_reaction_inputs == 2 .or. num_reaction_inputs == 3) then
1289 :
1290 1040094 : if (reaction_ye_rho_exponents(2,ir) == 0) then
1291 : ! treat as 1 body reaction
1292 0 : r = n% y(i1)
1293 0 : idr1 = i1
1294 0 : dr1 = 1
1295 0 : idr2 = i2
1296 0 : dr2 = 0
1297 : !write(*,*) 'get1_derivs rho=0: ' // trim(reaction_Name(ir))
1298 : !call mesa_error(__FILE__,__LINE__,'net_derivs')
1299 1040094 : else if ((cin1 == 1 .and. cin2 == 1) .or. reaction_ye_rho_exponents(2,ir) == 1) then
1300 : ! reaction_ye_rho_exponents(2,ir) == 1 for electron captures; treat as 2 body reaction
1301 1040094 : r = n% y(i1)*n% y(i2)
1302 1040094 : dr1 = n% y(i1)
1303 1040094 : idr1 = i2
1304 1040094 : dr2 = n% y(i2)
1305 1040094 : idr2 = i1
1306 0 : else if (cin1 == 2 .and. cin2 == 1) then
1307 0 : r = 0.5d0*n% y(i1)*n% y(i1)*n% y(i2)
1308 0 : dr1 = 0.5d0*n% y(i1)*n% y(i1)
1309 0 : idr1 = i2
1310 0 : dr2 = n% y(i1)*n% y(i2)
1311 0 : idr2 = i1
1312 0 : else if (cin1 == 1 .and. cin2 == 2) then
1313 : ! e.g., rhe4p, r_neut_he4_he4_to_be9, r_neut_h1_h1_to_h1_h2
1314 0 : r = n% y(i1)*0.5d0*n% y(i2)*n% y(i2)
1315 0 : dr1 = n% y(i1)*n% y(i2)
1316 0 : idr1 = i2
1317 0 : dr2 = 0.5d0*n% y(i2)*n% y(i2)
1318 0 : idr2 = i1
1319 0 : else if (cin1 == 2 .and. cin2 == 2) then
1320 : ! e.g., r_neut_neut_he4_he4_to_h3_li7, r_h1_h1_he4_he4_to_he3_be7
1321 0 : r = 0.5d0*n% y(i1)*n% y(i1)*0.5d0*n% y(i2)*n% y(i2)
1322 0 : dr1 = 0.5d0*n% y(i1)*n% y(i1)*n% y(i2)
1323 0 : idr1 = i2
1324 0 : dr2 = n% y(i1)*0.5d0*n% y(i2)*n% y(i2)
1325 0 : idr2 = i1
1326 : else
1327 0 : write(*,*) 'get1_derivs: ' // trim(reaction_Name(ir)) // ' invalid coefficient'
1328 0 : call mesa_error(__FILE__,__LINE__,'get1_derivs')
1329 : end if
1330 :
1331 1040094 : if (num_reaction_inputs == 3) then
1332 : ! we assume that the 3rd kind of input is just "along for the ride"
1333 : ! e.g., some compound reactions such as r34_pp2 are in this category.
1334 320004 : dr3 = 0
1335 320004 : idr3 = i3
1336 320004 : if (i3 == 0) then
1337 0 : write(*,*) trim(reaction_Name(ir))
1338 0 : call mesa_error(__FILE__,__LINE__,'get1_derivs: itab(in3) = 0')
1339 : end if
1340 : end if
1341 :
1342 : else
1343 :
1344 0 : write(*,*) 'get1_derivs: ' // trim(reaction_Name(ir)) // ' invalid specification'
1345 0 : call mesa_error(__FILE__,__LINE__,'get1_derivs')
1346 :
1347 : end if
1348 :
1349 :
1350 : ! after we've set these reaction values, pack them into arrays
1351 :
1352 1360150 : i_in(1) = i1; i_in(2) = i2; i_in(3) = i3
1353 1360150 : c_in(1) = din1; c_in(2) = din2; c_in(3) = din3
1354 :
1355 1360150 : i_out(1) = o1; i_out(2) = o2; i_out(3) = o3
1356 1360150 : c_out(1) = dout1; c_out(2) = dout2; c_out(3) = dout3
1357 :
1358 1360150 : dr(1) = dr1; dr(2) = dr2; dr(3) = dr3
1359 1360150 : idr(1) = idr1; idr(2) = idr2; idr(3) = idr3
1360 :
1361 :
1362 : ! for debugging
1363 :
1364 1360150 : if (num_reaction_inputs == 1) then
1365 :
1366 : if (num_reaction_outputs == 1) then
1367 : ! reaction of form din1 in1 -> dout1 out1
1368 : if (dbg) write(*,*) ' do_one_one din1', din1, trim(chem_isos% name(g% chem_id(i1)))
1369 : if (dbg) write(*,*) 'do_one_one dout1', dout1, trim(chem_isos% name(g% chem_id(o1)))
1370 :
1371 : else if (num_reaction_outputs == 2) then
1372 : ! reaction of form cin1 in1 -> dout1 out1 + dout2 out2
1373 : if (dbg) write(*,*) ' do_one_two din1', din1, trim(chem_isos% name(g% chem_id(i1)))
1374 : if (dbg) write(*,*) 'do_one_two dout1', dout1, trim(chem_isos% name(g% chem_id(o1)))
1375 : if (dbg) write(*,*) 'do_one_two dout2', dout2, trim(chem_isos% name(g% chem_id(o2)))
1376 :
1377 : if (.false. .and. reaction_Name(ir) == 'r_he4_he4_he4_to_h1_b11' .and. r > 0) then
1378 : write(*,'(3i6,3x,a,2x,99e20.10)') i, ir, &
1379 : reaction_ye_rho_exponents(2,ir), &
1380 : 'do_one_two ' // trim(reaction_Name(ir)) // ' ' // &
1381 : trim(chem_isos% name(g% chem_id(i1))) // ' => ' // &
1382 : trim(chem_isos% name(g% chem_id(o1))) // ' + ' // &
1383 : trim(chem_isos% name(g% chem_id(o2))), &
1384 : n% rate_screened(i), r, dr1, n% y(i1)
1385 : stop
1386 : end if
1387 :
1388 : else if (num_reaction_outputs == 3) then
1389 : ! reaction of form cin1 in1 -> dout1 out1 + dout2 out2 + dout3 out3
1390 : if (dbg) write(*,*) ' do_one_three din1', din1, trim(chem_isos% name(g% chem_id(i1)))
1391 : if (dbg) write(*,*) 'do_one_three dout1', dout1, trim(chem_isos% name(g% chem_id(o1)))
1392 : if (dbg) write(*,*) 'do_one_three dout2', dout2, trim(chem_isos% name(g% chem_id(o2)))
1393 : if (dbg) write(*,*) 'do_one_three dout3', dout3, trim(chem_isos% name(g% chem_id(o3)))
1394 :
1395 : else
1396 0 : write(*,*) trim(reaction_Name(ir))
1397 0 : write(*,*) 'too many reaction_outputs for num_reaction_inputs == 1'
1398 0 : call mesa_error(__FILE__, __LINE__)
1399 : end if
1400 :
1401 1040094 : else if (num_reaction_inputs == 2) then
1402 :
1403 : if (num_reaction_outputs == 1) then
1404 : ! reaction of form din1 in1 + din2 in2 -> dout1 out1
1405 : if (dbg) write(*,*) ' do_two_one din1', din1, trim(chem_isos% name(g% chem_id(i1)))
1406 : if (dbg) write(*,*) ' do_two_one din2', din2, trim(chem_isos% name(g% chem_id(i2)))
1407 : if (dbg) write(*,*) 'do_two_one dout1', dout1, trim(chem_isos% name(g% chem_id(o1)))
1408 :
1409 : if (.false. .and. reaction_Name(ir) == 'r_neut_he4_he4_to_be9' .and. r > 0 .and. &
1410 : abs(n% y(i1) - 7.7763751756339478D-05) < 1d-20) then
1411 : write(*,'(i3,3x,a,2x,99e20.10)') i, &
1412 : 'do_two_one ' // trim(reaction_Name(ir)) // ' ' // &
1413 : trim(chem_isos% name(g% chem_id(i1))) // ' + ' // &
1414 : trim(chem_isos% name(g% chem_id(i2))) // ' => ' // &
1415 : trim(chem_isos% name(g% chem_id(o1))), &
1416 : r, dr1, dr2, n% y(i1), n% y(i2)
1417 : !stop
1418 : end if
1419 :
1420 : if (.false. .and. reaction_Name(ir) == 'r_he4_si28_to_o16_o16') then
1421 : write(*,2) 'n% y(i1)', i1, n% y(i1)
1422 : write(*,2) 'n% y(i2)', i2, n% y(i2)
1423 : write(*,1) 'r', r
1424 : write(*,1) 'rate screened', n% rate_screened(i)
1425 : write(*,1) 'r*y1*y2', n% y(i1)*n% y(i2)*n% rate_screened(i)
1426 : !stop
1427 : end if
1428 :
1429 : else if (num_reaction_outputs == 2) then
1430 : ! reaction of form din1 in1 + din2 in2 -> dout1 out1 + dout2 out2
1431 : if (dbg) write(*,*) ' do_two_two din1', din1, trim(chem_isos% name(g% chem_id(i1)))
1432 : if (dbg) write(*,*) ' do_two_two din2', din2, trim(chem_isos% name(g% chem_id(i2)))
1433 : if (dbg) write(*,*) 'do_two_two dout1', dout1, trim(chem_isos% name(g% chem_id(o1)))
1434 : if (dbg) write(*,*) 'do_two_two dout2', dout2, trim(chem_isos% name(g% chem_id(o2)))
1435 :
1436 :
1437 : else if (num_reaction_outputs == 3) then
1438 : ! reaction of form din1 in1 + din2 in2 -> dout1 out1 + dout2 out2 + dout3 out3
1439 : if (dbg) write(*,*) ' do_two_three din1', din1, trim(chem_isos% name(g% chem_id(i1)))
1440 : if (dbg) write(*,*) ' do_two_three din2', din2, trim(chem_isos% name(g% chem_id(i2)))
1441 : if (dbg) write(*,*) 'do_two_three dout1', dout1, trim(chem_isos% name(g% chem_id(o1)))
1442 : if (dbg) write(*,*) 'do_two_three dout2', dout2, trim(chem_isos% name(g% chem_id(o2)))
1443 : if (dbg) write(*,*) 'do_two_three dout3', dout3, trim(chem_isos% name(g% chem_id(o3)))
1444 :
1445 : else
1446 0 : write(*,*) trim(reaction_Name(ir))
1447 0 : write(*,*) 'too many reaction_outputs for num_reaction_inputs == 2'
1448 0 : call mesa_error(__FILE__, __LINE__)
1449 : end if
1450 :
1451 320004 : else if (num_reaction_inputs == 3) then
1452 :
1453 320004 : if (num_reaction_outputs == 1) then
1454 : ! reaction of form din1 in1 + din2 in2 + din3 in3 -> dout1 out1
1455 : if (dbg) write(*,*) ' do_three_one din1', din1, trim(chem_isos% name(g% chem_id(i1)))
1456 : if (dbg) write(*,*) ' do_three_one din2', din2, trim(chem_isos% name(g% chem_id(i2)))
1457 : if (dbg) write(*,*) ' do_three_one din3', din3, trim(chem_isos% name(g% chem_id(i3)))
1458 : if (dbg) write(*,*) 'do_three_one dout1', dout1, trim(chem_isos% name(g% chem_id(o1)))
1459 :
1460 80000 : else if (num_reaction_outputs == 2) then
1461 : ! reaction of form din1 in1 + din2 in2 + din3 in3 -> dout1 out1 + dout2 out2
1462 : if (dbg) write(*,*) ' do_three_two din1', din1, trim(chem_isos% name(g% chem_id(i1)))
1463 : if (dbg) write(*,*) ' do_three_two din2', din2, trim(chem_isos% name(g% chem_id(i2)))
1464 : if (dbg) write(*,*) ' do_three_two din3', din3, trim(chem_isos% name(g% chem_id(i3)))
1465 : if (dbg) write(*,*) 'do_three_two dout1', dout1, trim(chem_isos% name(g% chem_id(o1)))
1466 : if (dbg) write(*,*) 'do_three_two dout2', dout2, trim(chem_isos% name(g% chem_id(o2)))
1467 :
1468 : else
1469 0 : write(*,*) trim(reaction_Name(ir))
1470 0 : write(*,*) 'too many reaction_outputs for num_reaction_inputs == 3'
1471 0 : call mesa_error(__FILE__, __LINE__)
1472 : end if
1473 :
1474 : else
1475 0 : write(*,*) 'too many reaction_inputs'
1476 0 : call mesa_error(__FILE__, __LINE__)
1477 : end if
1478 :
1479 :
1480 : ! all reactions are handled by do_in_out except for 1 -> reactions from weaklib
1481 : ! these must go through do_in_out_neu in order to use the weaklib Q and Qneu values
1482 :
1483 1360150 : if (num_reaction_inputs == 1 .and. num_reaction_outputs == 1) then
1484 :
1485 240026 : weak_id = g% weak_reaction_index(i)
1486 :
1487 240026 : done = .false.
1488 240026 : if (weak_id > 0) then
1489 160022 : if (weak_id > g% num_wk_reactions) then
1490 0 : write(*,2) 'i', i
1491 0 : write(*,2) 'ir', ir
1492 0 : write(*,2) 'weak_id', weak_id
1493 0 : write(*,2) 'g% num_wk_reactions', g% num_wk_reactions
1494 0 : write(*,*) trim(trim(reaction_Name(ir)))
1495 0 : call mesa_error(__FILE__,__LINE__,'derivs')
1496 : end if
1497 160022 : if (g% weaklib_ids(weak_id) > 0) then ! > 0 means included in weaklib
1498 :
1499 8 : n% rate_screened(i) = n% lambda(weak_id)
1500 8 : n% rate_screened_dT(i) = n % dlambda_dlnT(weak_id) / temp
1501 8 : n% rate_screened_dRho(i) = n% dlambda_dlnRho(weak_id) / den
1502 :
1503 : call do_in_out_neu( &
1504 : n, dydt, eps_nuc_MeV, i, r, &
1505 : num_reaction_inputs, i_in, c_in, &
1506 : num_reaction_outputs, i_out, c_out, &
1507 : idr, dr, n% Q(weak_id), &
1508 : n% Qneu(weak_id), n% dQneu_dlnT(weak_id)/temp, n% dQneu_dlnRho(weak_id)/den, &
1509 8 : deriv_flgs, symbolic, just_dydt)
1510 :
1511 : done = .true.
1512 : end if
1513 : end if
1514 : if (.not. done) then ! weak reaction not in weaklib
1515 :
1516 : if (.false. .and. trim(reaction_Name(ir)) == 'r_o15_wk_n15') then
1517 : write(*,'(A)')
1518 : write(*,'(2(i5,1x),a,2x,99e20.10)') i, ir, &
1519 : 'do_one_one_neu ' // trim(reaction_Name(ir)) // ' ' // &
1520 : trim(chem_isos% name(g% chem_id(i1))) // ' => ' // &
1521 : trim(chem_isos% name(g% chem_id(o1)))
1522 : call mesa_error(__FILE__,__LINE__,'weak reaction not in weaklib')
1523 : end if
1524 :
1525 : call do_in_out( &
1526 : n, dydt, eps_nuc_MeV, i, r, &
1527 : num_reaction_inputs, i_in, c_in, &
1528 : num_reaction_outputs, i_out, c_out, &
1529 : idr, dr, &
1530 240018 : deriv_flgs, symbolic, just_dydt)
1531 :
1532 : end if
1533 :
1534 : else ! all non 1->1 reactions
1535 :
1536 : call do_in_out( &
1537 : n, dydt, eps_nuc_MeV, i, r, &
1538 : num_reaction_inputs, i_in, c_in, &
1539 : num_reaction_outputs, i_out, c_out, &
1540 : idr, dr, &
1541 1120124 : deriv_flgs, symbolic, just_dydt)
1542 :
1543 : end if
1544 :
1545 2480160 : end subroutine get1_derivs
1546 :
1547 :
1548 9106 : subroutine eval_ni56_ec_rate( &
1549 : temp, den, ye, eta, zbar, weak_rate_factor, &
1550 : rate, d_rate_dlnT, d_rate_dlnRho, Q, Qneu, ierr)
1551 : real(dp), intent(in) :: temp, den, ye, eta, zbar, weak_rate_factor
1552 : real(dp), intent(out) :: rate, d_rate_dlnT, d_rate_dlnRho, Q, Qneu
1553 : integer, intent(out) :: ierr
1554 : real(dp) :: dQneu_dlnT, dQneu_dlnRho
1555 : include 'formats'
1556 : call eval1_weak_rate( &
1557 : weak_rate_id_for_ni56_ec, irni56ec_to_fe56, &
1558 : temp, ye, den, eta, zbar, &
1559 : weak_rate_factor, rate, d_rate_dlnT, d_rate_dlnRho, Q, &
1560 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
1561 9106 : ierr)
1562 9106 : end subroutine eval_ni56_ec_rate
1563 :
1564 :
1565 4 : subroutine eval_co56_ec_rate( &
1566 : temp, den, ye, eta, zbar, weak_rate_factor, &
1567 : rate, d_rate_dlnT, d_rate_dlnRho, Q, Qneu, ierr)
1568 : real(dp), intent(in) :: temp, den, ye, eta, zbar, weak_rate_factor
1569 : real(dp), intent(out) :: rate, d_rate_dlnT, d_rate_dlnRho, Q, Qneu
1570 : integer, intent(out) :: ierr
1571 : real(dp) :: dQneu_dlnT, dQneu_dlnRho
1572 : include 'formats'
1573 : call eval1_weak_rate( &
1574 : weak_rate_id_for_co56_ec, irco56ec_to_fe56, &
1575 : temp, ye, den, eta, zbar, &
1576 : weak_rate_factor, rate, d_rate_dlnT, d_rate_dlnRho, Q, &
1577 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
1578 4 : ierr)
1579 4 : end subroutine eval_co56_ec_rate
1580 :
1581 :
1582 9110 : subroutine eval1_weak_rate( &
1583 : id, ir, temp, ye, rho, eta, zbar, weak_rate_factor, &
1584 : rate_out, d_rate_dlnT, d_rate_dlnRho, Q_out, &
1585 : Qneu_out, dQneu_dlnT_out, dQneu_dlnRho_out, &
1586 : ierr)
1587 : use rates_def, only: Coulomb_Info
1588 : use rates_lib, only: eval_weak_reaction_info
1589 : integer, intent(in) :: id, ir
1590 : real(dp), intent(in) :: temp, ye, rho, eta, zbar, weak_rate_factor
1591 : real(dp), intent(out) :: &
1592 : rate_out, d_rate_dlnT, d_rate_dlnRho, Q_out, &
1593 : Qneu_out, dQneu_dlnT_out, dQneu_dlnRho_out
1594 : integer, intent(out) :: ierr
1595 :
1596 : integer :: ids(1), reaction_id_for_weak_reactions(1)
1597 : type(Coulomb_Info), pointer :: cc
1598 : type(Coulomb_Info), target :: cc_info
1599 18220 : real(dp) :: T9, YeRho, d_eta_dlnT, d_eta_dlnRho
1600 : ! lambda = combined rate (capture and decay)
1601 : ! Q and Qneu are for combined rate of beta decay and electron capture.
1602 : ! Q is total, so Q-Qneu is the actual thermal energy.
1603 : ! note: lambdas include Ye Rho factors for electron captures.
1604 : ! so treat the rates as if just beta decays
1605 : real(dp), dimension(1), target :: &
1606 45550 : lambda_a, dlambda_dlnT_a, dlambda_dlnRho_a, &
1607 45550 : Q_a, dQ_dlnT_a, dQ_dlnRho_a, &
1608 45550 : Qneu_a, dQneu_dlnT_a, dQneu_dlnRho_a
1609 : real(dp), dimension(:), pointer :: &
1610 9110 : lambda, dlambda_dlnT, dlambda_dlnRho, &
1611 9110 : Q, dQ_dlnT, dQ_dlnRho, &
1612 9110 : Qneu, dQneu_dlnT, dQneu_dlnRho
1613 :
1614 : include 'formats'
1615 :
1616 : ierr = 0
1617 9110 : if (id <= 0) then
1618 0 : ierr = -1
1619 0 : return
1620 : end if
1621 :
1622 9110 : lambda => lambda_a
1623 9110 : dlambda_dlnT => dlambda_dlnT_a
1624 9110 : dlambda_dlnRho => dlambda_dlnRho_a
1625 :
1626 9110 : Q => Q_a
1627 9110 : dQ_dlnT => dQ_dlnT_a
1628 9110 : dQ_dlnRho => dQ_dlnRho_a
1629 :
1630 9110 : Qneu => Qneu_a
1631 9110 : dQneu_dlnT => dQneu_dlnT_a
1632 9110 : dQneu_dlnRho => dQneu_dlnRho_a
1633 :
1634 9110 : ids(1) = id
1635 9110 : reaction_id_for_weak_reactions(1) = ir
1636 9110 : T9 = temp*1d-9
1637 9110 : YeRho = ye*rho
1638 9110 : d_eta_dlnT = 0
1639 9110 : d_eta_dlnRho = 0
1640 9110 : cc => cc_info
1641 : call eval_weak_reaction_info( &
1642 : 1, ids, reaction_id_for_weak_reactions, &
1643 : cc, T9, YeRho, &
1644 : eta, d_eta_dlnT, d_eta_dlnRho, &
1645 : lambda, dlambda_dlnT, dlambda_dlnRho, &
1646 : Q, dQ_dlnT, dQ_dlnRho, &
1647 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
1648 9110 : ierr)
1649 :
1650 9110 : if (ierr /= 0) then
1651 : return
1652 :
1653 : write(*,*) 'failed in eval_weak_reaction_info'
1654 : call mesa_error(__FILE__,__LINE__,'eval1_weak_rate')
1655 : end if
1656 :
1657 9110 : rate_out = lambda(1)*weak_rate_factor
1658 9110 : d_rate_dlnT = dlambda_dlnT(1)*weak_rate_factor
1659 9110 : d_rate_dlnRho = dlambda_dlnRho(1)*weak_rate_factor
1660 :
1661 :
1662 9110 : Q_out = Q(1)
1663 9110 : Qneu_out = Qneu(1)
1664 9110 : dQneu_dlnT_out = dQneu_dlnT(1)
1665 9110 : dQneu_dlnRho_out = dQneu_dlnRho(1)
1666 :
1667 18220 : end subroutine eval1_weak_rate
1668 :
1669 :
1670 2480160 : subroutine update_special_rates( &
1671 2480160 : n, dydt, eps_nuc_MeV, i, eta, ye, temp, den, abar, zbar, &
1672 2480160 : num_reactions, rate_factors, rtab, itab, &
1673 : deriv_flgs, symbolic, just_dydt, ierr)
1674 9110 : use rates_lib, only: eval_n14_electron_capture_rate
1675 : type (Net_Info) :: n
1676 : integer, intent(in) :: i, num_reactions
1677 : real(qp), intent(inout) :: dydt(:,:)
1678 : real(qp), intent(out) :: eps_nuc_MeV(num_rvs)
1679 : real(dp), intent(in) :: eta, ye, temp, den, abar, zbar, rate_factors(:)
1680 : integer, pointer, intent(in) :: rtab(:), itab(:)
1681 : logical, pointer :: deriv_flgs(:)
1682 : logical, intent(in) :: symbolic, just_dydt
1683 : integer, intent(out) :: ierr
1684 :
1685 : integer :: ir, j, prot, neut, h1, he4, c14, n14, ne20, ne22, &
1686 : mg21, mg22, mg23, mg24, al23, al24, si24, si25, si26, &
1687 : s28, s29, s30, cl31, ar32, ar33, ar34, k35, ca36, ca37, ca38, &
1688 : ti41, ti42, v43, cr44, cr45, cr46, cr56, &
1689 : fe48, fe49, fe50, fe51, fe52, fe54, fe56, fe58, fe60, fe62, fe64, &
1690 : ni52, ni53, ni54, ni55, ni56, ni58, ni60, ni62, &
1691 : zn57, zn58, zn59, zn60, ge62, ge63, ge64, &
1692 : se68, kr72, sr76, mo84, sn104
1693 : integer :: weak_id, num_reaction_inputs, in1, in2, in3, in4, in5
1694 : integer :: cin1, cin2, cin3, cin4, cin5
1695 9920640 : real(dp) :: din1, din2, din3, din4, din5
1696 : integer :: num_reaction_outputs, out1, out2, out3, out4, out5
1697 : integer :: cout1, cout2, cout3, cout4, cout5
1698 9920640 : real(dp) :: dout1, dout2, dout3, dout4, dout5
1699 : type (Net_General_Info), pointer :: g
1700 2480160 : integer, pointer :: reaction_id(:)
1701 : integer :: i1, i2, i3, idr1, idr2, idr3, o1, o2, o3
1702 2480160 : real(dp) :: r, dr1, dr2, dr3, rate, d_rate_dlnT, d_rate_dlnRho, Q, Qneu, rn14ec
1703 :
1704 : integer, dimension(3) :: i_in, i_out, idr
1705 24801600 : real(dp), dimension(3) :: c_in, c_out, dr
1706 :
1707 : logical :: done, has_prot, has_neut, has_h1, switch_to_prot
1708 : integer :: max_Z, Z_plus_N_for_max_Z
1709 : integer, parameter :: min_Z_for_switch_to_prot = 12
1710 :
1711 : logical, parameter :: dbg = .false.
1712 :
1713 : include 'formats'
1714 :
1715 2480160 : g => n% g
1716 2480160 : reaction_id => g% reaction_id
1717 :
1718 2480160 : ierr = 0
1719 2480160 : prot = itab(iprot)
1720 2480160 : neut = itab(ineut)
1721 2480160 : h1 = itab(ih1)
1722 2480160 : he4 = itab(ihe4)
1723 :
1724 2480160 : ir = reaction_id(i)
1725 :
1726 2560158 : if (reaction_outputs(1,ir) == 0) return ! skip aux reactions
1727 :
1728 : if (dbg) &
1729 : write(*,'(/,a,2i6)') ' reaction name <' // trim(reaction_Name(ir)) // '>', i, ir
1730 :
1731 1440148 : max_Z = g% reaction_max_Z(i)
1732 1440148 : Z_plus_N_for_max_Z = g% reaction_max_Z_plus_N_for_max_Z(i)
1733 :
1734 1440148 : din1 = 0d0
1735 1440148 : din2 = 0d0
1736 1440148 : din3 = 0d0
1737 1440148 : din4 = 0d0
1738 1440148 : din5 = 0d0
1739 1440148 : dout1 = 0d0
1740 1440148 : dout2 = 0d0
1741 1440148 : dout3 = 0d0
1742 1440148 : dout4 = 0d0
1743 1440148 : dout5 = 0d0
1744 :
1745 :
1746 2480160 : select case(ir)
1747 :
1748 : case(ir_he4_he4_he4_to_c12) ! triple alpha
1749 80002 : if (g% use_3a_fl87) then
1750 0 : call do_FL_3alf(i) ! Fushiki and Lamb, Apj, 317, 368-388, 1987
1751 0 : return
1752 : end if
1753 :
1754 :
1755 : case(irn14ag_lite) ! n14 + 1.5 alpha => ne20
1756 79998 : n14 = itab(in14)
1757 79998 : ne20 = itab(ine20)
1758 79998 : r = n% y(n14) * n% y(he4)
1759 79998 : dr1 = n% y(n14)
1760 79998 : dr2 = n% y(he4)
1761 :
1762 :
1763 79998 : i_in(1) = n14; i_in(2) = he4; i_in(3) = 0
1764 79998 : c_in(1) = 1d0; c_in(2) = 1.5d0; c_in(3) = 0d0
1765 :
1766 79998 : i_out(1) = ne20; i_out(2) = 0; i_out(3) = 0
1767 79998 : c_out(1) = 1d0; c_out(2) = 0d0; c_out(3) = 0d0
1768 :
1769 79998 : idr(1) = he4; idr(2) = n14; idr(3) = 0
1770 79998 : dr(1) = dr1; dr(2) = dr2; dr(3) = 0
1771 :
1772 : call do_in_out(n, dydt, eps_nuc_MeV, i, r, &
1773 : 2, i_in, c_in, &
1774 : 1, i_out, c_out, &
1775 : idr, dr, &
1776 79998 : deriv_flgs, symbolic, just_dydt)
1777 :
1778 1440148 : return
1779 :
1780 : end select
1781 :
1782 : contains
1783 :
1784 :
1785 0 : subroutine do_FL_3alf(i) ! Fushiki and Lamb, Apj, 317, 368-388, 1987
1786 : use rates_lib, only: eval_FL_epsnuc_3alf
1787 : integer, intent(in) :: i
1788 : integer :: he4, c12
1789 0 : real(dp) :: UE, XHe4, YHe4, &
1790 0 : FLeps_nuc, dFLeps_nuc_dT, dFLeps_nuc_dRho, r, drdT, drdRho, conv
1791 : include 'formats'
1792 0 : he4 = itab(ihe4)
1793 0 : c12 = itab(ic12)
1794 0 : UE = abar/zbar
1795 0 : YHe4 = n% y(he4)
1796 0 : XHe4 = 4d0*YHe4
1797 0 : if (YHe4 < 1d-50) then
1798 0 : n% rate_screened(i) = 0
1799 0 : n% rate_screened_dT(i) = 0
1800 0 : n% rate_screened_dRho(i) = 0
1801 0 : r = 0
1802 0 : dr1 = 0
1803 : else
1804 : call eval_FL_epsnuc_3alf( &
1805 0 : temp, den, XHe4, UE, FLeps_nuc, dFLeps_nuc_dT, dFLeps_nuc_dRho)
1806 0 : conv = Qconv*n% reaction_Qs(ir)
1807 0 : r = YHe4*YHe4*YHe4/6d0
1808 0 : dr1 = 0.5d0*YHe4*YHe4
1809 0 : n% rate_screened(i) = FLeps_nuc/r*rate_factors(i)/conv
1810 0 : n% rate_screened_dT(i) = dFLeps_nuc_dT/r*rate_factors(i)/conv
1811 0 : n% rate_screened_dRho(i) = dFLeps_nuc_dRho/r*rate_factors(i)/conv
1812 : end if
1813 :
1814 0 : i_in(1) = he4; i_in(2) = 0; i_in(3) = 0
1815 0 : c_in(1) = 3d0; c_in(2) = 0d0; c_in(3) = 0d0
1816 :
1817 0 : i_out(1) = c12; i_out(2) = 0; i_out(3) = 0
1818 0 : c_out(1) = 1d0; c_out(2) = 0d0; c_out(3) = 0d0
1819 :
1820 0 : idr(1) = he4; idr(2) = 0; idr(3) = 0
1821 0 : dr(1) = dr1; dr(2) = 0; dr(3) = 0
1822 :
1823 : call do_in_out(n, dydt, eps_nuc_MeV, i, r, &
1824 : 1, i_in, c_in, &
1825 : 1, i_out, c_out, &
1826 : idr, dr, &
1827 0 : deriv_flgs, symbolic, just_dydt)
1828 :
1829 0 : end subroutine do_FL_3alf
1830 :
1831 : end subroutine update_special_rates
1832 :
1833 : end module net_derivs
|
