c-----------------------------------------------------------------------
c # Update temperature in 3d case
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
subroutine TUpdate3(q,mx,my,mz,lb,ub,lbr,ubr,shaper,meqn)
c
implicit double precision (a-h,o-z)
c
include "ck.i"
c
integer meqn, mx, my, mz
double precision q(meqn,mx,my,mz)
integer lb(3), ub(3), lbr(3), ubr(3), shaper(3), getindx
integer i, j, k, m, ic, jc, kc, stride
c
if (meqn.le.6) return
c
stride = (ub(1) - lb(1))/(mx-1)
c
do 10 kc=lbr(3), ubr(3), stride
k = getindx(kc, lb(3), stride)
do 10 jc=lbr(2), ubr(2), stride
j = getindx(jc, lb(2), stride)
do 10 ic=lbr(1), ubr(1), stride
i = getindx(ic, lb(1), stride)
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
rhoek = 0.5d0*(q(Nsp+1,i,j,k)**2+
& q(Nsp+2,i,j,k)**2+q(Nsp+3,i,j,k)**2)/rho
rhoe = q(Nsp+4,i,j,k) - rhoek
call SolveTrhok(q(Nsp+5,i,j,K), rhoe, rhoW,
& q(1,i,j,K), Nsp, ifail)
if (ifail.ne.0) then
p = rhoW*RU*q(Nsp+5,i,j,k)
q(Nsp+4,i,j,k) = rhoek +
& avgtabip(T,q(1,i,j,k),hms,Nsp) - p
endif
c
10 continue
c
return
end