AMROC Main     Blockstructured Adaptive Mesh Refinement in object-oriented C++


Main Page   Class Hierarchy   Compound List   File List  

3d/integrator_ex/dimsp3ex.f

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,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)
         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)
         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)
         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


Quickstart     Users Guide     Programmers Reference     Installation      Examples     Download



AMROC Main      Home      Contact
last update: 06/01/04