c
c
c =====================================================
subroutine rpt2eurhok(ixy,maxm,meqn,mwaves,mbc,mx,
& ql,qr,maux,aux1,aux2,aux3,
& ilr,asdq,bmasdq,bpasdq)
c =====================================================
implicit double precision (a-h,o-z)
c
c # Riemann solver in the transverse direction for the Euler equations.
c # Split asdq (= A^* \Delta q, where * = + or -)
c # into down-going flux difference bmasdq (= B^- A^* \Delta q)
c # and up-going flux difference bpasdq (= B^+ A^* \Delta q)
c
c # Uses Roe averages and other quantities which were
c # computed in rpn2eu and stored in the common block comroe.
c
c # Copyright (C) 2002 Ralf Deiterding
c # Brandenburgische Universitaet Cottbus
c
include "ck.i"
c
dimension ql(1-mbc:maxm+mbc, meqn)
dimension qr(1-mbc:maxm+mbc, meqn)
dimension asdq(1-mbc:maxm+mbc, meqn)
dimension bmasdq(1-mbc:maxm+mbc, meqn)
dimension bpasdq(1-mbc:maxm+mbc, meqn)
c
parameter (maxm2 = 10005) !# assumes at most 10000x10000 grid with mbc=5
parameter (minm2 = -4) !# assumes at most mbc=5
common /comroe/ u(minm2:maxm2), v(minm2:maxm2), u2v2(minm2:maxm2),
& enth(minm2:maxm2), a(minm2:maxm2), g1a2(minm2:maxm2),
& dpY(minm2:maxm2), Y(LeNsp,minm2:maxm2), pk(LeNsp,minm2:maxm2)
dimension waveb(LeNsp+4,3),sb(3)
c
if (minm2.gt.1-mbc .or. maxm2 .lt. maxm+mbc) then
write(6,*) 'need to increase maxm2 in rpB'
stop
endif
c
if (ixy.eq.1) then
mu = Nsp+1
mv = Nsp+2
else
mu = Nsp+2
mv = Nsp+1
endif
mE = Nsp+3
mT = Nsp+4
c
c # compute the flux differences bmasdq and bpasdq
c
do 20 i=2-mbc,mx+mbc
dpdr = 0.d0
drho = 0.d0
do k = 1, Nsp
drho = drho + asdq(i,k)
dpdr = dpdr + pk(k,i) * asdq(i,k)
enddo
c
a2 = g1a2(i)*(dpdr - ( u(i)*asdq(i,mu) + v(i)*asdq(i,mv) )
& + asdq(i,mE) )
a3 = asdq(i,mu) - u(i)*drho
a4 = 0.5d0*( a2 - ( v(i)*drho - asdq(i,mv) )/a(i) )
a1 = a2 - a4
c
c # Compute the waves.
c # Note that the 1+k-waves, for 1 .le. k .le. Nsp travel at
c # the same speed and are lumped together in wave(.,.,2).
c # The 3-wave is then stored in wave(.,.,3).
c
do k = 1, Nsp
c # 1-wave
waveb(k,1) = a1*Y(k,i)
c # 2-wave
waveb(k,2) = asdq(i,k) - Y(k,i)*a2
c # 3-wave
waveb(k,3) = a4*Y(k,i)
enddo
c # 1-wave
waveb(mu,1) = a1*u(i)
waveb(mv,1) = a1*(v(i) - a(i))
waveb(mE,1) = a1*(enth(i) - v(i)*a(i))
waveb(mT,1) = 0.d0
sb(1) = v(i)-a(i)
c
c # 2-wave
waveb(mu,2) = (drho - a2)*u(i) + a3
waveb(mv,2) = (drho - a2)*v(i)
waveb(mE,2) = (drho - a2)*u2v2(i)
& - dpdr + dpY(i)*a2 + a3*u(i)
waveb(mT,2) = 0.d0
sb(2) = v(i)
c
c # 3-wave
waveb(mu,3) = a4*u(i)
waveb(mv,3) = a4*(v(i) + a(i))
waveb(mE,3) = a4*(enth(i) + v(i)*a(i))
waveb(mT,3) = 0.d0
sb(3) = v(i)+a(i)
c
do 10 m=1,meqn
bmasdq(i,m) = 0.d0
bpasdq(i,m) = 0.d0
do 10 mw=1,mwaves
bmasdq(i,m) = bmasdq(i,m)
& + dmin1(sb(mw), 0.d0) * waveb(m,mw)
bpasdq(i,m) = bpasdq(i,m)
& + dmax1(sb(mw), 0.d0) * waveb(m,mw)
10 continue
c
20 continue
c
return
end