c
c # CJ volume burn in the line of pages 316-317,
c # Charles L. Mader, Numerical Modelling of Detonations,
c # Los Alamos series in basic and applied sciences,
c # University of California Press, Berkeley and
c # Los Angeles, 1979.
c # Use this model only with Method(5)=1 or 3 and
c # do not use a Riemann solver for the fluid part.
c
c # Ynew satisfies
c # rho.lt.rho0 -> 1.d0.lt.Ynew
c # rho0.le.rho.le.rhoJ -> 0.d0.le.Ynew.lt.1.d0
c # rhoJ.lt.rho -> Ynew.lt.0.d0
c # which means burning for Ynew.le.1.d0
c
c Copyright (C) 2003-2007 California Institute of Technology
c Ralf Deiterding, ralf@amroc.net
c
c =========================================================
subroutine src(maxmx,maxmy,maxmz,meqn,mbc,ibx,iby,ibz,
& mx,my,mz,q,aux,maux,t,dt,ibnd)
c =========================================================
implicit double precision(a-h,o-z)
include "cuser.i"
c
dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc,1-ibz*mbc:maxmz+ibz*mbc)
dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc,1-ibz*mbc:maxmz+ibz*mbc)
c
do 10 k=1-ibnd*ibz*mbc,mz+ibnd*ibz*mbc
do 10 j=1-ibnd*iby*mbc,my+ibnd*iby*mbc
do 10 i=1-ibnd*ibx*mbc,mx+ibnd*ibx*mbc
c
if (aux(1,i,j,k).le.zt1.and.
& aux(2,i,j,k).gt.rf) then
q(1,i,j,k) = (1.d0-Yo)*rhoo
q(2,i,j,k) = Yo*rhoo
q(3,i,j,k) = 0.d0
q(4,i,j,k) = 0.d0
q(5,i,j,k) = uo*rhoo
q(6,i,J,k) = po/gamma1 + 0.5d0*rhoo*uo**2
goto 5
endif
c
rho=q(1,i,j,k)+q(2,i,j,k)
u=q(3,i,j,k)/rho
v=q(4,i,j,k)/rho
w=q(5,i,j,k)/rho
Yold=q(2,i,j,k)/rho
Ynew=1.d0-(1.d0/rho-V0)/(Vj-V0)
if (Ynew.gt.1.d0.or.Yold.lt.Yact)
& goto 5
if (Yold.lt.Ynew.and.Ynew.lt.0.5d0)
& Ynew = 0.d0
if (Ynew.lt.0.d0) Ynew = 0.d0
if (Yold.le.Ynew) goto 5
Ps = gamma1*(q(6,i,j,k)-q(2,i,j,k)*q0-
& 0.5d0*rho*(u**2+v**2+w**2))
if (Ynew.lt.0.9d0) Ps = (1.d0-Ynew)*Pj
q(2,i,j,k)=rho*Ynew
q(1,i,j,k)=rho-q(2,i,j,k)
q(6,i,j,k)=Ps/gamma1+q(2,i,j,k)*q0+
& 0.5d0*rho*(u**2+v**2+w**2)
5 continue
10 continue
c
return
end