c
c
c     ==================================================================
      subroutine flux3(ixyz,maxm,meqn,maux,mbc,mx,
     &                 q1d,dtdx1d,dtdy,dtdz,aux1,aux2,aux3,
     &                 faddm,faddp,gadd,hadd,cfl1d,
     &                 wave,s,amdq,apdq,cqxx,
     &                 bmamdq,bmapdq,bpamdq,bpapdq,
     &                 cmamdq,cmapdq,cpamdq,cpapdq,
     &                 cmamdq2,cmapdq2,cpamdq2,cpapdq2,
     &                 bmcqxxm,bmcqxxp,bpcqxxm,bpcqxxp,
     &                 cmcqxxm,cmcqxxp,cpcqxxm,cpcqxxp,
     &                 bmcmamdq,bmcmapdq,bpcmamdq,bpcmapdq,
     &                 bmcpamdq,bmcpapdq,bpcpamdq,bpcpapdq,
     &                 rpn3,rpt3, rptt3,mwaves,mcapa,method,mthlim,
     &                 use_fwaves)
c     ==================================================================
c
c     # clawpack routine ...  modified for AMRCLAW
c
c     # Compute the modification to fluxes f, g and h that are generated by
c     # all interfaces along a 1D slice of the 3D grid.
c     #    ixyz = 1  if it is a slice in x
c     #           2  if it is a slice in y
c     #           3  if it is a slice in z
c     # This value is passed into the Riemann solvers. The flux modifications
c     # go into the arrays fadd, gadd and hadd.  The notation is written
c     # assuming we are solving along a 1D slice in the x-direction.
c
c     # fadd(i,.) modifies F to the left of cell i
c     # gadd(i,.,1,slice) modifies G below cell i (in the z-direction)
c     # gadd(i,.,2,slice) modifies G above cell i
c     #                   The G flux in the surrounding slices may
c     #                   also be updated.
c     #                   slice  =  -1     The slice below in y-direction
c     #                   slice  =   0     The slice used in the 2D method
c     #                   slice  =   1     The slice above in y-direction
c     # hadd(i,.,1,slice) modifies H below cell i (in the y-direction)
c     # hadd(i,.,2,slice) modifies H above cell i
c     #                   The H flux in the surrounding slices may
c     #                   also be updated.
c     #                   slice  =  -1     The slice below in z-direction
c     #                   slice  =   0     The slice used in the 2D method
c     #                   slice  =   1     The slice above in z-direction
c     #
c     # The method used is specified by method(2) and method(3):
c
c        method(2) = 1 No correction waves
c                  = 2 if second order correction terms are to be added, with
c                      a flux limiter as specified by mthlim.  No transverse
c                      propagation of these waves.
c
c         method(3) specify how the transverse wave propagation
c         of the increment wave and the correction wave are performed.
c         Note that method(3) is given by a two digit number, in
c         contrast to what is the case for claw2. It is convenient
c         to define the scheme using the pair (method(2),method(3)).
c
c         method(3) <  0 Gives dimensional splitting using Godunov
c                        splitting, i.e. formally first order
c                        accurate.
c                      0 Gives the Donor cell method. No transverse
c                        propagation of neither the increment wave
c                        nor the correction wave.
c                   = 10 Transverse propagation of the increment wave
c                        as in 2D. Note that method (2,10) is
c                        unconditionally unstable.
c                   = 11 Corner transport upwind of the increment
c                        wave. Note that method (2,11) also is
c                        unconditionally unstable.
c                   = 20 Both the increment wave and the correction
c                        wave propagate as in the 2D case. Only to
c                        be used with method(2) = 2.
c                   = 21 Corner transport upwind of the increment wave,
c                        and the correction wave propagates as in 2D.
c                        Only to be used with method(2) = 2.
c                   = 22 3D propagation of both the increment wave and
c                        the correction wave. Only to be used with
c                        method(2) = 2.
c
c         Recommended settings:   First order schemes:
c                                       (1,10) Stable for CFL < 1/2
c                                       (1,11) Stable for CFL < 1
c                                 Second order schemes:
c                                        (2,20) Stable for CFL < 1/2
c                                        (2,22) Stable for CFL < 1
c
c         WARNING! The schemes (2,10), (2,11) are unconditionally
c                  unstable.
c
c                       ----------------------------------
c
c     Note that if method(6)=1 then the capa array comes into the second
c     order correction terms, and is already included in dtdx1d:
c     If ixyz = 1 then
c        dtdx1d(i) = dt/dx                      if method(6) = 0
c                  = dt/(dx*capa(i,jcom,kcom))  if method(6) = 1
c     If ixyz = 2 then
c        dtdx1d(j) = dt/dy                      if method(6) = 0
c                  = dt/(dy*capa(icom,j,kcom))  if method(6) = 1
c     If ixyz = 3 then
c        dtdx1d(k) = dt/dz                      if method(6) = 0
c                  = dt/(dz*capa(icom,jcom,k))  if method(6) = 1
c
c     Notation:
c        The jump in q (q1d(i,:)-q1d(i-1,:))  is split by rpn3 into
c            amdq =  the left-going flux difference  A^- Delta q
c            apdq = the right-going flux difference  A^+ Delta q
c        Each of these is split by rpt3 into
c            bmasdq = the down-going transverse flux difference B^- A^* Delta q
c            bpasdq =   the up-going transverse flux difference B^+ A^* Delta q
c        where A^* represents either A^- or A^+.
c
c        Finally, each bsasdq is split by rptt3 into :
c            cmbsasdq = C^- B^* A^* Dq
c            cpbsasdq = C^+ B^* A^* Dq
c
c
c      use amr_module
      implicit none

      external rpn3,rpt3, rptt3


      integer ixyz, maxm, meqn, maux, mbc, mx, mwaves, mcapa
      integer mthlim(mwaves), method(7)
      double precision cfl1d, dtdy, dtdz
      integer use_fwaves

      double precision     q1d(meqn,1-mbc:maxm+mbc)
      double precision    amdq(meqn,1-mbc:maxm+mbc)
      double precision    apdq(meqn,1-mbc:maxm+mbc)
      double precision  bmamdq(meqn,1-mbc:maxm+mbc)
      double precision  bmapdq(meqn,1-mbc:maxm+mbc)
      double precision  bpamdq(meqn,1-mbc:maxm+mbc)
      double precision  bpapdq(meqn,1-mbc:maxm+mbc)
      double precision    cqxx(meqn,1-mbc:maxm+mbc)
      double precision   faddm(meqn,1-mbc:maxm+mbc)
      double precision   faddp(meqn,1-mbc:maxm+mbc)
      double precision    gadd(meqn,1-mbc:maxm+mbc,2, -1:1)
      double precision    hadd(meqn,1-mbc:maxm+mbc,2, -1:1)
