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 # Integrate reaction terms in time in every cell.
c =========================================================
subroutine src(maxmx,meqn,mbc,ibx,mx,q,
& aux,maux,t,dt,ibnd)
c =========================================================
implicit double precision(a-h,o-z)
include "ck.i"
include "cuser.i"
include "call.i"
parameter ( lengrkaux = 2 )
c
dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc)
dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc)
dimension grkaux(lengrkaux), scale(LeNsp), qint(LeNsp)
external prod, Jprod, odeout
c
common /grkstat/ jsuc,jfail,jf,ja
common /output/ tout
c
logical mech
data mech /.true./ !# turn of mechanism
c
tol = 1.d-06
sclimit = 1.d-04
uround = 5.d-15
tend = dt
tstep = dt
tstart = 0.d0
tout = t
ichan = 0
ijac = 0
iscale = 1
c
do 10 i=1-ibnd*ibx*mbc,mx+ibnd*ibx*mbc
jsuc = 0
jfail = 0
jf = 0
ja = 0
c
ts = tstart
te = tend
hmax = tstep
c
c # Evaluate internal energy density and temperate
c # and pass them in array grkaux
rho = 0.d0
rhoW = 0.d0
do k=1,Nsp
rho = rho + q(k,i)
rhoW = rhoW + q(k,i)/Wk(k)
enddo
grkaux(1) = q(Nsp+2,i)-(q(Nsp+1,i)**2/2.d0)/rho
grkaux(2) = q(Nsp+3,i)
c
if (mech) then
do k=1,Nsp
scale(k) = dmax1(rho*sclimit,q(k,i))
enddo
Nint = Nsp
NTrys = 0
c
20 continue
do k=1,Nsp
qint(k) = q(k,i)
enddo
c
if (q(Nsp+4,i).gt.0.and.q(Nsp+4,i).lt.tstep/1.5d0) then
tst = q(Nsp+4,i)
else
tst = tstep
endif
do n=1,NTrys
tst = tst / 10.d0
if (tst.lt.uround) then
write(6,400) i,t,dt
goto 30
endif
enddo
call grk4a(Nint,prod,Jprod,ts,qint,te,tol,
& hmax,tst,lengrkaux,grkaux,ichan,ijac,odeout,
& iscale,scale,lengrkw,grkwork,lenint,igrkwork)
c
c !# Failure, try again with reduced step size
c
if (tst.eq.0.d0) then
NTrys = NTrys+1
goto 20
endif
c
c # Update partial densities and temperature from
c # calculated values
c
do k=1,Nsp
q(k,i) = qint(k)
enddo
q(Nsp+3,i) = grkaux(2)
q(Nsp+4,i) = tst
q(Nsp+5,i) = dble(jsuc)
q(Nsp+6,i) = dble(jfail)
c
30 continue
endif
c
rhoW = 0.d0
do k = 1, Nsp
rhoW = rhoW + q(k,i)/Wk(k)
enddo
call SolveTrhok(q(Nsp+3,i),grkaux(1),rhoW,q(1,i),Nsp,ifail)
c
10 continue
c
400 format(' Source integration failure in (',i5,
& ') t=',f20.12,' dt=',f20.12)
return
end
c
c
c # Production rates in partial densities
c # adiabatic, constant volume, iterative temperature
c =========================================================
subroutine prod(n,t,rhok,drhok,na,aux,ifcnfail)
c =========================================================
implicit double precision (a-h,o-z)
include "ck.i"
dimension rhok(n),drhok(n),aux(na),Ck(LeNsp)
c
ifcnfail = 0
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 prodckwc(aux(2),Ck,drhok,n)
do k = 1, Nsp
drhok(k) = drhok(k)*Wk(k)
enddo
c
return
end
c
c
c # Returns the molar production rates given the
c # temperature and molar concentrations
c =========================================================
subroutine prodckwc(T,Ck,dCk,n)
c =========================================================
implicit double precision (a-h,o-z)
include "ck.i"
include "cuser.i"
dimension Ck(n),dCk(n)
c
if (NCal.eq.1) then
call genckwcfs (T,Ck,dCk)
else
call ckwc(T,Ck,iwork,rwork,dCk)
endif
c
return
end