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