c
      double precision  cmamdq(meqn,1-mbc:maxm+mbc)
      double precision  cmapdq(meqn,1-mbc:maxm+mbc)
      double precision  cpamdq(meqn,1-mbc:maxm+mbc)
      double precision  cpapdq(meqn,1-mbc:maxm+mbc)
c
      double precision  cmamdq2(meqn,1-mbc:maxm+mbc)
      double precision  cmapdq2(meqn,1-mbc:maxm+mbc)
      double precision  cpamdq2(meqn,1-mbc:maxm+mbc)
      double precision  cpapdq2(meqn,1-mbc:maxm+mbc)
c
      double precision  bmcqxxm(meqn,1-mbc:maxm+mbc)
      double precision  bpcqxxm(meqn,1-mbc:maxm+mbc)
      double precision  bmcqxxp(meqn,1-mbc:maxm+mbc)
      double precision  bpcqxxp(meqn,1-mbc:maxm+mbc)
c
      double precision  cmcqxxm(meqn,1-mbc:maxm+mbc)
      double precision  cpcqxxm(meqn,1-mbc:maxm+mbc)
      double precision  cmcqxxp(meqn,1-mbc:maxm+mbc)
      double precision  cpcqxxp(meqn,1-mbc:maxm+mbc)

      double precision  bpcmamdq(meqn,1-mbc:maxm+mbc)
      double precision  bpcmapdq(meqn,1-mbc:maxm+mbc)
      double precision  bpcpamdq(meqn,1-mbc:maxm+mbc)
      double precision  bpcpapdq(meqn,1-mbc:maxm+mbc)
      double precision  bmcmamdq(meqn,1-mbc:maxm+mbc)
      double precision  bmcmapdq(meqn,1-mbc:maxm+mbc)
      double precision  bmcpamdq(meqn,1-mbc:maxm+mbc)
      double precision  bmcpapdq(meqn,1-mbc:maxm+mbc)
