c
c Copyright (C) 2002 Ralf Deiterding
c Brandenburgische Universitaet Cottbus
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"
include "call.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), qint(LeNsp)
external prodrhokAppTrho, Jprod, ownprod, odeoutput
c
common /grkstat/ jsuc,jfail,jf,ja,jsing
common /output/ tout
c
logical mech
data mech /.true./ !# turn off mechanism
c
tol = 1.d-05
sclimit = 1.d-02
uround = 5.d-15
tend = dt
tstep = dt
tstart = 0.d0
tout = t
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
hmax = tstep
c
rho = 0.d0
rhoW = 0.d0
do m=1,Nsp
rho = rho + q(m,i,j,k)
rhoW = rhoW + q(m,i,j,k)/Wk(m)
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
NTrys = 0
c
20 continue
do m=1,Nsp
qint(m) = q(m,i,j,k)
enddo
c
if (q(Nsp+6,i,j,k).gt.0.d0.and.
& q(Nsp+6,i,j,k).lt.tstep/1.5d0) then
tst = q(Nsp+6,i,j,k)
else
tst = tstep
endif
do n=1,NTrys
tst = tst / 10.d0
if (tst.lt.uround) then
write(6,400) i,j,k,t,dt
goto 30
endif
enddo
if (NCal.eq.1) then
call grk4a(Nint,ownprod,Jprod,ts,qint,te,tol,
& hmax,tst,lengrkaux,grkaux,ichan,ijac,
& odeoutput,iscale,scale,lengrkw,grkwork,
& lenint,igrkwork)
else
call grk4a(Nint,prodrhokAppTrho,Jprod,ts,qint,te,
& tol,hmax,tst,lengrkaux,grkaux,ichan,ijac,
& odeoutput,iscale,scale,lengrkw,grkwork,
& lenint,igrkwork)
endif
c
c !# Failure, try again with reduced step size
c
if (tst.eq.0.d0) then
NTrys = NTrys+1
goto 20
endif
do m=1,Nsp
q(m,i,j,k) = qint(m)
enddo
c
q(Nsp+5,i,j,k) = grkaux(2)
q(Nsp+6,i,j,k) = tst
c
30 continue
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)
if (ifail.ne.0) write (6,401) i,j,k,q(Nsp+5,i,j,k)
c
10 continue
c
400 format(' Source integration failure in (',i5,',',
& i5,',',i5,') t=',f20.12,' dt=',f20.12)
401 format('src (',i5,',',i5,',',i5,'): T out of range using: ',f16.8)
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
c **************************************************************
c Production rates in partial densities
c adiabatic, constant volume, iterative temperature
c Optimized function
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)
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)
if (ifcnfail.ne.0) write (6,600) aux(2)
call genckwcfs (aux(2),Ck,drhok)
c
do k = 1, Nsp
drhok(k) = drhok(k)*Wk(k)
enddo
c
600 format('T out of range in ownprod() using: ',f16.8)
return
end
c