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 van Leer's Flux-Vector-Splitting
c # following Shuen's approach
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)
dimension fvl(LeNsp+4), fvr(LeNsp+4), rkl(LeNsp), rkr(LeNsp)
double precision Ml, Mr
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
rkl(k) = qr(i-1,k)
rkr(k) = ql(i ,k)
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
rhoul = qr(i-1,mu)
rhour = ql(i ,mu)
ul = rhoul/rhol
ur = rhour/rhor
vl = qr(i-1,mv)/rhol
vr = ql(i ,mv)/rhor
c
c # left/right Temperatures already calculated
c
rhoel = qr(i-1,mE)-0.5d0*rhol*(ul**2+vl**2)
call SolveTrhok(qr(i-1,mT),rhoel,rhoWl,rkl,Nsp,ifail)
if (ifail.ne.0) write(6,600) i-1,qr(i-1,mT)
rhoer = ql(i ,mE)-0.5d0*rhor*(ur**2+vr**2)
call SolveTrhok(ql(i ,mT),rhoer,rhoWr,rkr,Nsp,ifail)
if (ifail.ne.0) write(6,601) i ,ql(i ,mT)
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
Ml = ul/al
Mr = ur/ar
c
el(1) = ul-al
el(2) = ul
el(3) = ul+al
er(1) = ur-ar
er(2) = ur
er(3) = ur+ar
c
if (Ml.gt.1.d0) then
do k = 1, Nsp
fvl(k) = Yl(k)*rhoul
enddo
fvl(mu) = rhoul*ul+pl
fvl(mv) = rhoul*vl
fvl(mE) = rhoul*Hl
else if (Ml.lt.-1.d0) then
do m = 1, Nsp+3
fvl(m) = 0.d0
enddo
else
fl = 0.25d0*rhol*al*(Ml+1.d0)**2
fhl = avgtabip(Tl,Yl,hms,Nsp)/al2
xl = fhl/(1.d0+2.d0*fhl)
do k = 1, Nsp
fvl(k) = fl*Yl(k)
enddo
fvl(mu) = fl*(ul-(ul-2.d0*al)/gammal)
fvl(mv) = fl* vl
fvl(mE) = fl*(Hl-xl*(ul-al)**2)
endif
fvl(mT) = 0.d0
c
if (Mr.lt.-1.d0) then
do k = 1, Nsp
fvr(k) = Yr(k)*rhour
enddo
fvr(mu) = rhour*ur+pr
fvr(mv) = rhour*vr
fvr(mE) = rhour*Hr
else if (Mr.gt.1.d0) then
do m = 1, Nsp+3
fvr(m) = 0.d0
enddo
else
fr = -0.25d0*rhor*ar*(Mr-1.d0)**2
fhr = avgtabip(Tr,Yr,hms,Nsp)/ar2
xr = fhr/(1.d0+2.d0*fhr)
do k = 1, Nsp
fvr(k) = fr*Yr(k)
enddo
fvr(mu) = fr*(ur-(ur+2.d0*ar)/gammar)
fvr(mv) = fr* vr
fvr(mE) = fr*(Hr-xr*(ur+ar)**2)
endif
fvr(mT) = 0.d0
c
do 20 m = 1,meqn
amdq(i,m) = fvl(m) + fvr(m)
apdq(i,m) = -amdq(i,m)
20 continue
c
if (dabs(Ml).lt.1.d0) then
fl = (gammal+3.d0)/(2.d0*gammal+dabs(Ml)*(3.d0-gammal))
else
fl = 1.d0
endif
if (dabs(Mr).lt.1.d0) then
fr = (gammar+3.d0)/(2.d0*gammar+dabs(Mr)*(3.d0-gammar))
else
fr = 1.d0
endif
c
do 10 mw=1,mwaves
s(i,mw) = dmax1(dabs(fl*el(mw)),dabs(fr*er(mw)))
do 10 m=1,meqn
wave(i,m,mw) = 0.d0
10 continue
c
600 format('rpn2eurhok qr(',i2,'): T out of range using: ',f16.8)
601 format('rpn2eurhok ql(',i2,'): T out of range using: ',f16.8)
return
end
c