program public_hotcno include 'implno.dek' include 'const.dek' include 'timers.dek' include 'burn_common.dek' include 'network.dek' c..this program exercises the hotcno network c..declare integer i,j,nok,nbad double precision tstart,tstep,conserv,tin,din,ein,vin,zin,xin(21), 1 tout,dout,eout,xout(21) c..initialize the network call init_hotcno c..keep coming back to here 20 continue call net_input(tstart,tstep,tin,din,vin,zin,ein,xin) c..start the clock call zsecond(timzer) c..a message write(6,*) write(6,*) 'starting integration' write(6,*) c..burn it call burner(tstep,tin,din,ein,xin, 1 tout,dout,eout,xout, 1 conserv,nok,nbad) c..output a summary of the integration call net_summary(tstep,tin,din,ein,tout,dout,eout,conserv, 1 nbad,nok,xout) c..back for another input point goto 20 end subroutine burner(tstep,tin,din,ein,xin,tout,dout,eout,xout, 1 conserv,nok,nbad) include 'implno.dek' include 'const.dek' include 'network.dek' c..given time step tstep, temperature tin, density din, thermal c..energy ein, and the composition xin, this routine returns the c..burned composition xout, final temperature tout, final density dout, c..and the final thermal energy eout. c..declare the pass integer nok,nbad double precision tstep,tin,din,ein,xin(1), 1 tout,dout,eout,xout(1),conserv c..local variables integer i,k,kk double precision abar,zbar c..for the integration driver integer kount,nbig double precision beg,stptry,stpmin,tend,ys2(abignet*nzmax), 1 odescal,tol parameter (tol = 1.0d-6, 1 odescal = 1.0d-6) external hotcno,shotcno,bhotcno,dhotcno c external forder_ma28 c external forder_umf c external forder_y12m c external forder_ludcmp c external forder_leqs c external forder_lapack c external forder_gift c external forder_biconj c external rosen_ma28 c external rosen_umf c external rosen_y12m c external rosen_ludcmp c external rosen_leqs c external rosen_lapack c external rosen_gift c external rosen_biconj external stifbs_ma28 c external stifbs_umf c external stifbs_y12m c external stifbs_ludcmp c external stifbs_leqs c external stifbs_lapack c external stifbs_gift c external stifbs_biconj c..load the mass fractions do i=1,ionmax xmass(i) = xin(i) enddo c..get abar, zbar and a few other composition variables call azbar(xmass,aion,zion,ionmax, 1 ymass,abar,zbar) c..stuff the initial conditions into ys2 do i=1,ionmax ys2(i) = ymass(i) enddo c..energy, temperature, density ys2(iener) = ein ys2(itemp) = tin ys2(iden) = din c..single step (tend=tstep), hydrostatic, or expansion ending times c..the variable tstep has two meanings here. tstep in single step mode c..is the size of the time step to try. tstep in hydrostatic or expansion c..mode is the ending integration time. the integration driver really c..gets some exercise if tstep is large in single step mode. beg = 0.0d0 tend = tstep if (one_step) then stptry = tstep stpmin = tstep * 1.0d-20 else stptry = max(beg * 1.0d-10,1.0d-16) stpmin = stptry * 1.0d-12 end if c..integrate the hotcno network call netint(beg,stptry,stpmin,tend,ys2, 1 tol,neqs,nok,nbad,kount,odescal, c 4 hotcno,shotcno,bhotcno,forder_ma28) c 4 hotcno,shotcno,bhotcno,forder_umf) c 4 hotcno,shotcno,bhotcno,forder_y12m) c 4 hotcno,dhotcno,bhotcno,forder_ludcmp) c 4 hotcno,dhotcno,bhotcno,forder_leqs) c 4 hotcno,dhotcno,bhotcno,forder_lapack) c 4 hotcno,dhotcno,bhotcno,forder_gift) c 4 hotcno,shotcno,bhotcno,forder_biconj) c 4 hotcno,shotcno,bhotcno,rosen_ma28) c 4 hotcno,shotcno,bhotcno,rosen_umf) c 4 hotcno,shotcno,bhotcno,rosen_y12m) c 4 hotcno,dhotcno,bhotcno,rosen_ludcmp) c 4 hotcno,dhotcno,bhotcno,rosen_leqs) c 4 hotcno,dhotcno,bhotcno,rosen_lapack) c 4 hotcno,dhotcno,bhotcno,rosen_gift) c 4 hotcno,shotcno,bhotcno,rosen_biconj) 4 hotcno,shotcno,bhotcno,stifbs_ma28) c 4 hotcno,shotcno,bhotcno,stifbs_umf) c 4 hotcno,shotcno,bhotcno,stifbs_y12m) c 4 hotcno,dhotcno,bhotcno,stifbs_ludcmp) c 4 hotcno,dhotcno,bhotcno,stifbs_leqs) c 4 hotcno,dhotcno,bhotcno,stifbs_lapack) c 4 hotcno,dhotcno,bhotcno,stifbs_gift) c 4 hotcno,shotcno,bhotcno,stifbs_biconj) c..set the output composition do i=1,ionmax xout(i) = ys2(i) * aion(i) enddo c..output temperature, density, and thermal energy tout = ys2(itemp) dout = ys2(iden) eout = ys2(iener) c..set the mass non-conservation conserv = 0.0d0 do i=1,ionmax conserv = conserv + xout(i) enddo conserv = 1.0d0 - conserv return end c..this file contains hotcno network c.. c..routine hotcno sets up the odes c..routine rhs evaluates the right hand sides c..routine dhotcno sets up the dense hotcno jacobian c..routine bhotcno build the nonzero locations for shotcno c..routine shotcno sets up the sparse hotcno jacobian c..routine hotcnorat generates the reaction rates for routine hotcno c..routine hotcnotab generates the raw rates from table interpolation c..routine screen_hotcno applies screening corrections to the raw rates c..routine init_hotcno initializes the hotcno network subroutine hotcno(tt,y,dydt) include 'implno.dek' include 'const.dek' include 'burn_common.dek' include 'network.dek' include 'vector_eos.dek' c..this routine sets up the system of ode's for 4 cno cycles, c..the 2 hot cno cycles, triple alfa, and an rp process approximation to ni56. c.. c..cycle1: h1, he4, c12, n13, c13, n14, o15, n15 c..cycle2: o16, f17, o17 c..cycle3: f18, o18 c..cycle4: f19 c..hot : o14 ne18 c..rp : ne19, ne20, mg22, s30 ni56 c..declare the pass double precision tt,y(1),dydt(1) c..local variables integer i double precision enuc,taud,taut,snucno,sum1,sum2,z, 2 zbarxx,ytot1,abar,zbar,ye,snuda,snudz double precision conv parameter (conv = ev2erg*1.0d6*avo) c..positive definite mass fractions do i=1,ionmax y(i) = min(1.0d0,max(y(i),1.0d-30)) enddo c..generate abar and zbar for this composition zbarxx = 0.0d0 ytot1 = 0.0d0 do i=1,ionmax ytot1 = ytot1 + y(i) zbarxx = zbarxx + zion(i) * y(i) enddo abar = 1.0d0/ytot1 zbar = zbarxx * abar ye = zbar * ytot1 c..positive definite temperatures and densities y(itemp) = min(1.0d11,max(y(itemp),1.0d4)) y(iden) = min(1.0d11,max(y(iden),1.0d-10)) c..for evolution equations with the network if (pure_network .eq. 0) then c..get the new temperature and density if need be if (trho_hist) call update2(tt,y(itemp),y(iden)) if (self_heat_const_pres) then jlo_eos = 1 jhi_eos = 1 ptot_row(1) = bpres temp_row(1) = y(itemp) abar_row(1) = abar zbar_row(1) = zbar den_row(1) = y(iden) call invert_helm_pt_quiet y(iden) = den_row(1) end if c..positive definite temperatures and densities y(itemp) = min(1.0d11,max(y(itemp),1.0d4)) y(iden) = min(1.0d11,max(y(iden),1.0d-10)) c..set the common block temperature and density btemp = y(itemp) bden = y(iden) c..for pure network else if (pure_network .eq. 1) then if (trho_hist) call update2(tt,btemp,bden) end if c..get the reaction rates if (use_tables .eq. 1) then call hotcnotab else call hotcnorat end if c..do the screening here because the corrections depend on the composition call screen_hotcno(y) c..get the right hand side of the odes call rhs(y,ratdum,dydt) c..if we are doing a pure network, we are done if (pure_network .eq. 1) return c..instantaneous energy generation rate enuc = 0.0d0 do i=1,ionmax enuc = enuc + dydt(i) * bion(i) enddo enuc = enuc * conv sdot = enuc c..get the neutrino losses call sneut5(btemp,bden,abar,zbar, 1 sneut,dsneutdt,dsneutdd,snuda,snudz) c..get the cno specific neutrino losses sneutcno = snucno(y(in13),bion(ic13),bion(in13), 1 y(io14),bion(in14),bion(io14), 2 y(io15),bion(in15),bion(io15), 3 y(if17),bion(io17),bion(if17), 4 y(if18),bion(io18),bion(if18)) c..sum 'em sneut = sneut + sneutcno c..append an energy equation dydt(iener) = enuc - sneut c..the type of temperature and density ode's depend c..on the burn mode: c..hydrostatic or single step cases if (hydrostatic .or. one_step .or. trho_hist) then dydt(itemp) = 0.0d0 dydt(iden) = 0.0d0 c..adiabatic expansion or contraction else if (expansion) then taud = 446.0d0/sqrt(den0) taut = 3.0d0 * taud dydt(itemp) = -psi * y(itemp)/taut dydt(iden) = -psi * y(iden)/taud c..self heating else if (self_heat_const_den) then c..call an eos temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos c..density equation dydt(iden) = 0.0d0 c..sum1 is d(ener)/d(yi)*d(yi)/dt = d(ener)/d(abar)*d(abar)/d(yi)*d(yi)/dt sum1 = 0.0d0 do i=1,ionmax sum1 = sum1 - dydt(i) enddo sum1 = sum1 * dea_row(1)*abar*abar c..sum2 is d(ener)/d(yi)*d(yi)/dt = d(ener)/d(zbar)*d(zbar)/d(yi)*d(yi)/dt sum2 = 0.0d0 do i=1,ionmax sum2 = sum2 + (zion(i) - zbar)*dydt(i) enddo sum2 = sum2 * dez_row(1)*abar c..temperature equation that is self-consistent with an eos z = 1.0d0/cv_row(1) dydt(itemp) = z*(dydt(iener) - ded_row(1)*dydt(iden) 1 - sum1 - sum2) end if return end subroutine rhs(y,rate,dydt) include 'implno.dek' include 'const.dek' include 'burn_common.dek' include 'network.dek' c..evaluates the right hand side of the achain odes c..declare the pass double precision y(1),rate(1),dydt(1) c..local variables integer i double precision sixth parameter (sixth = 1.0d0/6.0d0) c..zero the odes do i=1,neqs dydt(i) = 0.0d0 enddo c..add in the hot cno contributions c..h1 reactions dydt(ih1) = -y(ic12)*y(ih1)*rate(irc12pg) 1 + y(in13)*rate(irn13gp) 2 - y(ic13)*y(ih1)*rate(irc13pg) 3 + y(in14)*rate(irn14gp) 4 - y(in14)*y(ih1)*rate(irn14pg) 5 + y(io15)*rate(iro15gp) 6 - y(in15)*y(ih1)*rate(irn15pa) 7 + y(ic12)*y(ihe4)*rate(irc12ap) 8 - y(in15)*y(ih1)*rate(irn15pg) 9 + y(io16)*rate(iro16gp) & - y(io16)*y(ih1)*rate(iro16pg) 1 + y(if17)*rate(irf17gp) dydt(ih1) = dydt(ih1) 1 - y(io17)*y(ih1)*rate(iro17pa) 2 + y(in14)*y(ihe4)*rate(irn14ap) 3 - y(io17)*y(ih1)*rate(iro17pg) 4 + y(if18)*rate(irf18gp) 5 - y(io18)*y(ih1)*rate(iro18pa) 6 + y(in15)*y(ihe4)*rate(irn15ap) 7 - y(io18)*y(ih1)*rate(iro18pg) 8 + y(if19)*rate(irf19gp) 9 - y(if19)*y(ih1)*rate(irf19pa) & + y(io16)*y(ihe4)*rate(iro16ap) 1 - y(in13)*y(ih1)*rate(irn13pg) 2 + y(io14)*rate(iro14gp) dydt(ih1) = dydt(ih1) 1 + y(io14) * y(ihe4) * rate(iro14ap) 2 - y(if17) * y(ih1) * rate(irf17pa) 3 - y(if17) * y(ih1) * rate(irf17pg) 4 + y(ine18) * rate(irne18gp) 5 - y(if18) * y(ih1) * rate(irf18pa) 6 + y(io15) * y(ihe4) * rate(iro15ap) 7 - 8.0d0 * y(img22) * y(ih1) * 8 rate(iralam1)*(1.0d0 - ratdum(irdelta1)) 9 - 26.0d0 * y(is30) * y(ih1) * & rate(iralam2)*(1.0d0 - ratdum(irdelta2)) 1 - 3.0d0*y(ine19)*y(ih1)*rate(irne19pg) 2 - 2.0d0*y(ine20)*y(ih1)*rate(irne20pg) c..he4 reactions dydt(ihe4) = y(in15)*y(ih1)*rate(irn15pa) 1 - y(ic12)*y(ihe4)*rate(irc12ap) 2 + y(io17)*y(ih1)*rate(iro17pa) 3 - y(in14)*y(ihe4)*rate(irn14ap) 4 + y(io18)*y(ih1)*rate(iro18pa) 5 - y(in15)*y(ihe4)*rate(irn15ap) 6 + y(if19)*y(ih1)*rate(irf19pa) 7 - y(io16)*y(ihe4)*rate(iro16ap) 8 - y(io14) * y(ihe4) * rate(iro14ap) 9 + y(if17) * y(ih1) * rate(irf17pa) & + y(if18) * y(ih1) * rate(irf18pa) dydt(ihe4) = dydt(ihe4) 1 - y(io15) * y(ihe4) * rate(iro15ap) 2 - 0.5d0 * y(ihe4) * y(ihe4) * y(ihe4) * rate(ir3a) 3 + 3.0d0 * y(ic12) * rate(irg3a) 4 - y(io16) * y(ihe4) * rate(iroag) 5 - 2.0d0 * y(img22)* y(ihe4) * 6 rate(iralam1) * ratdum(irdelta1) 7 - 6.5d0 * y(is30) * y(ihe4) * 8 rate(iralam2) * ratdum(irdelta2) 9 + y(io16) * rate(iroga) & - y(ic12) * y(ihe4) * rate(ircag) 1 - y(ine18) * y(ihe4) * rate(irne18ap) 2 - y(io15) * y(ihe4) * rate(iro15ag) c..c12 reactions dydt(ic12) = -y(ic12)*y(ih1)*rate(irc12pg) 1 + y(in13)*rate(irn13gp) 2 + y(in15)*y(ih1)*rate(irn15pa) 3 - y(ic12)*y(ihe4)*rate(irc12ap) 4 + sixth * y(ihe4) * y(ihe4) * y(ihe4) * rate(ir3a) 5 - y(ic12) * rate(irg3a) 6 + y(io16) * rate(iroga) 7 - y(ic12) * y(ihe4) * rate(ircag) c..c13 reactions dydt(ic13) = y(in13)*rate(irn13enu) 1 - y(ic13)*y(ih1)*rate(irc13pg) 2 + y(in14)*rate(irn14gp) c..n13 reactions dydt(in13) = y(ic12)*y(ih1)*rate(irc12pg) 1 - y(in13)*rate(irn13gp) 2 - y(in13)*rate(irn13enu) 3 - y(in13)*y(ih1)*rate(irn13pg) 4 + y(io14)*rate(iro14gp) c..n14 reactions dydt(in14) = y(ic13)*y(ih1)*rate(irc13pg) 1 - y(in14)*rate(irn14gp) 2 - y(in14)*y(ih1)*rate(irn14pg) 3 + y(io15)*rate(iro15gp) 4 + y(io17)*y(ih1)*rate(iro17pa) 5 - y(in14)*y(ihe4)*rate(irn14ap) 6 + y(io14)*rate(iro14enu) c..n15 reactions dydt(in15) = y(io15)*rate(iro15enu) 1 - y(in15)*y(ih1)*rate(irn15pa) 2 + y(ic12)*y(ihe4)*rate(irc12ap) 3 - y(in15)*y(ih1)*rate(irn15pg) 4 + y(io16)*rate(iro16gp) 5 + y(io18)*y(ih1)*rate(iro18pa) 6 - y(in15)*y(ihe4)*rate(irn15ap) c..o14 reactions dydt(io14) = y(in13)*y(ih1)*rate(irn13pg) 1 - y(io14)*rate(iro14gp) 2 - y(io14)*rate(iro14enu) 3 - y(io14) * y(ihe4) * rate(iro14ap) 4 + y(if17) * y(ih1) * rate(irf17pa) c..o15 reactions dydt(io15) = y(in14)*y(ih1)*rate(irn14pg) 1 - y(io15)*rate(iro15gp) 2 - y(io15)*rate(iro15enu) 3 + y(if18) * y(ih1) * rate(irf18pa) 4 - y(io15) * y(ihe4) * rate(iro15ap) 5 - y(io15) * y(ihe4) * rate(iro15ag) c..o16 reactions dydt(io16) = y(in15)*y(ih1)*rate(irn15pg) 1 - y(io16)*rate(iro16gp) 2 - y(io16)*y(ih1)*rate(iro16pg) 3 + y(if17)*rate(irf17gp) 4 + y(if19)*y(ih1)*rate(irf19pa) 5 - y(io16)*y(ihe4)*rate(iro16ap) 6 - y(io16) * y(ihe4) * rate(iroag) 7 - y(io16) * rate(iroga) 8 + y(ic12) * y(ihe4) * rate(ircag) c..o17 reactions dydt(io17) = y(if17)*rate(irf17enu) 1 - y(io17)*y(ih1)*rate(iro17pa) 2 + y(in14)*y(ihe4)*rate(irn14ap) 3 - y(io17)*y(ih1)*rate(iro17pg) 4 + y(if18)*rate(irf18gp) c..o18 reactions dydt(io18) = y(if18)*rate(irf18enu) 1 - y(io18)*y(ih1)*rate(iro18pa) 2 + y(in15)*y(ihe4)*rate(irn15ap) 3 - y(io18)*y(ih1)*rate(iro18pg) 4 + y(if19)*rate(irf19gp) c..f17 reactions dydt(if17) = y(io16)*y(ih1)*rate(iro16pg) 1 - y(if17)*rate(irf17gp) 2 - y(if17)*rate(irf17enu) 3 + y(io14) * y(ihe4) * rate(iro14ap) 4 - y(if17) * y(ih1) * rate(irf17pa) 5 - y(if17) * y(ih1) * rate(irf17pg) 6 + y(ine18) * rate(irne18gp) c..f18 reactions dydt(if18) = y(io17)*y(ih1)*rate(iro17pg) 1 - y(if18)*rate(irf18gp) 2 - y(if18)*rate(irf18enu) 3 + y(ine18) * rate(irne18enu) 4 - y(if18) * y(ih1) * rate(irf18pa) 5 + y(io15) * y(ihe4) * rate(iro15ap) c..f19 reactions dydt(if19) = y(io18)*y(ih1)*rate(iro18pg) 1 - y(if19)*rate(irf19gp) 2 - y(if19)*y(ih1)*rate(irf19pa) 3 + y(io16)*y(ihe4)*rate(iro16ap) 4 + y(ine19) * rate(irne19enu) c..ne18 reactions dydt(ine18) = y(if17) * y(ih1) * rate(irf17pg) 1 - y(ine18) * rate(irne18gp) 2 - y(ine18) * rate(irne18enu) 3 - y(ine18) * y(ihe4) * rate(irne18ap) c..ne19 reactions dydt(ine19) = y(io15) * y(ihe4) * rate(iro15ag) 1 - y(ine19) * rate(irne19enu) 2 - y(ine19)*y(ih1)*rate(irne19pg) c..ne20 reactions dydt(ine20) = y(io16) * y(ihe4) * rate(iroag) 1 - y(ine20)*y(ih1)*rate(irne20pg) c..mg22 reactions c..for rp o16(a,g)ne20(p,g)na21(p,g)mg22 c.. f17(p,g)ne18(a,g)mg22 c.. o15(a,g)ne19(p,g)na20(p,g)mg21(e-nu)na21(p,g)mg22 dydt(img22) = -y(img22)*y(ih1)* 1 rate(iralam1)*(1.0d0 - ratdum(irdelta1)) 2 - y(img22)*y(ihe4)*rate(iralam1)*ratdum(irdelta1) 3 + y(ine18) * y(ihe4) * rate(irne18ap) 4 + y(ine19)*y(ih1)*rate(irne19pg) 5 + y(ine20)*y(ih1)*rate(irne20pg) c..s30 reactions dydt(is30) = y(img22)*y(ih1)* 1 rate(iralam1)*(1.0d0 - ratdum(irdelta1)) 2 + y(img22)*y(ihe4)*rate(iralam1)*ratdum(irdelta1) 3 - y(is30)*y(ih1)* 4 rate(iralam2)*(1.0d0 - ratdum(irdelta2)) 5 - y(is30)*y(ihe4)*rate(iralam2)*ratdum(irdelta2) c..ni56 reactions dydt(ini56) = y(is30)*y(ih1)* 1 rate(iralam2)*(1.0d0 - ratdum(irdelta2)) 2 + y(is30)*y(ihe4)*rate(iralam2)*ratdum(irdelta2) return end subroutine dhotcno(tt,y,dfdy,nlog,nphys) include 'implno.dek' include 'const.dek' include 'burn_common.dek' include 'network.dek' include 'vector_eos.dek' c..this routine sets up the dense hotcno jacobian c..declare the pass integer nlog,nphys double precision tt,y(1),dfdy(nphys,nphys) c..locals integer i,j double precision zbarxx,ytot1,abar,zbar,ye,taud,taut, 1 snuda,snudz,sum1,sum2,xx,zz double precision conv,sixth parameter (conv = ev2erg*1.0d6*avo) parameter (sixth = 1.0d0/6.0d0) c..zero the jacobian do j=1,nlog do i=1,nlog dfdy(i,j) = 0.0d0 enddo enddo c..positive definite mass fractions do i=1,ionmax y(i) = min(1.0d0,max(y(i),1.0d-30)) enddo c..generate abar and zbar for this composition zbarxx = 0.0d0 ytot1 = 0.0d0 do i=1,ionmax ytot1 = ytot1 + y(i) zbarxx = zbarxx + zion(i) * y(i) enddo abar = 1.0d0/ytot1 zbar = zbarxx * abar ye = zbar * ytot1 c..positive definite temperatures and densities y(itemp) = min(1.0d11,max(y(itemp),1.0d4)) y(iden) = min(1.0d11,max(y(iden),1.0d-10)) c..for evolution equations with the network if (pure_network .eq. 0) then c..get the new temperature and density if need be if (trho_hist) call update2(tt,y(itemp),y(iden)) if (self_heat_const_pres) then jlo_eos = 1 jhi_eos = 1 ptot_row(1) = bpres temp_row(1) = y(itemp) abar_row(1) = abar zbar_row(1) = zbar den_row(1) = y(iden) call invert_helm_pt_quiet y(iden) = den_row(1) end if c..positive definite temperatures and densities y(itemp) = min(1.0d11,max(y(itemp),1.0d4)) y(iden) = min(1.0d11,max(y(iden),1.0d-10)) c..set the common block temperature and density btemp = y(itemp) bden = y(iden) c..for pure network else if (pure_network .eq. 1) then if (trho_hist) call update2(tt,btemp,bden) end if c..get the reaction rates if (use_tables .eq. 1) then call hotcnotab else call hotcnorat end if c..do the screening here because the corrections depend on the composition call screen_hotcno(y) c..h1 jacobian elements dfdy(ih1,ih1) = -y(ic12)*ratdum(irc12pg) 1 - y(ic13)*ratdum(irc13pg) 2 - y(in14)*ratdum(irn14pg) 3 - y(in15)*ratdum(irn15pa) 4 - y(in15)*ratdum(irn15pg) 5 - y(io16)*ratdum(iro16pg) 6 - y(io17)*ratdum(iro17pa) 7 - y(io17)*ratdum(iro17pg) 8 - y(io18)*ratdum(iro18pa) 9 - y(io18)*ratdum(iro18pg) & - y(if19)*ratdum(irf19pa) 1 - y(in13)*ratdum(irn13pg) 2 - y(if17) * ratdum(irf17pa) 3 - y(if17) * ratdum(irf17pg) 4 - y(if18) * ratdum(irf18pa) dfdy(ih1,ih1) = dfdy(ih1,ih1) 1 - 3.0d0*y(ine19)*ratdum(irne19pg) 2 - 2.0d0*y(ine20)*ratdum(irne20pg) 3 - 8.0d0 * y(img22) * 4 ratdum(iralam1)*(1.0d0 - ratdum(irdelta1)) 5 - 26.0d0 * y(is30) * 6 ratdum(iralam2)*(1.0d0 - ratdum(irdelta2)) dfdy(ih1,ihe4) = dfdy(ih1,ihe4) 1 + y(ic12)*ratdum(irc12ap) 1 + y(in14)*ratdum(irn14ap) 2 + y(in15)*ratdum(irn15ap) 3 + y(io16)*ratdum(iro16ap) 4 + y(io14) * ratdum(iro14ap) 5 + y(io15) * ratdum(iro15ap) 6 - 8.0d0 * y(img22) * y(ih1) * 7 dratdumdy1(iralam1)*(1.0d0 - ratdum(irdelta1)) 8 - 26.0d0 * y(is30) * y(ih1) * 9 dratdumdy1(iralam2)*(1.0d0 - ratdum(irdelta2)) dfdy(ih1,ic12) = -y(ih1)*ratdum(irc12pg) 1 + y(ihe4)*ratdum(irc12ap) dfdy(ih1,ic13) = -y(ih1)*ratdum(irc13pg) dfdy(ih1,in13) = ratdum(irn13gp) 1 - y(ih1)*ratdum(irn13pg) dfdy(ih1,in14) = ratdum(irn14gp) 1 - y(ih1)*ratdum(irn14pg) 2 + y(ihe4)*ratdum(irn14ap) dfdy(ih1,in15) = -y(ih1)*ratdum(irn15pa) 1 - y(ih1)*ratdum(irn15pg) 2 + y(ihe4)*ratdum(irn15ap) dfdy(ih1,io14) = y(ihe4) * ratdum(iro14ap) 1 + ratdum(iro14gp) dfdy(ih1,io15) = ratdum(iro15gp) 1 + y(ihe4) * ratdum(iro15ap) dfdy(ih1,io16) = ratdum(iro16gp) 1 - y(ih1)*ratdum(iro16pg) 2 + y(ihe4)*ratdum(iro16ap) dfdy(ih1,io17) = -y(ih1)*ratdum(iro17pa) 1 - y(ih1)*ratdum(iro17pg) dfdy(ih1,io18) = -y(ih1)*ratdum(iro18pa) 1 - y(ih1)*ratdum(iro18pg) dfdy(ih1,if17) = ratdum(irf17gp) 1 - y(ih1) * ratdum(irf17pa) 2 - y(ih1) * ratdum(irf17pg) dfdy(ih1,if18) = ratdum(irf18gp) 1 - y(ih1) * ratdum(irf18pa) dfdy(ih1,if19) = ratdum(irf19gp) 1 - y(ih1)*ratdum(irf19pa) dfdy(ih1,ine18) = ratdum(irne18gp) dfdy(ih1,ine19) = -3.0d0*y(ih1)*ratdum(irne19pg) dfdy(ih1,ine20) = -2.0d0*y(ih1)*ratdum(irne20pg) dfdy(ih1,img22) = -8.0d0 * y(ih1) * 1 ratdum(iralam1)*(1.0d0 - ratdum(irdelta1)) dfdy(ih1,is30) = -26.0d0 * y(ih1) * 1 ratdum(iralam2)*(1.0d0 - ratdum(irdelta2)) c..he4 jacobian elements dfdy(ihe4,ih1) = y(in15)*ratdum(irn15pa) 1 + y(io17)*ratdum(iro17pa) 2 + y(io18)*ratdum(iro18pa) 3 + y(if19)*ratdum(irf19pa) 4 + y(if17) * ratdum(irf17pa) 5 + y(if18) * ratdum(irf18pa) dfdy(ihe4,ihe4) = dfdy(ihe4,ihe4) 1 - y(ic12)*ratdum(irc12ap) 1 - y(in14)*ratdum(irn14ap) 2 - y(in15)*ratdum(irn15ap) 3 - y(io16)*ratdum(iro16ap) 4 - y(io14) * ratdum(iro14ap) 5 - y(io15) * ratdum(iro15ap) 6 - 1.5d0 * y(ihe4) * y(ihe4) * ratdum(ir3a) 7 - y(io16) * ratdum(iroag) 8 - 2.0d0*y(img22)*ratdum(iralam1)*ratdum(irdelta1) 9 - 6.5d0*y(is30)*ratdum(iralam2)*ratdum(irdelta2) & - y(ic12) * ratdum(ircag) 1 - y(ine18) * ratdum(irne18ap) 2 - y(io15) * ratdum(iro15ag) dfdy(ihe4,ihe4) = dfdy(ihe4,ihe4) 1 - 2.0d0 * y(img22)* y(ihe4) * 2 dratdumdy1(iralam1) * ratdum(irdelta1) 3 - 6.5d0 * y(is30) * y(ihe4) * 4 dratdumdy1(iralam2) * ratdum(irdelta2) dfdy(ihe4,ic12) = -y(ihe4)*ratdum(irc12ap) 1 + 3.0d0 * ratdum(irg3a) 2 - y(ihe4) * ratdum(ircag) dfdy(ihe4,in14) = -y(ihe4)*ratdum(irn14ap) dfdy(ihe4,in15) = y(ih1)*ratdum(irn15pa) dfdy(ihe4,io14) = -y(ihe4) * ratdum(iro14ap) dfdy(ihe4,io15) = -y(ihe4) * ratdum(iro15ap) & - y(ihe4) * ratdum(iro15ag) dfdy(ihe4,io16) = -y(ihe4)*ratdum(iro16ap) 1 - y(ihe4) * ratdum(iroag) 2 + ratdum(iroga) dfdy(ihe4,io17) = y(ih1)*ratdum(iro17pa) dfdy(ihe4,io18) = y(ih1)*ratdum(iro18pa) dfdy(ihe4,if17) = y(ih1) * ratdum(irf17pa) dfdy(ihe4,if18) = y(ih1) * ratdum(irf18pa) dfdy(ihe4,if19) = y(ih1)*ratdum(irf19pa) dfdy(ihe4,ine18) = -y(ihe4) * ratdum(irne18ap) dfdy(ihe4,img22) = -2.0d0*y(ihe4)* 1 ratdum(iralam1)*ratdum(irdelta1) dfdy(ihe4,is30) = -6.5d0 * y(ihe4) * 1 ratdum(iralam2) * ratdum(irdelta2) c..c12 jacobian elements dfdy(ic12,ih1) = -y(ic12)*ratdum(irc12pg) 1 + y(in15)*ratdum(irn15pa) dfdy(ic12,ihe4) = -y(ic12)*ratdum(irc12ap) 1 + 0.5d0 * y(ihe4) * y(ihe4) * ratdum(ir3a) 2 - y(ic12) * ratdum(ircag) dfdy(ic12,ic12) = -y(ih1)*ratdum(irc12pg) 1 - y(ihe4)*ratdum(irc12ap) 2 - ratdum(irg3a) 3 - y(ihe4) * ratdum(ircag) dfdy(ic12,in13) = ratdum(irn13gp) dfdy(ic12,in15) = y(ih1)*ratdum(irn15pa) dfdy(ic12,io16) = ratdum(iroga) c..c13 jacobian elements dfdy(ic13,ih1) = -y(ic13)*ratdum(irc13pg) dfdy(ic13,ic13) = -y(ih1)*ratdum(irc13pg) dfdy(ic13,in13) = ratdum(irn13enu) dfdy(ic13,in14) = ratdum(irn14gp) c..n13 jacobian elements dfdy(in13,ih1) = y(ic12)*ratdum(irc12pg) 1 - y(in13)*ratdum(irn13pg) dfdy(in13,ic12) = y(ih1)*ratdum(irc12pg) dfdy(in13,in13) = -ratdum(irn13gp) 1 - ratdum(irn13enu) 2 - y(ih1)*ratdum(irn13pg) dfdy(in13,io14) = ratdum(iro14gp) c..n14 jacobian elements dfdy(in14,ih1) = y(ic13)*ratdum(irc13pg) 1 - y(in14)*ratdum(irn14pg) 2 + y(io17)*ratdum(iro17pa) dfdy(in14,ihe4) = -y(in14)*ratdum(irn14ap) dfdy(in14,ic13) = y(ih1)*ratdum(irc13pg) dfdy(in14,in14) = -ratdum(irn14gp) 1 - y(ih1)*ratdum(irn14pg) 2 - y(ihe4)*ratdum(irn14ap) dfdy(in14,io14) = ratdum(iro14enu) dfdy(in14,io15) = ratdum(iro15gp) dfdy(in14,io17) = y(ih1)*ratdum(iro17pa) c..n15 jacobian elements dfdy(in15,ih1) = -y(in15)*ratdum(irn15pa) 1 - y(in15)*ratdum(irn15pg) 2 + y(io18)*ratdum(iro18pa) dfdy(in15,ihe4) = y(ic12)*ratdum(irc12ap) 1 - y(in15)*ratdum(irn15ap) dfdy(in15,ic12) = y(ihe4)*ratdum(irc12ap) dfdy(in15,in15) = -y(ih1)*ratdum(irn15pa) 1 - y(ih1)*ratdum(irn15pg) 2 - y(ihe4)*ratdum(irn15ap) dfdy(in15,io15) = ratdum(iro15enu) dfdy(in15,io16) = ratdum(iro16gp) dfdy(in15,io18) = y(ih1)*ratdum(iro18pa) c..o14 jacobian elements dfdy(io14,ih1) = y(in13)*ratdum(irn13pg) 1 + y(if17) * ratdum(irf17pa) dfdy(io14,ihe4) = -y(io14) * ratdum(iro14ap) dfdy(io14,in13) = y(ih1)*ratdum(irn13pg) dfdy(io14,io14) = -ratdum(iro14gp) 1 - ratdum(iro14enu) 2 - y(ihe4) * ratdum(iro14ap) dfdy(io14,if17) = y(ih1) * ratdum(irf17pa) c..o15 jacobian elements dfdy(io15,ih1) = y(in14)*ratdum(irn14pg) 1 + y(if18) * ratdum(irf18pa) dfdy(io15,ihe4) = -y(io15) * ratdum(iro15ap) & - y(io15) * ratdum(iro15ag) dfdy(io15,in14) = y(ih1)*ratdum(irn14pg) dfdy(io15,io15) = -ratdum(iro15gp) 1 - ratdum(iro15enu) 2 - y(ihe4) * ratdum(iro15ap) 3 - y(ihe4) * ratdum(iro15ag) dfdy(io15,if18) = y(ih1) * ratdum(irf18pa) c..o16 jacobian elements dfdy(io16,ih1) = y(in15)*ratdum(irn15pg) 1 - y(io16)*ratdum(iro16pg) 2 + y(if19)*ratdum(irf19pa) dfdy(io16,ihe4) = -y(io16)*ratdum(iro16ap) 1 - y(io16) * ratdum(iroag) 2 + y(ic12) * ratdum(ircag) dfdy(io16,ic12) = y(ihe4) * ratdum(ircag) dfdy(io16,in15) = y(ih1)*ratdum(irn15pg) dfdy(io16,io16) = -ratdum(iro16gp) 1 - y(ih1)*ratdum(iro16pg) 2 - y(ihe4)*ratdum(iro16ap) 3 - y(ihe4) * ratdum(iroag) 4 - ratdum(iroga) dfdy(io16,if17) = ratdum(irf17gp) dfdy(io16,if19) = y(ih1)*ratdum(irf19pa) c..o17 jacobian elements dfdy(io17,ih1) = -y(io17)*ratdum(iro17pa) 1 - y(io17)*ratdum(iro17pg) dfdy(io17,ihe4) = y(in14)*ratdum(irn14ap) dfdy(io17,in14) = y(ihe4)*ratdum(irn14ap) dfdy(io17,io17) = -y(ih1)*ratdum(iro17pa) 1 - y(ih1)*ratdum(iro17pg) dfdy(io17,if17) = ratdum(irf17enu) dfdy(io17,if18) = ratdum(irf18gp) c..o18 jacobian elements dfdy(io18,ih1) = -y(io18)*ratdum(iro18pa) 1 - y(io18)*ratdum(iro18pg) dfdy(io18,ihe4) = y(in15)*ratdum(irn15ap) dfdy(io18,in15) = y(ihe4)*ratdum(irn15ap) dfdy(io18,io18) = -y(ih1)*ratdum(iro18pa) 1 - y(ih1)*ratdum(iro18pg) dfdy(io18,if18) = ratdum(irf18enu) dfdy(io18,if19) = ratdum(irf19gp) c..f17 jacobian elements dfdy(if17,ih1) = y(io16)*ratdum(iro16pg) 1 - y(if17) * ratdum(irf17pa) 2 - y(if17) * ratdum(irf17pg) dfdy(if17,ihe4) = y(io14) * ratdum(iro14ap) dfdy(if17,io14) = y(ihe4) * ratdum(iro14ap) dfdy(if17,io16) = y(ih1)*ratdum(iro16pg) dfdy(if17,if17) = -ratdum(irf17gp) 1 - ratdum(irf17enu) 2 - y(ih1) * ratdum(irf17pa) 3 - y(ih1) * ratdum(irf17pg) dfdy(if17,ine18) = ratdum(irne18gp) c..f18 jacobian elements dfdy(if18,ih1) = y(io17)*ratdum(iro17pg) 1 - y(if18) * ratdum(irf18pa) dfdy(if18,ihe4) = y(io15) * ratdum(iro15ap) dfdy(if18,io15) = y(ihe4) * ratdum(iro15ap) dfdy(if18,io17) = y(ih1)*ratdum(iro17pg) dfdy(if18,if18) = -ratdum(irf18gp) 1 - ratdum(irf18enu) 2 - y(ih1) * ratdum(irf18pa) dfdy(if18,ine18) = ratdum(irne18enu) c..f19 jacobian elements dfdy(if19,ih1) = y(io18)*ratdum(iro18pg) 1 - y(if19)*ratdum(irf19pa) dfdy(if19,ihe4) = y(io16)*ratdum(iro16ap) dfdy(if19,io16) = y(ihe4)*ratdum(iro16ap) dfdy(if19,io18) = y(ih1)*ratdum(iro18pg) dfdy(if19,if19) = -ratdum(irf19gp) 1 - y(ih1)*ratdum(irf19pa) dfdy(if19,ine19) = ratdum(irne19enu) c..ne18 jacobian elements dfdy(ine18,ih1) = y(if17) * ratdum(irf17pg) dfdy(ine18,ihe4) = -y(ine18) * ratdum(irne18ap) dfdy(ine18,if17) = y(ih1) * ratdum(irf17pg) dfdy(ine18,ine18) = -ratdum(irne18gp) 1 -ratdum(irne18enu) 2 - y(ihe4) * ratdum(irne18ap) c..ne19 jacobian elements dfdy(ine19,ih1) = -y(ine19)*ratdum(irne19pg) dfdy(ine19,ihe4) = y(io15) * ratdum(iro15ag) dfdy(ine19,io15) = y(ihe4) * ratdum(iro15ag) dfdy(ine19,ine19) = -ratdum(irne19enu) 1 - y(ih1)*ratdum(irne19pg) c..ne20 jacobian elements dfdy(ine20,ih1) = -y(ine20)*ratdum(irne20pg) dfdy(ine20,ihe4) = y(io16) * ratdum(iroag) dfdy(ine20,io16) = y(ihe4) * ratdum(iroag) dfdy(ine20,ine20) = -y(ih1)*ratdum(irne20pg) c..mg22 jacobian elements dfdy(img22,ih1) = y(ine19)*ratdum(irne19pg) 1 + y(ine20)*ratdum(irne20pg) 2 - y(img22)*ratdum(iralam1)* 3 (1.0d0 - ratdum(irdelta1)) dfdy(img22,ihe4) = -y(img22)*ratdum(iralam1)*ratdum(irdelta1) 1 + y(ine18) * ratdum(irne18ap) 2 - y(img22)*y(ih1)* 3 dratdumdy1(iralam1)*(1.0d0-ratdum(irdelta1)) 4 - y(img22)*y(ihe4)* 5 dratdumdy1(iralam1)*ratdum(irdelta1) dfdy(img22,ine18) = y(ihe4) * ratdum(irne18ap) dfdy(img22,ine19) = y(ih1)*ratdum(irne19pg) dfdy(img22,ine20) = y(ih1)*ratdum(irne20pg) dfdy(img22,img22) = -y(ih1)*ratdum(iralam1)* 1 (1.0d0 - ratdum(irdelta1)) 2 - y(ihe4)*ratdum(iralam1)*ratdum(irdelta1) c..s30 jacobian elements dfdy(is30,ih1) = y(img22)*ratdum(iralam1)* 1 (1.0d0 - ratdum(irdelta1)) 2 - y(is30)*ratdum(iralam2)* 3 (1.0d0 - ratdum(irdelta2)) dfdy(is30,ihe4) = y(img22)*ratdum(iralam1)*ratdum(irdelta1) 1 - y(is30)*ratdum(iralam2)*ratdum(irdelta2) 2 + y(img22)*y(ih1)* 3 dratdumdy1(iralam1)*(1.0d0 - ratdum(irdelta1)) 4 + y(img22)*y(ihe4)* 5 dratdumdy1(iralam1)*ratdum(irdelta1) 6 - y(is30)*y(ih1)* 7 dratdumdy1(iralam2)*(1.0d0 - ratdum(irdelta2)) 8 - y(is30)*y(ihe4)* 9 dratdumdy1(iralam2)*ratdum(irdelta2) dfdy(is30,img22) = y(ih1)*ratdum(iralam1)* 1 (1.0d0 - ratdum(irdelta1)) 2 + y(ihe4)*ratdum(iralam1)*ratdum(irdelta1) dfdy(is30,is30) = -y(ih1)*ratdum(iralam2)* 1 (1.0d0 - ratdum(irdelta2)) 2 - y(ihe4)*ratdum(iralam2)*ratdum(irdelta2) c..ni56 jacobian elements dfdy(ini56,ih1) = y(is30)*ratdum(iralam2)* 1 (1.0d0 - ratdum(irdelta2)) dfdy(ini56,ihe4) = y(is30)*ratdum(iralam2)*ratdum(irdelta2) 1 + y(is30)*y(ih1)* 2 dratdumdy1(iralam2)*(1.0d0-ratdum(irdelta2)) 3 + y(is30)*y(ihe4)* 4 dratdumdy1(iralam2)*ratdum(irdelta2) dfdy(ini56,is30) = y(ih1)*ratdum(iralam2)* 1 (1.0d0 - ratdum(irdelta2)) 2 + y(ihe4)*ratdum(iralam2)*ratdum(irdelta2) c..if we are doing a pure network, we are done if (pure_network .eq. 1) return c..append the temperature derivatives of the rate equations call rhs(y,dratdumdt,zwork1) do i=1,ionmax dfdy(i,itemp) = zwork1(i) enddo c..append the density derivatives of the rate equations call rhs(y,dratdumdd,zwork1) do i=1,ionmax dfdy(i,iden) = zwork1(i) enddo c..append the energy generation rate jacobian elements do j=1,ionmax do i=1,ionmax dfdy(iener,j) = dfdy(iener,j) + dfdy(i,j)*bion(i) enddo dfdy(iener,j) = dfdy(iener,j) * conv dfdy(iener,itemp) = dfdy(iener,itemp) + dfdy(j,itemp)*bion(j) dfdy(iener,iden) = dfdy(iener,iden) + dfdy(j,iden)*bion(j) enddo dfdy(iener,itemp) = dfdy(iener,itemp) * conv dfdy(iener,iden) = dfdy(iener,iden) * conv dsdotdt = dfdy(iener,itemp) dsdotdd = dfdy(iener,iden) c..account for the neutrino losses call sneut5(btemp,bden,abar,zbar, 1 sneut,dsneutdt,dsneutdd,snuda,snudz) do j=1,ionmax dfdy(iener,j) = dfdy(iener,j) 1 - (-abar*abar*snuda + (zion(j) - zbar)*abar*snudz) enddo dfdy(iener,itemp) = dfdy(iener,itemp) - dsneutdt dfdy(iener,iden) = dfdy(iener,iden) - dsneutdd c..for hydrostatic or one step or trho_hist burns c..all the temperature and density jacobian elements are zero, c..so there is nothing to do. c..adiabatic expansion if (expansion) then taud = 446.0d0/sqrt(den0) taut = 3.0d0 * taud dfdy(itemp,itemp) = -psi/taut dfdy(iden,iden) = -psi/taud c..for self-heating, we need the specific heat at constant volume else if (self_heat_const_den) then c..call an eos temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos zz = 1.0d0/cv_row(1) c..d(itemp)/d(yi) do j=1,ionmax dfdy(itemp,j) = zz*dfdy(iener,j) enddo xx = dea_row(1)*abar*abar*zz do j=1,ionmax do i=1,ionmax dfdy(itemp,j) = dfdy(itemp,j) - dfdy(i,j)*xx enddo enddo xx = dez_row(1)*abar*zz do j=1,ionmax do i=1,ionmax dfdy(itemp,j) = dfdy(itemp,j) - dfdy(i,j)*(zion(i)-zbar)*xx enddo enddo c..d(itemp)/d(temp) sum1 = 0.0d0 do i=1,ionmax sum1 = sum1 - dfdy(i,itemp) enddo sum1 = sum1 * dea_row(1)*abar*abar sum2 = 0.0d0 do i=1,ionmax sum2 = sum2 + (zion(i) - zbar)*dfdy(i,itemp) enddo sum2 = sum2 * dez_row(1)*abar dfdy(itemp,itemp) = zz*(dfdy(iener,itemp) - sum1 - sum2) c..d(itemp)/d(den) sum1 = 0.0d0 do i=1,ionmax sum1 = sum1 - dfdy(i,iden) enddo sum1 = sum1 * dea_row(1)*abar*abar sum2 = 0.0d0 do i=1,ionmax sum2 = sum2 + (zion(i) - zbar)*dfdy(i,iden) enddo sum2 = sum2 * dez_row(1)*abar dfdy(itemp,iden) = zz*(dfdy(iener,iden) - sum1 - sum2) end if return end subroutine bhotcno(iloc,jloc,nzo,np) include 'implno.dek' include 'network.dek' c.. c..this routine builds the nonzero matrix locations for shotcno c..input is the integer arrys iloc and jloc, both of dimension np, that c..on output contain nzo matrix element locations. c.. c..declare the pass integer np,iloc(np),jloc(np),nzo c..local variables integer i c..communicate with shotcno integer neloc parameter (neloc=250) integer eloc(neloc),nterms common /elchotcno/ eloc,nterms c..initialize nterms = 0 nzo = 0 do i=1,neloc eloc(i) = 0 enddo call tree_init(neqs) c..h1 jacobian elements call tree(ih1,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,ic12,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,ic13,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,in13,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,in14,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,in15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,io14,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,io15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,io17,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,io18,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,if17,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,if18,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,if19,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,ine18,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,ine19,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,ine20,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,img22,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ih1,is30,eloc,neloc,nterms,nzo,iloc,jloc,np) c..he4 jacobian elements call tree(ihe4,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,ic12,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,in14,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,in15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,io14,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,io15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,io17,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,io18,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,if17,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,if18,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,if19,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,ine18,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,img22,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,is30,eloc,neloc,nterms,nzo,iloc,jloc,np) c..c12 jacobian elements call tree(ic12,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ic12,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ic12,ic12,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ic12,in13,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ic12,in15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ic12,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) c..c13 jacobian elements call tree(ic13,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ic13,ic13,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ic13,in13,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ic13,in14,eloc,neloc,nterms,nzo,iloc,jloc,np) c..n13 jacobian elements call tree(in13,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in13,ic12,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in13,in13,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in13,io14,eloc,neloc,nterms,nzo,iloc,jloc,np) c..n14 jacobian elements call tree(in14,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in14,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in14,ic13,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in14,in14,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in14,io14,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in14,io15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in14,io17,eloc,neloc,nterms,nzo,iloc,jloc,np) c..n15 jacobian elements call tree(in15,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in15,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in15,ic12,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in15,in15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in15,io15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in15,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(in15,io18,eloc,neloc,nterms,nzo,iloc,jloc,np) c..o14 jacobian elements call tree(io14,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io14,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io14,in13,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io14,io14,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io14,if17,eloc,neloc,nterms,nzo,iloc,jloc,np) c..o15 jacobian elements call tree(io15,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io15,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io15,in14,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io15,io15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io15,if18,eloc,neloc,nterms,nzo,iloc,jloc,np) c..o16 jacobian elements call tree(io16,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io16,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io16,ic12,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io16,in15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io16,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io16,if17,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io16,if19,eloc,neloc,nterms,nzo,iloc,jloc,np) c..o17 jacobian elements call tree(io17,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io17,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io17,in14,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io17,io17,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io17,if17,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io17,if18,eloc,neloc,nterms,nzo,iloc,jloc,np) c..o18 jacobian elements call tree(io18,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io18,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io18,in15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io18,io18,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io18,if18,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io18,if19,eloc,neloc,nterms,nzo,iloc,jloc,np) c..f17 jacobian elements call tree(if17,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if17,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if17,io14,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if17,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if17,if17,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if17,ine18,eloc,neloc,nterms,nzo,iloc,jloc,np) c..f18 jacobian elements call tree(if18,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if18,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if18,io15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if18,io17,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if18,if18,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if18,ine18,eloc,neloc,nterms,nzo,iloc,jloc,np) c..f19 jacobian elements call tree(if19,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if19,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if19,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if19,io18,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if19,if19,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if19,ine19,eloc,neloc,nterms,nzo,iloc,jloc,np) c..ne18 jacobian elements call tree(ine18,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine18,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine18,if17,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine18,ine18,eloc,neloc,nterms,nzo,iloc,jloc,np) c..ne19 jacobian elements call tree(ine19,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine19,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine19,io15,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine19,ine19,eloc,neloc,nterms,nzo,iloc,jloc,np) c..ne20 jacobian elements call tree(ine20,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine20,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine20,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine20,ine20,eloc,neloc,nterms,nzo,iloc,jloc,np) c..mg22 jacobian elements call tree(img22,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(img22,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(img22,ine18,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(img22,ine19,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(img22,ine20,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(img22,img22,eloc,neloc,nterms,nzo,iloc,jloc,np) c..s30 jacobian elements call tree(is30,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(is30,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(is30,img22,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(is30,is30,eloc,neloc,nterms,nzo,iloc,jloc,np) c..ni56 jacobian elements call tree(ini56,ih1,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ini56,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ini56,is30,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ini56,ini56,eloc,neloc,nterms,nzo,iloc,jloc,np) c..if we are doing a pure network, we are done if (pure_network .eq. 1) return c..temperature contributions do i=1,ionmax call tree(i,itemp,eloc,neloc,nterms,nzo,iloc,jloc,np) end do c..density contributions do i=1,ionmax call tree(i,iden,eloc,neloc,nterms,nzo,iloc,jloc,np) end do c..energy equation jacobian elements do i=1,ionmax call tree(iener,i,eloc,neloc,nterms,nzo,iloc,jloc,np) enddo call tree(iener,iener,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(iener,itemp,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(iener,iden,eloc,neloc,nterms,nzo,iloc,jloc,np) c..neutrino losses do i=1,ionmax call tree(iener,i,eloc,neloc,nterms,nzo,iloc,jloc,np) enddo call tree(iener,itemp,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(iener,iden,eloc,neloc,nterms,nzo,iloc,jloc,np) c..the jacobian elements of the temperature and density equations c..depend on the burning mode c..hydrostatic or single step if (hydrostatic .or. one_step .or. trho_hist) then call tree(itemp,itemp,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(iden,iden,eloc,neloc,nterms,nzo,iloc,jloc,np) c..adiabatic expansion else if (expansion) then call tree(itemp,itemp,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(iden,iden,eloc,neloc,nterms,nzo,iloc,jloc,np) c..self heating else if (self_heat_const_den) then do i=1,ionmax call tree(itemp,i,eloc,neloc,nterms,nzo,iloc,jloc,np) enddo call tree(itemp,itemp,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(itemp,iden,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(iden,iden,eloc,neloc,nterms,nzo,iloc,jloc,np) end if c..store the number of non-zero matrix elements in common non_zero_elements = nzo c..write a diagnostic c write(6,*) ' ' c write(6,*) nzo,' matrix elements' c write(6,*) nterms,' jacobian contributions' c write(6,*) ' ' return end subroutine shotcno(tt,y,dfdy,nzo) include 'implno.dek' include 'const.dek' include 'burn_common.dek' include 'network.dek' include 'vector_eos.dek' c..this routine sets up the sparse hotcno jacobian. c..input is tt (irrelevant here) and the abundances y(1). c..output is the jacobian dfdy(nzo). c..declare the pass integer nzo double precision tt,y(1),dfdy(1) c..locals integer i,nt,iat double precision zbarxx,ytot1,abar,zbar,ye,taud,taut, 1 snuda,snudz,a1,a2,a3,a4,zz double precision conv,sixth parameter (conv = ev2erg*1.0d6*avo) parameter (sixth = 1.0d0/6.0d0) c..communicate with the jacobian builder integer neloc parameter (neloc=250) integer eloc(neloc),nterms common /elchotcno/ eloc,nterms c..initialize nt = 0 do i=1,nzo dfdy(i) = 0.0d0 enddo do i=1,ionmax xsum(i) = 0.0d0 ysum(i) = 0.0d0 zsum(i) = 0.0d0 enddo c..positive definite mass fractions do i=1,ionmax y(i) = min(1.0d0,max(y(i),1.0d-30)) enddo c..generate abar and zbar for this composition zbarxx = 0.0d0 ytot1 = 0.0d0 do i=1,ionmax ytot1 = ytot1 + y(i) zbarxx = zbarxx + zion(i) * y(i) enddo abar = 1.0d0/ytot1 zbar = zbarxx * abar ye = zbar * ytot1 c..positive definite temperatures and densities y(itemp) = min(1.0d11,max(y(itemp),1.0d4)) y(iden) = min(1.0d11,max(y(iden),1.0d-10)) c..for evolution equations with the network if (pure_network .eq. 0) then c..get the new temperature and density if need be if (trho_hist) call update2(tt,y(itemp),y(iden)) if (self_heat_const_pres) then jlo_eos = 1 jhi_eos = 1 ptot_row(1) = bpres temp_row(1) = y(itemp) abar_row(1) = abar zbar_row(1) = zbar den_row(1) = y(iden) call invert_helm_pt_quiet y(iden) = den_row(1) end if c..positive definite temperatures and densities y(itemp) = min(1.0d11,max(y(itemp),1.0d4)) y(iden) = min(1.0d11,max(y(iden),1.0d-10)) c..set the common block temperature and density btemp = y(itemp) bden = y(iden) c..for pure network else if (pure_network .eq. 1) then if (trho_hist) call update2(tt,btemp,bden) end if c..get the reaction rates if (use_tables .eq. 1) then call hotcnotab else call hotcnorat end if c..do the screening here because the corrections depend on the composition call screen_hotcno(y) c..h1 jacobian elements c..d(h1)/d(h1) a1 = -y(ic12)*ratdum(irc12pg) 1 - y(ic13)*ratdum(irc13pg) 2 - y(in14)*ratdum(irn14pg) 3 - y(in15)*ratdum(irn15pa) 4 - y(in15)*ratdum(irn15pg) 5 - y(io16)*ratdum(iro16pg) 6 - y(io17)*ratdum(iro17pa) 7 - y(io17)*ratdum(iro17pg) 8 - y(io18)*ratdum(iro18pa) 9 - y(io18)*ratdum(iro18pg) & - y(if19)*ratdum(irf19pa) 1 - y(in13)*ratdum(irn13pg) 2 - y(if17) * ratdum(irf17pa) 3 - y(if17) * ratdum(irf17pg) 4 - y(if18) * ratdum(irf18pa) a1 = a1 1 - 3.0d0*y(ine19)*ratdum(irne19pg) 2 - 2.0d0*y(ine20)*ratdum(irne20pg) 3 - 8.0d0 * y(img22) * 4 ratdum(iralam1)*(1.0d0 - ratdum(irdelta1)) 5 - 26.0d0 * y(is30) * 6 ratdum(iralam2)*(1.0d0 - ratdum(irdelta2)) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(ih1) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(ih1) - zbar) c..d(h1)/d(he4) a1 = y(ic12)*ratdum(irc12ap) 1 + y(in14)*ratdum(irn14ap) 2 + y(in15)*ratdum(irn15ap) 3 + y(io16)*ratdum(iro16ap) 4 + y(io14) * ratdum(iro14ap) 5 + y(io15) * ratdum(iro15ap) 6 - 8.0d0 * y(img22) * y(ih1) * 7 dratdumdy1(iralam1)*(1.0d0 - ratdum(irdelta1)) 8 - 26.0d0 * y(is30) * y(ih1) * 9 dratdumdy1(iralam2)*(1.0d0 - ratdum(irdelta2)) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ih1) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ih1) - zbar) c..d(h1)/d(c12) a1 = -y(ih1)*ratdum(irc12pg) 1 + y(ihe4)*ratdum(irc12ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(ih1) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(ih1) - zbar) c..d(h1)/d(c13) a1 = -y(ih1)*ratdum(irc13pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic13) = xsum(ic13) + a1 * bion(ih1) ysum(ic13) = ysum(ic13) - a1 zsum(ic13) = zsum(ic13) + a1 * (zion(ih1) - zbar) c..d(h1)/d(n13) a1 = ratdum(irn13gp) 1 - y(ih1)*ratdum(irn13pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in13) = xsum(in13) + a1 * bion(ih1) ysum(in13) = ysum(in13) - a1 zsum(in13) = zsum(in13) + a1 * (zion(ih1) - zbar) c..d(h1)/d(n14) a1 = ratdum(irn14gp) 1 - y(ih1)*ratdum(irn14pg) 2 + y(ihe4)*ratdum(irn14ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in14) = xsum(in14) + a1 * bion(ih1) ysum(in14) = ysum(in14) - a1 zsum(in14) = zsum(in14) + a1 * (zion(ih1) - zbar) c..d(h1)/d(n15) a1 = -y(ih1)*ratdum(irn15pa) 1 - y(ih1)*ratdum(irn15pg) 2 + y(ihe4)*ratdum(irn15ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in15) = xsum(in15) + a1 * bion(ih1) ysum(in15) = ysum(in15) - a1 zsum(in15) = zsum(in15) + a1 * (zion(ih1) - zbar) c..d(h1)/d(o14) a1 = y(ihe4) * ratdum(iro14ap) 1 + ratdum(iro14gp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io14) = xsum(io14) + a1 * bion(ih1) ysum(io14) = ysum(io14) - a1 zsum(io14) = zsum(io14) + a1 * (zion(ih1) - zbar) c..d(h1)/d(o15) a1 = ratdum(iro15gp) 1 + y(ihe4) * ratdum(iro15ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io15) = xsum(io15) + a1 * bion(ih1) ysum(io15) = ysum(io15) - a1 zsum(io15) = zsum(io15) + a1 * (zion(ih1) - zbar) c..d(h1)/d(o16) a1 = ratdum(iro16gp) 1 - y(ih1)*ratdum(iro16pg) 2 + y(ihe4)*ratdum(iro16ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(ih1) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(ih1) - zbar) c..d(h1)/d(o17) a1 = -y(ih1)*ratdum(iro17pa) 1 - y(ih1)*ratdum(iro17pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io17) = xsum(io17) + a1 * bion(ih1) ysum(io17) = ysum(io17) - a1 zsum(io17) = zsum(io17) + a1 * (zion(ih1) - zbar) c..d(h1)/d(o18) a1 = -y(ih1)*ratdum(iro18pa) 1 - y(ih1)*ratdum(iro18pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io18) = xsum(io18) + a1 * bion(ih1) ysum(io18) = ysum(io18) - a1 zsum(io18) = zsum(io18) + a1 * (zion(ih1) - zbar) c..d(h1)/d(f17) a1 = ratdum(irf17gp) 1 - y(ih1) * ratdum(irf17pa) 2 - y(ih1) * ratdum(irf17pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if17) = xsum(if17) + a1 * bion(ih1) ysum(if17) = ysum(if17) - a1 zsum(if17) = zsum(if17) + a1 * (zion(ih1) - zbar) c..d(h1)/d(f18) a1 = ratdum(irf18gp) 1 - y(ih1) * ratdum(irf18pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if18) = xsum(if18) + a1 * bion(ih1) ysum(if18) = ysum(if18) - a1 zsum(if18) = zsum(if18) + a1 * (zion(ih1) - zbar) c..d(h1)/d(f19) a1 = ratdum(irf19gp) 1 - y(ih1)*ratdum(irf19pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if19) = xsum(if19) + a1 * bion(ih1) ysum(if19) = ysum(if19) - a1 zsum(if19) = zsum(if19) + a1 * (zion(ih1) - zbar) c..d(h1)/d(ne18) a1 = ratdum(irne18gp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine18) = xsum(ine18) + a1 * bion(ih1) ysum(ine18) = ysum(ine18) - a1 zsum(ine18) = zsum(ine18) + a1 * (zion(ih1) - zbar) c..d(h1)/d(ne19) a1 = -3.0d0*y(ih1)*ratdum(irne19pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine19) = xsum(ine19) + a1 * bion(ih1) ysum(ine19) = ysum(ine19) - a1 zsum(ine19) = zsum(ine19) + a1 * (zion(ih1) - zbar) c..d(h1)/d(ne20) a1 = -2.0d0*y(ih1)*ratdum(irne20pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine20) = xsum(ine20) + a1 * bion(ih1) ysum(ine20) = ysum(ine20) - a1 zsum(ine20) = zsum(ine20) + a1 * (zion(ih1) - zbar) c..d(h1)/d(mg22) a1 = -8.0d0 * y(ih1)* 1 ratdum(iralam1)*(1.0d0 - ratdum(irdelta1)) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(img22) = xsum(img22) + a1 * bion(ih1) ysum(img22) = ysum(img22) - a1 zsum(img22) = zsum(img22) + a1 * (zion(ih1) - zbar) c..d(h1)/d(s30) a1 = -26.0d0 * y(ih1) * 1 ratdum(iralam2)*(1.0d0 - ratdum(irdelta2)) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(is30) = xsum(is30) + a1 * bion(ih1) ysum(is30) = ysum(is30) - a1 zsum(is30) = zsum(is30) + a1 * (zion(ih1) - zbar) c..he4 jacobian elements c..d(he4)/d(h1) a1 = y(in15)*ratdum(irn15pa) 1 + y(io17)*ratdum(iro17pa) 2 + y(io18)*ratdum(iro18pa) 3 + y(if19)*ratdum(irf19pa) 4 + y(if17) * ratdum(irf17pa) 5 + y(if18) * ratdum(irf18pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(ihe4) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(he4) a1 = -y(ic12)*ratdum(irc12ap) 1 - y(in14)*ratdum(irn14ap) 2 - y(in15)*ratdum(irn15ap) 3 - y(io16)*ratdum(iro16ap) 4 - y(io14) * ratdum(iro14ap) 5 - y(io15) * ratdum(iro15ap) 6 - 1.5d0 * y(ihe4) * y(ihe4) * ratdum(ir3a) 7 - y(io16) * ratdum(iroag) 8 - 2.0d0 * y(img22)*ratdum(iralam1)*ratdum(irdelta1) 9 - 6.5d0 * y(is30) * ratdum(iralam2) * ratdum(irdelta2) & - y(ic12) * ratdum(ircag) 1 - y(ine18) * ratdum(irne18ap) 2 - y(io15) * ratdum(iro15ag) a1 = a1 1 - 2.0d0 * y(img22)* y(ihe4) * 2 dratdumdy1(iralam1) * ratdum(irdelta1) 3 - 6.5d0 * y(is30) * y(ihe4) * 4 dratdumdy1(iralam2) * ratdum(irdelta2) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ihe4) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(c12) a1 = -y(ihe4)*ratdum(irc12ap) 1 + 3.0d0 * ratdum(irg3a) 2 - y(ihe4) * ratdum(ircag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(ihe4) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(n14) a1 = -y(ihe4)*ratdum(irn14ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in14) = xsum(in14) + a1 * bion(ihe4) ysum(in14) = ysum(in14) - a1 zsum(in14) = zsum(in14) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(n15) a1 = y(ih1)*ratdum(irn15pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in15) = xsum(in15) + a1 * bion(ihe4) ysum(in15) = ysum(in15) - a1 zsum(in15) = zsum(in15) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(o14) a1 = -y(ihe4) * ratdum(iro14ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io14) = xsum(io14) + a1 * bion(ihe4) ysum(io14) = ysum(io14) - a1 zsum(io14) = zsum(io14) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(o15) a1 = -y(ihe4) * ratdum(iro15ap) & - y(ihe4) * ratdum(iro15ag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io15) = xsum(io15) + a1 * bion(ihe4) ysum(io15) = ysum(io15) - a1 zsum(io15) = zsum(io15) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(o16) a1 = -y(ihe4)*ratdum(iro16ap) 1 - y(ihe4) * ratdum(iroag) 2 + ratdum(iroga) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(ihe4) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(o17) a1 = y(ih1)*ratdum(iro17pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io17) = xsum(io17) + a1 * bion(ihe4) ysum(io17) = ysum(io17) - a1 zsum(io17) = zsum(io17) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(o18) a1 = y(ih1)*ratdum(iro18pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io18) = xsum(io18) + a1 * bion(ihe4) ysum(io18) = ysum(io18) - a1 zsum(io18) = zsum(io18) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(f17) a1 = y(ih1) * ratdum(irf17pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if17) = xsum(if17) + a1 * bion(ihe4) ysum(if17) = ysum(if17) - a1 zsum(if17) = zsum(if17) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(f18) a1 = y(ih1) * ratdum(irf18pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if18) = xsum(if18) + a1 * bion(ihe4) ysum(if18) = ysum(if18) - a1 zsum(if18) = zsum(if18) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(f19) a1 = y(ih1)*ratdum(irf19pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if19) = xsum(if19) + a1 * bion(ihe4) ysum(if19) = ysum(if19) - a1 zsum(if19) = zsum(if19) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(ne18) a1 = -y(ihe4) * ratdum(irne18ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine18) = xsum(ine18) + a1 * bion(ihe4) ysum(ine18) = ysum(ine18) - a1 zsum(ine18) = zsum(ine18) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(mg22) a1 = -2.0d0*y(ihe4)*ratdum(iralam1)*ratdum(irdelta1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(img22) = xsum(img22) + a1 * bion(ihe4) ysum(img22) = ysum(img22) - a1 zsum(img22) = zsum(img22) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(s30) a1 = -6.5d0 * y(ihe4) * ratdum(iralam2) * ratdum(irdelta2) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(is30) = xsum(is30) + a1 * bion(ihe4) ysum(is30) = ysum(is30) - a1 zsum(is30) = zsum(is30) + a1 * (zion(ihe4) - zbar) c..c12 jacobian elements c..d(c12)/d(h1) a1 = -y(ic12)*ratdum(irc12pg) 1 + y(in15)*ratdum(irn15pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(ic12) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(ic12) - zbar) c..d(c12)/d(he4) a1 = -y(ic12)*ratdum(irc12ap) 1 + 0.5d0 * y(ihe4) * y(ihe4) * ratdum(ir3a) 2 - y(ic12) * ratdum(ircag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ic12) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ic12) - zbar) c..d(c12)/d(c12) a1 = -y(ih1)*ratdum(irc12pg) 1 - y(ihe4)*ratdum(irc12ap) 2 - ratdum(irg3a) 3 - y(ihe4) * ratdum(ircag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(ic12) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(ic12) - zbar) c..d(c12)/d(n13) a1 = ratdum(irn13gp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in13) = xsum(in13) + a1 * bion(ic12) ysum(in13) = ysum(in13) - a1 zsum(in13) = zsum(in13) + a1 * (zion(ic12) - zbar) c..d(c12)/d(n15) a1 = y(ih1)*ratdum(irn15pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in15) = xsum(in15) + a1 * bion(ic12) ysum(in15) = ysum(in15) - a1 zsum(in15) = zsum(in15) + a1 * (zion(ic12) - zbar) c..d(c12)/d(o16) a1 = ratdum(iroga) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(ic12) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(ic12) - zbar) c..c13 jacobian elements c..d(c13)/d(h1) a1 = -y(ic13)*ratdum(irc13pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(ic13) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(ic13) - zbar) c..d(c13)/d(c13) a1 = -y(ih1)*ratdum(irc13pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic13) = xsum(ic13) + a1 * bion(ic13) ysum(ic13) = ysum(ic13) - a1 zsum(ic13) = zsum(ic13) + a1 * (zion(ic13) - zbar) c..d(c13)/d(n13) a1 = ratdum(irn13enu) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in13) = xsum(in13) + a1 * bion(ic13) ysum(in13) = ysum(in13) - a1 zsum(in13) = zsum(in13) + a1 * (zion(ic13) - zbar) c..d(c13)/d(n14) a1 = ratdum(irn14gp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in14) = xsum(in14) + a1 * bion(ic13) ysum(in14) = ysum(in14) - a1 zsum(in14) = zsum(in14) + a1 * (zion(ic13) - zbar) c..n13 jacobian elements c..d(n13)/d(h1) a1 = y(ic12)*ratdum(irc12pg) 1 - y(in13)*ratdum(irn13pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(in13) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(in13) - zbar) c..d(n13)/d(c12) a1 = y(ih1)*ratdum(irc12pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(in13) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(in13) - zbar) c..d(n13)/d(n13) a1 = -ratdum(irn13gp) 1 - ratdum(irn13enu) 2 - y(ih1)*ratdum(irn13pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in13) = xsum(in13) + a1 * bion(in13) ysum(in13) = ysum(in13) - a1 zsum(in13) = zsum(in13) + a1 * (zion(in13) - zbar) c..d(n13)/d(o14) a1 = ratdum(iro14gp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io14) = xsum(io14) + a1 * bion(in13) ysum(io14) = ysum(io14) - a1 zsum(io14) = zsum(io14) + a1 * (zion(in13) - zbar) c..n14 jacobian elements c..d(n14)/d(h1) a1 = y(ic13)*ratdum(irc13pg) 1 - y(in14)*ratdum(irn14pg) 2 + y(io17)*ratdum(iro17pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(in14) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(in14) - zbar) c..d(n14)/d(he4) a1 = -y(in14)*ratdum(irn14ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(in14) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(in14) - zbar) c..d(n14)/d(c13) a1 = y(ih1)*ratdum(irc13pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic13) = xsum(ic13) + a1 * bion(in14) ysum(ic13) = ysum(ic13) - a1 zsum(ic13) = zsum(ic13) + a1 * (zion(in14) - zbar) c..d(n14)/d(n14) a1 = -ratdum(irn14gp) 1 - y(ih1)*ratdum(irn14pg) 2 - y(ihe4)*ratdum(irn14ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in14) = xsum(in14) + a1 * bion(in14) ysum(in14) = ysum(in14) - a1 zsum(in14) = zsum(in14) + a1 * (zion(in14) - zbar) c..d(n14)/d(o14) a1 = ratdum(iro14enu) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io14) = xsum(io14) + a1 * bion(in14) ysum(io14) = ysum(io14) - a1 zsum(io14) = zsum(io14) + a1 * (zion(in14) - zbar) c..d(n14)/d(o15) a1 = ratdum(iro15gp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io15) = xsum(io15) + a1 * bion(in14) ysum(io15) = ysum(io15) - a1 zsum(io15) = zsum(io15) + a1 * (zion(in14) - zbar) c..d(n14)/d(o17) a1 = y(ih1)*ratdum(iro17pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io17) = xsum(io17) + a1 * bion(in14) ysum(io17) = ysum(io17) - a1 zsum(io17) = zsum(io17) + a1 * (zion(in14) - zbar) c..n15 jacobian elements c..d(n15)/d(h1) a1 = -y(in15)*ratdum(irn15pa) 1 - y(in15)*ratdum(irn15pg) 2 + y(io18)*ratdum(iro18pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(in15) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(in15) - zbar) c..d(n15)/d(he4) a1 = y(ic12)*ratdum(irc12ap) 1 - y(in15)*ratdum(irn15ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(in15) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(in15) - zbar) c..d(n15)/d(c12) a1 = y(ihe4)*ratdum(irc12ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(in15) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(in15) - zbar) c..d(n15)/d(n15) a1 = -y(ih1)*ratdum(irn15pa) 1 - y(ih1)*ratdum(irn15pg) 2 - y(ihe4)*ratdum(irn15ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in15) = xsum(in15) + a1 * bion(in15) ysum(in15) = ysum(in15) - a1 zsum(in15) = zsum(in15) + a1 * (zion(in15) - zbar) c..d(n15)/d(o15) a1 = ratdum(iro15enu) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io15) = xsum(io15) + a1 * bion(in15) ysum(io15) = ysum(io15) - a1 zsum(io15) = zsum(io15) + a1 * (zion(in15) - zbar) c..d(n15)/d(o16) a1 = ratdum(iro16gp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(in15) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(in15) - zbar) c..d(n15)/d(o18) a1 = y(ih1)*ratdum(iro18pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io18) = xsum(io18) + a1 * bion(in15) ysum(io18) = ysum(io18) - a1 zsum(io18) = zsum(io18) + a1 * (zion(in15) - zbar) c..o14 jacobian elements c..d(o14)/d(h1) a1 = y(in13)*ratdum(irn13pg) 1 + y(if17) * ratdum(irf17pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(io14) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(io14) - zbar) c..d(o14)/d(he4) a1 = -y(io14) * ratdum(iro14ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(io14) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(io14) - zbar) c..d(o14)/d(n13) a1 = y(ih1)*ratdum(irn13pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in13) = xsum(in13) + a1 * bion(io14) ysum(in13) = ysum(in13) - a1 zsum(in13) = zsum(in13) + a1 * (zion(io14) - zbar) c..d(o14)/d(o14) a1 = -ratdum(iro14gp) 1 - ratdum(iro14enu) 2 - y(ihe4) * ratdum(iro14ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io14) = xsum(io14) + a1 * bion(io14) ysum(io14) = ysum(io14) - a1 zsum(io14) = zsum(io14) + a1 * (zion(io14) - zbar) c..d(o14)/d(f17) a1 = y(ih1) * ratdum(irf17pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if17) = xsum(if17) + a1 * bion(io14) ysum(if17) = ysum(if17) - a1 zsum(if17) = zsum(if17) + a1 * (zion(io14) - zbar) c..o15 jacobian elements c..d(o15)/d(h1) a1 = y(in14)*ratdum(irn14pg) 1 + y(if18) * ratdum(irf18pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(io15) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(io15) - zbar) c..d(o15)/d(he4) a1 = -y(io15) * ratdum(iro15ap) 1 - y(io15) * ratdum(iro15ag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(io15) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(io15) - zbar) c..d(o15)/d(n14) a1 = y(ih1)*ratdum(irn14pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in14) = xsum(in14) + a1 * bion(io15) ysum(in14) = ysum(in14) - a1 zsum(in14) = zsum(in14) + a1 * (zion(io15) - zbar) c..d(o15)/d(o15) a1 = -ratdum(iro15gp) 1 - ratdum(iro15enu) 2 - y(ihe4) * ratdum(iro15ap) 3 - y(ihe4) * ratdum(iro15ag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io15) = xsum(io15) + a1 * bion(io15) ysum(io15) = ysum(io15) - a1 zsum(io15) = zsum(io15) + a1 * (zion(io15) - zbar) c..d(o15)/d(f18) a1 = y(ih1) * ratdum(irf18pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if18) = xsum(if18) + a1 * bion(io15) ysum(if18) = ysum(if18) - a1 zsum(if18) = zsum(if18) + a1 * (zion(io15) - zbar) c..o16 jacobian elements c..d(o16)/d(h1) a1 = y(in15)*ratdum(irn15pg) 1 - y(io16)*ratdum(iro16pg) 2 + y(if19)*ratdum(irf19pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(io16) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(io16) - zbar) c..d(o16)/d(he4) a1 = -y(io16)*ratdum(iro16ap) 1 - y(io16) * ratdum(iroag) 2 + y(ic12) * ratdum(ircag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(io16) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(io16) - zbar) c..d(o16)/d(c12) a1 = y(ihe4) * ratdum(ircag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(io16) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(io16) - zbar) c..d(o16)/d(n15) a1 = y(ih1)*ratdum(irn15pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in15) = xsum(in15) + a1 * bion(io16) ysum(in15) = ysum(in15) - a1 zsum(in15) = zsum(in15) + a1 * (zion(io16) - zbar) c..d(o16)/d(o16) a1 = -ratdum(iro16gp) 1 - y(ih1)*ratdum(iro16pg) 2 - y(ihe4)*ratdum(iro16ap) 3 - y(ihe4) * ratdum(iroag) 4 - ratdum(iroga) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(io16) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(io16) - zbar) c..d(o16)/d(f17) a1 = ratdum(irf17gp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if17) = xsum(if17) + a1 * bion(io16) ysum(if17) = ysum(if17) - a1 zsum(if17) = zsum(if17) + a1 * (zion(io16) - zbar) c..d(o16)/d(f19) a1 = y(ih1)*ratdum(irf19pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if19) = xsum(if19) + a1 * bion(io16) ysum(if19) = ysum(if19) - a1 zsum(if19) = zsum(if19) + a1 * (zion(io16) - zbar) c..o17 jacobian elements c..d(o17)/d(h1) a1 = -y(io17)*ratdum(iro17pa) 1 - y(io17)*ratdum(iro17pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(io17) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(io17) - zbar) c..d(o17)/d(he4) a1 = y(in14)*ratdum(irn14ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(io17) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(io17) - zbar) c..d(o17)/d(n14) a1 = y(ihe4)*ratdum(irn14ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in14) = xsum(in14) + a1 * bion(io17) ysum(in14) = ysum(in14) - a1 zsum(in14) = zsum(in14) + a1 * (zion(io17) - zbar) c..d(o17)/d(o17) a1 = -y(ih1)*ratdum(iro17pa) 1 - y(ih1)*ratdum(iro17pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io17) = xsum(io17) + a1 * bion(io17) ysum(io17) = ysum(io17) - a1 zsum(io17) = zsum(io17) + a1 * (zion(io17) - zbar) c..d(o17)/d(f17) a1 = ratdum(irf17enu) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if17) = xsum(if17) + a1 * bion(io17) ysum(if17) = ysum(if17) - a1 zsum(if17) = zsum(if17) + a1 * (zion(io17) - zbar) c..d(o17)/d(f18) a1 = ratdum(irf18gp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if18) = xsum(if18) + a1 * bion(io17) ysum(if18) = ysum(if18) - a1 zsum(if18) = zsum(if18) + a1 * (zion(io17) - zbar) c..o18 jacobian elements c..d(o18)/d(h1) a1 = -y(io18)*ratdum(iro18pa) 1 - y(io18)*ratdum(iro18pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(io18) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(io18) - zbar) c..d(o18)/d(he4) a1 = y(in15)*ratdum(irn15ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(io18) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(io18) - zbar) c..d(o18)/d(n15) a1 = y(ihe4)*ratdum(irn15ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(in15) = xsum(in15) + a1 * bion(io18) ysum(in15) = ysum(in15) - a1 zsum(in15) = zsum(in15) + a1 * (zion(io18) - zbar) c..d(o18)/d(o18) a1 = -y(ih1)*ratdum(iro18pa) 1 - y(ih1)*ratdum(iro18pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io18) = xsum(io18) + a1 * bion(io18) ysum(io18) = ysum(io18) - a1 zsum(io18) = zsum(io18) + a1 * (zion(io18) - zbar) c..d(o18)/d(f18) a1 = ratdum(irf18enu) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if18) = xsum(if18) + a1 * bion(io18) ysum(if18) = ysum(if18) - a1 zsum(if18) = zsum(if18) + a1 * (zion(io18) - zbar) c..d(o18)/d(f19) a1 = ratdum(irf19gp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if19) = xsum(if19) + a1 * bion(io18) ysum(if19) = ysum(if19) - a1 zsum(if19) = zsum(if19) + a1 * (zion(io18) - zbar) c..f17 jacobian elements c..d(f17)/d(h1) a1 = y(io16)*ratdum(iro16pg) 1 - y(if17) * ratdum(irf17pa) 2 - y(if17) * ratdum(irf17pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(if17) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(if17) - zbar) c..d(f17)/d(he4) a1 = y(io14) * ratdum(iro14ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(if17) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(if17) - zbar) c..d(f17)/d(o14) a1 = y(ihe4) * ratdum(iro14ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io14) = xsum(io14) + a1 * bion(if17) ysum(io14) = ysum(io14) - a1 zsum(io14) = zsum(io14) + a1 * (zion(if17) - zbar) c..d(f17)/d(o16) a1 = y(ih1)*ratdum(iro16pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(if17) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(if17) - zbar) c..d(f17)/d(f17) a1 = -ratdum(irf17gp) 1 - ratdum(irf17enu) 2 - y(ih1) * ratdum(irf17pa) 3 - y(ih1) * ratdum(irf17pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if17) = xsum(if17) + a1 * bion(if17) ysum(if17) = ysum(if17) - a1 zsum(if17) = zsum(if17) + a1 * (zion(if17) - zbar) c..d(f17)/d(ne18) a1 = ratdum(irne18gp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine18) = xsum(ine18) + a1 * bion(if17) ysum(ine18) = ysum(ine18) - a1 zsum(ine18) = zsum(ine18) + a1 * (zion(if17) - zbar) c..f18 jacobian elements c..d(f18)/d(h1) a1 = y(io17)*ratdum(iro17pg) 1 - y(if18) * ratdum(irf18pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(if18) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(if18) - zbar) c..d(f18)/d(he4) a1 = y(io15) * ratdum(iro15ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(if18) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(if18) - zbar) c..d(f18)/d(o15) a1 = y(ihe4) * ratdum(iro15ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io15) = xsum(io15) + a1 * bion(if18) ysum(io15) = ysum(io15) - a1 zsum(io15) = zsum(io15) + a1 * (zion(if18) - zbar) c..d(f18)/d(o17) a1 = y(ih1)*ratdum(iro17pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io17) = xsum(io17) + a1 * bion(if18) ysum(io17) = ysum(io17) - a1 zsum(io17) = zsum(io17) + a1 * (zion(if18) - zbar) c..d(f18)/d(f18) a1 = -ratdum(irf18gp) 1 - ratdum(irf18enu) 2 - y(ih1) * ratdum(irf18pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if18) = xsum(if18) + a1 * bion(if18) ysum(if18) = ysum(if18) - a1 zsum(if18) = zsum(if18) + a1 * (zion(if18) - zbar) c..d(f18)/d(ne18) a1 = ratdum(irne18enu) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine22) = xsum(ine22) + a1 * bion(if18) ysum(ine22) = ysum(ine22) - a1 zsum(ine22) = zsum(ine22) + a1 * (zion(if18) - zbar) c..f19 jacobian elements c..d(f19)/d(h1) a1 = y(io18)*ratdum(iro18pg) 1 - y(if19)*ratdum(irf19pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(if19) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(if19) - zbar) c..d(f19)/d(he4) a1 = y(io16)*ratdum(iro16ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(if19) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(if19) - zbar) c..d(f19)/d(o16) a1 = y(ihe4)*ratdum(iro16ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(if19) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(if19) - zbar) c..d(f19)/d(o18) a1 = y(ih1)*ratdum(iro18pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io18) = xsum(io18) + a1 * bion(if19) ysum(io18) = ysum(io18) - a1 zsum(io18) = zsum(io18) + a1 * (zion(if19) - zbar) c..d(f19)/d(f19) a1 = -ratdum(irf19gp) 1 - y(ih1)*ratdum(irf19pa) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if19) = xsum(if19) + a1 * bion(if19) ysum(if19) = ysum(if19) - a1 zsum(if19) = zsum(if19) + a1 * (zion(if19) - zbar) c..d(f19)/d(ne19) a1 = ratdum(irne19enu) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine19) = xsum(ine19) + a1 * bion(if19) ysum(ine19) = ysum(ine19) - a1 zsum(ine19) = zsum(ine19) + a1 * (zion(if19) - zbar) c..ne18 jacobian elements c..d(ne18)/d(h1) a1 = y(if17) * ratdum(irf17pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(ine18) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(ine18) - zbar) c..d(ne18)/d(he4) a1 = -y(ine18) * ratdum(irne18ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ine18) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ine18) - zbar) c..d(ne18)/d(f17) a1 = y(ih1) * ratdum(irf17pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(if17) = xsum(if17) + a1 * bion(ine18) ysum(if17) = ysum(if17) - a1 zsum(if17) = zsum(if17) + a1 * (zion(ine18) - zbar) c..d(ne18)/d(ne18) a1 = -ratdum(irne18gp) 1 - ratdum(irne18enu) 2 - y(ihe4) * ratdum(irne18ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine18) = xsum(ine18) + a1 * bion(ine18) ysum(ine18) = ysum(ine18) - a1 zsum(ine18) = zsum(ine18) + a1 * (zion(ine18) - zbar) c..ne19 jacobian elements c..d(ne19)/d(h1) a1 = -y(ine19)*ratdum(irne19pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(ine19) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(ine19) - zbar) c..d(ne19)/d(he4) a1 = y(io15) * ratdum(iro15ag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ine19) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ine19) - zbar) c..d(ne19)/d(o15) a1 = y(ihe4) * ratdum(iro15ag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io15) = xsum(io15) + a1 * bion(ine19) ysum(io15) = ysum(io15) - a1 zsum(io15) = zsum(io15) + a1 * (zion(ine19) - zbar) c..d(ne19)/d(ne19) a1 = -ratdum(irne19enu) 1 - y(ih1)*ratdum(irne19pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine19) = xsum(ine19) + a1 * bion(ine19) ysum(ine19) = ysum(ine19) - a1 zsum(ine19) = zsum(ine19) + a1 * (zion(ine19) - zbar) c..ne20 jacobian elements c..d(ne20)/d(h1) a1 = -y(ine20)*ratdum(irne20pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(ine20) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(ine20) - zbar) c..d(ne20)/d(he4) a1 = y(io16) * ratdum(iroag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ine20) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ine20) - zbar) c..d(ne20)/d(o16) a1 = y(ihe4) * ratdum(iroag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(ine20) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(ine20) - zbar) c..d(ne20)/d(ne20) a1 = -y(ih1)*ratdum(irne20pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine20) = xsum(ine20) + a1 * bion(ine20) ysum(ine20) = ysum(ine20) - a1 zsum(ine20) = zsum(ine20) + a1 * (zion(ine20) - zbar) c..mg22 jacobian elements c..d(mg22)/d(h1) a1 = y(ine19)*ratdum(irne19pg) 1 + y(ine20)*ratdum(irne20pg) 2 - y(img22)*ratdum(iralam1)*(1.0d0 - ratdum(irdelta1)) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(img22) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(img22) - zbar) c..d(mg22)/d(he4) a1 = -y(img22)*ratdum(iralam1)*ratdum(irdelta1) 1 + y(ine18) * ratdum(irne18ap) 2 - y(img22)*y(ih1)* 3 dratdumdy1(iralam1)*(1.0d0-ratdum(irdelta1)) 4 - y(img22)*y(ihe4)* 5 dratdumdy1(iralam1)*ratdum(irdelta1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(img22) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(img22) - zbar) c..d(mg22)/d(ne18) a1 = y(ihe4) * ratdum(irne18ap) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine18) = xsum(ine18) + a1 * bion(img22) ysum(ine18) = ysum(ine18) - a1 zsum(ine18) = zsum(ine18) + a1 * (zion(img22) - zbar) c..d(mg22)/d(ne19) a1 = y(ih1)*ratdum(irne19pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine19) = xsum(ine19) + a1 * bion(img22) ysum(ine19) = ysum(ine19) - a1 zsum(ine19) = zsum(ine19) + a1 * (zion(img22) - zbar) c..d(mg22)/d(ne20) a1 = y(ih1)*ratdum(irne20pg) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine20) = xsum(ine20) + a1 * bion(img22) ysum(ine20) = ysum(ine20) - a1 zsum(ine20) = zsum(ine20) + a1 * (zion(img22) - zbar) c..d(mg22)/d(mg22) a1 = -y(ih1)*ratdum(iralam1)*(1.0d0 - ratdum(irdelta1)) 1 - y(ihe4)*ratdum(iralam1)*ratdum(irdelta1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(img22) = xsum(img22) + a1 * bion(img22) ysum(img22) = ysum(img22) - a1 zsum(img22) = zsum(img22) + a1 * (zion(img22) - zbar) c..s30 jacobian elements c..d(s30)/d(h1) a1 = y(img22)*ratdum(iralam1)*(1.0d0 - ratdum(irdelta1)) 1 - y(is30)*ratdum(iralam2)*(1.0d0 - ratdum(irdelta2)) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(is30) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(is30) - zbar) c..d(s30)/d(he4) a1 = y(img22)*ratdum(iralam1)*ratdum(irdelta1) 1 - y(is30)*ratdum(iralam2)*ratdum(irdelta2) 2 + y(img22)*y(ih1)* 3 dratdumdy1(iralam1)*(1.0d0 - ratdum(irdelta1)) 4 + y(img22)*y(ihe4)* 5 dratdumdy1(iralam1)*ratdum(irdelta1) 6 - y(is30)*y(ih1)* 7 dratdumdy1(iralam2)*(1.0d0 - ratdum(irdelta2)) 8 - y(is30)*y(ihe4)* 9 dratdumdy1(iralam2)*ratdum(irdelta2) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(is30) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(is30) - zbar) c..d(s30)/d(mg22) a1 = y(ih1)*ratdum(iralam1)*(1.0d0 - ratdum(irdelta1)) 1 + y(ihe4)*ratdum(iralam1)*ratdum(irdelta1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(img22) = xsum(img22) + a1 * bion(is30) ysum(img22) = ysum(img22) - a1 zsum(img22) = zsum(img22) + a1 * (zion(is30) - zbar) c..d(s30)/d(s30) a1 = -y(ih1)*ratdum(iralam2)*(1.0d0 - ratdum(irdelta2)) 1 - y(ihe4)*ratdum(iralam2)*ratdum(irdelta2) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(is30) = xsum(is30) + a1 * bion(is30) ysum(is30) = ysum(is30) - a1 zsum(is30) = zsum(is30) + a1 * (zion(is30) - zbar) c..ni56 jacobian elements c..d(ni56)/d(h1) a1 = y(is30)*ratdum(iralam2)*(1.0d0 - ratdum(irdelta2)) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ih1) = xsum(ih1) + a1 * bion(ini56) ysum(ih1) = ysum(ih1) - a1 zsum(ih1) = zsum(ih1) + a1 * (zion(ini56) - zbar) c..d(ni56)/d(he4) a1 = y(is30)*ratdum(iralam2)*ratdum(irdelta2) 1 + y(is30)*y(ih1)* 2 dratdumdy1(iralam2)*(1.0d0-ratdum(irdelta2)) 3 + y(is30)*y(ihe4)* 4 dratdumdy1(iralam2)*ratdum(irdelta2) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ini56) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ini56) - zbar) c..d(ni56)/d(s30) a1 = y(ih1)*ratdum(iralam2)*(1.0d0 - ratdum(irdelta2)) 1 + y(ihe4)*ratdum(iralam2)*ratdum(irdelta2) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(is30) = xsum(is30) + a1 * bion(ini56) ysum(is30) = ysum(is30) - a1 zsum(is30) = zsum(is30) + a1 * (zion(ini56) - zbar) c..d(ni56)/d(ni56) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ini56) = xsum(ini56) + a1 * bion(ini56) ysum(ini56) = ysum(ini56) - a1 zsum(ini56) = zsum(ini56) + a1 * (zion(ini56) - zbar) c..if we are doing a pure network, we are done if (pure_network .eq. 1) goto 678 c..append the temperature derivatives of the rate equations c..d(yi)/dtemp call rhs(y,dratdumdt,zwork1) do i=1,ionmax nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + zwork1(i) enddo c..append the density derivatives of the rate equations c..d(yi)/d(den) call rhs(y,dratdumdd,zwork2) do i=1,ionmax nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + zwork2(i) enddo c..energy jacobian elements c..d(iener)/d(ixx) do i=1,ionmax a1 = xsum(i) * conv nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 enddo c..d(iener)/d(iener) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(iener)/d(temp) a1 = 0.0d0 do i=1,ionmax a1 = a1 + zwork1(i)*bion(i) enddo a1 = a1 * conv nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 dsdotdt = dfdy(iat) c..d(iener)/d(den) a1 = 0.0d0 do i=1,ionmax a1 = a1 + zwork2(i)*bion(i) enddo a1 = a1 * conv nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 dsdotdd = dfdy(iat) c..account for the neutrino losses call sneut5(btemp,bden,abar,zbar, 1 sneut,dsneutdt,dsneutdd,snuda,snudz) c..d(ener)/d(yi) do i=1,ionmax a1 = -(-abar*abar*snuda + (zion(i) - zbar)*abar*snudz) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 enddo c..d(iener)/d(temp) a1 = -dsneutdt nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(iener)/d(den) a1 = -dsneutdd nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..the jacobian elements of the temperature and density equations c..depend on the burning mode c..hydrostatic if (hydrostatic .or. one_step .or. trho_hist) then c..d(itemp)/d(itemp) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(iden)/d(iden) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..adiabatic expansion else if (expansion) then taud = 446.0d0/sqrt(den0) taut = 3.0d0 * taud c..d(itemp)/d(itemp) a1 = -psi/taut nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(iden)/d(iden) a1 = -psi/taud nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..self heating else if (self_heat_const_den) then c..call an eos temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos zz = 1.0d0/cv_row(1) c..temperature jacobian elements c..d(itemp)/d(yi) do i=1,ionmax a1 = zz*(xsum(i) * conv 1 - ysum(i) * dea_row(1)*abar*abar 2 - zsum(i) * dez_row(1)*abar) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 enddo c..d(itemp)/d(itemp) a1 = 0.0d0 do i=1,ionmax a1 = a1 + zwork1(i)*bion(i) enddo a1 = a1*conv a2 = 0.0d0 do i=1,ionmax a2 = a2 - zwork1(i) enddo a2 = a2 * dea_row(1)*abar*abar a3 = 0.0d0 do i=1,ionmax a3 = a3 + zwork1(i)*(zion(i) - zbar) enddo a3 = a3 * dez_row(1) * abar a4 = (a1 - dsneutdt - a2 - a3) * zz nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a4 c..d(itemp)/d(iden) a1 = 0.0d0 do i=1,ionmax a1 = a1 + zwork2(i)*bion(i) enddo a1 = a1*conv a2 = 0.0d0 do i=1,ionmax a2 = a2 - zwork2(i) enddo a2 = a2 * dea_row(1)*abar*abar a3 = 0.0d0 do i=1,ionmax a3 = a3 + zwork2(i)*(zion(i) - zbar) enddo a3 = a3 * dez_row(1) * abar a4 = (a1 - dsneutdd - a2 - a3) * zz nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a4 c..d(iden)/d(iden) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..end of burning mode ifs end if c..bullet check the counting 678 if (nt .ne. nterms) then write(6,*) 'nt =',nt,' nterms =',nterms write(6,*) 'error in routine shotcno: nt .ne. nterms' stop 'error in routine shotcno' end if return end subroutine hotcnorat include 'implno.dek' include 'burn_common.dek' include 'network.dek' c.. c..this routine generates nuclear reaction rates for the hotcno network. c.. c..declare integer i double precision rrate,drratedt,drratedd c..zero the rates do i=1,nrat ratraw(i) = 0.0d0 enddo do i=1,nrat dratrawdt(i) = 0.0d0 enddo do i=1,nrat dratrawdd(i) = 0.0d0 enddo if (btemp .lt. 1.0e6) return c..get the temperature factors call tfactors(btemp) c..12c(p,g)13n call rate_c12pg(btemp,bden, 1 ratraw(irc12pg),dratrawdt(irc12pg),dratrawdd(irc12pg), 2 ratraw(irn13gp),dratrawdt(irn13gp),dratrawdd(irn13gp)) c..13n(e+nu)13c call rate_n13em(btemp,bden, 1 ratraw(irn13enu),dratrawdt(irn13enu),dratrawdd(irn13enu), 2 rrate,drratedt,drratedd) c..13c(p,g)14n call rate_c13pg(btemp,bden, 1 ratraw(irc13pg),dratrawdt(irc13pg),dratrawdd(irc13pg), 2 ratraw(irn14gp),dratrawdt(irc13pg),dratrawdd(irn14gp)) c..14n(p,g)15o call rate_n14pg(btemp,bden, 1 ratraw(irn14pg),dratrawdt(irn14pg),dratrawdd(irn14pg), 2 ratraw(iro15gp),dratrawdt(irn14pg),dratrawdd(iro15gp)) c..15o(e+nu)15n call rate_o15em(btemp,bden, 1 ratraw(iro15enu),dratrawdt(iro15enu),dratrawdd(iro15enu), 2 rrate,drratedt,drratedd) c..15n(p,a)12c fcz3 call rate_n15pa(btemp,bden, 1 ratraw(irn15pa),dratrawdt(irn15pa),dratrawdd(irn15pa), 2 ratraw(irc12ap),dratrawdt(irc12ap),dratrawdd(irc12ap)) c..15n(p,g)16o call rate_n15pg(btemp,bden, 1 ratraw(irn15pg),dratrawdt(irn15pg),dratrawdd(irn15pg), 2 ratraw(iro16gp),dratrawdt(iro16gp),dratrawdd(iro16gp)) c..16o(p,g)17f call rate_o16pg(btemp,bden, 1 ratraw(iro16pg),dratrawdt(iro16pg),dratrawdd(iro16pg), 2 ratraw(irf17gp),dratrawdt(irf17gp),dratrawdd(irf17gp)) c..17f(e+nu)17o call rate_f17em(btemp,bden, 1 ratraw(irf17enu),dratrawdt(irf17enu),dratrawdd(irf17enu), 2 rrate,drratedt,drratedd) c..17o(p,a)14n call rate_o17pa(btemp,bden, 1 ratraw(iro17pa),dratrawdt(iro17pa),dratrawdd(iro17pa), 2 ratraw(irn14ap),dratrawdt(irn14ap),dratrawdd(irn14ap)) c..17o(p,g)18f call rate_o17pg(btemp,bden, 1 ratraw(iro17pg),dratrawdt(iro17pg),dratrawdd(iro17pg), 2 ratraw(irf18gp),dratrawdt(irf18gp),dratrawdd(irf18gp)) c..18f(e+nu)18o call rate_f18em(btemp,bden, 1 ratraw(irf18enu),dratrawdt(irf18enu),dratrawdd(irf18enu), 2 rrate,drratedt,drratedd) c..18o(p,a)15n call rate_o18pa(btemp,bden, 1 ratraw(iro18pa),dratrawdt(iro18pa),dratrawdd(iro18pa), 2 ratraw(irn15ap),dratrawdt(irn15ap),dratrawdd(irn15ap)) c..18o(p,g)19f call rate_o18pg(btemp,bden, 1 ratraw(iro18pg),dratrawdt(iro18pg),dratrawdd(iro18pg), 2 ratraw(irf19gp),dratrawdt(irf19gp),dratrawdd(irf19gp)) c..19f(p,a)16o call rate_f19pa(btemp,bden, 1 ratraw(irf19pa),dratrawdt(irf19pa),dratrawdd(irf19pa), 2 ratraw(iro16ap),dratrawdt(iro16ap),dratrawdd(iro16ap)) c..add these for the hot cno cycles c..13n(p,g)14o call rate_n13pg(btemp,bden, 1 ratraw(irn13pg),dratrawdt(irn13pg),dratrawdd(irn13pg), 2 ratraw(iro14gp),dratrawdt(iro14gp),dratrawdd(iro14gp)) c..14o(e+nu)14n call rate_o14em(btemp,bden, 1 ratraw(iro14enu),dratrawdt(iro14enu),dratrawdd(iro14enu), 2 rrate,drratedt,drratedd) c..14o(a,p)17f cf88 q = 1.191 call rate_o14ap(btemp,bden, 1 ratraw(iro14ap),dratrawdt(iro14ap),dratrawdd(iro14ap), 2 ratraw(irf17pa),dratrawdt(irf17pa),dratrawdd(irf17pa)) c..17f(p,g)18ne call rate_f17pg(btemp,bden, 1 ratraw(irf17pg),dratrawdt(irf17pg),dratrawdd(irf17pg), 2 ratraw(irne18gp),dratrawdt(irne18gp),dratrawdd(irne18gp)) c..18ne(e+nu)18f call rate_ne18em(btemp,bden, 1 ratraw(irne18enu),dratrawdt(irne18enu),dratrawdd(irne18enu), 2 rrate,drratedt,drratedd) c..18f(p,a)15o call rate_f18pa(btemp,bden, 1 ratraw(irf18pa),dratrawdt(irf18pa),dratrawdd(irf18pa), 2 ratraw(iro15ap),dratrawdt(iro15ap),dratrawdd(iro15ap)) c..triple alpha to c12 call rate_tripalf(btemp,bden, 1 ratraw(ir3a),dratrawdt(ir3a),dratrawdd(ir3a), 2 ratraw(irg3a),dratrawdt(irg3a),dratrawdd(irg3a)) c..c12(a,g)o16 call rate_c12ag(btemp,bden, 1 ratraw(ircag),dratrawdt(ircag),dratrawdd(ircag), 2 ratraw(iroga),dratrawdt(iroga),dratrawdd(iroga)) c..16o(a,g)20ne call rate_o16ag(btemp,bden, 1 ratraw(iroag),dratrawdt(iroag),dratrawdd(iroag), 2 rrate,drratedt,drratedd) c..18ne(a,p)21na call rate_ne18ap(btemp,bden, 1 ratraw(irne18ap),dratrawdt(irne18ap),dratrawdd(irne18ap), 2 rrate,drratedt,drratedd) c..19ne(p,g)20na call rate_ne19pg(btemp,bden, 1 ratraw(irne19pg),dratrawdt(irne19pg),dratrawdd(irne19pg), 2 rrate,drratedt,drratedd) c..15o(a,g)19ne call rate_o15ag(btemp,bden, 1 ratraw(iro15ag),dratrawdt(iro15ag),dratrawdd(iro15ag), 2 rrate,drratedt,drratedd) c..44ti(a,p)47v call rate_ti44ap(btemp,bden, 1 ratraw(irtiap),dratrawdt(irtiap),dratrawdd(irtiap), 2 rrate,drratedt,drratedd) c g47v=1.0+exp(-1.206*t9m1+1.059+9.997d-2*t9) c ratraw(irtiap) = ratraw(irtiap)*g47v c..26si(a,p)29p call rate_si26ap(btemp,bden, 1 ratraw(irsi26ap),dratrawdt(irsi26ap),dratrawdd(irsi26ap), 2 rrate,drratedt,drratedd) c..19ne(e+nu)19f call rate_ne19em(btemp,bden, 1 ratraw(irne19enu),dratrawdt(irne19enu),dratrawdd(irne19enu), 2 rrate,drratedt,drratedd) c..20ne(p,g)21na call rate_ne20pg(btemp,bden, 1 ratraw(irne20pg),dratrawdt(irne20pg),dratrawdd(irne20pg), 2 rrate,drratedt,drratedd) return end subroutine hotcnotab include 'implno.dek' include 'burn_common.dek' include 'network.dek' c..uses tables instead of analytical expressions to evaluate the c..raw reaction rates. a cubic polynomial is hardwired for speed. integer i,j,imax,iat,mp,ifirst parameter (mp = 4) double precision tlo,thi,tstp,bden_sav,btemp_sav, 1 x,x1,x2,x3,x4,a,b,c,d,e,f,g,h,p,q, 2 alfa,beta,gama,delt data ifirst/0/ c..make the table if (ifirst .eq. 0) then ifirst = 1 c..set the log temperature loop limits thi = 10.0d0 tlo = 6.0d0 imax = int(thi-tlo)*120 + 1 if (imax .gt. nrattab) stop 'imax too small in aprox13tab' tstp = (thi - tlo)/float(imax-1) c..save the input btemp_sav = btemp bden_sav = bden c..form the table bden = 1.0d0 do i=1,imax btemp = tlo + float(i-1)*tstp btemp = 10.0d0**(btemp) call hotcnorat ttab(i) = btemp do j=1,nrat rattab(j,i) = ratraw(j) drattabdt(j,i) = dratrawdt(j) drattabdd(j,i) = dratrawdd(j) enddo enddo c..restore the input bden = bden_sav btemp = btemp_sav end if c..normal execution starts here c..set the density dependence vector dtab(irc12pg) = bden dtab(irn13gp) = 1.0d0 dtab(irn13enu) = 1.0d0 dtab(irc13pg) = bden dtab(irn14gp) = 1.0d0 dtab(irn14pg) = bden dtab(iro15gp) = 1.0d0 dtab(iro15enu) = 1.0d0 dtab(irn15pa) = bden dtab(irc12ap) = bden dtab(irn15pg) = bden dtab(iro16gp) = 1.0d0 dtab(iro16pg) = bden dtab(irf17gp) = 1.0d0 dtab(irf17enu) = 1.0d0 dtab(iro17pa) = bden dtab(irn14ap) = bden dtab(iro17pg) = bden dtab(irf18gp) = 1.0d0 dtab(irf18enu) = 1.0d0 dtab(iro18pa) = bden dtab(irn15ap) = bden dtab(iro18pg) = bden dtab(irf19gp) = 1.0d0 dtab(irf19pa) = bden dtab(iro16ap) = bden dtab(irn13pg) = bden dtab(iro14gp) = 1.0d0 dtab(iro14enu) = 1.0d0 dtab(iro14ap) = bden dtab(irf17pa) = bden dtab(irf17pg) = bden dtab(irne18gp) = 1.0d0 dtab(irne18enu)= 1.0d0 dtab(irf18pa) = bden dtab(iro15ap) = bden dtab(ir3a) = bden*bden dtab(irg3a) = 1.0d0 dtab(ircag) = bden dtab(iroga) = 1.0d0 dtab(iroag) = bden dtab(irne18ap) = bden dtab(irne19pg) = bden dtab(iro15ag) = bden dtab(irtiap) = bden dtab(irsi26ap) = bden dtab(irne19enu)= 1.0d0 dtab(irne20pg) = bden c..hash locate the temperature iat = int((log10(btemp) - tlo)/tstp) + 1 iat = max(1,min(iat - mp/2 + 1,imax - mp + 1)) c..setup the lagrange interpolation coefficients for a cubic x = btemp x1 = ttab(iat) x2 = ttab(iat+1) x3 = ttab(iat+2) x4 = ttab(iat+3) a = x - x1 b = x - x2 c = x - x3 d = x - x4 e = x1 - x2 f = x1 - x3 g = x1 - x4 h = x2 - x3 p = x2 - x4 q = x3 - x4 alfa = b*c*d/(e*f*g) beta = -a*c*d/(e*h*p) gama = a*b*d/(f*h*q) delt = -a*b*c/(g*p*q) c..crank off the raw reaction rates do j=1,nrat ratraw(j) = (alfa*rattab(j,iat) 1 + beta*rattab(j,iat+1) 2 + gama*rattab(j,iat+2) 3 + delt*rattab(j,iat+3) 4 ) * dtab(j) dratrawdt(j) = (alfa*drattabdt(j,iat) 1 + beta*drattabdt(j,iat+1) 2 + gama*drattabdt(j,iat+2) 3 + delt*drattabdt(j,iat+3) 4 ) * dtab(j) dratrawdd(j) = alfa*drattabdd(j,iat) 1 + beta*drattabdd(j,iat+1) 2 + gama*drattabdd(j,iat+2) 3 + delt*drattabdd(j,iat+3) enddo c..hand finish the three body reactions dratrawdd(ir3a) = bden * dratrawdd(ir3a) return end subroutine screen_hotcno(y) include 'implno.dek' include 'burn_common.dek' include 'network.dek' c..this routine computes the screening factors c..and applies them to the raw reaction rates, c..producing the final reaction rates used by the c..right hand sides and jacobian matrix elements c..this routine assumes screen_on = 1 or = 0 has been set at a higer level, c..presumably in the top level driver c..declare integer i,jscr double precision y(*),sc1a,sc1adt,sc1add,sc2a,sc2adt,sc2add, 1 sc3a,sc3adt,sc3add, 2 abar,zbar,z2bar,ytot1,zbarxx,z2barxx, 3 wp22mg,wp30s,alamt integer init data init/1/ c..initialize do i=1,nrat ratdum(i) = ratraw(i) dratdumdt(i) = dratrawdt(i) dratdumdd(i) = dratrawdd(i) scfac(i) = 1.0d0 dscfacdt(i) = 0.0d0 dscfacdt(i) = 0.0d0 end do c..if screening is on if (screen_on .ne. 0) then c..with the passed composition, compute abar,zbar and other variables zbarxx = 0.0d0 z2barxx = 0.0d0 ytot1 = 0.0d0 do i=1,ionmax ytot1 = ytot1 + y(i) zbarxx = zbarxx + zion(i) * y(i) z2barxx = z2barxx + zion(i) * zion(i) * y(i) enddo abar = 1.0d0/ytot1 zbar = zbarxx * abar z2bar = z2barxx * abar c..c12 + p jscr = 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ic12),aion(ic12),zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irc12pg) = ratraw(irc12pg) * sc1a dratdumdt(irc12pg) =dratrawdt(irc12pg)*sc1a+ratraw(irc12pg)*sc1adt dratdumdd(irc12pg) =dratrawdd(irc12pg)*sc1a+ratraw(irc12pg)*sc1add scfac(irc12pg) = sc1a dscfacdt(irc12pg) = sc1adt dscfacdd(irc12pg) = sc1add c..c13 + p jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ic13),aion(ic13),zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irc13pg) = ratraw(irc13pg) * sc1a dratdumdt(irc13pg) =dratrawdt(irc13pg)*sc1a+ratraw(irc13pg)*sc1adt dratdumdd(irc13pg) =dratrawdd(irc13pg)*sc1a+ratraw(irc13pg)*sc1add scfac(irc13pg) = sc1a dscfacdt(irc13pg) = sc1adt dscfacdd(irc13pg) = sc1add c..n14 + p jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(in14),aion(in14),zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irn14pg) = ratraw(irn14pg) * sc1a dratdumdt(irn14pg) =dratrawdt(irn14pg)*sc1a+ratraw(irn14pg)*sc1adt dratdumdd(irn14pg) =dratrawdd(irn14pg)*sc1a+ratraw(irn14pg)*sc1add scfac(irn14pg) = sc1a dscfacdt(irn14pg) = sc1adt dscfacdd(irn14pg) = sc1add c..n15 + p jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(in15),aion(in15),zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irn15pg) = ratraw(irn15pg) * sc1a dratdumdt(irn15pg) =dratrawdt(irn15pg)*sc1a+ratraw(irn15pg)*sc1adt dratdumdd(irn15pg) =dratrawdd(irn15pg)*sc1a+ratraw(irn15pg)*sc1add scfac(irn15pg) = sc1a dscfacdt(irn15pg) = sc1adt dscfacdd(irn15pg) = sc1add ratdum(irn15pa) = ratraw(irn15pa) * sc1a dratdumdt(irn15pa) =dratrawdt(irn15pa)*sc1a+ratraw(irn15pa)*sc1adt dratdumdd(irn15pa) =dratrawdd(irn15pa)*sc1a+ratraw(irn15pa)*sc1adt scfac(irn15pa) = sc1a dscfacdt(irn15pa) = sc1adt dscfacdd(irn15pa) = sc1add c..c12 + a jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ic12),aion(ic12),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1add,sc1add) ratdum(irc12ap) = ratraw(irc12ap) * sc1a dratdumdt(irc12ap) =dratrawdt(irc12ap)*sc1a+ratraw(irc12ap)*sc1adt dratdumdd(irc12ap) =dratrawdd(irc12ap)*sc1a+ratraw(irc12ap)*sc1add scfac(irc12ap) = sc1a dscfacdt(irc12ap) = sc1adt dscfacdd(irc12ap) = sc1add ratdum(ircag) = ratraw(ircag) * sc1a dratdumdt(ircag) = dratrawdt(ircag)*sc1a + ratraw(ircag)*sc1adt dratdumdd(ircag) = dratrawdd(ircag)*sc1a + ratraw(ircag)*sc1add scfac(ircag) = sc1a dscfacdt(ircag) = sc1adt dscfacdt(ircag) = sc1add c..o16 + p jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(io16),aion(io16),zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(iro16pg) = ratraw(iro16pg) * sc1a dratdumdt(iro16pg) =dratrawdt(iro16pg)*sc1a+ratraw(iro16pg)*sc1adt dratdumdd(iro16pg) =dratrawdd(iro16pg)*sc1a+ratraw(iro16pg)*sc1add scfac(iro16pg) = sc1a dscfacdt(iro16pg) = sc1adt dscfacdd(iro16pg) = sc1add c..o17 + p jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(io17),aion(io17),zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(iro17pg) = ratraw(iro17pg) * sc1a dratdumdt(iro17pg) =dratrawdt(iro17pg)*sc1a+ratraw(iro17pg)*sc1adt dratdumdd(iro17pg) =dratrawdd(iro17pg)*sc1a+ratraw(iro17pg)*sc1add scfac(iro17pg) = sc1a dscfacdt(iro17pg) = sc1adt dscfacdd(iro17pg) = sc1add ratdum(iro17pa) = ratraw(iro17pa) * sc1a dratdumdt(iro17pa) =dratrawdt(iro17pa)*sc1a+ratraw(iro17pa)*sc1adt dratdumdd(iro17pa) =dratrawdd(iro17pa)*sc1a+ratraw(iro17pa)*sc1add scfac(iro17pa) = sc1a dscfacdt(iro17pa) = sc1adt dscfacdd(iro17pa) = sc1add c..n14 + a jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(in14),aion(in14),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1add,sc1add) ratdum(irn14ap) = ratraw(irn14ap) * sc1a dratdumdt(irn14ap) =dratrawdt(irn14ap)*sc1a+ratraw(irn14ap)*sc1adt dratdumdd(irn14ap) =dratrawdd(irn14ap)*sc1a+ratraw(irn14ap)*sc1add scfac(irn14ap) = sc1a dscfacdt(irn14ap) = sc1adt dscfacdd(irn14ap) = sc1add c..o18 + p jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(io18),aion(io18),zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(iro18pg) = ratraw(iro18pg) * sc1a dratdumdt(iro18pg) =dratrawdt(iro18pg)*sc1a+ratraw(iro18pg)*sc1adt dratdumdd(iro18pg) =dratrawdd(iro18pg)*sc1a+ratraw(iro18pg)*sc1add scfac(iro18pg) = sc1a dscfacdt(iro18pg) = sc1adt dscfacdd(iro18pg) = sc1add ratdum(iro18pa) = ratraw(iro18pa) * sc1a dratdumdt(iro18pa) =dratrawdt(iro18pa)*sc1a+ratraw(iro18pa)*sc1adt dratdumdd(iro18pa) =dratrawdd(iro18pa)*sc1a+ratraw(iro18pa)*sc1add scfac(iro18pa) = sc1a dscfacdt(iro18pa) = sc1adt dscfacdd(iro18pa) = sc1add c..n15 + a jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(in15),aion(in15),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irn15ap) = ratraw(irn15ap) * sc1a dratdumdt(irn15ap) =dratrawdt(irn15ap)*sc1a+ratraw(irn15ap)*sc1adt dratdumdd(irn15ap) =dratrawdd(irn15ap)*sc1a+ratraw(irn15ap)*sc1add scfac(irn15ap) = sc1a dscfacdt(irn15ap) = sc1adt dscfacdd(irn15ap) = sc1add c..f19 + p jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(if19),aion(if19),zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irf19pa) = ratraw(irf19pa) * sc1a dratdumdt(irf19pa) =dratrawdt(irf19pa)*sc1a+ratraw(irf19pa)*sc1adt dratdumdd(irf19pa) =dratrawdd(irf19pa)*sc1a+ratraw(irf19pa)*sc1add scfac(irf19pa) = sc1a dscfacdt(irf19pa) = sc1adt dscfacdd(irf19pa) = sc1add c..o16 + a jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(io16),aion(io16),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(iro16ap) = ratraw(iro16ap) * sc1a dratdumdt(iro16ap) =dratrawdt(iro16ap)*sc1a+ratraw(iro16ap)*sc1adt dratdumdd(iro16ap) =dratrawdd(iro16ap)*sc1a+ratraw(iro16ap)*sc1add scfac(iro16ap) = sc1a dscfacdt(iro16ap) = sc1adt dscfacdd(iro16ap) = sc1add c..n13 + p jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(in13),aion(in13),zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irn13pg) = ratraw(irn13pg) * sc1a dratdumdt(irn13pg) =dratrawdt(irn13pg)*sc1a+ratraw(irn13pg)*sc1adt dratdumdd(irn13pg) =dratrawdd(irn13pg)*sc1a+ratraw(irn13pg)*sc1add scfac(irn13pg) = sc1a dscfacdt(irn13pg) = sc1adt dscfacdd(irn13pg) = sc1add c..f17 + p jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(if17),aion(if17),zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irf17pg) = ratraw(irf17pg) * sc1a dratdumdt(irf17pg) =dratrawdt(irf17pg)*sc1a+ratraw(irf17pg)*sc1adt dratdumdd(irf17pg) =dratrawdd(irf17pg)*sc1a+ratraw(irf17pg)*sc1add scfac(irf17pg) = sc1a dscfacdt(irf17pg) = sc1adt dscfacdd(irf17pg) = sc1add ratdum(irf17pa) = ratraw(irf17pa) * sc1a dratdumdt(irf17pa) =dratrawdt(irf17pa)*sc1a+ratraw(irf17pa)*sc1adt dratdumdd(irf17pa) =dratrawdd(irf17pa)*sc1a+ratraw(irf17pa)*sc1add scfac(irf17pa) = sc1a dscfacdt(irf17pa) = sc1adt dscfacdd(irf17pa) = sc1add c..o14 + a jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(io14),aion(io14),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(iro14ap) = ratraw(iro14ap) * sc1a dratdumdt(iro14ap) =dratrawdt(iro14ap)*sc1a+ratraw(iro14ap)*sc1adt dratdumdd(iro14ap) =dratrawdd(iro14ap)*sc1a+ratraw(iro14ap)*sc1add scfac(iro14ap) = sc1a dscfacdt(iro14ap) = sc1adt dscfacdd(iro14ap) = sc1add c..f18 + p jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(if18),aion(if18),zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irf18pa) = ratraw(irf18pa) * sc1a dratdumdt(irf18pa) =dratrawdt(irf18pa)*sc1a+ratraw(irf18pa)*sc1adt dratdumdd(irf18pa) =dratrawdd(irf18pa)*sc1a+ratraw(irf18pa)*sc1add scfac(irf18pa) = sc1a dscfacdt(irf18pa) = sc1adt dscfacdd(irf18pa) = sc1add c..o15 + a jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(io15),aion(io15),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(iro15ap) = ratraw(iro15ap) * sc1a dratdumdt(iro15ap) =dratrawdt(iro15ap)*sc1a+ratraw(iro15ap)*sc1adt dratdumdd(iro15ap) =dratrawdd(iro15ap)*sc1a+ratraw(iro15ap)*sc1add scfac(iro15ap) = sc1a dscfacdt(iro15ap) = sc1adt dscfacdd(iro15ap) = sc1add c..triple alpha jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ihe4),aion(ihe4),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ihe4),aion(ihe4),4.0d0,8.0d0, 2 jscr,init,sc2a,sc2adt,sc2add) sc3a = sc1a * sc2a sc3adt = sc1adt * sc2a + sc1a * sc2adt sc3add = sc1add * sc2a + sc1a * sc2add ratdum(ir3a) = ratraw(ir3a) * sc3a dratdumdt(ir3a) = dratrawdt(ir3a)*sc3a + ratraw(ir3a)*sc3adt dratdumdd(ir3a) = dratrawdd(ir3a)*sc3a + ratraw(ir3a)*sc3add scfac(ir3a) = sc3a dscfacdt(ir3a) = sc3adt dscfacdd(ir3a) = sc3add c..o16 to ne20 jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(io16),aion(io16),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(iroag) = ratraw(iroag) * sc1a dratdumdt(iroag) = dratrawdt(iroag)*sc1a + ratraw(iroag)*sc1adt dratdumdd(iroag) = dratrawdd(iroag)*sc1a + ratraw(iroag)*sc1add scfac(iroag) = sc1a dscfacdt(iroag) = sc1adt dscfacdd(iroag) = sc1add c..ne18 to na21 jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 10.0d0,18.0d0,zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irne18ap) = ratraw(irne18ap) * sc1a dratdumdt(irne18ap)= dratrawdt(irne18ap)*sc1a 1 + ratraw(irne18ap)*sc1adt dratdumdd(irne18ap)= dratrawdd(irne18ap)*sc1a 2 + ratraw(irne18ap)*sc1add scfac(irne18ap) = sc1a dscfacdt(irne18ap) = sc1adt dscfacdd(irne18ap) = sc1add c..o15 to ne19 jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(io15),aion(io15),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(iro15ag) = ratraw(iro15ag) * sc1a dratdumdt(iro15ag) =dratrawdt(iro15ag)*sc1a+ratraw(iro15ag)*sc1adt dratdumdd(iro15ag) =dratrawdd(iro15ag)*sc1a+ratraw(iro15ag)*sc1add scfac(iro15ag) = sc1a dscfacdt(iro15ag) = sc1adt dscfacdd(iro15ag) = sc1add c..ne19 + p jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 10.0d0,19.0d0,zion(ih1),aion(ih1), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irne19pg) = ratraw(irne19pg) * sc1a dratdumdt(irne19pg)= dratrawdt(irne19pg)*sc1a 1 + ratraw(irne19pg)*sc1adt dratdumdd(irne19pg)= dratrawdd(irne19pg)*sc1a 1 + ratraw(irne19pg)*sc1add scfac(irne19pg) = sc1a dscfacdt(irne19pg) = sc1adt dscfacdd(irne19pg) = sc1add c..si26 + a jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 14.0d0,26.0d0,zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irsi26ap) = ratraw(irsi26ap) * sc1a dratdumdt(irsi26ap)= dratrawdt(irsi26ap)*sc1a 1 + ratraw(irsi26ap)*sc1adt dratdumdd(irsi26ap)= dratrawdd(irsi26ap)*sc1a 1 + ratraw(irsi26ap)*sc1add scfac(irsi26ap) = sc1a dscfacdt(irsi26ap) = sc1adt dscfacdd(irsi26ap) = sc1add c..ti44 + alpha jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 22.0d0,44.0d0,zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irtiap) = ratraw(irtiap) * sc1a dratdumdt(irtiap) = dratrawdt(irtiap)*sc1a + ratraw(irtiap)*sc1adt dratdumdd(irtiap) = dratrawdd(irtiap)*sc1a + ratraw(irtiap)*sc1add scfac(irtiap) = sc1a dscfacdt(irtiap) = sc1adt dscfacdd(irtiap) = sc1add c..reset the screen initialization flag init = 0 c..end of screen block if end if c..test for 26si decay or 26si (a,p) wp22mg = 1.071d-1 ratdum(irdelta1) = 0.0d0 ratdum(iralam1) = wp22mg dratdumdt(iralam1) = 0.0d0 dratdumdd(iralam1) = 0.0d0 dratdumdy1(iralam1) = 0.0d0 alamt = y(ihe4)*ratdum(irsi26ap) if (ratdum(iralam1) .le. alamt) then ratdum(iralam1) = alamt dratdumdt(iralam1) = y(ihe4)*dratdumdt(irsi26ap) dratdumdd(iralam1) = y(ihe4)*dratdumdd(irsi26ap) dratdumdy1(iralam1) = ratdum(irsi26ap) ratdum(irdelta1) = 1.0d0 end if c..test for 44ti decay or 44ti (a,p) wp30s = 6.667d-2 ratdum(irdelta2) = 0.0d0 ratdum(i