c
c
c =========================================================
subroutine rp1eu(maxmx,meqn,mwaves,mbc,mx,ql,qr,maux,
& auxl,auxr,wave,s,fl,fr)
c =========================================================
c
c # solve Riemann problems for the 1D Euler equations using
c # van Leer's Flux-Vector-Splitting
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 # 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 routine step1, rp is called with ql = qr = q.
c
c # Copyright (C) 2002 Ralf Deiterding
c # Brandenburgische Universitaet Cottbus
c
implicit double precision (a-h,o-z)
dimension ql(1-mbc:maxmx+mbc, meqn)
dimension qr(1-mbc:maxmx+mbc, meqn)
dimension s(1-mbc:maxmx+mbc, mwaves)
dimension wave(1-mbc:maxmx+mbc, meqn, mwaves)
dimension fl(1-mbc:maxmx+mbc, meqn)
dimension fr(1-mbc:maxmx+mbc, meqn)
double precision Ml, Mr, sl(3), sr(3), fvl(3), fvr(3)
common /param/ gamma,gamma1
c
c # Method returns fluxes
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 1
c
gamma2 = gamma**2-1
do 10 i=2-mbc,mx+mbc
rhol = qr(i-1,1)
rhor = ql(i ,1)
ul = qr(i-1,2)/rhol
ur = ql(i ,2)/rhor
El = qr(i-1,3)/rhol
Er = ql(i ,3)/rhor
pl = gamma1*(qr(i-1,3) - 0.5d0*(qr(i-1,2)**2)/rhol)
pr = gamma1*(ql(i ,3) - 0.5d0*(ql(i ,2)**2)/rhor)
al = dsqrt(gamma*pl/rhol)
ar = dsqrt(gamma*pr/rhor)
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.gt.1d0) then
fvl(1) = rhol*ul
fvl(2) = fvl(1)*ul+pl
fvl(3) = ul*(rhol*El+pl)
else if (Ml.lt.-1.d0) then
do m = 1,meqn
fvl(m) = 0.d0
enddo
else
fvl(1) = 0.25d0*rhol*al*(Ml+1.d0)**2
fvl(2) = fvl(1)*2.d0*al/gamma*(0.5d0*gamma1*Ml+1.d0)
fvl(3) = fvl(1)*2.d0*al**2/gamma2*(0.5d0*gamma1*Ml+1.d0)**2
endif
c
if (Mr.lt.-1.d0) then
fvr(1) = rhor*ur
fvr(2) = fvr(1)*ur+pr
fvr(3) = ur*(rhor*Er+pr)
else if (Mr.gt.1.d0) then
do m = 1,meqn
fvr(m) = 0.d0
enddo
else
fvr(1) = -0.25d0*rhor*ar*(Mr-1.d0)**2
fvr(2) = fvr(1)*2.d0*ar/gamma*(0.5d0*gamma1*Mr-1.d0)
fvr(3) = fvr(1)*2.d0*ar**2/gamma2*(0.5d0*gamma1*Mr-1.d0)**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 mw=1,mwaves
s(i,mw) = dmax1(dabs(facl*sl(mw)),dabs(facr*sr(mw)))
do 10 m=1,meqn
wave(i,m,mw) = 0.d0
10 continue
c
return
end
c