c
c =========================================================
subroutine rp1eurhok(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 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: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 define local arrays
c
include "ck.i"
dimension Yl(LeNsp), Yr(LeNsp), el(3), er(3)
dimension fvl(LeNsp+3), fvr(LeNsp+3)
double precision Ml, Mr
c
c # Riemann solver returns fluxes
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 1
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
rhoul = qr(i-1,Nsp+1)
rhour = ql(i ,Nsp+1)
ul = rhoul/rhol
ur = rhour/rhor
c
c # left/right Temperatures already calculated
c
Tl = qr(i-1,Nsp+3)
Tr = ql(i ,Nsp+3)
c
pl = rhoWl*RU*Tl
pr = rhoWr*RU*Tr
Hl = (qr(i-1,Nsp+2)+pl)/rhol
Hr = (ql(i ,Nsp+2)+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(Nsp+1) = rhoul*ul+pl
fvl(Nsp+2) = rhoul*Hl
else if (Ml.lt.-1.d0) then
do m = 1, Nsp+2
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) = Yl(k)*fl
enddo
fvl(Nsp+1) = fl*(ul-(ul-2.d0*al)/gammal)
fvl(Nsp+2) = fl*(Hl-xl*(ul-al)**2)
endif
fvl(Nsp+3) = 0.d0
c
if (Mr.lt.-1.d0) then
do k = 1, Nsp
fvr(k) = Yr(k)*rhour
enddo
fvr(Nsp+1) = rhour*ur+pr
fvr(Nsp+2) = rhour*Hr
else if (Mr.gt.1.d0) then
do m = 1, Nsp+2
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) = Yr(k)*fr
enddo
fvr(Nsp+1) = fr*(ur-(ur+2.d0*ar)/gammar)
fvr(Nsp+2) = fr*(Hr-xr*(ur+ar)**2)
endif
fvr(Nsp+3) = 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
return
end
c