c
c ===================================================================
subroutine fmod1(maxmx,mvar,meqn,maux,mbc,mx,q,ql,qr,
& aux,dx,dt,method,cfl,fm,fp,q1,aux1,
& amdq,apdq)
c ===================================================================
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
implicit double precision (a-h,o-z)
include "ck.i"
include "cuser.i"
c
parameter (lenma=46, maxm2 = 10002) !# assumes at most 10000 grid points with mbc=2
dimension q(mvar, 1-mbc:maxmx+mbc)
dimension ql(mvar, 1-mbc:maxmx+mbc)
dimension qr(mvar, 1-mbc:maxmx+mbc)
dimension fm(mvar, 1-mbc:maxmx+mbc)
dimension fp(mvar, 1-mbc:maxmx+mbc)
dimension aux(maux, 1-mbc:maxmx+mbc)
dimension q1(1-mbc:maxmx+mbc, meqn)
dimension aux1(1-mbc:maxmx+mbc, maux)
dimension amdq(1-mbc:maxmx+mbc, meqn)
dimension apdq(1-mbc:maxmx+mbc, meqn)
dimension method(7)
dimension a(lenma, -1:maxm2)
dimension fdiv(leNsp+2), hR(leNsp), hL(leNsp)
c
logical tdiff
data tdiff /.false./ !# activate thermal diffusion
c
mu = Nsp+1
mE = Nsp+2
mT = Nsp+3
c
marho = 1
mau = marho+1
mapr = mau+1
maW = mapr+1
macp = maW+1
mavis = macp+1
macon = mavis+1
maX1 = macon+1
maXN = maX1+Nsp-1
maD1 = maXN+1
maDN = maD1+Nsp-1
maTD1 = maDN+1
maTDN = maTD1+Nsp-1
c
if (lenma.lt.maTDN) then
write (6,*) 'Set lenma in fnavs1g to ',maTDN,'!'
stop
endif
c
c # Compute derived quantities
do 10 i = 1-mbc, mx+mbc
c # Calculate total density
a(marho,i) = 0.d0
do k=1,Nsp
a(marho,i) = a(marho,i)+q(k,i)
enddo
a(mau,i) = q(mu,i)/a(marho,i)
c # Calculate mass fractions
do k=1,Nsp
a(maX1+k-1,i) = q(k,i)/a(marho,i)
enddo
call ckmmwy(a(maX1,i),iwork,rwork,a(maW,i))
a(macp,i) = avgtabip(q(mT,i),a(maX1,i),cpk,Nsp)
c
c # Calculate mole fractions
do k=1,Nsp
a(maX1+k-1,i) = a(maX1+k-1,i)*a(maW,i)/Wk(k)
enddo
c
a(mapr,i) = a(marho,i)/a(maW,i)*RU*q(mT,i)
c # Mixture viscosity
call mcavis(q(mT,i),a(maX1,i),rmcwork,a(mavis,i))
c # Mixture diffusion coefficients
call mcadif(a(mapr,i),q(mT,i),a(maX1,i),rmcwork,a(maD1,i))
c # Mixture thermal conductivity
call mcacon(q(mT,i),a(maX1,i),rmcwork,a(macon,i))
c # Thermal diffusion ratios for light species
if (tdiff)
& call mcatdr(q(mT,i),a(maX1,i),imcwork,rmcwork,a(maTD1,i))
c
10 continue
c
fac = 2.d0*dt/(dx**2)
do 110 i = 1, mx+1
c # Difference quotients
ux = (a(mau,i)-a(mau,i-1))/dx
Tx = (q(mT,i)-q(mT,i-1))/dx
px = (a(mapr,i)-a(mapr,i-1))/dx
c
c # Upwinding based on velocity u, if both sides have
c # the same velocity, otherwise use average
iR = i
iL = i-1
c if (a(mau,i).ge.0.d0.and.a(mau,i-1).ge.0.d0) iR = i-1
c if (a(mau,i).lt.0.d0.and.a(mau,i-1).lt.0.d0) iL = i
c
c # Construct intermediate values and material properties
rho = 0.5d0*(a(marho,iR)+a(marho,iL))
u = 0.5d0*(a(mau,iR)+a(mau,iL))
T = 0.5d0*(q(mT,iR)+q(mT,iL))
p = 0.5d0*(a(mapr,iR)+a(mapr,iL))
W = 0.5d0*(a(maW,iR)+a(maW,iL))
vis = 0.5d0*(a(mavis,iR)+a(mavis,iL))
con = 0.5d0*(a(macon,iR)+a(macon,iL))
cp = 0.5d0*(a(macp,iR)+a(macp,iL))
c
call tabintp(q(mT,iR),hR,hms,Nsp)
call tabintp(q(mT,iL),hL,hms,Nsp)
c
c # Diffusive momentum flux
fdiv(mu) = (4.d0/3.d0)*vis*ux
c # Thermal conductivity, energy change due to diffusive momentum
fdiv(mE) = con*Tx+u*fdiv(mu)
c
c # Stability condition from diffusive momentum and energy flux
cflvis = max(fac*dabs(vis/rho),fac*dabs(con/(rho*cp)))
c
do k=1, Nsp
c # Species-dependent difference quotient
Xx = (a(maX1+k-1,i)-a(maX1+k-1,i-1))/dx
c # Construct species-dependent intermediate values
Xk = 0.5d0*(a(maX1+k-1,iR)+a(maX1+k-1,iL))
rhok = 0.5d0*(q(k,iR)+q(k,iL))
Dmix = 0.5d0*(a(maD1+k-1,iR)+a(maD1+k-1,iL))
c Dmix = (a(maD1+k-1,iR)*a(maD1+k-1,iL))/
c & (0.5d0*(a(maD1+k-1,iR)+a(maD1+k-1,iL)))
TDmix = 0.5d0*(a(maTD1+k-1,iR)+a(maTD1+k-1,iL))
hk = 0.5d0*(hR(k)+hL(k))
c
c # Diffusion driving force
dk = Xx + (Xk-rhok/rho) * px/p
c
c # Diffusive multi-species flux
fdiv(k) = dk
c # add thermal diffusion
if (tdiff) fdiv(k) = fdiv(k) + TDmix*Tx/T
fdiv(k) = fdiv(k) * rho*Wk(k)*Dmix/W
c # incorporate multi-species flux into energy flux
fdiv(mE) = fdiv(mE) + hk*fdiv(k)
c # Dufour effect
if (Xk.gt.0.d0.and.tdiff)
& fdiv(mE) = fdiv(mE) + p*Dmix*TDmix*dk/Xk
cflvis = max(cflvis, fac*dabs(Dmix))
enddo
c
cfl = max(cfl,cflvis)
c
c # Add diffusive flux approximation
do k = 1, mE
fm(k,i) = fm(k,i) - fdiv(k)
fp(k,i) = fp(k,i) - fdiv(k)
enddo
c
110 continue
end