3d/integrator_ex/step3ex.f

c
c     ==================================================================
      subroutine step3(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     # Main entry from AMROC. Selects between the wave propagation method
c     # of Randall J. LeVeque and Godunov's dimensional splitting.
c
      implicit double precision (a-h,o-z)
      include "call.i"
      common /rpnflx/ mrpnflx
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
      if( method(3) .ge. 0 )then
c
c            # unsplit version
c             
         call 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
      else
c           # dimensional splitting (fractional steps)
c
         call 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
      endif
c
      if (mrpnflx.ne.0) then
         if (method(2).eq.2) then
            write(6,*) '*** Riemann solver returns fluxes.'
            write(6,*) '*** Wave limiting not possible.'
            write(6,*) '*** Set method(2)>=3 for slope limiting.'
            stop 
         endif
      else
         if (method(2).ge.3) then
            write(6,*) '*** Riemann solver returns flux differences.'
            write(6,*) '*** Slope limiting not possible.'
            write(6,*) '*** Set method(2)=2 for wave limiting.'
            stop 
         endif
      endif
c
      return
      end