c
c =========================================================
subroutine rp1eum(maxmx,meqn,mwaves,mbc,mx,ql,qr,maux,
& auxl,auxr,wave,s,amdq,apdq)
c =========================================================
c
c # solve Riemann problems for the 1D two-component Euler equations
c # using HLLC
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,
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 routine step1, rp is called with ql = qr = q.
c
c # Copyright (C) 2003-2007 California Institute of Technology
c # Ralf Deiterding, ralf@amroc.net
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 amdq(1-mbc:maxmx+mbc, meqn)
dimension apdq(1-mbc:maxmx+mbc, meqn)
c
c # local storage
c ---------------
parameter (max2 = 100002) !# assumes at most 100000 grid points with mbc=2
dimension qls(3), qrs(3)
c
c # Riemann solver returns flux differences
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 0
c
do 10 i = 2-mbc, mx+mbc
c
rl = qr(i-1,1)
ul = qr(i-1,2)/rl
pl = (qr(i-1,3) - 0.5d0*qr(i-1,2)**2/rl
& - qr(i-1,5) ) / qr(i-1,4)
c
rr = ql(i ,1)
ur = ql(i ,2)/rr
pr = (ql(i ,3) - 0.5d0*ql(i ,2)**2/rr
& - ql(i ,5) ) / ql(i ,4)
c
rhsqrtl = dsqrt(qr(i-1,1))
rhsqrtr = dsqrt(ql(i,1))
rhsq2 = rhsqrtl + rhsqrtr
gamma1 = rhsq2 / ( qr(i-1,4)*rhsqrtl + ql(i,4)*rhsqrtr )
ui = (qr(i-1,2)/rhsqrtl + ql(i,2)/rhsqrtr) / rhsq2
enthi = (((qr(i-1,3)+pl)/rhsqrtl
& + (ql(i ,3)+pr)/rhsqrtr)) / rhsq2
ai = dsqrt(gamma1*(enthi - .5d0*ui**2))
c
sl = ui-ai
sr = ui+ai
ss = (pr-pl+rl*ul*(sl-ul)-rr*ur*(sr-ur))/
& (rl*(sl-ul)-rr*(sr-ur))
c
qrs(1) = rr*(sr-ur)/(sr-ss)
qrs(2) = qrs(1)*ss
qrs(3) = qrs(1)*(ql(i ,3)/rr+
& (ss-ur)*(ss+pr/(rr*(sr-ur))))
c
qls(1) = rl*(sl-ul)/(sl-ss)
qls(2) = qls(1)*ss
qls(3) = qls(1)*(qr(i-1,3)/rl+
& (ss-ul)*(ss+pl/(rl*(sl-ul))))
c
do m=1,3
wave(i,m,1) = qls(m) - qr(i-1,m)
wave(i,m,2) = qrs(m) - qls(m)
wave(i,m,3) = ql(i,m) - qrs(m)
enddo
do m=4,5
wave(i,m,1) = 0.d0
wave(i,m,2) = ql(i,m) - qr(i-1,m)
wave(i,m,3) = 0.d0
enddo
c
s(i,1) = sl
s(i,2) = ss
s(i,3) = sr
c
do m=1,meqn
amdq(i,m) = 0.d0
apdq(i,m) = 0.d0
do mw=1,mwaves
if (s(i,mw) .lt. 0.d0) then
amdq(i,m) = amdq(i,m) + s(i,mw)*wave(i,m,mw)
else
apdq(i,m) = apdq(i,m) + s(i,mw)*wave(i,m,mw)
endif
enddo
enddo
10 continue
return
end
c