c
c Copyright (C) 2002 Ralf Deiterding
c Brandenburgische Universitaet Cottbus
c
c =====================================================
subroutine setaux(maxmx,maxmy,meqn,mbc,ibx,iby,mx,my,q,
& aux,maux,xc,yc,dx,dy,t,dt)
c =====================================================
c
implicit double precision (a-h, o-z)
include "cuser.i"
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)
dimension u(3), ul(3)
c
do 10 j = 1-iby*mbc, my+iby*mbc
do 10 i = 1-ibx*mbc, mx+ibx*mbc
call eigenval(meqn,q(1,i,j),u)
c
if (i.gt.1-ibx*mbc) then
call eigenval(meqn,q(1,i-1,j),ul)
aux(3,i,j) = dmax1(dmax1(dabs(u(2)-u(1)-(ul(2)-ul(1))),
& dabs(u(2)-ul(2))),
& dabs(u(2)+u(1)-(ul(2)+ul(1))))
endif
c
if (j.gt.1-iby*mbc) then
call eigenval(meqn,q(1,i,j-1),ul)
aux(4,i,j) = dmax1(dmax1(dabs(u(3)-u(1)-(ul(3)-ul(1))),
& dabs(u(3)-ul(3))),
& dabs(u(3)+u(1)-(ul(3)+ul(1))))
endif
10 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))
write (6,*) i,j,aux(1,i,j),aux(2,i,j)
20 continue
c
return
end
c
c =====================================================
subroutine eigenval(meqn,q,ul)
c =====================================================
c
implicit double precision (a-h, o-z)
include "cuser.i"
c
dimension q(meqn), ul(3)
c
p = gamma1*(q(4) - 0.5d0*(q(2)**2 + q(3)**2)/q(1))
ul(1) = dsqrt(gamma*p/q(1))
ul(2) = q(2)/q(1)
ul(3) = q(3)/q(1)
c
return
end
c