c
c Copyright (C) 2002 Ralf Deiterding
c Brandenburgische Universitaet Cottbus
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 "ck.i"
include "cuser.i"
parameter ( lengrkaux = 2 )
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)
dimension grkaux(lengrkaux), scale(LeNsp)
external prodrhokAppTrho, ownprod, Jprod, odeoutput
c
common /grkstat/ jsuc,jfail,jf,ja,jsing
common /output/ tout
c
logical mech, opt
data mech /.true./ !# turn of mechanism
data opt /.true./ !# use optimized function
c
tol = 1.d-04
sclimit = 1.d-03
tend = dt
tstep = dt
tstart = 0.d0
tout = t-dt
ichan = 0
ijac = 0
iscale = 1
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
jsuc = 0
jfail = 0
jf = 0
ja = 0
c
ts = tstart
te = tend
tst = tstep
hmax = tstep
c
rho = 0.d0
do m=1,Nsp
rho = rho+q(m,i,j,k)
enddo
grkaux(1) = q(Nsp+4,i,j,k)-(q(Nsp+1,i,j,k)**2+
& q(Nsp+2,i,j,k)**2+q(Nsp+3,i,j,k)**2)/(2.d0*rho)
grkaux(2) = q(Nsp+5,i,j,k)
c
if (mech) then
do m=1,Nsp
scale(m) = dmax1(rho*sclimit,q(m,i,j,k))
enddo
Nint = Nsp
if (opt) then
call grk4a(Nint,ownprod,Jprod,ts,q(1,i,j,k),te,
& tol,hmax,tst,lengrkaux,grkaux,ichan,ijac,
& odeoutput,iscale,scale,lengrkw,grkwork,
& lenint,igrkwork)
else
call grk4a(Nint,prodrhokAppTrho,Jprod,ts,
& q(1,i,j,k),te,tol,hmax,tst,lengrkaux,grkaux,
& ichan,ijac,odeoutput,iscale,scale,lengrkw,
& grkwork,lenint,igrkwork)
endif
q(Nsp+5,i,j,k) = grkaux(2)
endif
c
rhoW = 0.d0
do m = 1, Nsp
rhoW = rhoW + q(m,i,j,k)/Wk(m)
enddo
call SolveTrhok(q(Nsp+5,i,j,k),grkaux(1),rhoW,
& q(1,i,j,k),Nsp,ifail)
c
10 continue
return
end
c
c
c
c =========================================================
subroutine odeoutput(n,t,rhok,naux,aux,ind)
c =========================================================
implicit double precision(a-h,o-z)
include "ck.i"
include "cuser.i"
dimension rhok(n), aux(naux)
common /grkstat/ jsuc,jfail,jf,ja,jsing
common /output/ tout
return
end
c
c
c =========================================================
c Production rates in partial densities
c adiabatic, constant volume, iterative temperature
c Optimized function
c =========================================================
c
subroutine ownprod(n,t,rhok,drhok,na,aux,ifcnfail)
implicit double precision (a-h,o-z)
include "ck.i"
dimension rhok(n),drhok(n),aux(na),Ck(LeNsp)
external genckwcfs
c
ifcnfail = 0
c
rhoW = 0.d0
do k = 1, Nsp
Ck(k) = rhok(k)/Wk(k)
rhoW = rhoW + Ck(k)
enddo
call SolveTrhok(aux(2),aux(1),rhoW,rhok,Nsp,ifcnfail)
call genckwcfs (aux(2),Ck,drhok)
c
do k = 1, Nsp
drhok(k) = drhok(k)*Wk(k)
enddo
c
return
end
c