c
      double precision dtdx1d(1-mbc:maxm+mbc)
      double precision aux1(maux,1-mbc:maxm+mbc, 3)
      double precision aux2(maux,1-mbc:maxm+mbc, 3)
      double precision aux3(maux,1-mbc:maxm+mbc, 3)
c
      double precision    s(mwaves,1-mbc:maxm+mbc)
      double precision  wave(meqn,mwaves,1-mbc:maxm+mbc)
c
      logical limit

c      double precision dtcom, dxcom, dycom, dzcom, tcom
c      integer icom,jcom,kcom
c      common/comxyzt/dtcom,dxcom,dycom,dzcom,tcom,icom,jcom,kcom

      integer mw, m, i,m3, m4
      double precision dtdxave


      limit = .false.
      do mw = 1,mwaves
         if (mthlim(mw) .gt. 0) then 
            limit = .true.
         endif
      end do
c     
c     # initialize flux increments:
c     -----------------------------
c     
      do m = 1,meqn
         do i = 1-mbc, mx+mbc
            faddm(m,i) = 0.d0
            faddp(m,i) = 0.d0
            gadd(m,i,1,-1) = 0.d0
            gadd(m,i,1, 0) = 0.d0
            gadd(m,i,1, 1) = 0.d0
            hadd(m,i,1,-1) = 0.d0
            hadd(m,i,1, 0) = 0.d0
            hadd(m,i,1, 1) = 0.d0
            gadd(m,i,2,-1) = 0.d0
            gadd(m,i,2, 0) = 0.d0
            gadd(m,i,2, 1) = 0.d0
            hadd(m,i,2,-1) = 0.d0
            hadd(m,i,2, 0) = 0.d0
            hadd(m,i,2, 1) = 0.d0
         end do
      end do
     
c     # local method parameters
      if (method(3) .lt. 0) then
         m3 = -1
         m4 = 0
      else
         m3 = method(3)/10
         m4 = method(3) - 10*m3
      endif

c     -----------------------------------------------------------
c     # solve normal Riemann problem and compute Godunov updates
c     -----------------------------------------------------------
     
c     # aux2(1-mbc,1,2) is the start of a 1d array now used by rpn3
     
      if (maux > 0) then
         call rpn3(ixyz,maxm,meqn,mwaves,maux,mbc,mx,q1d,q1d,
     &           aux2(1,1-mbc,2),aux2(1,1-mbc,2),
     &           wave,s,amdq,apdq)
      else
         call rpn3(ixyz,maxm,meqn,mwaves,maux,mbc,mx,q1d,q1d,
     &           aux2,aux2,
     &           wave,s,amdq,apdq)
      endif
     
c     # Set fadd for the donor-cell upwind method (Godunov)
      do i = 1,mx+1
         do m = 1,meqn
            faddp(m,i) = faddp(m,i) - apdq(m,i)
            faddm(m,i) = faddm(m,i) + amdq(m,i)
         end do
      end do
     
c     # compute maximum wave speed for checking Courant number:
      cfl1d = 0.d0
      do mw = 1,mwaves
         do i = 1,mx+1
c          !  cfl1d = dmax1(cfl1d,dtdx1d(i)*dabs(s(mw,i))) OLD WAY
c          # if s>0 use dtdx1d(i) to compute CFL,
c          # if s<0 use dtdx1d(i-1) to compute CFL:
           cfl1d = dmax1(cfl1d,  dtdx1d(i)*s(mw,i), 
     &                    -dtdx1d(i-1)*s(mw,i))
         end do
      end do

     
      if (method(2).eq.1) then
        go to 130
      endif
    
