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@cacr.caltech.edu 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