************************************************************************
      function xup (f)
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, xup
!hpf$ align with TheCube:: f, xup

      c = 3./256.
      b = -25./256.
      a = .5-b-c

      xup(3:mx-3,:,:) = (
     &               a*(f(3:mx-3,:,:) + f(4:mx-2,:,:)) +
     &               b*(f(2:mx-4,:,:) + f(5:mx-1,:,:)) +
     &               c*(f(1:mx-5,:,:) + f(6:mx  ,:,:)))
      xup(mx-2,:,:) = (
     &               a*(f(mx-2,:,:) + f(mx-1,:,:)) +
     &               b*(f(mx-3,:,:) + f(mx,:,:)) +
     &               c*(f(mx-4,:,:) + f(1,:,:)))
      xup(mx-1,:,:) = (
     &               a*(f(mx-1,:,:) + f(mx,:,:)) +
     &               b*(f(mx-2,:,:) + f(1,:,:)) +
     &               c*(f(mx-3,:,:) + f(2,:,:)))
      xup(mx ,:,:) = (
     &               a*(f(mx  ,:,:) + f(1,:,:)) +
     &               b*(f(mx-1,:,:) + f(2,:,:)) +
     &               c*(f(mx-2,:,:) + f(3,:,:)))
      xup(1  ,:,:) = (
     &               a*(f(1   ,:,:) + f(2,:,:)) +
     &               b*(f(mx  ,:,:) + f(3,:,:)) +
     &               c*(f(mx-1,:,:) + f(4,:,:)))
      xup(2  ,:,:) = (
     &               a*(f(2   ,:,:) + f(3,:,:)) +
     &               b*(f(1   ,:,:) + f(4,:,:)) +
     &               c*(f(mx  ,:,:) + f(5,:,:)))

*
      end
************************************************************************
      function yup (f)
*
*  f is centered on (i-.5,j,k), fifth order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, yup
!hpf$ align with TheCube:: f, yup
*
*-----------------------------------------------------------------------
      c = 3./256.
      b = -25./256.
      a = .5-b-c
*

      yup(:,3:my-3,:) = (
     &                 a*(f(:,3:my-3,:) + f(:,4:my-2,:)) +
     &                 b*(f(:,2:my-4,:) + f(:,5:my-1,:)) +
     &                 c*(f(:,1:my-5,:) + f(:,6:my,:)))
      yup(:,2,:) = (
     &                 a*(f(:,2 ,:) + f(:,3,:)) +
     &                 b*(f(:,1 ,:) + f(:,4,:)) +
     &                 c*(f(:,my,:) + f(:,5,:)))
      yup(:,1,:) = (
     &                 a*(f(:,1 ,:) + f(:,2,:)) +
     &                 b*(f(:,my,:) + f(:,3,:)) +
     &                 c*(f(:,my-1,:) + f(:,4,:)))
      yup(:,my,:) = (
     &                 a*(f(:,my  ,:) + f(:,1,:)) +
     &                 b*(f(:,my-1,:) + f(:,2,:)) +
     &                 c*(f(:,my-2,:) + f(:,3,:)))
      yup(:,my-1,:) = (
     &                 a*(f(:,my-1,:) + f(:,my,:)) +
     &                 b*(f(:,my-2,:) + f(:,1,:)) +
     &                 c*(f(:,my-3,:) + f(:,2,:)))
      yup(:,my-2,:) = (
     &                 a*(f(:,my-2,:) + f(:,my-1,:)) +
     &                 b*(f(:,my-3,:) + f(:,my,:)) +
     &                 c*(f(:,my-4,:) + f(:,1,:)))
*
      end
************************************************************************
      function zup (f)
*
*  f is centered on (i-.5,j,k), fifth order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, zup
!hpf$ align with TheCube:: f, zup
*
*-----------------------------------------------------------------------
      c = 3./256.
      b = -25./256.
      a = .5-b-c
*
      zup = (a*(cshift(f,dim=3,shift=+1) + f)
     &   + b*(cshift(f,dim=3,shift=+2) + cshift(f,dim=3,shift=-1))
     &   + c*(cshift(f,dim=3,shift=+3) + cshift(f,dim=3,shift=-2))
     &     )

      end