c         # modify F fluxes for second order q_{xx} correction terms:
c         -----------------------------------------------------------
     
c     # apply limiter to waves:
      if (limit) then 
         call clawpack46_inlinelimiter(maxm,meqn,mwaves,
     &                  mbc,mx,wave,s,mthlim)
      endif
     
      if (use_fwaves .eq. 0) then
          do i = 1, mx+1
     
              dtdxave = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i))
     
              do m = 1,meqn
                  cqxx(m,i) = 0.d0
                  do mw = 1,mwaves
                      cqxx(m,i) = cqxx(m,i) + 0.5d0 * dabs(s(mw,i)) *
     &                     (1.d0 - dabs(s(mw,i))*dtdxave) * wave(m,mw,i)
                  end do
                  faddm(m,i) = faddm(m,i) + cqxx(m,i)
                  faddp(m,i) = faddp(m,i) + cqxx(m,i)
              enddo
          end do
      else
          do i = 1, mx+1

              dtdxave = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i)) 
c
              do m = 1,meqn
                  cqxx(m,i) = 0.d0
                  do mw = 1,mwaves
                      cqxx(m,i) = cqxx(m,i) + 0.5d0 * 
     &                    dsign(1.d0, s(mw,i)) *
     &              (1.d0 - dabs(s(mw,i))*dtdxave) * wave(m,mw,i)
                  end do
                  faddm(m,i) = faddm(m,i) + cqxx(m,i)
                  faddp(m,i) = faddp(m,i) + cqxx(m,i)
              end do
          end do
      endif
c     
 130  continue
c     
      if ( m3 .le. 0) then
          return
      endif
c     
c     --------------------------------------------
c     # TRANSVERSE PROPAGATION
c     --------------------------------------------
c     
c     # split the left-going flux difference into down-going and up-going
c     # flux differences (in the y-direction).
c     
      call rpt3(ixyz,2,1,maxm,meqn,mwaves,maux,mbc,mx,
     &             q1d,q1d,aux1,aux2,aux3,amdq,bmamdq,bpamdq)
c     
c     # split the right-going flux difference into down-going and up-going
c     # flux differences (in the y-direction).
c     
      call rpt3(ixyz,2,2,maxm,meqn,mwaves,maux,mbc,mx,
     &             q1d,q1d,aux1,aux2,aux3,apdq,bmapdq,bpapdq)
c     
c     # split the left-going flux difference into down-going and up-going
c     # flux differences (in the z-direction).
c     
      call rpt3(ixyz,3,1,maxm,meqn,mwaves,maux,mbc,mx,
     &             q1d,q1d,aux1,aux2,aux3,amdq,cmamdq,cpamdq)
c     
c     # split the right-going flux difference into down-going and up-going
c     # flux differences (in the y-direction).
c     
      call rpt3(ixyz,3,2,maxm,meqn,mwaves,maux,mbc,mx,
     &             q1d,q1d,aux1,aux2,aux3,apdq,cmapdq,cpapdq)
c     
c     # Split the correction wave into transverse propagating waves
c     # in the y-direction and z-direction.
c     
      if (m3 .eq. 2) then
          if (maux .gt. 0) then
c             # The corrections cqxx affect both cell i-1 to left and cell i
c             # to right of interface.  Transverse splitting will affect
c             # fluxes on both sides.
c             # If there are aux arrays, then we must split cqxx twice in
c             # each transverse direction, once with imp=1 and once with imp=2:

c             # imp = 1 or 2 is used to indicate whether we are propagating
c             # amdq or apdq, i.e. cqxxm or cqxxp


c             # in the y-like direction with imp=1
              call rpt3(ixyz,2,1,maxm,meqn,mwaves,maux,mbc,mx,
     &               q1d,q1d,aux1,aux2,aux3,cqxx,bmcqxxm,bpcqxxm)

