c
c
c ==================================================================
subroutine dimsp3(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 qold, using dimensional splitting.
c #
c # method(3) < 0 gives Godunov splitting:
c # time step dt in x-direction
c # time step dt in y-direction
c # time step dt in z-direction
c #
c # The common variable mpass>0 allows the update of a single direction.
c # This is necessary for method(3) = -2. This setting forces AMROC to
c # sychronize the ghost cells after each directional sweep.
c # While method(3) = -1 and method(3) = -2 give identical results on a
c # single node, method(3) = -2 should always be used for parallel execution.
c
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
mcapa = method(6)
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
cfl = 0.d0
dtdx = dt/dx
dtdy = dt/dy
dtdz = dt/dz
c
if (mpass.eq.0.or.mpass.eq.1) then
c
c # Take a full time step in x-direction
c
call step3ds(maxm,maxmx,maxmy,maxmz,mvar,meqn,
& maux,mwaves,mbc,mx,my,mz,
& qold,aux,dx,dy,dz,dt,method,mthlim,cflx,
& fm,fp,gm,gp,hm,hp,
& faddm,faddp,gaddm,gaddp,haddm,haddp,
& q1d,dtdx1d,dtdy1d,dtdz1d,
& aux1,aux2,aux3,work,mwork,rpn3,rpt3,1)
c
cfl = cflx
c
do 110 k=1,mz
do 110 j=1,my
do 110 i=1,mx
do 110 m=1,meqn
if (mcapa.gt.0) then
qold(m,i,j,k) = qold(m,i,j,k)
& - dtdx * (fm(m,i+1,j,k) - fp(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))
endif
110 continue
c
endif
c
if (mpass.eq.0.or.mpass.eq.2) then
c
c # Take full step in y-direction
c
call step3ds(maxm,maxmx,maxmy,maxmz,mvar,meqn,
& maux,mwaves,mbc,mx,my,mz,
& qold,aux,dx,dy,dz,dt,method,mthlim,cfly,
& fm,fp,gm,gp,hm,hp,
& faddm,faddp,gaddm,gaddp,haddm,haddp,
& q1d,dtdx1d,dtdy1d,dtdz1d,
& aux1,aux2,aux3,work,mwork,rpn3,rpt3,2)
c
cfl = cfly
c
do 120 k=1,mz
do 120 j=1,my
do 120 i=1,mx
do 120 m=1,meqn
if (mcapa.gt.0) then
qold(m,i,j,k) = qold(m,i,j,k)
& - dtdy * (gm(m,i,j+1,k) - gp(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)
& - dtdy * (gm(m,i,j+1,k) - gp(m,i,j,k))
endif
120 continue
c
endif
c
if (mpass.eq.0.or.mpass.eq.3) then
c
c # Take full step in z-direction
c
call step3ds(maxm,maxmx,maxmy,maxmz,mvar,meqn,
& maux,mwaves,mbc,mx,my,mz,
& qold,aux,dx,dy,dz,dt,method,mthlim,cflz,
& fm,fp,gm,gp,hm,hp,
& faddm,faddp,gaddm,gaddp,haddm,haddp,
& q1d,dtdx1d,dtdy1d,dtdz1d,
& aux1,aux2,aux3,work,mwork,rpn3,rpt3,3)
c
cfl = cflz
c
do 130 k=1,mz
do 130 j=1,my
do 130 i=1,mx
do 130 m=1,meqn
if (mcapa.gt.0) then
qold(m,i,j,k) = qold(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)
& - dtdz * (hm(m,i,j,k+1) - hp(m,i,j,k))
endif
130 continue
c
endif
c
if (mpass.eq.0) cfl = dmax1(cflx,cfly,cflz)
c
return
end