************************************************************************
      function xdn (f)
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, xdn
!hpf$ align with TheCube:: f, xdn

      c = 3./256.
      b = -25./256.
      a = .5-b-c

      xdn(4:mx-2  ,:,:) = (
     &               a*(f(3:mx-3,:,:) + f(4:mx-2,:,:)) +
     &               b*(f(2:mx-4,:,:) + f(5:mx-1,:,:)) +
     &               c*(f(1:mx-5,:,:) + f(6:mx  ,:,:)))
      xdn(mx-1,:,:) = (
     &               a*(f(mx-2,:,:) + f(mx-1,:,:)) +
     &               b*(f(mx-3,:,:) + f(mx  ,:,:)) +
     &               c*(f(mx-4,:,:) + f(1   ,:,:)))
      xdn(mx  ,:,:) = (
     &               a*(f(mx-1,:,:) + f(mx,:,:)) +
     &               b*(f(mx-2,:,:) + f(1,:,:)) +
     &               c*(f(mx-3,:,:) + f(2,:,:)))
      xdn(1   ,:,:) = (
     &               a*(f(mx  ,:,:) + f(1,:,:)) +
     &               b*(f(mx-1,:,:) + f(2,:,:)) +
     &               c*(f(mx-2,:,:) + f(3,:,:)))
      xdn(2  ,:,:) = (
     &               a*(f(1   ,:,:) + f(2,:,:)) +
     &               b*(f(mx  ,:,:) + f(3,:,:)) +
     &               c*(f(mx-1,:,:) + f(4,:,:)))
      xdn(3  ,:,:) = (
     &               a*(f(2   ,:,:) + f(3,:,:)) +
     &               b*(f(1   ,:,:) + f(4,:,:)) +
     &               c*(f(mx  ,:,:) + f(5,:,:)))
*
      end
************************************************************************
      function ydn (f)
*
*  f is centered on (i-.5,j,k), fifth order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, ydn
!hpf$ align with TheCube:: f, ydn
*
*-----------------------------------------------------------------------
      c = 3./256.
      b = -25./256.
      a = .5-b-c
*

      ydn(:,4:my-2,:) = (
     &                 a*(f(:,3:my-3,:) + f(:,4:my-2,:)) +
     &                 b*(f(:,2:my-4,:) + f(:,5:my-1,:)) +
     &                 c*(f(:,1:my-5,:) + f(:,6:my,:)))
      ydn(:,3,:) = (
     &                 a*(f(:,2 ,:) + f(:,3,:)) +
     &                 b*(f(:,1 ,:) + f(:,4,:)) +
     &                 c*(f(:,my,:) + f(:,5,:)))
      ydn(:,2,:) = (
     &                 a*(f(:,1 ,:) + f(:,2,:)) +
     &                 b*(f(:,my,:) + f(:,3,:)) +
     &                 c*(f(:,my-1,:) + f(:,4,:)))
      ydn(:,1,:) = (
     &                 a*(f(:,my  ,:) + f(:,1,:)) +
     &                 b*(f(:,my-1,:) + f(:,2,:)) +
     &                 c*(f(:,my-2,:) + f(:,3,:)))
      ydn(:,my,:) = (
     &                 a*(f(:,my-1,:) + f(:,my,:)) +
     &                 b*(f(:,my-2,:) + f(:,1,:)) +
     &                 c*(f(:,my-3,:) + f(:,2,:)))
      ydn(:,my-1,:) = (
     &                 a*(f(:,my-2,:) + f(:,my-1,:)) +
     &                 b*(f(:,my-3,:) + f(:,my,:)) +
     &                 c*(f(:,my-4,:) + f(:,1,:)))
*
      end
************************************************************************
      function zdn (f)
*
*  f is centered on (i-.5,j,k), fifth order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, zdn
!hpf$ align with TheCube:: f, zdn
*
*-----------------------------------------------------------------------
      c = 3./256.
      b = -25./256.
      a = .5-b-c
*
      zdn = (a*(cshift(f,dim=3,shift=-1) + f)
     &   + b*(cshift(f,dim=3,shift=-2) + cshift(f,dim=3,shift=1))
     &   + c*(cshift(f,dim=3,shift=-3) + cshift(f,dim=3,shift=2))
     &     )
*
      end
************************************************************************
      function ddxup (f)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddxup
