c
c
c ==================================================================
subroutine rpn3meu(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,
& maux,auxl,auxr,wave,s,amdq,apdq)
c ==================================================================
c
c # Solve Riemann problems for the 3D two-component Euler equations
c # using HLLC. Use flux difference splitting formulation for full
c # compatibility to Wave Propagation Method.
c
c # solve Riemann problems along one slice of data.
c # This data is along a slice in the x-direction if ixyz=1
c # the y-direction if ixyz=2.
c # the z-direction if ixyz=3.
c
c # On input, ql contains the state vector at the left edge of each cell
c # qr contains the state vector at the right edge of each cell
c # On output, wave contains the waves, s the speeds,
c # amdq and apdq the positive and negative flux.
c
c # Note that the i'th Riemann problem has left state qr(i-1,:)
c # and right state ql(i,:)
c # From the basic clawpack routines, this routine is called with ql = qr
c
c # Copyright (C) 2003-2007 California Institute of Technology
c # Ralf Deiterding, ralf@cacr.caltech.edu
c
implicit double precision(a-h,o-z)
c
dimension wave(1-mbc:maxm+mbc, meqn, mwaves)
dimension s(1-mbc:maxm+mbc, mwaves)
dimension ql(1-mbc:maxm+mbc, meqn)
dimension qr(1-mbc:maxm+mbc, meqn)
dimension amdq(1-mbc:maxm+mbc, meqn)
dimension apdq(1-mbc:maxm+mbc, meqn)
dimension auxl(1-mbc:maxm+mbc, maux, 3)
dimension auxr(1-mbc:maxm+mbc, maux, 3)
c
c local arrays -- common block comroe is passed to rpt3eum
c
c ------------
parameter (maxmrp = 1005) !# assumes atmost max(mx,my,mz) = 1000 with mbc=5
parameter (minmrp = -4) !# assumes at most mbc=5
dimension qls(5), qrs(5)
logical roespeed
c
common /comroe/ u2v2w2(minmrp:maxmrp),
& u(minmrp:maxmrp),v(minmrp:maxmrp),w(minmrp:maxmrp),
& enth(minmrp:maxmrp),a(minmrp:maxmrp),g1a2(minmrp:maxmrp),
& euv(minmrp:maxmrp),p(minmrp:maxmrp)
c
data roespeed /.false./ !# use Roe average for wave speed estimation
c
c # Riemann solver returns flux differences
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 0
c
if (minmrp.gt.1-mbc .or. maxmrp .lt. maxm+mbc) then
write(6,*) 'need to increase maxmrp in rpA'
stop
endif
c
c # set mu to point to the component of the system that corresponds
c # to momentum in the direction of this slice, mv and mw to the
c # orthogonal momentums:
c
if(ixyz .eq. 1)then
mu = 2
mv = 3
mw = 4
else if(ixyz .eq. 2)then
mu = 3
mv = 4
mw = 2
else
mu = 4
mv = 2
mw = 3
endif
c
c # note that notation for u,v, and w reflects assumption that the
c # Riemann problems are in the x-direction with u in the normal
c # direction and v and w in the orthogonal directions, but with the
c # above definitions of mu, mv, and mw the routine also works with
c # ixyz=2 and ixyz=3
c # and returns, for example, f0 as the Godunov flux g0 for the
c # Riemann problems u_t + g(u)_y = 0 in the y-direction.
c
c
c # Compute the Roe-averaged variables needed in the Roe solver.
c # These are stored in the common block comroe since they are
c # later used in routine rpt3eu to do the transverse wave
c # splitting.
c
do 10 i = 2-mbc, mx+mbc
if (qr(i-1,1).le.0.d0.or.ql(i,1).le.0.d0) then
write (6,*) 'Unrecoverable error in density',i
write (6,*) qr(i-1,1),ql(i,1)
stop
endif
c
rl = qr(i-1,1)
ul = qr(i-1,mu)/rl
vl = qr(i-1,mv)/rl
wl = qr(i-1,mw)/rl
El = qr(i-1,5)
gammal1 = 1.d0/qr(i-1,6)
gammal = gammal1 + 1.d0
pinfl = qr(i-1,7)*gammal1/gammal
pl = (El-0.5d0*(ul**2+vl**2+wl**2)*rl-qr(i-1,7))/qr(i-1,6)
if (pl+pinfl.le.0.d0.or.gammal.le.0.d0) then
write (6,*) 'Unrecoverable error in speed of sound l',i
write (6,*) pl,pinfl,pl+pinfl,gammal
stop
endif
al = dsqrt(gammal*(pl+pinfl)/rl)
c
rr = ql(i ,1)
ur = ql(i ,mu)/rr
vr = ql(i ,mv)/rr
wr = ql(i ,mw)/rr
Er = ql(i ,5)
gammar1 = 1.d0/ql(i ,6)
gammar = gammar1 + 1.d0
pinfr = ql(i ,7)*gammar1/gammar
pr = (Er-0.5d0*(ur**2+vr**2+wr**2)*rr-ql(i ,7))/ql(i ,6)
if (pr+pinfr.le.0.d0.or.gammar.le.0.d0) then
write (6,*) 'Unrecoverable error in speed of sound r',i
write (6,*) pr,pinfr,pr+pinfr,gammar
stop
endif
ar = dsqrt(gammar*(pr+pinfr)/rr)
c
rhsqrtl = dsqrt(qr(i-1,1))
rhsqrtr = dsqrt(ql(i,1))
rhsq2 = rhsqrtl + rhsqrtr
gamma1 = rhsq2 / ( qr(i-1,6)*rhsqrtl + ql(i,6)*rhsqrtr )
xjota = ( pl*qr(i-1,6)*rhsqrtl + pr*ql(i,6)*rhsqrtr ) / rhsq2
p(i) = xjota*gamma1
c
u(i) = (qr(i-1,mu)/rhsqrtl + ql(i,mu)/rhsqrtr) / rhsq2
v(i) = (qr(i-1,mv)/rhsqrtl + ql(i,mv)/rhsqrtr) / rhsq2
w(i) = (qr(i-1,mw)/rhsqrtl + ql(i,mw)/rhsqrtr) / rhsq2
enth(i) = (((qr(i-1,5)+pl)/rhsqrtl
& + (ql(i,5)+pr)/rhsqrtr)) / rhsq2
u2v2w2(i) = u(i)**2 + v(i)**2 + w(i)**2
a2 = gamma1*(enth(i) - .5d0*u2v2w2(i))
a(i) = dsqrt(a2)
if (a2.le.0.d0) a(i) = dmax1(al,ar)
g1a2(i) = gamma1 / a2
euv(i) = enth(i) - u2v2w2(i)
c
if (roespeed) then
sl = u(i)-a(i)
sr = u(i)+a(i)
else
sl = dmin1(ul-al,ur-ar)
sr = dmax1(ul+al,ur+ar)
endif
ss = (pr-pl+rl*ul*(sl-ul)-rr*ur*(sr-ur))/
& (rl*(sl-ul)-rr*(sr-ur))
c
qrs(1) = rr*(sr-ur)/(sr-ss)
qrs(mu) = qrs(1)*ss
qrs(mv) = qrs(1)*vr
qrs(mw) = qrs(1)*wr
qrs(5) = qrs(1)*(Er/rr+
& (ss-ur)*(ss+pr/(rr*(sr-ur))))
c
qls(1) = rl*(sl-ul)/(sl-ss)
qls(mu) = qls(1)*ss
qls(mv) = qls(1)*vl
qls(mw) = qls(1)*wl
qls(5) = qls(1)*(El/rl+
& (ss-ul)*(ss+pl/(rl*(sl-ul))))
c
do m=1,5
wave(i,m,1) = qls(m) - qr(i-1,m)
wave(i,m,2) = qrs(m) - qls(m)
wave(i,m,3) = ql(i,m) - qrs(m)
enddo
do m=6,7
wave(i,m,1) = 0.d0
wave(i,m,2) = ql(i,m) - qr(i-1,m)
wave(i,m,3) = 0.d0
enddo
c
s(i,1) = sl
s(i,2) = ss
s(i,3) = sr
c
do m=1,meqn
amdq(i,m) = 0.d0
apdq(i,m) = 0.d0
do mws=1,mwaves
if (s(i,mws) .lt. 0.d0) then
amdq(i,m) = amdq(i,m) + s(i,mws)*wave(i,m,mws)
else
apdq(i,m) = apdq(i,m) + s(i,mws)*wave(i,m,mws)
endif
enddo
enddo
10 continue
return
end