program public_cno include 'implno.dek' include 'const.dek' include 'timers.dek' include 'burn_common.dek' include 'network.dek' c..this program exercises the cno network c..declare integer i,j,nok,nbad double precision tstart,tstep,conserv,tin,din,ein,vin,zin,xin(14), 1 tout,dout,eout,xout(14) c..initialize the network call init_cno 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, 2 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 cno,scno,bcno,dcno 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 cno network call netint(beg,stptry,stpmin,tend,ys2, 1 tol,neqs,nok,nbad,kount,odescal, c 4 cno,scno,bcno,forder_ma28) c 4 cno,scno,bcno,forder_umf) c 4 cno,scno,bcno,forder_y12m) c 4 cno,dcno,bcno,forder_ludcmp) c 4 cno,dcno,bcno,forder_leqs) c 4 cno,dcno,bcno,forder_lapack) c 4 cno,dcno,bcno,forder_gift) c 4 cno,scno,bcno,forder_biconj) c 4 cno,scno,bcno,rosen_ma28) c 4 cno,scno,bcno,rosen_umf) c 4 cno,scno,bcno,rosen_y12m) c 4 cno,dcno,bcno,rosen_ludcmp) c 4 cno,dcno,bcno,rosen_leqs) c 4 cno,dcno,bcno,rosen_lapack) c 4 cno,dcno,bcno,rosen_gift) c 4 cno,scno,bcno,rosen_biconj) 4 cno,scno,bcno,stifbs_ma28) c 4 cno,scno,bcno,stifbs_umf) c 4 cno,scno,bcno,stifbs_y12m) c 4 cno,dcno,bcno,stifbs_ludcmp) c 4 cno,dcno,bcno,stifbs_leqs) c 4 cno,dcno,bcno,stifbs_lapack) c 4 cno,dcno,bcno,stifbs_gift) c 4 cno,scno,bcno,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 cno network c..routine cno sets up the odes c..routine rhs evaluates the right hand sides c..routine dcno sets up the dense cno jacobian c..routine bcno build the nonzero locations for scno c..routine scno sets up the sparse cno jacobian c..routine cnorat generates the reaction rates for routine cno c..routine cnotab generates the raw rates from table interpolation c..routine screen_cno applies screening corrections to the raw rates c..routine init_cno initializes the cno network subroutine cno(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 the 4 cno cycles. c.. c..isotopes cycle1: h1, he4, c12, n13, c13, n14, o15, n15 c cycle2: o16, f17, o17 c cycle3: f18, o18 c cycle4: f19 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, 1 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 cnotab else call cnorat end if c..do the screening here because the corrections depend on the composition call screen_cno(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; no o14 here sneutcno = snucno(y(in13),bion(ic13),bion(in13), 1 0.0d0,0.0d0,0.0d0, 2 y(io15),bion(in15),bion(io15), 3 y(if17),bion(io17),bion(if17), 4 y(if18),bion(io18),bion(if18)) c..sum 'em c 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 c..zero the odes do i=1,neqs dydt(i) = 0.0d0 enddo c..set up the system of odes: 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) 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) 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) 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) 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) 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..o15 reactions dydt(io15) = y(in14)*y(ih1)*rate(irn14pg) 1 - y(io15)*rate(iro15gp) 2 - y(io15)*rate(iro15enu) 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) 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) c..f18 reactions dydt(if18) = y(io17)*y(ih1)*rate(iro17pg) 1 - y(if18)*rate(irf18gp) 2 - y(if18)*rate(irf18enu) 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) return end subroutine dcno(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 cno jacobian c..declare the pass integer nlog,nphys double precision tt,y(1),dfdy(nphys,nphys) c..local variables integer i,j double precision zbarxx,ytot1,abar,zbar,ye,taud,taut, 1 snuda,snudz,sum1,sum2,zz,xx double precision conv parameter (conv = ev2erg*1.0d6*avo) 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 cnotab else call cnorat end if c..do the screening here because the corrections depend on the composition call screen_cno(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) dfdy(ih1,ihe4) = y(ic12)*ratdum(irc12ap) 1 + y(in14)*ratdum(irn14ap) 2 + y(in15)*ratdum(irn15ap) 3 + y(io16)*ratdum(iro16ap) dfdy(ih1,ic12) = -y(ih1)*ratdum(irc12pg) 1 + y(ihe4)*ratdum(irc12ap) dfdy(ih1,ic13) = -y(ih1)*ratdum(irc13pg) dfdy(ih1,in13) = ratdum(irn13gp) 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,io15) = ratdum(iro15gp) 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) dfdy(ih1,if18) = ratdum(irf18gp) dfdy(ih1,if19) = ratdum(irf19gp) 1 - y(ih1)*ratdum(irf19pa) 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) dfdy(ihe4,ihe4) = -y(ic12)*ratdum(irc12ap) 1 - y(in14)*ratdum(irn14ap) 2 - y(in15)*ratdum(irn15ap) 3 - y(io16)*ratdum(iro16ap) dfdy(ihe4,ic12) = -y(ihe4)*ratdum(irc12ap) dfdy(ihe4,in14) = -y(ihe4)*ratdum(irn14ap) dfdy(ihe4,in15) = y(ih1)*ratdum(irn15pa) dfdy(ihe4,io16) = -y(ihe4)*ratdum(iro16ap) dfdy(ihe4,io17) = y(ih1)*ratdum(iro17pa) dfdy(ihe4,io18) = y(ih1)*ratdum(iro18pa) dfdy(ihe4,if19) = y(ih1)*ratdum(irf19pa) c..c12 jacobian elements dfdy(ic12,ih1) = -y(ic12)*ratdum(irc12pg) 1 + y(in15)*ratdum(irn15pa) dfdy(ic12,ihe4) = -y(ic12)*ratdum(irc12ap) dfdy(ic12,ic12) = -y(ih1)*ratdum(irc12pg) 1 - y(ihe4)*ratdum(irc12ap) dfdy(ic12,in13) = ratdum(irn13gp) dfdy(ic12,in15) = y(ih1)*ratdum(irn15pa) 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) dfdy(in13,ic12) = y(ih1)*ratdum(irc12pg) dfdy(in13,in13) = -ratdum(irn13gp) 1 - ratdum(irn13enu) 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,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..o15 jacobian elements dfdy(io15,ih1) = y(in14)*ratdum(irn14pg) dfdy(io15,in14) = y(ih1)*ratdum(irn14pg) dfdy(io15,io15) = -ratdum(iro15gp) 1 - ratdum(iro15enu) 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) dfdy(io16,in15) = y(ih1)*ratdum(irn15pg) dfdy(io16,io16) = -ratdum(iro16gp) 1 - y(ih1)*ratdum(iro16pg) 2 - y(ihe4)*ratdum(iro16ap) 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) dfdy(if17,io16) = y(ih1)*ratdum(iro16pg) dfdy(if17,if17) = -ratdum(irf17gp) 1 - ratdum(irf17enu) c..f18 jacobian elements dfdy(if18,ih1) = y(io17)*ratdum(iro17pg) dfdy(if18,io17) = y(ih1)*ratdum(iro17pg) dfdy(if18,if18) = -ratdum(irf18gp) 1 - ratdum(irf18enu) 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) 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 bcno(iloc,jloc,nzo,np) include 'implno.dek' include 'network.dek' c..this routine builds the nonzero matrix locations for scno c..input is the integer arrys iloc and jloc, both of dimension np, that c..on output contain nzo matrix element locations. c..declare the pass integer np,iloc(np),jloc(np),nzo c..local variables integer i c..communicate with scno integer neloc parameter (neloc=158) integer eloc(neloc),nterms common /elccno/ eloc,nterms c..initialize nterms = 0 nzo = 0 do i=1,neloc eloc(i) = 0 enddo call tree_init(neqs) c..tag the nonzero locations 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,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) 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,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,if19,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) 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) 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,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..o15 jacobian elements call tree(io15,ih1,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) 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,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,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if17,if17,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,io17,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(if18,if18,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) c..if we are doing a pure network, we are done if (pure_network .eq. 1) return c..temperature elements do i=1,ionmax call tree(i,itemp,eloc,neloc,nterms,nzo,iloc,jloc,np) end do c..density elements 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 ddensity 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 scno(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 cno 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..local variables integer i,nt,iat double precision a1,a2,a3,a4,zz, 1 zbarxx,ytot1,abar,zbar,ye,taud,taut, 2 snuda,snudz double precision conv parameter (conv = ev2erg*1.0d6*avo) c..communicate with the jacobian builder integer neloc parameter (neloc=158) integer eloc(neloc),nterms common /elccno/ 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 cnotab else call cnorat end if c..do the screening here because the corrections depend on the composition call screen_cno(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) 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) 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) 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(o15) a1 = ratdum(iro15gp) 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) 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) 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..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) 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) 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) 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(o16) a1 = -y(ihe4)*ratdum(iro16ap) 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(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..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) 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) 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..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) 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) 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..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(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..o15 jacobian elements c..d(o15)/d(h1) a1 = y(in14)*ratdum(irn14pg) 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(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) 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..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) 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(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) 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) 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(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) 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..f18 jacobian elements c..d(f18)/d(h1) a1 = y(io17)*ratdum(iro17pg) 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(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) 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..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..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 ddensity 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 scno: nt .ne. nterms' stop 'error in routine scno' end if return end subroutine cnorat include 'implno.dek' include 'burn_common.dek' include 'network.dek' c.. c..this routine generates nuclear reaction rates for the cno 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..shut down various reactions c..these are for cno cycle 1 c ratraw(irc12pg) = 0.0d0 c ratraw(irn13gp) = 0.0d0 c ratraw(irn13enu) = 0.0d0 c ratraw(irc13pg) = 0.0d0 c ratraw(irn14gp) = 0.0d0 c ratraw(irn14pg) = 0.0d0 c ratraw(iro15gp) = 0.0d0 c ratraw(iro15enu) = 0.0d0 c ratraw(irn15pa) = 0.0d0 c ratraw(irc12ap) = 0.0d0 c..these are cno cycle 2 c ratraw(irn15pg) = 0.0d0 c ratraw(iro16gp) = 0.0d0 c ratraw(iro16pg) = 0.0d0 c ratraw(irf17gp) = 0.0d0 c ratraw(irf17enu) = 0.0d0 c ratraw(iro17pa) = 0.0d0 c ratraw(irn14ap) = 0.0d0 c..these are for cno cycle 3 c ratraw(iro17pg) = 0.0d0 c ratraw(irf18gp) = 0.0d0 c ratraw(irf18enu) = 0.0d0 c ratraw(iro18pa) = 0.0d0 c ratraw(irn15ap) = 0.0d0 c..these are for cno cyle 4 c ratraw(iro18pg) = 0.0d0 c ratraw(irf19gp) = 0.0d0 c ratraw(irf19pa) = 0.0d0 c ratraw(iro16ap) = 0.0d0 return end subroutine cnotab 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 c..use 120 points per decade imax = 481 if (imax .gt. nrattab) stop 'imax too small in cnotab' tlo = 6.0d0 thi = 10.0d0 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 cnorat 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 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 return end subroutine screen_cno(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,init double precision y(*),sc1a,sc1adt,sc1add, 1 abar,zbar,z2bar,ytot1,zbarxx,z2barxx 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 off if (screen_on .eq. 0) return c..screening is on 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)*sc1add 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,sc1adt,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 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,sc1adt,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..reset the screen initialization flag init = 0 return end subroutine init_cno include 'implno.dek' include 'network.dek' c.. c..this routine initializes stuff for the cno network c.. c..declare integer i c..for easy zeroing of the isotope pointers integer isotp(nisotp) equivalence (isotp(1),ih1) c..for easy zeroing of the rate pointers integer rts(numrates) equivalence (rts(1),ir3a) c..zero all the isotope pointers do i=1,nisotp isotp(i) = 0 enddo c..zero all the rate pointers do i=1,numrates rts(i) = 0 enddo c..set the size of the network and the number of rates idnet = idcno ionmax = 14 iener = ionmax + 1 itemp = ionmax + 2 iden = ionmax + 3 neqs = iden nrat = 26 netname = 'cno' c..set the id numbers of the elements ih1 = 1 ihe4 = 2 ic12 = 3 ic13 = 4 in13 = 5 in14 = 6 in15 = 7 io15 = 8 io16 = 9 io17 = 10 io18 = 11 if17 = 12 if18 = 13 if19 = 14 c..set the names of the elements ionam(ih1) = 'h1 ' ionam(ihe4) = 'he4 ' ionam(ic12) = 'c12 ' ionam(ic13) = 'c13 ' ionam(in13) = 'n13 ' ionam(in14) = 'n14 ' ionam(in15) = 'n15 ' ionam(io15) = 'o15 ' ionam(io16) = 'o16 ' ionam(io17) = 'o17 ' ionam(io18) = 'o18 ' ionam(if17) = 'f17 ' ionam(if18) = 'f18 ' ionam(if19) = 'f19 ' ionam(iener) = 'ener ' ionam(itemp) = 'temp ' ionam(iden) = 'den ' c..set the number of nucleons in the element aion(ih1) = 1.0d0 aion(ihe4) = 4.0d0 aion(ic12) = 12.0d0 aion(ic13) = 13.0d0 aion(in13) = 13.0d0 aion(in14) = 14.0d0 aion(in15) = 15.0d0 aion(io15) = 15.0d0 aion(io16) = 16.0d0 aion(io17) = 17.0d0 aion(io18) = 18.0d0 aion(if17) = 17.0d0 aion(if18) = 18.0d0 aion(if19) = 19.0d0 c..set the number of protons in the element zion(ih1) = 1.0d0 zion(ihe4) = 2.0d0 zion(ic12) = 6.0d0 zion(ic13) = 6.0d0 zion(in13) = 7.0d0 zion(in14) = 7.0d0 zion(in15) = 7.0d0 zion(io15) = 8.0d0 zion(io16) = 8.0d0 zion(io17) = 8.0d0 zion(io18) = 8.0d0 zion(if17) = 9.0d0 zion(if18) = 9.0d0 zion(if19) = 9.0d0 c..set the number of neutrons do i=1,ionmax nion(i) = aion(i) - zion(i) enddo c..set the binding energy of the element bion(ih1) = 0.0d0 bion(ihe4) = 28.2928d0 bion(ic12) = 92.1624d0 bion(ic13) = 97.1088d0 bion(in13) = 94.1064d0 bion(in14) = 104.6598d0 bion(in15) = 115.4932d0 bion(io15) = 111.9558d0 bion(io16) = 127.6202d0 bion(io17) = 131.7636d0 bion(io18) = 139.8080d0 bion(if17) = 128.2212d0 bion(if18) = 137.3706d0 bion(if19) = 147.8020d0 c..set the partition functions - statistical weights, ground-state only here do i=1,ionmax wpart(i) = 1.0d0 enddo c..set the id numbers of the reaction rates irc12pg = 1 irn13gp = 2 irn13enu = 3 irc13pg = 4 irn14gp = 5 irn14pg = 6 iro15gp = 7 iro15enu = 8 irn15pa = 9 irc12ap = 10 irn15pg = 11 iro16gp = 12 iro16pg = 13 irf17gp = 14 irf17enu = 15 iro17pa = 16 irn14ap = 17 iro17pg = 18 irf18gp = 19 irf18enu = 20 iro18pa = 21 irn15ap = 22 iro18pg = 23 irf19gp = 24 irf19pa = 25 iro16ap = 26 c..set the names of the reaction rates ratnam(irc12pg) = 'c12(p,g)n12' ratnam(irn13gp) = 'n13(g,p)c12' ratnam(irn13enu) = 'n13(p=>n)c13' ratnam(irc13pg) = 'c13(p,g)n14' ratnam(irn14gp) = 'n14(g,p)c13' ratnam(irn14pg) = 'n14(p,g)o15' ratnam(iro15gp) = 'o15(g,p)n14' ratnam(iro15enu) = 'o15(p=>n)n15' ratnam(irn15pa) = 'n15(p,a)c12' ratnam(irc12ap) = 'c12(a,p)n15' ratnam(irn15pg) = 'n15(p,g)o16' ratnam(iro16gp) = 'o16(g,p)n15' ratnam(iro16pg) = 'o16(p,g)f17' ratnam(irf17gp) = 'f17(g,p)o16' ratnam(irf17enu) = 'f17(p=>n)o17' ratnam(iro17pa) = 'o17(p,a)n14' ratnam(irn14ap) = 'n14(a,p)o17' ratnam(iro17pg) = 'o17(p,g)f18' ratnam(irf18gp) = 'f18(g,p)o17' ratnam(irf18enu) = 'f18(p=>n)o18' ratnam(iro18pa) = 'o18(p,a)n15' ratnam(irn15ap) = '15n(a,p)o18' ratnam(iro18pg) = 'o18(p,g)f19' ratnam(irf19gp) = 'f19(g,p)o18' ratnam(irf19pa) = 'f19(p,a)o16' ratnam(iro16ap) = 'o18(a,p)f19' return end c include '../src_method/netint.f' c-------------------------------------------------------------------- subroutine netint(start,stptry,stpmin,stopp,bc, 1 eps,ylogi,nok,nbad,kount,odescal, 2 derivs,jakob,bjakob,steper) include 'implno.dek' include 'burn_common.dek' include 'network.dek' c..ode integrator for stiff odes c..tuned for nnuclear reacton networks c..input: c..start = beginning integration point c..stptry = suggested first step size c..stpmin = minimum allowable step size c..stopp = ending integration point c..bc = initial conditions, array of physical dimension yphys c..eps = desired fraction error during the integration c..odescal = error scaling factor c..derivs = name of the routine that contains the odes c..jakob = name of the routine that contains the jacobian of the odes c..bjakob = name of the routine that sets the pointers of the sparse jacobian c..steper = name of the routine that will take a single step c..output: c..nok = number of succesful steps taken c..nbad = number of bad steps taken, bad but retried and then succesful c..kount = total number of steps taken c..declare the pass external derivs,jakob,bjakob,steper integer ylogi,nok,nbad,kount double precision start,stptry,stpmin,stopp,bc(ylogi),eps, 1 odescal c..for communicating a root find c..common block communication double precision nse_temp_switch common /nsetsw/ nse_temp_switch c..local variables character*5 cdtname integer nmax,stpmax,i,ii,nstp,idt parameter (nmax = abignet*nzmax, stpmax=200000) double precision yscal(nmax),y(nmax),dydx(nmax),xdum(nmax), 1 sum,cons,t9,tau_nse,tau_qse, 1 x,h,hdid,hnext,tiny parameter (tiny=1.0d-15) c..for a root find on the trajectory external time_switch double precision up_zbrent,time_switch,xb,tol_switch parameter (tol_switch = 1.0e-12) integer niter c..for smooth plot timesteps double precision ratio,xfloor,ychangemax,ynew,yold,yy,dtx c..for some more informative printouts double precision zbarxx,ytot1,abar,zbar,ye,ttz,ddz, 1 ttz1,ddz1,ttz2,ddz2,tlo,thi c..for nse character*3 cmode integer igues,nse_switch double precision xmun,xmup c..here are the format statements for printouts as we integrate 100 format(1x,i6,' ',a,a,1pe11.4,a,a,a,1pe11.4, 1 3(a6,1pe10.3),5(a5,1pe9.2)) 101 format(1x,1p12e10.2) c..initialize if (ylogi .gt. nmax) stop 'ylogi > nmax in routine netint' x = start h = sign(stptry,stopp-start) nok = 0 nbad = 0 kount = 0 idt = 0 nse_temp_switch = 10.0d9 thi = nse_temp_switch*(1.0d0 + tol_switch) tlo = nse_temp_switch*(1.0d0 - tol_switch) cmode = ' ' c..store the first step do i=1,ylogi y(i) = bc(i) enddo c..take at most stpmax steps do nstp=1,stpmax c..positive definite abundance fractions do i=1,ionmax y(i) = min(1.0d0, max(y(i),1.0d-30)) enddo c..form the mass fractions and nonconservation sum = 0.0d0 do i=1,ionmax xdum(i) = y(i) * aion(i) end do do i=1,ionmax sum = sum + xdum(i) end do cons = 1.0d0 - sum c..renorm the abundances c sum = 1.0d0/sum c do i=1,ionmax c xdum(i) = sum * xdum(i) c end do c do i=1,ionmax c y(i) = min(1.0d0,max(xdum(i)/aion(i),1.0d-30)) c end do c..get the right hand sides c if (nse_on .eq. 0) then call derivs(x,y,dydx) c..scaling vector used to monitor accuracy do i=1,ylogi yscal(i) = max(odescal,abs(y(i))) enddo c end if c..store intermediate results kount = kount+1 c..detailed file print if (iprint_files .eq. 1) call net_output(kount,x,y,derivs) c..screen print if (iprint_screen .eq. 1) then call azbar(xdum,aion,zion,ionmax, 1 zwork1,abar,zbar) ye = zbar/abar ttz = -1.0d0 ddz = -1.0d0 if (pure_network .eq. 0) then ttz = y(itemp) ddz = y(iden) else ttz = btemp ddz = bden end if if (trho_hist) call update2(x,ttz,ddz) if (idt .eq. 0) then cdtname = 'time ' else cdtname = ionam(idt) end if call indexx(ionmax,xdum,izwork1) write(6,100) kount,cmode,' time',x, 1 ' dt(',cdtname,')',hdid,(ionam(izwork1(ii)), 2 xdum(izwork1(ii)),ii=ionmax,ionmax-2,-1), 3 ' temp',ttz,' den ',ddz,' enuc',sdot-sneut, 4 ' ye ',ye,' sum ',cons c read(5,*) end if c call flush(6) c call flush_(6) c..if the step can overshoot the stop point cut it if ((x+h-stopp)*(x+h-start) .gt. 0.0d0) h = stopp - x c.. do an nse step - this should now be made a subroutine c..if nse evolution is allowed if (allow_nse_evol .eq. 1) then c..get the thermodymanic conditions if (pure_network .eq. 0) then ttz = y(itemp) ddz = y(iden) else ttz = btemp ddz = bden end if c if (trho_hist) call update2(x+h,ttz,ddz) c..if we are interpolating a trajectory c..get the values at the present time point and the suggested next time point if (trho_hist) then call update2(x,ttz1,ddz1) call update2(x+h,ttz2,ddz2) c..if both are above the nse point if (ttz1 .ge. nse_temp_switch .and. 1 ttz2 .ge. nse_temp_switch) then ttz = ttz2 ddz = ddz2 xb = 0.0d0 nse_switch = 0 c write(6,*) 'both above' c..if we are falling out of nse else if (ttz1 .ge. thi .and. 1 ttz2 .le. tlo) then xb = up_zbrent(time_switch,x,x+h,tol_switch,niter) h = max(xb - x,tol_switch) call update2(xb,ttz,ddz) nse_switch = 1 c write(6,*) 'ttz2 below',ttz1.ge.thi,ttz2.le.tlo c..if we are going into nse else if (ttz1 .le. tlo .and. 1 ttz2 .ge. thi) then xb = up_zbrent(time_switch,x,x+h,tol_switch,niter) h = max(xb - x,tol_switch) call update2(xb,ttz,ddz) nse_switch = 1 c write(6,*) 'ttz2 above',ttz1.le.tlo,ttz2.ge.thi c..if we are out of nse, then these values get reset c..this also applies if we are within tol_switch of nse_temp_switch else ttz = ttz2 ddz = ddz2 xb = 0.0d0 nse_switch = 0 c write(6,*) 'both below' end if end if c write(6,119) x,xb,h,ttz,ddz c write(6,119) ttz1,ttz2,ttz,tlo,thi c 119 format(1x,1p6e24.16) c read(5,*) c t9 = ttz * 1.0d-9 c tau_nse = ddz**(0.2d0) * exp(179.7d0/t9 - 40.5d0) c tau_qse = exp(149.7d0/t9 - 39.15d0) c..initialize for what type of integration cmode = 'int' nse_on = 0 c..check for nse conditions c if (ttz .ge. 10.0e9 .and. h.gt.tau_nse .and. x.gt.tau_nse) then if (ttz .ge. nse_temp_switch ) then cmode = 'nse' nse_on = 1 call azbar(xdum,aion,zion,ionmax, 1 zwork1,abar,zbar) ye = zbar/abar igues = 1 call nse(ttz,ddz,ye,igues,1,xdum,xmun,xmup,0) c..claim we did the requested time step x = x + h hdid = h c..estimate the next time step hnext = 1.0e30 ratio = 1.0d30 xfloor = 1.0e-5 ychangemax = 0.10d0 if (kount .ne. 1) then do i=1,ionmax ynew = xdum(i)/aion(i) yy = abs(ynew - y(i)) if (yy*ratio .gt. y(i) .and. y(i) .ge. xfloor) then ratio=y(i)/yy idt = i end if enddo end if hnext = min(ratio*h*ychangemax,2.0d0*h) if (nse_switch .eq. 1) then hnext = max(2.0d0*tol_switch,1.0d-2*hnext) end if if (hnext .eq. 2.0d0*h) idt = 0 c write(6,119) hnext c..update the abundance vector do i = 1,ionmax y(i) = xdum(i)/aion(i) end do c..end of nse if end if end if c..do an integration step if (nse_on .eq. 0) then call steper(y,dydx,ylogi,x,h,eps,yscal,hdid,hnext, 1 derivs,jakob,bjakob,nstp,idt) end if if (hdid.eq.h) then nok = nok+1 else nbad = nbad+1 end if c..this is the normal exit point, save the final step if ( (nstp .eq. stpmax) .or. 1 (x-stopp)*(stopp-start).ge. 0.0d0 .or. 3 (psi .ge. 1.0 .and. y(itemp) .lt. temp_stop) .or. 4 (psi .le. -1.0 .and. y(itemp) .gt. temp_stop) .or. 5 (detonation .and. y(iden) .lt. den_stop) .or. c 5 (detonation .and. abs(1.0d0-cs_cj/y(ivelx)).lt.1.0e-4) .or. 6 (y(id_stop)*aion(id_stop) .lt. xmass_stop) ) then c write(6,*) 'bailing' c write(6,*) id_stop,y(id_stop),aion(id_stop),xmass_stop c write(6,*) y(itemp),temp_stop c write(6,*) stopp do i=1,ylogi bc(i) = y(i) enddo kount = kount+1 c..detailed file print if (iprint_files .eq. 1) call net_output(kount,x,y,derivs) c..screen print if (iprint_screen .eq. 1) then call azbar(xdum,aion,zion,ionmax, 1 zwork1,abar,zbar) ye = zbar/abar ttz = -1.0d0 ddz = -1.0d0 if (pure_network .eq. 0) then ttz = y(itemp) ddz = y(iden) else ttz = btemp ddz = bden end if if (trho_hist) call update2(x,ttz,ddz) call indexx(ionmax,xdum,izwork1) write(6,100) kount,cmode,' time',x, 1 ' dt(',cdtname,')',hdid,(ionam(izwork1(ii)), 2 xdum(izwork1(ii)),ii=ionmax,ionmax-2,-1), 3 ' temp',ttz,' den ',ddz,' enuc',sdot-sneut, 4 ' ye ',ye,' sum ',cons end if c call flush(6) c call flush_(6) return end if c..set the step size for the next iteration; stay above stpmin c dtx = 1.0e30 c c if (kount .ne. 1) then c ratio = 1.0d30 c xfloor = 1.0e-5 c ychangemax = 0.30d0 c do i=1,ionmax c ynew = max(y(i),1.0e-20) c yold = yrk(i,kount-1) c yy = abs(ynew - yold) c if (yy*ratio .gt. yold .and. yold .ge. xfloor) ratio=yold/yy c enddo c dtx = min(ratio*h*ychangemax,2.0d0*h) c end if c c h = min(hnext,dtx) c..limit timestep changes to a factor of two c h = min(hnext,2.0d0*h) h = hnext if (abs(h).lt.stpmin) stop 'h < stpmin in netint' c..back for another iteration or death enddo stop 'more than stpmax steps required in netint' end c include '../src_net/net_input.f ' c-------------------------------------------------------------------- subroutine net_input(tstart,tstep,tin,din,vin,zin,ein,xin) include 'implno.dek' include 'const.dek' include 'vector_eos.dek' include 'burn_common.dek' include 'network.dek' include 'cjdet.dek' c..declare the pass double precision tstart,tstep,tin,din,vin,zin,ein,xin(*) c..local variables integer i,j,k,ibtype,ictype,igues,kkase,ians double precision xneut,xh1,xhe4,xc12,xn14,xo16,xne20,xne22,xsi28, 1 xfe52,xfe54,xni56,zye,sum,abar,zbar,ye_orig,xmup, 2 xmun,qdum,a,z,xelem,andgrev c..bigbang specifics double precision fac,f1,zeta3 parameter (zeta3 = 1.20205690315732d0) c..popular format statements 01 format(1x,a,a,a) 02 format(1x,a,'=',1pe10.3,' ',a,'=',1pe10.3,' ', 1 a,'=',1pe10.3,' ',a,'=',1pe10.3,' ', 2 a,'=',1pe10.3) 03 format(a) 04 format(1x,a,'=',i2,' ',a,'=',i2,' ', 1 a,'=',i2,' ',a,'=',i2,' ', 2 a,'=',i2) c..initialize the common block variables call net_initialize c..inititailize local variables ibtype = 0 ictype = 0 tstart = 0.0d0 tstep = 0.0d0 bpres = 0.0d0 tin = 0.0d0 din = 0.0d0 vin = 0.0d0 zin = 0.0d0 zye = 0.0d0 do i=1,ionmax xin(i) = 1.0d-30 end do c--------------------------------------------------------------------------- c..get the burn type 10 write(6,01) 'give burning mode:' write(6,01) ' ibtype = 0 = stop' write(6,01) ' 1 = onestep' write(6,01) ' 2 = hydrostatic' write(6,01) ' 3 = expansion' write(6,01) ' 4 = self heat constant density' write(6,01) ' 5 = self heat constant pressure' write(6,01) ' 6 = big bang ' write(6,01) ' 7 = detonation' write(6,01) ' 8 = temp-den history' read(5,*) ibtype if (ibtype .lt. 0 .or. ibtype .gt. 8) goto 10 c..general options 11 write(6,*) write(6,04) 'set general options:' write(6,04) 'screen_on',screen_on write(6,04) 'use_tables',use_tables write(6,04) 'weak_on',weak_on write(6,04) 'ffn_on',ffn_on write(6,04) 'pure_network',pure_network write(6,04) 'nse_analysis',nse_analysis write(6,04) 'allow_nse_evol',allow_nse_evol write(6,04) 'iprint_files',iprint_files write(6,04) 'iprint_screen',iprint_screen write(6,02) 'sthreshold',sthreshold,' set > 1 to disable' write(6,*) write(6,01) 'if these are ok, enter 1, otherwise enter 0 =>' read(5,*) ians if (ians .lt. 0 .or. ians .gt. 1) goto 11 if (ians .eq. 0) then 12 write(6,01) 'give the 9 integer and one real vector =>' read(5,*) screen_on, use_tables, weak_on, ffn_on, 1 pure_network, nse_analysis, allow_nse_evol, 2 iprint_files, iprint_screen, 3 sthreshold if (screen_on .lt. 0 .or. screen_on .gt. 1) goto 12 if (use_tables .lt. 0 .or. use_tables .gt. 1) goto 12 if (weak_on .lt. 0 .or. weak_on .gt. 1) goto 12 if (ffn_on .lt. 0 .or. ffn_on .gt. 1) goto 12 if (pure_network .lt. 0 .or. pure_network .gt. 1) goto 12 if (nse_analysis .lt. 0 .or. nse_analysis .gt. 1) goto 12 if (iprint_files .lt. 0 .or. iprint_files .gt. 1) goto 12 if (iprint_screen .lt. 0 .or. iprint_screen .gt. 1) goto 12 goto 11 end if c..get the thermodynamics if (ibtype .eq. 5) then write(6,01) 'give the ending time, temperature, pressure =>' read (5,*) tstep,tin,bpres else if (ibtype .eq.6) then write(6,01) 'give the ending time, temperature =>' read (5,*) tstep,tin,din else if (ibtype .ne. 8) then write(6,01) 'give the ending time, temperature, density =>' read (5,*) tstep,tin,din end if c..get the composition c if (ibtype .ne. 8) then 20 write(6,01) 'give initial composition:' write(6,01) ' ictype = 0 = leave alone; read from file' write(6,01) ' 1 = solar abundances' write(6,01) ' 2 = nse' write(6,01) ' 3 = specify initial composition' read(5,*) ictype if (ictype .lt. 0 .or. ictype .gt. 3) goto 20 if (ictype .eq. 3) then write(6,01) 1 'n h1 he4 c12 n14 o16 ne20 ne22 si28 fe52 fe54 ni56 =>' read(5,*) xneut,xh1,xhe4,xc12,xn14,xo16,xne20,xne22,xsi28, 1 xfe52,xfe54,xni56 end if c end if c..get the output root file name write(6,*) ' ' write(6,01) 'give output root name, for default "foo_"=>' read(5,03) hfile if (hfile(1:2) .eq. ' ') hfile = 'foo_' c--------------------------------------------------------------------------- c..crunch the users input c..set the burn type logical if (ibtype .eq. 0) then stop 'normal termination' else if (ibtype .eq. 1) then one_step = .true. else if (ibtype .eq. 2) then hydrostatic = .true. else if (ibtype .eq. 3) then expansion = .true. c..psi = 1 is an adiabatic expansion c..psi = -1 is an adiabatic implosion psi = 1.0d0 c psi = -1.0d0 den0 = din temp0 = tin temp_stop = 1.0d7 c temp_stop = 1.0d10 if ( (psi .ge. 1.0 .and. temp_stop .ge. tin) .or. 1 (psi .le. -1.0 .and. temp_stop .le. tin)) 2 stop 'bad adiabatic temp_stop in routine burner' else if (ibtype .eq. 4) then self_heat_const_den = .true. else if (ibtype .eq. 5) then self_heat_const_pres = .true. else if (ibtype .eq. 6) then bbang = .true. c..should probably read these parameters as input eta1 = 4.0e-10 xnnu = 3.0d0 hubble = 65.0d0 cmbtemp = 2.73d0 c..set the initial n and p abundances; equation 3 of wagoner et al 1967 fac = exp((mn - mp)*clight**2/(kerg*tin)) xneut = 1.0d0/(1.0d0 + fac) xh1 = 1.0d0 - xneut c..set the density from the temperature and eta1 f1 = 30.0d0 * zeta3/pi**4 * asol/(kerg*avo) din = f1 * eta1 * tin**3 c..detonation else if (ibtype .eq. 7) then detonation = .true. c..thermodynamic profile being given else if (ibtype .eq. 8) then trho_hist = .true. write(6,*) 'give the trho_history file =>' read(5,03) trho_file c..oops else stop 'unknown burn type' end if c--------------------------------------------------------------------------- c..read the thermodynamic trajectory and initial abundances c..transfer the info stored in xsum and zsum from the update2 call if (trho_hist) then call update2(tstart,tin,din) do i=1,ionmax xin(i) = xsum(i) enddo tstart = zsum(1) c tstart = 0.1685d0 tstep = zsum(2) zye = ysum(1) end if c..massage the input composition, includes possible changes to the c..the abundances read in from the trho_hist file c..solar abundances if (ictype .eq. 1) then do i=1,ionmax xin(i) = andgrev(ionam(i),z,a,xelem) enddo if (iprot .ne. 0) xin(iprot) = andgrev('h1 ',z,a,xelem) c..put it in nse else if (ictype .eq. 2) then if (zye .eq. 0.0) zye = 0.5d0 igues = 1 call nse(tin,din,zye,igues,1,xin,xmun,xmup,0) c..set the composition variables else if (ictype .eq. 3) then if (ineut .ne. 0) xin(ineut) = xneut if (ih1 .ne. 0) xin(ih1) = xh1 if (iprot .ne. 0) xin(iprot) = xh1 if (ih1 .ne. 0 .and. iprot .ne. 0) xin(iprot) = 0.0d0 if (ihe4 .ne. 0) xin(ihe4) = xhe4 if (ic12 .ne. 0) xin(ic12) = xc12 if (in14 .ne. 0) xin(in14) = xn14 if (io16 .ne. 0) xin(io16) = xo16 if (ine20 .ne. 0) xin(ine20) = xne20 if (ine22 .ne. 0) xin(ine22) = xne22 if (isi28 .ne. 0) xin(isi28) = xsi28 if (ife52 .ne. 0) xin(ife52) = xfe52 if (ife54 .ne. 0) xin(ife54) = xfe54 if (ini56 .ne. 0) xin(ini56) = xni56 end if c..write out the input composition so far c write(6,02) (ionam(i),xin(i), i=1,ionmax) c read(5,*) c..normalize the composition do i=1,ionmax xin(i) = min(1.0d0,max(xin(i),1.0d-30)) end do sum = 0.0d0 do i=1,ionmax sum = sum + xin(i) enddo sum = 1.0d0/sum do i=1,ionmax xin(i) = min(1.0d0,max(xin(i) * sum,1.0d-30)) enddo c--------------------------------------------------------------------------- c..get the ye of the initial compositon call azbar(xin,aion,zion,ionmax, 1 zwork1,abar,zbar) ye_orig = zbar/abar c..modify the composition if ye_orig is less than 0.55 c if (ye_orig .le. 0.55) then c c..set the mass fraction of fe58 to set the desired ye c ye_want = 0.495d0 c ye_want = 0.50d0 c if (ye_want .eq. 0.5) then c xin(ife58) = 0.0d0 c else c xin(ife58) = (ye_orig - ye_want) / c 1 (ye_orig - zion(ife58)/aion(ife58)) c end if c c..reset the mass fractions of everything else c sum = 1.0d0 - xin(ife58) c do i=1,ionmax c if (i .ne. ife58) xin(i) = xin(i) * sum c enddo c end if c--------------------------------------------------------------------------- c..modify for a detonation c..get the chapman-jouget solution if (detonation) then kkase = 1 mach = 0.0d0 do i=1,ionmax xmass_up(i) = xin(i) enddo temp_up = tin den_up = din call cjsolve(kkase,xmass_up,temp_up,den_up,mach, 1 qburn_cj,xmass_cj,ener_up,pres_up,cs_up, 2 vel_det,vel_cj,temp_cj,den_cj,ener_cj,pres_cj,cs_cj) write(6,*) ' ' write(6,63) 'cj state (should be sonic with vel_mat = cs_cj):' write(6,61) 'temp_cj',temp_cj,'den_cj ',den_cj, 1 'pres_cj',pres_cj write(6,61) 'cs_cj ',cs_cj, 2 'vel_mat',vel_cj,'vel_det',vel_det write(6,61) 'mach_cj',vel_cj/cs_cj,'qburn_cj',qburn_cj 63 format(1x,a) 61 format(1x,a7,'=',1pe10.3,' ',a7,'=',1pe10.3,' ', 1 a7,'=',1pe10.3,' ',a4,'=',1pe10.3) write(6,*) ' ' write(6,*) 'top 10 cj nse mass fractions:' call indexx(ionmax,xmass_cj,izwork1) write(6,02) (ionam(izwork1(i)), 1 xmass_cj(izwork1(i)), i=ionmax,ionmax-9,-1) c..get shock solution kkase = 4 mach_sh = vel_det/cs_up call cjsolve(kkase,xmass_up,temp_up,den_up,mach_sh, 1 qdum,xmass_up,ener_up,pres_up,cs_up, 2 vel_det,vel_sh,temp_sh,den_sh,ener_sh,pres_sh,cs_sh) c..reset the initial conditions for a znd detonation tin = temp_sh din = den_sh vin = vel_sh zin = 1.0e-16*vel_sh den_stop = 1.00d0 * den_cj write(6,*) write(6,*) 'resetting initial conditions for a detonation to:' write(6,64) 'tin=',tin,' din=',din,' vin=',vin,' zin=',zin 64 format(1x,4(a,1pe12.4) ) end if c--------------------------------------------------------------------------- c..get the abundance variables for the final mixture call azbar(xin,aion,zion,ionmax, 1 zwork1,abar,zbar) c..get the thermodynamic state temp_row(1) = tin den_row(1) = din ptot_row(1) = bpres abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 if (ibtype .eq. 5) then den_row(1) = bpres * abar/(avo * kerg * tin) call invert_helm_pt din = den_row(1) else call helmeos bpres = ptot_row(1) endif ein = etot_row(1) c--------------------------------------------------------------------------- c..write out the final input write(6,*) write(6,02) 'tstart',tstart,'tstep',tstep write(6,02) 'tin',tin,'din',din,'bpres',bpres,'ein',ein c..largest mass fractions call indexx(ionmax,xin,izwork1) j = min(20,ionmax) k = max(ionmax-19,1) write(6,*) j,' largest mass fractions' do i=ionmax,k,-1 if (xin(izwork1(i)) .gt. 1.0e-12) 1 write(6,02) ionam(izwork1(i)),xin(izwork1(i)) end do c..nonconservation, abar, zbar of the mixture sum = 0.0d0 do i=1,ionmax sum = sum + xin(i) enddo write(6,02) '1-sum',1.0d0 - sum write(6,02) 'abar',abar,'zbar',zbar,'ye',zbar/abar write(6,*) c read(5,*) c..there is probably a better place for this c..if requested, adjust the number of equations being solved if (pure_network .eq. 1) then neqs = ionmax btemp = tin bden = din end if return end c include '../src_net/net_auxillary.f' c-------------------------------------------------------------------- c.. c..this routine contains auxillary network routine c..routines for a tree construction to mark nonzero matrix locations c..routine screen5 computes screening factors c..routine snupp computes neutrino loss rates for the pp chain c..routine snucno computes neutrino loss rates for the cno cycles c..routine sneut5 computes neutrino loss rates c..routine ifermi12 does an inverse fermi integral of order 1/2 c..routine zfermim12 does an inverse fermi integral of order -1/2 c..routine ecapnuc02 computes electron capture rates c..routine ecapnuc computes electron capture rates c..routine mazurek computes ni56 electron capture rates c..interfaces to a balanced tree sort subroutine tree_init(n) implicit none common/locatdat/nmax integer nmax,n nmax=n+1 call avlinit(30*nmax+2048) return end subroutine tree(i,j,eloc,neloc,nterm,nzo,iloc,jloc,np) implicit none common/locatdat/nmax integer nmax integer i,j,neloc,eloc(neloc),nterm,nzo,np,iloc(np),jloc(np) integer iat,nzo_old nzo_old = nzo nterm = nterm + 1 if (nterm .gt. neloc) then write(6,10) 'nterm=',nterm write(6,10) 'neloc=',neloc 10 format(1x,a,' ',i6) stop 'nterm > neloc in routine tree' end if call avlinsert(i*nmax+j,iat,nzo) eloc(nterm) = iat if (nzo .gt. np) stop 'nzo > np in routine tree' if (nzo .ne. nzo_old) then eloc(nterm) = nzo iloc(nzo) = i jloc(nzo) = j end if return end subroutine tree_out(irow,icol,nzo,np) implicit none common/locatdat/nmax integer nmax,np,nzo,i,irow(np),icol(np) call avlgetlist(np,icol,nzo) call avlfree() do i=1,nzo irow(i)=icol(i)/nmax icol(i)=icol(i)-irow(i)*nmax enddo return end c.... AVL sort c.... c.... In 1960 two Russian mathematicians, Georgii Maksimovich c.... Adel'son-Vel'skii and Evgenii Mikhailovich Landis developed a c.... technique for keeping a binary search tree balanced as items are c.... inserted into it. called AVL trees. c.... c.... efficiently sort integers in N log N operations c.... c.... implemetation taken from c.... http://www.moorpark.cc.ca.us/~csalazar/cs20/nonlin.txt (buggy) c.... see also c.... http://swww.ee.uwa.edu.au/~plsd210/ds/AVL.html c.... http://www.purists.org/georg/avltree/ (my favorite site) c.... c.... implemented by Alexander Heger, 20010129 c.... avldelete by Alexander Heger, 20010205 c======================================================================= c======================================================================= MODULE avl_data implicit none c integer, parameter :: maxavldata = 65536 integer :: maxavldata integer, parameter :: maxavlindex = 4 integer, parameter :: NULL = 0 integer, parameter :: l_LEFT = 1 integer, parameter :: l_RIGHT = l_LEFT+1 ! do not change integer, parameter :: l_BAL = 3 integer, parameter :: l_KEY = 4 integer, parameter :: i_ROOT = 1 integer, parameter :: i_NODEOFFSET = 1 integer, parameter :: l_ROOT = l_RIGHT integer, parameter :: l_RIGHTHEAVY = 1 ! do not change integer, parameter :: l_BALANCED = 0 ! do not change integer, parameter :: l_LEFTHEAVY = -l_RIGHTHEAVY integer, parameter :: l_UNBALANCED = l_BALANCED+1 integer, parameter :: l_GARBAGE = l_LEFT c.... tree data integer :: maxel integer :: garbage c integer, dimension(maxavlindex,maxavldata) :: ichild integer, allocatable, dimension(:,:) :: ichild SAVE END MODULE avl_data c======================================================================= MODULE avl_stack implicit none integer, parameter :: maxdepth = 48 integer, parameter :: i_STACKBASE = 1 integer, dimension(maxdepth) :: istack,lrstack integer :: ipstack END MODULE avl_stack c======================================================================= subroutine avlinit(nmax) USE avl_data implicit none integer, intent(IN) :: nmax SELECT CASE (nmax) CASE (1:) maxavldata=nmax+1 CASE (0) maxavldata=1024 END SELECT IF (nmax >= 0) THEN CALL avlfree() ALLOCATE(ichild(maxavlindex,maxavldata)) ENDIF c.... initialize root pointer and zero number of elements ichild(l_ROOT,i_ROOT)=NULL maxel=0 c.... initialize garbage list garbage=NULL end c======================================================================= subroutine avlfree() USE avl_data implicit none IF (ALLOCATED(ichild)) DEALLOCATE(ichild) end c======================================================================= subroutine avlgetlist(nmax,list,n) USE avl_data USE avl_stack implicit none c.... some constants integer, intent(IN) :: nmax integer, dimension(nmax), intent(OUT) :: list integer, intent(OUT) :: n c.... running variables integer :: i, lr, ii n=0 i=ichild(l_ROOT,i_ROOT) IF (i == NULL) RETURN IF (nmax < maxel) THEN WRITE(*,"(' [AVL LIST] ERROR: list too small for data.')") n=-1 RETURN ENDIF c.... recursively traverse tree to get sorted list of key values ipstack=i_STACKBASE-1 lr=l_LEFT DO IF (lr <= l_LEFT) THEN c.... add left branch ii=ichild(l_LEFT,i) IF (ii /= NULL) THEN ipstack=ipstack+1 istack(ipstack)=i lrstack(ipstack)=l_RIGHT i=ii lr=l_LEFT CYCLE ENDIF ENDIF IF (lr <= l_RIGHT) THEN c.... add node n=n+1 list(n)=ichild(l_KEY,i) c.... add right branch ii=ichild(l_RIGHT,i) IF (ii /= NULL) THEN ipstack=ipstack+1 istack(ipstack)=i lrstack(ipstack)=l_RIGHT+1 i=ii lr=l_LEFT CYCLE ENDIF ENDIF IF (ipstack < i_STACKBASE) EXIT i=istack(ipstack) lr=lrstack(ipstack) ipstack=ipstack-1 ENDDO IF (n /= maxel) THEN WRITE(*,"(' [AVL LIST] ERROR in AVL data.')") n=-1 RETURN ENDIF END c======================================================================= subroutine avltree() USE avl_data USE avl_stack implicit none c.... some constants character*(*), parameter :: form = "(I5)" integer, parameter :: nwidth = 5 integer, parameter :: nmax = 1024 c.... running variables integer, dimension(nmax) :: level, index character*(nwidth*nmax),dimension(5,maxdepth+1) :: line character*(nwidth) :: item integer :: i, lr, ii integer :: n, maxlevel integer :: l1,l2 IF (nmax < maxel) THEN WRITE(*,"(' [AVL TREE] ERROR: too much data.')") RETURN ENDIF n=0 maxlevel=0 i=ichild(l_ROOT,i_ROOT) IF (i == NULL) RETURN c.... recursively traverse tree to get sorted list of key values ipstack=i_STACKBASE-1 lr=l_LEFT DO IF (lr <= l_LEFT) THEN c.... add left branch ii=ichild(l_LEFT,i) IF (ii /= NULL) THEN ipstack=ipstack+1 istack(ipstack)=i lrstack(ipstack)=l_RIGHT i=ii lr=l_LEFT CYCLE ENDIF ENDIF IF (lr <= l_RIGHT) THEN c.... add node n=n+1 index(n)=i level(n)=ipstack+1 maxlevel=MAX(maxlevel,ipstack+1) c.... add right branch ii=ichild(l_RIGHT,i) IF (ii /= NULL) THEN ipstack=ipstack+1 istack(ipstack)=i lrstack(ipstack)=l_RIGHT+1 i=ii lr=l_LEFT CYCLE ENDIF ENDIF IF (ipstack < i_STACKBASE) EXIT i=istack(ipstack) lr=lrstack(ipstack) ipstack=ipstack-1 ENDDO line=" " DO i=1,n l1=1+nwidth*(i-1) l2=nwidth*i write(item,form) index(i) line(1,level(i))(l1:l2)=item write(item,form) ichild(l_KEY,index(i)) line(2,level(i))(l1:l2)=item write(item,form) ichild(l_BAL,index(i)) line(3,level(i))(l1:l2)=item write(item,form) ichild(l_LEFT,index(i)) line(4,level(i))(l1:l2)=item write(item,form) ichild(l_RIGHT,index(i)) line(5,level(i))(l1:l2)=item line(1,maxlevel+1)(l1:l2)='------------' ENDDO WRITE(*,"(A)") line(1,maxlevel+1)(1:nwidth*n)//"|" DO ii=1,maxlevel DO i=1,5 WRITE(*,"(A)") line(i,ii)(1:nwidth*n)//"|" ENDDO ENDDO WRITE(*,"(A)") line(1,maxlevel+1)(1:nwidth*n)//"|" END c======================================================================= subroutine avlincrease USE avl_data implicit none integer, allocatable, dimension(:,:) :: ichild_temp c integer, dimension(:,:) :: ichild_temp ALLOCATE(ichild_temp(maxavlindex,maxavldata)) ichild_temp(:,:)=ichild(:,:) DEALLOCATE(ichild) maxavldata=maxavldata*2 ALLOCATE(ichild(maxavlindex,maxavldata)) ichild(:,:)=ichild_temp(:,:) DEALLOCATE(ichild_temp) WRITE (*,"(' [AVL INCREASE] INFO: now ',I8,' elements.')") & maxavldata end c======================================================================= subroutine avlinsert(key,ijk,nzo) USE avl_data USE avl_stack implicit none integer, intent(IN) :: key integer, intent(OUT):: ijk,nzo integer :: i, ii, lr, irevbal, ic, ip, lrs, lri ijk = 1 i=i_root lr=l_ROOT ii=ichild(lr,i) ipstack=i_STACKBASE istack(ipstack)=i lrstack(ipstack)=lr c.... find location and insert DO WHILE (ii /= NULL) i=ii ijk = i c write(6,*) ijk SELECT CASE (key - ichild(l_KEY,i)) CASE (0) c write(6,*) 'same as ',ijk-1 ijk = ijk - 1 RETURN ! element already present: done. CASE (:-1) lr=l_LEFT CASE DEFAULT lr=l_RIGHT END SELECT ipstack=ipstack+1 istack(ipstack)=i lrstack(ipstack)=lr ii=ichild(lr,i) ENDDO c.... initialize new element maxel=maxel+1 nzo = maxel IF (garbage /= NULL) THEN ii=garbage garbage=ichild(l_GARBAGE,garbage) ELSE IF (maxel == maxavldata-1) CALL avlincrease ii=maxel+i_NODEOFFSET ENDIF ichild(lr,i)=ii ichild(l_KEY,ii)=key ichild(l_BAL,ii)=l_BALANCED ichild(l_LEFT:l_RIGHT,ii)=NULL c.... balance tree irevbal=l_UNBALANCED DO WHILE ((ipstack > i_STACKBASE) .AND. (irevbal /= l_BALANCED)) i=istack(ipstack) lr=lrstack(ipstack) ipstack=ipstack-1 lri=(l_LEFT+l_RIGHT)-lr ! pointer to the opposite direction ! sign for balance determination lrs=2*lr-(l_LEFT+l_RIGHT) SELECT CASE (ichild(l_BAL,i)*lrs) CASE (l_LEFTHEAVY) ichild(l_BAL,i)=l_BALANCED irevbal=l_BALANCED CASE (l_BALANCED) ichild(l_BAL,i)=lrs CASE DEFAULT c.... update tree ic=ichild(lr,i) IF (ichild(l_BAL,ic) == lrs) THEN c.... single rotation ichild(l_BAL,i)=l_BALANCED ichild(l_BAL,ic)=l_BALANCED ichild(lr,i)=ichild(lri,ic) ichild(lri,ic)=i ichild(lrstack(ipstack),istack(ipstack))=ic ELSE IF (ichild(l_BAL,ic) == -lrs) THEN c.... double rotation ip=ichild(lri,ic) SELECT CASE (ichild(l_BAL,ip)*lrs) CASE (l_LEFTHEAVY) ichild(l_BAL,i)=l_BALANCED ichild(l_BAL,ic)=lrs CASE (l_BALANCED) ichild(l_BAL,i)=l_BALANCED ichild(l_BAL,ic)=l_BALANCED CASE DEFAULT ichild(l_BAL,i)=-lrs ichild(l_BAL,ic)=l_BALANCED END SELECT ichild(l_BAL,ip)=l_BALANCED ichild(lri,ic)=ichild(lr,ip) ichild(lr,ip)=ic ichild(lr,i)=ichild(lri,ip) ichild(lri,ip)=i ichild(lrstack(ipstack),istack(ipstack))=ip ENDIF irevbal=l_BALANCED END SELECT ENDDO END c======================================================================= subroutine avldelete(key) USE avl_data USE avl_stack implicit none integer, intent(IN) :: key integer :: i, ii, lr, irevbal, ic, ip, lrs, lri, i0, ipstack0 lr=l_ROOT ipstack=i_STACKBASE istack(ipstack)=i_ROOT lrstack(ipstack)=l_ROOT i=ichild(l_ROOT,i_ROOT) c.... find location to delete DO IF (i == NULL) RETURN ! element not present SELECT CASE (key - ichild(l_KEY,i)) CASE (0) EXIT ! element found CASE (:-1) lr=l_LEFT CASE DEFAULT lr=l_RIGHT END SELECT ipstack=ipstack+1 istack(ipstack)=i lrstack(ipstack)=lr i=ichild(lr,i) ENDDO i0=i c.... find closest element to replace it c.... decide whether to take left or right branch SELECT CASE (ichild(l_BAL,i)) CASE (l_LEFTHEAVY) lri=l_LEFT CASE (l_RIGHTHEAVY) lri=l_RIGHT CASE DEFAULT lri=(l_LEFT+l_RIGHT)-lr END SELECT c.... now search for it ii=ichild(lri,i) IF (ii /= NULL) THEN ipstack0=ipstack c.... go one step in lrx direction ipstack=ipstack+1 istack(ipstack)=i lrstack(ipstack)=lri c.... now seach element most in opposite direction lr=(l_LEFT+l_RIGHT)-lri i=ii ii=ichild(lr,i) DO WHILE (ii /= NULL) ipstack=ipstack+1 istack(ipstack)=i lrstack(ipstack)=lr i=ii ii=ichild(lr,i) ENDDO c.... found element to swap c.... do swap ic=ichild(lri,i) ichild(lri,i)=ichild(lri,i0) ichild(lr,i)=ichild(lr,i0) ichild(l_BAL,i)=ichild(l_BAL,i0) ichild(lrstack(ipstack0),istack(ipstack0))=i c.... CORRECT STACK istack(ipstack0+1)=i ichild(lrstack(ipstack),istack(ipstack))=ic ELSE c.... element last of chain c.... just remove it ic=NULL ENDIF c.... move rest of branch one level up ichild(lrstack(ipstack),istack(ipstack))=ic c.... (i): balance=balance - lrs c.... start regular re-balancing loop irevbal=l_UNBALANCED DO WHILE ((ipstack > i_STACKBASE) .AND. (irevbal /= l_BALANCED)) i=istack(ipstack) lr=lrstack(ipstack) ipstack=ipstack-1 c.... lri=(l_LEFT+l_RIGHT)-lr lrs=2*lr-(l_LEFT+l_RIGHT) SELECT CASE (ichild(l_BAL,i)*lrs) CASE (l_RIGHTHEAVY) ichild(l_BAL,i)=l_BALANCED CASE (l_BALANCED) ichild(l_BAL,i)=-lrs irevbal=l_BALANCED CASE DEFAULT c.... update tree ic=ichild(lri,i) IF (ichild(l_BAL,ic) == lrs) THEN c.... double rotation ip=ichild(lr,ic) SELECT CASE (ichild(l_BAL,ip)*lrs) CASE (l_RIGHTHEAVY) ichild(l_BAL,i)=l_BALANCED ichild(l_BAL,ic)=-lrs CASE (l_LEFTHEAVY) ichild(l_BAL,i)=lrs ichild(l_BAL,ic)=l_BALANCED CASE DEFAULT ichild(l_BAL,i)=l_BALANCED ichild(l_BAL,ic)=l_BALANCED END SELECT ichild(l_BAL,ip)=l_BALANCED ichild(lri,i)=ichild(lr,ip) ichild(lr,ip)=i ichild(lr,ic)=ichild(lri,ip) ichild(lri,ip)=ic ichild(lrstack(ipstack),istack(ipstack))=ip ELSE c.... single rotation IF (ichild(l_BAL,ic) == l_BALANCED) THEN ichild(l_BAL,i)=-lrs ichild(l_BAL,ic)=lrs irevbal=l_BALANCED ELSE ichild(l_BAL,i)=l_BALANCED ichild(l_BAL,ic)=l_BALANCED ENDIF ichild(lri,i)=ichild(lr,ic) ichild(lr,ic)=i ichild(lrstack(ipstack),istack(ipstack))=ic ENDIF END SELECT ENDDO c.... free element ichild(l_GARBAGE,i0)=garbage garbage=i0 maxel=maxel-1 END subroutine screen5(temp,den,zbar,abar,z2bar, 1 z1,a1,z2,a2,jscreen,init, 2 scor,scordt,scordd) include 'implno.dek' c..this subroutine calculates screening factors and their derivatives c..for nuclear reaction rates in the weak, intermediate and strong regimes. c..based on graboske, dewit, grossman and cooper apj 181 457 1973 for c..weak screening. based on alastuey and jancovici apj 226 1034 1978, c..with plasma parameters from itoh et al apj 234 1079 1979, for strong c..screening. c..input: c..temp = temperature c..den = density c..zbar = mean charge per nucleus c..abar = mean number of nucleons per nucleus c..z2bar = mean square charge per nucleus c..z1 a1 = charge and number in the entrance channel c..z2 a2 = charge and number in the exit channel c..jscreen = counter of which reaction is being calculated c..init = flag to compute the more expensive functions just once c..output: c..scor = screening correction c..scordt = derivative of screening correction with temperature c..scordd = d