c
c =========================================================
subroutine rp1eu(maxmx,meqn,mwaves,mbc,mx,ql,qr,maux,
& auxl,auxr,wave,s,amdq,apdq)
c =========================================================
c
c # solve Riemann problems for the 1D Euler equations using Roe's
c # approximate Riemann solver.
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
c # On output, wave contains the waves,
c # s the speeds,
c # amdq the left-going flux difference A^- \Delta q
c # apdq the right-going flux difference A^+ \Delta q
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 routines, this routine is called with ql = qr
c
c Author: Ralf Deiterding
c
implicit double precision (a-h,o-z)
c
dimension wave(1-mbc:maxmx+mbc, meqn, mwaves)
dimension s(1-mbc:maxmx+mbc, mwaves)
dimension ql(1-mbc:maxmx+mbc, meqn)
dimension qr(1-mbc:maxmx+mbc, meqn)
dimension auxl(1-mbc:maxmx+mbc, maux)
dimension auxr(1-mbc:maxmx+mbc, maux)
dimension apdq(1-mbc:maxmx+mbc, meqn)
dimension amdq(1-mbc:maxmx+mbc, meqn)
c
c # local storage
c ---------------
parameter (max2 = 100002) !# assumes at most 100000 grid points with mbc=2
dimension u(-1:max2), enth(-1:max2), a(-1:max2), smax(-1:max2)
dimension delta(3), fl(-1:max2,3), fr(-1:max2,3)
logical efix, hll, roe, hllfix
common /param/ gamma,gamma1
c
data efix /.true./ !# use entropy fix for transonic rarefactions
data hll /.true./ !# use HLL solver if unphysical values occur
data roe /.true./ !# use Roe solver
c
c # Riemann solver returns flux differences
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 0
c
if (-1.gt.1-mbc .or. max2 .lt. maxmx+mbc) then
write(6,*) 'need to increase max2 in rp'
stop
endif
c
c # Compute Roe-averaged quantities:
c
do 20 i=2-mbc,mx+mbc
rhsqrtl = dsqrt(qr(i-1,1))
rhsqrtr = dsqrt(ql(i,1))
ul = qr(i-1,2)/qr(i-1,1)
ur = ql(i ,2)/ql(i ,1)
pl = gamma1*(qr(i-1,3) - 0.5d0*(qr(i-1,2)**2)/qr(i-1,1))
pr = gamma1*(ql(i ,3) - 0.5d0*(ql(i ,2)**2)/ql(i ,1))
al = dsqrt(gamma*pl/qr(i-1,1))
ar = dsqrt(gamma*pr/ql(i ,1))
rhsq2 = rhsqrtl + rhsqrtr
u(i) = (qr(i-1,2)/rhsqrtl + ql(i,2)/rhsqrtr) / rhsq2
enth(i) = (((qr(i-1,3)+pl)/rhsqrtl
& + (ql(i,3)+pr)/rhsqrtr)) / rhsq2
a2 = gamma1*(enth(i) - .5d0*u(i)**2)
a(i) = dsqrt(a2)
smax(i) = dmax1(dmax1(dabs(ur-ar-(ul-al)),dabs(ur-ul)),
& dabs(ur+ar-(ul+al)))
c
20 continue
c
do 30 i=2-mbc,mx+mbc
c
c # find a1 thru a3, the coefficients of the 3 eigenvectors:
c
delta(1) = ql(i,1) - qr(i-1,1)
delta(2) = ql(i,2) - qr(i-1,2)
delta(3) = ql(i,3) - qr(i-1,3)
a2 = gamma1/a(i)**2 * ((enth(i)-u(i)**2)*delta(1)
& + u(i)*delta(2) - delta(3))
a3 = (delta(2) + (a(i)-u(i))*delta(1) - a(i)*a2) / (2.d0*a(i))
a1 = delta(1) - a2 - a3
c
c # Compute the waves.
c
wave(i,1,1) = a1
wave(i,2,1) = a1*(u(i)-a(i))
wave(i,3,1) = a1*(enth(i) - u(i)*a(i))
s(i,1) = u(i)-a(i)
c
wave(i,1,2) = a2
wave(i,2,2) = a2*u(i)
wave(i,3,2) = a2*0.5d0*u(i)**2
s(i,2) = u(i)
c
wave(i,1,3) = a3
wave(i,2,3) = a3*(u(i)+a(i))
wave(i,3,3) = a3*(enth(i)+u(i)*a(i))
s(i,3) = u(i)+a(i)
30 continue
c
c # compute flux differences as
c # (+/-)
c # A (Ur-Ul) = 0.5*( f(Ur)-f(Ul) +/- |A|(Ur-Ul) )
c --------------------------
c
call flx1(maxmx,meqn,mbc,mx,qr,maux,auxr,apdq)
call flx1(maxmx,meqn,mbc,mx,ql,maux,auxl,amdq)
c
do 35 i = 1-mbc, mx+mbc
do 35 m=1,meqn
fl(i,m) = amdq(i,m)
fr(i,m) = apdq(i,m)
35 continue
c
if (roe) then
do 40 i = 2-mbc, mx+mbc
do 40 m=1,meqn
amdq(i,m) = 0.5d0*(fl(i,m)-fr(i-1,m))
40 continue
c
do 50 i = 2-mbc, mx+mbc
do 50 m=1,meqn
sw = 0.d0
do 60 mw=1,mwaves
sl = dabs(s(i,mw))
if (efix.and.dabs(s(i,mw)).lt.smax(i).and.mw.ne.2)
& sl = s(i,mw)**2/(2.d0*smax(i))+0.5d0*smax(i)
sw = sw + sl*wave(i,m,mw)
60 continue
amdq(i,m) = amdq(i,m) - 0.5d0*sw
apdq(i,m) = amdq(i,m) + sw
50 continue
endif
c
if (hll) then
do 55 i = 2-mbc, mx+mbc
hllfix = .false.
if (.not.roe) hllfix = .true.
c
rhol = qr(i-1,1) + wave(i,1,1)
rhoul = qr(i-1,2) + wave(i,2,1)
El = qr(i-1,3) + wave(i,3,1)
pl = gamma1*(El - 0.5d0*rhoul**2/rhol)
if (rhol.le.0.d0.or.pl.le.0.d0) hllfix = .true.
c
rhor = ql(i,1) - wave(i,1,3)
rhour = ql(i,2) - wave(i,2,3)
Er = ql(i,3) - wave(i,3,3)
pr = gamma1*(Er - 0.5d0*rhour**2/rhor)
if (rhor.le.0.d0.or.pr.le.0.d0) hllfix = .true.
c
if (hllfix) then
c if (roe) write (6,*) 'Switching to HLL in',i
c
rl = qr(i-1,1)
ul = qr(i-1,2)/rl
pl = gamma1*(qr(i-1,3) - 0.5d0*qr(i-1,2)**2/rl)
al = dsqrt(gamma*pl/rl)
c
rr = ql(i ,1)
ur = ql(i ,2)/rr
pr = gamma1*(ql(i ,3) - 0.5d0*ql(i ,2)**2/rr)
ar = dsqrt(gamma*pr/rr)
c
sl = dmin1(ul-al,ur-ar)
sr = dmax1(ul+al,ur+ar)
c
do m=1,meqn
if (sl.ge.0.d0) fg = fr(i-1,m)
if (sr.le.0.d0) fg = fl(i,m)
if (sl.lt.0.d0.and.sr.gt.0.d0)
& fg = (sr*fr(i-1,m) - sl*fl(i,m) +
& sl*sr*(ql(i,m)-qr(i-1,m)))/ (sr-sl)
amdq(i,m) = fg-fr(i-1,m)
apdq(i,m) = -(fg-fl(i ,m))
enddo
s(i,1) = sl
s(i,2) = 0.d0
s(i,3) = sr
endif
55 continue
endif
c
return
end
c