!hpf$ align with TheCube:: f, ddxup

      c = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
      b = (-1.-120.*c)/24.
      a = (1.-3.*b-5.*c)/dx
      b = b/dx
      c = c/dx

      ddxup(3:mx-3,:,:) = (
     &               a*(f(4:mx-2,:,:) - f(3:mx-3,:,:)) +
     &               b*(f(5:mx-1,:,:) - f(2:mx-4,:,:)) +
     &               c*(f(6:mx  ,:,:) - f(1:mx-5,:,:)))
      ddxup(mx-2 ,:,:) = (
     &               a*(f(mx-1,:,:) - f(mx-2,:,:)) +
     &               b*(f(mx  ,:,:) - f(mx-3,:,:)) +
     &               c*(f(1   ,:,:) - f(mx-4,:,:)))
      ddxup(mx-1 ,:,:) = (
     &               a*(f(mx  ,:,:) - f(mx-1,:,:)) +
     &               b*(f(1   ,:,:) - f(mx-2,:,:)) +
     &               c*(f(2   ,:,:) - f(mx-3,:,:)))
      ddxup(mx  ,:,:) = (
     &               a*(f(1   ,:,:) - f(mx  ,:,:)) +
     &               b*(f(2   ,:,:) - f(mx-1,:,:)) +
     &               c*(f(3   ,:,:) - f(mx-2,:,:)))
      ddxup(1   ,:,:) = (
     &               a*(f(2   ,:,:) - f(1   ,:,:)) +
     &               b*(f(3   ,:,:) - f(mx  ,:,:)) +
     &               c*(f(4   ,:,:) - f(mx-1,:,:)))
      ddxup(2   ,:,:) = (
     &               a*(f(3   ,:,:) - f(2   ,:,:)) +
     &               b*(f(4   ,:,:) - f(1   ,:,:)) +
     &               c*(f(5   ,:,:) - f(mx  ,:,:)))

*
      end
************************************************************************
      function ddyup (f)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddyup
!hpf$ align with TheCube:: f, ddyup

      c = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
      b = (-1.-120.*c)/24.
      a = (1.-3.*b-5.*c)/dy
      b = b/dy
      c = c/dy

      ddyup(:,3:my-3,:) = (
     &               a*(f(:,4:my-2,:) - f(:,3:my-3,:)) +
     &               b*(f(:,5:my-1,:) - f(:,2:my-4,:)) +
     &               c*(f(:,6:my  ,:) - f(:,1:my-5,:)))
      ddyup(:,my-2,:) = (
     &               a*(f(:,my-1,:) - f(:,my-2,:)) +
     &               b*(f(:,my  ,:) - f(:,my-3,:)) +
     &               c*(f(:,1   ,:) - f(:,my-4,:)))
      ddyup(:,my-1,:) = (
     &               a*(f(:,my  ,:) - f(:,my-1,:)) +
     &               b*(f(:,1   ,:) - f(:,my-2,:)) +
     &               c*(f(:,2   ,:) - f(:,my-3,:)))
      ddyup(:,my  ,:) = (
     &               a*(f(:,1   ,:) - f(:,my  ,:)) +
     &               b*(f(:,2   ,:) - f(:,my-1,:)) +
     &               c*(f(:,3   ,:) - f(:,my-2,:)))
      ddyup(:,1  ,:) = (
     &               a*(f(:,2   ,:) - f(:,1   ,:)) +
     &               b*(f(:,3   ,:) - f(:,my  ,:)) +
     &               c*(f(:,4   ,:) - f(:,my-1,:)))
      ddyup(:,2  ,:) = (
     &               a*(f(:,3   ,:) - f(:,2   ,:)) +
     &               b*(f(:,4   ,:) - f(:,1   ,:)) +
     &               c*(f(:,5   ,:) - f(:,my  ,:)))

*
      end
************************************************************************
      function ddzup (f)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddzup
!hpf$ align with TheCube:: f, ddzup

      c = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
      b = (-1.-120.*c)/24.
      a = (1.-3.*b-5.*c)/dz
      b = b/dz
      c = c/dz

      ddzup = (a*(cshift(f,dim=3,shift=1) - f)
     &   + b*(cshift(f,dim=3,shift=2) - cshift(f,dim=3,shift=-1))
     &   + c*(cshift(f,dim=3,shift=3) - cshift(f,dim=3,shift=-2))
     &     )
