2d/integrator_ex/dimsp2ex.f

c
c
c     ==========================================================
      subroutine dimsp2(maxm,maxmx,maxmy,mvar,meqn,maux,mwaves,mbc,
     &                  mx,my,qold,aux,dx,dy,dt,method,mthlim,cfl,
     &		        fm,fp,gm,gp,faddm,faddp,gaddm,gaddp,
     &                  q1d,dtdx1d,dtdy1d,aux1,aux2,aux3,
     &                  work,mwork,rpn2,rpt2)
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
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 
c     # execution.
c
c     # See the flux2 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)
      dimension    fm(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
      dimension    fp(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
      dimension    gm(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
      dimension    gp(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+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)
      dimension gaddp(1-mbc:maxm+mbc, meqn, 2)
      dimension aux(maux, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
      dimension aux1(1-mbc:maxm+mbc, maux)
      dimension aux2(1-mbc:maxm+mbc, maux)
      dimension aux3(1-mbc:maxm+mbc, maux)
      dimension dtdx1d(1-mbc:maxm+mbc)
      dimension dtdy1d(1-mbc:maxm+mbc)
      dimension method(7),mthlim(mwaves)
      dimension work(mwork)
      external rpn2,rpt2
c
      mcapa = method(6)
c
      do 10 j=1-mbc,my+mbc
	 do 10 i=1-mbc,mx+mbc
            do 10 m=1,mvar
	       fm(m,i,j) = 0.d0
	       fp(m,i,j) = 0.d0
	       gm(m,i,j) = 0.d0
	       gp(m,i,j) = 0.d0
   10 continue
c
      cfl = 0.d0
      dtdx = dt/dx
      dtdy = dt/dy
c
      if (mpass.eq.0.or.mpass.eq.1) then
c
c     # Take a full time step in x-direction
c
         call step2ds(maxm,maxmx,maxmy,mvar,meqn,maux,mwaves,mbc,mx,my,
     &                qold,aux,dx,dy,dt,method,mthlim,cflx,
     &                fm,fp,gm,gp,faddm,faddp,gaddm,gaddp,
     &                q1d,dtdx1d,dtdy1d,aux1,aux2,aux3,
     &                work,mwork,rpn2,rpt2,1)
         cfl = cflx
c     
         do 110 j=1,my               
            do 110 i=1,mx
               do 110 m=1,meqn
                  if (mcapa.gt.0) then
                     qold(m,i,j) = qold(m,i,j) 
     &                    - dtdx * (fm(m,i+1,j) - fp(m,i,j)) / 
     &                    aux(mcapa,i,j)
c     # no capa array.  Standard flux differencing:
                  else
                     qold(m,i,j) = qold(m,i,j) 
     &                    - dtdx * (fm(m,i+1,j) - fp(m,i,j)) 
                  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 step2ds(maxm,maxmx,maxmy,mvar,meqn,maux,mwaves,mbc,mx,my,
     &                qold,aux,dx,dy,dt,method,mthlim,cfly,
     &                fm,fp,gm,gp,faddm,faddp,gaddm,gaddp,
     &                q1d,dtdx1d,dtdy1d,aux1,aux2,aux3,
     &                work,mwork,rpn2,rpt2,2)
         cfl = cfly
c
         do 120 j=1,my               
            do 120 i=1,mx
               do 120 m=1,meqn
                  if (mcapa.gt.0) then
                     qold(m,i,j) = qold(m,i,j) 
     &                    - dtdy * (gm(m,i,j+1) - gp(m,i,j)) / 
     &                    aux(mcapa,i,j)
c     # no capa array.  Standard flux differencing:
                  else
                     qold(m,i,j) = qold(m,i,j) 
     &                    - dtdy * (gm(m,i,j+1) - gp(m,i,j)) 
                  endif
 120     continue
c
      endif
c
      if (mpass.eq.0) cfl = dmax1(cflx,cfly)
c
      return
      end