c
c
c =====================================================
subroutine rpn3euznd(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,
& maux,auxl,auxr,wave,s,fl,fr)
c =====================================================
c
c # solve Riemann problems for the 3D ZND-Euler equations using
c # the Flux-Vector-Splitting of Vijayasundaram
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 # 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 output, wave contains the waves, s the speeds,
c # fl and fr 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) 2002 Ralf Deiterding
c # Brandenburgische Universitaet Cottbus
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 fl(1-mbc:maxm+mbc, meqn)
dimension fr(1-mbc:maxm+mbc, meqn)
dimension auxl(1-mbc:maxm+mbc, maux, 3)
dimension auxr(1-mbc:maxm+mbc, maux, 3)
double precision el(3), er(3)
common /param/ gamma,gamma1,q0
c
c # Riemann solver returns flux differences
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 1
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 = 3
mv = 4
mw = 5
else if(ixyz .eq. 2)then
mu = 4
mv = 5
mw = 3
else
mu = 5
mv = 3
mw = 4
endif
c
do 10 i=2-mbc,mx+mbc
rho1l = qr(i-1,1)
rho1r = ql(i ,1)
rho2l = qr(i-1,2)
rho2r = ql(i ,2)
rhoul = qr(i-1,mu)
rhour = ql(i ,mu)
rhovl = qr(i-1,mv)
rhovr = ql(i ,mv)
rhowl = qr(i-1,mw)
rhowr = ql(i ,mw)
rhoEl = qr(i-1,6)
rhoEr = ql(i ,6)
rhol = rho1l+rho2l
rhor = rho1r+rho2r
c
rho = 0.5d0*(rhol + rhor )
rho1 = 0.5d0*(rho1l + rho1r)
rho2 = 0.5d0*(rho2l + rho2r)
rhou = 0.5d0*(rhoul + rhour)
rhov = 0.5d0*(rhovl + rhovr)
rhoE = 0.5d0*(rhoEl + rhoEr)
c
Y1 = rho1/rho
Y2 = rho2/rho
u = rhou/rho
v = rhov/rho
w = rhov/rho
p = gamma1*(rhoE - rho2*q0 - 0.5d0*rho*(u**2+v**2+w**2))
H = (rhoE+p)/rho
if (p.le.0.d0.or.rho.le.0.d0)
& write (6,*) 'Error in middle state in',i,p,pl,pr,
& rho,rhol,rhor,a,al,ar
a = dsqrt(gamma*p/rho)
f = 0.5d0/a**2
c
el1 = 0.5d0*(u-a + dabs(u-a))
el2 = 0.5d0*(u + dabs(u) )
el3 = 0.5d0*(u+a + dabs(u+a))
er1 = 0.5d0*(u-a - dabs(u-a))
er2 = 0.5d0*(u - dabs(u) )
er3 = 0.5d0*(u+a - dabs(u+a))
c
zl = el1-el3
ol = el1-2.d0*el2+el3
zr = er1-er3
or = er1-2.d0*er2+er3
dul = a*(rhol*u-rhoul)
dur = a*(rhor*u-rhour)
dEl = gamma1*(rhoEl-rho2l*q0+0.5d0*rhol*(u**2+v**2+w**2)-
& rhoul*u-rhovl*v-rhowl*w)
dEr = gamma1*(rhoEr-rho2r*q0+0.5d0*rhor*(u**2+v**2+w**2)-
& rhour*u-rhovr*v-rhowr*w)
f1 = f*(zl*dul + ol*dEl + zr*dur + or*dEr)
f2 = a*f*(ol*dul + zl*dEl + or*dur + zr*dEr)
c
fl(i,1) = rho1l*el2 + rho1r*er2 + Y1*f1
fl(i,2) = rho2l*el2 + rho2r*er2 + Y2*f1
fl(i,mu) = rhoul*el2 + rhour*er2 + u*f1 - f2
fl(i,mv) = rhovl*el2 + rhovr*er2 + v*f1
fl(i,mw) = rhowl*el2 + rhowr*er2 + w*f1
fl(i,6) = rhoEl*el2 + rhoEr*er2 + H*f1 - u*f2
c
do 20 m = 1,meqn
fr(i,m) = -fl(i,m)
20 continue
c
s(i,1) = u-a
s(i,2) = u
s(i,3) = u+a
do 10 mws=1,mwaves
do 10 m=1,meqn
wave(i,m,mws) = 0.d0
10 continue
c
return
end
c