*
      end
************************************************************************
      function ddxdn (f)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddxdn
!hpf$ align with TheCube:: f, ddxdn

      c = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
      b = (-1.-120.*c)/24.
      a = (1.-3.*b-5.*c)/dx
      b = b/dx
      c = c/dx

      ddxdn(4:mx-2,:,:) = (
     &               a*(f(4:mx-2,:,:) - f(3:mx-3,:,:)) +
     &               b*(f(5:mx-1,:,:) - f(2:mx-4,:,:)) +
     &               c*(f(6:mx  ,:,:) - f(1:mx-5,:,:)))
      ddxdn(mx-1,:,:) = (
     &               a*(f(mx-1,:,:) - f(mx-2,:,:)) +
     &               b*(f(mx  ,:,:) - f(mx-3,:,:)) +
     &               c*(f(1   ,:,:) - f(mx-4,:,:)))
      ddxdn(mx  ,:,:) = (
     &               a*(f(mx  ,:,:) - f(mx-1,:,:)) +
     &               b*(f(1   ,:,:) - f(mx-2,:,:)) +
     &               c*(f(2   ,:,:) - f(mx-3,:,:)))
      ddxdn(1   ,:,:) = (
     &               a*(f(1   ,:,:) - f(mx  ,:,:)) +
     &               b*(f(2   ,:,:) - f(mx-1,:,:)) +
     &               c*(f(3   ,:,:) - f(mx-2,:,:)))
      ddxdn(2   ,:,:) = (
     &               a*(f(2   ,:,:) - f(1   ,:,:)) +
     &               b*(f(3   ,:,:) - f(mx  ,:,:)) +
     &               c*(f(4   ,:,:) - f(mx-1,:,:)))
      ddxdn(3   ,:,:) = (
     &               a*(f(3   ,:,:) - f(2   ,:,:)) +
     &               b*(f(4   ,:,:) - f(1   ,:,:)) +
     &               c*(f(5   ,:,:) - f(mx  ,:,:)))

*
      end
************************************************************************
      function ddydn (f)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddydn
!hpf$ align with TheCube:: f, ddydn

      c = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
      b = (-1.-120.*c)/24.
      a = (1.-3.*b-5.*c)/dy
      b = b/dy
      c = c/dy

      ddydn(:,4:my-2,:) = (
     &               a*(f(:,4:my-2,:) - f(:,3:my-3,:)) +
     &               b*(f(:,5:my-1,:) - f(:,2:my-4,:)) +
     &               c*(f(:,6:my  ,:) - f(:,1:my-5,:)))
      ddydn(:,my-1,:) = (
     &               a*(f(:,my-1,:) - f(:,my-2,:)) +
     &               b*(f(:,my  ,:) - f(:,my-3,:)) +
     &               c*(f(:,1   ,:) - f(:,my-4,:)))
      ddydn(:,my  ,:) = (
     &               a*(f(:,my  ,:) - f(:,my-1,:)) +
     &               b*(f(:,1   ,:) - f(:,my-2,:)) +
     &               c*(f(:,2   ,:) - f(:,my-3,:)))
      ddydn(:,1  ,:) = (
     &               a*(f(:,1   ,:) - f(:,my  ,:)) +
     &               b*(f(:,2   ,:) - f(:,my-1,:)) +
     &               c*(f(:,3   ,:) - f(:,my-2,:)))
      ddydn(:,2  ,:) = (
     &               a*(f(:,2   ,:) - f(:,1   ,:)) +
     &               b*(f(:,3   ,:) - f(:,my  ,:)) +
     &               c*(f(:,4   ,:) - f(:,my-1,:)))
      ddydn(:,3  ,:) = (
     &               a*(f(:,3   ,:) - f(:,2   ,:)) +
     &               b*(f(:,4   ,:) - f(:,1   ,:)) +
     &               c*(f(:,5   ,:) - f(:,my  ,:)))

*
      end
************************************************************************
      function ddzdn (f)
