c
c =====================================================
subroutine setaux(maxmx,maxmy,meqn,mbc,ibx,iby,mx,my,q,
& aux,maux,xc,yc,dx,dy,t,dt)
c =====================================================
c
c # Multi-dimensional computation of wave speeds, used in
c # multi-dimensional entropy correction in rpn2euefix.
c
c # Copyright (C) 2002 Ralf Deiterding
c # Brandenburgische Universitaet Cottbus
c
implicit double precision (a-h, o-z)
c
dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc)
dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc)
common /param/ gamma,gamma1,q0
c
do 10 j = 1-iby*mbc, my+iby*mbc
do 10 i = 1-ibx*mbc, mx+ibx*mbc
rho = q(1,i,j) + q(2,i,j)
p = gamma1*(q(5,i,j) - q(2,i,j)*q0 -
& 0.5d0*(q(3,i,j)**2+q(4,i,j)**2)/rho)
a = dsqrt(gamma*p/rho)
u = q(3,i,j)/rho
v = q(4,i,j)/rho
c
if (i.gt.1-ibx*mbc) then
rhol = q(1,i-1,j) + q(2,i-1,j)
pl = gamma1*(q(5,i-1,j) - q(2,i-1,j)*q0 -
& 0.5d0*(q(3,i-1,j)**2 + q(4,i-1,j)**2)/rhol)
al = dsqrt(gamma*pl/rhol)
ul = q(3,i-1,j)/rhol
aux(3,i,j) = dmax1(dmax1(dabs(u-a-(ul-al)),dabs(u-ul)),
& dabs(u+a-(ul+al)))
else
aux(3,i,j) = 0.d0
endif
c
if (j.gt.1-iby*mbc) then
rhod = q(1,i,j-1) + q(2,i,j-1)
pd = gamma1*(q(5,i,j-1) - q(2,i,j-1)*q0 -
& 0.5d0*(q(3,i,j-1)**2 + q(4,i,j-1)**2)/rhod)
ad = dsqrt(gamma*pd/rhod)
vd = q(4,i,j-1)/rhod
aux(4,i,j) = dmax1(dmax1(dabs(v-a-(vd-ad)),dabs(v-vd)),
& dabs(v+a-(vd+ad)))
else
aux(4,i,j) = 0.d0
endif
10 continue
c
do 11 j = 1-iby*mbc, my+iby*mbc
aux(1,1-ibx*mbc,j) = 0.d0
aux(2,1-ibx*mbc,j) = 0.d0
11 continue
c
do 12 i = 1-ibx*mbc, mx+ibx*mbc
aux(1,i,1-iby*mbc) = 0.d0
aux(2,i,1-iby*mbc) = 0.d0
12 continue
c
do 20 j = 2-iby*mbc, my+iby*mbc
do 20 i = 2-ibx*mbc, mx+ibx*mbc
aux(1,i,j) = dmax1(dmax1(aux(3,i,j),aux(4,i ,j )),
& aux(4,i-1,j ))
if (j.lt.my+iby*mbc)
& aux(1,i,j) = dmax1(dmax1(aux(1,i,j),aux(4,i ,j+1)),
& aux(4,i-1,j+1))
c
aux(2,i,j) = dmax1(dmax1(aux(4,i,j),aux(3,i ,j )),
& aux(3,i ,j-1))
if (i.lt.mx+ibx*mbc)
& aux(2,i,j) = dmax1(dmax1(aux(2,i,j),aux(3,i+1,j )),
& aux(3,i+1,j-1))
20 continue
c
return
end
c