c
c Copyright (C) 2002 Ralf Deiterding
c Brandenburgische Universitaet Cottbus
c
c =========================================================
subroutine src(maxmx,maxmy,meqn,mbc,ibx,iby,mx,my,q,
& aux,maux,t,dt,ibnd)
c =========================================================
implicit double precision(a-h,o-z)
include "cuser.i"
include "ck.i"
parameter ( lengrkaux = 2 )
c
dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc)
dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*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 i=1-ibnd*ibx*mbc,mx+ibnd*ibx*mbc
do 10 j=1-ibnd*iby*mbc,my+ibnd*iby*mbc
jsuc = 0
jfail = 0
jf = 0
ja = 0
c
ts = tstart
te = tend
tst = tstep
hmax = tstep
c
rho = 0.d0
do k=1,Nsp
rho = rho+q(k,i,j)
enddo
grkaux(1) = q(Nsp+3,i,j)-
& (q(Nsp+1,i,j)**2+q(Nsp+2,i,j)**2)/(2.d0*rho)
grkaux(2) = q(Nsp+4,i,j)
c
if (mech) then
do k=1,Nsp
scale(k) = dmax1(rho*sclimit,q(k,i,j))
enddo
Nint = Nsp
if (opt) then
call grk4a(Nint,ownprod,Jprod,ts,q(1,i,j),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),
& te,tol,hmax,tst,lengrkaux,grkaux,ichan,ijac,
& odeoutput,iscale,scale,lengrkw,grkwork,lenint,
& igrkwork)
endif
q(Nsp+4,i,j) = grkaux(2)
endif
c
rhoW = 0.d0
do k = 1, Nsp
rhoW = rhoW + q(k,i,j)/Wk(k)
enddo
call SolveTrhok(q(Nsp+4,i,j),grkaux(1),rhoW,q(1,i,j),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 "cuser.i"
include "ck.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