*
*  f is centered on (i-.5,j,k), fifth order stagger
*
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddzdn
!hpf$ align with TheCube:: f, ddzdn
*
*-----------------------------------------------------------------------
      c = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
      b = (-1.-120.*c)/24.
      a = (1.-3.*b-5.*c)/dz
      b = b/dz
      c = c/dz
*
      ddzdn = (a*(f - cshift(f,dim=3,shift=-1))
     &   + b*(cshift(f,dim=3,shift=1) - cshift(f,dim=3,shift=-2))
     &   + c*(cshift(f,dim=3,shift=2) - cshift(f,dim=3,shift=-3))
     &     )
*
      end
************************************************************************
      function xup1 (f)
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, xup1
!hpf$ align with TheCube:: f, xup1

      a = .5

      xup1(3:mx-3,:,:) = (
     &                 a*(f(3:mx-3,:,:) + f(4:mx-2,:,:)))
      xup1(mx-2,:,:) = (
     &               a*(f(mx-2,:,:) + f(mx-1,:,:)))
      xup1(mx-1,:,:) = (
     &               a*(f(mx-1,:,:) + f(mx,:,:)))
      xup1(mx ,:,:) = (
     &               a*(f(mx  ,:,:) + f(1,:,:)))
      xup1(1  ,:,:) = (
     &               a*(f(1   ,:,:) + f(2,:,:)))
      xup1(2  ,:,:) = (
     &               a*(f(2   ,:,:) + f(3,:,:)))

*
      end
************************************************************************
      function yup1 (f)
*
*  f is centered on (i-.5,j,k), fifth order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, yup1
!hpf$ align with TheCube:: f, yup1
*
*-----------------------------------------------------------------------
      a = .5
*

      yup1(:,3:my-3,:) = (
     &                 a*(f(:,3:my-3,:) + f(:,4:my-2,:)))
      yup1(:,2,:) = (
     &                 a*(f(:,2 ,:) + f(:,3,:)))
      yup1(:,1,:) = (
     &                 a*(f(:,1 ,:) + f(:,2,:)))
      yup1(:,my,:) = (
     &                 a*(f(:,my  ,:) + f(:,1,:)))
      yup1(:,my-1,:) = (
     &                 a*(f(:,my-1,:) + f(:,my,:)))
      yup1(:,my-2,:) = (
     &                 a*(f(:,my-2,:) + f(:,my-1,:)))
*
      end
************************************************************************
      function zup1 (f)
*
*  f is centered on (i-.5,j,k), fifth order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, zup1
!hpf$ align with TheCube:: f, zup1
*
*-----------------------------------------------------------------------
      a = .5
*
      zup1 = (a*(cshift(f,dim=3,shift=+1) + f)
     &     )

      end
************************************************************************
      function xdn1 (f)
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, xdn1
!hpf$ align with TheCube:: f, xdn1

      a = .5

      xdn1(4:mx-2  ,:,:) = (
     &                 a*(f(3:mx-3,:,:) + f(4:mx-2,:,:)))
      xdn1(mx-1,:,:) = (
     &               a*(f(mx-2,:,:) + f(mx-1,:,:)))
      xdn1(mx  ,:,:) = (
     &               a*(f(mx-1,:,:) + f(mx,:,:)))
      xdn1(1   ,:,:) = (
     &               a*(f(mx  ,:,:) + f(1,:,:)))
      xdn1(2  ,:,:) = (
     &               a*(f(1   ,:,:) + f(2,:,:)))
      xdn1(3  ,:,:) = (
     &               a*(f(2   ,:,:) + f(3,:,:)))

*
      end
************************************************************************
      function ydn1 (f)
*
*  f is centered on (i-.5,j,k), fifth order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, ydn1
!hpf$ align with TheCube:: f, ydn1
*
*-----------------------------------------------------------------------
      a = .5
*

      ydn1(:,4:my-2,:) = (
     &                 a*(f(:,3:my-3,:) + f(:,4:my-2,:)))
      ydn1(:,3,:) = (
     &                 a*(f(:,2 ,:) + f(:,3,:)))
      ydn1(:,2,:) = (
     &                 a*(f(:,1 ,:) + f(:,2,:)))
      ydn1(:,1,:) = (
     &                 a*(f(:,my  ,:) + f(:,1,:)))
      ydn1(:,my,:) = (
     &                 a*(f(:,my-1,:) + f(:,my,:)))
      ydn1(:,my-1,:) = (
     &                 a*(f(:,my-2,:) + f(:,my-1,:)))
