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 # Steger & Warming - 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 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,l,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 el(3), er(3)
common /param/ gamma,gamma1
c
c # Method returns fluxes
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 1
c
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
pl = gamma1*(qr(i-1,3) - 0.5d0*ul**2*rhol)
pr = gamma1*(ql(i ,3) - 0.5d0*ur**2*rhor)
Hl = (qr(i-1,3)+pl)/rhol
Hr = (ql(i ,3)+pr)/rhor
c
al2 = gamma*pl/rhol
al = dsqrt(al2)
ar2 = gamma*pr/rhor
ar = dsqrt(ar2)
c
el(1) = 0.5d0*(ul-al + dabs(ul-al))
el(2) = 0.5d0*(ul + dabs(ul) )
el(3) = 0.5d0*(ul+al + dabs(ul+al))
er(1) = 0.5d0*(ur-ar - dabs(ur-ar))
er(2) = 0.5d0*(ur - dabs(ur) )
er(3) = 0.5d0*(ur+ar - dabs(ur+ar))
c
facl = 0.5d0*rhol/gamma
facr = 0.5d0*rhor/gamma
c
taul = facl*(el(1) + 2.d0*gamma1*el(2) + el(3))
taur = facr*(er(1) + 2.d0*gamma1*er(2) + er(3))
zetal = al*facl*(el(1)-el(3))
zetar = ar*facr*(er(1)-er(3))
c
fl(i,1) = taul + taur
fl(i,2) = ul*taul - zetal + ur*taur - zetar
fl(i,3) = Hl*taul - ul*zetal - 2.d0*el(2)*facl*al2 +
& Hr*taur - ur*zetar - 2.d0*er(2)*facr*ar2
c
do 20 m = 1, meqn
fr(i,m) = -fl(i,m)
20 continue
c
do 10 mw=1,mwaves
s(i,mw) = dmax1(dabs(el(mw)),dabs(er(mw)))
do 10 m=1,meqn
wave(i,m,mw) = 0.d0
10 continue
c
return
end
c