c
c
c =====================================================
subroutine rpn3eu(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,
& maux,auxl,auxr,wave,s,fl,fr)
c =====================================================
c
c # FORCE scheme for the 3D Euler equations. The flux of the FORCE
c # scheme is the arithmetic mean of the fluxes of the finite difference
c # schemes of Richtmyer and Lax-Friedrichs. Use parameters
c # richtmyer, laxfriedrich to switch to the original schemes.
c
c # Eleuterio F. Toro, "Riemann solvers and numerical methods
c # for fluid dynamics", Springer-Verlag, Berlin 1997.
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 routine step1, 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)
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)
common /param/ gamma,gamma1
include "call.i"
c
c # local storage
c ---------------
parameter (maxmrp = 1005) !# assumes at most 1000 grid points with mbc=5
parameter (minmrp = -4) !# assumes at most mbc=5
dimension qint(minmrp:maxmrp2,5), fint(minmrp:maxmrp,5),
& auxint(minmrp:maxmrp,0,3)
logical richtmyer, laxfriedrich
c
data richtmyer /.true./
data laxfriedrich /.true./
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 momentum:
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
dxdt = 0.5d0*dxcom/dtcom
dtdx = 0.5d0*dtcom/dxcom
c
call flx3(ixyz,maxm,meqn,mbc,mx,ql,maux,auxl,fl)
call flx3(ixyz,maxm,meqn,mbc,mx,qr,maux,auxr,fr)
c
do 50 i = 2-mbc, mx+mbc
do 50 m=1,meqn
qint(i,m) = 0.5d0*(qr(i-1,m) + ql(i,m)) +
& dtdx*(fr(i-1,m) - fl(i,m))
50 continue
do 60 i = 2-mbc, mx+mbc
do 60 m=1,maux
auxint(i,m,2) = 0.5d0*(auxl(i,m,2) + auxr(i,m,2))
60 continue
call flx3(ixyz,max2,meqn,mbc,mx,qint,maux,auxint,fint)
c
do 100 i = 2-mbc, mx+mbc
ul = 0.5d0*qr(i-1,mu)/qr(i-1,1)
ur = 0.5d0*ql(i ,mu)/ql(i ,1)
pl = gamma1*(qr(i-1,5) - 0.5d0*(qr(i-1,mu)**2+
& qr(i-1,mv)**2+qr(i-1,mw)**2)/qr(i-1,1))
pr = gamma1*(ql(i ,5) - 0.5d0*(ql(i ,mu)**2+
& ql(i ,mv)**2+ql(i ,mw)**2)/ql(i ,1))
al = dsqrt(gamma*pl/qr(i-1,1))
ar = dsqrt(gamma*pr/ql(i ,1))
s(i,1) = dmax1(dabs(ul-al),dabs(ur-ar))
s(i,2) = dmax1(dabs(ul ),dabs(ur ))
s(i,3) = dmax1(dabs(ul+al),dabs(ur+ar))
do 110 mws=1,mwaves
do 110 m=1,meqn
wave(i,m,mws) = 0.d0
110 continue
do 100 m=1,meqn
if (richtmyer)
& fl(i,m) = fint(i,m)
if (laxfriedrich)
& fl(i,m) = dxdt*(qr(i-1,m) - ql(i,m)) +
& 0.5d0*(fr(i-1,m) + fl(i,m))
if (richtmyer.and.laxfriedrich)
& fl(i,m) = 0.5d0*(fl(i,m) + fint(i,m))
100 continue
c
do 120 i = 2-mbc, mx+mbc
do 120 m=1,meqn
fr(i,m) = -fl(i,m)
120 continue
c
return
end
c