c             # in the y-like direction with imp=2
              call rpt3(ixyz,2,2,maxm,meqn,mwaves,maux,mbc,mx,
     &               q1d,q1d,aux1,aux2,aux3,cqxx,bmcqxxp,bpcqxxp)

c             # in the z-like direction with imp=1
              call rpt3(ixyz,3,1,maxm,meqn,mwaves,maux,mbc,mx,
     &               q1d,q1d,aux1,aux2,aux3,cqxx,cmcqxxm,cpcqxxm)

c             # in the z-like direction with imp=2
              call rpt3(ixyz,3,2,maxm,meqn,mwaves,maux,mbc,mx,
     &               q1d,q1d,aux1,aux2,aux3,cqxx,cmcqxxp,cpcqxxp)
          else
c             # aux arrays aren't being used, so we only need to split
c             # cqxx once in each transverse direction and the same result can
c             # presumably be used to left and right.
c             # Set imp = 0 since this shouldn't be needed in rpt3 in this case.

c             # in the y-like direction
              call rpt3(ixyz,2,0,maxm,meqn,mwaves,maux,mbc,mx,
     &               q1d,q1d,aux1,aux2,aux3,cqxx,bmcqxxm,bpcqxxm)

c             # in the z-like direction
              call rpt3(ixyz,3,0,maxm,meqn,mwaves,maux,mbc,mx,
     &                   q1d,q1d,aux1,aux2,aux3,cqxx,cmcqxxm,cpcqxxm)

c             # use the same splitting to left and right:
              do i = 0,mx+2
                  do m = 1,meqn
                      bmcqxxp(m,i) = bmcqxxm(m,i)
                      bpcqxxp(m,i) = bpcqxxm(m,i)
                      cmcqxxp(m,i) = cmcqxxm(m,i)
                      cpcqxxp(m,i) = cpcqxxm(m,i)
                  enddo
              enddo
          endif
      endif

c     --------------------------------------------
c     # modify G fluxes in the y-like direction
c     --------------------------------------------
c     
c     # If the correction wave also propagates in a 3D sense, incorporate
c     # cpcqxx,... into cmamdq, cpamdq, ... so that it is split also.
c     
      if (m4 .eq. 1) then
          do i = 0, mx+2
              do m = 1, meqn
                  cpapdq2(m,i) = cpapdq(m,i)
                  cpamdq2(m,i) = cpamdq(m,i)
                  cmapdq2(m,i) = cmapdq(m,i)
                  cmamdq2(m,i) = cmamdq(m,i)
              end do
          end do
      else if (m4 .eq. 2) then
          do i = 0, mx+2
              do m = 1, meqn
                  cpapdq2(m,i) = cpapdq(m,i) - 3.d0*cpcqxxp(m,i)
                  cpamdq2(m,i) = cpamdq(m,i) + 3.d0*cpcqxxm(m,i)
                  cmapdq2(m,i) = cmapdq(m,i) - 3.d0*cmcqxxp(m,i)
                  cmamdq2(m,i) = cmamdq(m,i) + 3.d0*cmcqxxm(m,i)
              end do
          end do
      endif
c     
c     # The transverse flux differences in the z-direction are split
c     # into waves propagating in the y-direction. If m4 = 2,
c     # then the transverse propagating correction waves in the z-direction
c     # are also split. This yields terms of the form BCAu_{xzy} and
c     # BCAAu_{xxzy}.
c     
      if (m4 .gt. 0) then
          call rptt3(ixyz,2,2,2,maxm,meqn,mwaves,maux,mbc,mx,
     &                q1d,q1d,aux1,aux2,aux3,cpapdq2,
     &                bmcpapdq,bpcpapdq)
          call rptt3(ixyz,2,1,2,maxm,meqn,mwaves,maux,mbc,mx,
     &                q1d,q1d,aux1,aux2,aux3,cpamdq2,
     &                bmcpamdq,bpcpamdq)
          call rptt3(ixyz,2,2,1,maxm,meqn,mwaves,maux,mbc,mx,
     &                q1d,q1d,aux1,aux2,aux3,cmapdq2,
     &                bmcmapdq,bpcmapdq)
          call rptt3(ixyz,2,1,1,maxm,meqn,mwaves,maux,mbc,mx,
     &                q1d,q1d,aux1,aux2,aux3,cmamdq2,
     &                bmcmamdq,bpcmamdq)
      endif
