program drive_aprox13 include 'implno.dek' include 'const.dek' include 'vector_eos.dek' include 'burn_common.dek' include 'network.dek' include 'cjdet.dek' c..this program exercises the aprox13 network c..declare character*40 hfile logical lfile integer i,j,nok,nbad,nburn,iind(13),niter,kkase double precision tstep,conserv, 1 tin,din,ein,xin(13),tout,dout,eout,xout(13), 2 xneut,xh1,xhe4,xc12,xn14,xo16,xsi28,xfe52,xni56, 3 xdum,xdumnew,f,df,ratio,sum,z,a, 4 xelem,andgrev,ydum(13),abar,zbar, 5 qdum,vin,zin c..formats 01 format(1x,a,a,a) 02 format(1x,a4,'=',1pe10.3,' ',a4,'=',1pe10.3,' ', 1 a4,'=',1pe10.3,' ',a4,'=',1pe10.3,' ', 2 a4,'=',1pe10.3) 03 format(a) 04 format(1x,a,':',/, 1 1x,3(a,1pe20.12),/, 2 1x,3(a,1pe20.12),/, 3 1x,a,1pe11.3,2(a,i5)) 09 format(1x,a,i2,a) c..set the evolution histories write flag lfile = .true. hfile = 'aprox13_a' c..set the detonation run flag detonation = .true. c..set the upsteam temperature, denisty, and composition tstep = 1.0e30 tin = 5.0e7 din = 5.0e6 xhe4 = 1.0d0 c..set the peak expansion temperature if needed c..psi = 1 is an adiabatic expansion c..psi = -1 is an adiabatic implosion if (expansion) then psi = 1.0d0 c psi = -1.0d0 den0 = din temp0 = tin temp_stop = 1.0d7 c temp_stop = 1.0d10 if ( (psi .eq. 1.0 .and. temp_stop .ge. tin) .or. 1 (psi .eq. -1.0 .and. temp_stop .le. tin)) 2 stop 'bad adiabatic temp_stop in routine burner' else psi = 0.0d0 temp_stop = 1.0d30 end if c..screening and table options screen_on = 1 use_tables = 0 c..initialize the network call init_aprox13 c..form the input composition c..initialize do i=1,ionmax xin(i) = 1.0d-30 end do c..give it a solar seed if requested if (xh1 .lt. 0.0) 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..otherwise set the composition variables else if (idnet .eq. idaprox13 .or. idnet .eq. idaprox19 .or. 1 idnet .eq. idpp123 .or. idnet .eq. idcno .or. 2 idnet .eq. idhotcno .or. idnet .eq. idppcno .or. 3 idnet .eq. idpphotcno .or. idnet .eq. idwwrp) then if (ih1 .ne. 0) xin(ih1) = xh1 else if (iprot .ne. 0) xin(iprot) = xh1 end if if (ineut .ne. 0) xin(ineut) = xneut 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 (isi28 .ne. 0) xin(isi28) = xsi28 if (ife52 .ne. 0) xin(ife52) = xfe52 if (ini56 .ne. 0) xin(ini56) = xni56 endif 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..set the isotope name stop parameter and c..halt the integration when the name_stop mass fraction is below xmass_stop name_stop = 'he4 ' xmass_stop = -1.0d30 c..be sure the isotope is in the network do i=1,ionmax if (ionam(i) .eq. name_stop) then id_stop = i goto 11 end if enddo stop 'name_stop not in network in routine tmenet' 11 continue c..get the chapman-jouget detonation 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,iind) write(6,02) (ionam(iind(i)), 1 xmass_cj(iind(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) end if c..reset the initial conditions for a znd detonation if (detonation) then tin = temp_sh din = den_sh vin = vel_sh zin = 1.0d0 + 1.0e-16*vel_sh den_stop = 1.0d0 * 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) ) c..otherwise set the position and velocity to zero else vin = 0.0d0 zin = 0.0d0 end if c..get the internal energy corresponding to the c..starting temperature,density, and composition call azbar(xin,aion,zion,ionmax, 1 ydum,abar,zbar) c..call an eos temp_row(1) = tin den_row(1) = din abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos ein = etot_row(1) c..a message write(6,*) write(6,*) 'starting integration' write(6,*) c..burn it nburn = 1 do i=1,nburn call burner(tstep, 1 tin,din,vin,zin,ein,xin, 2 tout,dout,eout,xout, 3 conserv,nok,nbad, 4 lfile,hfile) enddo c..output write(6,*) ' ' write(6,04) netname, 1 ' tin =',tin,' din =',din,' ein =',ein, 2 ' tout=',tout,' dout=',dout,' eout=',eout, 3 ' sum =',conserv,' nbad=',nbad,' nok=',nok write(6,*) ' ' c..write out the biggest mass fractions call indexx(ionmax,xout,iind) j = min(20,ionmax) write(6,09) 'top ',j,' mass fractions:' j = max(ionmax-19,1) write(6,02) (ionam(iind(i)),xout(iind(i)), i=ionmax,j,-1) write(6,*) ' ' end subroutine burner(tstep, 1 tin,din,vin,zin,ein,xin, 2 tout,dout,eout,xout, 3 conserv,nok,nbad, 4 lfile,hfile) include 'implno.dek' include 'const.dek' include 'vector_eos.dek' include 'burn_common.dek' include 'network.dek' include 'cjdet.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 character*40 hfile logical lfile integer nok,nbad double precision tstep,tin,din,vin,zin,ein,xin(1), 1 tout,dout,eout,xout(1),conserv c..local variables integer i character*8 atim character*9 adat double precision xcons,yex,ydum(neqs), 1 dydt_dum(neqs),xdum(neqs),abar,zbar double precision conv parameter (conv = ev2erg*1.0d6*avo) external aprox13,saprox13,baprox13,daprox13 external stifbs_leqs c..tdim sets the size of the storage arrays c..keep tdim small if evolutions are not being written out c..and make tdim larger if evolution are being written out integer tdim,kount,iprint,kstore c parameter (tdim=10, iprint=0, kstore=0) c parameter (tdim=10000, iprint=0, kstore=tdim) parameter (tdim=1000, iprint=0, kstore=tdim) c..for the integration driver double precision beg,stptry,stpmin,tend,dxsav, 1 ys2(neqs),ttime(tdim),elem(neqs,tdim) parameter (dxsav = 0.0d0) c..for easy tolerance changes double precision odescal,tol parameter (tol = 1.0d-6, 1 odescal = 1.0d-6) c..for writing an evolution history character*80 string integer j,lop,ilop,jrem,kb,ke,k double precision sum c..popular format statements 01 format(1x,'*',t9,a,t20,a,t31,a,t42,a,t53,a,t64,a, 1 t75,a,t86,a,t97,a,t108,a,t119,a) 02 format(1x,a4,' =',1pe11.3,' ',a4,' =',1pe11.3,' ', 1 a4,' =',1pe11.3,' ',a4,' =',1pe11.3) 03 format(a30,i8,a4) 04 format(1x,1p15e11.3) 05 format(1x,i4,1p11e11.3) 06 format(1x,'* ',a,': ',2(a,i5)) 07 format(1x,'* ',a,4(a,1pe11.3)) 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 ys2(iener) = ein ys2(itemp) = tin ys2(iden) = din ys2(ivelx) = vin ys2(iposx) = zin 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 = 1.0d-16 stpmin = stptry * 1.0d-12 end if c..zero the output array do j=1,tdim do i=1,neqs elem(i,j) = 0.0d0 enddo enddo c..integrate the aprox13 network call netint(beg,stptry,stpmin,tend,ys2, 1 tol,dxsav,kstore, 2 ttime,elem,tdim,neqs,tdim,neqs, 3 nok,nbad,kount,odescal,iprint, 4 aprox13,daprox13,baprox13,stifbs_leqs) c..set the output composition do i=1,ionmax xout(i) = ys2(i) * aion(i) enddo c..output temperature, density, and thermal energy tout = elem(itemp,kount) dout = elem(iden,kount) eout = elem(iener,kount) c..also output the mass non-conservation conserv = 0.0d0 do i=1,ionmax conserv = conserv + xout(i) enddo conserv = 1.0d0 - conserv c..if the evolution history is to be written out if (lfile) then c..write the above quantities write(string,03) hfile,0,'.dat' call sqeeze(string) call today(adat,atim) open (unit=2, file=string, status='unknown') write(string,03) hfile,1,'.dat' call sqeeze(string) open (unit=3, file=string, status='unknown') write(2,01) adat,atim write(3,01) adat,atim if (one_step) then write(2,07) 'one_step:',' btemp=',btemp,' bden=',bden write(3,07) 'one_step:',' btemp=',btemp,' bden=',bden else if (hydrostatic) then write(2,07) 'hydrostatic:',' btemp=',btemp,' bden=',bden write(3,07) 'hydrostatic:',' btemp=',btemp,' bden=',bden else if (expansion) then write(2,07) 'expansion:',' temp0=',temp0,' den0=',den0, 1 ' temp_stop=',temp_stop write(3,07) 'expansion:',' temp0=',temp0,' den0=',den0, 1 ' temp_stop=',temp_stop else if (self_heat) then write(2,07) 'self_heat:',' temp0=',elem(itemp,1), 1 ' den0=',elem(iden,1) write(3,07) 'self_heat:',' temp0=',elem(itemp,1), 1 ' den0=',elem(iden,1) else if (detonation) then write(2,07) 'detonation:',' temp0=',temp_up, 1 ' den0=',den_up, 2 ' pres0=',pres_up, 3 ' mach=',mach_sh write(3,07) 'detonation:',' temp0=',temp_up, 1 ' den0=',den_up, 2 ' pres0=',pres_up, 3 ' mach=',mach_sh end if write(2,06) netname,' nbad=',nbad,' nok=',nok write(3,06) netname,' nbad=',nbad,' nok=',nok write(2,01) 'time','temp','den','ener','sdot','sneut', 1 's-snu','ye','1-sum' write(3,01) 'time','pos','vel','temp','den','pres','ener','cs' c..write the cj solution for detonation if (detonation) then i = kount write(3,05) i,ttime(i),elem(iposx,i),vel_cj, 1 temp_cj,den_cj,pres_cj, 2 ener_cj,cs_cj end if c..loop through the output points do i=1,kount c..form the mass fractions do j=1,ionmax xdum(j) = min(1.0d0,max(elem(j,i) * aion(j),1.0d-30)) enddo c..mass conservation sum = 0.0d0 do j=1,ionmax sum = sum + xdum(j) enddo sum = 1.0d0 - sum xcons = sum c..get ye using normalized mass fractions sum = 0.0d0 do j=1,ionmax sum = sum + xdum(j) enddo sum = 1.0d0/sum do j=1,ionmax xdum(j) = min(1.0d0,max(sum * xdum(j),1.0d-30)) enddo c..get abar, zbar and a few other composition variables call azbar(xdum,aion,zion,ionmax, 1 ydum,abar,zbar) yex = zbar/abar c..finish filling the dummy array do j=1,ionmax ydum(j) = elem(j,i) enddo ydum(iener) = elem(iener,i) ydum(itemp) = elem(itemp,i) ydum(iden) = elem(iden,i) c..get the right hand sides, exact energy generation rate and so on call aprox13(ttime(i),ydum,dydt_dum) c..and write what we found write(2,05) i,ttime(i),elem(itemp,i),elem(iden,i), 1 elem(iener,i),sdot,sneut,dydt_dum(iener), 2 yex,xcons c..call an eos temp_row(1) = elem(itemp,i) den_row(1) = elem(iden,i) abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos write(3,05) i,ttime(i),elem(iposx,i),elem(ivelx,i), 1 elem(itemp,i),elem(iden,i),ptot_row(1), 2 elem(iener,i),cs_row(1) c..end of kount loop enddo c..write the cj solution for detonation if (detonation) then i = kount write(3,05) i,ttime(i),elem(iposx,i),vel_cj, 1 temp_cj,den_cj,pres_cj, 2 ener_cj,cs_cj end if c..close up the files close(unit=2) close(unit=3) c..write out the isotopic mass fractions in blocks of 8 c..lop is how many groups of 8 exist; jrem is the remainder lop = ionmax/8 jrem = ionmax - 8*lop do ilop = 1,lop+1 kb = 1 + 8*(ilop-1) ke = 8 + 8*(ilop-1) if (ilop .eq. lop+1 .and. jrem .eq. 0) goto 60 if (ilop .eq. lop+1) ke = ionmax c..open the output file write(string,03) hfile,ilop+1,'.dat' call sqeeze(string) open (unit=2, file=string, status='unknown') write(2,01) adat,atim if (one_step) then write(2,07) 'one_step:',' btemp=',btemp,' bden=',bden else if (hydrostatic) then write(2,07) 'hydrostatic:',' btemp=',btemp,' bden=',bden else if (expansion) then write(2,07) 'expansion:',' temp0=',temp0,' den0=',den0, 1 ' temp_stop=',temp_stop else if (self_heat) then write(2,07) 'self_heat:',' temp0=',elem(itemp,1), 1 ' den0=',elem(iden,1) else if (detonation) then write(2,07) 'detonation:',' temp0=',temp_up, 1 ' den0=',den_up, 1 ' pres0=',pres_up, 2 ' mach=',mach_sh end if write(2,06) netname,' nbad=',nbad,' nok=',nok write(2,01) 'time',(ionam(k), k=kb,ke) do i=1,kount write(2,04) ttime(i),(elem(k,i)*aion(k), k=kb,ke) enddo close(unit=2) 60 continue enddo c..end of lfile if end if return end c..this file contains aprox13 network c..routine aprox13 sets up the odes c..routine faprox13 gets the flows from aprox13 c..routine daprox13 sets up the dense aprox13 jacobian c..routine baprox13 builds the nonzero locations for saprox13 c..routine saprox13 sets up the sparse aprox13 jacobian c..routine aprox13rat generates the reaction rates for routine aprox13 c..routine aprox13tab generates the raw rates using table interpolation c..routine screen_aprox13 applies screening corrections to the raw rates c..routine init_aprox13 initializes the aprox13 network subroutine aprox13(tt,y,dydt) include 'implno.dek' include 'const.dek' include 'burn_common.dek' include 'network.dek' include 'vector_eos.dek' include 'cjdet.dek' c..this routine sets up the system of ode's for the aprox13 nuclear reactions. c..this is an alpha chain + heavy ion network with (a,p)(p,g) links c.. c..isotopes: he4, c12, o16, ne20, mg24, si28, s32, c.. ar36, ca40, ti44, cr48, fe52, ni56 c..declare the pass double precision tt,y(1),dydt(1) c..local variables integer i double precision enuc,taud,taut,z,denom,suma,sumz,ww, 1 zbarxx,ytot1,abar,zbar,velx,posx,cs,dpde, 2 snudt,snudd,snuda,snudz,combo,phi,dtdp,conv parameter (conv = ev2erg*1.0d6*avo) c..positive definite mass fractions do i=1,ionmax y(i) = min(1.0d0,max(y(i),1.0d-30)) enddo c..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..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 c..get the reaction rates if (use_tables .eq. 1) then call aprox13tab else call aprox13rat end if c..do the screening here because the corrections depend on the composition call screen_aprox13(y) c..get the right hand side of the odes call rhs(y,ratdum,dydt) c..instantaneous energy generation rate enuc = 0.0d0 do i=1,ionmax enuc = enuc + dydt(i) * bion(i) enddo enuc = enuc * conv c..get the neutrino losses call sneut5(btemp,bden,abar,zbar, 1 sneut,snudt,snudd,snuda,snudz) c..append an energy equation sdot = enuc 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) then dydt(itemp) = 0.0d0 dydt(iden) = 0.0d0 dydt(ivelx) = 0.0d0 dydt(iposx) = 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 dydt(ivelx) = 0.0d0 dydt(iposx) = 0.0d0 c..self heating else if (self_heat) 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..for de/dy suma = 0.0d0 do i=1,ionmax suma = suma - dydt(i) enddo sumz = 0.0d0 do i=1,ionmax sumz = sumz + (zion(i) - zbar)*dydt(i) enddo c..temperature equation that is self-consistent with an eos z = 1.0d0/cv_row(1) ww = suma*dea_row(1)*abar*abar + sumz*dez_row(1)*abar dydt(itemp) = z*(dydt(iener) - ded_row(1)*dydt(iden) - ww) c..density, velocity, and position equations dydt(iden) = 0.0d0 dydt(ivelx) = 0.0d0 dydt(iposx) = 0.0d0 c..detonation else if (detonation) then c..map the rest of the input vector velx = y(ivelx) posx = y(iposx) 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..for de/dy and dp/dy suma = 0.0d0 do i=1,ionmax suma = suma - dydt(i) enddo sumz = 0.0d0 do i=1,ionmax sumz = sumz + (zion(i) - zbar)*dydt(i) enddo c..the possibly singular denominator cs = cs_row(1) denom = velx*velx - cs*cs c denom = velx*velx - cs_cj*cs_cj c..the function phi dpde = dpt_row(1)/det_row(1) z = suma*dpa_row(1)*abar*abar + sumz*dpz_row(1)*abar ww = suma*dea_row(1)*abar*abar + sumz*dez_row(1)*abar phi = dpde*(dydt(iener) - ww) - z c..a common combination if (denom .ne. 0.0) then combo = phi/denom else combo = 0.0d0 end if c..position equation dydt(iposx) = velx c..density equation dydt(iden) = combo c..velocity equations dydt(ivelx) = -velx/bden*dydt(iden) c..temperature equation dtdp = 1.0d0/dpt_row(1) ww = suma*dpa_row(1)*abar*abar + sumz*dpz_row(1)*abar dydt(itemp) = dtdp*((velx*velx - dpd_row(1))*dydt(iden) - ww) 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 aprox13 odes c..declare the pass double precision y(1),rate(1),dydt(1) c..local variables integer i double precision r1,s1,t1,u1,v1,w1,x1,y1,denom,sixth parameter (sixth = 1.0d0/6.0d0) c..zero the odes do i=1,ionmax dydt(i) = 0.0d0 enddo c..branching ratios for various dummy proton links r1 = 0.0d0 denom = ratdum(iralpa) + ratdum(iralpg) if (denom .ne. 0.0) r1 = ratdum(iralpa)/denom s1 = 0.0d0 denom = ratdum(irppa) + ratdum(irppg) if (denom .ne. 0.0) s1 = ratdum(irppa)/denom t1 = 0.0d0 denom = ratdum(irclpa) + ratdum(irclpg) if (denom .ne. 0.0) t1 = ratdum(irclpa)/denom u1 = 0.0d0 denom = ratdum(irkpa) + ratdum(irkpg) if (denom .ne. 0.0) u1 = ratdum(irkpa)/denom v1 = 0.0d0 denom = ratdum(irscpa) + ratdum(irscpg) if (denom .ne. 0.0) v1 = ratdum(irscpa)/denom w1 = 0.0d0 denom = ratdum(irvpa) + ratdum(irvpg) if (denom .ne. 0.0) w1 = ratdum(irvpa)/denom x1 = 0.0d0 denom = ratdum(irmnpa) + ratdum(irmnpg) if (denom .ne. 0.0) x1 = ratdum(irmnpa)/denom y1 = 0.0d0 denom = ratdum(ircopa) + ratdum(ircopg) if (denom .ne. 0.0) y1 = ratdum(ircopa)/denom c..set up the system of odes: c..he4 reactions c..heavy ion reactions dydt(ihe4) = 0.5d0 * y(ic12) * y(ic12) * rate(ir1212) 1 + 0.5d0 * y(ic12) * y(io16) * rate(ir1216) 2 + 0.56d0 * 0.5d0 * y(io16) * y(io16) * rate(ir1616) 3 + 0.34d0 * s1 * 0.5d0 * y(io16)*y(io16)*rate(ir1616) c..(a,g) and (g,a) reactions dydt(ihe4) = dydt(ihe4) 1 - 0.5d0 * y(ihe4) * y(ihe4) * y(ihe4) * rate(ir3a) 2 + 3.0d0 * y(ic12) * rate(irg3a) 3 - y(ihe4) * y(ic12) * rate(ircag) 4 + y(io16) * rate(iroga) 5 - y(ihe4) * y(io16) * rate(iroag) 6 + y(ine20) * rate(irnega) 7 - y(ihe4) * y(ine20) * rate(irneag) 8 + y(img24) * rate(irmgga) 9 - y(ihe4) * y(img24)* rate(irmgag) & + y(isi28) * rate(irsiga) 1 - y(ihe4) * y(isi28)*rate(irsiag) 2 + y(is32) * rate(irsga) dydt(ihe4) = dydt(ihe4) 1 - y(ihe4) * y(is32) * rate(irsag) 2 + y(iar36) * rate(irarga) 3 - y(ihe4) * y(iar36)*rate(irarag) 4 + y(ica40) * rate(ircaga) 5 - y(ihe4) * y(ica40)*rate(ircaag) 6 + y(iti44) * rate(irtiga) 7 - y(ihe4) * y(iti44)*rate(irtiag) 8 + y(icr48) * rate(ircrga) 9 - y(ihe4) * y(icr48)*rate(ircrag) & + y(ife52) * rate(irfega) 1 - y(ihe4) * y(ife52) * rate(irfeag) 2 + y(ini56) * rate(irniga) c..(a,p)(p,g) and (g,p)(p,a) reactions dydt(ihe4) = dydt(ihe4) 1 - y(ihe4) * y(img24) * rate(irmgap) * (1.0d0-r1) 2 + y(isi28) * rate(irsigp) * r1 3 - y(ihe4) * y(isi28) * rate(irsiap) * (1.0d0-s1) 4 + y(is32) * rate(irsgp) * s1 5 - y(ihe4) * y(is32) * rate(irsap) * (1.0d0-t1) 6 + y(iar36) * rate(irargp) * t1 7 - y(ihe4) * y(iar36) * rate(irarap) * (1.0d0-u1) 8 + y(ica40) * rate(ircagp) * u1 9 - y(ihe4) * y(ica40) * rate(ircaap) * (1.0d0-v1) & + y(iti44) * rate(irtigp) * v1 1 - y(ihe4) * y(iti44) * rate(irtiap) * (1.0d0-w1) 2 + y(icr48) * rate(ircrgp) * w1 3 - y(ihe4) * y(icr48) * rate(ircrap) * (1.0d0-x1) 4 + y(ife52) * rate(irfegp) * x1 5 - y(ihe4) * y(ife52) * rate(irfeap) * (1.0d0-y1) 6 + y(ini56) * rate(irnigp) * y1 c..c12 reactions dydt(ic12) = -y(ic12) * y(ic12) * rate(ir1212) 1 - y(ic12) * y(io16) * rate(ir1216) 2 + sixth * y(ihe4) * y(ihe4) * y(ihe4) * rate(ir3a) 3 - y(ic12) * rate(irg3a) 4 - y(ic12) * y(ihe4) * rate(ircag) 5 + y(io16) * rate(iroga) c..o16 reactions dydt(io16) = -y(ic12) * y(io16) * rate(ir1216) 1 - y(io16) * y(io16) * rate(ir1616) 2 + y(ic12) * y(ihe4) * rate(ircag) 3 - y(io16) * y(ihe4) * rate(iroag) 4 - y(io16) * rate(iroga) 5 + y(ine20) * rate(irnega) c..ne20 reactions dydt(ine20) = 0.5d0 * y(ic12) * y(ic12) * rate(ir1212) 1 + y(io16) * y(ihe4) * rate(iroag) 2 - y(ine20) * y(ihe4) * rate(irneag) 3 - y(ine20) * rate(irnega) 4 + y(img24) * rate(irmgga) c..mg24 reactions dydt(img24) = 0.5d0 * y(ic12) * y(io16) * rate(ir1216) 1 + y(ine20) * y(ihe4) * rate(irneag) 2 - y(img24) * y(ihe4) * rate(irmgag) 4 - y(img24) * rate(irmgga) 5 + y(isi28) * rate(irsiga) 6 - y(img24) * y(ihe4) * rate(irmgap) * (1.0d0-r1) 7 + y(isi28) * r1 * rate(irsigp) c..si28 reactions dydt(isi28) = 0.5d0 * y(ic12) * y(io16) * rate(ir1216) 1 + 0.56d0 * 0.5d0 * y(io16)*y(io16) * rate(ir1616) 2 + 0.34d0 * 0.5d0 * y(io16)*y(io16)*s1*rate(ir1616) 3 + y(img24) * y(ihe4) * rate(irmgag) 4 - y(isi28) * y(ihe4) * rate(irsiag) 5 - y(isi28) * rate(irsiga) 6 + y(is32) * rate(irsga) 7 + y(img24) * y(ihe4) * rate(irmgap) * (1.0d0-r1) 8 - y(isi28) * r1 * rate(irsigp) 9 - y(isi28) * y(ihe4) * rate(irsiap) * (1.0d0-s1) & + y(is32) * s1 * rate(irsgp) c..s32 reactions dydt(is32) = 0.34d0*0.5d0*y(io16)**2 * rate(ir1616)*(1.0d0-s1) 1 + 0.1d0 * 0.5d0 * y(io16) * y(io16) * rate(ir1616) 2 + y(isi28) * y(ihe4) * rate(irsiag) 3 - y(is32) * y(ihe4) * rate(irsag) 4 - y(is32) * rate(irsga) 5 + y(iar36) * rate(irarga) 6 + y(isi28) * y(ihe4) * rate(irsiap) * (1.0d0-s1) 7 - y(is32) * s1 * rate(irsgp) 8 - y(is32) * y(ihe4) * rate(irsap) * (1.0d0-t1) 9 + y(iar36) * t1 * rate(irargp) c..ar36 reactions dydt(iar36) = y(is32) * y(ihe4) * rate(irsag) 1 - y(iar36) * y(ihe4) * rate(irarag) 2 - y(iar36) * rate(irarga) 3 + y(ica40) * rate(ircaga) 4 + y(is32) * y(ihe4) * rate(irsap) * (1.0d0-t1) 5 - y(iar36) * t1 * rate(irargp) 6 - y(iar36) * y(ihe4) * rate(irarap) * (1.0d0-u1) 7 + y(ica40) * rate(ircagp) * u1 c..ca40 reactions dydt(ica40) = y(iar36) * y(ihe4) * rate(irarag) 1 - y(ica40) * y(ihe4) * rate(ircaag) 2 - y(ica40) * rate(ircaga) 3 + y(iti44) * rate(irtiga) 4 + y(iar36) * y(ihe4) * rate(irarap) * (1.0d0-u1) 5 - y(ica40) * rate(ircagp) * u1 6 - y(ica40) * y(ihe4) * rate(ircaap) * (1.0d0-v1) 7 + y(iti44) * rate(irtigp) * v1 c..ti44 reactions dydt(iti44) = y(ica40) * y(ihe4) * rate(ircaag) 1 - y(iti44) * y(ihe4) * rate(irtiag) 2 - y(iti44) * rate(irtiga) 3 + y(icr48) * rate(ircrga) 4 + y(ica40) * y(ihe4) * rate(ircaap)*(1.0d0-v1) 5 - y(iti44) * v1 * rate(irtigp) 6 - y(iti44) * y(ihe4) * rate(irtiap) * (1.0d0-w1) 7 + y(icr48) * w1 * rate(ircrgp) c..cr48 reactions dydt(icr48) = y(iti44) * y(ihe4) * rate(irtiag) 1 - y(icr48) * y(ihe4) * rate(ircrag) 2 - y(icr48) * rate(ircrga) 3 + y(ife52) * rate(irfega) 4 + y(iti44) * y(ihe4) * rate(irtiap)*(1.0d0-w1) 5 - y(icr48) * w1 * rate(ircrgp) 6 - y(icr48) * y(ihe4) * rate(ircrap) * (1.0d0-x1) 7 + y(ife52) * x1 * rate(irfegp) c..fe52 reactions dydt(ife52) = y(icr48) * y(ihe4) * rate(ircrag) 1 - y(ife52) * y(ihe4) * rate(irfeag) 2 - y(ife52) * rate(irfega) 3 + y(ini56) * rate(irniga) 4 + y(icr48) * y(ihe4) * rate(ircrap) * (1.0d0-x1) 5 - y(ife52) * x1 * rate(irfegp) 6 - y(ife52) * y(ihe4) * rate(irfeap) * (1.0d0-y1) 7 + y(ini56) * y1 * rate(irnigp) c..ni56 reactions dydt(ini56) = y(ife52) * y(ihe4) * rate(irfeag) 1 - y(ini56) * rate(irniga) 2 + y(ife52) * y(ihe4) * rate(irfeap) * (1.0d0-y1) 3 - y(ini56) * y1 * rate(irnigp) return end subroutine daprox13(tt,y,dfdy,nlog,nphys) include 'implno.dek' include 'const.dek' include 'burn_common.dek' include 'network.dek' include 'vector_eos.dek' include 'cjdet.dek' c..this routine sets up the dense aprox13 jacobian c..declare the pass integer nlog,nphys double precision tt,y(1),dfdy(nphys,nphys) c..local variables integer i,j double precision denom,denomdt,denomdd, 1 r1,r1dt,r1dd,s1,s1dt,s1dd,t1,t1dt,t1dd, 2 u1,u1dt,u1dd,v1,v1dt,v1dd,w1,w1dt,w1dd, 3 x1,x1dt,x1dd,y1,y1dt,y1dd,zz,xx,yy double precision zbarxx,ytot1,abar,zbar,taud,taut, 1 snudt,snudd,snuda,snudz, 2 dydt(neqs),enuc,velx,posx,suma,sumz,sum, 3 sum5,sum6,sum7,sum8,sum9,sum10,sum11,sum12, 4 cs,denombv,dpde,dpdebd,dpdebt,phi,phibd,phibt, 5 combo,combobv,combobd,combobt,z,zbd,zbt, 6 ww,wwbd,wwbt,dtdp,dtdpbd,dtdpbt, 7 dpdbd,dpdbt,foo(8),moo(8), 8 csbd,csbt,dptbd,dptbt,detbd,detbt,dpabd,dpabt, 9 dpzbd,dpzbt,deabd,deabt,dezbd,dezbt double precision conv parameter (conv = ev2erg*1.0d6*avo) c..zero the jacobian do j=1,nlog do i=1,nlog dfdy(i,j) = 0.0d0 enddo enddo c..positive definite mass fractions do i=1,ionmax y(i) = min(1.0d0,max(y(i),1.0d-30)) enddo c..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..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 c..get the reaction rates if (use_tables .eq. 1) then call aprox13tab else call aprox13rat end if c..do the screening here because the corrections depend on the composition call screen_aprox13(y) c..branching ratios for various dummy proton links r1 = 0.0d0 r1dt = 0.0d0 r1dd = 0.0d0 denom = ratdum(iralpa) + ratdum(iralpg) denomdt = dratdumdt(iralpa) + dratdumdt(iralpg) denomdd = dratdumdd(iralpa) + dratdumdd(iralpg) if (denom .ne. 0.0) then zz = 1.0d0/denom r1 = ratdum(iralpa)*zz r1dt = (dratdumdt(iralpa) - r1*denomdt)*zz r1dd = (dratdumdd(iralpa) - r1*denomdd)*zz end if s1 = 0.0d0 s1dt = 0.0d0 s1dd = 0.0d0 denom = ratdum(irppa) + ratdum(irppg) denomdt = dratdumdt(irppa) + dratdumdt(irppg) denomdd = dratdumdd(irppa) + dratdumdd(irppg) if (denom .ne. 0.0) then zz = 1.0d0/denom s1 = ratdum(irppa)*zz s1dt = (dratdumdt(irppa) - s1*denomdt)*zz s1dd = (dratdumdd(irppa) - s1*denomdd)*zz end if t1 = 0.0d0 t1dt = 0.0d0 t1dd = 0.0d0 denom = ratdum(irclpa) + ratdum(irclpg) denomdt = dratdumdt(irclpa) + dratdumdt(irclpg) denomdd = dratdumdd(irclpa) + dratdumdd(irclpg) if (denom .ne. 0.0) then zz = 1.0d0/denom t1 = ratdum(irclpa)*zz t1dt = (dratdumdt(irclpa) - t1*denomdt)*zz t1dd = (dratdumdd(irclpa) - t1*denomdd)*zz end if u1 = 0.0d0 u1dt = 0.0d0 u1dd = 0.0d0 denom = ratdum(irkpa) + ratdum(irkpg) denomdt = dratdumdt(irkpa) + dratdumdt(irkpg) denomdd = dratdumdd(irkpa) + dratdumdd(irkpg) if (denom .ne. 0.0) then zz = 1.0d0/denom u1 = ratdum(irkpa)*zz u1dt = (dratdumdt(irkpa) - u1*denomdt)*zz u1dd = (dratdumdd(irkpa) - u1*denomdd)*zz end if v1 = 0.0d0 v1dt = 0.0d0 v1dd = 0.0d0 denom = ratdum(irscpa) + ratdum(irscpg) denomdt = dratdumdt(irscpa) + dratdumdt(irscpg) denomdd = dratdumdd(irscpa) + dratdumdd(irscpg) if (denom .ne. 0.0) then zz = 1.0d0/denom v1 = ratdum(irscpa)*zz v1dt = (dratdumdt(irscpa) - v1*denomdt)*zz v1dd = (dratdumdd(irscpa) - v1*denomdd)*zz end if w1 = 0.0d0 w1dt = 0.0d0 w1dd = 0.0d0 denom = ratdum(irvpa) + ratdum(irvpg) denomdt = dratdumdt(irvpa) + dratdumdt(irvpg) denomdd = dratdumdd(irvpa) + dratdumdd(irvpg) if (denom .ne. 0.0) then zz = 1.0d0/denom w1 = ratdum(irvpa)*zz w1dt = (dratdumdt(irvpa) - w1*denomdt)*zz w1dd = (dratdumdd(irvpa) - w1*denomdd)*zz end if x1 = 0.0d0 x1dt = 0.0d0 x1dd = 0.0d0 denom = ratdum(irmnpa) + ratdum(irmnpg) denomdt = dratdumdt(irmnpa) + dratdumdt(irmnpg) denomdd = dratdumdd(irmnpa) + dratdumdd(irmnpg) if (denom .ne. 0.0) then zz = 1.0d0/denom x1 = ratdum(irmnpa)*zz x1dt = (dratdumdt(irmnpa) - x1*denomdt)*zz x1dd = (dratdumdd(irmnpa) - x1*denomdd)*zz endif y1 = 0.0d0 y1dt = 0.0d0 y1dd = 0.0d0 denom = ratdum(ircopa) + ratdum(ircopg) denomdt = dratdumdt(ircopa) + dratdumdt(ircopg) denomdd = dratdumdd(ircopa) + dratdumdd(ircopg) if (denom .ne. 0.0) then zz = 1.0d0/denom y1 = ratdum(ircopa)*zz y1dt = (dratdumdt(ircopa) - y1*denomdt)*zz y1dd = (dratdumdd(ircopa) - y1*denomdd)*zz end if c..he4 jacobian elements 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 - y(isi28) * ratdum(irsiag) 6 - y(is32) * ratdum(irsag) 7 - y(iar36) * ratdum(irarag) 8 - y(ica40) * ratdum(ircaag) 9 - y(iti44) * ratdum(irtiag) & - y(icr48) * ratdum(ircrag) 1 - y(ife52) * ratdum(irfeag) dfdy(ihe4,ihe4) = dfdy(ihe4,ihe4) 1 - y(img24) * ratdum(irmgap) * (1.0d0-r1) 2 - y(isi28) * ratdum(irsiap) * (1.0d0-s1) 3 - y(is32) * ratdum(irsap) * (1.0d0-t1) 4 - y(iar36) * ratdum(irarap) * (1.0d0-u1) 5 - y(ica40) * ratdum(ircaap) * (1.0d0-v1) 6 - y(iti44) * ratdum(irtiap) * (1.0d0-w1) 7 - y(icr48) * ratdum(ircrap) * (1.0d0-x1) 8 - y(ife52) * ratdum(irfeap) * (1.0d0-y1) dfdy(ihe4,ic12) = y(ic12) * ratdum(ir1212) 1 + 0.5d0 * y(io16) * ratdum(ir1216) 2 + 3.0d0 * ratdum(irg3a) 3 - y(ihe4) * ratdum(ircag) dfdy(ihe4,io16) = 0.5d0 * y(ic12) * ratdum(ir1216) 1 + 1.12d0 * 0.5d0 * y(io16) * ratdum(ir1616) 2 + 0.68d0 * s1 * 0.5d0*y(io16) * ratdum(ir1616) 3 + ratdum(iroga) 4 - y(ihe4) * ratdum(iroag) dfdy(ihe4,ine20) = ratdum(irnega) 1 - y(ihe4) * ratdum(irneag) dfdy(ihe4,img24) = ratdum(irmgga) 1 - y(ihe4) * ratdum(irmgag) 2 - y(ihe4) * ratdum(irmgap) * (1.0d0-r1) dfdy(ihe4,isi28) = ratdum(irsiga) 1 - y(ihe4) * ratdum(irsiag) 2 - y(ihe4) * ratdum(irsiap) * (1.0d0-s1) 3 + r1 * ratdum(irsigp) dfdy(ihe4,is32) = ratdum(irsga) 1 - y(ihe4) * ratdum(irsag) 2 - y(ihe4) * ratdum(irsap) * (1.0d0-t1) 3 + s1 * ratdum(irsgp) dfdy(ihe4,iar36) = ratdum(irarga) 1 - y(ihe4) * ratdum(irarag) 2 - y(ihe4) * ratdum(irarap) * (1.0d0-u1) 3 + t1 * ratdum(irargp) dfdy(ihe4,ica40) = ratdum(ircaga) 1 - y(ihe4) * ratdum(ircaag) 2 - y(ihe4) * ratdum(ircaap) * (1.0d0-v1) 3 + u1 * ratdum(ircagp) dfdy(ihe4,iti44) = ratdum(irtiga) 1 - y(ihe4) * ratdum(irtiag) 2 - y(ihe4) * ratdum(irtiap) * (1.0d0-w1) 3 + v1 * ratdum(irtigp) dfdy(ihe4,icr48) = ratdum(ircrga) 1 - y(ihe4) * ratdum(ircrag) 2 - y(ihe4) * ratdum(ircrap) * (1.0d0-x1) 3 + w1 * ratdum(ircrgp) dfdy(ihe4,ife52) = ratdum(irfega) 1 - y(ihe4) * ratdum(irfeag) 2 - y(ihe4) * ratdum(irfeap) * (1.0d0-y1) 3 + x1 * ratdum(irfegp) dfdy(ihe4,ini56) = ratdum(irniga) 1 + y1 * ratdum(irnigp) c..c12 jacobian elements dfdy(ic12,ihe4) = 0.5d0 * y(ihe4) * y(ihe4) * ratdum(ir3a) 1 - y(ic12) * ratdum(ircag) dfdy(ic12,ic12) = -2.0d0 * y(ic12) * ratdum(ir1212) 1 - y(io16) * ratdum(ir1216) 2 - ratdum(irg3a) 3 - y(ihe4) * ratdum(ircag) dfdy(ic12,io16) = -y(ic12) * ratdum(ir1216) 1 + ratdum(iroga) c..o16 jacobian elements dfdy(io16,ihe4) = y(ic12)*ratdum(ircag) 1 - y(io16)*ratdum(iroag) dfdy(io16,ic12) = -y(io16)*ratdum(ir1216) 1 + y(ihe4)*ratdum(ircag) dfdy(io16,io16) = - y(ic12) * ratdum(ir1216) 1 - 2.0d0 * y(io16) * ratdum(ir1616) 2 - y(ihe4) * ratdum(iroag) 3 - ratdum(iroga) dfdy(io16,ine20) = ratdum(irnega) c..ne20 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) = -y(ihe4) * ratdum(irneag) 1 - ratdum(irnega) dfdy(ine20,img24) = ratdum(irmgga) c..mg24 jacobian elements dfdy(img24,ihe4) = y(ine20) * ratdum(irneag) 1 -y(img24) * ratdum(irmgag) 2 -y(img24) * ratdum(irmgap) * (1.0d0-r1) 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) = -y(ihe4) * ratdum(irmgag) 1 - ratdum(irmgga) 2 - y(ihe4) * ratdum(irmgap) * (1.0d0-r1) dfdy(img24,isi28) = ratdum(irsiga) 1 + r1 * ratdum(irsigp) c..si28 jacobian elements dfdy(isi28,ihe4) = y(img24) * ratdum(irmgag) 1 - y(isi28) * ratdum(irsiag) 2 + y(img24) * ratdum(irmgap) * (1.0d0-r1) 3 - y(isi28) * ratdum(irsiap) * (1.0d0-s1) dfdy(isi28,ic12) = 0.5d0 * y(io16) * ratdum(ir1216) dfdy(isi28,io16) = 0.5d0 * y(ic12) * ratdum(ir1216) 1 + 1.12d0 * 0.5d0*y(io16) * ratdum(ir1616) 2 + 0.68d0 * 0.5d0*y(io16) * s1 * ratdum(ir1616) dfdy(isi28,img24) = y(ihe4) * ratdum(irmgag) 1 + y(ihe4) * ratdum(irmgap) * (1.0d0-r1) dfdy(isi28,isi28) = -y(ihe4) * ratdum(irsiag) 1 - ratdum(irsiga) 2 - r1 * ratdum(irsigp) 3 - y(ihe4) * ratdum(irsiap) * (1.0d0-s1) dfdy(isi28,is32) = ratdum(irsga) 1 + s1 * ratdum(irsgp) c..s32 jacobian elements dfdy(is32,ihe4) = y(isi28) * ratdum(irsiag) 1 - y(is32) * ratdum(irsag) 2 + y(isi28) * ratdum(irsiap) * (1.0d0-s1) 3 - y(is32) * ratdum(irsap) * (1.0d0-t1) dfdy(is32,io16) = 0.68d0*0.5d0*y(io16)*ratdum(ir1616)*(1.0d0-s1) 1 + 0.2d0 * 0.5d0*y(io16) * ratdum(ir1616) dfdy(is32,isi28) = y(ihe4) * ratdum(irsiag) 1 + y(ihe4) * ratdum(irsiap) * (1.0d0-s1) dfdy(is32,is32) = -y(ihe4) * ratdum(irsag) 1 - ratdum(irsga) 2 - s1 * ratdum(irsgp) 3 - y(ihe4) * ratdum(irsap) * (1.0d0-t1) dfdy(is32,iar36) = ratdum(irarga) 1 + t1 * ratdum(irargp) c..ar36 jacobian elements dfdy(iar36,ihe4) = y(is32) * ratdum(irsag) 1 - y(iar36) * ratdum(irarag) 2 + y(is32) * ratdum(irsap) * (1.0d0-t1) 3 - y(iar36) * ratdum(irarap) * (1.0d0-u1) dfdy(iar36,is32) = y(ihe4) * ratdum(irsag) 1 + y(ihe4) * ratdum(irsap) * (1.0d0-t1) dfdy(iar36,iar36) = -y(ihe4) * ratdum(irarag) 1 - ratdum(irarga) 2 - t1 * ratdum(irargp) 3 - y(ihe4) * ratdum(irarap) * (1.0d0-u1) dfdy(iar36,ica40) = ratdum(ircaga) 1 + ratdum(ircagp) * u1 c..ca40 jacobian elements dfdy(ica40,ihe4) = y(iar36) * ratdum(irarag) 1 - y(ica40) * ratdum(ircaag) 2 + y(iar36) * ratdum(irarap)*(1.0d0-u1) 3 - y(ica40) * ratdum(ircaap)*(1.0d0-v1) dfdy(ica40,iar36) = y(ihe4) * ratdum(irarag) 1 + y(ihe4) * ratdum(irarap)*(1.0d0-u1) dfdy(ica40,ica40) = -y(ihe4) * ratdum(ircaag) 1 - ratdum(ircaga) 2 - ratdum(ircagp) * u1 3 - y(ihe4) * ratdum(ircaap)*(1.0d0-v1) dfdy(ica40,iti44) = ratdum(irtiga) 1 + ratdum(irtigp) * v1 c..ti44 jacobian elements dfdy(iti44,ihe4) = y(ica40) * ratdum(ircaag) 1 - y(iti44) * ratdum(irtiag) 2 + y(ica40) * ratdum(ircaap)*(1.0d0-v1) 3 - y(iti44) * ratdum(irtiap)*(1.0d0-w1) dfdy(iti44,ica40) = y(ihe4) * ratdum(ircaag) 1 + y(ihe4) * ratdum(ircaap)*(1.0d0-v1) dfdy(iti44,iti44) = -y(ihe4) * ratdum(irtiag) 1 - ratdum(irtiga) 2 - v1 * ratdum(irtigp) 3 - y(ihe4) * ratdum(irtiap)*(1.0d0-w1) dfdy(iti44,icr48) = ratdum(ircrga) 1 + w1 * ratdum(ircrgp) c..cr48 jacobian elements dfdy(icr48,ihe4) = y(iti44) * ratdum(irtiag) 1 - y(icr48) * ratdum(ircrag) 2 + y(iti44) * ratdum(irtiap)*(1.0d0-w1) 3 - y(icr48) * ratdum(ircrap)*(1.0d0-x1) dfdy(icr48,iti44) = y(ihe4) * ratdum(irtiag) 1 + y(ihe4) * ratdum(irtiap)*(1.0d0-w1) dfdy(icr48,icr48) = -y(ihe4) * ratdum(ircrag) 1 - ratdum(ircrga) 2 - w1 * ratdum(ircrgp) 3 - y(ihe4) * ratdum(ircrap)*(1.0d0-x1) dfdy(icr48,ife52) = ratdum(irfega) 1 + x1 * ratdum(irfegp) c..fe52 jacobian elements dfdy(ife52,ihe4) = y(icr48) * ratdum(ircrag) 1 - y(ife52) * ratdum(irfeag) 2 + y(icr48) * ratdum(ircrap) * (1.0d0-x1) 3 - y(ife52) * ratdum(irfeap) * (1.0d0-y1) dfdy(ife52,icr48) = y(ihe4) * ratdum(ircrag) 1 + y(ihe4) * ratdum(ircrap) * (1.0d0-x1) dfdy(ife52,ife52) = - y(ihe4) * ratdum(irfeag) 1 - ratdum(irfega) 2 - x1 * ratdum(irfegp) 3 - y(ihe4) * ratdum(irfeap) * (1.0d0-y1) dfdy(ife52,ini56) = ratdum(irniga) 1 + y1 * ratdum(irnigp) c..ni56 jacobian elements dfdy(ini56,ihe4) = y(ife52) * ratdum(irfeag) 1 + y(ife52) * ratdum(irfeap) * (1.0d0-y1) dfdy(ini56,ife52) = y(ihe4) * ratdum(irfeag) 1 + y(ihe4) * ratdum(irfeap) * (1.0d0-y1) dfdy(ini56,ini56) = -ratdum(irniga) 1 - y1 * ratdum(irnigp) c..append the temperature derivatives of the rate equations call rhs(y,dratdumdt,dydt) c..add in the parts from the dummy proton links dydt(ihe4) = dydt(ihe4) 1 + 0.34d0 *s1dt*0.5d0*y(io16)*y(io16)*ratdum(ir1616) 2 + y(ihe4) * y(img24) * ratdum(irmgap)* r1dt 3 + y(isi28) * ratdum(irsigp) * r1dt 4 + y(ihe4) * y(isi28) * ratdum(irsiap) * s1dt 5 + y(is32) * ratdum(irsgp) * s1dt 6 + y(ihe4) * y(is32) * ratdum(irsap) * t1dt 7 + y(iar36) * ratdum(irargp) * t1dt 8 + y(ihe4) * y(iar36) * ratdum(irarap) * u1dt 9 + y(ica40) * ratdum(ircagp) * u1dt & + y(ihe4) * y(ica40) * ratdum(ircaap) * v1dt 1 + y(iti44) * ratdum(irtigp) * v1dt 2 + y(ihe4) * y(iti44) * ratdum(irtiap) * w1dt 3 + y(icr48) * ratdum(ircrgp) * w1dt 4 + y(ihe4) * y(icr48) * ratdum(ircrap) * x1dt 5 + y(ife52) * ratdum(irfegp) * x1dt dydt(ihe4) = dydt(ihe4) 1 + y(ihe4) * y(ife52) * ratdum(irfeap) * y1dt 2 + y(ini56) * ratdum(irnigp) * y1dt dydt(img24) = dydt(img24) 1 + y(img24) * y(ihe4) *ratdum(irmgap)*r1dt 2 + y(isi28) * r1dt * ratdum(irsigp) dydt(isi28) = dydt(isi28) 1 + 0.34d0*0.5d0*y(io16)*y(io16)*s1dt*ratdum(ir1616) 2 - y(img24) * y(ihe4) * ratdum(irmgap)*r1dt 3 - y(isi28) * r1dt * ratdum(irsigp) 4 + y(isi28) * y(ihe4) * ratdum(irsiap)*s1dt 5 + y(is32) * s1dt * ratdum(irsgp) dydt(is32) = dydt(is32) 1 - 0.34d0*0.5d0*y(io16)**2 * ratdum(ir1616)*s1dt 2 - y(isi28) * y(ihe4) * ratdum(irsiap) * s1dt 3 - y(is32) * s1dt * ratdum(irsgp) 4 - y(is32) * y(ihe4) * ratdum(irsap) * t1dt 5 + y(iar36) * t1dt * ratdum(irargp) dydt(iar36) = dydt(iar36) 1 - y(is32) * y(ihe4) * ratdum(irsap) * t1dt 2 - y(iar36) * t1dt * ratdum(irargp) 3 + y(iar36) * y(ihe4) * ratdum(irarap) * u1dt 4 + y(ica40) * ratdum(ircagp) * u1dt dydt(ica40) = dydt(ica40) 1 - y(iar36) * y(ihe4) * ratdum(irarap) * u1dt 2 - y(ica40) * ratdum(ircagp) * u1dt 3 + y(ica40) * y(ihe4) * ratdum(ircaap) * v1dt 4 + y(iti44) * ratdum(irtigp) * v1dt dydt(iti44) = dydt(iti44) 1 - y(ica40) * y(ihe4) * ratdum(ircaap)*v1dt 2 - y(iti44) * v1dt * ratdum(irtigp) 3 + y(iti44) * y(ihe4) * ratdum(irtiap) * w1dt 4 + y(icr48) * w1dt * ratdum(ircrgp) dydt(icr48) = dydt(icr48) 1 - y(iti44) * y(ihe4) * ratdum(irtiap)*w1dt 2 - y(icr48) * w1dt * ratdum(ircrgp) 3 + y(icr48) * y(ihe4) * ratdum(ircrap) * x1dt 4 + y(ife52) * x1dt * ratdum(irfegp) dydt(ife52) = dydt(ife52) 1 - y(icr48) * y(ihe4) * ratdum(ircrap) * x1dt 2 - y(ife52) * x1dt * ratdum(irfegp) 3 + y(ife52) * y(ihe4) * ratdum(irfeap) * y1dt 4 + y(ini56) * y1dt * ratdum(irnigp) dydt(ini56) = dydt(ini56) 1 - y(ife52) * y(ihe4) * ratdum(irfeap) * y1dt 2 - y(ini56) * y1dt * ratdum(irnigp) do i=1,ionmax dfdy(i,itemp) = dydt(i) enddo c..append the density derivatives of the rate equations call rhs(y,dratdumdd,dydt) c..add in the parts from the dummy proton links dydt(ihe4) = dydt(ihe4) 1 + 0.34d0 *s1dd*0.5d0*y(io16)*y(io16)*ratdum(ir1616) 2 + y(ihe4) * y(img24) * ratdum(irmgap)* r1dd 3 + y(isi28) * ratdum(irsigp) * r1dd 4 + y(ihe4) * y(isi28) * ratdum(irsiap) * s1dd 5 + y(is32) * ratdum(irsgp) * s1dd 6 + y(ihe4) * y(is32) * ratdum(irsap) * t1dd 7 + y(iar36) * ratdum(irargp) * t1dd 8 + y(ihe4) * y(iar36) * ratdum(irarap) * u1dd 9 + y(ica40) * ratdum(ircagp) * u1dd & + y(ihe4) * y(ica40) * ratdum(ircaap) * v1dd 1 + y(iti44) * ratdum(irtigp) * v1dd 2 + y(ihe4) * y(iti44) * ratdum(irtiap) * w1dd 3 + y(icr48) * ratdum(ircrgp) * w1dd 4 + y(ihe4) * y(icr48) * ratdum(ircrap) * x1dd 5 + y(ife52) * ratdum(irfegp) * x1dd dydt(ihe4) = dydt(ihe4) 1 + y(ihe4) * y(ife52) * ratdum(irfeap) * y1dd 2 + y(ini56) * ratdum(irnigp) * y1dd dydt(img24) = dydt(img24) 1 + y(img24) * y(ihe4) *ratdum(irmgap)*r1dd 2 + y(isi28) * r1dd * ratdum(irsigp) dydt(isi28) = dydt(isi28) 1 + 0.34d0*0.5d0*y(io16)*y(io16)*s1dd*ratdum(ir1616) 2 - y(img24) * y(ihe4) * ratdum(irmgap)*r1dd 3 - y(isi28) * r1dd * ratdum(irsigp) 4 + y(isi28) * y(ihe4) * ratdum(irsiap)*s1dd 5 + y(is32) * s1dd * ratdum(irsgp) dydt(is32) = dydt(is32) 1 - 0.34d0*0.5d0*y(io16)**2 * ratdum(ir1616)*s1dd 2 - y(isi28) * y(ihe4) * ratdum(irsiap) * s1dd 3 - y(is32) * s1dd * ratdum(irsgp) 4 - y(is32) * y(ihe4) * ratdum(irsap) * t1dd 5 + y(iar36) * t1dd * ratdum(irargp) dydt(iar36) = dydt(iar36) 1 - y(is32) * y(ihe4) * ratdum(irsap) * t1dd 2 - y(iar36) * t1dd * ratdum(irargp) 3 + y(iar36) * y(ihe4) * ratdum(irarap) * u1dd 4 + y(ica40) * ratdum(ircagp) * u1dd dydt(ica40) = dydt(ica40) 1 - y(iar36) * y(ihe4) * ratdum(irarap) * u1dd 2 - y(ica40) * ratdum(ircagp) * u1dd 3 + y(ica40) * y(ihe4) * ratdum(ircaap) * v1dd 4 + y(iti44) * ratdum(irtigp) * v1dd dydt(iti44) = dydt(iti44) 1 - y(ica40) * y(ihe4) * ratdum(ircaap)*v1dd 2 - y(iti44) * v1dd * ratdum(irtigp) 3 + y(iti44) * y(ihe4) * ratdum(irtiap) * w1dd 4 + y(icr48) * w1dd * ratdum(ircrgp) dydt(icr48) = dydt(icr48) 1 - y(iti44) * y(ihe4) * ratdum(irtiap)*w1dd 2 - y(icr48) * w1dd * ratdum(ircrgp) 3 + y(icr48) * y(ihe4) * ratdum(ircrap) * x1dd 4 + y(ife52) * x1dd * ratdum(irfegp) dydt(ife52) = dydt(ife52) 1 - y(icr48) * y(ihe4) * ratdum(ircrap) * x1dd 2 - y(ife52) * x1dd * ratdum(irfegp) 3 + y(ife52) * y(ihe4) * ratdum(irfeap) * y1dd 4 + y(ini56) * y1dd * ratdum(irnigp) dydt(ini56) = dydt(ini56) 1 - y(ife52) * y(ihe4) * ratdum(irfeap) * y1dd 2 - y(ini56) * y1dd * ratdum(irnigp) do i=1,ionmax dfdy(i,iden) = dydt(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 c..account for the neutrino losses call sneut5(btemp,bden,abar,zbar, 1 sneut,snudt,snudd,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) - snudt dfdy(iener,iden) = dfdy(iener,iden) - snudd c..for hydrostatic or one step burns all the temperature and density c..jacobian elements are zero, 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 else if (self_heat) 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) suma = 0.0d0 do i=1,ionmax suma = suma - dfdy(i,itemp) enddo sumz = 0.0d0 do i=1,ionmax sumz = sumz + (zion(i) - zbar)*dfdy(i,itemp) enddo ww = suma*dea_row(1)*abar*abar + sumz*dez_row(1)*abar dfdy(itemp,itemp) = zz*(dfdy(iener,itemp) - ww) c..d(itemp)/d(den) suma = 0.0d0 do i=1,ionmax suma = suma - dfdy(i,iden) enddo sumz = 0.0d0 do i=1,ionmax sumz = sumz + (zion(i) - zbar)*dfdy(i,iden) enddo ww = suma*dea_row(1)*abar*abar + sumz*dez_row(1)*abar dfdy(itemp,iden) = zz*(dfdy(iener,iden) - ww) c..for detonations else if (detonation) then c..get the right hand sides call rhs(y,ratdum,dydt) c..instantaneous energy generation rate enuc = 0.0d0 do i=1,ionmax enuc = enuc + dydt(i) * bion(i) enddo enuc = enuc * conv dydt(iener) = enuc - sneut c..map the rest of the input vector velx = y(ivelx) posx = y(iposx) c..it appears as if we need the derivatives of derivatibe based eos quantities. grrr. z = bden xx = 0.01d0*z bden = z + xx temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos foo(1) = cs_row(1) foo(2) = dpt_row(1) foo(3) = dpt_row(1)/det_row(1) foo(4) = dpd_row(1) foo(5) = dpa_row(1) foo(6) = dpz_row(1) foo(7) = dea_row(1) foo(8) = dez_row(1) bden = z - xx temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos moo(1) = cs_row(1) moo(2) = dpt_row(1) moo(3) = dpt_row(1)/det_row(1) moo(4) = dpd_row(1) moo(5) = dpa_row(1) moo(6) = dpz_row(1) moo(7) = dea_row(1) moo(8) = dez_row(1) bden = z z = 0.5d0/xx csbd = (foo(1) - moo(1))*z dptbd = (foo(2) - moo(2))*z dpdebd = (foo(3) - moo(3))*z dpdbd = (foo(4) - moo(4))*z dpabd = (foo(5) - moo(5))*z dpzbd = (foo(6) - moo(6))*z deabd = (foo(7) - moo(7))*z dezbd = (foo(8) - moo(8))*z z = btemp xx = 0.01d0*z btemp = z + xx temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos foo(1) = cs_row(1) foo(2) = dpt_row(1) foo(3) = dpt_row(1)/det_row(1) foo(4) = dpd_row(1) foo(5) = dpa_row(1) foo(6) = dpz_row(1) foo(7) = dea_row(1) foo(8) = dez_row(1) btemp = z - xx temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos moo(1) = cs_row(1) moo(2) = dpt_row(1) moo(3) = dpt_row(1)/det_row(1) moo(4) = dpd_row(1) moo(5) = dpa_row(1) moo(6) = dpz_row(1) moo(7) = dea_row(1) moo(8) = dez_row(1) btemp = z z = 0.5d0/xx csbt = (foo(1) - moo(1))*z dptbt = (foo(2) - moo(2))*z dpdebt = (foo(3) - moo(3))*z dpdbt = (foo(4) - moo(4))*z dpabt = (foo(5) - moo(5))*z dpzbt = (foo(6) - moo(6))*z deabt = (foo(7) - moo(7))*z dezbt = (foo(8) - moo(8))*z 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..for de/dy and dp/dy suma = 0.0d0 do i=1,ionmax suma = suma - dydt(i) enddo sumz = 0.0d0 do i=1,ionmax sumz = sumz + (zion(i) - zbar)*dydt(i) enddo sum = 0.0d0 do i=1,ionmax sum = sum - dfdy(i,iden) enddo sum5 = sum * dea_row(1)*abar*abar sum9 = sum * dpa_row(1)*abar*abar sum = 0.0d0 do i=1,ionmax sum = sum + (zion(i) - zbar)*dfdy(i,iden) enddo sum6 = sum * dez_row(1)*abar sum10 = sum * dpz_row(1)*abar sum = 0.0d0 do i=1,ionmax sum = sum - dfdy(i,itemp) enddo sum7 = sum * dea_row(1)*abar*abar sum11 = sum * dpa_row(1)*abar*abar sum = 0.0d0 do i=1,ionmax sum = sum + (zion(i) - zbar)*dfdy(i,itemp) enddo sum8 = sum * dez_row(1)*abar sum12 = sum * dpz_row(1)*abar c..the possibly singular denominator cs = cs_row(1) denom = velx*velx - cs*cs denombv = 2.0d0*velx denomdd = -2.0d0*cs*csbd denomdt = -2.0d0*cs*csbt c..the function phi dpde = dpt_row(1)/det_row(1) z = suma*dpa_row(1)*abar*abar + sumz*dpz_row(1)*abar zbd = suma*dpabd*abar*abar + sumz*dpzbd*abar + sum9 + sum10 zbt = suma*dpabt*abar*abar + sumz*dpzbt*abar + sum11 + sum12 ww = suma*dea_row(1)*abar*abar + sumz*dez_row(1)*abar wwbd = suma*deabd*abar*abar + sumz*dezbd*abar + sum5 + sum6 wwbt = suma*deabt*abar*abar + sumz*dezbt*abar + sum7 + sum8 phi = dpde*(dydt(iener) - ww) - z phibd = dpdebd*(dydt(iener) - ww) - zbd 1 + dpde*(dfdy(iener,iden) - wwbd) phibt = dpdebt*(dydt(iener) - ww) - zbt 1 + dpde*(dfdy(iener,itemp) - wwbt) c..a common combination if (denom .ne. 0.0) then combo = phi/denom combobv = -combo/denom*denombv combobd = -combo/denom*denomdd + phibd/denom combobt = -combo/denom*denomdt + phibt/denom else combo = 0.0d0 combobv = 0.0d0 combobd = 0.0d0 combobt = 0.0d0 end if c..position equation dydt(iposx) = velx dfdy(iposx,ivelx) = 1.0d0 c..density equation dydt(iden) = combo dfdy(iden,ivelx) = combobv dfdy(iden,iden) = combobd dfdy(iden,itemp) = combobt c..d(iden)/d(yi) yy = 1.0d0/denom zz = dpde*yy do j=1,ionmax xx = 0.0d0 ww = 0.0d0 do i=1,ionmax xx = xx + dfdy(i,j)*(-dea_row(1)*abar*abar 1 + dez_row(1)*(zion(i)-zbar)*abar) ww = ww + dfdy(i,j)*(-dpa_row(1)*abar*abar 1 + dpz_row(1)*(zion(i)-zbar)*abar) enddo dfdy(iden,j) = zz*(dfdy(iener,j) - xx) - ww*yy enddo c..velocity equation z = velx/bden dydt(ivelx) = -z*dydt(iden) dfdy(ivelx,ivelx) = -dydt(iden)/bden - z*dfdy(iden,ivelx) dfdy(ivelx,iden) = z/bden*dydt(iden) - z*dfdy(iden,iden) dfdy(ivelx,itemp) = -z*dfdy(iden,itemp) c..d(ivelx)/d(yi) do j=1,ionmax dfdy(ivelx,j) = -z*dfdy(iden,j) enddo c..temperature equation dtdp = 1.0d0/dpt_row(1) dtdpbd = -dtdp*dtdp*dptbd dtdpbt = -dtdp*dtdp*dptbt ww = suma*dpa_row(1)*abar*abar + sumz*dpz_row(1)*abar wwbd = suma*dpabd*abar*abar + sumz*dpzbd*abar + sum9 + sum10 wwbt = suma*dpabt*abar*abar + sumz*dpzbt*abar + sum11 + sum12 dydt(itemp) = dtdp*((velx*velx - dpd_row(1))*dydt(iden) - ww) dfdy(itemp,ivelx) = dtdp*(2.0d0*velx*dydt(iden) 1 + (velx*velx - dpd_row(1))*dfdy(iden,ivelx)) dfdy(itemp,iden) = dtdpbd*((velx*velx-dpd_row(1))*dydt(iden)-ww) 1 + dtdp*((velx*velx-dpd_row(1))*dfdy(iden,iden) 2 - dpdbd*dydt(iden) - wwbd) dfdy(itemp,itemp) = dtdpbt*((velx*velx-dpd_row(1))*dydt(iden)-ww) 1 + dtdp*((velx*velx-dpd_row(1))*dfdy(iden,itemp) 2 - dpdbt*dydt(iden) - wwbt) c..d(itemp)/d(yi) do j=1,ionmax ww = 0.0d0 c do i=1,ionmax c ww = ww + dfdy(i,j)*(-dpa_row(1)*abar*abar c 1 + dpz_row(1)*(zion(i)-zbar)*abar) c enddo dfdy(itemp,j) = dtdp*((velx*velx - dpd_row(1))*dfdy(iden,j) -ww) enddo c..end of burning mode ifs end if return end subroutine baprox13(iloc,jloc,nzo,np) include 'implno.dek' include 'network.dek' c.. c..this routine builds the nonzero matrix locations for saprox13 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 integer iloc(1),jloc(1),nzo,np,i c..communicate with saprox13 integer neloc parameter (neloc=172) integer eloc(neloc),nterms common /elca13/ eloc,nterms c..initialize nterms = 0 nzo = 0 do i=1,neloc eloc(i) = 0 enddo c..tag the nonzero locations c..he4 jacobian elements call mcord(ihe4,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,ic12,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,io16,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,ine20,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,img24,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,isi28,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,is32,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,iar36,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,ica40,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,iti44,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,icr48,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,ife52,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ihe4,ini56,iloc,jloc,nzo,np,eloc,nterms,neloc) c..c12 jacobian elements call mcord(ic12,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ic12,ic12,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ic12,io16,iloc,jloc,nzo,np,eloc,nterms,neloc) c..o16 jacobian elements call mcord(io16,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(io16,ic12,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(io16,io16,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(io16,ine20,iloc,jloc,nzo,np,eloc,nterms,neloc) c..ne20 jacobian elements call mcord(ine20,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ine20,ic12,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ine20,io16,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ine20,ine20,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ine20,img24,iloc,jloc,nzo,np,eloc,nterms,neloc) c..mg24 jacobian elements call mcord(img24,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(img24,ic12,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(img24,io16,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(img24,ine20,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(img24,img24,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(img24,isi28,iloc,jloc,nzo,np,eloc,nterms,neloc) c..si28 jacobian elements call mcord(isi28,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(isi28,ic12,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(isi28,io16,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(isi28,img24,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(isi28,isi28,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(isi28,is32,iloc,jloc,nzo,np,eloc,nterms,neloc) c..s32 jacobian elements call mcord(is32,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(is32,io16,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(is32,isi28,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(is32,is32,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(is32,iar36,iloc,jloc,nzo,np,eloc,nterms,neloc) c..ar36 jacobian elements call mcord(iar36,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iar36,is32,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iar36,iar36,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iar36,ica40,iloc,jloc,nzo,np,eloc,nterms,neloc) c..ca40 jacobian elements call mcord(ica40,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ica40,iar36,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ica40,ica40,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ica40,iti44,iloc,jloc,nzo,np,eloc,nterms,neloc) c..ti44 jacobian elements call mcord(iti44,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iti44,ica40,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iti44,iti44,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iti44,icr48,iloc,jloc,nzo,np,eloc,nterms,neloc) c..cr48 jacobian elements call mcord(icr48,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(icr48,iti44,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(icr48,icr48,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(icr48,ife52,iloc,jloc,nzo,np,eloc,nterms,neloc) c..fe52 jacobian elements call mcord(ife52,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ife52,icr48,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ife52,ife52,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ife52,ini56,iloc,jloc,nzo,np,eloc,nterms,neloc) c..ni56 jacobian elements call mcord(ini56,ihe4,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ini56,ife52,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ini56,ini56,iloc,jloc,nzo,np,eloc,nterms,neloc) c..temperature contributions do i=1,ionmax call mcord(i,itemp,iloc,jloc,nzo,np,eloc,nterms,neloc) end do c..density contributions do i=1,ionmax call mcord(i,iden,iloc,jloc,nzo,np,eloc,nterms,neloc) end do c..energy equation jacobian elements do i=1,ionmax call mcord(iener,i,iloc,jloc,nzo,np,eloc,nterms,neloc) enddo call mcord(iener,iener,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iener,itemp,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iener,iden,iloc,jloc,nzo,np,eloc,nterms,neloc) c..neutrino losses do i=1,ionmax call mcord(iener,i,iloc,jloc,nzo,np,eloc,nterms,neloc) enddo call mcord(iener,itemp,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iener,iden,iloc,jloc,nzo,np,eloc,nterms,neloc) c..the jacobian elements depend on the burning mode c..hydrostatic or single step if (hydrostatic .or. one_step) then call mcord(itemp,itemp,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iden,iden,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ivelx,ivelx,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iposx,iposx,iloc,jloc,nzo,np,eloc,nterms,neloc) c..adiabatic expansion else if (expansion) then call mcord(itemp,itemp,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iden,iden,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ivelx,ivelx,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iposx,iposx,iloc,jloc,nzo,np,eloc,nterms,neloc) c..self heating else if (self_heat) then do i=1,ionmax call mcord(itemp,i,iloc,jloc,nzo,np,eloc,nterms,neloc) enddo call mcord(itemp,itemp,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(itemp,iden,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iden,iden,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ivelx,ivelx,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iposx,iposx,iloc,jloc,nzo,np,eloc,nterms,neloc) c..detonation else if (detonation) then call mcord(iposx,iposx,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iposx,ivelx,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iden,ivelx,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iden,iden,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(iden,itemp,iloc,jloc,nzo,np,eloc,nterms,neloc) do i=1,ionmax call mcord(iden,i,iloc,jloc,nzo,np,eloc,nterms,neloc) enddo call mcord(ivelx,ivelx,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ivelx,iden,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(ivelx,itemp,iloc,jloc,nzo,np,eloc,nterms,neloc) do i=1,ionmax call mcord(ivelx,i,iloc,jloc,nzo,np,eloc,nterms,neloc) enddo call mcord(itemp,ivelx,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(itemp,iden,iloc,jloc,nzo,np,eloc,nterms,neloc) call mcord(itemp,itemp,iloc,jloc,nzo,np,eloc,nterms,neloc) do i=1,ionmax call mcord(itemp,i,iloc,jloc,nzo,np,eloc,nterms,neloc) enddo end if c..write a diagnostic c write(6,*) ' ' c write(6,*) nzo,' matrix elements' c write(6,*) nterms,' jacobian contributions' c write(6,*) ' ' c read(5,*) return end subroutine saprox13(tt,y,dfdy,nzo) include 'implno.dek' include 'const.dek' include 'burn_common.dek' include 'network.dek' include 'vector_eos.dek' include 'cjdet.dek' c..this routine sets up the sparse aprox13 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 denom,denomdt,denomdd,a1,a2,a3,a4, 1 r1,r1dt,r1dd,s1,s1dt,s1dd,t1,t1dt,t1dd, 2 u1,u1dt,u1dd,v1,v1dt,v1dd,w1,w1dt,w1dd, 3 x1,x1dt,x1dd,y1,y1dt,y1dd,zz,xx, 4 teqn(ionmax),deqn(ionmax), 5 xsum(ionmax),ysum(ionmax),zsum(ionmax) double precision zbarxx,ytot1,abar,zbar,taud,taut, 1 snudt,snudd,snuda,snudz, 2 dydt(neqs),enuc,velx,posx,suma,sumz,sum1,sum, 3 sum5,sum6,sum7,sum8,sum9,sum10,sum11,sum12, 4 cs,denombv,dpde,dpdebd,dpdebt,phi,phibd,phibt, 5 combo,combobv,combobd,combobt,z,zbd,zbt, 6 ww,wwbd,wwbt,dtdp,dtdpbd,dtdpbt, 7 dpdbd,dpdbt,foo(8),moo(8),enucbd,enucbt, 8 csbd,csbt,dptbd,dptbt,detbd,detbt,dpabd,dpabt, 9 dpzbd,dpzbt,deabd,deabt,dezbd,dezbt double precision conv parameter (conv = ev2erg*1.0d6*avo) c..communicate with the jacobian builder integer neloc parameter (neloc=172) integer eloc(neloc),nterms common /elca13/ 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..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..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 c..get the reaction rates if (use_tables .eq. 1) then call aprox13tab else call aprox13rat end if c..do the screening here because the corrections depend on the composition call screen_aprox13(y) c..branching ratios for various dummy proton links r1 = 0.0d0 r1dt = 0.0d0 r1dd = 0.0d0 denom = ratdum(iralpa) + ratdum(iralpg) denomdt = dratdumdt(iralpa) + dratdumdt(iralpg) denomdd = dratdumdd(iralpa) + dratdumdd(iralpg) if (denom .ne. 0.0) then zz = 1.0d0/denom r1 = ratdum(iralpa)*zz r1dt = (dratdumdt(iralpa) - r1*denomdt)*zz r1dd = (dratdumdd(iralpa) - r1*denomdd)*zz end if s1 = 0.0d0 s1dt = 0.0d0 s1dd = 0.0d0 denom = ratdum(irppa) + ratdum(irppg) denomdt = dratdumdt(irppa) + dratdumdt(irppg) denomdd = dratdumdd(irppa) + dratdumdd(irppg) if (denom .ne. 0.0) then zz = 1.0d0/denom s1 = ratdum(irppa)*zz s1dt = (dratdumdt(irppa) - s1*denomdt)*zz s1dd = (dratdumdd(irppa) - s1*denomdd)*zz end if t1 = 0.0d0 t1dt = 0.0d0 t1dd = 0.0d0 denom = ratdum(irclpa) + ratdum(irclpg) denomdt = dratdumdt(irclpa) + dratdumdt(irclpg) denomdd = dratdumdd(irclpa) + dratdumdd(irclpg) if (denom .ne. 0.0) then zz = 1.0d0/denom t1 = ratdum(irclpa)*zz t1dt = (dratdumdt(irclpa) - t1*denomdt)*zz t1dd = (dratdumdd(irclpa) - t1*denomdd)*zz end if u1 = 0.0d0 u1dt = 0.0d0 u1dd = 0.0d0 denom = ratdum(irkpa) + ratdum(irkpg) denomdt = dratdumdt(irkpa) + dratdumdt(irkpg) denomdd = dratdumdd(irkpa) + dratdumdd(irkpg) if (denom .ne. 0.0) then zz = 1.0d0/denom u1 = ratdum(irkpa)*zz u1dt = (dratdumdt(irkpa) - u1*denomdt)*zz u1dd = (dratdumdd(irkpa) - u1*denomdd)*zz end if v1 = 0.0d0 v1dt = 0.0d0 v1dd = 0.0d0 denom = ratdum(irscpa) + ratdum(irscpg) denomdt = dratdumdt(irscpa) + dratdumdt(irscpg) denomdd = dratdumdd(irscpa) + dratdumdd(irscpg) if (denom .ne. 0.0) then zz = 1.0d0/denom v1 = ratdum(irscpa)*zz v1dt = (dratdumdt(irscpa) - v1*denomdt)*zz v1dd = (dratdumdd(irscpa) - v1*denomdd)*zz end if w1 = 0.0d0 w1dt = 0.0d0 w1dd = 0.0d0 denom = ratdum(irvpa) + ratdum(irvpg) denomdt = dratdumdt(irvpa) + dratdumdt(irvpg) denomdd = dratdumdd(irvpa) + dratdumdd(irvpg) if (denom .ne. 0.0) then zz = 1.0d0/denom w1 = ratdum(irvpa)*zz w1dt = (dratdumdt(irvpa) - w1*denomdt)*zz w1dd = (dratdumdd(irvpa) - w1*denomdd)*zz end if x1 = 0.0d0 x1dt = 0.0d0 x1dd = 0.0d0 denom = ratdum(irmnpa) + ratdum(irmnpg) denomdt = dratdumdt(irmnpa) + dratdumdt(irmnpg) denomdd = dratdumdd(irmnpa) + dratdumdd(irmnpg) if (denom .ne. 0.0) then zz = 1.0d0/denom x1 = ratdum(irmnpa)*zz x1dt = (dratdumdt(irmnpa) - x1*denomdt)*zz x1dd = (dratdumdd(irmnpa) - x1*denomdd)*zz endif y1 = 0.0d0 y1dt = 0.0d0 y1dd = 0.0d0 denom = ratdum(ircopa) + ratdum(ircopg) denomdt = dratdumdt(ircopa) + dratdumdt(ircopg) denomdd = dratdumdd(ircopa) + dratdumdd(ircopg) if (denom .ne. 0.0) then zz = 1.0d0/denom y1 = ratdum(ircopa)*zz y1dt = (dratdumdt(ircopa) - y1*denomdt)*zz y1dd = (dratdumdd(ircopa) - y1*denomdd)*zz end if c..he4 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 - y(isi28) * ratdum(irsiag) 6 - y(is32) * ratdum(irsag) 7 - y(iar36) * ratdum(irarag) 8 - y(ica40) * ratdum(ircaag) 9 - y(iti44) * ratdum(irtiag) & - y(icr48) * ratdum(ircrag) 1 - y(ife52) * ratdum(irfeag) a1 = a1 1 - y(img24) * ratdum(irmgap) * (1.0d0-r1) 2 - y(isi28) * ratdum(irsiap) * (1.0d0-s1) 3 - y(is32) * ratdum(irsap) * (1.0d0-t1) 4 - y(iar36) * ratdum(irarap) * (1.0d0-u1) 5 - y(ica40) * ratdum(ircaap) * (1.0d0-v1) 6 - y(iti44) * ratdum(irtiap) * (1.0d0-w1) 7 - y(icr48) * ratdum(ircrap) * (1.0d0-x1) 8 - y(ife52) * ratdum(irfeap) * (1.0d0-y1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ihe4) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(c12) a1 = y(ic12) * ratdum(ir1212) 1 + 0.5d0 * y(io16) * ratdum(ir1216) 2 + 3.0d0 * ratdum(irg3a) 3 - y(ihe4) * ratdum(ircag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(ihe4) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(o16) a1 = 0.5d0 * y(ic12) * ratdum(ir1216) 1 + 1.12d0 * 0.5d0*y(io16) * ratdum(ir1616) 2 + 0.68d0 * s1 * 0.5d0*y(io16) * ratdum(ir1616) 3 + ratdum(iroga) 4 - 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) 2 - y(ihe4) * ratdum(irmgap) * (1.0d0-r1) 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 - y(ihe4) * ratdum(irsiag) 2 - y(ihe4) * ratdum(irsiap) * (1.0d0-s1) 3 + r1 * ratdum(irsigp) 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(s32) a1 = ratdum(irsga) 1 - y(ihe4) * ratdum(irsag) 2 - y(ihe4) * ratdum(irsap) * (1.0d0-t1) 3 + s1 * ratdum(irsgp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(is32) = xsum(is32) + a1 * bion(ihe4) ysum(is32) = ysum(is32) - a1 zsum(is32) = zsum(is32) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(ar36) a1 = ratdum(irarga) 1 - y(ihe4) * ratdum(irarag) 2 - y(ihe4) * ratdum(irarap) * (1.0d0-u1) 3 + t1 * ratdum(irargp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(iar36) = xsum(iar36) + a1 * bion(ihe4) ysum(iar36) = ysum(iar36) - a1 zsum(iar36) = zsum(iar36) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(ca40) a1 = ratdum(ircaga) 1 - y(ihe4) * ratdum(ircaag) 2 - y(ihe4) * ratdum(ircaap) * (1.0d0-v1) 3 + u1 * ratdum(ircagp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ica40) = xsum(ica40) + a1 * bion(ihe4) ysum(ica40) = ysum(ica40) - a1 zsum(ica40) = zsum(ica40) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(ti44) a1 = ratdum(irtiga) 1 - y(ihe4) * ratdum(irtiag) 2 - y(ihe4) * ratdum(irtiap) * (1.0d0-w1) 3 + v1 * ratdum(irtigp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(iti44) = xsum(iti44) + a1 * bion(ihe4) ysum(iti44) = ysum(iti44) - a1 zsum(iti44) = zsum(iti44) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(cr48) a1 = ratdum(ircrga) 1 - y(ihe4) * ratdum(ircrag) 2 - y(ihe4) * ratdum(ircrap) * (1.0d0-x1) 3 + w1 * ratdum(ircrgp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(icr48) = xsum(icr48) + a1 * bion(ihe4) ysum(icr48) = ysum(icr48) - a1 zsum(icr48) = zsum(icr48) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(fe52) a1 = ratdum(irfega) 1 - y(ihe4) * ratdum(irfeag) 2 - y(ihe4) * ratdum(irfeap) * (1.0d0-y1) 3 + x1 * ratdum(irfegp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ife52) = xsum(ife52) + a1 * bion(ihe4) ysum(ife52) = ysum(ife52) - a1 zsum(ife52) = zsum(ife52) + a1 * (zion(ihe4) - zbar) c..d(he4)/d(ni56) a1 = ratdum(irniga) 1 + y1 * ratdum(irnigp) 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..c12 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 = -2.0d0 * y(ic12) * ratdum(ir1212) 1 - y(io16) * ratdum(ir1216) 2 - ratdum(irg3a) 3 - y(ihe4) * ratdum(ircag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(ic12) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(ic12) - zbar) c..d(c12)/d(o16) a1 = -y(ic12) * ratdum(ir1216) 1 + ratdum(iroga) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(ic12) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(ic12) - zbar) c..o16 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(io16)*ratdum(ir1216) 1 + y(ihe4)*ratdum(ircag) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ic12) = xsum(ic12) + a1 * bion(io16) ysum(ic12) = ysum(ic12) - a1 zsum(ic12) = zsum(ic12) + a1 * (zion(io16) - zbar) c..d(o16)/d(o16) a1 = - y(ic12) * ratdum(ir1216) 1 - 2.0d0 * y(io16) * ratdum(ir1616) 2 - y(ihe4) * ratdum(iroag) 3 - ratdum(iroga) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(io16) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(io16) - zbar) c..d(o16)/d(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..ne20 jacobian elements c..d(ne20)/d(he4) a1 = y(io16) * ratdum(iroag) 1 - 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 = -y(ihe4) * ratdum(irneag) 1 - ratdum(irnega) 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..mg24 jacobian elements c..d(mg24)/d(he4) a1 = y(ine20) * ratdum(irneag) 1 -y(img24) * ratdum(irmgag) 2 -y(img24) * ratdum(irmgap) * (1.0d0-r1) 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 = -y(ihe4) * ratdum(irmgag) 1 - ratdum(irmgga) 2 - y(ihe4) * ratdum(irmgap) * (1.0d0-r1) 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) 1 + r1 * ratdum(irsigp) 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..si28 jacobian elements c..d(si28)/d(he4) a1 = y(img24) * ratdum(irmgag) 1 - y(isi28) * ratdum(irsiag) 2 + y(img24) * ratdum(irmgap) * (1.0d0-r1) 3 - y(isi28) * ratdum(irsiap) * (1.0d0-s1) 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 = 0.5d0 * y(ic12) * ratdum(ir1216) 1 + 1.12d0 * 0.5d0*y(io16) * ratdum(ir1616) 2 + 0.68d0 * 0.5d0*y(io16) * s1 * ratdum(ir1616) 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) 1 + y(ihe4) * ratdum(irmgap) * (1.0d0-r1) 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 = -y(ihe4) * ratdum(irsiag) 1 - ratdum(irsiga) 2 - r1 * ratdum(irsigp) 3 - y(ihe4) * ratdum(irsiap) * (1.0d0-s1) 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(s32) a1 = ratdum(irsga) 1 + s1 * ratdum(irsgp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(is32) = xsum(is32) + a1 * bion(isi28) ysum(is32) = ysum(is32) - a1 zsum(is32) = zsum(is32) + a1 * (zion(isi28) - zbar) c..s32 jacobian elements c..d(s32)/d(he4) a1 = y(isi28) * ratdum(irsiag) 1 - y(is32) * ratdum(irsag) 2 + y(isi28) * ratdum(irsiap) * (1.0d0-s1) 3 - y(is32) * ratdum(irsap) * (1.0d0-t1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(is32) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(is32) - zbar) c..d(s32)/d(o16) a1 = 0.68d0 * 0.5d0*y(io16) * ratdum(ir1616) * (1.0d0-s1) 1 + 0.2d0 * 0.5d0*y(io16) * ratdum(ir1616) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(io16) = xsum(io16) + a1 * bion(is32) ysum(io16) = ysum(io16) - a1 zsum(io16) = zsum(io16) + a1 * (zion(is32) - zbar) c..d(s32)/d(si28) a1 = y(ihe4) * ratdum(irsiag) 1 + y(ihe4) * ratdum(irsiap) * (1.0d0-s1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(isi28) = xsum(isi28) + a1 * bion(is32) ysum(isi28) = ysum(isi28) - a1 zsum(isi28) = zsum(isi28) + a1 * (zion(is32) - zbar) c..d(s32)/d(s32) a1 = -y(ihe4) * ratdum(irsag) 1 - ratdum(irsga) 2 - s1 * ratdum(irsgp) 3 - y(ihe4) * ratdum(irsap) * (1.0d0-t1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(is32) = xsum(is32) + a1 * bion(is32) ysum(is32) = ysum(is32) - a1 zsum(is32) = zsum(is32) + a1 * (zion(is32) - zbar) c..d(s32)/d(ar36) a1 = ratdum(irarga) 1 + t1 * ratdum(irargp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(iar36) = xsum(iar36) + a1 * bion(is32) ysum(iar36) = ysum(iar36) - a1 zsum(iar36) = zsum(iar36) + a1 * (zion(is32) - zbar) c..ar36 jacobian elements c..d(ar36)/d(he4) a1 = y(is32) * ratdum(irsag) 1 - y(iar36) * ratdum(irarag) 2 + y(is32) * ratdum(irsap) * (1.0d0-t1) 3 - y(iar36) * ratdum(irarap) * (1.0d0-u1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(iar36) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(iar36) - zbar) c..d(ar36)/d(s32) a1 = y(ihe4) * ratdum(irsag) 1 + y(ihe4) * ratdum(irsap) * (1.0d0-t1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(is32) = xsum(is32) + a1 * bion(iar36) ysum(is32) = ysum(is32) - a1 zsum(is32) = zsum(is32) + a1 * (zion(iar36) - zbar) c..d(ar36)/d(ar36) a1 = -y(ihe4) * ratdum(irarag) 1 - ratdum(irarga) 2 - t1 * ratdum(irargp) 3 - y(ihe4) * ratdum(irarap) * (1.0d0-u1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(iar36) = xsum(iar36) + a1 * bion(iar36) ysum(iar36) = ysum(iar36) - a1 zsum(iar36) = zsum(iar36) + a1 * (zion(iar36) - zbar) c..d(ar36)/d(ca40) a1 = ratdum(ircaga) 1 + ratdum(ircagp) * u1 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ica40) = xsum(ica40) + a1 * bion(iar36) ysum(ica40) = ysum(ica40) - a1 zsum(ica40) = zsum(ica40) + a1 * (zion(iar36) - zbar) c..ca40 jacobian elements c..d(ca40)/d(he4) a1 = y(iar36) * ratdum(irarag) 1 - y(ica40) * ratdum(ircaag) 2 + y(iar36) * ratdum(irarap)*(1.0d0-u1) 3 - y(ica40) * ratdum(ircaap)*(1.0d0-v1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ica40) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ica40) - zbar) c..d(ca40)/d(ar36) a1 = y(ihe4) * ratdum(irarag) 1 + y(ihe4) * ratdum(irarap)*(1.0d0-u1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(iar36) = xsum(iar36) + a1 * bion(ica40) ysum(iar36) = ysum(iar36) - a1 zsum(iar36) = zsum(iar36) + a1 * (zion(ica40) - zbar) c..d(ca40)/d(ca40) a1 = -y(ihe4) * ratdum(ircaag) 1 - ratdum(ircaga) 2 - ratdum(ircagp) * u1 3 - y(ihe4) * ratdum(ircaap)*(1.0d0-v1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ica40) = xsum(ica40) + a1 * bion(ica40) ysum(ica40) = ysum(ica40) - a1 zsum(ica40) = zsum(ica40) + a1 * (zion(ica40) - zbar) c..d(ca40)/d(ti44) a1 = ratdum(irtiga) 1 + ratdum(irtigp) * v1 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(iti44) = xsum(iti44) + a1 * bion(ica40) ysum(iti44) = ysum(iti44) - a1 zsum(iti44) = zsum(iti44) + a1 * (zion(ica40) - zbar) c..ti44 jacobian elements c..d(ti44)/d(he4) a1 = y(ica40) * ratdum(ircaag) 1 - y(iti44) * ratdum(irtiag) 2 + y(ica40) * ratdum(ircaap)*(1.0d0-v1) 3 - y(iti44) * ratdum(irtiap)*(1.0d0-w1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(iti44) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(iti44) - zbar) c..d(ti44)/d(ca40) a1 = y(ihe4) * ratdum(ircaag) 1 + y(ihe4) * ratdum(ircaap)*(1.0d0-v1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ica40) = xsum(ica40) + a1 * bion(iti44) ysum(ica40) = ysum(ica40) - a1 zsum(ica40) = zsum(ica40) + a1 * (zion(iti44) - zbar) c..d(ti44)/d(ti44) a1 = -y(ihe4) * ratdum(irtiag) 1 - ratdum(irtiga) 2 - v1 * ratdum(irtigp) 3 - y(ihe4) * ratdum(irtiap)*(1.0d0-w1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(iti44) = xsum(iti44) + a1 * bion(iti44) ysum(iti44) = ysum(iti44) - a1 zsum(iti44) = zsum(iti44) + a1 * (zion(iti44) - zbar) c..d(ti44)/d(cr48) a1 = ratdum(ircrga) 1 + w1 * ratdum(ircrgp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(icr48) = xsum(icr48) + a1 * bion(iti44) ysum(icr48) = ysum(icr48) - a1 zsum(icr48) = zsum(icr48) + a1 * (zion(iti44) - zbar) c..cr48 jacobian elements c..d(cr48)/d(he4) a1 = y(iti44) * ratdum(irtiag) 1 - y(icr48) * ratdum(ircrag) 2 + y(iti44) * ratdum(irtiap)*(1.0d0-w1) 3 - y(icr48) * ratdum(ircrap)*(1.0d0-x1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(icr48) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(icr48) - zbar) c..d(cr48)/d(ti44) a1 = y(ihe4) * ratdum(irtiag) 1 + y(ihe4) * ratdum(irtiap)*(1.0d0-w1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(iti44) = xsum(iti44) + a1 * bion(icr48) ysum(iti44) = ysum(iti44) - a1 zsum(iti44) = zsum(iti44) + a1 * (zion(icr48) - zbar) c..d(cr48)/d(cr48) a1 = -y(ihe4) * ratdum(ircrag) 1 - ratdum(ircrga) 2 - w1 * ratdum(ircrgp) 3 - y(ihe4) * ratdum(ircrap)*(1.0d0-x1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(icr48) = xsum(icr48) + a1 * bion(icr48) ysum(icr48) = ysum(icr48) - a1 zsum(icr48) = zsum(icr48) + a1 * (zion(icr48) - zbar) c..d(cr48)/d(fe52) a1 = ratdum(irfega) 1 + x1 * ratdum(irfegp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ife52) = xsum(ife52) + a1 * bion(icr48) ysum(ife52) = ysum(ife52) - a1 zsum(ife52) = zsum(ife52) + a1 * (zion(icr48) - zbar) c..fe52 jacobian elements c..d(fe52)/d(he4) a1 = y(icr48) * ratdum(ircrag) 1 - y(ife52) * ratdum(irfeag) 2 + y(icr48) * ratdum(ircrap) * (1.0d0-x1) 3 - y(ife52) * ratdum(irfeap) * (1.0d0-y1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ihe4) = xsum(ihe4) + a1 * bion(ife52) ysum(ihe4) = ysum(ihe4) - a1 zsum(ihe4) = zsum(ihe4) + a1 * (zion(ife52) - zbar) c..d(fe52)/d(cr48) a1 = y(ihe4) * ratdum(ircrag) 1 + y(ihe4) * ratdum(ircrap) * (1.0d0-x1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(icr48) = xsum(icr48) + a1 * bion(ife52) ysum(icr48) = ysum(icr48) - a1 zsum(icr48) = zsum(icr48) + a1 * (zion(ife52) - zbar) c..d(fe52)/d(fe52) a1 = -y(ihe4) * ratdum(irfeag) 1 - ratdum(irfega) 2 - x1 * ratdum(irfegp) 3 - y(ihe4) * ratdum(irfeap) * (1.0d0-y1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ife52) = xsum(ife52) + a1 * bion(ife52) ysum(ife52) = ysum(ife52) - a1 zsum(ife52) = zsum(ife52) + a1 * (zion(ife52) - zbar) c..d(fe52)/d(ni56) a1 = ratdum(irniga) 1 + y1 * ratdum(irnigp) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ini56) = xsum(ini56) + a1 * bion(ife52) ysum(ini56) = ysum(ini56) - a1 zsum(ini56) = zsum(ini56) + a1 * (zion(ife52) - zbar) c..ni56 jacobian elements c..d(ni56)/d(he4) a1 = y(ife52) * ratdum(irfeag) 1 + y(ife52) * ratdum(irfeap) * (1.0d0-y1) 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(fe52) a1 = y(ihe4) * ratdum(irfeag) 1 + y(ihe4) * ratdum(irfeap) * (1.0d0-y1) nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 xsum(ife52) = xsum(ife52) + a1 * bion(ini56) ysum(ife52) = ysum(ife52) - a1 zsum(ife52) = zsum(ife52) + a1 * (zion(ini56) - zbar) c..d(ni56)/d(ni56) a1 = -ratdum(irniga) 1 - y1 * ratdum(irnigp) 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..append the temperature derivatives of the rate equations call rhs(y,dratdumdt,teqn) c..add in the parts from the dummy proton links teqn(ihe4) = teqn(ihe4) 1 + 0.34d0 *s1dt*0.5d0*y(io16)*y(io16)*ratdum(ir1616) 2 + y(ihe4) * y(img24) * ratdum(irmgap)* r1dt 3 + y(isi28) * ratdum(irsigp) * r1dt 4 + y(ihe4) * y(isi28) * ratdum(irsiap) * s1dt 5 + y(is32) * ratdum(irsgp) * s1dt 6 + y(ihe4) * y(is32) * ratdum(irsap) * t1dt 7 + y(iar36) * ratdum(irargp) * t1dt 8 + y(ihe4) * y(iar36) * ratdum(irarap) * u1dt 9 + y(ica40) * ratdum(ircagp) * u1dt & + y(ihe4) * y(ica40) * ratdum(ircaap) * v1dt 1 + y(iti44) * ratdum(irtigp) * v1dt 2 + y(ihe4) * y(iti44) * ratdum(irtiap) * w1dt 3 + y(icr48) * ratdum(ircrgp) * w1dt 4 + y(ihe4) * y(icr48) * ratdum(ircrap) * x1dt 5 + y(ife52) * ratdum(irfegp) * x1dt teqn(ihe4) = teqn(ihe4) 1 + y(ihe4) * y(ife52) * ratdum(irfeap) * y1dt 2 + y(ini56) * ratdum(irnigp) * y1dt teqn(img24) = teqn(img24) 1 + y(img24) * y(ihe4) *ratdum(irmgap)*r1dt 2 + y(isi28) * r1dt * ratdum(irsigp) teqn(isi28) = teqn(isi28) 1 + 0.34d0*0.5d0*y(io16)*y(io16)*s1dt*ratdum(ir1616) 2 - y(img24) * y(ihe4) * ratdum(irmgap)*r1dt 3 - y(isi28) * r1dt * ratdum(irsigp) 4 + y(isi28) * y(ihe4) * ratdum(irsiap)*s1dt 5 + y(is32) * s1dt * ratdum(irsgp) teqn(is32) = teqn(is32) 1 - 0.34d0*0.5d0*y(io16)**2 * ratdum(ir1616)*s1dt 2 - y(isi28) * y(ihe4) * ratdum(irsiap) * s1dt 3 - y(is32) * s1dt * ratdum(irsgp) 4 - y(is32) * y(ihe4) * ratdum(irsap) * t1dt 5 + y(iar36) * t1dt * ratdum(irargp) teqn(iar36) = teqn(iar36) 1 - y(is32) * y(ihe4) * ratdum(irsap) * t1dt 2 - y(iar36) * t1dt * ratdum(irargp) 3 + y(iar36) * y(ihe4) * ratdum(irarap) * u1dt 4 + y(ica40) * ratdum(ircagp) * u1dt teqn(ica40) = teqn(ica40) 1 - y(iar36) * y(ihe4) * ratdum(irarap) * u1dt 2 - y(ica40) * ratdum(ircagp) * u1dt 3 + y(ica40) * y(ihe4) * ratdum(ircaap) * v1dt 4 + y(iti44) * ratdum(irtigp) * v1dt teqn(iti44) = teqn(iti44) 1 - y(ica40) * y(ihe4) * ratdum(ircaap)*v1dt 2 - y(iti44) * v1dt * ratdum(irtigp) 3 + y(iti44) * y(ihe4) * ratdum(irtiap) * w1dt 4 + y(icr48) * w1dt * ratdum(ircrgp) teqn(icr48) = teqn(icr48) 1 - y(iti44) * y(ihe4) * ratdum(irtiap)*w1dt 2 - y(icr48) * w1dt * ratdum(ircrgp) 3 + y(icr48) * y(ihe4) * ratdum(ircrap) * x1dt 4 + y(ife52) * x1dt * ratdum(irfegp) teqn(ife52) = teqn(ife52) 1 - y(icr48) * y(ihe4) * ratdum(ircrap) * x1dt 2 - y(ife52) * x1dt * ratdum(irfegp) 3 + y(ife52) * y(ihe4) * ratdum(irfeap) * y1dt 4 + y(ini56) * y1dt * ratdum(irnigp) teqn(ini56) = teqn(ini56) 1 - y(ife52) * y(ihe4) * ratdum(irfeap) * y1dt 2 - y(ini56) * y1dt * ratdum(irnigp) c..append the density derivatives of the rate equations call rhs(y,dratdumdd,deqn) c..add in the parts from the dummy proton links deqn(ihe4) = deqn(ihe4) 1 + 0.34d0 *s1dd*0.5d0*y(io16)*y(io16)*ratdum(ir1616) 2 + y(ihe4) * y(img24) * ratdum(irmgap)* r1dd 3 + y(isi28) * ratdum(irsigp) * r1dd 4 + y(ihe4) * y(isi28) * ratdum(irsiap) * s1dd 5 + y(is32) * ratdum(irsgp) * s1dd 6 + y(ihe4) * y(is32) * ratdum(irsap) * t1dd 7 + y(iar36) * ratdum(irargp) * t1dd 8 + y(ihe4) * y(iar36) * ratdum(irarap) * u1dd 9 + y(ica40) * ratdum(ircagp) * u1dd & + y(ihe4) * y(ica40) * ratdum(ircaap) * v1dd 1 + y(iti44) * ratdum(irtigp) * v1dd 2 + y(ihe4) * y(iti44) * ratdum(irtiap) * w1dd 3 + y(icr48) * ratdum(ircrgp) * w1dd 4 + y(ihe4) * y(icr48) * ratdum(ircrap) * x1dd 5 + y(ife52) * ratdum(irfegp) * x1dd deqn(ihe4) = deqn(ihe4) 1 + y(ihe4) * y(ife52) * ratdum(irfeap) * y1dd 2 + y(ini56) * ratdum(irnigp) * y1dd deqn(img24) = deqn(img24) 1 + y(img24) * y(ihe4) *ratdum(irmgap)*r1dd 2 + y(isi28) * r1dd * ratdum(irsigp) deqn(isi28) = deqn(isi28) 1 + 0.34d0*0.5d0*y(io16)*y(io16)*s1dd*ratdum(ir1616) 2 - y(img24) * y(ihe4) * ratdum(irmgap)*r1dd 3 - y(isi28) * r1dd * ratdum(irsigp) 4 + y(isi28) * y(ihe4) * ratdum(irsiap)*s1dd 5 + y(is32) * s1dd * ratdum(irsgp) deqn(is32) = deqn(is32) 1 - 0.34d0*0.5d0*y(io16)**2 * ratdum(ir1616)*s1dd 2 - y(isi28) * y(ihe4) * ratdum(irsiap) * s1dd 3 - y(is32) * s1dd * ratdum(irsgp) 4 - y(is32) * y(ihe4) * ratdum(irsap) * t1dd 5 + y(iar36) * t1dd * ratdum(irargp) deqn(iar36) = deqn(iar36) 1 - y(is32) * y(ihe4) * ratdum(irsap) * t1dd 2 - y(iar36) * t1dd * ratdum(irargp) 3 + y(iar36) * y(ihe4) * ratdum(irarap) * u1dd 4 + y(ica40) * ratdum(ircagp) * u1dd deqn(ica40) = deqn(ica40) 1 - y(iar36) * y(ihe4) * ratdum(irarap) * u1dd 2 - y(ica40) * ratdum(ircagp) * u1dd 3 + y(ica40) * y(ihe4) * ratdum(ircaap) * v1dd 4 + y(iti44) * ratdum(irtigp) * v1dd deqn(iti44) = deqn(iti44) 1 - y(ica40) * y(ihe4) * ratdum(ircaap)*v1dd 2 - y(iti44) * v1dd * ratdum(irtigp) 3 + y(iti44) * y(ihe4) * ratdum(irtiap) * w1dd 4 + y(icr48) * w1dd * ratdum(ircrgp) deqn(icr48) = deqn(icr48) 1 - y(iti44) * y(ihe4) * ratdum(irtiap)*w1dd 2 - y(icr48) * w1dd * ratdum(ircrgp) 3 + y(icr48) * y(ihe4) * ratdum(ircrap) * x1dd 4 + y(ife52) * x1dd * ratdum(irfegp) deqn(ife52) = deqn(ife52) 1 - y(icr48) * y(ihe4) * ratdum(ircrap) * x1dd 2 - y(ife52) * x1dd * ratdum(irfegp) 3 + y(ife52) * y(ihe4) * ratdum(irfeap) * y1dd 4 + y(ini56) * y1dd * ratdum(irnigp) deqn(ini56) = deqn(ini56) 1 - y(ife52) * y(ihe4) * ratdum(irfeap) * y1dd 2 - y(ini56) * y1dd * ratdum(irnigp) c..d(yi)/dtemp do i=1,ionmax nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + teqn(i) enddo c..d(yi)/d(den) do i=1,ionmax nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + deqn(i) enddo c..energy jacobian elements c..d(ener)/d(yi) 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 + teqn(i)*bion(i) enddo a1 = a1 * conv nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 enucbt = a1 c..d(iener)/d(den) a1 = 0.0d0 do i=1,ionmax a1 = a1 + deqn(i)*bion(i) enddo a1 = a1 * conv nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 enucbd = a1 c..account for the neutrino losses call sneut5(btemp,bden,abar,zbar, 1 sneut,snudt,snudd,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 = -snudt nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 enucbt = enucbt + a1 c..d(iener)/d(den) a1 = -snudd nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 enucbd = enucbd + a1 c..the jacobian elements of the temperature and density equations c..depend on the burning mode c..hydrostatic if (hydrostatic .or. one_step) 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..d(ivelx)/d(ivelx) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(iposx)/d(iposx) 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..d(ivelx)/d(ivelx) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(iposx)/d(iposx) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..self heating else if (self_heat) 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 + teqn(i)*bion(i) enddo a1 = a1*conv a2 = 0.0d0 do i=1,ionmax a2 = a2 - teqn(i) enddo a2 = a2 * dea_row(1)*abar*abar a3 = 0.0d0 do i=1,ionmax a3 = a3 + teqn(i)*(zion(i) - zbar) enddo a3 = a3 * dez_row(1) * abar a4 = (a1 - snudt - 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 + deqn(i)*bion(i) enddo a1 = a1*conv a2 = 0.0d0 do i=1,ionmax a2 = a2 - deqn(i) enddo a2 = a2 * dea_row(1)*abar*abar a3 = 0.0d0 do i=1,ionmax a3 = a3 + deqn(i)*(zion(i) - zbar) enddo a3 = a3 * dez_row(1) * abar a4 = (a1 - snudd - 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..d(ivelx)/d(ivelx) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(iposx)/d(iposx) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..detonation else if (detonation) then c..get the right hand sides call rhs(y,ratdum,dydt) c do i=1,neqs c write(6,123) i,ionam(i),y(i),dydt(i) c 123 format(1x,i4,' ',a,' ',1p2e11.3) c enddo c read(5,*) c..instantaneous energy generation rate enuc = 0.0d0 do i=1,ionmax enuc = enuc + dydt(i) * bion(i) enddo enuc = enuc * conv dydt(iener) = enuc - sneut c..map the rest of the input vector velx = y(ivelx) posx = y(iposx) c..it appears as if we need the derivatives of derivative based c..eos quantities. grrr. z = bden xx = 0.01d0*z bden = z + xx temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos foo(1) = cs_row(1) foo(2) = dpt_row(1) foo(3) = dpt_row(1)/det_row(1) foo(4) = dpd_row(1) foo(5) = dpa_row(1) foo(6) = dpz_row(1) foo(7) = dea_row(1) foo(8) = dez_row(1) bden = z - xx temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos moo(1) = cs_row(1) moo(2) = dpt_row(1) moo(3) = dpt_row(1)/det_row(1) moo(4) = dpd_row(1) moo(5) = dpa_row(1) moo(6) = dpz_row(1) moo(7) = dea_row(1) moo(8) = dez_row(1) bden = z z = 0.5d0/xx csbd = (foo(1) - moo(1))*z dptbd = (foo(2) - moo(2))*z dpdebd = (foo(3) - moo(3))*z dpdbd = (foo(4) - moo(4))*z dpabd = (foo(5) - moo(5))*z dpzbd = (foo(6) - moo(6))*z deabd = (foo(7) - moo(7))*z dezbd = (foo(8) - moo(8))*z z = btemp xx = 0.01d0*z btemp = z + xx temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos foo(1) = cs_row(1) foo(2) = dpt_row(1) foo(3) = dpt_row(1)/det_row(1) foo(4) = dpd_row(1) foo(5) = dpa_row(1) foo(6) = dpz_row(1) foo(7) = dea_row(1) foo(8) = dez_row(1) btemp = z - xx temp_row(1) = btemp den_row(1) = bden abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos moo(1) = cs_row(1) moo(2) = dpt_row(1) moo(3) = dpt_row(1)/det_row(1) moo(4) = dpd_row(1) moo(5) = dpa_row(1) moo(6) = dpz_row(1) moo(7) = dea_row(1) moo(8) = dez_row(1) btemp = z z = 0.5d0/xx csbt = (foo(1) - moo(1))*z dptbt = (foo(2) - moo(2))*z dpdebt = (foo(3) - moo(3))*z dpdbt = (foo(4) - moo(4))*z dpabt = (foo(5) - moo(5))*z dpzbt = (foo(6) - moo(6))*z deabt = (foo(7) - moo(7))*z dezbt = (foo(8) - moo(8))*z 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..for de/dy and dp/dy suma = 0.0d0 do i=1,ionmax suma = suma - dydt(i) enddo sumz = 0.0d0 do i=1,ionmax sumz = sumz + (zion(i) - zbar)*dydt(i) enddo sum = 0.0d0 do i=1,ionmax sum = sum - deqn(i) enddo sum5 = sum*dea_row(1)*abar*abar sum9 = sum*dpa_row(1)*abar*abar sum = 0.0d0 do i=1,ionmax sum = sum + deqn(i)*(zion(i)-zbar) enddo sum6 = sum*dez_row(1)*abar sum10 = sum*dpz_row(1)*abar sum = 0.0d0 do i=1,ionmax sum = sum - teqn(i) enddo sum7 = sum*dea_row(1)*abar*abar sum11 = sum*dpa_row(1)*abar*abar sum = 0.0d0 do i=1,ionmax sum = sum + teqn(i)*(zion(i)-zbar) enddo sum8 = sum*dez_row(1)*abar sum12 = sum*dpz_row(1)*abar c..the possibly singular denominator cs = cs_row(1) denom = velx*velx - cs*cs denombv = 2.0d0*velx denomdd = -2.0d0*cs*csbd denomdt = -2.0d0*cs*csbt c denom = velx*velx - cs_cj*cs_cj c denomdd = 0.0d0 c denomdt = 0.0d0 c..the function phi dpde = dpt_row(1)/det_row(1) z = suma*dpa_row(1)*abar*abar + sumz*dpz_row(1)*abar zbd = suma*dpabd*abar*abar + sumz*dpzbd*abar + sum9 + sum10 zbt = suma*dpabt*abar*abar + sumz*dpzbt*abar + sum11 + sum12 ww = suma*dea_row(1)*abar*abar + sumz*dez_row(1)*abar wwbd = suma*deabd*abar*abar + sumz*dezbd*abar + sum5 + sum6 wwbt = suma*deabt*abar*abar + sumz*dezbt*abar + sum7 + sum8 phi = dpde*(dydt(iener) - ww) - z phibd = dpdebd*(dydt(iener) - ww) - zbd 1 + dpde*(enucbd - wwbd) phibt = dpdebt*(dydt(iener) - ww) - zbt 1 + dpde*(enucbt - wwbt) c..a common combination if (denom .ne. 0.0) then combo = phi/denom combobv = -combo/denom*denombv combobd = -combo/denom*denomdd + phibd/denom combobt = -combo/denom*denomdt + phibt/denom else combo = 0.0d0 combobv = 0.0d0 combobd = 0.0d0 combobt = 0.0d0 end if c..position equation dydt(iposx) = velx c..d(iposx)/d(iposx) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(iposx)/d(ivelx) a1 = 1.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..density equation dydt(iden) = combo c..d(iden)/d(ivelx) a1 = combobv c a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(iden)/d(iden) a1 = combobd c a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(iden)/d(itemp) a1 = combobt c a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(iden)/d(yi) xx = 1.0d0/denom zz = dpde*xx do i=1,ionmax a1 = zz*(xsum(i) * conv 1 - ysum(i) * dea_row(1)*abar*abar 2 - zsum(i) * dez_row(1)*abar) 3 - xx*(ysum(i) * dpa_row(1)*abar*abar 4 + zsum(i) * dpz_row(1)*abar) c a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 enddo c..velocity equation z = velx/bden dydt(ivelx) = -z*dydt(iden) c..d(ivelx)/d(ivelx) a1 = -dydt(iden)/bden - z*combobv a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(ivelx)/d(iden) a1 = z/bden*dydt(iden) - z*combobd a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(ivelx)/d(itemp) a1 = -z*combobt a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(ivelx)/d(yi) zz = -z * dpde/denom do i=1,ionmax a1 = zz*(xsum(i) * conv 1 - ysum(i) * dea_row(1)*abar*abar 2 - zsum(i) * dez_row(1)*abar) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 enddo c..temperature equation dtdp = 1.0d0/dpt_row(1) dtdpbd = -dtdp*dtdp*dptbd dtdpbt = -dtdp*dtdp*dptbt ww = suma*dpa_row(1)*abar*abar + sumz*dpz_row(1)*abar wwbd = suma*dpabd*abar*abar + sumz*dpzbd*abar + sum9 + sum10 wwbt = suma*dpabt*abar*abar + sumz*dpzbt*abar + sum11 + sum12 dydt(itemp) = dtdp*((velx*velx - dpd_row(1))*dydt(iden) - ww) c..d(itemp)/d(ivelx) a1 = dtdp*(2.0d0*velx*dydt(iden) 1 + (velx*velx - dpd_row(1))*combobv) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(itemp)/d(iden) a1 = dtdpbd*((velx*velx-dpd_row(1))*dydt(iden)-ww) 1 + dtdp*((velx*velx-dpd_row(1))*combobd 2 - dpdbd*dydt(iden) 3 - wwbd) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(itemp)/d(itemp) a1 = dtdpbt*((velx*velx-dpd_row(1))*dydt(iden)-ww) 1 + dtdp*((velx*velx-dpd_row(1))*combobt 2 - dpdbt*dydt(iden) 3 - wwbt) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 c..d(itemp)/d(yi) xx = 1.0d0/denom zz = dpde*xx do i=1,ionmax a3 = zz*(xsum(i) * conv 1 - ysum(i) * dea_row(1)*abar*abar 2 - zsum(i) * dez_row(1)*abar) 3 - xx*(ysum(i) * dpa_row(1)*abar*abar 4 + zsum(i) * dpz_row(1)*abar) a2 = 0.0d0 c a2 = ysum(i) * dpa_row(1)*abar*abar c 1 + zsum(i) * dpz_row(1)*abar a1 = dtdp * ((velx*velx - dpd_row(1)) * a3 - a2) a1 = 0.0d0 nt = nt + 1 iat = eloc(nt) dfdy(iat) = dfdy(iat) + a1 enddo c..end of burning mode ifs end if c..bullet check the counting if (nt .ne. nterms) then write(6,*) 'nt =',nt,' nterms =',nterms write(6,*) 'error in routine saprox13: nt .ne. nterms' stop 'error in routine saprox13' end if return end subroutine aprox13rat include 'implno.dek' include 'burn_common.dek' include 'network.dek' c..this routine generates nuclear reaction rates for the aprox13 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..o16 + o16 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..mg24(a,p)al27 call rate_mg24ap(btemp,bden, 1 ratraw(irmgap),dratrawdt(irmgap),dratrawdd(irmgap), 2 ratraw(iralpa),dratrawdt(iralpa),dratrawdd(iralpa)) c..al27(p,g)si28 call rate_al27pg(btemp,bden, 1 ratraw(iralpg),dratrawdt(iralpg),dratrawdd(iralpg), 2 ratraw(irsigp),dratrawdt(irsigp),dratrawdd(irsigp)) c..si28(a,g)s32 call rate_si28ag(btemp,bden, 1 ratraw(irsiag),dratrawdt(irsiag),dratrawdd(irsiag), 2 ratraw(irsga),dratrawdt(irsga),dratrawdd(irsga)) c..si28(a,p)p31 call rate_si28ap(btemp,bden, 1 ratraw(irsiap),dratrawdt(irsiap),dratrawdd(irsiap), 2 ratraw(irppa),dratrawdt(irppa),dratrawdd(irppa)) c..p31(p,g)s32 call rate_p31pg(btemp,bden, 1 ratraw(irppg),dratrawdt(irppg),dratrawdd(irppg), 2 ratraw(irsgp),dratrawdt(irsgp),dratrawdd(irsgp)) c..s32(a,g)ar36 call rate_s32ag(btemp,bden, 1 ratraw(irsag),dratrawdt(irsag),dratrawdd(irsag), 2 ratraw(irarga),dratrawdt(irarga),dratrawdd(irarga)) c..s32(a,p)cl35 call rate_s32ap(btemp,bden, 1 ratraw(irsap),dratrawdt(irsap),dratrawdd(irsap), 2 ratraw(irclpa),dratrawdt(irclpa),dratrawdd(irclpa)) c..cl35(p,g)ar36 call rate_cl35pg(btemp,bden, 1 ratraw(irclpg),dratrawdt(irclpg),dratrawdd(irclpg), 2 ratraw(irargp),dratrawdt(irargp),dratrawdd(irargp)) c..ar36(a,g)ca40 call rate_ar36ag(btemp,bden, 1 ratraw(irarag),dratrawdt(irarag),dratrawdd(irarag), 2 ratraw(ircaga),dratrawdt(ircaga),dratrawdd(ircaga)) c..ar36(a,p)k39 call rate_ar36ap(btemp,bden, 1 ratraw(irarap),dratrawdt(irarap),dratrawdd(irarap), 2 ratraw(irkpa),dratrawdt(irkpa),dratrawdd(irkpa)) c..k39(p,g)ca40 call rate_k39pg(btemp,bden, 1 ratraw(irkpg),dratrawdt(irkpg),dratrawdd(irkpg), 2 ratraw(ircagp),dratrawdt(ircagp),dratrawdd(ircagp)) 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..ca40(a,p)sc43 call rate_ca40ap(btemp,bden, 1 ratraw(ircaap),dratrawdt(ircaap),dratrawdd(ircaap), 2 ratraw(irscpa),dratrawdt(irscpa),dratrawdd(irscpa)) c..sc43(p,g)ti44 call rate_sc43pg(btemp,bden, 1 ratraw(irscpg),dratrawdt(irscpg),dratrawdd(irscpg), 2 ratraw(irtigp),dratrawdt(irtigp),dratrawdd(irtigp)) c..ti44(a,g)cr48 call rate_ti44ag(btemp,bden, 1 ratraw(irtiag),dratrawdt(irtiag),dratrawdd(irtiag), 2 ratraw(ircrga),dratrawdt