*
      end
************************************************************************
      function zdn1 (f)
*
*  f is centered on (i-.5,j,k), fifth order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, zdn1
!hpf$ align with TheCube:: f, zdn1
*
*-----------------------------------------------------------------------
      a = .5
*
      zdn1 = (a*(cshift(f,dim=3,shift=-1) + f)
     &     )
*
      end
************************************************************************
      function ddxup1 (f)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddxup1
!hpf$ align with TheCube:: f, ddxup1

      adx1=1./dx
      ddxup1(3:mx-3,:,:) = (
     &               adx1*(f(4:mx-2,:,:) - f(3:mx-3,:,:)))
      ddxup1(mx-2,:,:) = (
     &               adx1*(f(mx-1,:,:) - f(mx-2,:,:)))
      ddxup1(mx-1 ,:,:) = (
     &               adx1*(f(mx  ,:,:) - f(mx-1,:,:)))
      ddxup1(mx  ,:,:) = (
     &               adx1*(f(1   ,:,:) - f(mx  ,:,:)))
      ddxup1(1   ,:,:) = (
     &               adx1*(f(2   ,:,:) - f(1   ,:,:)))
      ddxup1(2   ,:,:) = (
     &               adx1*(f(3   ,:,:) - f(2   ,:,:)))

*
      end
************************************************************************
      function ddyup1 (f)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddyup1
!hpf$ align with TheCube:: f, ddyup1

      ady1=1./dy
      ddyup1(:,3:my-3,:) = (
     &               ady1*(f(:,4:my-2,:) - f(:,3:my-3,:)))
      ddyup1(:,my-2,:) = (
     &               ady1*(f(:,my-1,:) - f(:,my-2,:)))
      ddyup1(:,my-1,:) = (
     &               ady1*(f(:,my  ,:) - f(:,my-1,:)))
      ddyup1(:,my  ,:) = (
     &               ady1*(f(:,1   ,:) - f(:,my  ,:)))
      ddyup1(:,1  ,:) = (
     &               ady1*(f(:,2   ,:) - f(:,1   ,:)))
      ddyup1(:,2  ,:) = (
     &               ady1*(f(:,3   ,:) - f(:,2   ,:)))

*
      end
************************************************************************
      function ddzup1 (f)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddzup1
!hpf$ align with TheCube:: f, ddzup1

      adz1=1./dz
      ddzup1 = (adz1*(cshift(f,dim=3,shift=1) - f)
     &     )
*
      end
************************************************************************
      function ddxdn1 (f)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddxdn1
!hpf$ align with TheCube:: f, ddxdn1

      adx1=1./dx
      ddxdn1(4:mx-2,:,:) = (
     &               adx1*(f(4:mx-2,:,:) - f(3:mx-3,:,:)))
      ddxdn1(mx-1,:,:) = (
     &               adx1*(f(mx-1,:,:) - f(mx-2,:,:)))
      ddxdn1(mx  ,:,:) = (
     &               adx1*(f(mx  ,:,:) - f(mx-1,:,:)))
      ddxdn1(1   ,:,:) = (
     &               adx1*(f(1   ,:,:) - f(mx  ,:,:)))
      ddxdn1(2   ,:,:) = (
     &               adx1*(f(2   ,:,:) - f(1   ,:,:)))
      ddxdn1(3   ,:,:) = (
     &               adx1*(f(3   ,:,:) - f(2   ,:,:)))

*
      end
************************************************************************
      function ddydn1 (f)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddydn1
!hpf$ align with TheCube:: f, ddydn1

      ady1=1./dy
      ddydn1(:,4:my-2,:) = (
     &               ady1*(f(:,4:my-2,:) - f(:,3:my-3,:)))
      ddydn1(:,my-1,:) = (
     &               ady1*(f(:,my-1,:) - f(:,my-2,:)))
      ddydn1(:,my  ,:) = (
     &               ady1*(f(:,my  ,:) - f(:,my-1,:)))
      ddydn1(:,1  ,:) = (
     &               ady1*(f(:,1   ,:) - f(:,my  ,:)))
      ddydn1(:,2  ,:) = (
     &               ady1*(f(:,2   ,:) - f(:,1   ,:)))
      ddydn1(:,3  ,:) = (
     &               ady1*(f(:,3   ,:) - f(:,2   ,:)))

