program public_iso7 include 'implno.dek' include 'const.dek' include 'timers.dek' include 'burn_common.dek' include 'network.dek' c..this program exercises the iso7 network c..declare integer i,j,nok,nbad double precision tstart,tstep,conserv,tin,din,ein,vin,zin,xin(7), 1 tout,dout,eout,xout(7) c..initialize the network call init_iso7 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 iso7,siso7,biso7,diso7 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 iso7 network call netint(beg,stptry,stpmin,tend,ys2, 1 tol,neqs,nok,nbad,kount,odescal, c 4 iso7,siso7,biso7,forder_ma28) c 4 iso7,siso7,biso7,forder_umf) c 4 iso7,siso7,biso7,forder_y12m) c 4 iso7,diso7,biso7,forder_ludcmp) c 4 iso7,diso7,biso7,forder_leqs) c 4 iso7,diso7,biso7,forder_lapack) c 4 iso7,diso7,biso7,forder_gift) c 4 iso7,siso7,biso7,forder_biconj) c 4 iso7,siso7,biso7,rosen_ma28) c 4 iso7,siso7,biso7,rosen_umf) c 4 iso7,siso7,biso7,rosen_y12m) c 4 iso7,diso7,biso7,rosen_ludcmp) c 4 iso7,diso7,biso7,rosen_leqs) c 4 iso7,diso7,biso7,rosen_lapack) c 4 iso7,diso7,biso7,rosen_gift) c 4 iso7,siso7,biso7,rosen_biconj) 4 iso7,siso7,biso7,stifbs_ma28) c 4 iso7,siso7,biso7,stifbs_umf) c 4 iso7,siso7,biso7,stifbs_y12m) c 4 iso7,diso7,biso7,stifbs_ludcmp) c 4 iso7,diso7,biso7,stifbs_leqs) c 4 iso7,diso7,biso7,stifbs_lapack) c 4 iso7,diso7,biso7,stifbs_gift) c 4 iso7,siso7,biso7,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 iso7 network c.. c..routine iso7 sets up the odes c..routine rhs evaluates the right hand sides c..routine diso7 sets up the dense iso7 jacobian c..routine biso7 build the nonzero locations for siso7 c..routine siso7 sets up the sparse iso7 jacobian c..routine iso7rat generates the reaction rates for routine iso7 c..routine iso7tab generates the raw rates by table interpolation c..routine screen_iso7 applies screening corrections to the raw rates c..routine init_iso7 initializes the iso7 network subroutine iso7(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 iso7 nuclear reactions. c.. c..isotopes: he4, c12, o16, ne20, mg25, si28, ni56 c..declare the pass double precision tt,y(1),dydt(1) c..local variables integer i double precision enuc,taud,taut,sum1,sum2,z, 1 zbarxx,ytot1,abar,zbar,ye,snuda,snudz double precision conv,sixth parameter (conv = ev2erg*1.0d6*avo, 1 sixth = 1.0d0/6.0d0) 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 iso7tab else call iso7rat end if c..do the screening here because the corrections depend on the composition call screen_iso7(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..append an energy equation dydt(iener) = enuc - sneut c..the type of temperature and density ode's depend c..on the burn mode: c..hydrostatic or single step cases if (hydrostatic .or. one_step .or. trho_hist) then dydt(itemp) = 0.0d0 dydt(iden) = 0.0d0 c..adiabatic expansion or contraction else if (expansion) then taud = 446.0d0/sqrt(den0) taut = 3.0d0 * taud dydt(itemp) = -psi * y(itemp)/taut dydt(iden) = -psi * y(iden)/taud c..self heating else if (self_heat_const_den) then c..call an eos temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos c..density equation dydt(iden) = 0.0d0 c..sum1 is d(ener)/d(yi)*d(yi)/dt = d(ener)/d(abar)*d(abar)/d(yi)*d(yi)/dt sum1 = 0.0d0 do i=1,ionmax sum1 = sum1 - dydt(i) enddo sum1 = sum1 * dea_row(1)*abar*abar c..sum2 is d(ener)/d(yi)*d(yi)/dt = d(ener)/d(zbar)*d(zbar)/d(yi)*d(yi)/dt sum2 = 0.0d0 do i=1,ionmax sum2 = sum2 + (zion(i) - zbar)*dydt(i) enddo sum2 = sum2 * dez_row(1)*abar c..temperature equation that is self-consistent with an eos z = 1.0d0/cv_row(1) dydt(itemp) = z*(dydt(iener) - ded_row(1)*dydt(iden) 1 - sum1 - sum2) end if return end subroutine rhs(y,rate,dydt) include 'implno.dek' include 'const.dek' include 'burn_common.dek' include 'network.dek' c..evaluates the right hand side of the achain odes c..declare the pass double precision y(1),rate(1),dydt(1) c..local variables integer i double precision sixth parameter (sixth = 1.0d0/6.0d0) c..zero the odes do i=1,ionmax dydt(i) = 0.0d0 enddo c..set up the system of ode's : c..4he reactions dydt(ihe4) = 3.0d0 * y(ic12) * rate(irg3a) 1 - 0.5d0 * y(ihe4) * y(ihe4) * y(ihe4) * rate(ir3a) 2 + y(io16) * rate(iroga) 3 - y(ic12) * y(ihe4) * rate(ircag) 4 + 0.5d0 * y(ic12) * y(ic12) * rate(ir1212) 5 + 0.5d0 * y(ic12) * y(io16) * rate(ir1216) 6 + 0.5d0 * y(io16) * y(io16) * rate(ir1616) 7 - y(io16) * y(ihe4) * rate(iroag) 8 + y(ine20) * rate(irnega) 9 + y(img24) * rate(irmgga) & - y(ine20) * y(ihe4) * rate(irneag) 1 + y(isi28) * rate(irsiga) 2 - y(img24) * y(ihe4) * rate(irmgag) c..12c reactions dydt(ic12) = sixth * y(ihe4) * y(ihe4) * y(ihe4) * rate(ir3a) 1 - y(ic12) * rate(irg3a) 2 + y(io16) * rate(iroga) 3 - y(ic12) * y(ihe4) * rate(ircag) 4 - y(ic12) * y(ic12) * rate(ir1212) 5 - y(ic12) * y(io16) * rate(ir1216) c..16o reactions dydt(io16) = -y(io16) * rate(iroga) 1 + y(ic12) * y(ihe4) * rate(ircag) 2 - y(ic12) * y(io16) * rate(ir1216) 3 - y(io16) * y(io16) * rate(ir1616) 4 - y(io16) * y(ihe4) * rate(iroag) 5 + y(ine20) * rate(irnega) c..20ne reactions dydt(ine20) = 0.5d0 * y(ic12) * y(ic12) * rate(ir1212) 1 + y(io16) * y(ihe4) * rate(iroag) 2 - y(ine20) * rate(irnega) 3 + y(img24) * rate(irmgga) 4 - y(ine20) * y(ihe4) * rate(irneag) c..24mg reactions dydt(img24) = 0.5d0 * y(ic12) * y(io16) * rate(ir1216) 1 - y(img24) * rate(irmgga) 2 + y(ine20) * y(ihe4) * rate(irneag) 3 + y(isi28) * rate(irsiga) 4 - y(img24) * y(ihe4) * rate(irmgag) c..28si reactions dydt(isi28) = 0.5d0 * y(ic12) * y(io16) * rate(ir1216) 1 + 0.5d0 * y(io16) * y(io16) * rate(ir1616) 2 - y(isi28) * rate(irsiga) 3 + y(img24) * y(ihe4) * rate(irmgag) 4 - rate(irsi2ni) * y(ihe4) 5 + rate(irni2si) * y(ini56) c..ni56 reactions dydt(ini56) = rate(irsi2ni) * y(ihe4) 1 - rate(irni2si) * y(ini56) return end subroutine diso7(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 iso7 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,xx,zz double precision conv,sixth parameter (conv = ev2erg*1.0d6*avo, 1 sixth = 1.0d0/6.0d0) c..zero the jacobian do j=1,nlog do i=1,nlog dfdy(i,j) = 0.0d0 enddo enddo c..positive definite mass fractions do i=1,ionmax y(i) = min(1.0d0,max(y(i),1.0d-30)) enddo c..generate abar and zbar for this composition zbarxx = 0.0d0 ytot1 = 0.0d0 do i=1,ionmax ytot1 = ytot1 + y(i) zbarxx = zbarxx + zion(i) * y(i) enddo abar = 1.0d0/ytot1 zbar = zbarxx * abar ye = zbar * ytot1 c..positive definite temperatures and densities y(itemp) = min(1.0d11,max(y(itemp),1.0d4)) y(iden) = min(1.0d11,max(y(iden),1.0d-10)) c..for evolution equations with the network if (pure_network .eq. 0) then c..get the new temperature and density if need be if (trho_hist) call update2(tt,y(itemp),y(iden)) if (self_heat_const_pres) then jlo_eos = 1 jhi_eos = 1 ptot_row(1) = bpres temp_row(1) = y(itemp) abar_row(1) = abar zbar_row(1) = zbar den_row(1) = y(iden) call invert_helm_pt_quiet y(iden) = den_row(1) end if c..positive definite temperatures and densities y(itemp) = min(1.0d11,max(y(itemp),1.0d4)) y(iden) = min(1.0d11,max(y(iden),1.0d-10)) c..set the common block temperature and density btemp = y(itemp) bden = y(iden) c..for pure network else if (pure_network .eq. 1) then if (trho_hist) call update2(tt,btemp,bden) end if c..get the reaction rates if (use_tables .eq. 1) then call iso7tab else call iso7rat end if c..do the screening here because the corrections depend on the composition call screen_iso7(y) c..set up the jacobian c..4he jacobian elementss dfdy(ihe4,ihe4) = -1.5d0 * y(ihe4) * y(ihe4) * ratdum(ir3a) 1 - y(ic12) * ratdum(ircag) 2 - y(io16) * ratdum(iroag) 3 - y(ine20) * ratdum(irneag) 4 - y(img24) * ratdum(irmgag) 5 - 7.0d0 * ratdum(irsi2ni) 6 - 7.0d0 * dratdumdy1(irsi2ni) * y(ihe4) 7 + 7.0d0 * dratdumdy1(irni2si) * y(ini56) dfdy(ihe4,ic12) = 3.0d0 * ratdum(irg3a) 1 - y(ihe4) * ratdum(ircag) 2 + y(ic12) * ratdum(ir1212) 3 + 0.5d0 * y(io16) * ratdum(ir1216) dfdy(ihe4,io16) = ratdum(iroga) 1 + 0.5d0 * y(ic12) * ratdum(ir1216) 2 + y(io16) * ratdum(ir1616) 3 - y(ihe4) * ratdum(iroag) dfdy(ihe4,ine20) = ratdum(irnega) 1 - y(ihe4) * ratdum(irneag) dfdy(ihe4,img24) = ratdum(irmgga) 1 - y(ihe4) * ratdum(irmgag) dfdy(ihe4,isi28) = ratdum(irsiga) 1 - 7.0d0 * dratdumdy2(irsi2ni) * y(ihe4) dfdy(ihe4,ini56) = 7.0d0 * ratdum(irni2si) c..12c jacobian elements dfdy(ic12,ihe4) = 0.5d0 * y(ihe4) * y(ihe4) * ratdum(ir3a) 1 - y(ic12) * ratdum(ircag) dfdy(ic12,ic12) = -ratdum(irg3a) 1 - y(ihe4) * ratdum(ircag) 2 - 2.0d0 * y(ic12) * ratdum(ir1212) 3 - y(io16) * ratdum(ir1216) dfdy(ic12,io16) = ratdum(iroga) 1 - y(ic12) * ratdum(ir1216) c..16o jacobian elements dfdy(io16,ihe4) = y(ic12) * ratdum(ircag) 1 - y(io16) * ratdum(iroag) dfdy(io16,ic12) = y(ihe4) * ratdum(ircag) 1 - y(io16) * ratdum(ir1216) dfdy(io16,io16) = -ratdum(iroga) 1 - y(ic12) * ratdum(ir1216) 2 - 2.0d0 * y(io16) * ratdum(ir1616) 3 - y(ihe4) * ratdum(iroag) dfdy(io16,ine20) = ratdum(irnega) c..20ne jacobian elements dfdy(ine20,ihe4) = y(io16) * ratdum(iroag) 1 - y(ine20) * ratdum(irneag) dfdy(ine20,ic12) = y(ic12) * ratdum(ir1212) dfdy(ine20,io16) = y(ihe4) * ratdum(iroag) dfdy(ine20,ine20) = -ratdum(irnega) 1 - y(ihe4) * ratdum(irneag) dfdy(ine20,img24) = ratdum(irmgga) c..24mg jacobian elements dfdy(img24,ihe4) = y(ine20) * ratdum(irneag) 1 - y(img24) * ratdum(irmgag) dfdy(img24,ic12) = 0.5d0 * y(io16) * ratdum(ir1216) dfdy(img24,io16) = 0.5d0 * y(ic12) * ratdum(ir1216) dfdy(img24,ine20) = y(ihe4) * ratdum(irneag) dfdy(img24,img24) = -ratdum(irmgga) 1 - y(ihe4) * ratdum(irmgag) dfdy(img24,isi28) = ratdum(irsiga) c..28si jacobian elements dfdy(isi28,ihe4) = y(img24) * ratdum(irmgag) 1 - ratdum(irsi2ni) 2 - dratdumdy1(irsi2ni) * y(ihe4) 3 + dratdumdy1(irni2si) * y(ini56) dfdy(isi28,ic12) = 0.5d0 * y(io16) * ratdum(ir1216) dfdy(isi28,io16) = y(io16) * ratdum(ir1616) 1 + 0.5d0 * y(ic12) * ratdum(ir1216) dfdy(isi28,img24) = y(ihe4) * ratdum(irmgag) dfdy(isi28,isi28) = -ratdum(irsiga) 1 - dratdumdy2(irsi2ni) * y(ihe4) dfdy(isi28,ini56) = ratdum(irni2si) c..ni56 jacobian elements dfdy(ini56,ihe4) = ratdum(irsi2ni) 1 + dratdumdy1(irsi2ni) * y(ihe4) 2 - dratdumdy1(irni2si) * y(ini56) dfdy(ini56,isi28) = dratdumdy2(irsi2ni) * y(ihe4) dfdy(ini56,ini56) = -ratdum(irni2si) 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 c..shut down the temperature and density derivatives c do i=1,ionmax c dfdy(i,itemp) = 0.0d0 c dfdy(i,iden) = 0.0d0 c enddo c dfdy(iener,itemp) = 0.0d0 c dfdy(iener,iden) = 0.0d0 return end subroutine biso7(iloc,jloc,nzo,np) include 'implno.dek' include 'network.dek' c.. c..this routine builds the nonzero matrix locations for siso7 c..input is the integer arrys iloc and jloc, both of dimension np, that c..on output contain nzo matrix element locations. c.. c..declare the pass integer np,iloc(np),jloc(np),nzo c..local variables integer i c..communicate with siso7 integer neloc parameter (neloc=77) integer eloc(neloc),nterms common /elcpp/ 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..4he jacobian elementss 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,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,ine20,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,img24,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,isi28,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ihe4,ini56,eloc,neloc,nterms,nzo,iloc,jloc,np) c..12c jacobian elements 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,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) c..16o jacobian elements call tree(io16,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io16,ic12,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io16,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(io16,ine20,eloc,neloc,nterms,nzo,iloc,jloc,np) c..20ne jacobian elements call tree(ine20,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine20,ic12,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine20,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine20,ine20,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ine20,img24,eloc,neloc,nterms,nzo,iloc,jloc,np) c..24mg jacobian elements call tree(img24,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(img24,ic12,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(img24,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(img24,ine20,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(img24,img24,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(img24,isi28,eloc,neloc,nterms,nzo,iloc,jloc,np) c..28si jacobian elements call tree(isi28,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(isi28,ic12,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(isi28,io16,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(isi28,img24,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(isi28,isi28,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(isi28,ini56,eloc,neloc,nterms,nzo,iloc,jloc,np) c..ni56 jacobian elements call tree(ini56,ihe4,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ini56,isi28,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(ini56,ini56,eloc,neloc,nterms,nzo,iloc,jloc,np) c..if we are doing a pure network, we are done if (pure_network .eq. 1) return c..temperature contributions do i=1,ionmax call tree(i,itemp,eloc,neloc,nterms,nzo,iloc,jloc,np) end do c..density contributions do i=1,ionmax call tree(i,iden,eloc,neloc,nterms,nzo,iloc,jloc,np) end do c..energy equation jacobian elements do i=1,ionmax call tree(iener,i,eloc,neloc,nterms,nzo,iloc,jloc,np) enddo call tree(iener,iener,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(iener,itemp,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(iener,iden,eloc,neloc,nterms,nzo,iloc,jloc,np) c..neutrino losses do i=1,ionmax call tree(iener,i,eloc,neloc,nterms,nzo,iloc,jloc,np) enddo call tree(iener,itemp,eloc,neloc,nterms,nzo,iloc,jloc,np) call tree(iener,iden,eloc,neloc,nterms,nzo,iloc,jloc,np) c..the jacobian elements 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 siso7(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 iso7 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 zbarxx,ytot1,abar,zbar,ye,taud,taut, 1 snuda,snudz,a1,a2,a3,a4,zz double precision conv,sixth parameter (conv = ev2erg*1.0d6*avo, 1 sixth = 1.0d0/6.0d0) c..communicate with the jacobian builder integer neloc parameter (neloc=77) integer eloc(neloc),nterms common /elcpp/ 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 iso7tab else call iso7rat end if c..do the screening here because the corrections depend on the composition call screen_iso7(y) c..set up the jacobian c..4he jacobian elements c..d(he4)/d(he4) a1 = -1.5d0 * y(ihe4) * y(ihe4) * ratdum(ir3a) 1 - y(ic12) * ratdum(ircag) 2 - y(io16) * ratdum(iroag) 3 - y(ine20) * ratdum(irneag) 4 - y(img24) * ratdum(irmgag) 5 - 7.0d0 * ratdum(irsi2ni) 6 - 7.0d0 * dratdumdy1(irsi2ni) * y(ihe4) 7 + 7.0d0 * dratdumdy1(irni2si) * y(ini56) 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 = 3.0d0 * ratdum(irg3a) 1 - y(ihe4) * ratdum(ircag) 2 + y(ic12) * ratdum(ir1212) 3 + 0.5d0 * y(io16) * ratdum(ir1216) 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(o16) a1 = ratdum(iroga) 1 + 0.5d0 * y(ic12) * ratdum(ir1216) 2 + y(io16) * ratdum(ir1616) 3 - y(ihe4) * ratdum(iroag) 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(ne20) a1 = ratdum(irnega) 1 - y(ihe4) * ratdum(irneag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine20) = xsum(ine20) + a1 * bion(ihe4) ysum(ine20) = ysum(ine20) - a1 zsum(ine20) = zsum(ine20) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(mg24) a1 = ratdum(irmgga) 1 - y(ihe4) * ratdum(irmgag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(img24) = xsum(img24) + a1 * bion(ihe4) ysum(img24) = ysum(img24) - a1 zsum(img24) = zsum(img24) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(si28) a1 = ratdum(irsiga) 1 - 7.0d0 * dratdumdy2(irsi2ni) * y(ihe4) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(isi28) = xsum(isi28) + a1 * bion(ihe4) ysum(isi28) = ysum(isi28) - a1 zsum(isi28) = zsum(isi28) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(ni56) a1 = 7.0d0 * ratdum(irni2si) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ini56) = xsum(ini56) + a1 * bion(ihe4) ysum(ini56) = ysum(ini56) - a1 zsum(ini56) = zsum(ini56) + a1 * (zion(ihe4) - zbar) c..12c jacobian elements c..d(c12)/d(he4) a1 = 0.5d0 * y(ihe4) * y(ihe4) * ratdum(ir3a) 1 - y(ic12) * ratdum(ircag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ic12) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ic12) - zbar) c..d(c12)/d(c12) a1 = -ratdum(irg3a) 1 - y(ihe4) * ratdum(ircag) 2 - 2.0d0 * y(ic12) * ratdum(ir1212) 3 - y(io16) * ratdum(ir1216) 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(o16) a1 = ratdum(iroga) 1 - y(ic12) * ratdum(ir1216) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(ic12) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(ic12) - zbar) c..16o jacobian elements c..d(o16)/d(he4) a1 = y(ic12) * ratdum(ircag) 1 - y(io16) * ratdum(iroag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(io16) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(io16) - zbar) c..d(o16)/d(c12) a1 = y(ihe4) * ratdum(ircag) 1 - y(io16) * ratdum(ir1216) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(io16) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(io16) - zbar) c..d(o16)/d(o16) a1 = -ratdum(iroga) 1 - y(ic12) * ratdum(ir1216) 2 - 2.0d0 * y(io16) * ratdum(ir1616) 3 - y(ihe4) * ratdum(iroag) 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(ne20) a1 = ratdum(irnega) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine20) = xsum(ine20) + a1 * bion(io16) ysum(ine20) = ysum(ine20) - a1 zsum(ine20) = zsum(ine20) + a1 * (zion(io16) - zbar) c..20ne jacobian elements c..d(ne20)/d(he4) a1 = y(io16) * ratdum(iroag) - y(ine20) * ratdum(irneag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ine20) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ine20) - zbar) c..d(ne20)/d(c12) a1 = y(ic12) * ratdum(ir1212) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(ine20) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(ine20) - zbar) c..d(ne20)/d(o16) a1 = y(ihe4) * ratdum(iroag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(ine20) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(ine20) - zbar) c..d(ne20)/d(ne20) a1 = -ratdum(irnega) - y(ihe4) * ratdum(irneag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine20) = xsum(ine20) + a1 * bion(ine20) ysum(ine20) = ysum(ine20) - a1 zsum(ine20) = zsum(ine20) + a1 * (zion(ine20) - zbar) c..d(ne20)/d(mg24) a1 = ratdum(irmgga) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(img24) = xsum(img24) + a1 * bion(ine20) ysum(img24) = ysum(img24) - a1 zsum(img24) = zsum(img24) + a1 * (zion(ine20) - zbar) c..24mg jacobian elements c..d(mg24)/d(he4) a1 = y(ine20) * ratdum(irneag) 1 - y(img24) * ratdum(irmgag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(img24) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(img24) - zbar) c..d(mg24)/d(c12) a1 = 0.5d0 * y(io16) * ratdum(ir1216) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(img24) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(img24) - zbar) c..d(mg24)/d(o16) a1 = 0.5d0 * y(ic12) * ratdum(ir1216) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(img24) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(img24) - zbar) c..d(mg24)/d(ne20) a1 = y(ihe4) * ratdum(irneag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ine20) = xsum(ine20) + a1 * bion(img24) ysum(ine20) = ysum(ine20) - a1 zsum(ine20) = zsum(ine20) + a1 * (zion(img24) - zbar) c..d(mg24)/d(mg24) a1 = -ratdum(irmgga) 1 - y(ihe4) * ratdum(irmgag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(img24) = xsum(img24) + a1 * bion(img24) ysum(img24) = ysum(img24) - a1 zsum(img24) = zsum(img24) + a1 * (zion(img24) - zbar) c..d(mg24)/d(si28) a1 = ratdum(irsiga) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(isi28) = xsum(isi28) + a1 * bion(img24) ysum(isi28) = ysum(isi28) - a1 zsum(isi28) = zsum(isi28) + a1 * (zion(img24) - zbar) c..28si jacobian elements c..d(si28)/d(he4) a1 = y(img24) * ratdum(irmgag) 1 - ratdum(irsi2ni) 2 - dratdumdy1(irsi2ni) * y(ihe4) 3 + dratdumdy1(irni2si) * y(ini56) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(isi28) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(isi28) - zbar) c..d(si28)/d(c12) a1 = 0.5d0 * y(io16) * ratdum(ir1216) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(isi28) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(isi28) - zbar) c..d(si28)/d(o16) a1 = y(io16) * ratdum(ir1616) 1 + 0.5d0 * y(ic12) * ratdum(ir1216) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(isi28) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(isi28) - zbar) c..d(si28)/d(mg24) a1 = y(ihe4) * ratdum(irmgag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(img24) = xsum(img24) + a1 * bion(isi28) ysum(img24) = ysum(img24) - a1 zsum(img24) = zsum(img24) + a1 * (zion(isi28) - zbar) c..d(si28)/d(si28) a1 = -ratdum(irsiga) 1 - dratdumdy2(irsi2ni) * y(ihe4) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(isi28) = xsum(isi28) + a1 * bion(isi28) ysum(isi28) = ysum(isi28) - a1 zsum(isi28) = zsum(isi28) + a1 * (zion(isi28) - zbar) c..d(si28)/d(ni56) a1 = ratdum(irni2si) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ini56) = xsum(ini56) + a1 * bion(isi28) ysum(ini56) = ysum(ini56) - a1 zsum(ini56) = zsum(ini56) + a1 * (zion(isi28) - zbar) c..ni56 jacobian elements c..d(ni56)/d(he4) a1 = ratdum(irsi2ni) 1 + dratdumdy1(irsi2ni) * y(ihe4) 2 - dratdumdy1(irni2si) * y(ini56) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ini56) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ini56) - zbar) c..d(ni56)/d(si28) a1 = dratdumdy2(irsi2ni) * y(ihe4) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(isi28) = xsum(isi28) + a1 * bion(ini56) ysum(isi28) = ysum(isi28) - a1 zsum(isi28) = zsum(isi28) + a1 * (zion(ini56) - zbar) c..d(ni56)/d(ni56) a1 = -ratdum(irni2si) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ini56) = xsum(ini56) + a1 * bion(ini56) ysum(ini56) = ysum(ini56) - a1 zsum(ini56) = zsum(ini56) + a1 * (zion(ini56) - zbar) c..if we are doing a pure network, we are done, head to the error check 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 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 siso7: nt .ne. nterms' stop 'error in routine siso7' end if return end subroutine iso7rat include 'implno.dek' include 'burn_common.dek' include 'network.dek' c..this routine generates nuclear reaction rates for the iso7 network. 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..c12(a,g)o16 call rate_c12ag(btemp,bden, 1 ratraw(ircag),dratrawdt(ircag),dratrawdd(ircag), 2 ratraw(iroga),dratrawdt(iroga),dratrawdd(iroga)) c..triple alpha to c12 call rate_tripalf(btemp,bden, 1 ratraw(ir3a),dratrawdt(ir3a),dratrawdd(ir3a), 2 ratraw(irg3a),dratrawdt(irg3a),dratrawdd(irg3a)) c..c12 + c12 call rate_c12c12(btemp,bden, 1 ratraw(ir1212),dratrawdt(ir1212),dratrawdd(ir1212), 2 rrate,drratedt,drratedd) c..c12 + o16 call rate_c12o16(btemp,bden, 1 ratraw(ir1216),dratrawdt(ir1216),dratrawdd(ir1216), 2 rrate,drratedt,drratedd) c..16o + 16o call rate_o16o16(btemp,bden, 1 ratraw(ir1616),dratrawdt(ir1616),dratrawdd(ir1616), 2 rrate,drratedt,drratedd) c..o16(a,g)ne20 call rate_o16ag(btemp,bden, 1 ratraw(iroag),dratrawdt(iroag),dratrawdd(iroag), 2 ratraw(irnega),dratrawdt(irnega),dratrawdd(irnega)) c..ne20(a,g)mg24 call rate_ne20ag(btemp,bden, 1 ratraw(irneag),dratrawdt(irneag),dratrawdd(irneag), 2 ratraw(irmgga),dratrawdt(irmgga),dratrawdd(irmgga)) c..mg24(a,g)si28 call rate_mg24ag(btemp,bden, 1 ratraw(irmgag),dratrawdt(irmgag),dratrawdd(irmgag), 2 ratraw(irsiga),dratrawdt(irsiga),dratrawdd(irsiga)) c..ca40(a,g)ti44 call rate_ca40ag(btemp,bden, 1 ratraw(ircaag),dratrawdt(ircaag),dratrawdd(ircaag), 2 ratraw(irtiga),dratrawdt(irtiga),dratrawdd(irtiga)) c..temp depend of partition functions c gca40=1.0d0 +exp((-4.150E+01 + 1.636E+00*t9 + 1.483E-01*t92)*t9ri) c gti44=1.0d0 +exp((-1.111E+01 + 6.293E-01*t9 + 1.732E-01*t92)*t9ri) c ratraw(irtiga) = ratraw(irtiga) * gca40/gti44 c ratraw(ircag) = 0.0d0 c ratraw(iroga) = 0.0d0 c ratraw(ir3a) = 0.0d0 c ratraw(irg3a) = 0.0d0 c ratraw(ir1212) = 0.0d0 c ratraw(ir1216) = 0.0d0 c ratraw(ir1616) = 0.0d0 c ratraw(iroag) = 0.0d0 c ratraw(irnega) = 0.0d0 c ratraw(irneag) = 0.0d0 c ratraw(irmgga) = 0.0d0 c ratraw(irmgag) = 0.0d0 c ratraw(irsiga) = 0.0 c ratraw(ircaag) = 1.0d-40 c ratraw(irtiga) = 1.0d-40 return end subroutine iso7tab include 'implno.dek' include 'burn_common.dek' include 'network.dek' c..uses tables instead of analytical expressions to evaluate the c..raw reaction rates. a cubic polynomial is hardwired for speed. integer i,j,imax,iat,mp,ifirst parameter (mp = 4) double precision tlo,thi,tstp,bden_sav,btemp_sav, 1 x,x1,x2,x3,x4,a,b,c,d,e,f,g,h,p,q, 2 alfa,beta,gama,delt data ifirst/0/ c..make the table if (ifirst .eq. 0) then ifirst = 1 c..set the log temperature loop limits thi = 10.0d0 tlo = 6.0d0 imax = int(thi-tlo)*120 + 1 if (imax .gt. nrattab) stop 'imax too small in aprox13tab' tstp = (thi - tlo)/float(imax-1) c..save the input btemp_sav = btemp bden_sav = bden c..form the table bden = 1.0d0 do i=1,imax btemp = tlo + float(i-1)*tstp btemp = 10.0d0**(btemp) call iso7rat 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(ircag) = bden dtab(iroga) = 1.0d0 dtab(ir3a) = bden*bden dtab(irg3a) = 1.0d0 dtab(ir1212) = bden dtab(ir1216) = bden dtab(ir1616) = bden dtab(iroag) = bden dtab(irnega) = 1.0d0 dtab(irneag) = bden dtab(irmgga) = 1.0d0 dtab(irmgag) = bden dtab(irsiga) = 1.0d0 dtab(ircaag) = bden dtab(irtiga) = 1.0d0 c..hash locate the temperature iat = int((log10(btemp) - tlo)/tstp) + 1 iat = max(1,min(iat - mp/2 + 1,imax - mp + 1)) c..setup the lagrange interpolation coefficients for a cubic x = btemp x1 = ttab(iat) x2 = ttab(iat+1) x3 = ttab(iat+2) x4 = ttab(iat+3) a = x - x1 b = x - x2 c = x - x3 d = x - x4 e = x1 - x2 f = x1 - x3 g = x1 - x4 h = x2 - x3 p = x2 - x4 q = x3 - x4 alfa = b*c*d/(e*f*g) beta = -a*c*d/(e*h*p) gama = a*b*d/(f*h*q) delt = -a*b*c/(g*p*q) c..crank off the raw reaction rates do j=1,nrat ratraw(j) = (alfa*rattab(j,iat) 1 + beta*rattab(j,iat+1) 2 + gama*rattab(j,iat+2) 3 + delt*rattab(j,iat+3) 4 ) * dtab(j) dratrawdt(j) = (alfa*drattabdt(j,iat) 1 + beta*drattabdt(j,iat+1) 2 + gama*drattabdt(j,iat+2) 3 + delt*drattabdt(j,iat+3) 4 ) * dtab(j) dratrawdd(j) = alfa*drattabdd(j,iat) 1 + beta*drattabdd(j,iat+1) 2 + gama*drattabdd(j,iat+2) 3 + delt*drattabdd(j,iat+3) enddo c..hand finish the three body reactions dratrawdd(ir3a) = bden * dratrawdd(ir3a) return end subroutine screen_iso7(y) include 'implno.dek' include 'burn_common.dek' include 'tfactors.dek' include 'network.dek' c..this routine computes the screening factors c..and applies them to the raw reaction rates, c..producing the final reaction rates used by the c..right hand sides and jacobian matrix elements c..this routine assumes screen_on = 1 or = 0 has been set at a higer level, c..presumably in the top level driver c..declare integer i,jscr double precision y(*),sc1a,sc1adt,sc1add,sc2a,sc2adt,sc2add, 1 sc3a,sc3adt,sc3add, 2 abar,zbar,z2bar,ytot1,zbarxx,z2barxx, 3 t992,t9i92,denom,zz, 4 yeff_ca40,yeff_ca40dt,yeff_ti44,yeff_ti44dt integer init data init/1/ c..roll all of them 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..form those lovely dummy links c..rsi2ni is the rate for silicon to nickel c..rni2si is the rate for nickel to silicon ratdum(irsi2ni) = 0.0d0 dratdumdy1(irsi2ni) = 0.0d0 dratdumdy2(irsi2ni) = 0.0d0 dratdumdt(irsi2ni) = 0.0d0 dratdumdd(irsi2ni) = 0.0d0 ratdum(irni2si) = 0.0d0 dratdumdy1(irni2si) = 0.0d0 dratdumdt(irni2si) = 0.0d0 dratdumdd(irni2si) = 0.0d0 c..if screening is on if (screen_on .ne. 0) then c..with the passed composition, compute abar,zbar and other variables zbarxx = 0.0d0 z2barxx = 0.0d0 ytot1 = 0.0d0 do i=1,ionmax ytot1 = ytot1 + y(i) zbarxx = zbarxx + zion(i) * y(i) z2barxx = z2barxx + zion(i) * zion(i) * y(i) enddo abar = 1.0d0/ytot1 zbar = zbarxx * abar z2bar = z2barxx * abar c..first the always fun triple alpha and its inverse jscr = 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ihe4),aion(ihe4),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ihe4),aion(ihe4),4.0d0,8.0d0, 2 jscr,init,sc2a,sc2adt,sc2add) sc3a = sc1a * sc2a sc3adt = sc1adt*sc2a + sc1a*sc2adt sc3add = sc1add*sc2a + sc1a*sc2add ratdum(ir3a) = ratraw(ir3a) * sc3a dratdumdt(ir3a) = dratrawdt(ir3a)*sc3a + ratraw(ir3a)*sc3adt dratdumdd(ir3a) = dratrawdd(ir3a)*sc3a + ratraw(ir3a)*sc3add scfac(ir3a) = sc3a dscfacdt(ir3a) = sc3adt dscfacdd(ir3a) = sc3add c..c12 to o16 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(ircag) = ratraw(ircag) * sc1a dratdumdt(ircag) = dratrawdt(ircag)*sc1a + ratraw(ircag)*sc1adt dratdumdd(ircag) = dratrawdd(ircag)*sc1a + ratraw(ircag)*sc1add scfac(ircag) = sc1a dscfacdt(ircag) = sc1adt dscfacdt(ircag) = sc1add c..c12 + c12 jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ic12),aion(ic12),zion(ic12),aion(ic12), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(ir1212) = ratraw(ir1212) * sc1a dratdumdt(ir1212) = dratrawdt(ir1212)*sc1a + ratraw(ir1212)*sc1adt dratdumdd(ir1212) = dratrawdd(ir1212)*sc1a + ratraw(ir1212)*sc1add scfac(ir1212) = sc1a dscfacdt(ir1212) = sc1adt dscfacdd(ir1212) = sc1add c..c12 + o16 jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ic12),aion(ic12),zion(io16),aion(io16), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(ir1216) = ratraw(ir1216) * sc1a dratdumdt(ir1216) = dratrawdt(ir1216)*sc1a + ratraw(ir1216)*sc1adt dratdumdd(ir1216) = dratrawdd(ir1216)*sc1a + ratraw(ir1216)*sc1add scfac(ir1216) = sc1a dscfacdt(ir1216) = sc1adt dscfacdd(ir1216) = sc1add c..o16 + o16 jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ic12),aion(ic12),zion(io16),aion(io16), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(ir1216) = ratraw(ir1216) * sc1a dratdumdt(ir1216) = dratrawdt(ir1216)*sc1a + ratraw(ir1216)*sc1adt dratdumdd(ir1216) = dratrawdd(ir1216)*sc1a + ratraw(ir1216)*sc1add scfac(ir1216) = sc1a dscfacdt(ir1216) = sc1adt dscfacdd(ir1216) = sc1add c..o16 to ne20 jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(io16),aion(io16),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(iroag) = ratraw(iroag) * sc1a dratdumdt(iroag) = dratrawdt(iroag)*sc1a + ratraw(iroag)*sc1adt dratdumdd(iroag) = dratrawdd(iroag)*sc1a + ratraw(iroag)*sc1add scfac(iroag) = sc1a dscfacdt(iroag) = sc1adt dscfacdd(iroag) = sc1add c..o16 to mg24 jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ine20),aion(ine20),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irneag) = ratraw(irneag) * sc1a dratdumdt(irneag) = dratrawdt(irneag)*sc1a + ratraw(irneag)*sc1adt dratdumdd(irneag) = dratrawdd(irneag)*sc1a + ratraw(irneag)*sc1add scfac(irneag) = sc1a dscfacdt(irneag) = sc1adt dscfacdd(irneag) = sc1add c..mg24 to si28 jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(img24),aion(img24),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(irmgag) = ratraw(irmgag) * sc1a dratdumdt(irmgag) = dratrawdt(irmgag)*sc1a + ratraw(irmgag)*sc1adt dratdumdd(irmgag) = dratrawdd(irmgag)*sc1a + ratraw(irmgag)*sc1add scfac(irmgag) = sc1a dscfacdt(irmgag) = sc1adt dscfacdd(irmgag) = sc1add c..ca40 to ti44 jscr = jscr + 1 call screen5(btemp,bden,zbar,abar,z2bar, 1 zion(ica40),aion(ica40),zion(ihe4),aion(ihe4), 2 jscr,init,sc1a,sc1adt,sc1add) ratdum(ircaag) = ratraw(ircaag) * sc1a dratdumdt(ircaag) = dratrawdt(ircaag)*sc1a + ratraw(ircaag)*sc1adt dratdumdd(ircaag) = dratrawdd(ircaag)*sc1a + ratraw(ircaag)*sc1add scfac(ircaag) = sc1a dscfacdt(ircaag) = sc1adt dscfacdd(ircaag) = sc1add c..reset the screen initialization flag init = 0 c..end of the screening if block end if if (t9 .gt. 2.5 .and. y(ic12)+y(io16) .le. 4.0d-3) then t992 = t972 * t9 t9i92 = 1.0d0/t992 yeff_ca40 = t9i92 * exp(239.42*t9i - 74.741) yeff_ca40dt = -yeff_ca40*(239.42*t9i2 + 4.5d0*t9i) yeff_ti44 = t992 * exp(-274.12*t9i + 74.914) yeff_ti44dt = yeff_ti44*(274.12*t9i2 + 4.5d0*t9i) denom = (bden * y(ihe4))**3 ratdum(irsi2ni) = yeff_ca40*denom*ratdum(ircaag)*y(isi28) dratdumdy1(irsi2ni) = 3.0d0 * ratdum(irsi2ni)/y(ihe4) dratdumdy2(irsi2ni) = yeff_ca40*denom*ratdum(ircaag) dratdumdt(irsi2ni) = (yeff_ca40dt*ratdum(ircaag) 1 + yeff_ca40*dratdumdt(ircaag))*denom*y(isi28)*1.0d-9 dratdumdd(irsi2ni) = 3.0d0*ratdum(irsi2ni)/bden 1 + yeff_ca40*denom*dratdumdd(ircaag)*y(isi28) if (denom .ne. 0.0) then zz = 1.0d0/denom ratdum(irni2si) = min(1.0d10,yeff_ti44*ratdum(irtiga)*zz) if (ratdum(irni2si) .eq. 1.0d10) then dratdumdy1(irni2si) = 0.0d0 dratdumdt(irni2si) = 0.0d0 dratdumdd(irni2si) = 0.0d0 else dratdumdy1(irni2si) = -3.0d0 * ratdum(irni2si)/y(ihe4) dratdumdt(irni2si) = (yeff_ti44dt*ratdum(irtiga) 1 + yeff_ti44*dratdumdt(irtiga))*zz*1.0d-9 dratdumdd(irni2si) = -3.0d0 * ratdum(irni2si)/bden 1 + yeff_ti44*dratdumdd(irtiga)*zz end if endif end if return end subroutine init_iso7 include 'implno.dek' include 'network.dek' c..this routine initializes stuff for the iso7 network 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 = idiso7 ionmax = 7 iener = ionmax + 1 itemp = ionmax + 2 iden = ionmax + 3 neqs = iden nrat = 17 netname = 'iso7' c..set the id numbers of the elements ihe4 = 1 ic12 = 2 io16 = 3 ine20 = 4 img24 = 5 isi28 = 6 ini56 = 7 c..set the names of the elements ionam(ihe4) = 'he4 ' ionam(ic12) = 'c12 ' ionam(io16) = 'o16 ' ionam(ine20) = 'ne20' ionam(img24) = 'mg24' ionam(isi28) = 'si28' ionam(ini56) = 'ni56' ionam(iener) = 'ener ' ionam(itemp) = 'temp ' ionam(iden) = 'den ' c..set the number of nucleons in the element aion(ihe4) = 4.0d0 aion(ic12) = 12.0d0 aion(io16) = 16.0d0 aion(ine20) = 20.0d0 aion(img24) = 24.0d0 aion(isi28) = 28.0d0 aion(ini56) = 56.0d0 c..set the number of protons in the element zion(ihe4) = 2.0d0 zion(ic12) = 6.0d0 zion(io16) = 8.0d0 zion(ine20) = 10.0d0 zion(img24) = 12.0d0 zion(isi28) = 14.0d0 zion(ini56) = 28.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(ihe4) = 28.29603d0 bion(ic12) = 92.16294d0 bion(io16) = 127.62093d0 bion(ine20) = 160.64788d0 bion(img24) = 198.25790d0 bion(isi28) = 236.53790d0 bion(ini56) = 484.00300d0 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 ircag = 1 iroga = 2 ir3a = 3 irg3a = 4 ir1212 = 5 ir1216 = 6 ir1616 = 7 iroag = 8 irnega = 9 irneag = 10 irmgga = 11 irmgag = 12 irsiga = 13 ircaag = 14 irtiga = 15 irsi2ni = 16 irni2si = 17 c..set the names of the reaction rates ratnam(ircag) = 'rcag ' ratnam(iroga) = 'roga ' ratnam(ir3a) = 'r3a ' ratnam(irg3a) = 'rg3a ' ratnam(ir1212) = 'r1212' ratnam(ir1216) = 'r1216' ratnam(ir1616) = 'r1616' ratnam(iroag) = 'roag ' ratnam(irnega) = 'rnega' ratnam(irneag) = 'rneag' ratnam(irmgga) = 'rmgga' ratnam(irmgag) = 'rmgag' ratnam(irsiga) = 'rsiga' ratnam(ircaag) = 'rcaag' ratnam(irtiga) = 'rtiga' ratnam(irsi2ni) = 'rsi2ni' ratnam(irni2si) = 'rni2si' 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 = derivative of screening correction with density c..declare the pass integer jscreen,init double precision temp,den,zbar,abar,z2bar,z1,a1,z2,a2, 1 scor,scordt,scordd c..local variables double precision aa,daadt,daadd,bb,cc,dccdt,dccdd, 1 pp,dppdt,dppdd,qq,dqqdt,dqqdd,rr,drrdt,drrdd, 2 ss,dssdt,dssdd,tt,dttdt,dttdd,uu,duudt,duudd, 3 vv,dvvdt,dvvdd,a3,da3,tempi,dtempi,deni, 2 qlam0z,qlam0zdt,qlam0zdd, 3 h12w,dh12wdt,dh12wdd,h12,dh12dt,dh12dd, 4 taufac,taufacdt,gamp,gampdt,gampdd, 5 gamef,gamefdt,gamefdd, 6 tau12,tau12dt,alph12,alph12dt,alph12dd, 7 xlgfac,dxlgfacdt,dxlgfacdd, 8 gamp14,gamp14dt,gamp14dd, 9 xni,dxnidd,ytot, & temp_old,den_old,zbar_old,abar_old c..screening variables c..zs13 = (z1+z2)**(1./3.) c..zhat = combination of z1 and z2 raised to the 5/3 power c..zhat2 = combination of z1 and z2 raised to the 5/12 power c..lzav = log of effective charge c..aznut = combination of a1,z1,a2,z2 raised to 1/3 power integer abigrat parameter (abigrat = 8000) double precision zs13(abigrat),zhat(abigrat), 1 zhat2(abigrat),lzav(abigrat), 2 aznut(abigrat),zs13inv(abigrat) c..parameter fact is the cube root of 2 double precision x13,x14,x53,x532,x512,fact,co2 parameter (x13 = 1.0d0/3.0d0, 1 x14 = 1.0d0/4.0d0, 3 x53 = 5.0d0/3.0d0, 4 x532 = 5.0d0/32.0d0, 5 x512 = 5.0d0/12.0d0, 6 fact = 1.25992104989487d0, 7 co2 = x13 * 4.248719d3) data temp_old/-1.0d0/, den_old/-1.0d0/, 1 zbar_old/-1.0d0/, abar_old/-1.0d0/ c..compute and store the more expensive screening factors if (init .eq. 1) then if (jscreen .gt. abigrat) stop 'jscreen > abigrat in screen5' zs13(jscreen) = (z1 + z2)**x13 zs13inv(jscreen) = 1.0d0/zs13(jscreen) zhat(jscreen) = (z1 + z2)**x53 - z1**x53 - z2**x53 zhat2(jscreen) = (z1 + z2)**x512 - z1**x512 -z2**x512 lzav(jscreen) = x53 * log(z1*z2/(z1 + z2)) aznut(jscreen) = (z1**2 * z2**2 * a1*a2 / (a1 + a2))**x13 endif c..calculate average plasma, if need be if (temp_old .ne. temp .or. 1 den_old .ne. den .or. 2 zbar_old .ne. zbar .or. 3 abar_old .ne. abar ) then temp_old = temp den_old = den zbar_old = zbar abar_old = abar ytot = 1.0d0/abar rr = den * ytot tempi = 1.0d0/temp dtempi = -tempi*tempi deni = 1.0d0/den pp = sqrt(rr*tempi*(z2bar + zbar)) qq = 0.5d0/pp *(z2bar + zbar) dppdt = qq*rr*dtempi dppdd = qq*ytot*tempi qlam0z = 1.88d8 * tempi * pp qlam0zdt = 1.88d8 * (dtempi*pp + tempi*dppdt) qlam0zdd = 1.88d8 * tempi * dppdd taufac = co2 * tempi**x13 taufacdt = -x13*taufac*tempi qq = rr*zbar xni = qq**x13 dxnidd = x13 * xni * deni aa = 2.27493d5 * tempi * xni daadt = 2.27493d5 * dtempi * xni daadd = 2.27493d5 * tempi * dxnidd end if c..calculate individual screening factors bb = z1 * z2 gamp = aa gampdt = daadt gampdd = daadd qq = fact * bb * zs13inv(jscreen) gamef = qq * gamp gamefdt = qq * gampdt gamefdd = qq * gampdd tau12 = taufac * aznut(jscreen) tau12dt = taufacdt * aznut(jscreen) qq = 1.0d0/tau12 alph12 = gamef * qq alph12dt = (gamefdt - alph12*tau12dt) * qq alph12dd = gamefdd * qq c..limit alph12 to 1.6 to prevent unphysical behavior. c..this should really be replaced by a pycnonuclear reaction rate formula if (alph12 .gt. 1.6) then alph12 = 1.6d0 alph12dt = 0.0d0 alph12dd = 0.0d0 gamef = 1.6d0 * tau12 gamefdt = 1.6d0 * tau12dt gamefdd = 0.0d0 qq = zs13(jscreen)/(fact * bb) gamp = gamef * qq gampdt = gamefdt * qq gampdd = 0.0d0 end if c..weak screening regime h12w = bb * qlam0z dh12wdt = bb * qlam0zdt dh12wdd = bb * qlam0zdd h12 = h12w dh12dt = dh12wdt dh12dd = dh12wdd c..intermediate and strong sceening regime if (gamef .gt. 0.3) then gamp14 = gamp**x14 rr = 1.0d0/gamp qq = 0.25d0*gamp14*rr gamp14dt = qq * gampdt gamp14dd = qq * gampdd cc = 0.896434d0 * gamp * zhat(jscreen) 1 - 3.44740d0 * gamp14 * zhat2(jscreen) 2 - 0.5551d0 * (log(gamp) + lzav(jscreen)) 3 - 2.996d0 dccdt = 0.896434d0 * gampdt * zhat(jscreen) 1 - 3.44740d0 * gamp14dt * zhat2(jscreen) 2 - 0.5551d0*rr*gampdt dccdd = 0.896434d0 * gampdd * zhat(jscreen) 1 - 3.44740d0 * gamp14dd * zhat2(jscreen) 2 - 0.5551d0*rr*gampdd a3 = alph12 * alph12 * alph12 da3 = 3.0d0 * alph12 * alph12 qq = 0.014d0 + 0.0128d0*alph12 dqqdt = 0.0128d0*alph12dt dqqdd = 0.0128d0*alph12dd rr = x532 - alph12*qq drrdt = -(alph12dt*qq + alph12*dqqdt) drrdd = -(alph12dd*qq + alph12*dqqdd) ss = tau12*rr dssdt = tau12dt*rr + tau12*drrdt dssdd = tau12*drrdd tt = -0.0098d0 + 0.0048d0*alph12 dttdt = 0.0048d0*alph12dt dttdd = 0.0048d0*alph12dd uu = 0.0055d0 + alph12*tt duudt = alph12dt*tt + alph12*dttdt duudd = alph12dd*tt + alph12*dttdd vv = gamef * alph12 * uu dvvdt= gamefdt*alph12*uu + gamef*alph12dt*uu + gamef*alph12*duudt dvvdd= gamefdd*alph12*uu + gamef*alph12dd*uu + gamef*alph12*duudd h12 = cc - a3 * (ss + vv) rr = da3 * (ss + vv) dh12dt = dccdt - rr*alph12dt - a3*(dssdt + dvvdt) dh12dd = dccdd - rr*alph12dd - a3*(dssdd + dvvdd) rr = 1.0d0 - 0.0562d0*a3 ss = -0.0562d0*da3 drrdt = ss*alph12dt drrdd = ss*alph12dd if (rr .ge. 0.77d0) then xlgfac = rr dxlgfacdt = drrdt dxlgfacdd = drrdd else xlgfac = 0.77d0 dxlgfacdt = 0.0d0 dxlgfacdd = 0.0d0 end if h12 = log(xlgfac) + h12 rr = 1.0d0/xlgfac dh12dt = rr*dxlgfacdt + dh12dt dh12dd = rr*dxlgfacdd + dh12dd if (gamef .le. 0.8) then rr = 2.0d0*(0.8d0-gamef) drrdt = -2.0d0*gamefdt drrdd = -2.0d0*gamefdd ss = 2.0d0*(gamef-0.3d0) dssdt = 2.0d0*gamefdt dssdd = 2.0d0*gamefdd vv = h12 h12 = h12w*rr + vv*ss dh12dt = dh12wdt*rr + h12w*drrdt + dh12dt*ss + vv*dssdt dh12dd = dh12wdd*rr + h12w*drrdd + dh12dd*ss + vv*dssdd end if c..end of intermediate and strong screening if end if c..machine limit the output h12 = max(min(h12,300.0d0),0.0d0) scor = exp(h12) if (h12 .eq. 300.0d0) then scordt = 0.0d0 scordd = 0.0d0 else scordt = scor * dh12dt scordd = scor * dh12dd end if c write(6,111) 'weak =',h12w,' total =',h12, c 1 ' 1-ratio =',1.0d0-h12w/h12,' correction',scor c 111 format(1x,4(a,1pe13.6)) c read(5,*) return end double precision function snupp(yp,ratepp,ybe7,ratebeec, 1 yb8,rateb8epnu) include 'implno.dek' include 'const.dek' c..computes approximate neutrino losses from pp chain reactions c..see page 142 of astro 289j notes for these loss formulas c..input: c..yp = proton molar abbundance c..ratepp = pp reaction rate c..ybe7 = be7 molar abundance c..ratebeec = be7 electron capture reaction rate c..yb8 = b8 molar abundance c..rateb8epnu = b8 decay reaction rate c..declare the pass double precision yp,ratepp,ybe7,ratebeec,yb8,rateb8epnu c..local variables double precision pp1nu,pp2nu,pp3nu,conv parameter (conv = ev2erg*1.0d6*avo) c..nu losses from p(p,e-nu)h2 pp1nu = yp*yp*ratepp * 0.5d0 * 0.263d0 c..nu losses from be7(n=>p)li7 pp2nu = ybe7 * ratebeec * 0.81d0 c..nu losses from b8(p=>n)be8=>2a pp3nu = yb8 * rateb8epnu * 7.73d0 c..sum the pp-chain neutrino losses and convert to erg/g/s snupp = (pp1nu + pp2nu + pp3nu) * conv return end double precision function snucno(yn13,bc13,bn13,yo14,bn14,bo14, 1 yo15,bn15,bo15,yf17,bo17,bf17, 2 yf18,bo18,bf18) include 'implno.dek' include 'const.dek' c..computes approximate neutrino losses from cno cycle reactions c..see page 142 of astro 289j notes for these loss formulas c..input: c..yn13 = n13 molar abundance c..bc13 = c13 binding energy in mev c..bn13 = n13 binding energy in mev c..yo14 = o14 molar abundance c..bn14 = n14 binding energy in mev c..bo14 = o14 binding energy in mev c..yo15 = o15 molar abundance c..bn15 = n15 binding energy in mev c..bo15 = o15 binding energy in mev c..yf17 = f17 molar abundance c..bo17 = o17 binding energy in mev c..bf17 = f17 binding energy in mev c..yf18 = f18 molar abundance c..bo18 = o18 binding energy in mev c..bf18 = f18 binding energy in mev c..declare the pass double precision yn13,bc13,bn13,yo14,bn14,bo14, 1 yo15,bn15,bo15,yf17,bo17,bf17, 2 yf18,bo18,bf18 c..local variables double precision sum,sum2,enu13n,enu14o,enu15o,enu17f,enu18f, 1 conv,lntwo,tm1,tm2,tm3,tm4,tm5 parameter (conv = ev2erg*1.0d6*avo, 1 lntwo = 0.693147181d0, 2 tm1 = lntwo/597.9d0, 3 tm2 = lntwo/70.606d0, 4 tm3 = lntwo/124.0, 5 tm4 = lntwo/64.49, 6 tm5 = lntwo/6586.2) c..13n(e+nu)13c sum = bc13 - bn13 - 0.782d0 - 1.022d0 sum = 1.0d0 + sum/0.511d0 sum2 = sum*sum enu13n = 0.5d0 * sum * 0.511d0 * (1.0d0 - 1.0d0/sum2) 1 * (1.0d0 - 1.0d0/(4.0d0*sum) - 1.0d0/(9.0d0*sum2)) enu13n = yn13 * enu13n * tm1 c..hot cno cycle 14o(e+nu)14n sum = bn14 - bo14 - 0.782d0 - 1.022d0 sum = 1.0d0 + sum/0.511d0 sum2 = sum*sum enu14o = 0.5d0 * sum * 0.511d0 * (1.0d0 - 1.0d0/sum2) 1 * (1.0d0 - 1.0d0/(4.0d0*sum) - 1.0d0/(9.0d0*sum2)) enu14o = yo14 * enu14o * tm2 c..15o(e+nu)15n sum = bn15 - bo15 - 0.782d0 - 1.022d0 sum = 1.0d0 + sum/0.511d0 sum2 = sum*sum enu15o = 0.5d0 * sum * 0.511d0 * (1.0d0 - 1.0d0/sum2) 1 * (1.0d0 - 1.0d0/(4.0d0*sum) - 1.0d0/(9.0d0*sum2)) enu15o = yo15 * enu15o * tm3 c..17f(e+nu)17o sum = bo17 - bf17 - 0.782d0 - 1.022d0 sum = 1.0d0 + sum/0.511d0 sum2 = sum*sum enu17f = 0.5d0 * sum * 0.511d0 * (1.0d0 - 1.0d0/sum2) 1 * (1.0d0 - 1.0d0/(4.0d0*sum) - 1.0d0/(9.0d0*sum2)) enu17f = yf17 * enu17f * tm4 c..18f(e+nu)18o sum = bo18 - bf18 - 0.782d0 - 1.022d0 sum = 1.0d0 + sum/0.511d0 sum2 = sum*sum enu18f = 0.5d0 * sum * 0.511d0 * (1.0d0 - 1.0d0/sum2) 1 * (1.0d0 - 1.0d0/(4.0d0*sum) - 1.0d0/(9.0d0*sum2)) enu18f = yf18 * enu18f * tm5 c..sum the cno cycle losses and convert to erg/g/s snucno = (enu13n + enu14o + enu15o + enu17f + enu18f) * conv return end subroutine sneut5(temp,den,abar,zbar, 1 snu,dsnudt,dsnudd,dsnuda,dsnudz) include 'implno.dek' include 'const.dek' c..this routine computes neutrino losses from the analytic fits of c..itoh et al. apjs 102, 411, 1996, and also returns their derivatives. c..input: c..temp = temperature c..den = density c..abar = mean atomic weight c..zbar = mean charge c..output: c..snu = total neutrino loss rate in erg/g/sec c..dsnudt = derivative of snu with temperature c..dsnudd = derivative of snu with density c..dsnuda = derivative of snu with abar c..dsnudz = derivative of snu with zbar c..declare the pass double precision temp,den,abar,zbar, 1 snu,dsnudt,dsnudd,dsnuda,dsnudz c..local variables double precision spair,spairdt,spairdd,spairda,spairdz, 1 splas,splasdt,splasdd,splasda,splasdz, 2 sphot,sphotdt,sphotdd,sphotda,sphotdz, 3 sbrem,sbremdt,sbremdd,sbremda,sbremdz, 4 sreco,srecodt,srecodd,srecoda,srecodz double precision t9,xl,xldt,xlp5,xl2,xl3,xl4,xl5,xl6,xl7,xl8,xl9, 1 xlmp5,xlm1,xlm2,xlm3,xlm4,xlnt,cc,den6,tfermi, 2 a0,a1,a2,a3,b1,b2,c00,c01,c02,c03,c04,c05,c06, 3 c10,c11,c12,c13,c14,c15,c16,c20,c21,c22,c23,c24, 4 c25,c26,dd00,dd01,dd02,dd03,dd04,dd05,dd11,dd12, 5 dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,b,c,d,f0, 6 f1,deni,tempi,abari,zbari,f2,f3,z,xmue,ye, 7 dum,dumdt,dumdd,dumda,dumdz, 8 gum,gumdt,gumdd,gumda,gumdz c..pair production double precision rm,rmdd,rmda,rmdz,rmi,gl,gldt, 1 zeta,zetadt,zetadd,zetada,zetadz,zeta2,zeta3, 2 xnum,xnumdt,xnumdd,xnumda,xnumdz, 3 xden,xdendt,xdendd,xdenda,xdendz, 4 fpair,fpairdt,fpairdd,fpairda,fpairdz, 5 qpair,qpairdt,qpairdd,qpairda,qpairdz c..plasma double precision gl2,gl2dt,gl2dd,gl2da,gl2dz,gl12,gl32,gl72,gl6, 1 ft,ftdt,ftdd,ftda,ftdz,fl,fldt,fldd,flda,fldz, 2 fxy,fxydt,fxydd,fxyda,fxydz c..photo double precision tau,taudt,cos1,cos2,cos3,cos4,cos5,sin1,sin2, 1 sin3,sin4,sin5,last,xast, 2 fphot,fphotdt,fphotdd,fphotda,fphotdz, 3 qphot,qphotdt,qphotdd,qphotda,qphotdz c..brem double precision t8,t812,t832,t82,t83,t85,t86,t8m1,t8m2,t8m3,t8m5, 1 t8m6, 2 eta,etadt,etadd,etada,etadz,etam1,etam2,etam3, 3 fbrem,fbremdt,fbremdd,fbremda,fbremdz, 4 gbrem,gbremdt,gbremdd,gbremda,gbremdz, 5 u,gm1,gm2,gm13,gm23,gm43,gm53,v,w,fb,gt,gb, 6 fliq,fliqdt,fliqdd,fliqda,fliqdz, 7 gliq,gliqdt,gliqdd,gliqda,gliqdz c..recomb double precision ifermi12,zfermim12,nu,nudt,nudd,nuda,nudz, 1 nu2,nu3,bigj,bigjdt,bigjdd,bigjda,bigjdz c..numerical constants double precision fac1,fac2,fac3,oneth,twoth,con1,sixth,iln10 parameter (fac1 = 5.0d0 * pi / 3.0d0, 2 fac2 = 10.0d0 * pi, 3 fac3 = pi / 5.0d0, 4 oneth = 1.0d0/3.0d0, 5 twoth = 2.0d0/3.0d0, 6 con1 = 1.0d0/5.9302d0, 7 sixth = 1.0d0/6.0d0, 8 iln10 = 4.342944819032518d-1) c..theta is sin**2(theta_weinberg) = 0.2319 plus/minus 0.00005 (1996) c..xnufam is the number of neutrino flavors = 3.02 plus/minus 0.005 (1998) c..change theta and xnufam if need be, and the changes will automatically c..propagate through the routine. cv and ca are the vector and axial currents. double precision theta,xnufam,cv,ca,cvp,cap,tfac1,tfac2,tfac3, 1 tfac4,tfac5,tfac6 parameter (theta = 0.2319d0, 1 xnufam = 3.0d0, 2 cv = 0.5d0 + 2.0d0 * theta, 3 cvp = 1.0d0 - cv, 4 ca = 0.5d0, 5 cap = 1.0d0 - ca, 6 tfac1 = cv*cv + ca*ca + 7 (xnufam-1.0d0) * (cvp*cvp+cap*cap), 8 tfac2 = cv*cv - ca*ca + 9 (xnufam-1.0d0) * (cvp*cvp - cap*cap), & tfac3 = tfac2/tfac1, 1 tfac4 = 0.5d0 * tfac1, 2 tfac5 = 0.5d0 * tfac2, 3 tfac6 = cv*cv + 1.5d0*ca*ca + (xnufam - 1.0d0)* 4 (cvp*cvp + 1.5d0*cap*cap)) c..initialize spair = 0.0d0 spairdt = 0.0d0 spairdd = 0.0d0 spairda = 0.0d0 spairdz = 0.0d0 splas = 0.0d0 splasdt = 0.0d0 splasdd = 0.0d0 splasda = 0.0d0 splasdz = 0.0d0 sphot = 0.0d0 sphotdt = 0.0d0 sphotdd = 0.0d0 sphotda = 0.0d0 sphotdz = 0.0d0 sbrem = 0.0d0 sbremdt = 0.0d0 sbremdd = 0.0d0 sbremda = 0.0d0 sbremdz = 0.0d0 sreco = 0.0d0 srecodt = 0.0d0 srecodd = 0.0d0 srecoda = 0.0d0 srecodz = 0.0d0 snu = 0.0d0 dsnudt = 0.0d0 dsnudd = 0.0d0 dsnuda = 0.0d0 dsnudz = 0.0d0 if (temp .lt. 1.0e7) return c..to avoid lots of divisions deni = 1.0d0/den tempi = 1.0d0/temp abari = 1.0d0/abar zbari = 1.0d0/zbar c..some composition variables ye = zbar*abari xmue = abar*zbari c..some frequent factors t9 = temp * 1.0d-9 xl = t9 * con1 xldt = 1.0d-9 * con1 xlp5 = sqrt(xl) xl2 = xl*xl xl3 = xl2*xl xl4 = xl3*xl xl5 = xl4*xl xl6 = xl5*xl xl7 = xl6*xl xl8 = xl7*xl xl9 = xl8*xl xlmp5 = 1.0d0/xlp5 xlm1 = 1.0d0/xl xlm2 = xlm1*xlm1 xlm3 = xlm1*xlm2 xlm4 = xlm1*xlm3 rm = den*ye rmdd = ye rmda = -rm*abari rmdz = den*abari rmi = 1.0d0/rm a0 = rm * 1.0d-9 a1 = a0**oneth zeta = a1 * xlm1 zetadt = -a1 * xlm2 * xldt a2 = oneth * a1*rmi * xlm1 zetadd = a2 * rmdd zetada = a2 * rmda zetadz = a2 * rmdz zeta2 = zeta * zeta zeta3 = zeta2 * zeta c..pair neutrino section c..for reactions like e+ + e- => nu_e + nubar_e c..equation 2.8 gl = 1.0d0 - 13.04d0*xl2 +133.5d0*xl4 +1534.0d0*xl6 +918.6d0*xl8 gldt = xldt*(-26.08d0*xl +534.0d0*xl3 +9204.0d0*xl5 +7348.8d0*xl7) c..equation 2.7 a1 = 6.002d19 + 2.084d20*zeta + 1.872d21*zeta2 a2 = 2.084d20 + 2.0d0*1.872d21*zeta if (t9 .lt. 10.0) then b1 = exp(-5.5924d0*zeta) b2 = -b1*5.5924d0 else b1 = exp(-4.9924d0*zeta) b2 = -b1*4.9924d0 end if xnum = a1 * b1 c = a2*b1 + a1*b2 xnumdt = c*zetadt xnumdd = c*zetadd xnumda = c*zetada xnumdz = c*zetadz if (t9 .lt. 10.0) then a1 = 9.383d-1*xlm1 - 4.141d-1*xlm2 + 5.829d-2*xlm3 a2 = -9.383d-1*xlm2 + 2.0d0*4.141d-1*xlm3 - 3.0d0*5.829d-2*xlm4 else a1 = 1.2383d0*xlm1 - 8.141d-1*xlm2 a2 = -1.2383d0*xlm2 + 2.0d0*8.141d-1*xlm3 end if b1 = 3.0d0*zeta2 xden = zeta3 + a1 xdendt = b1*zetadt + a2*xldt xdendd = b1*zetadd xdenda = b1*zetada xdendz = b1*zetadz a1 = 1.0d0/xden fpair = xnum*a1 fpairdt = (xnumdt - fpair*xdendt)*a1 fpairdd = (xnumdd - fpair*xdendd)*a1 fpairda = (xnumda - fpair*xdenda)*a1 fpairdz = (xnumdz - fpair*xdendz)*a1 c..equation 2.6 a1 = 10.7480d0*xl2 + 0.3967d0*xlp5 + 1.005d0 a2 = xldt*(2.0d0*10.7480d0*xl + 0.5d0*0.3967d0*xlmp5) xnum = 1.0d0/a1 xnumdt = -xnum*xnum*a2 a1 = 7.692d7*xl3 + 9.715d6*xlp5 a2 = xldt*(3.0d0*7.692d7*xl2 + 0.5d0*9.715d6*xlmp5) c = 1.0d0/a1 b1 = 1.0d0 + rm*c xden = b1**(-0.3d0) d = -0.3d0*xden/b1 xdendt = -d*rm*c*c*a2 xdendd = d*rmdd*c xdenda = d*rmda*c xdendz = d*rmdz*c qpair = xnum*xden qpairdt = xnumdt*xden + xnum*xdendt qpairdd = xnum*xdendd qpairda = xnum*xdenda qpairdz = xnum*xdendz c..equation 2.5 a1 = exp(-2.0d0*xlm1) a2 = a1*2.0d0*xlm2*xldt spair = a1*fpair spairdt = a2*fpair + a1*fpairdt spairdd = a1*fpairdd spairda = a1*fpairda spairdz = a1*fpairdz a1 = spair spair = gl*a1 spairdt = gl*spairdt + gldt*a1 spairdd = gl*spairdd spairda = gl*spairda spairdz = gl*spairdz a1 = tfac4*(1.0d0 + tfac3 * qpair) a2 = tfac4*tfac3 a3 = spair spair = a1*a3 spairdt = a1*spairdt + a2*qpairdt*a3 spairdd = a1*spairdd + a2*qpairdd*a3 spairda = a1*spairda + a2*qpairda*a3 spairdz = a1*spairdz + a2*qpairdz*a3 c..plasma neutrino section c..for collective reactions like gamma_plasmon => nu_e + nubar_e c..equation 4.6 a1 = 1.019d-6*rm a2 = a1**twoth a3 = twoth*a2/a1 b1 = sqrt(1.0d0 + a2) b2 = 1.0d0/b1 c00 = 1.0d0/(temp*temp*b1) gl2 = 1.1095d11 * rm * c00 gl2dt = -2.0d0*gl2*tempi d = rm*c00*b2*0.5d0*b2*a3*1.019d-6 gl2dd = 1.1095d11 * (rmdd*c00 - d*rmdd) gl2da = 1.1095d11 * (rmda*c00 - d*rmda) gl2dz = 1.1095d11 * (rmdz*c00 - d*rmdz) gl = sqrt(gl2) gl12 = sqrt(gl) gl32 = gl * gl12 gl72 = gl2 * gl32 gl6 = gl2 * gl2 * gl2 c..equation 4.7 ft = 2.4d0 + 0.6d0*gl12 + 0.51d0*gl + 1.25d0*gl32 gum = 1.0d0/gl2 a1 =(0.25d0*0.6d0*gl12 +0.5d0*0.51d0*gl +0.75d0*1.25d0*gl32)*gum ftdt = a1*gl2dt ftdd = a1*gl2dd ftda = a1*gl2da ftdz = a1*gl2dz c..equation 4.8 a1 = 8.6d0*gl2 + 1.35d0*gl72 a2 = 8.6d0 + 1.75d0*1.35d0*gl72*gum b1 = 225.0d0 - 17.0d0*gl + gl2 b2 = -0.5d0*17.0d0*gl*gum + 1.0d0 c = 1.0d0/b1 fl = a1*c d = (a2 - fl*b2)*c fldt = d*gl2dt fldd = d*gl2dd flda = d*gl2da fldz = d*gl2dz c..equation 4.9 and 4.10 cc = log10(2.0d0*rm) xlnt = log10(temp) xnum = sixth * (17.5d0 + cc - 3.0d0*xlnt) xnumdt = -iln10*0.5d0*tempi a2 = iln10*sixth*rmi xnumdd = a2*rmdd xnumda = a2*rmda xnumdz = a2*rmdz xden = sixth * (-24.5d0 + cc + 3.0d0*xlnt) xdendt = iln10*0.5d0*tempi xdendd = a2*rmdd xdenda = a2*rmda xdendz = a2*rmdz c..equation 4.11 if (abs(xnum) .gt. 0.7d0 .or. xden .lt. 0.0d0) then fxy = 1.0d0 fxydt = 0.0d0 fxydd = 0.0d0 fxydz = 0.0d0 fxyda = 0.0d0 else a1 = 0.39d0 - 1.25d0*xnum - 0.35d0*sin(4.5d0*xnum)