c
c =========================================================
subroutine rpn2eurhok(ixy,maxm,meqn,mwaves,mbc,mx,ql,qr,maux,
& auxl,auxr,wave,s,amdq,apdq)
c =========================================================
c
c # solve Riemann problems for the 2D Euler equations of multiple
c # thermally perfect gases using the Flux-Vector-Splitting
c # of Steger & Warming
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 positive flux
c # apdq the negative flux
c # (the fluxes are stored twice to be consistent with the
c # flux-difference splitting formulation)
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:maxm+mbc, meqn)
dimension qr(1-mbc:maxm+mbc, meqn)
dimension s(1-mbc:maxm+mbc, mwaves)
dimension wave(1-mbc:maxm+mbc, meqn, mwaves)
dimension amdq(1-mbc:maxm+mbc, meqn)
dimension apdq(1-mbc:maxm+mbc, meqn)
c
c define local arrays
c
include "ck.i"
dimension Yl(LeNsp), Yr(LeNsp), el(3), er(3)
c
c # Riemann solver 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 to the orthogonal
c # momentum:
c
if (ixy.eq.1) then
mu = Nsp+1
mv = Nsp+2
else
mu = Nsp+2
mv = Nsp+1
endif
mE = Nsp+3
mT = Nsp+4
c
do 10 i=2-mbc,mx+mbc
rhol = 0.d0
rhor = 0.d0
rhoWl = 0.d0
rhoWr = 0.d0
do k = 1, Nsp
rhol = rhol + qr(i-1,k)
rhor = rhor + ql(i ,k)
rhoWl = rhoWl + qr(i-1,k)/Wk(k)
rhoWr = rhoWr + ql(i ,k)/Wk(k)
enddo
do k = 1, Nsp
Yl(k) = qr(i-1,k)/rhol
Yr(k) = ql(i ,k)/rhor
enddo
ul = qr(i-1,mu)/rhol
ur = ql(i ,mu)/rhor
vl = qr(i-1,mv)/rhol
vr = ql(i ,mv)/rhor
c
c # left/right Temperatures already calculated
c
Tl = qr(i-1,mT)
Tr = ql(i ,mT)
c
pl = rhoWl*RU*Tl
pr = rhoWr*RU*Tr
Hl = (qr(i-1,mE)+pl)/rhol
Hr = (ql(i ,mE)+pr)/rhor
c
c # Evaluate temperature depended gamma for left/right mixture
c
Cpl = avgtabip(Tl,Yl,cpk,Nsp)
gamma1l = RU / ( rhol/rhoWl*Cpl - RU )
gammal = gamma1l + 1.d0
Cpr = avgtabip(Tr,Yr,cpk,Nsp)
gamma1r = RU / ( rhor/rhoWr*Cpr - RU )
gammar = gamma1r + 1.d0
c
al2 = gammal*pl/rhol
al = dsqrt(al2)
ar2 = gammar*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
fl = 0.5d0*rhol/gammal
fr = 0.5d0*rhor/gammar
c
taul = fl*(el(1) + 2.d0*gamma1l*el(2) + el(3))
taur = fr*(er(1) + 2.d0*gamma1r*er(2) + er(3))
zetal = al*fl*(el(1)-el(3))
zetar = ar*fr*(er(1)-er(3))
c
do k = 1, Nsp
amdq(i,k) = Yl(k)*taul + Yr(k)*taur
enddo
amdq(i,mu) = ul*taul - zetal + ur*taur - zetar
amdq(i,mv) = vl*taul + vr*taur
amdq(i,mE) = Hl*taul - ul*zetal - 2.d0*el(2)*fl*al2 +
& Hr*taur - ur*zetar - 2.d0*er(2)*fr*ar2
amdq(i,mT) = 0.d0
c
do 20 m = 1, meqn
apdq(i,m) = -amdq(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