c
c ==================================================================
subroutine unsp3(maxm,maxmx,maxmy,maxmz,mvar,meqn,
& maux,mwaves,mbc,mx,my,mz,
& qold,aux,dx,dy,dz,dt,method,mthlim,cfl,
& fm,fp,gm,gp,hm,hp,
& faddm,faddp,gaddm,gaddp,haddm,haddp,
& q1d,dtdx1d,dtdy1d,dtdz1d,
& aux1,aux2,aux3,work,mwork,rpn3,rpt3)
c ==================================================================
c
c # Take one time step, updating q with the wave propagation method
c # of Randall J. LeVeque.
c
c # fm, fp are fluxes to left and right of single cell edge,
c # gm, gp are fluxes at bottom and top,
c # hm, hp are fluxes at back and front.
c # fadd, gadd and hadd are used to return flux increments from flux3.
c # See the flux3 documentation for more information.
c
implicit double precision (a-h,o-z)
include "call.i"
c
dimension qold(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension fm(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension fp(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension gm(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension gp(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension hm(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension hp(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension q1d(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)
dimension aux(maux, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+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)
dimension dtdx1d(1-mbc:maxm+mbc)
dimension dtdy1d(1-mbc:maxm+mbc)
dimension dtdz1d(1-mbc:maxm+mbc)
dimension work(mwork)
dimension method(7),mthlim(mwaves)
external rpn3, rpt3
c
c # partition work array into pieces needed for local storage in
c # flux3 routine. Find starting index of each piece:
c
i0wave = 1
i0s = i0wave + (maxm+2*mbc)*meqn*mwaves
i0amdq = i0s + (maxm+2*mbc)*mwaves
i0apdq = i0amdq + (maxm+2*mbc)*meqn
i0cqxx = i0apdq + (maxm+2*mbc)*meqn
i0bmamdq = i0cqxx + (maxm+2*mbc)*meqn
i0bmapdq = i0bmamdq + (maxm+2*mbc)*meqn
i0bpamdq = i0bmapdq + (maxm+2*mbc)*meqn
i0bpapdq = i0bpamdq + (maxm+2*mbc)*meqn
i0cmamdq = i0bpapdq + (maxm+2*mbc)*meqn
i0cmapdq = i0cmamdq + (maxm+2*mbc)*meqn
i0cpamdq = i0cmapdq + (maxm+2*mbc)*meqn
i0cpapdq = i0cpamdq + (maxm+2*mbc)*meqn
i0cmamdq2 = i0cpapdq + (maxm+2*mbc)*meqn
i0cmapdq2 = i0cmamdq2 + (maxm+2*mbc)*meqn
i0cpamdq2 = i0cmapdq2 + (maxm+2*mbc)*meqn
i0cpapdq2 = i0cpamdq2 + (maxm+2*mbc)*meqn
i0bmcqxx = i0cpapdq2 + (maxm+2*mbc)*meqn
i0bpcqxx = i0bmcqxx + (maxm+2*mbc)*meqn
i0cmcqxx = i0bpcqxx + (maxm+2*mbc)*meqn
i0cpcqxx = i0cmcqxx + (maxm+2*mbc)*meqn
i0bmcmamdq = i0cpcqxx + (maxm+2*mbc)*meqn
i0bmcmapdq = i0bmcmamdq + (maxm+2*mbc)*meqn
i0bpcmamdq = i0bmcmapdq + (maxm+2*mbc)*meqn
i0bpcmapdq = i0bpcmamdq + (maxm+2*mbc)*meqn
i0bmcpamdq = i0bpcmapdq + (maxm+2*mbc)*meqn
i0bmcpapdq = i0bmcpamdq + (maxm+2*mbc)*meqn
i0bpcpamdq = i0bmcpapdq + (maxm+2*mbc)*meqn
i0bpcpapdq = i0bpcpamdq + (maxm+2*mbc)*meqn
if (method(1).ge.1) then
i0qls = i0bpcpapdq + (maxm+2*mbc)*meqn
i0qrs = i0qls + (maxmx+2*mbc)*(maxmy+2*mbc)*(maxmz+2*mbc)*mvar
i0qbs = i0qrs + (maxmx+2*mbc)*(maxmy+2*mbc)*(maxmz+2*mbc)*mvar
i0qts = i0qbs + (maxmx+2*mbc)*(maxmy+2*mbc)*(maxmz+2*mbc)*mvar
i0qfs = i0qts + (maxmx+2*mbc)*(maxmy+2*mbc)*(maxmz+2*mbc)*mvar
i0qks = i0qfs + (maxmx+2*mbc)*(maxmy+2*mbc)*(maxmz+2*mbc)*mvar
i0ql = i0qks + (maxmx+2*mbc)*(maxmy+2*mbc)*(maxmz+2*mbc)*mvar
i0qr = i0ql + (maxm+2*mbc)*meqn
i0slope = i0ql
else
i0slope = i0bpcpapdq + (maxm+2*mbc)*meqn
endif
iused = i0slope + (maxm+2*mbc)*meqn*4 - 1
mslope = iused-i0slope+1
c
if (iused.gt.mwork) then
write(6,*) '*** not enough work space in unsp3'
write(6,*) '*** iused = ', iused, ' mwork =',mwork
stop
endif
c
mcapa = method(6)
c
cfl = 0.d0
dtdx = dt/dx
dtdy = dt/dy
dtdz = dt/dz
c
do 10 k=1-mbc,mz+mbc
do 10 j=1-mbc,my+mbc
do 10 i=1-mbc,mx+mbc
do 10 m=1,mvar
fm(m,i,j,k) = 0.d0
fp(m,i,j,k) = 0.d0
gm(m,i,j,k) = 0.d0
gp(m,i,j,k) = 0.d0
hm(m,i,j,k) = 0.d0
hp(m,i,j,k) = 0.d0
10 continue
c
if (mcapa.eq.0) then
c # no capa array:
do 5 i=1-mbc,maxm+mbc
dtdx1d(i) = dtdx
dtdy1d(i) = dtdy
dtdz1d(i) = dtdz
5 continue
endif
c
c # method parameters
if (method(3).gt.0) then
m3 = method(3)/10
m4 = method(3) - 10*m3
c
else
m3 = 0
m4 = 0
endif
c
c
c # perform x-sweeps
c ==================
c
do 50 k = 0,mz+1
do 50 j = 0,my+1
c
do 20 i = 1-mbc, mx+mbc
do 20 m=1,meqn
c # copy data along a slice into 1d array:
q1d(i,m) = qold(m,i,j,k)
20 continue
c
if (mcapa.gt.0) then
do 23 i = 1-mbc, mx+mbc
dtdx1d(i) = dtdx / aux(mcapa,i,j,k)
23 continue
endif
c
if (maux .gt. 0) then
do 22 i = 1-mbc, mx+mbc
do 22 ma=1,maux
aux1(i,ma,1) = aux(ma,i,j-1,k-1)
aux1(i,ma,2) = aux(ma,i,j-1,k)
aux1(i,ma,3) = aux(ma,i,j-1,k+1)
aux2(i,ma,1) = aux(ma,i,j,k-1)
aux2(i,ma,2) = aux(ma,i,j,k)
aux2(i,ma,3) = aux(ma,i,j,k+1)
aux3(i,ma,1) = aux(ma,i,j+1,k-1)
aux3(i,ma,2) = aux(ma,i,j+1,k)
aux3(i,ma,3) = aux(ma,i,j+1,k+1)
22 continue
endif
c
c # Store the value of j and k along this slice in the common block
c # comxyt in case it is needed in the Riemann solver (for
c # variable coefficient problems)
c
jcom = j
kcom = k
c
c # compute modifications fadd, gadd and hadd to fluxes along
c # this slice:
c
call flux3(1,maxm,meqn,maux,mwaves,mbc,mx,
& q1d,dtdx1d,dtdy,dtdz,aux1,aux2,aux3,
& method,mthlim,
& faddm,faddp,gaddm,gaddp,haddm,haddp,cfl1d,
& work(i0wave),work(i0s),work(i0amdq),
& work(i0apdq),work(i0cqxx),
& work(i0bmamdq),work(i0bmapdq),
& work(i0bpamdq),work(i0bpapdq),
& work(i0cmamdq),work(i0cmapdq),
& work(i0cpamdq),work(i0cpapdq),
& work(i0cmamdq2),work(i0cmapdq2),
& work(i0cpamdq2),work(i0cpapdq2),
& work(i0bmcqxx),work(i0bpcqxx),
& work(i0cmcqxx),work(i0cpcqxx),
& work(i0bmcmamdq),work(i0bmcmapdq),
& work(i0bpcmamdq),work(i0bpcmapdq),
& work(i0bmcpamdq),work(i0bmcpapdq),
& work(i0bpcpamdq),work(i0bpcpapdq),
& work(i0slope),mslope,rpn3,rpt3)
c
cfl = dmax1(cfl,cfl1d)
c
c # update arrays f, g and h
c
do 30 i=1,mx+1
do 30 m=1,meqn
fm(m,i,j,k) = fm(m,i,j,k) + faddm(i,m)
fp(m,i,j,k) = fp(m,i,j,k) + faddp(i,m)
30 continue
c
if (m3.gt.0) then
do 32 i=1,mx+1
do 32 m=1,meqn
gm(m,i,j,k ) = gm(m,i,j,k ) + gaddm(i,m,1,0)
gp(m,i,j,k ) = gp(m,i,j,k ) + gaddp(i,m,1,0)
gm(m,i,j+1,k ) = gm(m,i,j+1,k ) + gaddm(i,m,2,0)
gp(m,i,j+1,k ) = gp(m,i,j+1,k ) + gaddp(i,m,2,0)
if (m4.gt.0) then
gm(m,i,j,k-1) = gm(m,i,j,k-1) + gaddm(i,m,1,-1)
gp(m,i,j,k-1) = gp(m,i,j,k-1) + gaddp(i,m,1,-1)
gm(m,i,j,k+1) = gm(m,i,j,k+1) + gaddm(i,m,1,1)
gp(m,i,j,k+1) = gp(m,i,j,k+1) + gaddp(i,m,1,1)
gm(m,i,j+1,k-1)=gm(m,i,j+1,k-1)+gaddm(i,m,2,-1)
gp(m,i,j+1,k-1)=gp(m,i,j+1,k-1)+gaddp(i,m,2,-1)
gm(m,i,j+1,k+1)=gm(m,i,j+1,k+1)+gaddm(i,m,2,1)
gp(m,i,j+1,k+1)=gp(m,i,j+1,k+1)+gaddp(i,m,2,1)
endif
32 continue
c
do 34 i=1,mx+1
do 34 m=1,meqn
hm(m,i,j, k) = hm(m,i,j, k) + haddm(i,m,1,0)
hp(m,i,j, k) = hp(m,i,j, k) + haddp(i,m,1,0)
hm(m,i,j, k+1) = hm(m,i,j, k+1) + haddm(i,m,2,0)
hp(m,i,j, k+1) = hp(m,i,j, k+1) + haddp(i,m,2,0)
if (m4.gt.0) then
hm(m,i,j-1,k) = hm(m,i,j-1,k) + haddm(i,m,1,-1)
hp(m,i,j-1,k) = hp(m,i,j-1,k) + haddp(i,m,1,-1)
hm(m,i,j+1,k) = hm(m,i,j+1,k) + haddm(i,m,1,1)
hp(m,i,j+1,k) = hp(m,i,j+1,k) + haddp(i,m,1,1)
hm(m,i,j-1,k+1)=hm(m,i,j-1,k+1)+haddm(i,m,2,-1)
hp(m,i,j-1,k+1)=hp(m,i,j-1,k+1)+haddp(i,m,2,-1)
hm(m,i,j+1,k+1)=hm(m,i,j+1,k+1)+haddm(i,m,2,1)
hp(m,i,j+1,k+1)=hp(m,i,j+1,k+1)+haddp(i,m,2,1)
endif
34 continue
endif
c
if (method(1).ge.1.and.method(2).ge.3)
& call saverec3(1,maxm,maxmx,maxmy,maxmz,mvar,meqn,
& mbc,mx,my,mz,qold,work(i0qls),work(i0qrs),
& work(i0qbs),work(i0qts),work(i0qfs),work(i0qks),
& work(i0ql),work(i0qr))
c
50 continue
c
c
c # perform y sweeps
c ==================
c
c
do 100 k = 0, mz+1
do 100 i = 0, mx+1
c
do 70 j = 1-mbc, my+mbc
do 70 m=1,meqn
c # copy data along a slice into 1d array:
q1d(j,m) = qold(m,i,j,k)
70 continue
c
if (mcapa.gt.0) then
do 71 j = 1-mbc, my+mbc
dtdy1d(j) = dtdy / aux(mcapa,i,j,k)
71 continue
endif
c
if (maux .gt. 0) then
do 72 j = 1-mbc, my+mbc
do 72 ma=1,maux
aux1(j,ma,1) = aux(ma,i-1,j,k-1)
aux1(j,ma,2) = aux(ma,i ,j,k-1)
aux1(j,ma,3) = aux(ma,i+1,j,k-1)
aux2(j,ma,1) = aux(ma,i-1,j,k)
aux2(j,ma,2) = aux(ma,i ,j,k)
aux2(j,ma,3) = aux(ma,i+1,j,k)
aux3(j,ma,1) = aux(ma,i-1,j,k+1)
aux3(j,ma,2) = aux(ma,i ,j,k+1)
aux3(j,ma,3) = aux(ma,i+1,j,k+1)
72 continue
endif
c
c # Store the value of i and k along this slice in the common block
c # comxyzt in case it is needed in the Riemann solver (for
c # variable coefficient problems)
c
icom = i
kcom = k
c
c # compute modifications fadd, gadd and hadd to fluxes along this
c # slice:
c
call flux3(2,maxm,meqn,maux,mwaves,mbc,my,
& q1d,dtdy1d,dtdz,dtdx,aux1,aux2,aux3,
& method,mthlim,
& faddm,faddp,gaddm,gaddp,haddm,haddp,cfl1d,
& work(i0wave),work(i0s),work(i0amdq),
& work(i0apdq),work(i0cqxx),
& work(i0bmamdq),work(i0bmapdq),
& work(i0bpamdq),work(i0bpapdq),
& work(i0cmamdq),work(i0cmapdq),
& work(i0cpamdq),work(i0cpapdq),
& work(i0cmamdq2),work(i0cmapdq2),
& work(i0cpamdq2),work(i0cpapdq2),
& work(i0bmcqxx),work(i0bpcqxx),
& work(i0cmcqxx),work(i0cpcqxx),
& work(i0bmcmamdq),work(i0bmcmapdq),
& work(i0bpcmamdq),work(i0bpcmapdq),
& work(i0bmcpamdq),work(i0bmcpapdq),
& work(i0bpcpamdq),work(i0bpcpapdq),
& work(i0slope),mslope,rpn3,rpt3)
c
cfl = dmax1(cfl,cfl1d)
c
c # Note that the roles of the flux updates are changed.
c # fadd - modifies the g-fluxes
c # gadd - modifies the h-fluxes
c # hadd - modifies the f-fluxes
do 80 j=1,my+1
do 80 m=1,meqn
gm(m,i,j,k) = gm(m,i,j,k) + faddm(j,m)
gp(m,i,j,k) = gp(m,i,j,k) + faddp(j,m)
80 continue
c
if (m3.gt.0) then
do 82 j=1,my+1
do 82 m=1,meqn
fm(m,i,j,k ) = fm(m,i,j,k ) + haddm(j,m,1,0)
fp(m,i,j,k ) = fp(m,i,j,k ) + haddp(j,m,1,0)
fm(m,i+1,j,k ) = fm(m,i+1,j,k ) + haddm(j,m,2,0)
fp(m,i+1,j,k ) = fp(m,i+1,j,k ) + haddp(j,m,2,0)
if (m4.gt.0) then
fm(m,i,j,k-1) = fm(m,i,j,k-1) + haddm(j,m,1,-1)
fp(m,i,j,k-1) = fp(m,i,j,k-1) + haddp(j,m,1,-1)
fm(m,i,j,k+1) = fm(m,i,j,k+1) + haddm(j,m,1,1)
fp(m,i,j,k+1) = fp(m,i,j,k+1) + haddp(j,m,1,1)
fm(m,i+1,j,k-1)=fm(m,i+1,j,k-1)+haddm(j,m,2,-1)
fp(m,i+1,j,k-1)=fp(m,i+1,j,k-1)+haddp(j,m,2,-1)
fm(m,i+1,j,k+1)=fm(m,i+1,j,k+1)+haddm(j,m,2,1)
fp(m,i+1,j,k+1)=fp(m,i+1,j,k+1)+haddp(j,m,2,1)
endif
82 continue
c
do 84 j=1,my+1
do 84 m=1,meqn
hm(m,i, j,k) = hm(m,i, j,k) + gaddm(j,m,1,0)
hp(m,i, j,k) = hp(m,i, j,k) + gaddp(j,m,1,0)
hm(m,i, j,k+1) = hm(m,i, j,k+1) + gaddm(j,m,2,0)
hp(m,i, j,k+1) = hp(m,i, j,k+1) + gaddp(j,m,2,0)
if (m4.gt.0) then
hm(m,i-1,j,k) = hm(m,i-1,j,k) + gaddm(j,m,1,-1)
hp(m,i-1,j,k) = hp(m,i-1,j,k) + gaddp(j,m,1,-1)
hm(m,i+1,j,k) = hm(m,i+1,j,k) + gaddm(j,m,1,1)
hp(m,i+1,j,k) = hp(m,i+1,j,k) + gaddp(j,m,1,1)
hm(m,i-1,j,k+1)=hm(m,i-1,j,k+1)+gaddm(j,m,2,-1)
hp(m,i-1,j,k+1)=hp(m,i-1,j,k+1)+gaddp(j,m,2,-1)
hm(m,i+1,j,k+1)=hm(m,i+1,j,k+1)+gaddm(j,m,2,1)
hp(m,i+1,j,k+1)=hp(m,i+1,j,k+1)+gaddp(j,m,2,1)
endif
84 continue
endif
c
if (method(1).ge.1.and.method(2).ge.3)
& call saverec3(2,maxm,maxmx,maxmy,maxmz,mvar,meqn,
& mbc,mx,my,mz,qold,work(i0qls),work(i0qrs),
& work(i0qbs),work(i0qts),work(i0qfs),work(i0qks),
& work(i0ql),work(i0qr))
c
100 continue
c
c
c # perform z sweeps
c ==================
c
c
do 150 j = 0, my+1
do 150 i = 0, mx+1
c
do 110 k = 1-mbc, mz+mbc
do 110 m=1,meqn
c # copy data along a slice into 1d array:
q1d(k,m) = qold(m,i,j,k)
110 continue
c
if (mcapa.gt.0) then
do 130 k = 1-mbc, mz+mbc
dtdz1d(k) = dtdz / aux(mcapa,i,j,k)
130 continue
endif
c
if (maux .gt. 0) then
do 131 k = 1-mbc, mz+mbc
do 131 ma=1,maux
aux1(k,ma,1) = aux(ma,i-1,j-1,k)
aux1(k,ma,2) = aux(ma,i-1,j ,k)
aux1(k,ma,3) = aux(ma,i-1,j+1,k)
aux2(k,ma,1) = aux(ma,i ,j-1,k)
aux2(k,ma,2) = aux(ma,i ,j ,k)
aux2(k,ma,3) = aux(ma,i ,j+1,k)
aux3(k,ma,1) = aux(ma,i+1,j-1,k)
aux3(k,ma,2) = aux(ma,i+1,j ,k)
aux3(k,ma,3) = aux(ma,i+1,j+1,k)
131 continue
endif
c
c # Store the value of i and j along this slice in the common block
c # comxyzt in case it is needed in the Riemann solver (for
c # variable coefficient problems)
c
icom = i
jcom = j
c
c # compute modifications fadd, gadd and hadd to fluxes along this
c # slice:
c
call flux3(3,maxm,meqn,maux,mwaves,mbc,mz,
& q1d,dtdz1d,dtdx,dtdy,aux1,aux2,aux3,
& method,mthlim,
& faddm,faddp,gaddm,gaddp,haddm,haddp,cfl1d,
& work(i0wave),work(i0s),work(i0amdq),
& work(i0apdq),work(i0cqxx),
& work(i0bmamdq),work(i0bmapdq),
& work(i0bpamdq),work(i0bpapdq),
& work(i0cmamdq),work(i0cmapdq),
& work(i0cpamdq),work(i0cpapdq),
& work(i0cmamdq2),work(i0cmapdq2),
& work(i0cpamdq2),work(i0cpapdq2),
& work(i0bmcqxx),work(i0bpcqxx),
& work(i0cmcqxx),work(i0cpcqxx),
& work(i0bmcmamdq),work(i0bmcmapdq),
& work(i0bpcmamdq),work(i0bpcmapdq),
& work(i0bmcpamdq),work(i0bmcpapdq),
& work(i0bpcpamdq),work(i0bpcpapdq),
& work(i0slope),mslope,rpn3,rpt3)
c
cfl = dmax1(cfl,cfl1d)
c
c # Note that the roles of the flux updates are changed.
c # fadd - modifies the h-fluxes
c # gadd - modifies the f-fluxes
c # hadd - modifies the g-fluxes
c
do 120 k=1,mz+1
do 120 m=1,meqn
hm(m,i,j,k) = hm(m,i,j,k) + faddm(k,m)
hp(m,i,j,k) = hp(m,i,j,k) + faddp(k,m)
120 continue
c
if (m3.gt.0) then
do 122 k=1,mz+1
do 122 m=1,meqn
fm(m,i,j, k) = fm(m,i,j, k) + gaddm(k,m,1,0)
fp(m,i,j, k) = fp(m,i,j, k) + gaddp(k,m,1,0)
fm(m,i+1,j, k) = fm(m,i+1,j, k) + gaddm(k,m,2,0)
fp(m,i+1,j, k) = fp(m,i+1,j, k) + gaddp(k,m,2,0)
if (m4.gt.0) then
fm(m,i,j-1,k) = fm(m,i,j-1,k) + gaddm(k,m,1,-1)
fp(m,i,j-1,k) = fp(m,i,j-1,k) + gaddp(k,m,1,-1)
fm(m,i,j+1,k) = fm(m,i,j+1,k) + gaddm(k,m,1,1)
fp(m,i,j+1,k) = fp(m,i,j+1,k) + gaddp(k,m,1,1)
fm(m,i+1,j-1,k)=fm(m,i+1,j-1,k)+gaddm(k,m,2,-1)
fp(m,i+1,j-1,k)=fp(m,i+1,j-1,k)+gaddp(k,m,2,-1)
fm(m,i+1,j+1,k)=fm(m,i+1,j+1,k)+gaddm(k,m,2,1)
fp(m,i+1,j+1,k)=fp(m,i+1,j+1,k)+gaddp(k,m,2,1)
endif
122 continue
c
do 124 k=1,mz+1
do 124 m=1,meqn
gm(m,i, j,k) = gm(m,i, j,k) + haddm(k,m,1,0)
gp(m,i, j,k) = gp(m,i, j,k) + haddp(k,m,1,0)
gm(m,i, j+1,k) = gm(m,i, j+1,k) + haddm(k,m,2,0)
gp(m,i, j+1,k) = gp(m,i, j+1,k) + haddp(k,m,2,0)
if (m4.gt.0) then
gm(m,i-1,j,k) = gm(m,i-1,j,k) + haddm(k,m,1,-1)
gp(m,i-1,j,k) = gp(m,i-1,j,k) + haddp(k,m,1,-1)
gm(m,i+1,j,k) = gm(m,i+1,j,k) + haddm(k,m,1,1)
gp(m,i+1,j,k) = gp(m,i+1,j,k) + haddp(k,m,1,1)
gm(m,i-1,j+1,k)=gm(m,i-1,j+1,k)+haddm(k,m,2,-1)
gp(m,i-1,j+1,k)=gp(m,i-1,j+1,k)+haddp(k,m,2,-1)
gm(m,i+1,j+1,k)=gm(m,i+1,j+1,k)+haddm(k,m,2,1)
gp(m,i+1,j+1,k)=gp(m,i+1,j+1,k)+haddp(k,m,2,1)
endif
124 continue
endif
c
if (method(1).ge.1.and.method(2).ge.3)
& call saverec3(3,maxm,maxmx,maxmy,maxmz,mvar,meqn,
& mbc,mx,my,mz,qold,work(i0qls),work(i0qrs),
& work(i0qbs),work(i0qts),work(i0qfs),work(i0qks),
& work(i0ql),work(i0qr))
c
150 continue
c
if (method(1).ge.1) then
if (method(2).le.2) then
call fmod3(maxm,maxmx,maxmy,maxmz,mvar,meqn,maux,mbc,mx,
& my,mz,qold,qold,qold,qold,qold,qold,qold,aux,
& dx,dy,dz,dt,method,cfl,fm,fp,gm,gp,hm,hp,
& q1d,aux2,faddm,faddp,0)
else
call fmod3(maxm,maxmx,maxmy,maxmz,mvar,meqn,maux,mbc,mx,
& my,mz,qold,work(i0qls),work(i0qrs),work(i0qbs),
& work(i0qts),work(i0qfs),work(i0qks),aux,
& dx,dy,dz,dt,method,cfl,fm,fp,gm,gp,hm,hp,
& q1d,aux2,faddm,faddp,0)
endif
endif
c
do 200 k=1,mz
do 200 j=1,my
do 200 i=1,mx
do 200 m=1,meqn
if (mcapa.gt.0) then
c # with capa array.
qold(m,i,j,k) = qold(m,i,j,k)
& - (dtdx * (fm(m,i+1,j,k) - fp(m,i,j,k))
& + dtdy * (gm(m,i,j+1,k) - gp(m,i,j,k))
& + dtdz * (hm(m,i,j,k+1) - hp(m,i,j,k))) /
& aux(mcapa,i,j,k)
c # no capa array. Standard flux differencing:
else
qold(m,i,j,k) = qold(m,i,j,k)
& - dtdx * (fm(m,i+1,j,k) - fp(m,i,j,k))
& - dtdy * (gm(m,i,j+1,k) - gp(m,i,j,k))
& - dtdz * (hm(m,i,j,k+1) - hp(m,i,j,k))
endif
c
200 continue
c
return
end