c
c
c ==================================================================
subroutine flux3(ixyz,maxm,meqn,maux,mwaves,mbc,mx,
& q1d,dtdx1d,dtdy,dtdz,aux1,aux2,aux3,
& method,mthlim,
& faddm,faddp,gaddm,gaddp,haddm,haddp,cfl1d,
& wave,s,amdq,apdq,cqxx,
& bmamdq,bmapdq,bpamdq,bpapdq,
& cmamdq,cmapdq,cpamdq,cpapdq,
& cmamdq2,cmapdq2,cpamdq2,cpapdq2,
& bmcqxx,bpcqxx,cmcqxx,cpcqxx,
& bmcmamdq,bmcmapdq,bpcmamdq,bpcmapdq,
& bmcpamdq,bmcpapdq,bpcpamdq,bpcpapdq,
& work,mwork,rpn3,rpt3)
c ==================================================================
c
c Author: Randall J. LeVeque
c Modified for AMROC: Ralf Deiterding
c
c # Compute the modification to fluxes f, g and h that are generated by
c # all interfaces along a 1D slice of the 3D grid.
c # ixyz = 1 if it is a slice in x
c # 2 if it is a slice in y
c # 3 if it is a slice in z
c # This value is passed into the Riemann solvers. The flux modifications
c # go into the arrays fadd, gadd and hadd. The notation is written
c # assuming we are solving along a 1D slice in the x-direction.
c
c # fadd(i,.) modifies F to the left of cell i
c # gadd(i,.,1,slice) modifies G below cell i (in the z-direction)
c # gadd(i,.,2,slice) modifies G above cell i
c # The G flux in the surrounding slices may
c # also be updated.
c # slice = -1 The slice below in y-direction
c # slice = 0 The slice used in the 2D method
c # slice = 1 The slice above in y-direction
c # hadd(i,.,1,slice) modifies H below cell i (in the y-direction)
c # hadd(i,.,2,slice) modifies H above cell i
c # The H flux in the surrounding slices may
c # also be updated.
c # slice = -1 The slice below in z-direction
c # slice = 0 The slice used in the 2D method
c # slice = 1 The slice above in z-direction
c #
c # The method used is specified by method(2) and method(3):
c
c method(2) = 1 if only first order increment waves are to be used.
c = 2 if second order correction terms are to be added, with
c a flux limiter as specified by mthlim.
c = 3 ==> Slope limiting of the conserved variables, with a
c slope limiter as specified by mthlim(1).
c > 3 ==> User defined slope limiter method.
c Slope-limiting is intended to be used with dimensional
c splitting, but schemes that return the fluctuations might
c also be used in combination with transverse wave
c propagation. Note, that the application of MUSCL before
c transverse propagation corresponds to method(3)=2
c although internally the algorithm of method(3)=1 is
c used.
c
c method(3) specify how the transverse wave propagation
c of the increment wave and the correction wave are performed.
c Note that method(3) is given by a two digit number, in
c contrast to what is the case for claw2. It is convenient
c to define the scheme using the pair (method(2),method(3)).
c
c method(3) = -1 Gives dimensional splitting using Godunov
c splitting, i.e. formally first order accurate.
c = -2 Dimensional splitting using Godunov splitting with
c boundary update after each directional step.
c The necessary ghost cell synchronization is done by
c the surrounding AMROC framework.
c This selection ensures that the solution of the
c splitting method is independent of the number of
c computational nodes.
c 0 Gives the Donor cell method. No transverse
c propagation of neither the increment wave
c nor the correction wave.
c = 10 Transverse propagation of the increment wave
c as in 2D. Note that method (2,10) is
c unconditionally unstable.
c = 11 Corner transport upwind of the increment
c wave. Note that method (2,11) also is
c unconditionally unstable.
c = 20 Both the increment wave and the correction
c wave propagate as in the 2D case. Only to
c be used with method(2) = 2.
c = 21 Corner transport upwind of the increment wave,
c and the correction wave propagates as in 2D.
c Only to be used with method(2) = 2.
c = 22 3D propagation of both the increment wave and
c the correction wave. Only to be used with
c method(2) = 2.
c
c Recommended settings: First order schemes:
c (1,10) Stable for CFL < 1/2
c (1,11) Stable for CFL < 1
c Second order schemes:
c (2,20) Stable for CFL < 1/2
c (2,22) Stable for CFL < 1
c (3,10) Stable for CFL < 1/2
c (3,11) Stable for CFL < 1
c
c WARNING! The schemes (2,10), (2,11) are unconditionally
c unstable.
c
c ----------------------------------
c
c Note that if method(6)=1 then the capa array comes into the second
c order correction terms, and is already included in dtdx1d:
c If ixyz = 1 then
c dtdx1d(i) = dt/dx if method(6) = 0
c = dt/(dx*capa(i,jcom,kcom)) if method(6) = 1
c If ixyz = 2 then
c dtdx1d(j) = dt/dy if method(6) = 0
c = dt/(dy*capa(icom,j,kcom)) if method(6) = 1
c If ixyz = 3 then
c dtdx1d(k) = dt/dz if method(6) = 0
c = dt/(dz*capa(icom,jcom,k)) if method(6) = 1
c
c Notation:
c The jump in q (q1d(i,:)-q1d(i-1,:)) is split by rpn3 into
c amdq = the left-going flux difference A^- Delta q
c apdq = the right-going flux difference A^+ Delta q
c Each of these is split by rpt3 into
c bmasdq = the down-going transverse flux difference B^- A^* Delta q
c bpasdq = the up-going transverse flux difference B^+ A^* Delta q
c where A^* represents either A^- or A^+.
c
c
implicit double precision (a-h,o-z)
include "call.i"
dimension q1d(1-mbc:maxm+mbc, meqn)
dimension amdq(1-mbc:maxm+mbc, meqn)
dimension apdq(1-mbc:maxm+mbc, meqn)
dimension bmamdq(1-mbc:maxm+mbc, meqn)
dimension bmapdq(1-mbc:maxm+mbc, meqn)
dimension bpamdq(1-mbc:maxm+mbc, meqn)
dimension bpapdq(1-mbc:maxm+mbc, meqn)
dimension cqxx(1-mbc:maxm+mbc, meqn)
dimension faddm(1-mbc:maxm+mbc, meqn)
dimension faddp(1-mbc:maxm+mbc, meqn)
dimension gaddm(1-mbc:maxm+mbc, meqn, 2, -1:1)
dimension gaddp(1-mbc:maxm+mbc, meqn, 2, -1:1)
dimension haddm(1-mbc:maxm+mbc, meqn, 2, -1:1)
dimension haddp(1-mbc:maxm+mbc, meqn, 2, -1:1)
c
dimension cmamdq(1-mbc:maxm+mbc, meqn)
dimension cmapdq(1-mbc:maxm+mbc, meqn)
dimension cpamdq(1-mbc:maxm+mbc, meqn)
dimension cpapdq(1-mbc:maxm+mbc, meqn)
c
dimension cmamdq2(1-mbc:maxm+mbc, meqn)
dimension cmapdq2(1-mbc:maxm+mbc, meqn)
dimension cpamdq2(1-mbc:maxm+mbc, meqn)
dimension cpapdq2(1-mbc:maxm+mbc, meqn)
c
dimension bmcqxx(1-mbc:maxm+mbc, meqn)
dimension bpcqxx(1-mbc:maxm+mbc, meqn)
dimension cmcqxx(1-mbc:maxm+mbc, meqn)
dimension cpcqxx(1-mbc:maxm+mbc, meqn)
c
dimension bpcmamdq(1-mbc:maxm+mbc, meqn)
dimension bpcmapdq(1-mbc:maxm+mbc, meqn)
dimension bpcpamdq(1-mbc:maxm+mbc, meqn)
dimension bpcpapdq(1-mbc:maxm+mbc, meqn)
dimension bmcmamdq(1-mbc:maxm+mbc, meqn)
dimension bmcmapdq(1-mbc:maxm+mbc, meqn)
dimension bmcpamdq(1-mbc:maxm+mbc, meqn)
dimension bmcpapdq(1-mbc:maxm+mbc, meqn)
c
dimension dtdx1d(1-mbc:maxm+mbc)
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
dimension s(1-mbc:maxm+mbc, mwaves)
dimension wave(1-mbc:maxm+mbc, meqn, mwaves)
dimension method(7),mthlim(mwaves)
dimension work(mwork)
external rpn3, rpt3
c
logical limit
common /rpnflx/ mrpnflx
c
i0ql = 1
i0qr = i0ql + (maxm+2*mbc)*meqn
i0fl = i0qr + (maxm+2*mbc)*meqn
i0fr = i0fl + (maxm+2*mbc)*meqn
iused = i0fr + (maxm+2*mbc)*meqn - 1
c
if (iused.gt.mwork) then
write(6,*) '*** not enough work space in flux3ex'
write(6,*) '*** iused = ', iused, ' mwork =',mwork
stop
endif
c
limit = .false.
do 5 mw=1,mwaves
if (mthlim(mw) .gt. 0) limit = .true.
5 continue
c
c # initialize flux increments:
c -----------------------------
c
do 11 m=1,meqn
do 11 i = 1-mbc, mx+mbc
faddm(i,m) = 0.d0
faddp(i,m) = 0.d0
11 continue
c
c # local method parameters
if (method(3).gt.0) then
m3 = method(3)/10
m4 = method(3) - 10*m3
c
do 10 jside=1,2
do 10 m=1,meqn
do 10 i = 1-mbc, mx+mbc
gaddm(i,m,jside,-1) = 0.d0
gaddm(i,m,jside, 0) = 0.d0
gaddm(i,m,jside, 1) = 0.d0
gaddp(i,m,jside,-1) = 0.d0
gaddp(i,m,jside, 0) = 0.d0
gaddp(i,m,jside, 1) = 0.d0
haddm(i,m,jside,-1) = 0.d0
haddm(i,m,jside, 0) = 0.d0
haddm(i,m,jside, 1) = 0.d0
haddp(i,m,jside,-1) = 0.d0
haddp(i,m,jside, 0) = 0.d0
haddp(i,m,jside, 1) = 0.d0
10 continue
c
else
m3 = 0
m4 = 0
endif
c
c # solve Riemann problem at each interface and compute Godunov flux F0
c ---------------------------------------------------------------------
c
if (method(2).le.2) then
call rpn3(ixyz,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,
& aux2,aux2,wave,s,amdq,apdq)
else
call slope3(ixyz,maxm,meqn,maux,mwaves,mbc,mx,q1d,aux2,dx,dt,
& method,mthlim,wave,s,amdq,apdq,dtdx1d,rpn3,
& work(i0ql),work(i0qr),work(i0fl),work(i0fr))
endif
c
c # Set fadd for the donor-cell upwind method (Godunov)
do 40 m=1,meqn
do 40 i=1,mx+1
faddp(i,m) = faddp(i,m) - apdq(i,m)
faddm(i,m) = faddm(i,m) + amdq(i,m)
40 continue
c
c # compute maximum wave speed for checking Courant number:
cfl1d = 0.d0
do 50 mw=1,mwaves
do 50 i=1,mx+1
cfl1d = dmax1(cfl1d,dtdx1d(i)*dabs(s(i,mw)))
50 continue
c
if (method(2).eq.1 .and. m3.eq.0 .and. m4.eq.0) go to 999
if (method(2).eq.1.or.method(2).ge.3) go to 130
c
if (mrpnflx.ne.0) then
write(6,*) '*** Riemann solver returns fluxes.'
write(6,*) '*** Wave limiting not possible.'
write(6,*) '*** Set method(2)>=3 for slope limiting.'
stop
endif
c
c # modify F fluxes for second order q_{xx} correction terms:
c -----------------------------------------------------------
c
c # apply limiter to waves:
if (limit) call limiter(maxm,meqn,mwaves,mbc,mx,wave,s,mthlim)
c
do 115 m = 1,meqn
do 115 i = 1-mbc, mx+mbc
cqxx(i,m) = 0.d0
115 continue
c
do 120 i = 1, mx+1
c
c # For correction terms below, need average of dtdx in cell
c # i-1 and i. Compute these and overwrite dtdx1d:
c
dtdx1d(i-1) = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i))
c
do 120 m=1,meqn
do 119 mw=1,mwaves
cqxx(i,m) = cqxx(i,m) + 0.5d0 * dabs(s(i,mw))
& * (1.d0 - dabs(s(i,mw))*dtdx1d(i-1)) * wave(i,m,mw)
119 continue
faddm(i,m) = faddm(i,m) + cqxx(i,m)
faddp(i,m) = faddp(i,m) + cqxx(i,m)
120 continue
c
130 continue
c
if( m3.eq.0 .and. m4.eq.0 ) go to 999
c
if (mrpnflx.ne.0) then
write(6,*) '*** Riemann solver returns fluxes.'
write(6,*) '*** Transverse wave propagation not possible.'
write(6,*) '*** Set method(3)<0 for dimensional splitting.'
stop
endif
c
c # modify G fluxes for transverse propagation
c --------------------------------------------
c
c # split the left-going flux difference into down-going and up-going
c # flux differences (in the y-direction).
c
call rpt3(ixyz,2,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,aux1,aux2,
& aux3,1,amdq,bmamdq,bpamdq)
c
c # split the right-going flux difference into down-going and up-going
c # flux differences (in the y-direction).
c
call rpt3(ixyz,2,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,aux1,aux2,
& aux3,2,apdq,bmapdq,bpapdq)
c
c # split the left-going flux difference into down-going and up-going
c # flux differences (in the z-direction).
c
call rpt3(ixyz,3,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,aux1,aux2,
& aux3,1,amdq,cmamdq,cpamdq)
c
c # split the right-going flux difference into down-going and up-going
c # flux differences (in the y-direction).
c
call rpt3(ixyz,3,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,aux1,aux2,
& aux3,2,apdq,cmapdq,cpapdq)
c
c # Split the correction wave into transverse propagating waves
c # in the y-direction and z-direction.
c
if (m3.eq.2.and.method(2).le.2) then
call rpt3(ixyz,2,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,
& aux1,aux2,aux3,0,cqxx,bmcqxx,bpcqxx)
call rpt3(ixyz,3,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,
& aux1,aux2,aux3,0,cqxx,cmcqxx,cpcqxx)
endif
c
c # If the correction wave also propagates in a 3D sense, incorporate
c # cpcqxx,... into cmamdq, cpamdq, ... so that it is split also.
c
if(m4 .eq. 1)then
do 145 m = 1, meqn
do 145 i = 0, mx+2
cpapdq2(i,m) = cpapdq(i,m)
cpamdq2(i,m) = cpamdq(i,m)
cmapdq2(i,m) = cmapdq(i,m)
cmamdq2(i,m) = cmamdq(i,m)
145 continue
else if(m4.eq.2.and.method(2).le.2)then
do 146 m = 1, meqn
do 146 i = 0, mx+2
cpapdq2(i,m) = cpapdq(i,m) - 3.d0*cpcqxx(i,m)
cpamdq2(i,m) = cpamdq(i,m) + 3.d0*cpcqxx(i,m)
cmapdq2(i,m) = cmapdq(i,m) - 3.d0*cmcqxx(i,m)
cmamdq2(i,m) = cmamdq(i,m) + 3.d0*cmcqxx(i,m)
146 continue
endif
c
c # The transverse flux differences in the z-direction are split
c # into waves propagating in the y-direction. If m4 = 2,
c # then the transverse propagating correction waves in the z-direction
c # are also split. This yields terms of the form BCAu_{xzy} and
c # BCAAu_{xxzy}.
c
if( m4.gt.0 )then
call rpt3(ixyz,2,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,
& aux1,aux2,aux3,0,cpapdq2,bmcpapdq,bpcpapdq)
call rpt3(ixyz,2,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,
& aux1,aux2,aux3,0,cpamdq2,bmcpamdq,bpcpamdq)
call rpt3(ixyz,2,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,
& aux1,aux2,aux3,0,cmapdq2,bmcmapdq,bpcmapdq)
call rpt3(ixyz,2,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,
& aux1,aux2,aux3,0,cmamdq2,bmcmamdq,bpcmamdq)
endif
c
c # The updates--------------------------------------------------
c
do 180 m=1,meqn
do 180 i = 1, mx+1
c
c # Transverse propagation of the increment waves
c # between cells sharing interfaces, i.e. the 2D approach.
c # Yields BAu_{xy}.
c
gupdate = 0.5d0*dtdx1d(i-1)*bmamdq(i,m)
gaddm(i-1,m,1,0) = gaddm(i-1,m,1,0) - gupdate
gaddp(i-1,m,1,0) = gaddp(i-1,m,1,0) - gupdate
gupdate = 0.5d0*dtdx1d(i-1)*bpamdq(i,m)
gaddm(i-1,m,2,0) = gaddm(i-1,m,2,0) - gupdate
gaddp(i-1,m,2,0) = gaddp(i-1,m,2,0) - gupdate
gupdate = 0.5d0*dtdx1d(i-1)*bmapdq(i,m)
gaddm(i,m,1,0) = gaddm(i,m,1,0) - gupdate
gaddp(i,m,1,0) = gaddp(i,m,1,0) - gupdate
gupdate = 0.5d0*dtdx1d(i-1)*bpapdq(i,m)
gaddm(i,m,2,0) = gaddm(i,m,2,0) - gupdate
gaddp(i,m,2,0) = gaddp(i,m,2,0) - gupdate
c
c # Transverse propagation of the increment wave (and the
c # correction wave if m4=2) between cells
c # only having a corner or edge in common. Yields terms of the
c # BCAu_{xzy} and BCAAu_{xxzy}.
c
if( m4.gt.0 )then
c
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz
& * (bpcpapdq(i,m) - bpcmapdq(i,m))
gaddm(i,m,2,0) = gaddm(i,m,2,0) + gupdate
gaddp(i,m,2,0) = gaddp(i,m,2,0) + gupdate
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz
& * (bmcpapdq(i,m) - bmcmapdq(i,m))
gaddm(i,m,1,0) = gaddm(i,m,1,0) + gupdate
gaddp(i,m,1,0) = gaddp(i,m,1,0) + gupdate
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz*bpcpapdq(i,m)
gaddm(i,m,2,1) = gaddm(i,m,2,1) - gupdate
gaddp(i,m,2,1) = gaddp(i,m,2,1) - gupdate
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz*bmcpapdq(i,m)
gaddm(i,m,1,1) = gaddm(i,m,1,1) - gupdate
gaddp(i,m,1,1) = gaddp(i,m,1,1) - gupdate
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz*bpcmapdq(i,m)
gaddm(i,m,2,-1) = gaddm(i,m,2,-1) + gupdate
gaddp(i,m,2,-1) = gaddp(i,m,2,-1) + gupdate
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz*bmcmapdq(i,m)
gaddm(i,m,1,-1) = gaddm(i,m,1,-1) + gupdate
gaddp(i,m,1,-1) = gaddp(i,m,1,-1) + gupdate
c
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz
& * (bpcpamdq(i,m) - bpcmamdq(i,m))
gaddm(i-1,m,2,0) = gaddm(i-1,m,2,0) + gupdate
gaddp(i-1,m,2,0) = gaddp(i-1,m,2,0) + gupdate
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz
& * (bmcpamdq(i,m) - bmcmamdq(i,m))
gaddm(i-1,m,1,0) = gaddm(i-1,m,1,0) + gupdate
gaddp(i-1,m,1,0) = gaddp(i-1,m,1,0) + gupdate
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz*bpcpamdq(i,m)
gaddm(i-1,m,2,1) = gaddm(i-1,m,2,1) - gupdate
gaddp(i-1,m,2,1) = gaddp(i-1,m,2,1) - gupdate
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz*bmcpamdq(i,m)
gaddm(i-1,m,1,1) = gaddm(i-1,m,1,1) - gupdate
gaddp(i-1,m,1,1) = gaddp(i-1,m,1,1) - gupdate
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz*bpcmamdq(i,m)
gaddm(i-1,m,2,-1) = gaddm(i-1,m,2,-1) + gupdate
gaddp(i-1,m,2,-1) = gaddp(i-1,m,2,-1) + gupdate
gupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdz*bmcmamdq(i,m)
gaddm(i-1,m,1,-1) = gaddm(i-1,m,1,-1) + gupdate
gaddp(i-1,m,1,-1) = gaddp(i-1,m,1,-1) + gupdate
c
endif
c
c # Transverse propagation of the correction wave between
c # cells sharing faces. This gives BAAu_{xxy}.
c
if(m3.lt.2.or.method(2).ge.3) go to 180
gupdate = dtdx1d(i-1)*bpcqxx(i,m)
gaddm(i,m,2,0) = gaddm(i,m,2,0) + gupdate
gaddp(i,m,2,0) = gaddp(i,m,2,0) + gupdate
gupdate = dtdx1d(i-1)*bmcqxx(i,m)
gaddm(i,m,1,0) = gaddm(i,m,1,0) + gupdate
gaddp(i,m,1,0) = gaddp(i,m,1,0) + gupdate
gupdate = dtdx1d(i-1)*bpcqxx(i,m)
gaddm(i-1,m,2,0) = gaddm(i-1,m,2,0) - gupdate
gaddp(i-1,m,2,0) = gaddp(i-1,m,2,0) - gupdate
gupdate = dtdx1d(i-1)*bmcqxx(i,m)
gaddm(i-1,m,1,0) = gaddm(i-1,m,1,0) - gupdate
gaddp(i-1,m,1,0) = gaddp(i-1,m,1,0) - gupdate
c
180 continue
c
c
c # modify H fluxes
c --------------------------------------------
c
c # If the correction wave also propagates in a 3D sense, incorporate
c # cqxx into bmamdq, bpamdq, ... so that is is split also.
c
if(m4.eq.2.and.method(2).le.2)then
do 155 m = 1, meqn
do 155 i = 0, mx+2
bpapdq(i,m) = bpapdq(i,m) - 3.d0*bpcqxx(i,m)
bpamdq(i,m) = bpamdq(i,m) + 3.d0*bpcqxx(i,m)
bmapdq(i,m) = bmapdq(i,m) - 3.d0*bmcqxx(i,m)
bmamdq(i,m) = bmamdq(i,m) + 3.d0*bmcqxx(i,m)
155 continue
endif
c
c # The transverse flux differences in the y-direction are split
c # into waves propagating in the z-direction. If m4 = 2,
c # then the transverse propagating correction waves in the y-direction
c # are also split. This yields terms of the form BCAu_{xzy} and
c # BCAAu_{xxzy}.
c
if( m4.gt.0 )then
call rpt3(ixyz,3,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,
& aux1,aux2,aux3,0,bpapdq,bmcpapdq,bpcpapdq)
call rpt3(ixyz,3,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,
& aux1,aux2,aux3,0,bpamdq,bmcpamdq,bpcpamdq)
call rpt3(ixyz,3,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,
& aux1,aux2,aux3,0,bmapdq,bmcmapdq,bpcmapdq)
call rpt3(ixyz,3,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,
& aux1,aux2,aux3,0,bmamdq,bmcmamdq,bpcmamdq)
endif
c
c # The updates--------------------------------------------------
c
do 200 m=1,meqn
do 200 i = 1, mx+1
c
c # Transverse propagation of the increment waves
c # between cells sharing interfaces, i.e. the 2D approach.
c # Yields CAu_{xy}.
c
hupdate = 0.5d0*dtdx1d(i-1)*cmamdq(i,m)
haddm(i-1,m,1,0) = haddm(i-1,m,1,0) - hupdate
haddp(i-1,m,1,0) = haddp(i-1,m,1,0) - hupdate
hupdate = 0.5d0*dtdx1d(i-1)*cpamdq(i,m)
haddm(i-1,m,2,0) = haddm(i-1,m,2,0) - hupdate
haddp(i-1,m,2,0) = haddp(i-1,m,2,0) - hupdate
hupdate = 0.5d0*dtdx1d(i-1)*cmapdq(i,m)
haddm(i,m,1,0) = haddm(i,m,1,0) - hupdate
haddp(i,m,1,0) = haddp(i,m,1,0) - hupdate
hupdate = 0.5d0*dtdx1d(i-1)*cpapdq(i,m)
haddm(i,m,2,0) = haddm(i,m,2,0) - hupdate
haddp(i,m,2,0) = haddp(i,m,2,0) - hupdate
c
c # Transverse propagation of the increment wave (and the
c # correction wave if m4=2) between cells
c # only having a corner or edge in common. Yields terms of the
c # CBAu_{xzy} and CBAAu_{xxzy}.
c
if( m4.gt.0 )then
c
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy
& * (bpcpapdq(i,m) - bpcmapdq(i,m))
haddm(i,m,2,0) = haddm(i,m,2,0) + hupdate
haddp(i,m,2,0) = haddp(i,m,2,0) + hupdate
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy
& * (bmcpapdq(i,m) - bmcmapdq(i,m))
haddm(i,m,1,0) = haddm(i,m,1,0) + hupdate
haddp(i,m,1,0) = haddp(i,m,1,0) + hupdate
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy*bpcpapdq(i,m)
haddm(i,m,2,1) = haddm(i,m,2,1) - hupdate
haddp(i,m,2,1) = haddp(i,m,2,1) - hupdate
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy*bmcpapdq(i,m)
haddm(i,m,1,1) = haddm(i,m,1,1) - hupdate
haddp(i,m,1,1) = haddp(i,m,1,1) - hupdate
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy*bpcmapdq(i,m)
haddm(i,m,2,-1) = haddm(i,m,2,-1) + hupdate
haddp(i,m,2,-1) = haddp(i,m,2,-1) + hupdate
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy*bmcmapdq(i,m)
haddm(i,m,1,-1) = haddm(i,m,1,-1) + hupdate
haddp(i,m,1,-1) = haddp(i,m,1,-1) + hupdate
c
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy
& * (bpcpamdq(i,m) - bpcmamdq(i,m))
haddm(i-1,m,2,0) = haddm(i-1,m,2,0) + hupdate
haddp(i-1,m,2,0) = haddp(i-1,m,2,0) + hupdate
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy
& * (bmcpamdq(i,m) - bmcmamdq(i,m))
haddm(i-1,m,1,0) = haddm(i-1,m,1,0) + hupdate
haddp(i-1,m,1,0) = haddp(i-1,m,1,0) + hupdate
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy*bpcpamdq(i,m)
haddm(i-1,m,2,1) = haddm(i-1,m,2,1) - hupdate
haddp(i-1,m,2,1) = haddp(i-1,m,2,1) - hupdate
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy*bmcpamdq(i,m)
haddm(i-1,m,1,1) = haddm(i-1,m,1,1) - hupdate
haddp(i-1,m,1,1) = haddp(i-1,m,1,1) - hupdate
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy*bpcmamdq(i,m)
haddm(i-1,m,2,-1) = haddm(i-1,m,2,-1) + hupdate
haddp(i-1,m,2,-1) = haddp(i-1,m,2,-1) + hupdate
hupdate = (1.d0/6.d0)*dtdx1d(i-1)*dtdy*bmcmamdq(i,m)
haddm(i-1,m,1,-1) = haddm(i-1,m,1,-1) + hupdate
haddp(i-1,m,1,-1) = haddp(i-1,m,1,-1) + hupdate
c
endif
c
c # Transverse propagation of the correction wave between
c # cells sharing faces. This gives CAAu_{xxy}.
c
if(m3.lt.2.or.method(2).ge.3) go to 200
hupdate = dtdx1d(i-1)*cpcqxx(i,m)
haddm(i,m,2,0) = haddm(i,m,2,0) + hupdate
haddp(i,m,2,0) = haddp(i,m,2,0) + hupdate
hupdate = dtdx1d(i-1)*cmcqxx(i,m)
haddm(i,m,1,0) = haddm(i,m,1,0) + hupdate
haddp(i,m,1,0) = haddp(i,m,1,0) + hupdate
hupdate = dtdx1d(i-1)*cpcqxx(i,m)
haddm(i-1,m,2,0) = haddm(i-1,m,2,0) - hupdate
haddp(i-1,m,2,0) = haddp(i-1,m,2,0) - hupdate
hupdate = dtdx1d(i-1)*cmcqxx(i,m)
haddm(i-1,m,1,0) = haddm(i-1,m,1,0) - hupdate
haddp(i-1,m,1,0) = haddp(i-1,m,1,0) - hupdate
c
200 continue
c
999 continue
return
end
c
c
c ===================================================================
subroutine slope3(ixyz,maxm,meqn,maux,mwaves,mbc,mx,
& q,aux,dx,dt,method,mthlim,
& wave,s,amdq,apdq,dtdx,rpn3,
& ql,qr,fl,fr)
c ===================================================================
c
c # Implements the standard MUSCL-Hancock method. MUSCL must
c # be used to obtain 2nd order accuracy with conventional
c # finite-volume schemes that return the numerical fluxes instead
c # of fluctuations. Schemes returning the fluctuations can use MUSCL
c # slope-limiting or wave-limiting. Slope-limiting is intended to
c # be used with dimensional splitting, but wave propagation schemes
c # also can apply it in combination with transverse wave propagation.
c
c # Author: Ralf Deiterding, ralf@cacr.caltech.edu
c
implicit double precision (a-h,o-z)
include "call.i"
common /rpnflx/ mrpnflx
c
dimension q(1-mbc:maxm+mbc, meqn)
dimension ql(1-mbc:maxm+mbc, meqn)
dimension qr(1-mbc:maxm+mbc, meqn)
dimension aux(1-mbc:maxm+mbc, maux, 3)
dimension fl(1-mbc:maxm+mbc, meqn)
dimension fr(1-mbc:maxm+mbc, meqn)
dimension wave(1-mbc:maxm+mbc, meqn, mwaves)
dimension s(1-mbc:maxm+mbc, mwaves)
dimension amdq(1-mbc:maxm+mbc, meqn)
dimension apdq(1-mbc:maxm+mbc, meqn)
dimension dtdx(1-mbc:maxm+mbc)
dimension method(7),mthlim(mwaves)
external rpn3
c
mlim = 0
do 90 mw=1,mwaves
if (mthlim(mw) .gt. 0) then
mlim = mthlim(mw)
goto 95
endif
90 continue
95 continue
c
do 100 m=1,meqn
ql(1-mbc,m) = q(1-mbc,m)
qr(1-mbc,m) = q(1-mbc,m)
ql(mx+mbc,m) = q(mx+mbc,m)
qr(mx+mbc,m) = q(mx+mbc,m)
100 continue
c
if (method(2).gt.3) then
call rec3(ixyz,maxm,meqn,mwaves,mbc,mx,q,method,mthlim,ql,qr)
c
c # MUSCL reconstruction with slope-limiting for conserved
c # variables. Linear reconstruction: om=0.d0
c # Quadratic spatial reconstuction: om!=0.d0
c # 2nd order spatial reconstruction: om=1.d0/3.d0
c # Reconstructions with om!=0.d0 are not conservative!
c
else
om = 0.d0
do 110 i=2-mbc,mx+mbc-1
do 110 m=1,meqn
call reclim(q(i,m),q(i-1,m),q(i+1,m),
& mlim,om,ql(i,m),qr(i,m))
110 continue
endif
c
call flx3(ixyz,maxm,meqn,mbc,mx,ql,maux,aux,fl)
call flx3(ixyz,maxm,meqn,mbc,mx,qr,maux,aux,fr)
c
do 200 i=2-mbc,mx+mbc-1
do 200 m=1,meqn
ql(i,m) = ql(i,m) + 0.5d0*dtdx(i)*(fl(i,m)-fr(i,m))
qr(i,m) = qr(i,m) + 0.5d0*dtdx(i)*(fl(i,m)-fr(i,m))
200 continue
c
call rpn3(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,maux,
& aux,aux,wave,s,amdq,apdq)
c
c # Add differences of fluxes between original and reconstructed values
c # to fluctuations, if a wave propagation scheme is applied.
c
if (mrpnflx.eq.0) then
call flx3(ixyz,maxm,meqn,mbc,mx,ql,maux,aux,fl)
call flx3(ixyz,maxm,meqn,mbc,mx,qr,maux,aux,fr)
do 300 i=2-mbc,mx+mbc-1
do 300 m=1,meqn
amdq(i,m) = amdq(i,m) + fr(i-1,m)
apdq(i,m) = apdq(i,m) - fl(i ,m)
300 continue
call flx3(ixyz,maxm,meqn,mbc,mx,q,maux,aux,fl)
do 310 i=2-mbc,mx+mbc-1
do 310 m=1,meqn
amdq(i,m) = amdq(i,m) - fl(i-1,m)
apdq(i,m) = apdq(i,m) + fl(i ,m)
310 continue
endif
c
return
end
c
c
c ===================================================================
subroutine saverec3(ixyz,maxm,maxmx,maxmy,maxmz,mvar,meqn,mbc,
& mx,my,mz,q,qls,qrs,qbs,qts,qfs,qks,ql,qr)
c ===================================================================
c
c # Store reconstructed values for later use.
c
implicit double precision (a-h,o-z)
include "call.i"
c
dimension q(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension qls(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension qrs(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension qbs(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension qts(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension qfs(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension qks(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension ql(1-mbc:maxm+mbc, meqn)
dimension qr(1-mbc:maxm+mbc, meqn)
c
if (ixyz.eq.1) then
do 10 i = 1-mbc, mx+mbc
do m=1,meqn
qls(m,i,jcom,kcom) = ql(i,m)
qrs(m,i,jcom,kcom) = qr(i,m)
enddo
do m=meqn+1,mvar
qls(m,i,jcom,kcom) = q(m,i,jcom,kcom)
qrs(m,i,jcom,kcom) = q(m,i,jcom,kcom)
enddo
10 continue
endif
c
if (ixyz.eq.2) then
do 20 j = 1-mbc, my+mbc
do m=1,meqn
qbs(m,icom,j,kcom) = ql(j,m)
qts(m,icom,j,kcom) = qr(j,m)
enddo
do m=meqn+1,mvar
qbs(m,icom,j,kcom) = q(m,icom,j,kcom)
qts(m,icom,j,kcom) = q(m,icom,j,kcom)
enddo
20 continue
endif
c
if (ixyz.eq.3) then
do 30 k = 1-mbc, mz+mbc
do m=1,meqn
qfs(m,icom,jcom,k) = ql(k,m)
qks(m,icom,jcom,k) = qr(k,m)
enddo
do m=meqn+1,mvar
qfs(m,icom,jcom,k) = q(m,icom,jcom,k)
qks(m,icom,jcom,k) = q(m,icom,jcom,k)
enddo
30 continue
endif
c
return
end
c