c     
c     -----------------------------
c     # The updates for G fluxes :
c     -----------------------------
c     
      ig_loop : do i = 1, mx+1
          meqng_loop : do m = 1,meqn
     
c             # Transverse propagation of the increment waves
c             # between cells sharing interfaces, i.e. the 2D approach.
c             # Yields BAu_{xy}.
     
              gadd(m,i-1,1,0) = gadd(m,i-1,1,0)
     &                           - 0.5d0*dtdx1d(i-1)*bmamdq(m,i)
              gadd(m,i-1,2,0) = gadd(m,i-1,2,0)
     &                           - 0.5d0*dtdx1d(i-1)*bpamdq(m,i)
              gadd(m,i,1,0)   = gadd(m,i,1,0)
     &                           - 0.5d0*dtdx1d(i)*bmapdq(m,i)
              gadd(m,i,2,0)   = gadd(m,i,2,0)
     &                           - 0.5d0*dtdx1d(i)*bpapdq(m,i)
c     
c             # Transverse propagation of the increment wave (and the
c             # correction wave if m4=2) between cells
c             # only having a corner or edge in common. Yields terms of the
c             # BCAu_{xzy} and BCAAu_{xxzy}.
c     
              if (m4 .gt. 0) then
c     
                  gadd(m,i,2,0) = gadd(m,i,2,0)
     &                            + (1.d0/6.d0)*dtdx1d(i-1)*dtdz
     &                            * (bpcpapdq(m,i) - bpcmapdq(m,i))
                  gadd(m,i,1,0) = gadd(m,i,1,0)
     &                            + (1.d0/6.d0)*dtdx1d(i)*dtdz
     &                            * (bmcpapdq(m,i) - bmcmapdq(m,i))
                  gadd(m,i,2,1) = gadd(m,i,2,1)
     &                            - (1.d0/6.d0)*dtdx1d(i)*dtdz
     &                            * bpcpapdq(m,i)
                  gadd(m,i,1,1) = gadd(m,i,1,1)
     &                             - (1.d0/6.d0)*dtdx1d(i)*dtdz
     &                             * bmcpapdq(m,i)
                  gadd(m,i,2,-1) = gadd(m,i,2,-1)
     &                             + (1.d0/6.d0)*dtdx1d(i)*dtdz
     &                             * bpcmapdq(m,i)
                  gadd(m,i,1,-1) = gadd(m,i,1,-1)
     &                             + (1.d0/6.d0)*dtdx1d(i)*dtdz
     &                             * bmcmapdq(m,i)     
                  gadd(m,i-1,2,0) = gadd(m,i-1,2,0)
     &                             + (1.d0/6.d0)*dtdx1d(i-1)*dtdz
     &                             * (bpcpamdq(m,i) - bpcmamdq(m,i))
                  gadd(m,i-1,1,0) = gadd(m,i-1,1,0)
     &                             + (1.d0/6.d0)*dtdx1d(i-1)*dtdz
     &                             * (bmcpamdq(m,i) - bmcmamdq(m,i))
                  gadd(m,i-1,2,1) = gadd(m,i-1,2,1)
     &                             - (1.d0/6.d0)*dtdx1d(i-1)*dtdz
     &                             * bpcpamdq(m,i)
                  gadd(m,i-1,1,1) = gadd(m,i-1,1,1)
     &                             - (1.d0/6.d0)*dtdx1d(i-1)*dtdz
     &                             * bmcpamdq(m,i)
                  gadd(m,i-1,2,-1) = gadd(m,i-1,2,-1)
     &                             + (1.d0/6.d0)*dtdx1d(i-1)*dtdz
     &                             * bpcmamdq(m,i)
                  gadd(m,i-1,1,-1) = gadd(m,i-1,1,-1)
     &                             + (1.d0/6.d0)*dtdx1d(i-1)*dtdz
     &                             * bmcmamdq(m,i)
