c
c =========================================================
subroutine flx3(ixyz,maxmx,meqn,mbc,mx,q,maux,aux,f)
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@cacr.caltech.edu
c
implicit double precision (a-h,o-z)
include "ck.i"
c
dimension q(1-mbc:maxmx+mbc, meqn)
dimension aux(1-mbc:maxmx+mbc, maux)
dimension f(1-mbc:maxmx+mbc, meqn)
dimension rk(LeNsp)
c
if (ixyz .eq. 1) then
mu = Nsp+1
mv = Nsp+2
mw = Nsp+3
else if(ixyz .eq. 2)then
mu = Nsp+2
mv = Nsp+3
mw = Nsp+1
else
mu = Nsp+3
mv = Nsp+1
mw = Nsp+2
endif
mE = Nsp+4
mT = Nsp+5
c
do i=1-mbc,mx+mbc
rho = 0.d0
rhoW = 0.d0
do m = 1, Nsp
rk(m) = q(i,m)
rho = rho + q(i,m)
rhoW = rhoW + q(i,m)/Wk(m)
enddo
u = q(i,mu)/rho
v = q(i,mv)/rho
w = q(i,mw)/rho
rhoe = q(i,mE)-0.5d0*rho*(u**2+v**2+w**2)
call SolveTrhok(q(i,mT),rhoe,rhoW,rk,Nsp,ifail)
if (ifail.ne.0) write(6,100) i,q(i,mT)
p = rhoW*RU*q(i,mT)
do m = 1, Nsp
f(i,m) = u*q(i,m)
enddo
f(i,mu) = u*q(i,mu)+p
f(i,mv) = u*q(i,mv)
f(i,mw) = u*q(i,mw)
f(i,mE) = u*(q(i,mE)+p)
f(i,mT) = 0.d0
enddo
c
100 format('flx3eurhok q(',i2,'): T out of range using: ',f16.8)
return
end
c