c
c =====================================================
subroutine setaux(maxmx,meqn,mbc,ibx,mx,q,
& aux,maux,xc,dx,t,dt)
c =====================================================
c
c # Copyright (C) 2002 Ralf Deiterding
c # Brandenburgische Universitaet Cottbus
c
implicit double precision (a-h, o-z)
c
include "ck.i"
c
dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc)
dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc)
dimension rk(LeNsp), rkl(LeNsp)
c
mu = Nsp+1
mE = Nsp+2
mT = Nsp+3
c
do 5 i = 1-ibx*mbc, mx+ibx*mbc
rho = 0.d0
rhoW = 0.d0
do k = 1, Nsp
rk(k) = q(k,i)
rho = rho + q(k,i)
rhoW = rhoW + q(k,i)/Wk(k)
enddo
rhoe = q(mE,i)-0.5d0*q(mu,i)**2/rho
call SolveTrhok(q(mT,i),rhoe,rhoW,rk,Nsp,ifail)
5 continue
c
aux(1,1-ibx*mbc) = 0.d0
do 10 i = 2-ibx*mbc, mx+ibx*mbc
rho = 0.d0
rhoW = 0.d0
do k = 1, Nsp
rk(k) = q(k,i)
rho = rho + q(k,i)
rhoW = rhoW + q(k,i)/Wk(k)
enddo
u = q(mu,i)/rho
T0 = q(mT,i)
p = rhoW*RU*T0
rhoCp = avgtabip( T0, rk, cpk, Nsp )
gamma = RU / ( rhoCp/rhoW - RU ) + 1.d0
a = dsqrt(gamma*p/rho)
c
rhol = 0.d0
rhoWl = 0.d0
do k = 1, Nsp
rkl(k) = q(k,i-1)
rhol = rhol + q(k,i-1)
rhoWl = rhoWl + q(k,i-1)/Wk(k)
enddo
ul = q(mu,i-1)/rhol
Tl = q(mT,i-1)
pl = rhoWl*RU*Tl
rhoCpl = avgtabip( Tl, rkl, cpk, Nsp )
gammal = RU / ( rhoCpl/rhoWl - RU ) + 1.d0
al = dsqrt(gammal*pl/rhol)
aux(1,i) = dmax1(dabs(u-a-(ul-al)),dabs(u-ul),
& dabs(u+a-(ul+al)))
10 continue
c
return
end
c