c     
              endif
c     
c             # Transverse propagation of the correction wave between
c             # cells sharing faces. This gives BAAu_{xxy}.
c     
              if (m3 .ge. 2) then
                  gadd(m,i,2,0)   = gadd(m,i,2,0)
     &                              + dtdx1d(i)*bpcqxxp(m,i)
                  gadd(m,i,1,0)   = gadd(m,i,1,0)
     &                              + dtdx1d(i)*bmcqxxp(m,i)
                  gadd(m,i-1,2,0) = gadd(m,i-1,2,0)
     &                              - dtdx1d(i-1)*bpcqxxm(m,i)
                  gadd(m,i-1,1,0) = gadd(m,i-1,1,0)
     &                              - dtdx1d(i-1)*bmcqxxm(m,i)
              endif
          end do meqng_loop
      end do ig_loop
c     
c     
c     --------------------------------------------
c     # modify H fluxes in the z-like direction
c     --------------------------------------------
c     
c     # If the correction wave also propagates in a 3D sense, incorporate
c     # cqxx into bmamdq, bpamdq, ... so that is is split also.
c     
      if (m4 .eq. 2) then
          do i = 0, mx+2
              do m = 1, meqn
                  bpapdq(m,i) = bpapdq(m,i) - 3.d0*bpcqxxp(m,i)
                  bpamdq(m,i) = bpamdq(m,i) + 3.d0*bpcqxxm(m,i)
                  bmapdq(m,i) = bmapdq(m,i) - 3.d0*bmcqxxp(m,i)
                  bmamdq(m,i) = bmamdq(m,i) + 3.d0*bmcqxxm(m,i)
              end do
         end do
      endif
c     
c     # The transverse flux differences in the y-direction are split
c     # into waves propagating in the z-direction. If m4 = 2,
c     # then the transverse propagating correction waves in the y-direction
c     # are also split. This yields terms of the form BCAu_{xzy} and
c     # BCAAu_{xxzy}.
c     
c     # Note that the output to rptt3 below should logically be named
c     # cmbsasdq and cpbsasdq rather than bmcsasdq and bpcsasdq, but
c     # we are re-using the previous storage rather than requiring new arrays.
c     
      if (m4 .gt. 0) then
           call rptt3(ixyz,3,2,2,maxm,meqn,mwaves,maux,mbc,mx,
     &                q1d,q1d,aux1,aux2,aux3,bpapdq,bmcpapdq,bpcpapdq)
           call rptt3(ixyz,3,1,2,maxm,meqn,mwaves,maux,mbc,mx,
     &                q1d,q1d,aux1,aux2,aux3,bpamdq,bmcpamdq,bpcpamdq)
           call rptt3(ixyz,3,2,1,maxm,meqn,mwaves,maux,mbc,mx,
     &                q1d,q1d,aux1,aux2,aux3,bmapdq,bmcmapdq,bpcmapdq)
           call rptt3(ixyz,3,1,1,maxm,meqn,mwaves,maux,mbc,mx,
     &                q1d,q1d,aux1,aux2,aux3,bmamdq,bmcmamdq,bpcmamdq)
      endif
c     
c     -----------------------------
c     # The updates for H fluxes :
c     -----------------------------
c     
      ih_loop : do i = 1, mx+1
          meqnh_loop : do m = 1,meqn
c     
c             # Transverse propagation of the increment waves
c             # between cells sharing interfaces, i.e. the 2D approach.
c             # Yields CAu_{xy}.
c     
              hadd(m,i-1,1,0) = hadd(m,i-1,1,0)
     &                         - 0.5d0*dtdx1d(i-1)*cmamdq(m,i)
              hadd(m,i-1,2,0) = hadd(m,i-1,2,0)
     &                         - 0.5d0*dtdx1d(i-1)*cpamdq(m,i)
              hadd(m,i,1,0)   = hadd(m,i,1,0)
     &                         - 0.5d0*dtdx1d(i)*cmapdq(m,i)
              hadd(m,i,2,0)   = hadd(m,i,2,0)
     &                         - 0.5d0*dtdx1d(i)*cpapdq(m,i)
