Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2014 Frank Timmes & 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_approx21
21 : use const_def, only: dp, qp, avo, clight
22 : use utils_lib, only: is_bad, mesa_error
23 : use math_lib
24 :
25 :
26 : implicit none
27 :
28 :
29 :
30 : !logical :: plus_co56 ! Must now be passed as an argument
31 :
32 : logical, parameter :: reduced_net_for_testing = .true.
33 : !logical, parameter :: reduced_net_for_testing = .false.
34 :
35 : integer, parameter :: species_21 = 21, species_co56 = 22
36 :
37 :
38 : integer :: iso_cid(species_co56) ! these are corresponding chem ids for the isos
39 : ! e.g., iso_cid(ife52) is = the iso number for fe52 as defined in mesa/chem
40 : ! Define as largest possible array
41 : integer :: & ! these are indices in y vector
42 : ih1, &
43 : ihe3, &
44 : ihe4, &
45 : ic12, &
46 : in14, &
47 : io16, &
48 : ine20, &
49 : img24, &
50 : isi28, &
51 : is32, &
52 : iar36, &
53 : ica40, &
54 : iti44, &
55 : icr48, &
56 : icrx, &
57 : ife52, &
58 : ife53, &
59 : ife54, &
60 : ife55, &
61 : ife56, &
62 : ico56, &
63 : ini56, &
64 : ineut, &
65 : iprot
66 :
67 : integer, parameter :: approx21_num_mesa_reactions_21 = 93, approx21_nrat = 116
68 : integer, parameter :: approx21_num_mesa_reactions_co56 = approx21_num_mesa_reactions_21+1, &
69 : approx21_plus_co56_nrat = approx21_nrat+1
70 :
71 : ! integer :: num_mesa_reactions
72 : ! integer :: num_reactions
73 :
74 : integer :: rate_id(approx21_num_mesa_reactions_co56) ! rate ids for the mesa reactions
75 : ! e.g., rate_id(ir3a) is reaction id for triple alpha as defined in mesa/rates
76 : ! Define as largest possible array
77 : integer :: & ! these are indices in rates arrays
78 : ir3a, &
79 : irg3a, &
80 : ircag, &
81 : ir1212, &
82 : ir1216, &
83 : ir1616, &
84 : iroga, &
85 : iroag, &
86 : irnega, &
87 : irneag, &
88 : irmgga, &
89 : irmgag, &
90 : irsiga, &
91 : irmgap, &
92 : iralpa, &
93 : iralpg, &
94 : irsigp, &
95 : irsiag, &
96 : irsga, &
97 : irsiap, &
98 : irppa, &
99 : irppg, &
100 : irsgp, &
101 : irsag, &
102 : irarga, &
103 : irsap, &
104 : irclpa, &
105 : irclpg, &
106 : irargp, &
107 : irarag, &
108 : ircaga, &
109 : irarap, &
110 : irkpa, &
111 : irkpg, &
112 : ircagp, &
113 : ircaag, &
114 : irtiga, &
115 : ircaap, &
116 : irscpa, &
117 : irscpg, &
118 : irtigp, &
119 : irtiag, &
120 : ircrga, &
121 : irtiap, &
122 : irvpa , &
123 : irvpg, &
124 : ircrgp, &
125 : ircrag, &
126 : irfega, &
127 : ircrap, &
128 : irmnpa, &
129 : irmnpg, &
130 : irfegp, &
131 : irfeag, &
132 : irniga, &
133 : irfeap, &
134 : ircopa, &
135 : ircopg, &
136 : irnigp, &
137 :
138 : ! for fe54 photodisintegration
139 : ir52ng, &
140 : ir53gn, &
141 : ir53ng, &
142 : ir54gn, &
143 : irfepg, &
144 : ircogp, &
145 :
146 : ! for he4 photodisintegration
147 : irheng, &
148 : irhegn, &
149 : irhng, &
150 : irdgn, &
151 : irdpg, &
152 : irhegp, &
153 :
154 : ! weak reactions
155 : irpen, &
156 : irnep, &
157 : irn56ec, &
158 : irco56ec, &
159 :
160 : ! ppchain
161 : irpp, &
162 : ir33, &
163 : irhe3ag, &
164 : ir_be7_wk_li7, &
165 : ir_be7_pg_b8, &
166 :
167 : ! cno cycles
168 : ircpg, &
169 : irnpg, &
170 : iropg, &
171 : irnag, &
172 :
173 : ! for reactions to fe56
174 : ir54ng, &
175 : ir55gn, &
176 : ir55ng, &
177 : ir56gn, &
178 : irfe54ap, &
179 : irco57pa, &
180 : irfe56pg, &
181 : irco57gp, &
182 :
183 : ! for n15 branching
184 : irn15pa, &
185 : irn15pg, &
186 :
187 : ! the equilibrium links
188 : ifa, &
189 : ifg, &
190 :
191 : irr1, &
192 : irs1, &
193 : irt1, &
194 : iru1, &
195 : irv1, &
196 : irw1, &
197 : irx1, &
198 :
199 : ir1f54, &
200 : ir2f54, &
201 : ir3f54, &
202 : ir4f54, &
203 : ir5f54, &
204 : ir6f54, &
205 : ir7f54, &
206 : ir8f54, &
207 :
208 : iralf1, &
209 : iralf2, &
210 :
211 : irfe56_aux1, &
212 : irfe56_aux2, &
213 : irfe56_aux3, &
214 : irfe56_aux4
215 :
216 : ! names
217 : character (len=40) :: ratnam(approx21_plus_co56_nrat)
218 : ! Define as largest possible array
219 :
220 : contains
221 :
222 :
223 :
224 : ! call this after get raw rates
225 9106 : subroutine approx21_pa_pg_fractions( &
226 9106 : ratraw,dratrawdt,dratrawdd,ierr)
227 : real(dp), dimension(:) :: ratraw,dratrawdt,dratrawdd
228 : integer, intent(out) :: ierr
229 :
230 : include 'formats'
231 :
232 9106 : ierr = 0
233 :
234 9106 : call set1(ifa,irn15pg,irn15pa)
235 9106 : ratraw(ifg) = 1.0d0 - ratraw(ifa)
236 9106 : dratrawdt(ifg) = -dratrawdt(ifa)
237 9106 : dratrawdd(ifg) = -dratrawdd(ifa)
238 :
239 9106 : call set1(irr1,iralpg,iralpa) ! al27
240 9106 : call set1(irs1,irppg,irppa) ! p31
241 9106 : call set1(irt1,irclpg,irclpa) ! cl35
242 9106 : call set1(iru1,irkpg,irkpa) ! k39
243 9106 : call set1(irv1,irscpg,irscpa) ! sc43
244 9106 : call set1(irw1,irvpg,irvpa) ! v47
245 9106 : call set1(irx1,irmnpg,irmnpa) ! mn51
246 :
247 :
248 : contains
249 :
250 72848 : subroutine set1(ifa,irn15pg,irn15pa)
251 : integer, intent(in) :: ifa,irn15pg,irn15pa
252 72848 : real(dp) :: ff1,dff1dt,dff1dd,ff2,dff2dt,dff2dd, &
253 72848 : tot,dtotdt,dtotdd,invtot
254 :
255 72848 : ff1 = ratraw(irn15pg)
256 72848 : dff1dt = dratrawdt(irn15pg)
257 72848 : dff1dd = dratrawdd(irn15pg)
258 :
259 72848 : ff2 = ratraw(irn15pa)
260 72848 : dff2dt = dratrawdt(irn15pa)
261 72848 : dff2dd = dratrawdd(irn15pa)
262 :
263 72848 : tot = ff1 + ff2
264 72848 : dtotdt = dff1dt + dff2dt
265 72848 : dtotdd = dff1dd + dff2dd
266 :
267 72848 : if (tot > 1d-30) then
268 72832 : invtot = 1.0d0/tot
269 72832 : ratraw(ifa) = ff2 * invtot
270 72832 : dratrawdt(ifa) = dff2dt * invtot - ff2 * invtot*invtot * dtotdt
271 72832 : dratrawdd(ifa) = dff2dd * invtot - ff2 * invtot*invtot * dtotdd
272 : else
273 16 : ratraw(ifa) = 0.0d0
274 16 : dratrawdt(ifa) = 0.0d0
275 16 : dratrawdd(ifa) = 0.0d0
276 : end if
277 :
278 72848 : end subroutine set1
279 :
280 : end subroutine approx21_pa_pg_fractions
281 :
282 :
283 : ! call this before screening
284 9106 : subroutine approx21_weak_rates( &
285 9106 : y, ratraw, dratrawdt, dratrawdd, &
286 : temp, den, ye, eta, zbar, &
287 : weak_rate_factor, plus_co56, ierr)
288 : use rates_lib, only: eval_ecapnuc_rate
289 : use net_derivs, only: eval_ni56_ec_rate, eval_co56_ec_rate
290 :
291 : real(dp), dimension(:) :: y, ratraw, dratrawdt, dratrawdd
292 : real(dp), intent(in) :: temp, den, ye, eta, zbar, weak_rate_factor
293 : logical, intent(in) :: plus_co56
294 : integer, intent(out) :: ierr
295 :
296 : real(dp) :: rpen, rnep, spen, snep, &
297 : rate, d_rate_dlnT, d_rate_dlnRho, Q, Qneu
298 : include 'formats'
299 :
300 9106 : ierr = 0
301 :
302 9106 : call eval_ecapnuc_rate(eta, temp, den, rpen, rnep, spen, snep)
303 :
304 9106 : ratraw(irpen) = rpen
305 9106 : dratrawdt(irpen) = 0
306 9106 : dratrawdd(irpen) = 0
307 : if (rpen > 0) then
308 : Qneu = spen/rpen
309 : else
310 : Qneu = 0
311 : end if
312 :
313 9106 : ratraw(irnep) = rnep
314 9106 : dratrawdt(irnep) = 0
315 9106 : dratrawdd(irnep) = 0
316 : if (rnep > 0) then
317 : Qneu = snep/rnep
318 : else
319 : Qneu = 0
320 : end if
321 :
322 : call eval_ni56_ec_rate( &
323 : temp, den, ye, eta, zbar, weak_rate_factor, &
324 : rate, d_rate_dlnT, d_rate_dlnRho, Q, Qneu, &
325 9106 : ierr)
326 9106 : if (ierr /= 0) then
327 : !write(*,*) 'failed in eval_ni56_ec_rate'
328 : return
329 : end if
330 9106 : ratraw(irn56ec) = rate
331 9106 : dratrawdt(irn56ec) = 0
332 9106 : dratrawdd(irn56ec) = 0
333 :
334 9106 : if (plus_co56) then
335 : call eval_co56_ec_rate( &
336 : temp, den, ye, eta, zbar, weak_rate_factor, &
337 : rate, d_rate_dlnT, d_rate_dlnRho, Q, Qneu, &
338 4 : ierr)
339 4 : if (ierr /= 0) then
340 : !write(*,*) 'failed in eval_co56_ec_rate'
341 : return
342 : end if
343 4 : ratraw(irco56ec) = rate
344 4 : dratrawdt(irco56ec) = 0
345 4 : dratrawdd(irco56ec) = 0
346 : end if
347 :
348 9106 : end subroutine approx21_weak_rates
349 :
350 :
351 : ! call this after screening -- depends on y, so don't reuse results.
352 9106 : subroutine approx21_special_reactions( &
353 9106 : btemp, bden, abar, zbar, y, &
354 : use_3a_FL, conv_eps_3a, &
355 9106 : ratdum, dratdumdt, dratdumdd, dratdumdy1, dratdumdy2, &
356 : plus_co56, ierr)
357 9106 : use math_lib
358 : use utils_lib, only: mesa_error
359 : use rates_lib, only: eval_FL_epsnuc_3alf
360 : real(dp), intent(in) :: &
361 : btemp, bden, abar, zbar, y(:), conv_eps_3a
362 : logical, intent(in) :: use_3a_fl
363 : real(dp), dimension(:) :: &
364 : ratdum,dratdumdt,dratdumdd,dratdumdy1,dratdumdy2
365 : logical, intent(in) :: plus_co56
366 : integer, intent(out) :: ierr
367 :
368 9106 : real(dp) :: denom, denomdt, denomdd, zz, xx, eps, deps_dT, deps_dRho
369 : real(dp), parameter :: tiny_denom = 1d-50, tiny_y = 1d-30
370 : integer :: i
371 : logical :: okay
372 : include 'formats'
373 :
374 9106 : ierr = 0
375 :
376 9106 : if (use_3a_FL) then
377 : ! Fushiki and Lamb, Apj, 317, 368-388, 1987
378 0 : if (y(ihe4) < tiny_y) then
379 0 : ratdum(ir3a) = 0.0d0
380 0 : dratdumdt(ir3a) = 0.0d0
381 0 : dratdumdd(ir3a) = 0.0d0
382 : else
383 : call eval_FL_epsnuc_3alf( &
384 0 : btemp, bden, 4*y(ihe4), abar/zbar, eps, deps_dT, deps_dRho)
385 : ! convert from eps back to rate
386 0 : xx = conv_eps_3a*y(ihe4)*y(ihe4)*y(ihe4)/6d0
387 0 : ratdum(ir3a) = eps/xx
388 0 : dratdumdt(ir3a) = deps_dT/xx
389 0 : dratdumdd(ir3a) = deps_dRho/xx
390 : end if
391 : end if
392 :
393 :
394 9106 : okay = .true.
395 855968 : do i=1,num_mesa_reactions(plus_co56)
396 855968 : if (ratdum(i) < 0d0) then
397 0 : write(*,2) 'approx21 missing rate for ' // ratnam(i), i, ratdum(i), &
398 0 : btemp, log10(btemp), bden, log10(bden)
399 0 : okay = .false.
400 : end if
401 : end do
402 9106 : if (.not. okay) call mesa_error(__FILE__,__LINE__)
403 :
404 : ! for debugging: sum(cat)/eps_nuc
405 :
406 : if (reduced_net_for_testing) then
407 :
408 : !if (.true.) then
409 :
410 : !end if
411 :
412 : if (.false.) then
413 :
414 : ! turn off PP
415 : call turn_off_reaction(irpp)
416 : call turn_off_reaction(ir33)
417 : call turn_off_reaction(irhe3ag)
418 :
419 : !turn off CNO
420 : call turn_off_reaction(ircpg)
421 : call turn_off_reaction(irnpg)
422 : call turn_off_reaction(iropg)
423 :
424 : ! turn off n14ag + 3alpha
425 : call turn_off_reaction(irnag)
426 : call turn_off_reaction(ir3a)
427 : call turn_off_reaction(irg3a)
428 :
429 : !c12+c12, c12ag, o16ag
430 : call turn_off_reaction(ir1212)
431 : call turn_off_reaction(iroag)
432 : call turn_off_reaction(irnega)
433 : call turn_off_reaction(ircag)
434 : call turn_off_reaction(iroga)
435 :
436 : !Ne/O burn
437 : call turn_off_reaction(ir1216)
438 : call turn_off_reaction(ir1616)
439 : call turn_off_reaction(irneag)
440 : call turn_off_reaction(irmgga)
441 :
442 : !alpha links
443 : call turn_off_reaction(irmgag)
444 : call turn_off_reaction(irsiga)
445 : call turn_off_reaction(irmgap)
446 : call turn_off_reaction(iralpa)
447 : call turn_off_reaction(iralpg)
448 : call turn_off_reaction(irsigp)
449 : call turn_off_reaction(irsiag)
450 : call turn_off_reaction(irsga)
451 :
452 : call turn_off_reaction(irppa)
453 : call turn_off_reaction(irppg)
454 : call turn_off_reaction(irsiap)
455 : call turn_off_reaction(irsgp)
456 :
457 : call turn_off_reaction(irsag)
458 : call turn_off_reaction(irarga)
459 : call turn_off_reaction(irsap)
460 : call turn_off_reaction(irclpa)
461 : call turn_off_reaction(irclpg)
462 : call turn_off_reaction(irargp)
463 :
464 : call turn_off_reaction(irarag)
465 : call turn_off_reaction(ircaga)
466 : call turn_off_reaction(irarap)
467 : call turn_off_reaction(irkpa)
468 : call turn_off_reaction(irkpg)
469 : call turn_off_reaction(ircagp)
470 :
471 : call turn_off_reaction(ircaag)
472 : call turn_off_reaction(irtiga)
473 : call turn_off_reaction(ircaap)
474 : call turn_off_reaction(irscpa)
475 : call turn_off_reaction(irscpg)
476 : call turn_off_reaction(irtigp)
477 : call turn_off_reaction(irtiag)
478 :
479 : call turn_off_reaction(ircrga)
480 : call turn_off_reaction(irtiap)
481 : call turn_off_reaction(irvpa )
482 : call turn_off_reaction(irvpg)
483 : call turn_off_reaction(ircrgp)
484 : call turn_off_reaction(ircrag)
485 : call turn_off_reaction(irfega)
486 : call turn_off_reaction(ircrap)
487 : call turn_off_reaction(irmnpa)
488 : call turn_off_reaction(irmnpg)
489 : call turn_off_reaction(irfegp)
490 : call turn_off_reaction(irfeag)
491 :
492 : !iron group
493 : call turn_off_reaction(irniga)
494 : call turn_off_reaction(irfeap)
495 : call turn_off_reaction(ircopa)
496 :
497 : call turn_off_reaction(irnigp)
498 : call turn_off_reaction(irfepg)
499 : call turn_off_reaction(ircogp)
500 :
501 : call turn_off_reaction(irheng)
502 : call turn_off_reaction(irhegn)
503 :
504 : call turn_off_reaction(irhng)
505 : call turn_off_reaction(irdgn)
506 :
507 : call turn_off_reaction(irdpg)
508 : call turn_off_reaction(irhegp)
509 :
510 : call turn_off_reaction(irpen)
511 : call turn_off_reaction(irnep)
512 :
513 : call turn_off_reaction(ircopg)
514 :
515 : call turn_off_reaction(ir54ng)
516 : call turn_off_reaction(ir55gn)
517 : call turn_off_reaction(ir55ng)
518 : call turn_off_reaction(ir56gn)
519 :
520 : call turn_off_reaction(ir52ng)
521 : call turn_off_reaction(ir53gn)
522 : call turn_off_reaction(ir53ng)
523 : call turn_off_reaction(ir54gn)
524 :
525 : call turn_off_reaction(irfe56pg)
526 : call turn_off_reaction(irco57gp)
527 :
528 : call turn_off_reaction(irfe54ap)
529 : call turn_off_reaction(irco57pa)
530 :
531 : call turn_off_reaction(irco56ec)
532 : call turn_off_reaction(irn56ec)
533 :
534 : end if
535 :
536 : !if (.true.) then
537 :
538 :
539 : !end if
540 :
541 : end if
542 :
543 : ! fe52(n,g)fe53(n,g)fe54 equilibrium links
544 9106 : ratdum(ir1f54) = 0.0d0
545 9106 : dratdumdy1(ir1f54) = 0.0d0
546 9106 : dratdumdt(ir1f54) = 0.0d0
547 9106 : dratdumdd(ir1f54) = 0.0d0
548 :
549 9106 : ratdum(ir2f54) = 0.0d0
550 9106 : dratdumdy1(ir2f54) = 0.0d0
551 9106 : dratdumdt(ir2f54) = 0.0d0
552 9106 : dratdumdd(ir2f54) = 0.0d0
553 :
554 9106 : denom = ratdum(ir53gn) + y(ineut)*ratdum(ir53ng)
555 9106 : denomdt = dratdumdt(ir53gn) + y(ineut)*dratdumdt(ir53ng)
556 9106 : denomdd = dratdumdd(ir53gn) + y(ineut)*dratdumdd(ir53ng)
557 :
558 9106 : if (denom > tiny_denom .and. btemp > 1.5d9) then
559 5934 : zz = 1.0d0/denom
560 :
561 5934 : ratdum(ir1f54) = ratdum(ir54gn)*ratdum(ir53gn)*zz
562 5934 : dratdumdy1(ir1f54) = -ratdum(ir1f54)*zz * ratdum(ir53ng)
563 : dratdumdt(ir1f54) = dratdumdt(ir54gn)*ratdum(ir53gn)*zz &
564 : + ratdum(ir54gn)*dratdumdt(ir53gn)*zz &
565 5934 : - ratdum(ir1f54)*zz*denomdt
566 : dratdumdd(ir1f54) = dratdumdd(ir54gn)*ratdum(ir53gn)*zz &
567 : + ratdum(ir54gn)*dratdumdd(ir53gn)*zz &
568 5934 : - ratdum(ir1f54)*zz*denomdd
569 :
570 5934 : ratdum(ir2f54) = ratdum(ir52ng)*ratdum(ir53ng)*zz
571 5934 : dratdumdy1(ir2f54) = -ratdum(ir2f54)*zz * ratdum(ir53ng)
572 : dratdumdt(ir2f54) = dratdumdt(ir52ng)*ratdum(ir53ng)*zz &
573 : + ratdum(ir52ng)*dratdumdt(ir53ng)*zz &
574 5934 : - ratdum(ir2f54)*zz*denomdt
575 : dratdumdd(ir2f54) = dratdumdd(ir52ng)*ratdum(ir53ng)*zz &
576 : + ratdum(ir52ng)*dratdumdd(ir53ng)*zz &
577 5934 : - ratdum(ir2f54)*zz*denomdd
578 : end if
579 :
580 : ! fe54(n,g)fe55(n,g)fe56 equilibrium links
581 9106 : ratdum(irfe56_aux1) = 0.0d0
582 9106 : dratdumdy1(irfe56_aux1) = 0.0d0
583 9106 : dratdumdt(irfe56_aux1) = 0.0d0
584 9106 : dratdumdd(irfe56_aux1) = 0.0d0
585 :
586 9106 : ratdum(irfe56_aux2) = 0.0d0
587 9106 : dratdumdy1(irfe56_aux2) = 0.0d0
588 9106 : dratdumdt(irfe56_aux2) = 0.0d0
589 9106 : dratdumdd(irfe56_aux2) = 0.0d0
590 :
591 9106 : denom = ratdum(ir55gn) + y(ineut)*ratdum(ir55ng)
592 9106 : denomdt = dratdumdt(ir55gn) + y(ineut)*dratdumdt(ir55ng)
593 9106 : denomdd = dratdumdd(ir55gn) + y(ineut)*dratdumdd(ir55ng)
594 :
595 9106 : if (denom > tiny_denom .and. btemp > 1.5d9) then
596 5934 : zz = 1.0d0/denom
597 :
598 5934 : ratdum(irfe56_aux1) = ratdum(ir56gn)*ratdum(ir55gn)*zz
599 5934 : dratdumdy1(irfe56_aux1) = -ratdum(irfe56_aux1)*zz * ratdum(ir55ng)
600 : dratdumdt(irfe56_aux1) = dratdumdt(ir56gn)*ratdum(ir55gn)*zz &
601 : + ratdum(ir56gn)*dratdumdt(ir55gn)*zz &
602 5934 : - ratdum(irfe56_aux1)*zz*denomdt
603 : dratdumdd(irfe56_aux1) = dratdumdd(ir56gn)*ratdum(ir55gn)*zz &
604 : + ratdum(ir56gn)*dratdumdd(ir55gn)*zz &
605 5934 : - ratdum(irfe56_aux1)*zz*denomdd
606 :
607 5934 : ratdum(irfe56_aux2) = ratdum(ir54ng)*ratdum(ir55ng)*zz
608 5934 : dratdumdy1(irfe56_aux2) = -ratdum(irfe56_aux2)*zz * ratdum(ir55ng)
609 : dratdumdt(irfe56_aux2) = dratdumdt(ir54ng)*ratdum(ir55ng)*zz &
610 : + ratdum(ir54ng)*dratdumdt(ir55ng)*zz &
611 5934 : - ratdum(irfe56_aux2)*zz*denomdt
612 : dratdumdd(irfe56_aux2) = dratdumdd(ir54ng)*ratdum(ir55ng)*zz &
613 : + ratdum(ir54ng)*dratdumdd(ir55ng)*zz &
614 5934 : - ratdum(irfe56_aux2)*zz*denomdd
615 :
616 : end if
617 :
618 : ! fe54(a,p)co57(g,p)fe56 equilibrium links
619 :
620 9106 : ratdum(irfe56_aux3) = 0.0d0
621 9106 : dratdumdy1(irfe56_aux3) = 0.0d0
622 9106 : dratdumdt(irfe56_aux3) = 0.0d0
623 9106 : dratdumdd(irfe56_aux3) = 0.0d0
624 :
625 9106 : ratdum(irfe56_aux4) = 0.0d0
626 9106 : dratdumdy1(irfe56_aux4) = 0.0d0
627 9106 : dratdumdt(irfe56_aux4) = 0.0d0
628 9106 : dratdumdd(irfe56_aux4) = 0.0d0
629 :
630 9106 : denom = ratdum(irco57gp) + y(iprot)*ratdum(irco57pa)
631 9106 : denomdt = dratdumdt(irco57gp) + y(iprot)*dratdumdt(irco57pa)
632 9106 : denomdd = dratdumdd(irco57gp) + y(iprot)*dratdumdd(irco57pa)
633 :
634 9106 : if (denom > tiny_denom .and. btemp > 1.5d9) then
635 5934 : zz = 1.0d0/denom
636 :
637 5934 : ratdum(irfe56_aux3) = ratdum(irfe56pg) * ratdum(irco57pa) * zz
638 5934 : dratdumdy1(irfe56_aux3) = -ratdum(irfe56_aux3) * zz * ratdum(irco57pa)
639 : dratdumdt(irfe56_aux3) = dratdumdt(irfe56pg) * ratdum(irco57pa) * zz &
640 : + ratdum(irfe56pg) * dratdumdt(irco57pa) * zz &
641 5934 : - ratdum(irfe56_aux3) * zz * denomdt
642 : dratdumdd(irfe56_aux3) = dratdumdd(irfe56pg) * ratdum(irco57pa) * zz &
643 : + ratdum(irfe56pg) * dratdumdd(irco57pa) * zz &
644 5934 : - ratdum(irfe56_aux3) * zz * denomdd
645 :
646 5934 : ratdum(irfe56_aux4) = ratdum(irfe54ap) * ratdum(irco57gp) * zz
647 5934 : dratdumdy1(irfe56_aux4) = -ratdum(irfe56_aux4) * zz * ratdum(irco57pa)
648 : dratdumdt(irfe56_aux4) = dratdumdt(irfe54ap) * ratdum(irco57gp) * zz &
649 : + ratdum(irfe54ap) * dratdumdt(irco57gp) * zz &
650 5934 : - ratdum(irfe56_aux4) * zz * denomdt
651 : dratdumdd(irfe56_aux4) = dratdumdd(irfe54ap) * ratdum(irco57gp) * zz &
652 : + ratdum(irfe54ap) * dratdumdd(irco57gp) * zz &
653 5934 : - ratdum(irfe56_aux4) * zz * denomdd
654 : end if
655 :
656 :
657 : ! fe54(p,g)co55(p,g)ni56 equilibrium links r3f54 r4f54
658 : ! fe52(a,p)co55(g,p)fe54 equilibrium links r5f54 r6f54
659 : ! fe52(a,p)co55(p,g)ni56 equilibrium links r7f54 r8f54
660 :
661 9106 : ratdum(ir3f54) = 0.0d0
662 9106 : dratdumdy1(ir3f54) = 0.0d0
663 9106 : dratdumdt(ir3f54) = 0.0d0
664 9106 : dratdumdd(ir3f54) = 0.0d0
665 :
666 9106 : ratdum(ir4f54) = 0.0d0
667 9106 : dratdumdy1(ir4f54) = 0.0d0
668 9106 : dratdumdt(ir4f54) = 0.0d0
669 9106 : dratdumdd(ir4f54) = 0.0d0
670 :
671 9106 : ratdum(ir5f54) = 0.0d0
672 9106 : dratdumdy1(ir5f54) = 0.0d0
673 9106 : dratdumdt(ir5f54) = 0.0d0
674 9106 : dratdumdd(ir5f54) = 0.0d0
675 :
676 9106 : ratdum(ir6f54) = 0.0d0
677 9106 : dratdumdy1(ir6f54) = 0.0d0
678 9106 : dratdumdt(ir6f54) = 0.0d0
679 9106 : dratdumdd(ir6f54) = 0.0d0
680 :
681 9106 : ratdum(ir7f54) = 0.0d0
682 9106 : dratdumdy1(ir7f54) = 0.0d0
683 9106 : dratdumdt(ir7f54) = 0.0d0
684 9106 : dratdumdd(ir7f54) = 0.0d0
685 :
686 9106 : ratdum(ir8f54) = 0.0d0
687 9106 : dratdumdy1(ir8f54) = 0.0d0
688 9106 : dratdumdt(ir8f54) = 0.0d0
689 9106 : dratdumdd(ir8f54) = 0.0d0
690 :
691 9106 : denom = ratdum(ircogp)+y(iprot)*(ratdum(ircopg)+ratdum(ircopa))
692 :
693 9106 : if (denom > tiny_denom .and. btemp > 1.5d9) then
694 :
695 : denomdt = dratdumdt(ircogp) &
696 5934 : + y(iprot)*(dratdumdt(ircopg) + dratdumdt(ircopa))
697 : denomdd = dratdumdd(ircogp) &
698 5934 : + y(iprot)*(dratdumdd(ircopg) + dratdumdd(ircopa))
699 :
700 5934 : zz = 1.0d0/denom
701 :
702 5934 : ratdum(ir3f54) = ratdum(irfepg) * ratdum(ircopg) * zz
703 : dratdumdy1(ir3f54) = -ratdum(ir3f54) * zz * &
704 5934 : (ratdum(ircopg) + ratdum(ircopa))
705 : dratdumdt(ir3f54) = dratdumdt(irfepg) * ratdum(ircopg) * zz &
706 : + ratdum(irfepg) * dratdumdt(ircopg) * zz &
707 5934 : - ratdum(ir3f54)*zz*denomdt
708 : dratdumdd(ir3f54) = dratdumdd(irfepg) * ratdum(ircopg) * zz &
709 : + ratdum(irfepg) * dratdumdd(ircopg) * zz &
710 5934 : - ratdum(ir3f54)*zz*denomdd
711 :
712 5934 : ratdum(ir4f54) = ratdum(irnigp) * ratdum(ircogp) * zz
713 : dratdumdy1(ir4f54) = -ratdum(ir4f54) * zz * &
714 5934 : (ratdum(ircopg)+ratdum(ircopa))
715 : dratdumdt(ir4f54) = dratdumdt(irnigp) * ratdum(ircogp) * zz &
716 : + ratdum(irnigp) * dratdumdt(ircogp) * zz &
717 5934 : - ratdum(ir4f54)*zz*denomdt
718 : dratdumdd(ir4f54) = dratdumdd(irnigp) * ratdum(ircogp) * zz &
719 : + ratdum(irnigp) * dratdumdd(ircogp) * zz &
720 5934 : - ratdum(ir4f54)*zz*denomdd
721 :
722 5934 : ratdum(ir5f54) = ratdum(irfepg) * ratdum(ircopa) * zz
723 : dratdumdy1(ir5f54) = -ratdum(ir5f54) * zz * &
724 5934 : (ratdum(ircopg)+ratdum(ircopa))
725 : dratdumdt(ir5f54) = dratdumdt(irfepg) * ratdum(ircopa) * zz &
726 : + ratdum(irfepg) * dratdumdt(ircopa) * zz &
727 5934 : - ratdum(ir5f54) * zz * denomdt
728 : dratdumdd(ir5f54) = dratdumdd(irfepg) * ratdum(ircopa) * zz &
729 : + ratdum(irfepg) * dratdumdd(ircopa) * zz &
730 5934 : - ratdum(ir5f54) * zz * denomdd
731 :
732 5934 : ratdum(ir6f54) = ratdum(irfeap) * ratdum(ircogp) * zz
733 : dratdumdy1(ir6f54) = -ratdum(ir6f54) * zz * &
734 5934 : (ratdum(ircopg)+ratdum(ircopa))
735 : dratdumdt(ir6f54) = dratdumdt(irfeap) * ratdum(ircogp) * zz &
736 : + ratdum(irfeap) * dratdumdt(ircogp) * zz &
737 5934 : - ratdum(ir6f54) * zz * denomdt
738 : dratdumdd(ir6f54) = dratdumdd(irfeap) * ratdum(ircogp) * zz &
739 : + ratdum(irfeap) * dratdumdd(ircogp) * zz &
740 5934 : - ratdum(ir6f54) * zz * denomdd
741 :
742 5934 : ratdum(ir7f54) = ratdum(irfeap) * ratdum(ircopg) * zz
743 :
744 : dratdumdy1(ir7f54) = -ratdum(ir7f54) * zz * &
745 5934 : (ratdum(ircopg)+ratdum(ircopa))
746 : dratdumdt(ir7f54) = dratdumdt(irfeap) * ratdum(ircopg) * zz &
747 : + ratdum(irfeap) * dratdumdt(ircopg) * zz &
748 5934 : - ratdum(ir7f54) * zz * denomdt
749 : dratdumdd(ir7f54) = dratdumdd(irfeap) * ratdum(ircopg) * zz &
750 : + ratdum(irfeap) * dratdumdd(ircopg) * zz &
751 5934 : - ratdum(ir7f54) * zz * denomdd
752 :
753 5934 : ratdum(ir8f54) = ratdum(irnigp) * ratdum(ircopa) * zz
754 :
755 : dratdumdy1(ir8f54) = -ratdum(ir8f54) * zz * &
756 5934 : (ratdum(ircopg)+ratdum(ircopa))
757 : dratdumdt(ir8f54) = dratdumdt(irnigp) * ratdum(ircopa) * zz &
758 : + ratdum(irnigp) * dratdumdt(ircopa) * zz &
759 5934 : - ratdum(ir8f54) * zz * denomdt
760 : dratdumdd(ir8f54) = dratdumdd(irnigp) * ratdum(ircopa) * zz &
761 : + ratdum(irnigp) * dratdumdd(ircopa) * zz &
762 5934 : - ratdum(ir8f54) * zz * denomdd
763 :
764 :
765 : end if
766 :
767 :
768 : ! p(n,g)h2(n,g)3h(p,g)he4 photodisintegrated n and p back to he4 equilibrium links
769 : ! p(n,g)h2(p,g)he3(n,g)he4
770 :
771 9106 : ratdum(iralf1) = 0.0d0
772 9106 : dratdumdy1(iralf1) = 0.0d0
773 9106 : dratdumdy2(iralf1) = 0.0d0
774 9106 : dratdumdt(iralf1) = 0.0d0
775 9106 : dratdumdd(iralf1) = 0.0d0
776 :
777 9106 : ratdum(iralf2) = 0.0d0
778 9106 : dratdumdy1(iralf2) = 0.0d0
779 9106 : dratdumdy2(iralf2) = 0.0d0
780 9106 : dratdumdt(iralf2) = 0.0d0
781 9106 : dratdumdd(iralf2) = 0.0d0
782 :
783 : denom = ratdum(irhegp)*ratdum(irdgn) + &
784 : y(ineut)*ratdum(irheng)*ratdum(irdgn) + &
785 9106 : y(ineut)*y(iprot)*ratdum(irheng)*ratdum(irdpg)
786 :
787 9106 : if (denom > tiny_denom .and. btemp > 1.5d9) then
788 :
789 : denomdt = dratdumdt(irhegp)*ratdum(irdgn) &
790 : + ratdum(irhegp)*dratdumdt(irdgn) &
791 : + y(ineut) * (dratdumdt(irheng)*ratdum(irdgn) &
792 : + ratdum(irheng)*dratdumdt(irdgn)) &
793 : + y(ineut)*y(iprot) * (dratdumdt(irheng)*ratdum(irdpg) &
794 5934 : + ratdum(irheng)*dratdumdt(irdpg))
795 :
796 : denomdd = dratdumdd(irhegp)*ratdum(irdgn) &
797 : + ratdum(irhegp)*dratdumdd(irdgn) &
798 : + y(ineut) * (dratdumdd(irheng)*ratdum(irdgn) &
799 : + ratdum(irheng)*dratdumdd(irdgn)) &
800 : + y(ineut)*y(iprot) * (dratdumdd(irheng)*ratdum(irdpg) &
801 5934 : + ratdum(irheng)*dratdumdd(irdpg))
802 :
803 5934 : zz = 1.0d0/denom
804 :
805 : ratdum(iralf1) = ratdum(irhegn) * ratdum(irhegp)* &
806 5934 : ratdum(irdgn) * zz
807 : dratdumdy1(iralf1) = -ratdum(iralf1) * zz * &
808 : (ratdum(irheng)*ratdum(irdgn) + &
809 5934 : y(iprot)*ratdum(irheng)*ratdum(irdpg))
810 : dratdumdy2(iralf1) = -ratdum(iralf1) * zz * y(ineut) * &
811 5934 : ratdum(irheng) * ratdum(irdpg)
812 : dratdumdt(iralf1) = dratdumdt(irhegn)*ratdum(irhegp)* &
813 : ratdum(irdgn) * zz &
814 : + ratdum(irhegn)*dratdumdt(irhegp)*ratdum(irdgn)*zz &
815 : + ratdum(irhegn)*ratdum(irhegp)*dratdumdt(irdgn)*zz &
816 5934 : - ratdum(iralf1)*zz*denomdt
817 : dratdumdd(iralf1) = dratdumdd(irhegn) * ratdum(irhegp)* &
818 : ratdum(irdgn) * zz &
819 : + ratdum(irhegn)*dratdumdd(irhegp)*ratdum(irdgn)*zz &
820 : + ratdum(irhegn)*ratdum(irhegp)*dratdumdd(irdgn)*zz &
821 5934 : - ratdum(iralf1)*zz*denomdd
822 :
823 :
824 : ratdum(iralf2) = ratdum(irheng)*ratdum(irdpg)* &
825 5934 : ratdum(irhng)*zz
826 : dratdumdy1(iralf2) = -ratdum(iralf2) * zz * &
827 : (ratdum(irheng)*ratdum(irdgn) + &
828 5934 : y(iprot)*ratdum(irheng)*ratdum(irdpg))
829 :
830 :
831 : denom = ratdum(irhegp)*ratdum(irdgn) + &
832 : y(ineut)*ratdum(irheng)*ratdum(irdgn) + &
833 5934 : y(ineut)*y(iprot)*ratdum(irheng)*ratdum(irdpg)
834 :
835 5934 : if (is_bad(dratdumdy1(iralf2))) then
836 0 : write(*,1) 'denom', denom
837 0 : write(*,1) 'zz', zz
838 0 : write(*,1) 'dratdumdy1(iralf2)', dratdumdy1(iralf2)
839 0 : write(*,1) 'ratdum(irhegp)*ratdum(irdgn)', ratdum(irhegp)*ratdum(irdgn)
840 0 : write(*,1) 'y(ineut)*ratdum(irheng)*ratdum(irdgn)', y(ineut)*ratdum(irheng)*ratdum(irdgn)
841 0 : write(*,1) 'y(ineut)*y(iprot)*ratdum(irheng)*ratdum(irdpg)', &
842 0 : y(ineut)*y(iprot)*ratdum(irheng)*ratdum(irdpg)
843 0 : write(*,1) 'ratdum(irhegp)', ratdum(irhegp)
844 0 : write(*,1) 'ratdum(irdgn)', ratdum(irdgn)
845 0 : write(*,1) 'ratdum(irheng)', ratdum(irheng)
846 0 : write(*,1) 'ratdum(irdgn)', ratdum(irdgn)
847 0 : write(*,1) 'y(ineut)', y(ineut)
848 0 : write(*,1) 'y(iprot)', y(iprot)
849 0 : stop
850 : end if
851 :
852 :
853 : dratdumdy2(iralf2) = -ratdum(iralf2) * zz * y(ineut)* &
854 5934 : ratdum(irheng) * ratdum(irdpg)
855 : dratdumdt(iralf2) = dratdumdt(irheng)*ratdum(irdpg) * &
856 : ratdum(irhng) * zz &
857 : + ratdum(irheng)*dratdumdt(irdpg)*ratdum(irhng)*zz &
858 : + ratdum(irheng)*ratdum(irdpg)*dratdumdt(irhng)*zz &
859 5934 : - ratdum(iralf2)*zz*denomdt
860 : dratdumdd(iralf2) = dratdumdd(irheng)*ratdum(irdpg)* &
861 : ratdum(irhng)*zz &
862 : + ratdum(irheng)*dratdumdd(irdpg)*ratdum(irhng)*zz &
863 : + ratdum(irheng)*ratdum(irdpg)*dratdumdd(irhng)*zz &
864 5934 : - ratdum(iralf2)*zz*denomdd
865 :
866 : end if
867 :
868 :
869 :
870 : ! he3(a,g)be7(p,g)8b(e+nu)8be(2a)
871 : ! beta limit he3+he4 by the 8B decay half life
872 9106 : if (y(ihe4) > tiny_y) then
873 9102 : xx = 0.896d0/y(ihe4)
874 9102 : ratdum(irhe3ag) = min(ratdum(irhe3ag),xx)
875 9102 : if (ratdum(irhe3ag) == xx) then
876 9100 : dratdumdy1(irhe3ag) = -xx/y(ihe4)
877 9100 : dratdumdt(irhe3ag) = 0.0d0
878 9100 : dratdumdd(irhe3ag) = 0.0d0
879 : else
880 2 : dratdumdy1(irhe3ag) = 0.0d0
881 : end if
882 : end if
883 :
884 :
885 : ! beta limit n14(p,g)o15(enu)o16 and o16(p,g)f17(e+nu)17o(p,a)n14
886 11726 : if (y(ih1) > tiny_y) then
887 :
888 2620 : xx = 5.68d-03/(y(ih1)*1.57d0)
889 2620 : ratdum(irnpg) = min(ratdum(irnpg),xx)
890 2620 : if (ratdum(irnpg) == xx) then
891 0 : dratdumdy1(irnpg) = -xx/y(ih1)
892 0 : dratdumdt(irnpg) = 0.0d0
893 0 : dratdumdd(irnpg) = 0.0d0
894 : else
895 2620 : dratdumdy1(irnpg) = 0.0d0
896 : end if
897 :
898 2620 : xx = 0.0105d0/y(ih1)
899 2620 : ratdum(iropg) = min(ratdum(iropg),xx)
900 2620 : if (ratdum(iropg) == xx) then
901 0 : dratdumdy1(iropg) = -xx/y(ih1)
902 0 : dratdumdt(iropg) = 0.0d0
903 0 : dratdumdd(iropg) = 0.0d0
904 : else
905 2620 : dratdumdy1(iropg) = 0.0d0
906 : end if
907 :
908 : end if
909 :
910 :
911 : contains
912 :
913 : subroutine turn_off_reaction(i)
914 : integer, intent(in) :: i
915 : if (i == 0) return
916 : ratdum(i) = 0
917 : dratdumdt(i) = 0
918 : dratdumdd(i) = 0
919 : dratdumdy1(i) = 0
920 : dratdumdy2(i) = 0
921 9106 : end subroutine turn_off_reaction
922 :
923 : end subroutine approx21_special_reactions
924 :
925 :
926 21194 : subroutine approx21_dydt( &
927 21194 : y, rate, ratdum, dydt, deriva, &
928 : fe56ec_fake_factor_in, min_T, fe56ec_n_neut, temp, den, plus_co56, ierr)
929 : logical, intent(in) :: deriva ! false for dydt, true for partials wrt T, Rho
930 : real(dp), dimension(:), intent(in) :: y, rate, ratdum
931 : integer, intent(in) :: fe56ec_n_neut
932 : real(dp), dimension(:), intent(out) :: dydt
933 : real(dp), intent(in) :: fe56ec_fake_factor_in, temp, den
934 21194 : real(dp) :: fe56ec_fake_factor, min_T
935 : logical, intent(in) :: plus_co56
936 : integer, intent(out) :: ierr
937 :
938 : integer :: i
939 :
940 : ! quad precision dydt sums
941 21194 : real(qp) :: a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,&
942 21194 : a11,a12,a13,a14,a15,a16,a17,a18,a19,a20
943 466280 : real(qp) :: qray(species(plus_co56))
944 :
945 : logical :: okay
946 :
947 : include 'formats'
948 :
949 21194 : ierr = 0
950 :
951 : ! Turn on special fe56ec rate above some temperature
952 21194 : fe56ec_fake_factor = 0d0
953 21194 : if(.not.deriva) then
954 21194 : fe56ec_fake_factor = eval_fe56ec_fake_factor(fe56ec_fake_factor_in, min_T, temp)
955 : end if
956 :
957 466280 : dydt(1:species(plus_co56)) = 0.0d0
958 466280 : qray(1:species(plus_co56)) = 0.0_qp
959 :
960 : ! hydrogen reactions
961 21194 : a1 = -1.5d0 * y(ih1) * y(ih1) * rate(irpp)
962 21194 : a2 = y(ihe3) * y(ihe3) * rate(ir33)
963 21194 : a3 = -y(ihe3) * y(ihe4) * rate(irhe3ag)
964 21194 : a4 = -2.0d0 * y(ic12) * y(ih1) * rate(ircpg)
965 21194 : a5 = -2.0d0 * y(in14) * y(ih1) * rate(irnpg)
966 21194 : a6 = -2.0d0 * y(io16) * y(ih1) * rate(iropg)
967 21194 : a7 = -3.0d0 * y(ih1) * rate(irpen)
968 :
969 21194 : qray(ih1) = qray(ih1) + a1 + a2 + a3 + a4 + a5 + a6 + a7
970 :
971 : ! he3 reactions
972 :
973 21194 : a1 = 0.5d0 * y(ih1) * y(ih1) * rate(irpp)
974 21194 : a2 = -y(ihe3) * y(ihe3) * rate(ir33)
975 21194 : a3 = -y(ihe3) * y(ihe4) * rate(irhe3ag)
976 21194 : a4 = y(ih1) * rate(irpen)
977 :
978 21194 : qray(ihe3) = qray(ihe3) + a1 + a2 + a3 + a4
979 :
980 :
981 : ! he4 reactions
982 : ! heavy ion reactions
983 21194 : a1 = 0.5d0 * y(ic12) * y(ic12) * rate(ir1212)
984 21194 : a2 = 0.5d0 * y(ic12) * y(io16) * rate(ir1216)
985 21194 : a3 = 0.56d0 * 0.5d0 * y(io16) * y(io16) * rate(ir1616)
986 21194 : a4 = -y(ihe4) * y(in14) * rate(irnag) * 1.5d0 ! n14 + 1.5 alpha => ne20
987 21194 : qray(ihe4) = qray(ihe4) + a1 + a2 + a3 + a4
988 :
989 :
990 : ! (a,g) and (g,a) reactions
991 :
992 21194 : a1 = -0.5d0 * y(ihe4) * y(ihe4) * y(ihe4) * rate(ir3a)
993 21194 : a2 = 3.0d0 * y(ic12) * rate(irg3a)
994 21194 : a3 = -y(ihe4) * y(ic12) * rate(ircag)
995 21194 : a4 = y(io16) * rate(iroga)
996 21194 : a5 = -y(ihe4) * y(io16) * rate(iroag)
997 21194 : a6 = y(ine20) * rate(irnega)
998 21194 : a7 = -y(ihe4) * y(ine20) * rate(irneag)
999 21194 : a8 = y(img24) * rate(irmgga)
1000 21194 : a9 = -y(ihe4) * y(img24)* rate(irmgag)
1001 21194 : a10 = y(isi28) * rate(irsiga)
1002 21194 : a11 = -y(ihe4) * y(isi28)*rate(irsiag)
1003 21194 : a12 = y(is32) * rate(irsga)
1004 :
1005 : qray(ihe4) = qray(ihe4) + &
1006 21194 : a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + a10 + a11 + a12
1007 :
1008 21194 : a1 = -y(ihe4) * y(is32) * rate(irsag)
1009 21194 : a2 = y(iar36) * rate(irarga)
1010 21194 : a3 = -y(ihe4) * y(iar36)*rate(irarag)
1011 21194 : a4 = y(ica40) * rate(ircaga)
1012 21194 : a5 = -y(ihe4) * y(ica40)*rate(ircaag)
1013 21194 : a6 = y(iti44) * rate(irtiga)
1014 21194 : a7 = -y(ihe4) * y(iti44)*rate(irtiag)
1015 21194 : a8 = y(icr48) * rate(ircrga)
1016 21194 : a9 = -y(ihe4) * y(icr48)*rate(ircrag)
1017 21194 : a10 = y(ife52) * rate(irfega)
1018 21194 : a11 = -y(ihe4) * y(ife52) * rate(irfeag)
1019 21194 : a12 = y(ini56) * rate(irniga)
1020 :
1021 : qray(ihe4) = qray(ihe4) + &
1022 21194 : a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + a10 + a11 + a12
1023 :
1024 :
1025 : ! (a,p)(p,g) and (g,p)(p,a) reactions
1026 :
1027 21194 : if (.not.deriva) then
1028 9106 : a1 = 0.34d0*0.5d0*y(io16)*y(io16)*rate(irs1)*rate(ir1616)
1029 9106 : a2 = -y(ihe4) * y(img24) * rate(irmgap)*(1.0d0-rate(irr1))
1030 9106 : a3 = y(isi28) * rate(irsigp) * rate(irr1)
1031 9106 : a4 = -y(ihe4) * y(isi28) * rate(irsiap)*(1.0d0-rate(irs1))
1032 9106 : a5 = y(is32) * rate(irsgp) * rate(irs1)
1033 9106 : a6 = -y(ihe4) * y(is32) * rate(irsap)*(1.0d0-rate(irt1))
1034 9106 : a7 = y(iar36) * rate(irargp) * rate(irt1)
1035 9106 : a8 = -y(ihe4) * y(iar36) * rate(irarap)*(1.0d0-rate(iru1))
1036 9106 : a9 = y(ica40) * rate(ircagp) * rate(iru1)
1037 9106 : a10 = -y(ihe4) * y(ica40) * rate(ircaap)*(1.0d0-rate(irv1))
1038 9106 : a11 = y(iti44) * rate(irtigp) * rate(irv1)
1039 9106 : a12 = -y(ihe4) * y(iti44) * rate(irtiap)*(1.0d0-rate(irw1))
1040 9106 : a13 = y(icr48) * rate(ircrgp) * rate(irw1)
1041 9106 : a14 = -y(ihe4) * y(icr48) * rate(ircrap)*(1.0d0-rate(irx1))
1042 9106 : a15 = y(ife52) * rate(irfegp) * rate(irx1)
1043 :
1044 : qray(ihe4) = qray(ihe4) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + &
1045 9106 : a8 + a9 + a10 + a11 + a12 + a13 + a14 + a15
1046 :
1047 : else
1048 12088 : a1 = 0.34d0*0.5d0*y(io16)*y(io16) * ratdum(irs1)*rate(ir1616)
1049 12088 : a2 = 0.34d0*0.5d0*y(io16)*y(io16) * rate(irs1) * ratdum(ir1616)
1050 12088 : a3 = -y(ihe4)*y(img24) * rate(irmgap)*(1.0d0 - ratdum(irr1))
1051 12088 : a4 = y(ihe4)*y(img24) * ratdum(irmgap)*rate(irr1)
1052 12088 : a5 = y(isi28) * ratdum(irsigp) * rate(irr1)
1053 12088 : a6 = y(isi28) * rate(irsigp) * ratdum(irr1)
1054 12088 : a7 = -y(ihe4)*y(isi28) * rate(irsiap)*(1.0d0 - ratdum(irs1))
1055 12088 : a8 = y(ihe4)*y(isi28) * ratdum(irsiap) * rate(irs1)
1056 12088 : a9 = y(is32) * ratdum(irsgp) * rate(irs1)
1057 12088 : a10 = y(is32) * rate(irsgp) * ratdum(irs1)
1058 :
1059 12088 : qray(ihe4) = qray(ihe4) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + a10
1060 :
1061 12088 : a1 = -y(ihe4)*y(is32) * rate(irsap)*(1.0d0 - ratdum(irt1))
1062 12088 : a2 = y(ihe4)*y(is32) * ratdum(irsap)*rate(irt1)
1063 12088 : a3 = y(iar36) * ratdum(irargp) * rate(irt1)
1064 12088 : a4 = y(iar36) * rate(irargp) * ratdum(irt1)
1065 12088 : a5 = -y(ihe4)*y(iar36) * rate(irarap)*(1.0d0 - ratdum(iru1))
1066 12088 : a6 = y(ihe4)*y(iar36) * ratdum(irarap)*rate(iru1)
1067 12088 : a7 = y(ica40) * ratdum(ircagp) * rate(iru1)
1068 12088 : a8 = y(ica40) * rate(ircagp) * ratdum(iru1)
1069 12088 : a9 = -y(ihe4)*y(ica40) * rate(ircaap)*(1.0d0-ratdum (irv1))
1070 12088 : a10 = y(ihe4)*y(ica40) * ratdum(ircaap)*rate(irv1)
1071 12088 : a11 = y(iti44) * ratdum(irtigp) * rate(irv1)
1072 12088 : a12 = y(iti44) * rate(irtigp) * ratdum(irv1)
1073 :
1074 : qray(ihe4) = qray(ihe4) + &
1075 12088 : a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + a10 + a11 + a12
1076 :
1077 12088 : a1 = -y(ihe4)*y(iti44) * rate(irtiap)*(1.0d0 - ratdum(irw1))
1078 12088 : a2 = y(ihe4)*y(iti44) * ratdum(irtiap)*rate(irw1)
1079 12088 : a3 = y(icr48) * ratdum(ircrgp) * rate(irw1)
1080 12088 : a4 = y(icr48) * rate(ircrgp) * ratdum(irw1)
1081 12088 : a5 = -y(ihe4)*y(icr48) * rate(ircrap)*(1.0d0 - ratdum(irx1))
1082 12088 : a6 = y(ihe4)*y(icr48) * ratdum(ircrap)*rate(irx1)
1083 12088 : a7 = y(ife52) * ratdum(irfegp) * rate(irx1)
1084 12088 : a8 = y(ife52) * rate(irfegp) * ratdum(irx1)
1085 :
1086 12088 : qray(ihe4) = qray(ihe4) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8
1087 : end if
1088 :
1089 :
1090 : ! photodisintegration reactions
1091 21194 : a1 = y(ife54) * y(iprot) * y(iprot) * rate(ir5f54)
1092 21194 : a2 = -y(ife52) * y(ihe4) * rate(ir6f54)
1093 21194 : a3 = -y(ife52) * y(ihe4) * y(iprot) * rate(ir7f54)
1094 21194 : a4 = y(ini56) * y(iprot) * rate(ir8f54)
1095 21194 : a5 = -y(ihe4) * rate(iralf1)
1096 21194 : a6 = y(ineut)*y(ineut) * y(iprot)*y(iprot) * rate(iralf2)
1097 21194 : a7 = y(ife56) * y(iprot) * y(iprot) * rate(irfe56_aux3)
1098 21194 : a8 = -y(ife54) * y(ihe4) * rate(irfe56_aux4)
1099 :
1100 21194 : qray(ihe4) = qray(ihe4) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8
1101 :
1102 :
1103 : ! ppchain
1104 21194 : a1 = 0.5d0 * y(ihe3) * y(ihe3) * rate(ir33)
1105 21194 : a2 = y(ihe3) * y(ihe4) * rate(irhe3ag)
1106 :
1107 21194 : qray(ihe4) = qray(ihe4) + a1 + a2
1108 :
1109 :
1110 : ! cno cycles
1111 21194 : a1 = y(io16) * y(ih1) * rate(iropg)
1112 :
1113 21194 : qray(ihe4) = qray(ihe4) + a1 + a2
1114 :
1115 21194 : if (.not. deriva) then
1116 9106 : a1 = y(in14) * y(ih1) * rate(ifa) * rate(irnpg)
1117 9106 : qray(ihe4) = qray(ihe4) + a1
1118 : else
1119 12088 : a1 = y(in14) * y(ih1) * rate(ifa) * ratdum(irnpg)
1120 12088 : a2 = y(in14) * y(ih1) * ratdum(ifa) * rate(irnpg)
1121 12088 : qray(ihe4) = qray(ihe4) + a1 + a2
1122 : end if
1123 :
1124 :
1125 : ! c12 reactions
1126 21194 : a1 = -y(ic12) * y(ic12) * rate(ir1212)
1127 21194 : a2 = -y(ic12) * y(io16) * rate(ir1216)
1128 21194 : a3 = (1d0/6d0) * y(ihe4) * y(ihe4) * y(ihe4) * rate(ir3a)
1129 21194 : a4 = -y(ic12) * rate(irg3a)
1130 21194 : a5 = -y(ic12) * y(ihe4) * rate(ircag)
1131 21194 : a6 = y(io16) * rate(iroga)
1132 21194 : a7 = -y(ic12) * y(ih1) * rate(ircpg)
1133 :
1134 21194 : qray(ic12) = qray(ic12) + a1 + a2 + a3 + a4 + a5 + a6 + a7
1135 :
1136 21194 : if (.not. deriva) then
1137 9106 : a1 = y(in14) * y(ih1) * rate(ifa) * rate(irnpg)
1138 9106 : qray(ic12) = qray(ic12) + a1
1139 : else
1140 12088 : a1 = y(in14) * y(ih1) * rate(ifa) * ratdum(irnpg)
1141 12088 : a2 = y(in14) * y(ih1) * ratdum(ifa) * rate(irnpg)
1142 12088 : qray(ic12) = qray(ic12) + a1 + a2
1143 : end if
1144 :
1145 :
1146 : ! n14 reactions
1147 21194 : a1 = y(ic12) * y(ih1) * rate(ircpg)
1148 21194 : a2 = -y(in14) * y(ih1) * rate(irnpg)
1149 21194 : a3 = y(io16) * y(ih1) * rate(iropg)
1150 21194 : a4 = -y(ihe4) * y(in14) * rate(irnag) ! n14 + 1.5 alpha => ne20
1151 :
1152 21194 : qray(in14) = qray(in14) + a1 + a2 + a3 + a4
1153 :
1154 :
1155 : ! o16 reactions
1156 21194 : a1 = -y(ic12) * y(io16) * rate(ir1216)
1157 21194 : a2 = -y(io16) * y(io16) * rate(ir1616)
1158 21194 : a3 = y(ic12) * y(ihe4) * rate(ircag)
1159 21194 : a4 = -y(io16) * y(ihe4) * rate(iroag)
1160 21194 : a5 = -y(io16) * rate(iroga)
1161 21194 : a6 = y(ine20) * rate(irnega)
1162 21194 : a7 = -y(io16) * y(ih1) * rate(iropg)
1163 :
1164 21194 : qray(io16) = qray(io16) + a1 + a2 + a3 + a4 + a5 + a6 + a7
1165 :
1166 21194 : if (.not. deriva) then
1167 9106 : a1 = y(in14) * y(ih1) * rate(ifg) * rate(irnpg)
1168 9106 : qray(io16) = qray(io16) + a1
1169 : else
1170 12088 : a1 = y(in14) * y(ih1) * rate(ifg) * ratdum(irnpg)
1171 12088 : a2 = y(in14) * y(ih1) * ratdum(ifg) * rate(irnpg)
1172 12088 : qray(io16) = qray(io16) + a1 + a2
1173 : end if
1174 :
1175 :
1176 : ! ne20 reactions
1177 21194 : a1 = 0.5d0 * y(ic12) * y(ic12) * rate(ir1212)
1178 21194 : a2 = y(io16) * y(ihe4) * rate(iroag)
1179 21194 : a3 = -y(ine20) * y(ihe4) * rate(irneag)
1180 21194 : a4 = -y(ine20) * rate(irnega)
1181 21194 : a5 = y(img24) * rate(irmgga)
1182 21194 : a6 = y(in14) * y(ihe4) * rate(irnag) ! n14 + 1.5 alpha => ne20
1183 :
1184 21194 : qray(ine20) = qray(ine20) + a1 + a2 + a3 + a4 + a5 + a6
1185 :
1186 :
1187 : ! mg24 reactions
1188 21194 : a1 = 0.5d0 * y(ic12) * y(io16) * rate(ir1216)
1189 21194 : a2 = y(ine20) * y(ihe4) * rate(irneag)
1190 21194 : a3 = -y(img24) * y(ihe4) * rate(irmgag)
1191 21194 : a4 = -y(img24) * rate(irmgga)
1192 21194 : a5 = y(isi28) * rate(irsiga)
1193 :
1194 21194 : qray(img24) = qray(img24) + a1 + a2 + a3 + a4 + a5
1195 :
1196 21194 : if (.not.deriva) then
1197 9106 : a1 = -y(img24) * y(ihe4) * rate(irmgap)*(1.0d0-rate(irr1))
1198 9106 : a2 = y(isi28) * rate(irr1) * rate(irsigp)
1199 :
1200 9106 : qray(img24) = qray(img24) + a1 + a2
1201 :
1202 : else
1203 12088 : a1 = -y(img24)*y(ihe4) * rate(irmgap)*(1.0d0 - ratdum(irr1))
1204 12088 : a2 = y(img24)*y(ihe4) * ratdum(irmgap)*rate(irr1)
1205 12088 : a3 = y(isi28) * ratdum(irr1) * rate(irsigp)
1206 12088 : a4 = y(isi28) * rate(irr1) * ratdum(irsigp)
1207 :
1208 12088 : qray(img24) = qray(img24) + a1 + a2 + a3 + a4
1209 : end if
1210 :
1211 :
1212 : ! si28 reactions
1213 21194 : a1 = 0.5d0 * y(ic12) * y(io16) * rate(ir1216)
1214 21194 : a2 = 0.56d0 * 0.5d0*y(io16) * y(io16) * rate(ir1616)
1215 21194 : a3 = y(img24) * y(ihe4) * rate(irmgag)
1216 21194 : a4 = -y(isi28) * y(ihe4) * rate(irsiag)
1217 21194 : a5 = -y(isi28) * rate(irsiga)
1218 21194 : a6 = y(is32) * rate(irsga)
1219 :
1220 21194 : qray(isi28) = qray(isi28) + a1 + a2 + a3 + a4 + a5 + a6
1221 :
1222 21194 : if (.not.deriva) then
1223 :
1224 9106 : a1 = 0.34d0*0.5d0*y(io16)*y(io16)*rate(irs1)*rate(ir1616)
1225 9106 : a2 = y(img24) * y(ihe4) * rate(irmgap)*(1.0d0-rate(irr1))
1226 9106 : a3 = -y(isi28) * rate(irr1) * rate(irsigp)
1227 9106 : a4 = -y(isi28) * y(ihe4) * rate(irsiap)*(1.0d0-rate(irs1))
1228 9106 : a5 = y(is32) * rate(irs1) * rate(irsgp)
1229 :
1230 9106 : qray(isi28) = qray(isi28) + a1 + a2 + a3 + a4 + a5
1231 :
1232 : else
1233 12088 : a1 = 0.34d0*0.5d0*y(io16)*y(io16) * ratdum(irs1)*rate(ir1616)
1234 12088 : a2 = 0.34d0*0.5d0*y(io16)*y(io16) * rate(irs1)*ratdum(ir1616)
1235 12088 : a3 = y(img24)*y(ihe4) * rate(irmgap)*(1.0d0 - ratdum(irr1))
1236 12088 : a4 = -y(img24)*y(ihe4) * ratdum(irmgap)*rate(irr1)
1237 12088 : a5 = -y(isi28) * ratdum(irr1) * rate(irsigp)
1238 12088 : a6 = -y(isi28) * rate(irr1) * ratdum(irsigp)
1239 12088 : a7 = -y(isi28)*y(ihe4) * rate(irsiap)*(1.0d0 - ratdum(irs1))
1240 12088 : a8 = y(isi28)*y(ihe4) * ratdum(irsiap)*rate(irs1)
1241 12088 : a9 = y(is32) * ratdum(irs1) * rate(irsgp)
1242 12088 : a10 = y(is32) * rate(irs1) * ratdum(irsgp)
1243 :
1244 : qray(isi28) = qray(isi28) + &
1245 12088 : a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + a10
1246 : end if
1247 :
1248 :
1249 : ! s32 reactions
1250 21194 : a1 = 0.1d0 * 0.5d0*y(io16) * y(io16) * rate(ir1616)
1251 21194 : a2 = y(isi28) * y(ihe4) * rate(irsiag)
1252 21194 : a3 = -y(is32) * y(ihe4) * rate(irsag)
1253 21194 : a4 = -y(is32) * rate(irsga)
1254 21194 : a5 = y(iar36) * rate(irarga)
1255 :
1256 21194 : qray(is32) = qray(is32) + a1 + a2 + a3 + a4 + a5
1257 :
1258 21194 : if (.not.deriva) then
1259 :
1260 9106 : a1 = 0.34d0*0.5d0*y(io16)*y(io16)* rate(ir1616)*(1.0d0-rate(irs1))
1261 9106 : a2 = y(isi28) * y(ihe4) * rate(irsiap)*(1.0d0-rate(irs1))
1262 9106 : a3 = -y(is32) * rate(irs1) * rate(irsgp)
1263 9106 : a4 = -y(is32) * y(ihe4) * rate(irsap)*(1.0d0-rate(irt1))
1264 9106 : a5 = y(iar36) * rate(irt1) * rate(irargp)
1265 :
1266 9106 : qray(is32) = qray(is32) + a1 + a2 + a3 + a4 + a5
1267 :
1268 : else
1269 12088 : a1 = 0.34d0*0.5d0*y(io16)*y(io16) * rate(ir1616)*(1.0d0-ratdum(irs1))
1270 12088 : a2 = -0.34d0*0.5d0*y(io16)*y(io16) * ratdum(ir1616)*rate(irs1)
1271 12088 : a3 = y(isi28)*y(ihe4) * rate(irsiap)*(1.0d0-ratdum(irs1))
1272 12088 : a4 = -y(isi28)*y(ihe4) * ratdum(irsiap)*rate(irs1)
1273 12088 : a5 = -y(is32) * ratdum(irs1) * rate(irsgp)
1274 12088 : a6 = -y(is32) * rate(irs1) * ratdum(irsgp)
1275 12088 : a7 = -y(is32)*y(ihe4) * rate(irsap)*(1.0d0-ratdum(irt1))
1276 12088 : a8 = y(is32)*y(ihe4) * ratdum(irsap)*rate(irt1)
1277 12088 : a9 = y(iar36) * ratdum(irt1) * rate(irargp)
1278 12088 : a10 = y(iar36) * rate(irt1) * ratdum(irargp)
1279 :
1280 : qray(is32) = qray(is32) + &
1281 12088 : a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + a10
1282 : end if
1283 :
1284 :
1285 : ! ar36 reactions
1286 21194 : a1 = y(is32) * y(ihe4) * rate(irsag)
1287 21194 : a2 = -y(iar36) * y(ihe4) * rate(irarag)
1288 21194 : a3 = -y(iar36) * rate(irarga)
1289 21194 : a4 = y(ica40) * rate(ircaga)
1290 :
1291 21194 : qray(iar36) = qray(iar36) + a1 + a2 + a3 + a4
1292 :
1293 21194 : if (.not.deriva) then
1294 9106 : a1 = y(is32) * y(ihe4) * rate(irsap)*(1.0d0-rate(irt1))
1295 9106 : a2 = -y(iar36) * rate(irt1) * rate(irargp)
1296 9106 : a3 = -y(iar36) * y(ihe4) * rate(irarap)*(1.0d0-rate(iru1))
1297 9106 : a4 = y(ica40) * rate(ircagp) * rate(iru1)
1298 :
1299 9106 : qray(iar36) = qray(iar36) + a1 + a2 + a3 + a4
1300 :
1301 : else
1302 12088 : a1 = y(is32)*y(ihe4) * rate(irsap)*(1.0d0 - ratdum(irt1))
1303 12088 : a2 = -y(is32)*y(ihe4) * ratdum(irsap)*rate(irt1)
1304 12088 : a3 = -y(iar36) * ratdum(irt1) * rate(irargp)
1305 12088 : a4 = -y(iar36) * rate(irt1) * ratdum(irargp)
1306 12088 : a5 = -y(iar36)*y(ihe4) * rate(irarap)*(1.0d0-ratdum(iru1))
1307 12088 : a6 = y(iar36)*y(ihe4) * ratdum(irarap)*rate(iru1)
1308 12088 : a7 = y(ica40) * ratdum(ircagp) * rate(iru1)
1309 12088 : a8 = y(ica40) * rate(ircagp) * ratdum(iru1)
1310 :
1311 12088 : qray(iar36) = qray(iar36) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8
1312 : end if
1313 :
1314 :
1315 : ! ca40 reactions
1316 21194 : a1 = y(iar36) * y(ihe4) * rate(irarag)
1317 21194 : a2 = -y(ica40) * y(ihe4) * rate(ircaag)
1318 21194 : a3 = -y(ica40) * rate(ircaga)
1319 21194 : a4 = y(iti44) * rate(irtiga)
1320 :
1321 21194 : qray(ica40) = qray(ica40) + a1 + a2 + a3 + a4
1322 :
1323 21194 : if (.not.deriva) then
1324 :
1325 9106 : a1 = y(iar36) * y(ihe4) * rate(irarap)*(1.0d0-rate(iru1))
1326 9106 : a2 = -y(ica40) * rate(ircagp) * rate(iru1)
1327 9106 : a3 = -y(ica40) * y(ihe4) * rate(ircaap)*(1.0d0-rate(irv1))
1328 9106 : a4 = y(iti44) * rate(irtigp) * rate(irv1)
1329 :
1330 9106 : qray(ica40) = qray(ica40) + a1 + a2 + a3 + a4
1331 :
1332 : else
1333 12088 : a1 = y(iar36)*y(ihe4) * rate(irarap)*(1.0d0-ratdum(iru1))
1334 12088 : a2 = -y(iar36)*y(ihe4) * ratdum(irarap)*rate(iru1)
1335 12088 : a3 = -y(ica40) * ratdum(ircagp) * rate(iru1)
1336 12088 : a4 = -y(ica40) * rate(ircagp) * ratdum(iru1)
1337 12088 : a5 = -y(ica40)*y(ihe4) * rate(ircaap)*(1.0d0-ratdum(irv1))
1338 12088 : a6 = y(ica40)*y(ihe4) * ratdum(ircaap)*rate(irv1)
1339 12088 : a7 = y(iti44) * ratdum(irtigp) * rate(irv1)
1340 12088 : a8 = y(iti44) * rate(irtigp) * ratdum(irv1)
1341 :
1342 12088 : qray(ica40) = qray(ica40) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8
1343 : end if
1344 :
1345 :
1346 : ! ti44 reactions
1347 21194 : a1 = y(ica40) * y(ihe4) * rate(ircaag)
1348 21194 : a2 = -y(iti44) * y(ihe4) * rate(irtiag)
1349 21194 : a3 = -y(iti44) * rate(irtiga)
1350 21194 : a4 = y(icr48) * rate(ircrga)
1351 :
1352 21194 : qray(iti44) = qray(iti44) + a1 + a2 + a3 + a4
1353 :
1354 21194 : if (.not.deriva) then
1355 9106 : a1 = y(ica40) * y(ihe4) * rate(ircaap)*(1.0d0-rate(irv1))
1356 9106 : a2 = -y(iti44) * rate(irv1) * rate(irtigp)
1357 9106 : a3 = -y(iti44) * y(ihe4) * rate(irtiap)*(1.0d0-rate(irw1))
1358 9106 : a4 = y(icr48) * rate(irw1) * rate(ircrgp)
1359 :
1360 9106 : qray(iti44) = qray(iti44) + a1 + a2 + a3 + a4
1361 :
1362 : else
1363 12088 : a1 = y(ica40)*y(ihe4) * rate(ircaap)*(1.0d0-ratdum(irv1))
1364 12088 : a2 = -y(ica40)*y(ihe4) * ratdum(ircaap)*rate(irv1)
1365 12088 : a3 = -y(iti44) * ratdum(irv1) * rate(irtigp)
1366 12088 : a4 = -y(iti44) * rate(irv1) * ratdum(irtigp)
1367 12088 : a5 = -y(iti44)*y(ihe4) * rate(irtiap)*(1.0d0-ratdum(irw1))
1368 12088 : a6 = y(iti44)*y(ihe4) * ratdum(irtiap)*rate(irw1)
1369 12088 : a7 = y(icr48) * ratdum(irw1) * rate(ircrgp)
1370 12088 : a8 = y(icr48) * rate(irw1) * ratdum(ircrgp)
1371 :
1372 12088 : qray(iti44) = qray(iti44) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8
1373 : end if
1374 :
1375 :
1376 : ! cr48 reactions
1377 21194 : a1 = y(iti44) * y(ihe4) * rate(irtiag)
1378 21194 : a2 = -y(icr48) * y(ihe4) * rate(ircrag)
1379 21194 : a3 = -y(icr48) * rate(ircrga)
1380 21194 : a4 = y(ife52) * rate(irfega)
1381 :
1382 21194 : qray(icr48) = qray(icr48) + a1 + a2 + a3 + a4
1383 :
1384 21194 : if (.not.deriva) then
1385 9106 : a1 = y(iti44) * y(ihe4) * rate(irtiap)*(1.0d0-rate(irw1))
1386 9106 : a2 = -y(icr48) * rate(irw1) * rate(ircrgp)
1387 9106 : a3 = -y(icr48) * y(ihe4) * rate(ircrap)*(1.0d0-rate(irx1))
1388 9106 : a4 = y(ife52) * rate(irx1) * rate(irfegp)
1389 :
1390 9106 : qray(icr48) = qray(icr48) + a1 + a2 + a3 + a4
1391 :
1392 : else
1393 12088 : a1 = y(iti44)*y(ihe4) * rate(irtiap)*(1.0d0-ratdum(irw1))
1394 12088 : a2 = -y(iti44)*y(ihe4) * ratdum(irtiap)*rate(irw1)
1395 12088 : a3 = -y(icr48) * ratdum(irw1) * rate(ircrgp)
1396 12088 : a4 = -y(icr48) * rate(irw1) * ratdum(ircrgp)
1397 12088 : a5 = -y(icr48)*y(ihe4) * rate(ircrap)*(1.0d0-ratdum(irx1))
1398 12088 : a6 = y(icr48)*y(ihe4) * ratdum(ircrap)*rate(irx1)
1399 12088 : a7 = y(ife52) * ratdum(irx1) * rate(irfegp)
1400 12088 : a8 = y(ife52) * rate(irx1) * ratdum(irfegp)
1401 :
1402 12088 : qray(icr48) = qray(icr48) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8
1403 : end if
1404 :
1405 :
1406 : ! crx reactions
1407 21194 : a1 = y(ife56) * fe56ec_fake_factor * rate(irn56ec)
1408 :
1409 21194 : qray(icrx) = qray(icrx) + a1
1410 :
1411 : ! fe52 reactions
1412 21194 : a1 = y(icr48) * y(ihe4) * rate(ircrag)
1413 21194 : a2 = -y(ife52) * y(ihe4) * rate(irfeag)
1414 21194 : a3 = -y(ife52) * rate(irfega)
1415 21194 : a4 = y(ini56) * rate(irniga)
1416 :
1417 21194 : qray(ife52) = qray(ife52) + a1 + a2 + a3 + a4
1418 :
1419 21194 : if (.not.deriva) then
1420 9106 : a1 = y(icr48) * y(ihe4) * rate(ircrap)*(1.0d0-rate(irx1))
1421 9106 : a2 = -y(ife52) * rate(irx1) * rate(irfegp)
1422 :
1423 9106 : qray(ife52) = qray(ife52) + a1 + a2
1424 :
1425 : else
1426 12088 : a1 = y(icr48)*y(ihe4) * rate(ircrap)*(1.0d0-ratdum(irx1))
1427 12088 : a2 = -y(icr48)*y(ihe4) * ratdum(ircrap)*rate(irx1)
1428 12088 : a3 = -y(ife52) * ratdum(irx1) * rate(irfegp)
1429 12088 : a4 = -y(ife52) * rate(irx1) * ratdum(irfegp)
1430 :
1431 12088 : qray(ife52) = qray(ife52) + a1 + a2 + a3 + a4
1432 : end if
1433 :
1434 21194 : a1 = y(ife54) * rate(ir1f54)
1435 21194 : a2 = -y(ife52) * y(ineut) * y(ineut) * rate(ir2f54)
1436 21194 : a3 = y(ife54) * y(iprot) * y(iprot) * rate(ir5f54)
1437 21194 : a4 = -y(ife52) * y(ihe4) * rate(ir6f54)
1438 21194 : a5 = -y(ife52) * y(ihe4) * y(iprot) * rate(ir7f54)
1439 21194 : a6 = y(ini56) * y(iprot) * rate(ir8f54)
1440 :
1441 21194 : qray(ife52) = qray(ife52) + a1 + a2 + a3 + a4 + a5 + a6
1442 :
1443 :
1444 : ! fe54 reactions
1445 21194 : a1 = -y(ife54) * rate(ir1f54)
1446 21194 : a2 = y(ife52) * y(ineut) * y(ineut) * rate(ir2f54)
1447 21194 : a3 = -y(ife54) * y(iprot) * y(iprot) * rate(ir3f54)
1448 21194 : a4 = y(ini56) * rate(ir4f54)
1449 21194 : a5 = -y(ife54) * y(iprot) * y(iprot) * rate(ir5f54)
1450 21194 : a6 = y(ife52) * y(ihe4) * rate(ir6f54)
1451 21194 : a7 = y(ife56) * rate(irfe56_aux1)
1452 21194 : a8 = -y(ife54) * y(ineut) * y(ineut) * rate(irfe56_aux2)
1453 21194 : a9 = y(ife56) * y(iprot) * y(iprot) * rate(irfe56_aux3)
1454 21194 : a10 = -y(ife54) * y(ihe4) * rate(irfe56_aux4)
1455 :
1456 : qray(ife54) = qray(ife54) + &
1457 21194 : a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + a10
1458 :
1459 :
1460 : ! fe56 reactions
1461 21194 : if (plus_co56) then
1462 12 : a1 = y(ico56) * rate(irco56ec)
1463 : else
1464 21182 : a1 = y(ini56) * rate(irn56ec)
1465 : end if
1466 21194 : a2 = -y(ife56) * fe56ec_fake_factor * rate(irn56ec)
1467 21194 : a3 = -y(ife56) * rate(irfe56_aux1)
1468 21194 : a4 = y(ife54) * y(ineut) * y(ineut) * rate(irfe56_aux2)
1469 21194 : a5 = -y(ife56) * y(iprot) * y(iprot) * rate(irfe56_aux3)
1470 21194 : a6 = y(ife54) * y(ihe4) * rate(irfe56_aux4)
1471 :
1472 21194 : qray(ife56) = qray(ife56) + a1 + a2 + a3 + a4 + a5 + a6
1473 :
1474 21194 : if (plus_co56) then
1475 : ! co56 reactions
1476 12 : a1 = y(ini56) * rate(irn56ec)
1477 12 : a2 = -y(ico56) * rate(irco56ec)
1478 :
1479 12 : qray(ico56) = qray(ico56) + a1 + a2
1480 : end if
1481 :
1482 : ! ni56 reactions
1483 21194 : a1 = y(ife52) * y(ihe4) * rate(irfeag)
1484 21194 : a2 = -y(ini56) * rate(irniga)
1485 21194 : a3 = -y(ini56) * rate(irn56ec)
1486 :
1487 21194 : qray(ini56) = qray(ini56) + a1 + a2 + a3
1488 :
1489 21194 : a1 = y(ife54) * y(iprot) * y(iprot) * rate(ir3f54)
1490 21194 : a2 = -y(ini56) * rate(ir4f54)
1491 21194 : a3 = y(ife52) * y(ihe4)* y(iprot) * rate(ir7f54)
1492 21194 : a4 = -y(ini56) * y(iprot) * rate(ir8f54)
1493 :
1494 21194 : qray(ini56) = qray(ini56) + a1 + a2 + a3 + a4
1495 :
1496 : ! neutrons
1497 21194 : a1 = 2.0d0 * y(ife54) * rate(ir1f54)
1498 21194 : a2 = -2.0d0 * y(ife52) * y(ineut) * y(ineut) * rate(ir2f54)
1499 21194 : a3 = 2.0d0 * y(ihe4) * rate(iralf1)
1500 21194 : a4 = -2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * rate(iralf2)
1501 21194 : a5 = y(iprot) * rate(irpen)
1502 21194 : a6 = -y(ineut) * rate(irnep)
1503 21194 : a7 = 2.0d0 * y(ife56) * rate(irfe56_aux1)
1504 21194 : a8 = -2.0d0 * y(ife54) * y(ineut) * y(ineut) * rate(irfe56_aux2)
1505 21194 : a9 = -fe56ec_n_neut * y(ife56) * fe56ec_fake_factor * rate(irn56ec)
1506 :
1507 21194 : qray(ineut) = qray(ineut) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9
1508 :
1509 : ! photodisintegration protons
1510 21194 : a1 = -2.0d0 * y(ife54) * y(iprot) * y(iprot) * rate(ir3f54)
1511 21194 : a2 = 2.0d0 * y(ini56) * rate(ir4f54)
1512 21194 : a3 = -2.0d0 * y(ife54) * y(iprot) * y(iprot) * rate(ir5f54)
1513 21194 : a4 = 2.0d0 * y(ife52) * y(ihe4) * rate(ir6f54)
1514 21194 : a5 = 2.0d0 * y(ihe4) * rate(iralf1)
1515 21194 : a6 = -2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * rate(iralf2)
1516 21194 : a7 = -y(iprot) * rate(irpen)
1517 21194 : a8 = y(ineut) * rate(irnep)
1518 21194 : a9 = -2.0d0 * y(ife56) * y(iprot) * y(iprot) * rate(irfe56_aux3)
1519 21194 : a10 = 2.0d0 * y(ife54) * y(ihe4) * rate(irfe56_aux4)
1520 :
1521 : qray(iprot) = qray(iprot) + &
1522 21194 : a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + a10
1523 :
1524 : ! now set the real(dp) return argument dydt
1525 21194 : okay = .true.
1526 466280 : do i=1,species(plus_co56)
1527 445086 : dydt(i) = qray(i)
1528 466280 : if (is_bad(dydt(i))) then
1529 0 : write(*,*) 'dydt(i)', i, dydt(i), y(i)
1530 0 : okay = .false.
1531 : end if
1532 : end do
1533 21194 : if (.not. okay) then
1534 :
1535 0 : write(*,*) 'log10(temp) = ',safe_log10(temp)
1536 0 : write(*,*) 'log10(rho) = ',safe_log10(den)
1537 :
1538 0 : do i=1,num_reactions(plus_co56)
1539 0 : write(*,*) trim(ratnam(i)), i, rate(i)
1540 : end do
1541 0 : call mesa_error(__FILE__,__LINE__,'approx21_dydt')
1542 : end if
1543 :
1544 21194 : end subroutine approx21_dydt
1545 :
1546 :
1547 9106 : real(dp) function approx21_eval_PPII_fraction(y, rate) result(fII)
1548 : real(dp), dimension(:), intent(in) :: y, rate
1549 9106 : real(dp) :: rateII, rateIII, rsum
1550 : include 'formats'
1551 9106 : rateII = rate(ir_be7_wk_li7)
1552 9106 : rateIII = y(ih1) * rate(ir_be7_pg_b8)
1553 9106 : rsum = rateII + rateIII
1554 9106 : if (rsum < 1d-50) then
1555 : fII = 0.5d0
1556 : else
1557 8362 : fII = rateII / rsum
1558 : end if
1559 9106 : end function approx21_eval_PPII_fraction
1560 :
1561 :
1562 21194 : subroutine approx21_eps_info( &
1563 21194 : n, y, mion, dydt, rate, fII, &
1564 : Qtotal_rpp, Qneu_rpp, &
1565 : Qr33, &
1566 : Qtotal_rpp2, Qneu_rpp2, &
1567 : Qtotal_rpp3, Qneu_rpp3, &
1568 : Qtotal_rcpg, Qneu_rcpg, &
1569 : Qtotal_rnpg, Qneu_rnpg, &
1570 : Qtotal_ropg, Qneu_ropg, &
1571 : Qrn14_to_o16, &
1572 : Qtotal_rpen, Qneu_rpen, &
1573 : Qtotal_rnep, Qneu_rnep, &
1574 : Qtotal_rni56ec, Qneu_rni56ec, &
1575 : Qtotal_rco56ec, Qneu_rco56ec, &
1576 : Qtotal_rfe56ec, Qneu_rfe56ec, &
1577 : Qrn14ag, &
1578 : Qr3alf, &
1579 : Qrc12ag, Qro16ag,&
1580 : Qr1212, &
1581 : Qr1216_to_mg24, Qr1216_to_si28, &
1582 : Qr1616a, Qr1616g, &
1583 : Qrne20ag, &
1584 : Qrmg24ag, &
1585 : Qrsi28ag, &
1586 : Qrs32ag, &
1587 : Qrar36ag, &
1588 : Qrca40ag, &
1589 : Qrti44ag, &
1590 : Qrcr48ag, &
1591 : Qrfe52ag, &
1592 : Qrfe52ng, &
1593 : Qrfe53ng, &
1594 : Qrfe54ng, &
1595 : Qrfe55ng, &
1596 : Qrfe52neut_to_fe54, &
1597 : Qrfe52aprot_to_fe54, &
1598 : Qrfe54ng_to_fe56, &
1599 : Qrfe54aprot_to_fe56, &
1600 : Qrfe52aprot_to_ni56, &
1601 : Qrfe54prot_to_ni56, &
1602 : Qrhe4_breakup, &
1603 : Qrhe4_rebuild, &
1604 : eps_total, eps_neu, &
1605 21194 : do_eps_nuc_categories, eps_nuc_categories, &
1606 : dbg, &
1607 : plus_co56, &
1608 : ierr)
1609 : use const_def, only: Qconv
1610 : use chem_def, only: num_categories, category_name, chem_isos, &
1611 : ipp, icno, i3alf, i_burn_c, i_burn_n, i_burn_o, i_burn_ne, i_burn_na, &
1612 : i_burn_mg, i_burn_si, i_burn_s, i_burn_ar, i_burn_ca, i_burn_ti, i_burn_cr, &
1613 : i_burn_fe, icc, ico, ioo, ipnhe4, iphoto, i_ni56_co56, i_co56_fe56, iother
1614 : use net_def, only: Net_Info
1615 : type (Net_Info) :: n
1616 : real(dp), dimension(:), intent(in) :: y, mion, dydt, rate
1617 : real(dp), intent(in) :: fII, &
1618 : Qtotal_rpp, Qneu_rpp, Qr33, &
1619 : Qtotal_rpp2, Qneu_rpp2, &
1620 : Qtotal_rpp3, Qneu_rpp3, &
1621 : Qtotal_rcpg, Qneu_rcpg, &
1622 : Qtotal_rnpg, Qneu_rnpg, &
1623 : Qtotal_ropg, Qneu_ropg, &
1624 : Qrn14_to_o16, Qtotal_rpen, Qneu_rpen, &
1625 : Qtotal_rnep, Qneu_rnep, &
1626 : Qtotal_rni56ec, Qneu_rni56ec, &
1627 : Qtotal_rco56ec, Qneu_rco56ec, &
1628 : Qtotal_rfe56ec, Qneu_rfe56ec, &
1629 : Qrn14ag, &
1630 : Qr3alf, &
1631 : Qrc12ag, Qro16ag, &
1632 : Qr1212, &
1633 : Qr1216_to_mg24, Qr1216_to_si28, &
1634 : Qr1616a, Qr1616g, &
1635 : Qrne20ag, &
1636 : Qrmg24ag, &
1637 : Qrsi28ag, &
1638 : Qrs32ag, &
1639 : Qrar36ag, &
1640 : Qrca40ag, &
1641 : Qrti44ag, &
1642 : Qrcr48ag, &
1643 : Qrfe52ag, &
1644 : Qrfe52ng, &
1645 : Qrfe53ng, &
1646 : Qrfe54ng, &
1647 : Qrfe55ng, &
1648 : Qrfe52neut_to_fe54, &
1649 : Qrfe52aprot_to_fe54, &
1650 : Qrfe54ng_to_fe56, &
1651 : Qrfe54aprot_to_fe56, &
1652 : Qrfe52aprot_to_ni56, &
1653 : Qrfe54prot_to_ni56, &
1654 : Qrhe4_breakup, &
1655 : Qrhe4_rebuild
1656 : logical, intent(in) :: do_eps_nuc_categories
1657 : real(dp), intent(out) :: eps_total, eps_neu, eps_nuc_categories(:)
1658 : logical, intent(in) :: dbg
1659 : integer, intent(out) :: ierr
1660 :
1661 : integer :: i, fe56ec_n_neut
1662 529850 : real(qp) :: a1, a2, xx, eps_neu_q, eps_nuc_cat(num_categories), eps_total_q, &
1663 21194 : eps_nuc_q, sum_categories_q
1664 21194 : real(dp) :: enuc_conv2, sum_categories, eps_nuc, fe56ec_fake_factor
1665 : logical, intent(in) :: plus_co56
1666 :
1667 : include 'formats'
1668 :
1669 : !write(*,1) 'reaction_Qs(irn14_to_o16) Qrn14_to_o16*Qconv', Qrn14_to_o16*Qconv
1670 :
1671 21194 : ierr = 0
1672 :
1673 21194 : xx = 0.0_qp
1674 466280 : do i=1,species(plus_co56)
1675 445086 : a1 = dydt(i)
1676 445086 : a2 = mion(i)
1677 466280 : xx = xx + a1*a2
1678 : end do
1679 21194 : eps_total_q = -m3(avo,clight,clight) * xx
1680 21194 : eps_total = eps_total_q
1681 :
1682 : fe56ec_fake_factor = eval_fe56ec_fake_factor( &
1683 42388 : n% g% fe56ec_fake_factor, n% g% min_T_for_fe56ec_fake_factor, n% temp)
1684 21194 : fe56ec_n_neut = n% g% fe56ec_n_neut
1685 :
1686 : eps_neu_q = &
1687 : m5(Qneu_rpp, 0.5d0, y(ih1), y(ih1), rate(irpp)) + &
1688 : m5(Qneu_rpp2, y(ihe3), y(ihe4), rate(irhe3ag), fII) + &
1689 : m5(Qneu_rpp3, y(ihe3), y(ihe4), rate(irhe3ag), (1d0-fII)) + &
1690 : m4(Qneu_rcpg, y(ic12), y(ih1), rate(ircpg)) + &
1691 : m5(Qneu_rnpg, y(in14), y(ih1), rate(irnpg), rate(ifa)) + &
1692 : m4(Qneu_ropg, y(io16), y(ih1), rate(iropg)) + &
1693 : m3(Qneu_rpen, y(ih1), rate(irpen)) + &
1694 : m3(Qneu_rpen, y(iprot), rate(irpen)) + &
1695 : m3(Qneu_rnep, y(ineut), rate(irnep)) + &
1696 21194 : m4(Qneu_rfe56ec, y(ife56), fe56ec_fake_factor, rate(irn56ec))
1697 :
1698 21194 : if (plus_co56) then
1699 : eps_neu_q = eps_neu_q + &
1700 : m3(Qneu_rni56ec, y(ini56), rate(irn56ec)) + &
1701 12 : m3(Qneu_rco56ec, y(ico56), rate(irco56ec))
1702 : else
1703 : eps_neu_q = eps_neu_q + &
1704 21182 : m3(Qneu_rni56ec + Qneu_rco56ec, y(ini56), rate(irn56ec))
1705 : end if
1706 21194 : eps_neu_q = eps_neu_q * Qconv
1707 21194 : eps_neu = eps_neu_q
1708 :
1709 21194 : eps_nuc_q = eps_total_q - eps_neu_q
1710 21194 : eps_nuc = eps_nuc_q
1711 :
1712 21194 : if (.not. do_eps_nuc_categories) return
1713 :
1714 227650 : do i=1,num_categories
1715 227650 : eps_nuc_cat(i) = 0d0
1716 : end do
1717 :
1718 : ! Set the rates
1719 1065406 : n% raw_rate = 0d0
1720 1065406 : n% screened_rate = 0d0
1721 1065406 : n% eps_nuc_rate = 0d0
1722 1065406 : n% eps_neu_rate = 0d0
1723 :
1724 : ! Set neutrinos first
1725 9106 : n% eps_neu_rate(irpp) = m5(Qneu_rpp, 0.5d0, y(ih1), y(ih1), rate(irpp))
1726 : n% eps_neu_rate(irhe3ag) = m5(Qneu_rpp2, y(ihe3), y(ihe4), rate(irhe3ag), fII) + &
1727 9106 : m5(Qneu_rpp3, y(ihe3), y(ihe4), rate(irhe3ag), (1d0-fII))
1728 9106 : n% eps_neu_rate(ircpg) = m4(Qneu_rcpg, y(ic12), y(ih1), rate(ircpg))
1729 9106 : n% eps_neu_rate(irnpg) = m5(Qneu_rnpg, y(in14), y(ih1), rate(irnpg), rate(ifa))
1730 9106 : n% eps_neu_rate(iropg) = m4(Qneu_ropg, y(io16), y(ih1), rate(iropg))
1731 : n% eps_neu_rate(irpen) = m3(Qneu_rpen, y(ih1), rate(irpen)) + &
1732 9106 : m3(Qneu_rpen, y(iprot), rate(irpen))
1733 9106 : n% eps_neu_rate(irnep) = m3(Qneu_rnep, y(ineut), rate(irnep))
1734 9106 : n% eps_neu_rate(irn56ec) = m4(Qneu_rfe56ec, y(ife56), fe56ec_fake_factor, rate(irn56ec))
1735 9106 : if (plus_co56) then
1736 : n% eps_neu_rate(irn56ec) = n% eps_neu_rate(irn56ec) + &
1737 : m3(Qneu_rni56ec, y(ini56), rate(irn56ec)) + &
1738 4 : m3(Qneu_rco56ec, y(ico56), rate(irco56ec))
1739 : else
1740 : n% eps_neu_rate(irn56ec) = n% eps_neu_rate(irn56ec) + &
1741 9102 : m3(Qneu_rni56ec + Qneu_rco56ec, y(ini56), rate(irn56ec))
1742 : end if
1743 :
1744 1065406 : n% eps_neu_rate = n% eps_neu_rate * Qconv
1745 :
1746 36424 : call set_eps_nuc(Qtotal_rpp - Qneu_rpp,[0.5d0, y(ih1), y(ih1)],irpp, ipp)
1747 36424 : call set_eps_nuc(Qr33, [0.5d0, y(ihe3), y(ihe3)], ir33, ipp)
1748 : call set_eps_nuc(( &
1749 : (Qtotal_rpp2 - Qneu_rpp2)*fII + &
1750 : (Qtotal_rpp3 - Qneu_rpp3)*(1d0 - fII)), &
1751 27318 : [y(ihe3), y(ihe4)],irhe3ag, ipp)
1752 :
1753 27318 : call set_eps_nuc(Qtotal_rcpg - Qneu_rcpg, [y(ic12), y(ih1)],ircpg,icno)
1754 27318 : call set_eps_nuc(Qtotal_rcpg - Qneu_rnpg, [y(in14), y(ih1)],irnpg,icno)
1755 36424 : call set_eps_nuc(Qtotal_ropg - Qneu_ropg, [y(io16), y(ih1),rate(ifa)],iropg,icno)
1756 :
1757 45530 : call set_eps_nuc(Qr3alf, [1d0/6d0,y(ihe4), y(ihe4), y(ihe4)],ir3a,i3alf)
1758 :
1759 27318 : call set_eps_nuc(Qrc12ag, [y(ic12), y(ihe4)],ircag,i_burn_c)
1760 :
1761 27318 : call set_eps_nuc(Qrn14ag, [ y(ihe4), y(in14)],irnag,i_burn_n)
1762 36424 : call set_eps_nuc(Qrn14_to_o16, [y(in14), y(ih1),rate(ifg)],irnpg,i_burn_n)
1763 :
1764 27318 : call set_eps_nuc(Qro16ag, [y(io16), y(ihe4)], iroag, i_burn_o)
1765 :
1766 36424 : call set_eps_nuc(Qr1212, [0.5d0,y(ic12), y(ic12)],ir1212,icc)
1767 :
1768 27318 : call set_eps_nuc(0.5d0*(Qr1216_to_mg24 + Qr1216_to_si28), [y(ic12), y(io16)], ir1216, ico )
1769 :
1770 : ! these make he4 + si28
1771 36424 : call set_eps_nuc( Qr1616a * (0.56d0 + 0.34d0*rate(irs1)), [0.5d0,y(io16), y(io16)], ir1616, ioo)
1772 : ! these make s32
1773 36424 : call set_eps_nuc( Qr1616g * (0.1d0 + 0.34d0*(1d0 - rate(irs1))) , [0.5d0,y(io16), y(io16)], ir1616, ioo )
1774 :
1775 27318 : call set_eps_nuc(Qrne20ag, [y(ihe4), y(ine20)], irneag, i_burn_ne)
1776 :
1777 27318 : call set_eps_nuc(Qrmg24ag, [y(ihe4), y(img24)],irmgag,i_burn_mg)
1778 36424 : call set_eps_nuc(Qrmg24ag, [y(ihe4), y(img24),1.0d0-rate(irr1)],irmgap,i_burn_mg)
1779 :
1780 27318 : call set_eps_nuc(Qrsi28ag, [y(ihe4), y(isi28)],irsiag,i_burn_si)
1781 36424 : call set_eps_nuc(Qrsi28ag, [y(ihe4), y(isi28),(1.0d0-rate(irs1))],irsiap,i_burn_si)
1782 :
1783 27318 : call set_eps_nuc(Qrs32ag, [y(ihe4), y(is32)], irsag, i_burn_s)
1784 36424 : call set_eps_nuc(Qrs32ag, [y(ihe4), y(is32),(1.0d0-rate(irt1))], irsap, i_burn_s)
1785 :
1786 27318 : call set_eps_nuc(Qrar36ag, [y(ihe4), y(iar36)], irarag, i_burn_ar)
1787 36424 : call set_eps_nuc(Qrar36ag, [y(ihe4), y(iar36),(1.0d0-rate(iru1))], irarap, i_burn_ar)
1788 :
1789 27318 : call set_eps_nuc(Qrca40ag, [y(ihe4), y(ica40)], ircaag, i_burn_ca)
1790 36424 : call set_eps_nuc(Qrca40ag, [y(ihe4), y(ica40), (1.0d0-rate(irv1))], ircaap, i_burn_ca)
1791 :
1792 27318 : call set_eps_nuc(Qrti44ag, [y(ihe4), y(iti44)], irtiag, i_burn_ti)
1793 36424 : call set_eps_nuc(Qrti44ag, [y(ihe4), y(iti44),(1.0d0-rate(irw1))], irtiap, i_burn_ti)
1794 :
1795 27318 : call set_eps_nuc(Qrcr48ag, [y(ihe4), y(icr48)], ircrag, i_burn_cr)
1796 36424 : call set_eps_nuc(Qrcr48ag, [y(ihe4), y(icr48),(1.0d0-rate(irx1))], ircrap, i_burn_cr)
1797 :
1798 27318 : call set_eps_nuc(Qrfe52ag, [y(ihe4), y(ife52)], irfeag, i_burn_fe)
1799 36424 : call set_eps_nuc(Qrfe52aprot_to_ni56, [y(ife52), y(ihe4), y(iprot)], ir7f54, i_burn_fe)
1800 36424 : call set_eps_nuc(Qrfe52neut_to_fe54, [ y(ife52), y(ineut), y(ineut)], ir2f54, i_burn_fe)
1801 36424 : call set_eps_nuc(Qrfe54ng_to_fe56, [ y(ife54), y(ineut), y(ineut)], irfe56_aux2, i_burn_fe)
1802 27318 : call set_eps_nuc(Qrfe52aprot_to_fe54, [y(ife52), y(ihe4)], ir6f54, i_burn_fe)
1803 27318 : call set_eps_nuc(Qrfe54aprot_to_fe56, [y(ife54), y(ihe4)], irfe56_aux4, i_burn_fe)
1804 36424 : call set_eps_nuc(Qrfe54prot_to_ni56, [y(ife54), y(iprot),y(iprot)], ir3f54, i_burn_fe)
1805 27318 : call set_eps_nuc(Qtotal_rfe56ec - Qneu_rfe56ec, [y(ife56),fe56ec_fake_factor], irn56ec, i_burn_fe)
1806 :
1807 45530 : call set_eps_nuc(Qrhe4_rebuild, [y(ineut),y(ineut), y(iprot),y(iprot)],iralf2,ipnhe4)
1808 :
1809 18212 : call set_eps_nuc(Qtotal_rpen, [y(ih1)],irpen,iother)
1810 18212 : call set_eps_nuc(Qtotal_rpen, [y(iprot)],irpen,iother)
1811 18212 : call set_eps_nuc(Qtotal_rnep, [y(ineut)],irnep,iother)
1812 :
1813 : ! m4(Qrfe52aprot_to_ni56, y(ini56), y(iprot), rate(ir8f54)) + &
1814 : ! m5(Qrfe52aprot_to_fe54, y(ife54), y(iprot), y(iprot), rate(ir5f54)) + &
1815 : ! m3(Qrfe52ag, y(ini56), rate(irniga)) + &
1816 : ! m3(Qrfe52neut_to_fe54, y(ife54), rate(ir1f54)) + &
1817 : ! m3(Qrfe54ng_to_fe56, y(ife56), rate(irfe56_aux1)) + &
1818 : ! m5(Qrfe54aprot_to_fe56, y(ife56), y(iprot), y(iprot), rate(irfe56_aux3)) + &
1819 : ! m3(Qrfe54prot_to_ni56, y(ini56), rate(ir4f54)))
1820 :
1821 18212 : call set_eps_nuc(Qrhe4_breakup,[ y(ihe4)],iralf1,iphoto) ! note: Qrhe4_breakup < 0
1822 18212 : call set_eps_nuc(-Qrc12ag,[ y(io16)],iroga, iphoto) ! all the rest are > 0 Q's for forward reactions
1823 18212 : call set_eps_nuc(-Qr3alf,[ y(ic12)],irg3a, iphoto)
1824 18212 : call set_eps_nuc(-Qro16ag,[ y(ine20)],irnega, iphoto)
1825 18212 : call set_eps_nuc(-Qrne20ag,[ y(img24)],irmgga, iphoto)
1826 :
1827 18212 : call set_eps_nuc(-Qrmg24ag,[ y(isi28)],irsiga, iphoto)
1828 27318 : call set_eps_nuc(-Qrmg24ag,[ y(isi28),rate(irr1)],irsigp, iphoto)
1829 :
1830 18212 : call set_eps_nuc(-Qrsi28ag,[ y(is32)],irsga, iphoto)
1831 27318 : call set_eps_nuc(-Qrsi28ag,[ y(is32),rate(irs1)],irsgp, iphoto)
1832 :
1833 18212 : call set_eps_nuc(-Qrs32ag,[ y(iar36)],irarga, iphoto)
1834 27318 : call set_eps_nuc(-Qrs32ag,[ y(iar36),rate(irt1)],irargp, iphoto)
1835 :
1836 18212 : call set_eps_nuc(-Qrar36ag,[ y(ica40)],ircaga, iphoto)
1837 27318 : call set_eps_nuc(-Qrar36ag,[ y(ica40),rate(iru1)],ircagp, iphoto)
1838 :
1839 18212 : call set_eps_nuc(-Qrca40ag,[ y(iti44)],irtiga, iphoto)
1840 27318 : call set_eps_nuc(-Qrca40ag,[ y(iti44),rate(irv1)],irtigp, iphoto)
1841 :
1842 18212 : call set_eps_nuc(-Qrti44ag,[ y(icr48)],ircrga, iphoto)
1843 27318 : call set_eps_nuc(-Qrti44ag,[ y(icr48),rate(irw1)],ircrgp, iphoto)
1844 :
1845 18212 : call set_eps_nuc(-Qrcr48ag,[ y(ife52)],irfega, iphoto)
1846 27318 : call set_eps_nuc(-Qrcr48ag,[ y(ife52),rate(irx1)],irfegp, iphoto)
1847 :
1848 27318 : call set_eps_nuc(-Qrfe52aprot_to_ni56,[ y(ini56), y(iprot)],ir8f54, iphoto)
1849 36424 : call set_eps_nuc(-Qrfe52aprot_to_fe54,[ y(ife54), y(iprot), y(iprot)],ir5f54, iphoto)
1850 18212 : call set_eps_nuc(-Qrfe52ag,[ y(ini56)],irniga, iphoto)
1851 18212 : call set_eps_nuc(-Qrfe52neut_to_fe54,[ y(ife54)],ir1f54, iphoto)
1852 18212 : call set_eps_nuc(-Qrfe54ng_to_fe56,[ y(ife56)],irfe56_aux1, iphoto)
1853 36424 : call set_eps_nuc(-Qrfe54aprot_to_fe56,[ y(ife56), y(iprot), y(iprot)],irfe56_aux3, iphoto)
1854 18212 : call set_eps_nuc(-Qrfe54prot_to_ni56,[ y(ini56)],ir4f54, iphoto)
1855 :
1856 :
1857 :
1858 18212 : call set_eps_nuc(Qtotal_rni56ec - Qneu_rni56ec, [y(ini56)], irn56ec, i_ni56_co56)
1859 :
1860 9106 : if (plus_co56) then
1861 8 : call set_eps_nuc(Qtotal_rco56ec - Qneu_rco56ec, [y(ico56)], irco56ec, i_co56_fe56)
1862 : else
1863 18204 : call set_eps_nuc(Qtotal_rco56ec - Qneu_rco56ec, [y(ini56)], irn56ec, i_co56_fe56)
1864 :
1865 : end if
1866 :
1867 227650 : eps_nuc_cat = eps_nuc_cat * Qconv
1868 1065406 : n% eps_nuc_rate = n% eps_nuc_rate * Qconv
1869 :
1870 227650 : do i=1,num_categories
1871 227650 : eps_nuc_categories(i) = eps_nuc_cat(i)
1872 : end do
1873 :
1874 : ! check eps_nuc vs sum(eps_nuc_cat)
1875 :
1876 9106 : sum_categories_q = sum(eps_nuc_cat)
1877 9106 : sum_categories = sum_categories_q
1878 :
1879 : if (.false. .and. &
1880 : abs(eps_nuc) > 1d-10*abs(eps_nuc_cat(iphoto)) .and. abs(eps_nuc) > 1d0 .and. &
1881 9106 : abs(sum_categories - eps_nuc) > 1d-2*min(abs(sum_categories),abs(eps_nuc))) then
1882 : !$OMP critical (net21_crit1)
1883 : !write(*,*) '>>>>> problem in net_approx21_procs.inc, approx21_eps_info <<<<<<'
1884 : !write(*,*) ' please report it. can edit the file in eos/private to remove this test. '
1885 : write(*,1) 'logT', n% logT
1886 : write(*,1) 'approx21_eps_info'
1887 : write(*,'(A)')
1888 : do i=1,num_categories
1889 : if (abs(eps_nuc_categories(i)) > 1d-6) then
1890 : write(*,1) trim(category_name(i)), eps_nuc_categories(i)
1891 : end if
1892 : end do
1893 : write(*,*)
1894 : write(*,1) 'eps_total', eps_total
1895 : write(*,1) 'eps_neu', eps_neu
1896 : write(*,1) 'eps_nuc', eps_nuc
1897 : write(*,1) 'sum(cat)', sum_categories_q
1898 : write(*,1) 'sum(cat) - eps_nuc', sum_categories_q - eps_nuc_q
1899 : write(*,1) 'sum(cat)/eps_nuc - 1', (sum_categories_q - eps_nuc_q)/eps_nuc_q
1900 : write(*,'(A)')
1901 : call mesa_error(__FILE__,__LINE__)
1902 : !$OMP end critical (net21_crit1)
1903 : end if
1904 :
1905 : ! for debugging use reduced_net_for_testing
1906 :
1907 : if (.false. .and. n% logT >= 9.220336900d0 .and. n% logT <= 9.2203369009d0 .and. &
1908 : abs(eps_nuc) > 1d-10*abs(eps_nuc_cat(iphoto)) .and. abs(eps_nuc) > 1d0 .and. &
1909 21194 : abs(sum_categories - eps_nuc) > 1d-2*min(abs(sum_categories),abs(eps_nuc))) then
1910 : write(*,1) 'logT', n% logT
1911 : write(*,2) 'icrx ' // trim(chem_isos% name(n% g% approx21_ye_iso)), icrx
1912 : write(*,2) 'fe56ec_n_neut', fe56ec_n_neut
1913 : write(*,1) 'fe56ec_fake_factor', fe56ec_fake_factor
1914 : write(*,'(A)')
1915 : do i=1,num_categories
1916 : write(*,1) trim(category_name(i)), eps_nuc_categories(i)
1917 : end do
1918 : write(*,*)
1919 : write(*,1) 'eps_total', eps_total
1920 : write(*,1) 'eps_neu', eps_neu
1921 : write(*,1) 'eps_nuc', eps_nuc
1922 : write(*,1) 'sum(cat)', sum_categories
1923 : write(*,1) 'sum(cat) - eps_nuc', sum_categories_q - eps_nuc_q
1924 : write(*,1) 'sum(cat)/eps_nuc - 1', (sum_categories_q - eps_nuc_q)/eps_nuc_q
1925 : call mesa_error(__FILE__,__LINE__,'approx21_eps_info')
1926 : end if
1927 :
1928 :
1929 : contains
1930 :
1931 : real(qp) function m2(a1,a2)
1932 : real(dp), intent(in) :: a1, a2
1933 : real(qp) :: q1, q2
1934 : q1 = a1; q2 = a2; m2 = q1*q2
1935 21194 : end function m2
1936 :
1937 60600 : real(qp) function m3(a1,a2,a3)
1938 : real(dp), intent(in) :: a1, a2, a3
1939 60600 : real(qp) :: q1, q2, q3
1940 30300 : q1 = a1; q2 = a2; q3 = a3; m3 = q1*q2*q3
1941 : end function m3
1942 :
1943 30300 : real(qp) function m4(a1,a2,a3,a4)
1944 : real(dp), intent(in) :: a1, a2, a3, a4
1945 30300 : real(qp) :: q1, q2, q3, q4
1946 30300 : q1 = a1; q2 = a2; q3 = a3; q4 = a4; m4 = q1*q2*q3*q4
1947 : end function m4
1948 :
1949 30300 : real(qp) function m5(a1,a2,a3,a4,a5)
1950 : real(dp), intent(in) :: a1, a2, a3, a4, a5
1951 30300 : real(qp) :: q1, q2, q3, q4, q5
1952 30300 : q1 = a1; q2 = a2; q3 = a3; q4 = a4; q5 = a5; m5 = q1*q2*q3*q4*q5
1953 : end function m5
1954 :
1955 :
1956 637420 : subroutine set_eps_nuc(q, ys, rat_id, eps_id)
1957 : real(dp),intent(in) :: q
1958 : real(dp),intent(in) :: ys(:)
1959 : integer :: rat_id, eps_id
1960 637420 : real(qp) :: val1,val2,val3
1961 :
1962 1939578 : val1 = product(real(ys,kind=qp))
1963 637420 : val2 = val1 * real(rate(rat_id),kind=qp)
1964 :
1965 637420 : val3 = real(q,kind=qp) * val2
1966 :
1967 637420 : eps_nuc_cat(eps_id) = eps_nuc_cat(eps_id) + val3
1968 :
1969 637420 : if(rat_id < ifa) then
1970 509936 : n% raw_rate(rat_id) = n% raw_rate(rat_id) + val1 * real(n% rate_raw(rat_id),kind=qp)
1971 509936 : n% screened_rate(rat_id) = n% screened_rate(rat_id) + val2
1972 509936 : n% eps_nuc_rate(rat_id) = n% eps_nuc_rate(rat_id) + val3
1973 : else
1974 127484 : n% raw_rate(rat_id) = 0d0
1975 127484 : n% screened_rate(rat_id) = 0d0
1976 127484 : n% eps_nuc_rate(rat_id) = 0d0
1977 : end if
1978 :
1979 637420 : end subroutine set_eps_nuc
1980 :
1981 :
1982 : end subroutine approx21_eps_info
1983 :
1984 :
1985 6044 : subroutine approx21_d_epsneu_dy( &
1986 6044 : y, rate, &
1987 : Qneu_rpp, Qneu_rpp2, Qneu_rpp3, &
1988 : Qneu_rcpg, Qneu_rnpg, Qneu_ropg, &
1989 6044 : d_epsneu_dy, plus_co56, ierr)
1990 : use const_def, only: Qconv
1991 : real(dp), dimension(:), intent(in) :: y, rate
1992 : real(dp), intent(in) :: &
1993 : Qneu_rpp, Qneu_rpp2, Qneu_rpp3, &
1994 : Qneu_rcpg, Qneu_rnpg, Qneu_ropg
1995 : real(dp), intent(inout) :: d_epsneu_dy(:)
1996 : logical, intent(in) :: plus_co56
1997 : integer, intent(out) :: ierr
1998 :
1999 6044 : real(dp) :: fII
2000 :
2001 6044 : ierr = 0
2002 :
2003 6044 : fII = 0.5d0 ! fix this
2004 :
2005 132972 : d_epsneu_dy(1:species(plus_co56)) = 0d0
2006 :
2007 : d_epsneu_dy(ih1) = Qconv*( &
2008 : Qneu_rpp * y(ih1) * rate(irpp) + & ! rpp_to_he3
2009 : Qneu_rcpg * y(ic12) * rate(ircpg) + & ! C of CNO
2010 : Qneu_rnpg * y(in14) * rate(irnpg) + & ! N of CNO
2011 6044 : Qneu_ropg * y(io16) * rate(iropg)) ! O of CNO
2012 :
2013 : d_epsneu_dy(ihe3) = Qconv*( &
2014 : Qneu_rpp2 * y(ihe4) * rate(irhe3ag) * fII + & ! r34_pp2
2015 6044 : Qneu_rpp3 * y(ihe4) * rate(irhe3ag) * (1d0-fII)) ! r34_pp3
2016 :
2017 : d_epsneu_dy(ihe4) = Qconv*( &
2018 : Qneu_rpp2 * y(ihe3) * rate(irhe3ag) * fII + & ! r34_pp2
2019 6044 : Qneu_rpp3 * y(ihe3) * rate(irhe3ag) * (1d0-fII)) ! r34_pp3
2020 :
2021 : d_epsneu_dy(ic12) = Qconv* &
2022 6044 : Qneu_rcpg * y(ih1) * rate(ircpg) ! C of CNO
2023 :
2024 : d_epsneu_dy(in14) = Qconv* &
2025 6044 : Qneu_rnpg * y(ih1) * rate(irnpg) ! N of CNO
2026 :
2027 : d_epsneu_dy(io16) = Qconv* &
2028 6044 : Qneu_ropg * y(ih1) * rate(iropg) ! O of CNO
2029 :
2030 6044 : end subroutine approx21_d_epsneu_dy
2031 :
2032 :
2033 6044 : subroutine approx21_dfdy( &
2034 6044 : y, dfdy, fe56ec_fake_factor_in, min_T, fe56ec_n_neut, &
2035 6044 : ratdum, dratdumdt, dratdumdd, dratdumdy1, dratdumdy2, btemp, plus_co56, ierr)
2036 : real(dp), intent(in) :: fe56ec_fake_factor_in, min_T, btemp
2037 : integer, intent(in) :: fe56ec_n_neut
2038 : real(dp), intent(in), dimension(:) :: &
2039 : y, ratdum, dratdumdt, dratdumdd, dratdumdy1, dratdumdy2
2040 : real(dp), intent(inout) :: dfdy(:,:)
2041 : logical, intent(in) :: plus_co56
2042 : integer, intent(out) :: ierr
2043 :
2044 : integer :: i,j
2045 : real(dp) :: abar,zbar,ye,taud,taut, b1, &
2046 : snuda,snudz,enuc,velx,posx,zz
2047 6044 : real(dp) :: fe56ec_fake_factor
2048 :
2049 6044 : ierr = 0
2050 :
2051 : ! Turn on special fe56ec rate above some temperature
2052 6044 : fe56ec_fake_factor=eval_fe56ec_fake_factor(fe56ec_fake_factor_in,min_T,btemp)
2053 :
2054 : ! NOTE: use of quad precision for dfdy doesn't make a difference.
2055 :
2056 2798548 : dfdy(1:species(plus_co56),1:species(plus_co56)) = 0.0d0
2057 :
2058 : ! h1 jacobian elements
2059 : dfdy(ih1,ih1) = -3.0d0 * y(ih1) * ratdum(irpp) &
2060 : - 2.0d0 * y(ic12) * ratdum(ircpg) &
2061 : - 2.0d0 * y(in14) * ratdum(irnpg) &
2062 : - 2.0d0 * y(in14) * y(ih1) * dratdumdy1(irnpg) &
2063 : - 2.0d0 * y(io16) * ratdum(iropg) &
2064 : - 2.0d0 * y(io16) * y(ih1) * dratdumdy1(iropg) &
2065 6044 : - 3.0d0 * ratdum(irpen)
2066 :
2067 : dfdy(ih1,ihe3) = 2.0d0 * y(ihe3) * ratdum(ir33) &
2068 6044 : - y(ihe4) * ratdum(irhe3ag)
2069 :
2070 : dfdy(ih1,ihe4) = -y(ihe3) * ratdum(irhe3ag) &
2071 6044 : - y(ihe3) * y(ihe4) * dratdumdy1(irhe3ag)
2072 :
2073 6044 : dfdy(ih1,ic12) = -2.0d0 * y(ih1) * ratdum(ircpg)
2074 :
2075 6044 : dfdy(ih1,in14) = -2.0d0 * y(ih1) * ratdum(irnpg)
2076 :
2077 6044 : dfdy(ih1,io16) = -2.0d0 * y(ih1) * ratdum(iropg)
2078 :
2079 :
2080 : ! he3 jacobian elements
2081 : dfdy(ihe3,ih1) = y(ih1) * ratdum(irpp) &
2082 6044 : + ratdum(irpen)
2083 :
2084 : dfdy(ihe3,ihe3) = -2.0d0 * y(ihe3) * ratdum(ir33) &
2085 6044 : - y(ihe4) * ratdum(irhe3ag)
2086 :
2087 : dfdy(ihe3,ihe4) = -y(ihe3) * ratdum(irhe3ag) &
2088 6044 : - y(ihe3) * y(ihe4) * dratdumdy1(irhe3ag)
2089 :
2090 :
2091 : ! he4 jacobian elements
2092 : dfdy(ihe4,ih1) = y(in14) * ratdum(ifa) * ratdum(irnpg) &
2093 : + y(in14) * y(ih1) * ratdum(ifa) * dratdumdy1(irnpg) &
2094 : + y(io16) * ratdum(iropg) &
2095 6044 : + y(io16) * y(ih1) * dratdumdy1(iropg)
2096 :
2097 : dfdy(ihe4,ihe3) = y(ihe3) * ratdum(ir33) &
2098 6044 : + y(ihe4) * ratdum(irhe3ag)
2099 :
2100 :
2101 : dfdy(ihe4,ihe4) = -1.5d0 * y(ihe4) * y(ihe4) * ratdum(ir3a) &
2102 : - y(ic12) * ratdum(ircag) &
2103 : - y(io16) * ratdum(iroag) &
2104 : - y(ine20) * ratdum(irneag) &
2105 : - y(img24) * ratdum(irmgag) &
2106 : - y(isi28) * ratdum(irsiag) &
2107 : - y(is32) * ratdum(irsag) &
2108 : - y(iar36) * ratdum(irarag) &
2109 : - y(ica40) * ratdum(ircaag) &
2110 : - y(iti44) * ratdum(irtiag) &
2111 : - y(icr48) * ratdum(ircrag) &
2112 6044 : - y(ife52) * ratdum(irfeag)
2113 :
2114 : dfdy(ihe4,ihe4) = dfdy(ihe4,ihe4) &
2115 : - y(img24) * ratdum(irmgap) * (1.0d0-ratdum(irr1)) &
2116 : - y(isi28) * ratdum(irsiap) * (1.0d0-ratdum(irs1)) &
2117 : - y(is32) * ratdum(irsap) * (1.0d0-ratdum(irt1)) &
2118 : - y(iar36) * ratdum(irarap) * (1.0d0-ratdum(iru1)) &
2119 : - y(ica40) * ratdum(ircaap) * (1.0d0-ratdum(irv1)) &
2120 : - y(iti44) * ratdum(irtiap) * (1.0d0-ratdum(irw1)) &
2121 6044 : - y(icr48) * ratdum(ircrap) * (1.0d0-ratdum(irx1))
2122 :
2123 : dfdy(ihe4,ihe4) = dfdy(ihe4,ihe4) &
2124 : - y(ife52) * ratdum(ir6f54) &
2125 : - y(ife52) * y(iprot) * ratdum(ir7f54) &
2126 : - ratdum(iralf1) &
2127 6044 : - y(ife54) * ratdum(irfe56_aux4)
2128 :
2129 :
2130 : dfdy(ihe4,ihe4) = dfdy(ihe4,ihe4) &
2131 : + y(ihe3) * ratdum(irhe3ag) &
2132 : + y(ihe3) * y(ihe4) * dratdumdy1(irhe3ag) &
2133 6044 : - y(in14) * ratdum(irnag) * 1.5d0
2134 :
2135 :
2136 : dfdy(ihe4,ic12) = y(ic12) * ratdum(ir1212) &
2137 : + 0.5d0 * y(io16) * ratdum(ir1216) &
2138 : + 3.0d0 * ratdum(irg3a) &
2139 6044 : - y(ihe4) * ratdum(ircag)
2140 :
2141 :
2142 : dfdy(ihe4,in14) = y(ih1) * ratdum(ifa) * ratdum(irnpg) &
2143 6044 : - y(ihe4) * ratdum(irnag) * 1.5d0
2144 :
2145 :
2146 : dfdy(ihe4,io16) = 0.5d0 * y(ic12) * ratdum(ir1216) &
2147 : + 1.12d0 * 0.5d0*y(io16) * ratdum(ir1616) &
2148 : + 0.68d0 * ratdum(irs1) * 0.5d0*y(io16) * ratdum(ir1616) &
2149 : + ratdum(iroga) &
2150 : - y(ihe4) * ratdum(iroag) &
2151 6044 : + y(ih1) * ratdum(iropg)
2152 :
2153 : dfdy(ihe4,ine20) = ratdum(irnega) &
2154 6044 : - y(ihe4) * ratdum(irneag)
2155 :
2156 : dfdy(ihe4,img24) = ratdum(irmgga) &
2157 : - y(ihe4) * ratdum(irmgag) &
2158 6044 : - y(ihe4) * ratdum(irmgap) * (1.0d0-ratdum(irr1))
2159 :
2160 : dfdy(ihe4,isi28) = ratdum(irsiga) &
2161 : - y(ihe4) * ratdum(irsiag) &
2162 : - y(ihe4) * ratdum(irsiap) * (1.0d0-ratdum(irs1)) &
2163 6044 : + ratdum(irr1) * ratdum(irsigp)
2164 :
2165 : dfdy(ihe4,is32) = ratdum(irsga) &
2166 : - y(ihe4) * ratdum(irsag) &
2167 : - y(ihe4) * ratdum(irsap) * (1.0d0-ratdum(irt1)) &
2168 6044 : + ratdum(irs1) * ratdum(irsgp)
2169 :
2170 : dfdy(ihe4,iar36) = ratdum(irarga) &
2171 : - y(ihe4) * ratdum(irarag) &
2172 : - y(ihe4) * ratdum(irarap) * (1.0d0-ratdum(iru1)) &
2173 6044 : + ratdum(irt1) * ratdum(irargp)
2174 :
2175 : dfdy(ihe4,ica40) = ratdum(ircaga) &
2176 : - y(ihe4) * ratdum(ircaag) &
2177 : - y(ihe4) * ratdum(ircaap) * (1.0d0-ratdum(irv1)) &
2178 6044 : + ratdum(iru1) * ratdum(ircagp)
2179 :
2180 : dfdy(ihe4,iti44) = ratdum(irtiga) &
2181 : - y(ihe4) * ratdum(irtiag) &
2182 : - y(ihe4) * ratdum(irtiap) * (1.0d0-ratdum(irw1)) &
2183 6044 : + ratdum(irv1) * ratdum(irtigp)
2184 :
2185 : dfdy(ihe4,icr48) = ratdum(ircrga) &
2186 : - y(ihe4) * ratdum(ircrag) &
2187 : - y(ihe4) * ratdum(ircrap) * (1.0d0-ratdum(irx1)) &
2188 6044 : + ratdum(irw1) * ratdum(ircrgp)
2189 :
2190 : dfdy(ihe4,ife52) = ratdum(irfega) &
2191 : - y(ihe4) * ratdum(irfeag) &
2192 : + ratdum(irx1) * ratdum(irfegp) &
2193 : - y(ihe4) * ratdum(ir6f54) &
2194 6044 : - y(ihe4) * y(iprot) * ratdum(ir7f54)
2195 :
2196 : dfdy(ihe4,ife54) = y(iprot) * y(iprot) * ratdum(ir5f54) &
2197 6044 : - y(ihe4) * ratdum(irfe56_aux4)
2198 :
2199 6044 : dfdy(ihe4,ife56) = y(iprot) * y(iprot) * ratdum(irfe56_aux3)
2200 :
2201 : dfdy(ihe4,ini56) = ratdum(irniga) &
2202 6044 : + y(iprot) * ratdum(ir8f54)
2203 :
2204 :
2205 : dfdy(ihe4,ineut) = -y(ihe4) * dratdumdy1(iralf1) &
2206 : + 2.0d0 * y(ineut) * y(iprot)*y(iprot) * ratdum(iralf2) &
2207 6044 : + y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy1(iralf2)
2208 :
2209 : include 'formats'
2210 :
2211 : dfdy(ihe4,iprot) = 2.0d0 * y(ife54) * y(iprot) * ratdum(ir5f54) &
2212 : + y(ife54) * y(iprot) * y(iprot) * dratdumdy1(ir5f54) &
2213 : - y(ihe4) * y(ife52) * dratdumdy1(ir6f54) &
2214 : - y(ife52) * y(ihe4) * ratdum(ir7f54) &
2215 : - y(ife52) * y(ihe4) * y(iprot) * dratdumdy1(ir7f54) &
2216 : + y(ini56) * ratdum(ir8f54) &
2217 : + y(ini56) * y(iprot) * dratdumdy1(ir8f54) &
2218 : - y(ihe4) * dratdumdy2(iralf1) &
2219 : + 2.0d0 * y(ineut)*y(ineut) * y(iprot) * ratdum(iralf2) &
2220 : + y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy2(iralf2) &
2221 : + 2.0d0 * y(ife56) * y(iprot) * ratdum(irfe56_aux3) &
2222 : + y(ife56) * y(iprot) * y(iprot) * dratdumdy1(irfe56_aux3) &
2223 6044 : - y(ihe4) * y(ife54) * dratdumdy1(irfe56_aux4)
2224 :
2225 :
2226 :
2227 : ! c12 jacobian elements
2228 : dfdy(ic12,ih1) = -y(ic12) * ratdum(ircpg) &
2229 : + y(in14) * ratdum(ifa) * ratdum(irnpg) &
2230 6044 : + y(in14) * y(ih1) * ratdum(ifa) * dratdumdy1(irnpg)
2231 :
2232 : dfdy(ic12,ihe4) = 0.5d0 * y(ihe4) * y(ihe4) * ratdum(ir3a) &
2233 6044 : - y(ic12) * ratdum(ircag)
2234 :
2235 : dfdy(ic12,ic12) = -2.0d0 * y(ic12) * ratdum(ir1212) &
2236 : - y(io16) * ratdum(ir1216) &
2237 : - ratdum(irg3a) &
2238 : - y(ihe4) * ratdum(ircag) &
2239 6044 : - y(ih1) * ratdum(ircpg)
2240 :
2241 6044 : dfdy(ic12,in14) = y(ih1) * ratdum(ifa) * ratdum(irnpg)
2242 :
2243 : dfdy(ic12,io16) = -y(ic12) * ratdum(ir1216) &
2244 6044 : + ratdum(iroga)
2245 :
2246 :
2247 : ! n14 jacobian elements
2248 : dfdy(in14,ih1) = y(ic12) * ratdum(ircpg) &
2249 : - y(in14) * ratdum(irnpg) &
2250 : - y(in14) * y(ih1) * dratdumdy1(irnpg) &
2251 : + y(io16) * ratdum(iropg) &
2252 6044 : + y(io16) * y(ih1) * dratdumdy1(iropg)
2253 :
2254 6044 : dfdy(in14,ihe4) = -y(in14) * ratdum(irnag)
2255 :
2256 6044 : dfdy(in14,ic12) = y(ih1) * ratdum(ircpg)
2257 :
2258 : dfdy(in14,in14) = -y(ih1) * ratdum(irnpg) &
2259 6044 : - y(ihe4) * ratdum(irnag)
2260 :
2261 6044 : dfdy(in14,io16) = y(ih1) * ratdum(iropg)
2262 :
2263 :
2264 :
2265 : ! o16 jacobian elements
2266 : dfdy(io16,ih1) = y(in14) * ratdum(ifg) * ratdum(irnpg) &
2267 : + y(in14) * y(ih1) * ratdum(ifg) * dratdumdy1(irnpg) &
2268 : - y(io16) * ratdum(iropg) &
2269 6044 : - y(io16) * y(ih1) * dratdumdy1(iropg)
2270 :
2271 : dfdy(io16,ihe4) = y(ic12)*ratdum(ircag) &
2272 6044 : - y(io16)*ratdum(iroag)
2273 :
2274 : dfdy(io16,ic12) = -y(io16)*ratdum(ir1216) &
2275 6044 : + y(ihe4)*ratdum(ircag)
2276 :
2277 6044 : dfdy(io16,in14) = y(ih1) * ratdum(ifg) * ratdum(irnpg)
2278 :
2279 : dfdy(io16,io16) = - y(ic12) * ratdum(ir1216) &
2280 : - 2.0d0 * y(io16) * ratdum(ir1616) &
2281 : - y(ihe4) * ratdum(iroag) &
2282 : - ratdum(iroga) &
2283 6044 : - y(ih1) * ratdum(iropg)
2284 :
2285 6044 : dfdy(io16,ine20) = ratdum(irnega)
2286 :
2287 :
2288 : ! ne20 jacobian elements
2289 : dfdy(ine20,ihe4) = y(io16) * ratdum(iroag) &
2290 : - y(ine20) * ratdum(irneag) &
2291 6044 : + y(in14) * ratdum(irnag)
2292 :
2293 6044 : dfdy(ine20,ic12) = y(ic12) * ratdum(ir1212)
2294 :
2295 6044 : dfdy(ine20,in14) = y(ihe4) * ratdum(irnag)
2296 :
2297 6044 : dfdy(ine20,io16) = y(ihe4) * ratdum(iroag)
2298 :
2299 : dfdy(ine20,ine20) = -y(ihe4) * ratdum(irneag) &
2300 6044 : - ratdum(irnega)
2301 :
2302 6044 : dfdy(ine20,img24) = ratdum(irmgga)
2303 :
2304 :
2305 :
2306 : ! mg24 jacobian elements
2307 : dfdy(img24,ihe4) = y(ine20) * ratdum(irneag) &
2308 : -y(img24) * ratdum(irmgag) &
2309 6044 : -y(img24) * ratdum(irmgap) * (1.0d0-ratdum(irr1))
2310 :
2311 6044 : dfdy(img24,ic12) = 0.5d0 * y(io16) * ratdum(ir1216)
2312 :
2313 6044 : dfdy(img24,io16) = 0.5d0 * y(ic12) * ratdum(ir1216)
2314 :
2315 6044 : dfdy(img24,ine20) = y(ihe4) * ratdum(irneag)
2316 :
2317 : dfdy(img24,img24) = -y(ihe4) * ratdum(irmgag) &
2318 : - ratdum(irmgga) &
2319 6044 : - y(ihe4)*ratdum(irmgap)*(1.0d0-ratdum(irr1))
2320 :
2321 : dfdy(img24,isi28) = ratdum(irsiga) &
2322 6044 : + ratdum(irr1) * ratdum(irsigp)
2323 :
2324 :
2325 :
2326 : ! si28 jacobian elements
2327 : dfdy(isi28,ihe4) = y(img24) * ratdum(irmgag) &
2328 : - y(isi28) * ratdum(irsiag) &
2329 : + y(img24) * ratdum(irmgap) * (1.0d0-ratdum(irr1)) &
2330 6044 : - y(isi28) * ratdum(irsiap) * (1.0d0-ratdum(irs1))
2331 :
2332 6044 : dfdy(isi28,ic12) = 0.5d0 * y(io16) * ratdum(ir1216)
2333 :
2334 : dfdy(isi28,io16) = 0.5d0 * y(ic12) * ratdum(ir1216) &
2335 : + 1.12d0 * 0.5d0*y(io16) * ratdum(ir1616) &
2336 6044 : + 0.68d0 * 0.5d0*y(io16) * ratdum(irs1) * ratdum(ir1616)
2337 :
2338 : dfdy(isi28,img24) = y(ihe4) * ratdum(irmgag) &
2339 6044 : + y(ihe4) * ratdum(irmgap) * (1.0d0-ratdum(irr1))
2340 :
2341 : dfdy(isi28,isi28) = -y(ihe4) * ratdum(irsiag) &
2342 : - ratdum(irsiga) &
2343 : - ratdum(irr1) * ratdum(irsigp) &
2344 6044 : - y(ihe4) * ratdum(irsiap) * (1.0d0-ratdum(irs1))
2345 :
2346 : dfdy(isi28,is32) = ratdum(irsga) &
2347 6044 : + ratdum(irs1) * ratdum(irsgp)
2348 :
2349 :
2350 :
2351 : ! s32 jacobian elements
2352 : dfdy(is32,ihe4) = y(isi28) * ratdum(irsiag) &
2353 : - y(is32) * ratdum(irsag) &
2354 : + y(isi28) * ratdum(irsiap) * (1.0d0-ratdum(irs1)) &
2355 6044 : - y(is32) * ratdum(irsap) * (1.0d0-ratdum(irt1))
2356 :
2357 : dfdy(is32,io16) = &
2358 : + 0.68d0*0.5d0*y(io16)*ratdum(ir1616)*(1.0d0-ratdum(irs1)) &
2359 6044 : + 0.2d0 * 0.5d0*y(io16) * ratdum(ir1616)
2360 :
2361 : dfdy(is32,isi28) = y(ihe4) * ratdum(irsiag) &
2362 6044 : + y(ihe4) * ratdum(irsiap) * (1.0d0-ratdum(irs1))
2363 :
2364 : dfdy(is32,is32) = -y(ihe4) * ratdum(irsag) &
2365 : - ratdum(irsga) &
2366 : - ratdum(irs1) * ratdum(irsgp) &
2367 6044 : - y(ihe4) * ratdum(irsap) * (1.0d0-ratdum(irt1))
2368 :
2369 : dfdy(is32,iar36) = ratdum(irarga) &
2370 6044 : + ratdum(irt1) * ratdum(irargp)
2371 :
2372 :
2373 :
2374 : ! ar36 jacobian elements
2375 : dfdy(iar36,ihe4) = y(is32) * ratdum(irsag) &
2376 : - y(iar36) * ratdum(irarag) &
2377 : + y(is32) * ratdum(irsap) * (1.0d0-ratdum(irt1)) &
2378 6044 : - y(iar36) * ratdum(irarap) * (1.0d0-ratdum(iru1))
2379 :
2380 : dfdy(iar36,is32) = y(ihe4) * ratdum(irsag) &
2381 6044 : + y(ihe4) * ratdum(irsap) * (1.0d0-ratdum(irt1))
2382 :
2383 : dfdy(iar36,iar36) = -y(ihe4) * ratdum(irarag) &
2384 : - ratdum(irarga) &
2385 : - ratdum(irt1) * ratdum(irargp) &
2386 6044 : - y(ihe4) * ratdum(irarap) * (1.0d0-ratdum(iru1))
2387 :
2388 : dfdy(iar36,ica40) = ratdum(ircaga) &
2389 6044 : + ratdum(ircagp) * ratdum(iru1)
2390 :
2391 :
2392 :
2393 : ! ca40 jacobian elements
2394 : dfdy(ica40,ihe4) = y(iar36) * ratdum(irarag) &
2395 : - y(ica40) * ratdum(ircaag) &
2396 : + y(iar36) * ratdum(irarap)*(1.0d0-ratdum(iru1)) &
2397 6044 : - y(ica40) * ratdum(ircaap)*(1.0d0-ratdum(irv1))
2398 :
2399 : dfdy(ica40,iar36) = y(ihe4) * ratdum(irarag) &
2400 6044 : + y(ihe4) * ratdum(irarap)*(1.0d0-ratdum(iru1))
2401 :
2402 : dfdy(ica40,ica40) = -y(ihe4) * ratdum(ircaag) &
2403 : - ratdum(ircaga) &
2404 : - ratdum(ircagp) * ratdum(iru1) &
2405 6044 : - y(ihe4) * ratdum(ircaap)*(1.0d0-ratdum(irv1))
2406 :
2407 : dfdy(ica40,iti44) = ratdum(irtiga) &
2408 6044 : + ratdum(irtigp) * ratdum(irv1)
2409 :
2410 :
2411 :
2412 : ! ti44 jacobian elements
2413 : dfdy(iti44,ihe4) = y(ica40) * ratdum(ircaag) &
2414 : - y(iti44) * ratdum(irtiag) &
2415 : + y(ica40) * ratdum(ircaap)*(1.0d0-ratdum(irv1)) &
2416 6044 : - y(iti44) * ratdum(irtiap)*(1.0d0-ratdum(irw1))
2417 :
2418 : dfdy(iti44,ica40) = y(ihe4) * ratdum(ircaag) &
2419 6044 : + y(ihe4) * ratdum(ircaap)*(1.0d0-ratdum(irv1))
2420 :
2421 : dfdy(iti44,iti44) = -y(ihe4) * ratdum(irtiag) &
2422 : - ratdum(irtiga) &
2423 : - ratdum(irv1) * ratdum(irtigp) &
2424 6044 : - y(ihe4) * ratdum(irtiap)*(1.0d0-ratdum(irw1))
2425 :
2426 : dfdy(iti44,icr48) = ratdum(ircrga) &
2427 6044 : + ratdum(irw1) * ratdum(ircrgp)
2428 :
2429 :
2430 :
2431 : ! cr48 jacobian elements
2432 : dfdy(icr48,ihe4) = y(iti44) * ratdum(irtiag) &
2433 : - y(icr48) * ratdum(ircrag) &
2434 : + y(iti44) * ratdum(irtiap)*(1.0d0-ratdum(irw1)) &
2435 6044 : - y(icr48) * ratdum(ircrap)*(1.0d0-ratdum(irx1))
2436 :
2437 : dfdy(icr48,iti44) = y(ihe4) * ratdum(irtiag) &
2438 6044 : + y(ihe4) * ratdum(irtiap)*(1.0d0-ratdum(irw1))
2439 :
2440 : dfdy(icr48,icr48) = -y(ihe4) * ratdum(ircrag) &
2441 : - ratdum(ircrga) &
2442 : - ratdum(irw1) * ratdum(ircrgp) &
2443 6044 : - y(ihe4) * ratdum(ircrap)*(1.0d0-ratdum(irx1))
2444 :
2445 : dfdy(icr48,ife52) = ratdum(irfega) &
2446 6044 : + ratdum(irx1) * ratdum(irfegp)
2447 :
2448 :
2449 : ! crx jacobian elements
2450 6044 : dfdy(icrx,ife56) = fe56ec_fake_factor * ratdum(irn56ec)
2451 :
2452 :
2453 : ! fe52 jacobian elements
2454 : dfdy(ife52,ihe4) = y(icr48) * ratdum(ircrag) &
2455 : - y(ife52) * ratdum(irfeag) &
2456 : + y(icr48) * ratdum(ircrap) * (1.0d0-ratdum(irx1)) &
2457 : - y(ife52) * ratdum(ir6f54) &
2458 6044 : - y(ife52) * y(iprot) * ratdum(ir7f54)
2459 :
2460 : dfdy(ife52,icr48) = y(ihe4) * ratdum(ircrag) &
2461 6044 : + y(ihe4) * ratdum(ircrap) * (1.0d0-ratdum(irx1))
2462 :
2463 : dfdy(ife52,ife52) = - y(ihe4) * ratdum(irfeag) &
2464 : - ratdum(irfega) &
2465 : - ratdum(irx1) * ratdum(irfegp) &
2466 : - y(ineut) * y(ineut) * ratdum(ir2f54) &
2467 : - y(ihe4) * ratdum(ir6f54) &
2468 6044 : - y(ihe4) * y(iprot) * ratdum(ir7f54)
2469 :
2470 : dfdy(ife52,ife54) = ratdum(ir1f54) + &
2471 6044 : y(iprot) * y(iprot) * ratdum(ir5f54)
2472 :
2473 : dfdy(ife52,ini56) = ratdum(irniga) &
2474 6044 : + y(iprot) * ratdum(ir8f54)
2475 :
2476 : dfdy(ife52,ineut) = &
2477 : y(ife54) * dratdumdy1(ir1f54) &
2478 : - 2.0d0 * y(ife52) * y(ineut) * ratdum(ir2f54) &
2479 6044 : - y(ife52) * y(ineut) * y(ineut) * dratdumdy1(ir2f54)
2480 :
2481 : dfdy(ife52,iprot) = 2.0d0 * y(ife54) * y(iprot) * ratdum(ir5f54) &
2482 : + y(ife54) * y(iprot) * y(iprot) * dratdumdy1(ir5f54) &
2483 : - y(ihe4) * y(ife52) * dratdumdy1(ir6f54) &
2484 : - y(ife52) * y(ihe4) * ratdum(ir7f54) &
2485 : - y(ife52) * y(ihe4) * y(iprot) * dratdumdy1(ir7f54) &
2486 : + y(ini56) * ratdum(ir8f54) &
2487 6044 : + y(ini56) * y(iprot) * dratdumdy1(ir8f54)
2488 :
2489 :
2490 : ! fe54 jacobian elements
2491 : dfdy(ife54,ihe4) = y(ife52) * ratdum(ir6f54) &
2492 6044 : - y(ife54) * ratdum(irfe56_aux4)
2493 :
2494 : dfdy(ife54,ife52) = &
2495 : y(ineut) * y(ineut) * ratdum(ir2f54) + &
2496 6044 : y(ihe4) * ratdum(ir6f54)
2497 :
2498 : dfdy(ife54,ife54) = &
2499 : - ratdum(ir1f54) &
2500 : - y(ineut) * y(ineut) * ratdum(irfe56_aux2) &
2501 : - y(iprot) * y(iprot) * ratdum(ir3f54) &
2502 : - y(iprot) * y(iprot) * ratdum(ir5f54) &
2503 6044 : - y(ihe4) * ratdum(irfe56_aux4)
2504 :
2505 : dfdy(ife54,ife56) = &
2506 : ratdum(irfe56_aux1) + &
2507 6044 : y(iprot) * y(iprot) * ratdum(irfe56_aux3)
2508 :
2509 6044 : dfdy(ife54,ini56) = ratdum(ir4f54)
2510 :
2511 : dfdy(ife54,ineut) = &
2512 : - y(ife54) * dratdumdy1(ir1f54) &
2513 : + 2.0d0 * y(ife52) * y(ineut) * ratdum(ir2f54) &
2514 : + y(ife52) * y(ineut) * y(ineut) * dratdumdy1(ir2f54) &
2515 : + y(ife56) * dratdumdy1(irfe56_aux1) &
2516 : - 2.0d0 * y(ife54) * y(ineut) * ratdum(irfe56_aux2) &
2517 6044 : - y(ife54) * y(ineut) * y(ineut) * dratdumdy1(irfe56_aux2)
2518 :
2519 : dfdy(ife54,iprot) = -2.0d0 * y(ife54) * y(iprot) * ratdum(ir3f54) &
2520 : - y(ife54) * y(iprot) * y(iprot) * dratdumdy1(ir3f54) &
2521 : + y(ini56) * dratdumdy1(ir4f54) &
2522 : - 2.0d0 * y(ife54) * y(iprot) * ratdum(ir5f54) &
2523 : - y(ife54) * y(iprot) * y(iprot) * dratdumdy1(ir5f54) &
2524 : + y(ihe4) * y(ife52) * dratdumdy1(ir6f54) &
2525 : + 2.0d0 * y(ife56) * y(iprot) * ratdum(irfe56_aux3) &
2526 : + y(ife56) * y(iprot) * y(iprot) * dratdumdy1(irfe56_aux3) &
2527 6044 : - y(ihe4) * y(ife54) * dratdumdy1(irfe56_aux4)
2528 :
2529 :
2530 : ! fe56 jacobian elements
2531 :
2532 6044 : dfdy(ife56,ihe4) = y(ife54) * ratdum(irfe56_aux4)
2533 :
2534 :
2535 : dfdy(ife56,ife54) = &
2536 : y(ineut) * y(ineut) * ratdum(irfe56_aux2) + &
2537 6044 : y(ihe4) * ratdum(irfe56_aux4)
2538 :
2539 : dfdy(ife56,ife56) = - fe56ec_fake_factor * ratdum(irn56ec) &
2540 : - ratdum(irfe56_aux1) &
2541 6044 : - y(iprot) * y(iprot) * ratdum(irfe56_aux3)
2542 :
2543 6044 : if (plus_co56) then
2544 4 : dfdy(ife56,ico56) = ratdum(irco56ec)
2545 : else
2546 6040 : dfdy(ife56,ini56) = ratdum(irn56ec)
2547 : end if
2548 :
2549 :
2550 : dfdy(ife56,ineut) = &
2551 : -y(ife56) * dratdumdy1(irfe56_aux1) &
2552 : + 2.0d0 * y(ife54) * y(ineut) * ratdum(irfe56_aux2) &
2553 6044 : + y(ife54) * y(ineut) * y(ineut) * dratdumdy1(irfe56_aux2)
2554 :
2555 :
2556 : dfdy(ife56,iprot) = -2.0d0 * y(ife56) * y(iprot) * ratdum(irfe56_aux3) &
2557 : - y(ife56) * y(iprot) * y(iprot) * dratdumdy1(irfe56_aux3) &
2558 6044 : + y(ihe4) * y(ife54) * dratdumdy1(irfe56_aux4)
2559 :
2560 6044 : if (plus_co56) then
2561 : ! co56 jacobian elements
2562 4 : dfdy(ico56,ini56) = ratdum(irn56ec)
2563 4 : dfdy(ico56,ico56) = -ratdum(irco56ec)
2564 : end if
2565 :
2566 :
2567 : ! ni56 jacobian elements
2568 : dfdy(ini56,ihe4) = y(ife52) * ratdum(irfeag) &
2569 6044 : + y(ife52) * y(iprot) * ratdum(ir7f54)
2570 :
2571 : dfdy(ini56,ife52) = y(ihe4) * ratdum(irfeag) &
2572 6044 : + y(ihe4)* y(iprot) * ratdum(ir7f54)
2573 :
2574 6044 : dfdy(ini56,ife54) = y(iprot) * y(iprot) * ratdum(ir3f54)
2575 :
2576 : dfdy(ini56,ini56) = -ratdum(irniga) &
2577 : - ratdum(ir4f54) &
2578 : - y(iprot) * ratdum(ir8f54) &
2579 6044 : - ratdum(irn56ec)
2580 :
2581 : dfdy(ini56,iprot) = 2.0d0 * y(ife54) * y(iprot) * ratdum(ir3f54) &
2582 : + y(ife54) * y(iprot) * y(iprot) * dratdumdy1(ir3f54) &
2583 : - y(ini56) * dratdumdy1(ir4f54) &
2584 : + y(ife52) * y(ihe4)* ratdum(ir7f54) &
2585 : + y(ife52) * y(ihe4)* y(iprot) * dratdumdy1(ir7f54) &
2586 : - y(ini56) * ratdum(ir8f54) &
2587 6044 : - y(ini56) * y(iprot) * dratdumdy1(ir8f54)
2588 :
2589 :
2590 : ! photodisintegration neutrons jacobian elements
2591 6044 : dfdy(ineut,ihe4) = 2.0d0 * ratdum(iralf1)
2592 :
2593 6044 : dfdy(ineut,ife52) = -2.0d0 * y(ineut) * y(ineut) * ratdum(ir2f54)
2594 :
2595 : dfdy(ineut,ife54) = 2.0d0 * ratdum(ir1f54) &
2596 6044 : - 2.0d0 * y(ineut) * y(ineut) * ratdum(irfe56_aux2)
2597 :
2598 : dfdy(ineut,ife56) = 2.0d0 * ratdum(irfe56_aux1) &
2599 6044 : - fe56ec_n_neut * fe56ec_fake_factor * ratdum(irn56ec)
2600 :
2601 : dfdy(ineut,ineut) = &
2602 : 2.0d0 * y(ife54) * dratdumdy1(ir1f54) &
2603 : - 4.0d0 * y(ife52) * y(ineut) * ratdum(ir2f54) &
2604 : - 2.0d0 * y(ife52) * y(ineut) * y(ineut) * dratdumdy1(ir2f54) &
2605 : + 2.0d0 * y(ihe4) * dratdumdy1(iralf1) &
2606 : - 4.0d0 * y(ineut) * y(iprot)*y(iprot) * ratdum(iralf2) &
2607 : - 2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy1(iralf2) &
2608 : - ratdum(irnep) &
2609 : + 2.0d0 * y(ife56) * dratdumdy1(irfe56_aux1) &
2610 : - 4.0d0 * y(ife54) * y(ineut) * ratdum(irfe56_aux2) &
2611 6044 : - 2.0d0 * y(ife54) * y(ineut) * y(ineut) * dratdumdy1(irfe56_aux2)
2612 :
2613 : dfdy(ineut,iprot) = 2.0d0 * y(ihe4) * dratdumdy2(iralf1) &
2614 : - 4.0d0 * y(ineut)*y(ineut) * y(iprot) * ratdum(iralf2) &
2615 : - 2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy2(iralf2) &
2616 6044 : + ratdum(irpen)
2617 :
2618 : ! photodisintegration protons jacobian elements
2619 : dfdy(iprot,ihe4) = 2.0d0 * y(ife52) * ratdum(ir6f54) &
2620 : + 2.0d0 * ratdum(iralf1) &
2621 6044 : + 2.0d0 * y(ife54) * ratdum(irfe56_aux4)
2622 :
2623 6044 : dfdy(iprot,ife52) = 2.0d0 * y(ihe4) * ratdum(ir6f54)
2624 :
2625 : dfdy(iprot,ife54) = -2.0d0 * y(iprot) * y(iprot) * ratdum(ir3f54) &
2626 : - 2.0d0 * y(iprot) * y(iprot) * ratdum(ir5f54) &
2627 6044 : + 2.0d0 * y(ihe4) * ratdum(irfe56_aux4)
2628 :
2629 6044 : dfdy(iprot,ife56) = -2.0d0 * y(iprot) * y(iprot) * ratdum(irfe56_aux3)
2630 :
2631 6044 : dfdy(iprot,ini56) = 2.0d0 * ratdum(ir4f54)
2632 :
2633 : dfdy(iprot,ineut) = 2.0d0 * y(ihe4) * dratdumdy1(iralf1) &
2634 : - 4.0d0 * y(ineut) * y(iprot)*y(iprot) * ratdum(iralf2) &
2635 : - 2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy1(iralf2) &
2636 6044 : + ratdum(irnep)
2637 :
2638 : dfdy(iprot,iprot) = -4.0d0 * y(ife54) * y(iprot) * ratdum(ir3f54) &
2639 : - 2.0d0 * y(ife54) * y(iprot)*y(iprot)*dratdumdy1(ir3f54) &
2640 : + 2.0d0 * y(ini56) * dratdumdy1(ir4f54) &
2641 : - 4.0d0 * y(ife54) * y(iprot) * ratdum(ir5f54) &
2642 : - 2.0d0 * y(ife54) * y(iprot)*y(iprot)*dratdumdy1(ir5f54) &
2643 : + 2.0d0 * y(ihe4) * y(ife52) * dratdumdy1(ir6f54) &
2644 : + 2.0d0 * y(ihe4) * dratdumdy2(iralf1) &
2645 : - 4.0d0 * y(ineut)*y(ineut) * y(iprot) * ratdum(iralf2) &
2646 : - 2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy2(iralf2) &
2647 : - ratdum(irpen) &
2648 : - 4.0d0 * y(ife56) * y(iprot) * ratdum(irfe56_aux3) &
2649 : - 2.0d0 * y(ife56) * y(iprot) * y(iprot) * dratdumdy1(irfe56_aux3) &
2650 6044 : + 2.0d0 * y(ihe4) * y(ife54) * dratdumdy1(irfe56_aux4)
2651 :
2652 6044 : end subroutine approx21_dfdy
2653 :
2654 :
2655 6044 : subroutine approx21_dfdT_dfdRho( & ! epstotal includes neutrinos
2656 12088 : y, mion, dfdy, ratdum, dratdumdt, dratdumdd, &
2657 : fe56ec_fake_factor, min_T, fe56ec_n_neut, temp, den, &
2658 12088 : dfdT, dfdRho, d_epstotal_dy, plus_co56, ierr)
2659 : real(dp), intent(in), dimension(:) :: &
2660 : y, mion, ratdum, dratdumdt, dratdumdd
2661 : real(dp), intent(in) :: fe56ec_fake_factor, min_T, temp, den, dfdy(:,:)
2662 : integer, intent(in) :: fe56ec_n_neut
2663 : real(dp), intent(inout), dimension(:) :: d_epstotal_dy, dfdT, dfdRho
2664 : logical, intent(in) :: plus_co56
2665 : integer, intent(out) :: ierr
2666 :
2667 : integer :: i, j
2668 6044 : real(dp) :: enuc_conv2
2669 : logical, parameter :: deriva = .true.
2670 :
2671 : ! temperature dependence of the rate equations
2672 132972 : dfdT(1:species(plus_co56)) = 0d0
2673 : call approx21_dydt( &
2674 : y,dratdumdt,ratdum,dfdT,deriva,&
2675 6044 : fe56ec_fake_factor,min_T,fe56ec_n_neut,temp,den,plus_co56,ierr)
2676 6044 : if (ierr /= 0) return
2677 :
2678 : ! density dependence of the rate equations
2679 132972 : dfdRho(1:species(plus_co56)) = 0d0
2680 : call approx21_dydt( &
2681 : y,dratdumdd,ratdum,dfdRho,deriva,&
2682 6044 : fe56ec_fake_factor,min_T,fe56ec_n_neut,temp,den,plus_co56,ierr)
2683 6044 : if (ierr /= 0) return
2684 :
2685 : ! energy generation rate partials (total energy; do neutrinos elsewhere)
2686 132972 : enuc_conv2 = -avo*clight*clight
2687 132972 : d_epstotal_dy(1:species(plus_co56)) = 0d0
2688 132972 : do j=1,species(plus_co56)
2689 2792504 : do i=1,species(plus_co56)
2690 2792504 : d_epstotal_dy(j) = d_epstotal_dy(j) + dfdy(i,j)*mion(i)
2691 : end do
2692 132972 : d_epstotal_dy(j) = d_epstotal_dy(j) * enuc_conv2
2693 : end do
2694 :
2695 : end subroutine approx21_dfdT_dfdRho
2696 :
2697 :
2698 10 : subroutine mark_approx21(handle, ierr)
2699 : use net_def, only: Net_General_Info, get_net_ptr
2700 : use chem_def, only: chem_isos
2701 : integer, intent(in) :: handle
2702 : integer, intent(out) :: ierr
2703 : type (Net_General_Info), pointer :: g
2704 : include 'formats'
2705 10 : call get_net_ptr(handle, g, ierr)
2706 10 : if (ierr /= 0) then
2707 0 : write(*,*) 'invalid handle for do_mark_approx21_on_coprocessor'
2708 0 : return
2709 : end if
2710 : call mark_approx21_isos( &
2711 10 : g% net_iso, chem_isos% name(g% approx21_ye_iso),g% add_co56_to_approx21, ierr)
2712 10 : if (ierr /= 0) return
2713 10 : call mark_approx21_reactions(g% net_reaction,g% add_co56_to_approx21, ierr)
2714 10 : if (ierr /= 0) return
2715 10 : end subroutine mark_approx21
2716 :
2717 :
2718 10 : subroutine set_approx21(handle, ierr)
2719 10 : use net_def, only: Net_General_Info, get_net_ptr
2720 : use chem_def, only: chem_isos
2721 : integer, intent(in) :: handle
2722 : integer, intent(out) :: ierr
2723 : type (Net_General_Info), pointer :: g
2724 : include 'formats'
2725 10 : call get_net_ptr(handle, g, ierr)
2726 10 : if (ierr /= 0) then
2727 0 : write(*,*) 'invalid handle for do_mark_approx21_on_coprocessor'
2728 0 : return
2729 : end if
2730 : call set_approx21_isos( &
2731 10 : g% net_iso, chem_isos% name(g% approx21_ye_iso),g% add_co56_to_approx21,ierr)
2732 10 : if (ierr /= 0) return
2733 10 : call set_approx21_reactions(g% net_reaction,g% add_co56_to_approx21,ierr)
2734 10 : if (ierr /= 0) return
2735 10 : end subroutine set_approx21
2736 :
2737 :
2738 10 : subroutine mark_approx21_isos(itab, ye_iso_name,plus_co56, ierr)
2739 10 : use chem_lib, only: chem_get_iso_id
2740 : integer :: itab(:)
2741 : character (len=*), intent(in) :: ye_iso_name
2742 : logical, intent(in) :: plus_co56
2743 : integer, intent(out) :: ierr
2744 : integer :: i, cid
2745 10 : ierr = 0
2746 :
2747 10 : call do1('h1')
2748 10 : call do1('he3')
2749 10 : call do1('he4')
2750 10 : call do1('c12')
2751 10 : call do1('n14')
2752 10 : call do1('o16')
2753 10 : call do1('ne20')
2754 10 : call do1('mg24')
2755 10 : call do1('si28')
2756 10 : call do1('s32')
2757 10 : call do1('ar36')
2758 10 : call do1('ca40')
2759 10 : call do1('ti44')
2760 10 : call do1('cr48')
2761 10 : call do1('fe52')
2762 10 : call do1('fe54')
2763 10 : call do1('fe56')
2764 10 : if (plus_co56) call do1('co56')
2765 10 : call do1('ni56')
2766 10 : call do1('neut')
2767 10 : call do1('prot')
2768 10 : call do1(ye_iso_name)
2769 :
2770 : contains
2771 :
2772 214 : subroutine do1(str)
2773 10 : use utils_lib, only: mesa_error
2774 : character (len=*), intent(in) :: str
2775 : integer :: cid
2776 214 : cid = chem_get_iso_id(str)
2777 214 : if (cid <= 0) then
2778 0 : ierr = -1
2779 0 : write(*,*) 'mark_approx21_isos failed for ' // trim(str)
2780 0 : call mesa_error(__FILE__,__LINE__)
2781 : end if
2782 214 : itab(cid) = 1
2783 214 : end subroutine do1
2784 :
2785 : end subroutine mark_approx21_isos
2786 :
2787 :
2788 10 : subroutine set_approx21_isos(itab, ye_iso_name, plus_co56, ierr)
2789 : use chem_lib, only: chem_get_iso_id
2790 : use const_def, only: ev2erg, clight
2791 : integer :: itab(:)
2792 : character (len=*), intent(in) :: ye_iso_name
2793 : logical, intent(in) :: plus_co56
2794 : integer, intent(out) :: ierr
2795 : integer :: i, cid
2796 10 : ierr = 0
2797 :
2798 10 : ih1 = do1('h1')
2799 10 : ihe3 = do1('he3')
2800 10 : ihe4 = do1('he4')
2801 10 : ic12 = do1('c12')
2802 10 : in14 = do1('n14')
2803 10 : io16 = do1('o16')
2804 10 : ine20 = do1('ne20')
2805 10 : img24 = do1('mg24')
2806 10 : isi28 = do1('si28')
2807 10 : is32 = do1('s32')
2808 10 : iar36 = do1('ar36')
2809 10 : ica40 = do1('ca40')
2810 10 : iti44 = do1('ti44')
2811 10 : icr48 = do1('cr48')
2812 10 : ife52 = do1('fe52')
2813 10 : ife54 = do1('fe54')
2814 10 : ife56 = do1('fe56')
2815 10 : if (plus_co56) ico56 = do1('co56')
2816 10 : ini56 = do1('ni56')
2817 10 : ineut = do1('neut')
2818 10 : iprot = do1('prot')
2819 10 : icrx = do1(ye_iso_name)
2820 10 : iso_cid(icrx) = -1 ! different for different approx21 nets
2821 :
2822 : contains
2823 :
2824 214 : integer function do1(str)
2825 10 : use chem_def, only: chem_isos
2826 : use utils_lib, only: mesa_error
2827 : character (len=*), intent(in) :: str
2828 : integer :: cid
2829 214 : cid = chem_get_iso_id(str)
2830 214 : if (cid <= 0) then
2831 0 : write(*,*) 'set_approx21_isos failed for ' // trim(str)
2832 0 : call mesa_error(__FILE__,__LINE__)
2833 : end if
2834 214 : do1 = itab(cid)
2835 214 : iso_cid(do1) = cid
2836 214 : end function do1
2837 :
2838 : end subroutine set_approx21_isos
2839 :
2840 :
2841 10 : subroutine mark_approx21_reactions(rtab, plus_co56, ierr)
2842 : use rates_lib, only: rates_reaction_id
2843 : integer :: rtab(:)
2844 : logical, intent(in) :: plus_co56
2845 : integer, intent(out) :: ierr
2846 : integer :: i, ir
2847 : include 'formats'
2848 10 : ierr = 0
2849 :
2850 10 : call do1('r_he4_he4_he4_to_c12')
2851 10 : call do1('r_c12_to_he4_he4_he4')
2852 10 : call do1('r_c12_ag_o16')
2853 10 : call do1('r1212')
2854 10 : call do1('r1216')
2855 10 : call do1('r1616')
2856 10 : call do1('r_o16_ga_c12')
2857 10 : call do1('r_o16_ag_ne20')
2858 10 : call do1('r_ne20_ga_o16')
2859 10 : call do1('r_ne20_ag_mg24')
2860 10 : call do1('r_mg24_ga_ne20')
2861 :
2862 10 : call do1('r_mg24_ag_si28')
2863 10 : call do1('r_si28_ga_mg24')
2864 10 : call do1('r_mg24_ap_al27')
2865 10 : call do1('r_al27_pa_mg24')
2866 10 : call do1('r_al27_pg_si28')
2867 10 : call do1('r_si28_gp_al27')
2868 10 : call do1('r_si28_ag_s32')
2869 10 : call do1('r_s32_ga_si28')
2870 10 : call do1('r_si28_ap_p31')
2871 10 : call do1('r_p31_pa_si28')
2872 10 : call do1('r_p31_pg_s32')
2873 10 : call do1('r_s32_gp_p31')
2874 10 : call do1('r_s32_ag_ar36')
2875 10 : call do1('r_ar36_ga_s32')
2876 10 : call do1('r_s32_ap_cl35')
2877 10 : call do1('r_cl35_pa_s32')
2878 10 : call do1('r_cl35_pg_ar36')
2879 10 : call do1('r_ar36_gp_cl35')
2880 10 : call do1('r_ar36_ag_ca40')
2881 10 : call do1('r_ca40_ga_ar36')
2882 10 : call do1('r_ar36_ap_k39')
2883 10 : call do1('r_k39_pa_ar36')
2884 10 : call do1('r_k39_pg_ca40')
2885 10 : call do1('r_ca40_gp_k39')
2886 10 : call do1('r_ca40_ag_ti44')
2887 10 : call do1('r_ti44_ga_ca40')
2888 10 : call do1('r_ca40_ap_sc43')
2889 10 : call do1('r_sc43_pa_ca40')
2890 10 : call do1('r_sc43_pg_ti44')
2891 10 : call do1('r_ti44_gp_sc43')
2892 10 : call do1('r_ti44_ag_cr48')
2893 10 : call do1('r_cr48_ga_ti44')
2894 10 : call do1('r_ti44_ap_v47')
2895 10 : call do1('r_v47_pa_ti44')
2896 10 : call do1('r_v47_pg_cr48')
2897 10 : call do1('r_cr48_gp_v47')
2898 10 : call do1('r_cr48_ag_fe52')
2899 10 : call do1('r_fe52_ga_cr48')
2900 10 : call do1('r_cr48_ap_mn51')
2901 10 : call do1('r_mn51_pa_cr48')
2902 10 : call do1('r_mn51_pg_fe52')
2903 10 : call do1('r_fe52_gp_mn51')
2904 10 : call do1('r_fe52_ag_ni56')
2905 10 : call do1('r_ni56_ga_fe52')
2906 10 : call do1('r_fe52_ap_co55')
2907 10 : call do1('r_co55_pa_fe52')
2908 10 : call do1('r_co55_pg_ni56')
2909 10 : call do1('r_ni56_gp_co55')
2910 :
2911 : ! for fe54 photodisintegration
2912 10 : call do1('r_fe52_ng_fe53')
2913 10 : call do1('r_fe53_gn_fe52')
2914 10 : call do1('r_fe53_ng_fe54')
2915 10 : call do1('r_fe54_gn_fe53')
2916 10 : call do1('r_fe54_pg_co55')
2917 10 : call do1('r_co55_gp_fe54')
2918 :
2919 : ! for he4 photodisintegration
2920 10 : call do1('r_he3_ng_he4')
2921 10 : call do1('r_he4_gn_he3')
2922 10 : call do1('r_h1_ng_h2')
2923 10 : call do1('r_h2_gn_h1')
2924 10 : call do1('r_h2_pg_he3')
2925 10 : call do1('r_he3_gp_h2')
2926 :
2927 : ! for weak reactions
2928 10 : call do1('rprot_to_neut')
2929 10 : call do1('rneut_to_prot')
2930 :
2931 10 : if (plus_co56) then
2932 4 : call do1('rni56ec_to_co56')
2933 4 : call do1('rco56ec_to_fe56')
2934 : else
2935 6 : call do1('rni56ec_to_fe56')
2936 : end if
2937 :
2938 : ! ppchain
2939 10 : call do1('rpp_to_he3')
2940 10 : call do1('r_he3_he3_to_h1_h1_he4')
2941 10 : call do1('r_he3_ag_be7')
2942 10 : call do1('r_be7_wk_li7')
2943 10 : call do1('r_be7_pg_b8')
2944 :
2945 : ! cno cycles
2946 10 : call do1('r_c12_pg_n13')
2947 10 : call do1('r_n14_pg_o15')
2948 10 : call do1('r_o16_pg_f17')
2949 10 : call do1('r_n15_pg_o16')
2950 10 : call do1('r_n15_pa_c12')
2951 10 : call do1('r_n14_ag_f18')
2952 :
2953 : ! for reactions to fe56
2954 10 : call do1('r_fe54_ng_fe55')
2955 10 : call do1('r_fe55_gn_fe54')
2956 10 : call do1('r_fe55_ng_fe56')
2957 10 : call do1('r_fe56_gn_fe55')
2958 10 : call do1('r_fe54_ap_co57')
2959 10 : call do1('r_co57_pa_fe54')
2960 10 : call do1('r_fe56_pg_co57')
2961 10 : call do1('r_co57_gp_fe56')
2962 :
2963 : contains
2964 :
2965 934 : subroutine do1(str)
2966 : use utils_lib, only: mesa_error
2967 : character (len=*), intent(in) :: str
2968 : integer :: ir
2969 934 : ir = rates_reaction_id(str)
2970 934 : if (ir <= 0) then
2971 0 : ierr = -1
2972 0 : write(*,*) 'mark_approx21_reactions failed for ' // trim(str)
2973 0 : call mesa_error(__FILE__,__LINE__)
2974 : end if
2975 934 : rtab(ir) = 1
2976 934 : end subroutine do1
2977 :
2978 : end subroutine mark_approx21_reactions
2979 :
2980 :
2981 10 : subroutine set_approx21_reactions(rtab, plus_co56, ierr)
2982 : use rates_lib, only: rates_reaction_id
2983 : use utils_lib, only: mesa_error
2984 : integer :: rtab(:)
2985 : logical, intent(in) :: plus_co56
2986 : integer, intent(out) :: ierr
2987 10 : ierr = 0
2988 :
2989 10 : ir3a = do1('r_he4_he4_he4_to_c12')
2990 10 : irg3a = do1('r_c12_to_he4_he4_he4')
2991 10 : ircag = do1('r_c12_ag_o16')
2992 10 : ir1212 = do1('r1212')
2993 10 : ir1216 = do1('r1216')
2994 10 : ir1616 = do1('r1616')
2995 10 : iroga = do1('r_o16_ga_c12')
2996 10 : iroag = do1('r_o16_ag_ne20')
2997 10 : irnega = do1('r_ne20_ga_o16')
2998 10 : irneag = do1('r_ne20_ag_mg24')
2999 10 : irmgga = do1('r_mg24_ga_ne20')
3000 :
3001 10 : irmgag = do1('r_mg24_ag_si28')
3002 10 : irsiga = do1('r_si28_ga_mg24')
3003 10 : irmgap = do1('r_mg24_ap_al27')
3004 10 : iralpa = do1('r_al27_pa_mg24')
3005 10 : iralpg = do1('r_al27_pg_si28')
3006 10 : irsigp = do1('r_si28_gp_al27')
3007 10 : irsiag = do1('r_si28_ag_s32')
3008 10 : irsga = do1('r_s32_ga_si28')
3009 10 : irsiap = do1('r_si28_ap_p31')
3010 10 : irppa = do1('r_p31_pa_si28')
3011 10 : irppg = do1('r_p31_pg_s32')
3012 10 : irsgp = do1('r_s32_gp_p31')
3013 10 : irsag = do1('r_s32_ag_ar36')
3014 10 : irarga = do1('r_ar36_ga_s32')
3015 10 : irsap = do1('r_s32_ap_cl35')
3016 10 : irclpa = do1('r_cl35_pa_s32')
3017 10 : irclpg = do1('r_cl35_pg_ar36')
3018 10 : irargp = do1('r_ar36_gp_cl35')
3019 10 : irarag = do1('r_ar36_ag_ca40')
3020 10 : ircaga = do1('r_ca40_ga_ar36')
3021 10 : irarap = do1('r_ar36_ap_k39')
3022 10 : irkpa = do1('r_k39_pa_ar36')
3023 10 : irkpg = do1('r_k39_pg_ca40')
3024 10 : ircagp = do1('r_ca40_gp_k39')
3025 10 : ircaag = do1('r_ca40_ag_ti44')
3026 10 : irtiga = do1('r_ti44_ga_ca40')
3027 10 : ircaap = do1('r_ca40_ap_sc43')
3028 10 : irscpa = do1('r_sc43_pa_ca40')
3029 10 : irscpg = do1('r_sc43_pg_ti44')
3030 10 : irtigp = do1('r_ti44_gp_sc43')
3031 10 : irtiag = do1('r_ti44_ag_cr48')
3032 10 : ircrga = do1('r_cr48_ga_ti44')
3033 10 : irtiap = do1('r_ti44_ap_v47')
3034 10 : irvpa = do1('r_v47_pa_ti44')
3035 10 : irvpg = do1('r_v47_pg_cr48')
3036 10 : ircrgp = do1('r_cr48_gp_v47')
3037 10 : ircrag = do1('r_cr48_ag_fe52')
3038 10 : irfega = do1('r_fe52_ga_cr48')
3039 10 : ircrap = do1('r_cr48_ap_mn51')
3040 10 : irmnpa = do1('r_mn51_pa_cr48')
3041 10 : irmnpg = do1('r_mn51_pg_fe52')
3042 10 : irfegp = do1('r_fe52_gp_mn51')
3043 10 : irfeag = do1('r_fe52_ag_ni56')
3044 10 : irniga = do1('r_ni56_ga_fe52')
3045 10 : irfeap = do1('r_fe52_ap_co55')
3046 10 : ircopa = do1('r_co55_pa_fe52')
3047 10 : ircopg = do1('r_co55_pg_ni56')
3048 10 : irnigp = do1('r_ni56_gp_co55')
3049 :
3050 : ! for fe54 photodisintegration
3051 10 : ir52ng = do1('r_fe52_ng_fe53')
3052 10 : ir53gn = do1('r_fe53_gn_fe52')
3053 10 : ir53ng = do1('r_fe53_ng_fe54')
3054 10 : ir54gn = do1('r_fe54_gn_fe53')
3055 10 : irfepg = do1('r_fe54_pg_co55')
3056 10 : ircogp = do1('r_co55_gp_fe54')
3057 :
3058 : ! for he4 photodisintegration
3059 10 : irheng = do1('r_he3_ng_he4')
3060 10 : irhegn = do1('r_he4_gn_he3')
3061 10 : irhng = do1('r_h1_ng_h2')
3062 10 : irdgn = do1('r_h2_gn_h1')
3063 10 : irdpg = do1('r_h2_pg_he3')
3064 10 : irhegp = do1('r_he3_gp_h2')
3065 :
3066 : ! for weak reactions
3067 10 : irpen = do1('rprot_to_neut')
3068 10 : irnep = do1('rneut_to_prot')
3069 :
3070 10 : if (plus_co56) then
3071 4 : irn56ec = do1('rni56ec_to_co56')
3072 4 : irco56ec = do1('rco56ec_to_fe56')
3073 : else
3074 6 : irn56ec = do1('rni56ec_to_fe56')
3075 : end if
3076 :
3077 : ! ppchain
3078 10 : irpp = do1('rpp_to_he3')
3079 10 : ir33 = do1('r_he3_he3_to_h1_h1_he4')
3080 10 : irhe3ag = do1('r_he3_ag_be7')
3081 10 : ir_be7_wk_li7 = do1('r_be7_wk_li7')
3082 10 : ir_be7_pg_b8 = do1('r_be7_pg_b8')
3083 :
3084 : ! cno cycles
3085 10 : ircpg = do1('r_c12_pg_n13')
3086 10 : irnpg = do1('r_n14_pg_o15')
3087 10 : iropg = do1('r_o16_pg_f17')
3088 10 : irn15pg = do1('r_n15_pg_o16')
3089 10 : irn15pa = do1('r_n15_pa_c12')
3090 10 : irnag = do1('r_n14_ag_f18')
3091 :
3092 : ! for reactions to fe56
3093 10 : ir54ng = do1('r_fe54_ng_fe55')
3094 10 : ir55gn = do1('r_fe55_gn_fe54')
3095 10 : ir55ng = do1('r_fe55_ng_fe56')
3096 10 : ir56gn = do1('r_fe56_gn_fe55')
3097 10 : irfe54ap = do1('r_fe54_ap_co57')
3098 10 : irco57pa = do1('r_co57_pa_fe54')
3099 10 : irfe56pg = do1('r_fe56_pg_co57')
3100 10 : irco57gp = do1('r_co57_gp_fe56')
3101 :
3102 : ! the equilibrium links come after the mesa reactions
3103 10 : ifa = num_mesa_reactions(plus_co56)+1
3104 10 : ifg = ifa+1
3105 :
3106 10 : irr1 = ifg+1
3107 10 : irs1 = irr1+1
3108 10 : irt1 = irs1+1
3109 10 : iru1 = irt1+1
3110 10 : irv1 = iru1+1
3111 10 : irw1 = irv1+1
3112 10 : irx1 = irw1+1
3113 :
3114 10 : ir1f54 = irx1+1
3115 10 : ir2f54 = ir1f54+1
3116 10 : ir3f54 = ir2f54+1
3117 10 : ir4f54 = ir3f54+1
3118 10 : ir5f54 = ir4f54+1
3119 10 : ir6f54 = ir5f54+1
3120 10 : ir7f54 = ir6f54+1
3121 10 : ir8f54 = ir7f54+1
3122 :
3123 10 : iralf1 = ir8f54+1
3124 10 : iralf2 = iralf1+1
3125 :
3126 10 : irfe56_aux1 = iralf2+1
3127 10 : irfe56_aux2 = irfe56_aux1+1
3128 10 : irfe56_aux3 = irfe56_aux2+1
3129 10 : irfe56_aux4 = irfe56_aux3+1
3130 :
3131 10 : if( (plus_co56 .and. irfe56_aux4 /= num_reactions(plus_co56)) .or. &
3132 10 : (.not.plus_co56 .and. irfe56_aux4 /= num_reactions(plus_co56))) then
3133 0 : write(*,*) 'set_approx21_reactions found bad num_reactions'
3134 0 : write(*,*) plus_co56,irfe56_aux4,num_reactions(plus_co56)
3135 0 : call mesa_error(__FILE__,__LINE__)
3136 : end if
3137 :
3138 10 : call init_approx21(plus_co56)
3139 :
3140 : contains
3141 :
3142 934 : integer function do1(str)
3143 : use utils_lib, only: mesa_error
3144 : character (len=*), intent(in) :: str
3145 : integer :: ir
3146 934 : ir = rates_reaction_id(str)
3147 934 : if (ir <= 0) then
3148 0 : write(*,*) 'set_approx21_reactions failed for ' // trim(str)
3149 0 : call mesa_error(__FILE__,__LINE__)
3150 : end if
3151 934 : do1 = rtab(ir)
3152 934 : if (do1 <= 0) then
3153 0 : write(*,*) 'set_approx21_reactions failed to find rate for ' // trim(str)
3154 0 : call mesa_error(__FILE__,__LINE__)
3155 : end if
3156 934 : rate_id(do1) = ir
3157 934 : end function do1
3158 :
3159 : end subroutine set_approx21_reactions
3160 :
3161 :
3162 : ! call this after have set rate numbers
3163 10 : subroutine init_approx21(plus_co56)
3164 : integer :: i
3165 : logical, intent(in) :: plus_co56
3166 : include 'formats'
3167 :
3168 : ! set the names of the reaction rates (use mesa standard names)
3169 :
3170 10 : ratnam(ir3a) = 'r_he4_he4_he4_to_c12'
3171 10 : ratnam(irg3a) = 'r_c12_to_he4_he4_he4'
3172 10 : ratnam(ircag) = 'r_c12_ag_o16'
3173 10 : ratnam(ir1212) = 'r1212'
3174 10 : ratnam(ir1216) = 'r1216'
3175 10 : ratnam(ir1616) = 'r1616'
3176 10 : ratnam(iroga) = 'r_o16_ga_c12'
3177 10 : ratnam(iroag) = 'r_o16_ag_ne20'
3178 10 : ratnam(irnega) = 'r_ne20_ga_o16'
3179 10 : ratnam(irneag) = 'r_ne20_ag_mg24'
3180 10 : ratnam(irmgga) = 'r_mg24_ga_ne20'
3181 :
3182 10 : ratnam(irmgag) = 'r_mg24_ag_si28'
3183 10 : ratnam(irsiga) = 'r_si28_ga_mg24'
3184 10 : ratnam(irmgap) = 'r_mg24_ap_al27'
3185 10 : ratnam(iralpa) = 'r_al27_pa_mg24'
3186 10 : ratnam(iralpg) = 'r_al27_pg_si28'
3187 10 : ratnam(irsigp) = 'r_si28_gp_al27'
3188 10 : ratnam(irsiag) = 'r_si28_ag_s32'
3189 10 : ratnam(irsga) = 'r_s32_ga_si28'
3190 10 : ratnam(irsiap) = 'r_si28_ap_p31'
3191 10 : ratnam(irppa) = 'r_p31_pa_si28'
3192 10 : ratnam(irppg) = 'r_p31_pg_s32'
3193 10 : ratnam(irsgp) = 'r_s32_gp_p31'
3194 10 : ratnam(irsag) = 'r_s32_ag_ar36'
3195 10 : ratnam(irarga) = 'r_ar36_ga_s32'
3196 10 : ratnam(irsap) = 'r_s32_ap_cl35'
3197 10 : ratnam(irclpa) = 'r_cl35_pa_s32'
3198 10 : ratnam(irclpg) = 'r_cl35_pg_ar36'
3199 10 : ratnam(irargp) = 'r_ar36_gp_cl35'
3200 10 : ratnam(irarag) = 'r_ar36_ag_ca40'
3201 10 : ratnam(ircaga) = 'r_ca40_ga_ar36'
3202 10 : ratnam(irarap) = 'r_ar36_ap_k39'
3203 10 : ratnam(irkpa) = 'r_k39_pa_ar36'
3204 10 : ratnam(irkpg) = 'r_k39_pg_ca40'
3205 10 : ratnam(ircagp) = 'r_ca40_gp_k39'
3206 10 : ratnam(ircaag) = 'r_ca40_ag_ti44'
3207 10 : ratnam(irtiga) = 'r_ti44_ga_ca40'
3208 10 : ratnam(ircaap) = 'r_ca40_ap_sc43'
3209 10 : ratnam(irscpa) = 'r_sc43_pa_ca40'
3210 10 : ratnam(irscpg) = 'r_sc43_pg_ti44'
3211 10 : ratnam(irtigp) = 'r_ti44_gp_sc43'
3212 10 : ratnam(irtiag) = 'r_ti44_ag_cr48'
3213 10 : ratnam(ircrga) = 'r_cr48_ga_ti44'
3214 10 : ratnam(irtiap) = 'r_ti44_ap_v47'
3215 10 : ratnam(irvpa) = 'r_v47_pa_ti44'
3216 10 : ratnam(irvpg) = 'r_v47_pg_cr48'
3217 10 : ratnam(ircrgp) = 'r_cr48_gp_v47'
3218 10 : ratnam(ircrag) = 'r_cr48_ag_fe52'
3219 10 : ratnam(irfega) = 'r_fe52_ga_cr48'
3220 10 : ratnam(ircrap) = 'r_cr48_ap_mn51'
3221 10 : ratnam(irmnpa) = 'r_mn51_pa_cr48'
3222 10 : ratnam(irmnpg) = 'r_mn51_pg_fe52'
3223 10 : ratnam(irfegp) = 'r_fe52_gp_mn51'
3224 10 : ratnam(irfeag) = 'r_fe52_ag_ni56'
3225 10 : ratnam(irniga) = 'r_ni56_ga_fe52'
3226 10 : ratnam(irfeap) = 'r_fe52_ap_co55'
3227 10 : ratnam(ircopa) = 'r_co55_pa_fe52'
3228 10 : ratnam(ircopg) = 'r_co55_pg_ni56'
3229 10 : ratnam(irnigp) = 'r_ni56_gp_co55'
3230 :
3231 : ! for fe54 photodisintegration
3232 10 : ratnam(ir52ng) = 'r_fe52_ng_fe53'
3233 10 : ratnam(ir53gn) = 'r_fe53_gn_fe52'
3234 10 : ratnam(ir53ng) = 'r_fe53_ng_fe54'
3235 10 : ratnam(ir54gn) = 'r_fe54_gn_fe53'
3236 10 : ratnam(irfepg) = 'r_fe54_pg_co55'
3237 10 : ratnam(ircogp) = 'r_co55_gp_fe54'
3238 :
3239 : ! for he4 photodisintegration
3240 10 : ratnam(irheng) = 'r_he3_ng_he4'
3241 10 : ratnam(irhegn) = 'r_he4_gn_he3'
3242 10 : ratnam(irhng) = 'r_h1_ng_h2'
3243 10 : ratnam(irdgn) = 'r_h2_gn_h1'
3244 10 : ratnam(irdpg) = 'r_h2_pg_he3'
3245 10 : ratnam(irhegp) = 'r_he3_gp_h2'
3246 :
3247 : ! for weak reactions
3248 10 : ratnam(irpen) = 'rprot_to_neut'
3249 10 : ratnam(irnep) = 'rneut_to_prot'
3250 :
3251 10 : if (plus_co56) then
3252 4 : ratnam(irn56ec) = 'r_ni56_wk_co56'
3253 4 : ratnam(irco56ec) = 'r_co56_wk_fe56'
3254 : else
3255 6 : ratnam(irn56ec) = 'rni56ec_to_fe56'
3256 : end if
3257 :
3258 : ! ppchain
3259 10 : ratnam(irpp) = 'rpp_to_he3'
3260 10 : ratnam(ir33) = 'r_he3_he3_to_h1_h1_he4'
3261 10 : ratnam(irhe3ag) = 'r_he3_ag_be7'
3262 10 : ratnam(ir_be7_wk_li7) = 'r_be7_wk_li7'
3263 10 : ratnam(ir_be7_pg_b8) = 'r_be7_pg_b8'
3264 :
3265 : ! cno cycles
3266 10 : ratnam(ircpg) = 'r_c12_pg_n13'
3267 10 : ratnam(irnpg) = 'r_n14_pg_o15'
3268 10 : ratnam(iropg) = 'r_o16_pg_f17'
3269 :
3270 10 : ratnam(irn15pg) = 'r_n15_pg_o16'
3271 10 : ratnam(irn15pa) = 'r_n15_pa_c12'
3272 :
3273 10 : ratnam(irnag) = 'r_n14_ag_f18'
3274 :
3275 10 : ratnam(ir54ng) = 'r_fe54_ng_fe55'
3276 10 : ratnam(ir55gn) = 'r_fe55_gn_fe54'
3277 10 : ratnam(ir55ng) = 'r_fe55_ng_fe56'
3278 10 : ratnam(ir56gn) = 'r_fe56_gn_fe55'
3279 10 : ratnam(irfe54ap) = 'r_fe54_ap_co57'
3280 10 : ratnam(irco57pa) = 'r_co57_pa_fe54'
3281 10 : ratnam(irfe56pg) = 'r_fe56_pg_co57'
3282 10 : ratnam(irco57gp) = 'r_co57_gp_fe56'
3283 :
3284 :
3285 : ! the combo links
3286 :
3287 10 : ratnam(ifa) ='fa' ! this is fraction of n15 that goes to c12 by pa
3288 10 : ratnam(ifg) ='fg' ! this is fraction of n15 that goes to o16 by pg
3289 :
3290 10 : ratnam(irr1) = 'r1'
3291 10 : ratnam(irs1) = 's1'
3292 10 : ratnam(irt1) = 't1'
3293 10 : ratnam(iru1) = 'u1'
3294 10 : ratnam(irv1) = 'v1'
3295 10 : ratnam(irw1) = 'w1'
3296 10 : ratnam(irx1) = 'x1'
3297 :
3298 10 : ratnam(ir1f54) = 'r1f54'
3299 10 : ratnam(ir2f54) = 'r2f54'
3300 10 : ratnam(ir3f54) = 'r3f54' ! rfe54prot_to_ni56
3301 10 : ratnam(ir4f54) = 'r4f54' ! rni56gprot_to_fe54
3302 10 : ratnam(ir5f54) = 'r5f54'
3303 10 : ratnam(ir6f54) = 'r6f54'
3304 10 : ratnam(ir7f54) = 'r7f54' ! rfe52aprot_to_ni56
3305 10 : ratnam(ir8f54) = 'r8f54' ! rni56gprot_to_fe52
3306 :
3307 10 : ratnam(iralf1) = 'ralf1'
3308 10 : ratnam(iralf2) = 'ralf2'
3309 :
3310 10 : ratnam(irfe56_aux1) = 'rfe56aux1'
3311 10 : ratnam(irfe56_aux2) = 'rfe56aux2'
3312 10 : ratnam(irfe56_aux3) = 'rfe56aux3'
3313 10 : ratnam(irfe56_aux4) = 'rfe56aux4'
3314 :
3315 : return
3316 :
3317 : do i=1,num_mesa_reactions(plus_co56)
3318 : write(*,2) trim(ratnam(i)), i
3319 : end do
3320 : write(*,*) ''
3321 : do i=num_mesa_reactions(plus_co56)+1,num_reactions(plus_co56)
3322 : write(*,2) 'extra ' // trim(ratnam(i)), i
3323 : end do
3324 : call mesa_error(__FILE__,__LINE__,'init_approx21')
3325 :
3326 : end subroutine init_approx21
3327 :
3328 36344 : real(dp) function eval_fe56ec_fake_factor(fe56ec_fake_factor,min_T,temp)
3329 : real(dp), intent(in) :: fe56ec_fake_factor,min_T,temp
3330 :
3331 36344 : eval_fe56ec_fake_factor = 0.d0
3332 36344 : if(temp >= min_T)then
3333 29660 : eval_fe56ec_fake_factor = fe56ec_fake_factor
3334 : end if
3335 :
3336 0 : end function eval_fe56ec_fake_factor
3337 :
3338 :
3339 27342 : pure integer function num_reactions(plus_co56)
3340 : logical, intent(in) :: plus_co56
3341 :
3342 27342 : if(plus_co56) then
3343 : num_reactions = approx21_plus_co56_nrat
3344 : else
3345 27322 : num_reactions = approx21_nrat
3346 : end if
3347 :
3348 27332 : end function num_reactions
3349 :
3350 :
3351 9116 : pure integer function num_mesa_reactions(plus_co56)
3352 : logical, intent(in) :: plus_co56
3353 :
3354 9116 : if(plus_co56) then
3355 : num_mesa_reactions = approx21_num_mesa_reactions_co56
3356 : else
3357 9108 : num_mesa_reactions = approx21_num_mesa_reactions_21
3358 : end if
3359 :
3360 0 : end function num_mesa_reactions
3361 :
3362 21194 : pure integer function species(plus_co56)
3363 : logical, intent(in) :: plus_co56
3364 :
3365 60520 : if(plus_co56) then
3366 : species = species_co56
3367 : else
3368 60484 : species = species_21
3369 : end if
3370 :
3371 0 : end function species
3372 :
3373 :
3374 : end module net_approx21
|
