c
c
c ==================================================================
subroutine rpt3eurhok(ixyz,icoor,maxm,meqn,mwaves,mbc,mx,
& ql,qr,maux,aux1,aux2,aux3,ilr,asdq,
& bmasdq,bpasdq)
c ==================================================================
c
c # solve Riemann problems in the transverse direction for the 3D Euler
c # equations of multiple thermally perfect gases using Roe's approximate
c # Riemann solver.
c #
c # On input,
c
c # ql,qr is the data along some one-dimensional slice, as in rpn3
c # This slice is
c # in the x-direction if ixyz=1,
c # in the y-direction if ixyz=2, or
c # in the z-direction if ixyz=3.
c # asdq is an array of flux differences (A^* \Delta q).
c # asdq(i,:) is the flux difference propagating away from
c # the interface between cells i-1 and i.
c # imp = 1 if asdq = A^- \Delta q, the left-going flux difference
c # 2 if asdq = A^+ \Delta q, the right-going flux difference
c
c # aux2 is the auxiliary array (if method(7)=maux>0) along
c # the plane where this slice lies, say at j=J if ixyz=1.
c # aux2(:,:,1) contains data along j=J, k=k-1
c # aux2(:,:,2) contains data along j=J, k=k
c # aux2(:,:,3) contains data along j=J, k=k+1
c # aux1 is the auxiliary array along the plane with j=J-1
c # aux3 is the auxiliary array along the plane with j=J+1
c
c # if ixyz=2 then aux2 is in some plane k=K, and
c # aux2(:,:,1) contains data along i=I-1, k=K, etc.
c
c # if ixyz=3 then aux2 is in some plane i=I, and
c # aux2(:,:,1) contains data along j=j-1, i=I, etc.
c
c # On output,
c # If data is in x-direction (ixyz=1) then this routine does the
c # splitting of asdq (= A^* \Delta q, where * = + or -) ...
c
c # into down-going flux difference bmasdq (= B^- A^* \Delta q)
c # and up-going flux difference bpasdq (= B^+ A^* \Delta q)
c # when icoor = 2,
c
c # or
c
c # into down-going flux difference bmasdq (= C^- A^* \Delta q)
c # and up-going flux difference bpasdq (= C^+ A^* \Delta q)
c # when icoor = 3.
c #
c
c # More generally, ixyz specifies what direction the slice of data is
c # in, and icoor tells which transverse direction to do the splitting in:
c
c # If ixyz = 1, data is in x-direction and then
c # icoor = 2 => split in the y-direction
c # icoor = 3 => split in the z-direction
c
c # If ixyz = 2, data is in y-direction and then
c # icoor = 2 => split in the z-direction
c # icoor = 3 => split in the x-direction
c
c # If ixyz = 3, data is in z-direction and then
c # icoor = 2 => split in the x-direction
c # icoor = 3 => split in the y-direction
c
c #
c # Uses Roe averages and other quantities which were
c # computed in rpn3eurhok and stored in the common block comroe.
c
c # Copyright (C) 2002 Ralf Deiterding
c # Brandenburgische Universitaet Cottbus
c
implicit double precision (a-h,o-z)
c
include "ck.i"
c
dimension ql(1-mbc:maxm+mbc, meqn)
dimension qr(1-mbc:maxm+mbc, meqn)
dimension asdq(1-mbc:maxm+mbc, meqn)
dimension bmasdq(1-mbc:maxm+mbc, meqn)
dimension bpasdq(1-mbc:maxm+mbc, meqn)
dimension aux1(1-mbc:maxm+mbc, maux, 3)
dimension aux2(1-mbc:maxm+mbc, maux, 3)
dimension aux3(1-mbc:maxm+mbc, maux, 3)
c
c local arrays -- common block comroe is passed to rpt3eurhok
c ------------
parameter (maxmrp = 1005) !# assumes atmost max(mx,my,mz) = 1000 with mbc=5
parameter (minmrp = -4) !# assumes at most mbc=5
common /comroe/ u(minmrp:maxmrp), v(minmrp:maxmrp),
& w(minmrp:maxmrp), u2v2w2(minmrp:maxmrp),
& enth(minmrp:maxmrp), a(minmrp:maxmrp),
& g1a2(minmrp:maxmrp), dpY(minmrp:maxmrp),
& Y(LeNsp,minmrp:maxmrp), pk(LeNsp,-1:maxmrp)
c define local arrays
dimension waveb(LeNsp+5,3),sb(3)
dimension rkl(LeNsp), rkr(LeNsp)
dimension hkl(LeNsp), hkr(LeNsp)
c
if (minmrp.gt.1-mbc .or. maxmrp .lt. maxm+mbc) then
write(6,*) 'need to increase maxmrp in rp3t'
stop
endif
c
if(ixyz .eq. 1)then
mu = Nsp+1
mv = Nsp+2
mw = Nsp+3
else if(ixyz .eq. 2)then
mu = Nsp+2
mv = Nsp+3
mw = Nsp+1
else
mu = Nsp+3
mv = Nsp+1
mw = Nsp+2
endif
mE = Nsp+4
mT = Nsp+5
c
c # Solve Riemann problem in the second coordinate direction
c
if( icoor .eq. 2 )then
do 10 i=2-mbc,mx+mbc
dpdr = 0.d0
drho = 0.d0
do k = 1, Nsp
drho = drho + asdq(i,k)
dpdr = dpdr + pk(k,i) * asdq(i,k)
enddo
c
a2 = g1a2(i)*(dpdr - ( u(i)*asdq(i,mu) + v(i)*asdq(i,mv) +
& w(i)*asdq(i,mw) ) + asdq(i,mE) )
a3 = asdq(i,mu) - u(i)*drho
a4 = asdq(i,mw) - w(i)*drho
a5 = 0.5d0*( a2 - ( v(i)*drho - asdq(i,mv) )/a(i) )
a1 = a2 - a5
c
c # Compute the waves.
c # Note that the 1+k-waves, for 1 .le. k .le. Nsp travel at
c # the same speed and are lumped together in wave(.,.,2).
c # The 3-wave is then stored in wave(.,.,3).
c
do k = 1, Nsp
c # 1-wave
waveb(k,1) = a1*Y(k,i)
c # 2-wave
waveb(k,2) = asdq(i,k) - Y(k,i)*a2
c # 3-wave
waveb(k,3) = a5*Y(k,i)
enddo
c # 1-wave
waveb(mu,1) = a1*u(i)
waveb(mv,1) = a1*(v(i) - a(i))
waveb(mw,1) = a1*w(i)
waveb(mE,1) = a1*(enth(i) - v(i)*a(i))
waveb(mT,1) = 0.d0
sb(1) = v(i)-a(i)
c
c # 2-wave
waveb(mu,2) = (drho - a2)*u(i) + a3
waveb(mv,2) = (drho - a2)*v(i)
waveb(mw,2) = (drho - a2)*w(i) + a4
waveb(mE,2) = (drho - a2)*u2v2w2(i)
& - dpdr + dpY(i)*a2 + a3*u(i) + a4*w(i)
waveb(mT,2) = 0.d0
sb(2) = v(i)
c
c # 3-wave
waveb(mu,3) = a5*u(i)
waveb(mv,3) = a5*(v(i) + a(i))
waveb(mw,3) = a5*w(i)
waveb(mE,3) = a5*(enth(i) + v(i)*a(i))
waveb(mT,3) = 0.d0
sb(3) = v(i)+a(i)
c
do 20 m=1,meqn
bmasdq(i,m) = 0.d0
bpasdq(i,m) = 0.d0
do 20 mws=1,mwaves
bmasdq(i,m) = bmasdq(i,m)
& + dmin1(sb(mws), 0.d0) * waveb(m,mws)
bpasdq(i,m) = bpasdq(i,m)
& + dmax1(sb(mws), 0.d0) * waveb(m,mws)
20 continue
c
10 continue
c
else
c
c # Solve Riemann problem in the third coordinate direction
c
do 30 i = 2-mbc, mx+mbc
dpdr = 0.d0
drho = 0.d0
do k = 1, Nsp
drho = drho + asdq(i,k)
dpdr = dpdr + pk(k,i) * asdq(i,k)
enddo
c
a2 = g1a2(i)*(dpdr - ( u(i)*asdq(i,mu) + v(i)*asdq(i,mv) +
& w(i)*asdq(i,mw) ) + asdq(i,mE) )
a3 = asdq(i,mu) - u(i)*drho
a4 = asdq(i,mv) - v(i)*drho
a5 = 0.5d0*( a2 - ( w(i)*drho - asdq(i,mw) )/a(i) )
a1 = a2 - a5
c
c # Compute the waves.
c # Note that the 1+k-waves, for 1 .le. k .le. Nsp travel at
c # the same speed and are lumped together in wave(.,.,2).
c # The 3-wave is then stored in wave(.,.,3).
c
do k = 1, Nsp
c # 1-wave
waveb(k,1) = a1*Y(k,i)
c # 2-wave
waveb(k,2) = asdq(i,k) - Y(k,i)*a2
c # 3-wave
waveb(k,3) = a5*Y(k,i)
enddo
c # 1-wave
waveb(mu,1) = a1*u(i)
waveb(mv,1) = a1*v(i)
waveb(mw,1) = a1*(w(i) - a(i))
waveb(mE,1) = a1*(enth(i) - w(i)*a(i))
waveb(mT,1) = 0.d0
sb(1) = w(i)-a(i)
c
c # 2-wave
waveb(mu,2) = (drho - a2)*u(i) + a3
waveb(mv,2) = (drho - a2)*v(i) + a4
waveb(mw,2) = (drho - a2)*w(i)
waveb(mE,2) = (drho - a2)*u2v2w2(i)
& - dpdr + dpY(i)*a2 + a3*u(i) + a4*v(i)
waveb(mT,2) = 0.d0
sb(2) = w(i)
c
c # 3-wave
waveb(mu,3) = a5*u(i)
waveb(mv,3) = a5*v(i)
waveb(mw,3) = a5*(w(i) + a(i))
waveb(mE,3) = a5*(enth(i) + w(i)*a(i))
waveb(mT,3) = 0.d0
sb(3) = w(i)+a(i)
c
do 40 m=1,meqn
bmasdq(i,m) = 0.d0
bpasdq(i,m) = 0.d0
do 40 mws=1,mwaves
bmasdq(i,m) = bmasdq(i,m)
& + dmin1(sb(mws), 0.d0) * waveb(m,mws)
bpasdq(i,m) = bpasdq(i,m)
& + dmax1(sb(mws), 0.d0) * waveb(m,mws)
40 continue
c
30 continue
c
endif
c
return
end