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 # van Leer's Flux Vector Splitting following Shuen's approach for
c # multicomponent gas mixtures
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 Ml, Mr
dimension fvl(6), fvr(6), sl(3), sr(3)
common /param/ gamma,gamma1,q0
c
c # Method returns fluxes
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
c # Van Leer's Flux Vector Splitting
c
gamma2 = gamma**2-1
do 10 i=2-mbc,mx+mbc
rhol = qr(i-1,1)+qr(i-1,2)
rhor = ql(i ,1)+ql(i ,2)
Y1l = qr(i-1,1)/rhol
Y2l = qr(i-1,2)/rhol
Y1r = ql(i ,1)/rhor
Y2r = ql(i ,2)/rhor
ul = qr(i-1,mu)/rhol
ur = ql(i ,mu)/rhor
vl = qr(i-1,mv)/rhol
vr = ql(i ,mv)/rhor
wl = qr(i-1,mw)/rhol
wr = ql(i ,mw)/rhor
El = qr(i-1,6)/rhol
Er = ql(i ,6)/rhor
pl = gamma1*(qr(i-1,6) - qr(i-1,2)*q0 -
& 0.5d0*(qr(i-1,mu)**2+qr(i-1,mv)**2+qr(i-1,mw)**2)/rhol)
pr = gamma1*(ql(i ,6) - ql(i ,2)*q0 -
& 0.5d0*(ql(i ,mu)**2+ql(i ,mv)**2+ql(i ,mw)**2)/rhor)
Hl = El+pl/rhol
Hr = Er+pr/rhor
c
al2 = gamma*pl/rhol
al = dsqrt(al2)
ar2 = gamma*pr/rhor
ar = dsqrt(ar2)
c
Ml = ul/al
Mr = ur/ar
c
sl(1) = ul-al
sl(2) = ul
sl(3) = ul+al
sr(1) = ur-ar
sr(2) = ur
sr(3) = ur+ar
c
if (Ml.ge.1.d0) then
facl = rhol*ul
fvl(1) = facl*Y1l
fvl(2) = facl*Y2l
fvl(mu) = facl*ul+pl
fvl(mv) = facl*vl
fvl(mw) = facl*wl
fvl(6) = facl*El+ul*pl
else if (Ml.le.-1.d0) then
do m = 1,meqn
fvl(m) = 0.d0
enddo
else
fhl = gamma*(El - Y2l*q0 - 0.5d0*(ul**2+vl**2+wl**2))/al2
xl = fhl/(1.d0+2.d0*fhl)
facl = 0.25d0*rhol*al*(Ml+1.d0)**2
fvl(1) = facl*Y1l
fvl(2) = facl*Y2l
fvl(mu) = facl*2.d0*al/gamma*(0.5d0*gamma1*Ml+1.d0)
fvl(mv) = facl*vl
fvl(mw) = facl*wl
fvl(6) = facl*(Hl-xl*(ul-al)**2)
endif
c
if (Mr.le.-1.d0) then
facr = rhor*ur
fvr(1) = facr*Y1r
fvr(2) = facr*Y2r
fvr(mu) = facr*ur+pr
fvr(mv) = facr*vr
fvr(mw) = facr*wr
fvr(6) = facr*Er+ur*pr
else if (Mr.ge.1.d0) then
do m = 1,meqn
fvr(m) = 0.d0
enddo
else
fhr = gamma*(Er - Y2r*q0 - 0.5d0*(ur**2+vr**2+wr**2))/ar2
xr = fhr/(1.d0+2.d0*fhr)
facr = -0.25d0*rhor*ar*(Mr-1.d0)**2
fvr(1) = facr*Y1r
fvr(2) = facr*Y2r
fvr(mu) = facr*2.d0*ar/gamma*(0.5d0*gamma1*Mr-1.d0)
fvr(mv) = facr*vr
fvr(mw) = facr*wr
fvr(6) = facr*(Hr-xr*(ur+ar)**2)
endif
c
do 20 m = 1,meqn
fl(i,m) = fvl(m) + fvr(m)
fr(i,m) = -fl(i,m)
20 continue
c
if (dabs(Ml).lt.1.d0) then
facl = (gamma+3.d0)/(2.d0*gamma+dabs(Ml)*(3.d0-gamma))
else
facl = 1.d0
endif
if (dabs(Mr).lt.1.d0) then
facr = (gamma+3.d0)/(2.d0*gamma+dabs(Mr)*(3.d0-gamma))
else
facr = 1.d0
endif
c
do 10 mws=1,mwaves
s(i,mws) = dmax1(dabs(facl*sl(mws)),dabs(facr*sr(mws)))
do 10 m=1,meqn
wave(i,m,mws) = 0.d0
10 continue
c
return
end
c