*
      end
************************************************************************
      function ddzdn1 (f)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddzdn1
!hpf$ align with TheCube:: f, ddzdn1

      adz1=1./dz
      ddzdn1 = (adz1*(f - cshift(f,dim=3,shift=-1))
     &     )
*
      end
*-----------------------------------------------------------------------
      function difx1(f)
*
*  Finite difference one zone left to two zones right.
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, difx1
!hpf$ align with TheCube:: f, difx1
*-----------------------------------------------------------------------
      difx1 = f - cshift(f,dim=1,shift=1)
      nflop=nflop+1
      END
*-----------------------------------------------------------------------
      function dify1(f)
*
*  Finite difference one zone left to two zones right.
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, dify1
!hpf$ align with TheCube:: f, dify1
*-----------------------------------------------------------------------
      if (my .eq. 1) then
	dify1 = 0.0
	return
      endif
      dify1 = f - cshift(f,dim=2,shift=1)
      nflop=nflop+1
      END
*-----------------------------------------------------------------------
      function difz1(f)
*
*  Finite difference one zone left to two zones right.
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, difz1
!hpf$ align with TheCube:: f, difz1
*-----------------------------------------------------------------------
      if (mz .eq. 1) then
	difz1 = 0.0
	return
      endif
      difz1 = f - cshift(f,dim=3,shift=1)
      nflop=nflop+1
      END
************************************************************************
      function max3 (f)
      include 'cparam.inc'
      include 'cdata.inc'
*
      real, dimension (mx,my,mz) :: f, max3
!hpf$ align with TheCube:: f, max3

      max3 = amax1(f,
     &      cshift(f,dim=1,shift=+1),
     &      cshift(f,dim=1,shift=-1))
      if (my .gt. 1) then
        max3 = amax1(max3,
     &        cshift(max3,dim=2,shift=+1),
     &        cshift(max3,dim=2,shift=-1))
      endif
      max3 = amax1(max3,
     &      cshift(max3,dim=3,shift=+1),
     &      cshift(max3,dim=3,shift=-1))

      nflop=nflop+9
      end
*----------------------------------------------------------------------
      function smooth3 (f)
      include 'cparam.inc'
      include 'cdata.inc'
*
      real, dimension (mx,my,mz) :: f, smooth3
!hpf$ align with TheCube:: f, smooth3

      smooth3=(1./3.)*(f + 
     & cshift(f,dim=1,shift=+1) + cshift(f,dim=1,shift=-1))
      if (my .gt. 1) then
        smooth3=(1./3.)*(smooth3 + 
     & cshift(smooth3,dim=2,shift=+1) + cshift(smooth3,dim=2,shift=-1))
      endif
      smooth3=(1./3.)*(smooth3 + 
     & cshift(smooth3,dim=3,shift=+1) + cshift(smooth3,dim=3,shift=-1))

      nflop=nflop+9
      end

************************************************************************
      subroutine init_derivs ()
*
*  Initialize the derivatives setup
*
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cstagger.inc'
*-----------------------------------------------------------------------
      c = 3./256.
      b = -25./256.
      a = .5-b-c
      cdx = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
      bdx = (-1.-120.*cdx)/24.
      adx = (1.-3.*bdx-5.*cdx)/dx
      bdx = bdx/dx
      cdx = cdx/dx
      cdy = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
      bdy = (-1.-120.*cdy)/24.
      ady = (1.-3.*bdy-5.*cdy)/dy
      bdy = bdy/dy
      cdy = cdy/dy
      cdz = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
      bdz = (-1.-120.*cdz)/24.
      adz = (1.-3.*bdz-5.*cdz)/dz
      bdz = bdz/dz
      cdz = cdz/dz

      adx1 = 1./dx
      ady1 = 1./dy
      adz1 = 1./dz

      end
*-----------------------------------------------------------------------
      subroutine init_stagger ()
*
*  Initialize the derivatives setup
*
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cstagger.inc'
*-----------------------------------------------------------------------
      c = 3./256.
      b = -25./256.
      a = .5-b-c

      end