c     
c             # Transverse propagation of the increment wave (and the
c             # correction wave if m4=2) between cells
c             # only having a corner or edge in common. Yields terms of the
c             # CBAu_{xzy} and CBAAu_{xxzy}.
c     
              if (m4 .gt. 0) then
c     
                  hadd(m,i,2,0)  = hadd(m,i,2,0)
     &                        + (1.d0/6.d0)*dtdx1d(i)*dtdy
     &                        * (bpcpapdq(m,i) - bpcmapdq(m,i))
                  hadd(m,i,1,0)  = hadd(m,i,1,0)
     &                        + (1.d0/6.d0)*dtdx1d(i)*dtdy
     &                        * (bmcpapdq(m,i) - bmcmapdq(m,i))
                  hadd(m,i,2,1)  = hadd(m,i,2,1)
     &                        - (1.d0/6.d0)*dtdx1d(i)*dtdy
     &                        * bpcpapdq(m,i)
                  hadd(m,i,1,1)  = hadd(m,i,1,1)
     &                        - (1.d0/6.d0)*dtdx1d(i)*dtdy
     &                        * bmcpapdq(m,i)
                  hadd(m,i,2,-1) = hadd(m,i,2,-1)
     &                        + (1.d0/6.d0)*dtdx1d(i)*dtdy
     &                        * bpcmapdq(m,i)
                  hadd(m,i,1,-1) = hadd(m,i,1,-1)
     &                        + (1.d0/6.d0)*dtdx1d(i)*dtdy
     &                        * bmcmapdq(m,i)     
                  hadd(m,i-1,2,0)  = hadd(m,i-1,2,0)
     &                        + (1.d0/6.d0)*dtdx1d(i-1)*dtdy
     &                        * (bpcpamdq(m,i) - bpcmamdq(m,i))
                  hadd(m,i-1,1,0)  = hadd(m,i-1,1,0)
     &                        + (1.d0/6.d0)*dtdx1d(i-1)*dtdy
     &                        * (bmcpamdq(m,i) - bmcmamdq(m,i))
                  hadd(m,i-1,2,1)  = hadd(m,i-1,2,1)
     &                        - (1.d0/6.d0)*dtdx1d(i-1)*dtdy
     &                        * bpcpamdq(m,i)
                  hadd(m,i-1,1,1)  = hadd(m,i-1,1,1)
     &                        - (1.d0/6.d0)*dtdx1d(i-1)*dtdy
     &                        * bmcpamdq(m,i)
                  hadd(m,i-1,2,-1) = hadd(m,i-1,2,-1)
     &                        + (1.d0/6.d0)*dtdx1d(i-1)*dtdy
     &                        * bpcmamdq(m,i)
                  hadd(m,i-1,1,-1) = hadd(m,i-1,1,-1)
     &                        + (1.d0/6.d0)*dtdx1d(i-1)*dtdy
     &                        * bmcmamdq(m,i)
c     
              endif
c     
c             # Transverse propagation of the correction wave between
c             # cells sharing faces. This gives CAAu_{xxy}.
c     
              if (m3 .ge. 2) then
                  hadd(m,i,2,0)   = hadd(m,i,2,0)
     &                               + dtdx1d(i)*cpcqxxp(m,i)
                  hadd(m,i,1,0)   = hadd(m,i,1,0)
     &                               + dtdx1d(i)*cmcqxxp(m,i)
                  hadd(m,i-1,2,0) = hadd(m,i-1,2,0)
     &                               - dtdx1d(i-1)*cpcqxxm(m,i)
                  hadd(m,i-1,1,0) = hadd(m,i-1,1,0)
     &                               - dtdx1d(i-1)*cmcqxxm(m,i)
              endif
          end do meqnh_loop
      end do ih_loop
      return
      end
