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

      if (mx .eq. 1) then
        xup=f
        return
      endif

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

C$ do across local(i,j,k)
      do k=1,mz
       do j=1,my
        xup(mx-2,j,k) = (
     &                 a*(f(mx-2,j,k)+f(mx-1,j,k))
     &               + b*(f(mx-3,j,k)+f(mx  ,j,k))
     &               + c*(f(mx-4,j,k)+f(1   ,j,k)))
        xup(mx-1,j,k) = (
     &                 a*(f(mx-1,j,k)+f(mx  ,j,k))
     &               + b*(f(mx-2,j,k)+f(1   ,j,k))
     &               + c*(f(mx-3,j,k)+f(2   ,j,k)))
        xup(mx  ,j,k) = (
     &                 a*(f(mx  ,j,k)+f(1   ,j,k))
     &               + b*(f(mx-1,j,k)+f(2   ,j,k))
     &               + c*(f(mx-2,j,k)+f(3   ,j,k)))
        xup(1   ,j,k) = (
     &                 a*(f(1   ,j,k)+f(2   ,j,k))
     &               + b*(f(mx  ,j,k)+f(3   ,j,k))
     &               + c*(f(mx-1,j,k)+f(4   ,j,k)))
        xup(2   ,j,k) = (
     &                 a*(f(2   ,j,k)+f(3   ,j,k))
     &               + b*(f(1   ,j,k)+f(4   ,j,k))
     &               + c*(f(mx  ,j,k)+f(5   ,j,k)))
       end do
       do j=1,my
        do i=3,mx-3
          xup(i  ,j,k) = (
     &                  a*(f(i  ,j,k)+f(i+1 ,j,k))
     &                + b*(f(i-1,j,k)+f(i+2 ,j,k))
     &                + c*(f(i-2,j,k)+f(i+3 ,j,k)))
        end do
       end do
      end do

*
      nflop=nflop+8
      end
************************************************************************
      function yup (f)
*
*  f is centered on (i,j-.5,k)
*  yup(f) is centered on (i,j,k), fifth order stagger
*
*-----------------------------------------------------------------------
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, yup
!hpf$ align with TheCube:: f, yup
*
      if (my .eq. 1) then
        yup=f
        return
      endif

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

      if (my .eq. 1) then
	do k=1,mz
	do i=1,mx
	  yup(i,1,k)=f(i,1,k)
	enddo
	enddo
	return
      endif

C$ do across local(i,j,k)
      do k=1,mz
        do i=1,mx
          yup(i,my-2,k) = (
     &                   a*(f(i,my-2,k)+f(i,my-1,k)) +
     &                   b*(f(i,my-3,k)+f(i,my  ,k)) + 
     &                   c*(f(i,my-4,k)+f(i,1   ,k)))
          yup(i,my-1,k) = (
     &                   a*(f(i,my-1,k)+f(i,my  ,k)) +
     &                   b*(f(i,my-2,k)+f(i,1   ,k)) + 
     &                   c*(f(i,my-3,k)+f(i,2   ,k)))
          yup(i,my  ,k) = (
     &                   a*(f(i,my  ,k)+f(i,1   ,k)) +
     &                   b*(f(i,my-1,k)+f(i,2   ,k)) + 
     &                   c*(f(i,my-2,k)+f(i,3   ,k)))
          yup(i,1   ,k) = (
     &                   a*(f(i,1   ,k)+f(i,2   ,k)) +
     &                   b*(f(i,my  ,k)+f(i,3   ,k)) + 
     &                   c*(f(i,my-1,k)+f(i,4   ,k)))
          yup(i,2   ,k) = (
     &                   a*(f(i,2   ,k)+f(i,3   ,k)) +
     &                   b*(f(i,1   ,k)+f(i,4   ,k)) + 
     &                   c*(f(i,my  ,k)+f(i,5   ,k)))
        end do
        do j=3,my-3
         do i=1,mx
          yup(i,j  ,k) = (
     &                   a*(f(i,j   ,k)+f(i,j+1 ,k)) +
     &                   b*(f(i,j-1 ,k)+f(i,j+2 ,k)) + 
     &                   c*(f(i,j-2 ,k)+f(i,j+3 ,k)))
         end do
        end do
      end do
*
      nflop=nflop+8
      end
************************************************************************
      function zup (f)
*
*  f is centered on (i,j,k-.5)
*  zup(f) is centered on (i,j,k), fifth order stagger
*
*-----------------------------------------------------------------------
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, zup
!hpf$ align with TheCube:: f, zup
*
      if (mz .eq. 1) then
        zup=f
        return
      endif

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

C*$* ASSERT DO (SERIAL)
      do j=1,my
        do i=1,mx
          zup(i,j,mz-2) = (
     &                   a*(f(i,j,mz-2)+f(i,j,mz-1)) +
     &                   b*(f(i,j,mz-3)+f(i,j,mz  )) +
     &                   c*(f(i,j,mz-4)+f(i,j,1   )))
          zup(i,j,mz-1) = (
     &                   a*(f(i,j,mz-1)+f(i,j,mz  )) +
     &                   b*(f(i,j,mz-2)+f(i,j,1   )) +
     &                   c*(f(i,j,mz-3)+f(i,j,2   )))
          zup(i,j,mz  ) = (
     &                   a*(f(i,j,mz  )+f(i,j,1   )) +
     &                   b*(f(i,j,mz-1)+f(i,j,2   )) +
     &                   c*(f(i,j,mz-2)+f(i,j,3   )))
          zup(i,j,1   ) = (
     &                   a*(f(i,j,1   )+f(i,j,2   )) +
     &                   b*(f(i,j,mz  )+f(i,j,3   )) +
     &                   c*(f(i,j,mz-1)+f(i,j,4   )))
          zup(i,j,2   ) = (
     &                   a*(f(i,j,2   )+f(i,j,3   )) +
     &                   b*(f(i,j,1   )+f(i,j,4   )) +
     &                   c*(f(i,j,mz  )+f(i,j,5   )))
        end do
      end do

C$ do across local(i,j,k)
      do k=3,mz-3
        do j=1,my
          do i=1,mx
            zup(i,j,k) = (
     &                   a*(f(i,j,k   )+f(i,j,k+1 )) +
     &                   b*(f(i,j,k-1 )+f(i,j,k+2 )) + 
     &                   c*(f(i,j,k-2 )+f(i,j,k+3 )))
          end do
        end do
      end do
      nflop=nflop+8
      end
************************************************************************
      function xdn (f)
*
*  f is centered on (i,j,k)
*  xdn(f) is centered on (i-.5,j,k), fifth order stagger
*
*-----------------------------------------------------------------------
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, xdn
!hpf$ align with TheCube:: f, xdn

      if (mx .eq. 1) then
        xdn=f
        return
      endif

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

C$ do across local(i,j,k)
      do k=1,mz
      do j=1,my
        xdn(mx-1,j,k) = (
     &                 a*(f(mx-2,j,k)+f(mx-1,j,k))
     &               + b*(f(mx-3,j,k)+f(mx  ,j,k))
     &               + c*(f(mx-4,j,k)+f(1   ,j,k)))
        xdn(mx  ,j,k) = (
     &                 a*(f(mx-1,j,k)+f(mx  ,j,k))
     &               + b*(f(mx-2,j,k)+f(1   ,j,k))
     &               + c*(f(mx-3,j,k)+f(2   ,j,k)))
        xdn(1  ,j,k) = (
     &                 a*(f(mx  ,j,k)+f(1   ,j,k))
     &               + b*(f(mx-1,j,k)+f(2   ,j,k))
     &               + c*(f(mx-2,j,k)+f(3   ,j,k)))
        xdn(2  ,j,k) = (
     &                 a*(f(1   ,j,k)+f(2   ,j,k))
     &               + b*(f(mx  ,j,k)+f(3   ,j,k))
     &               + c*(f(mx-1,j,k)+f(4   ,j,k)))
        xdn(3  ,j,k) = (
     &                 a*(f(2   ,j,k)+f(3   ,j,k))
     &               + b*(f(1   ,j,k)+f(4   ,j,k))
     &               + c*(f(mx  ,j,k)+f(5   ,j,k)))
      end do
      do j=1,my
        do i=3,mx-3
          xdn(i+1,j,k) = (
     &                 a*(f(i   ,j,k)+f(i+1 ,j,k))
     &               + b*(f(i-1 ,j,k)+f(i+2 ,j,k))
     &               + c*(f(i-2 ,j,k)+f(i+3 ,j,k)))
        end do
      end do
      do j=1,my
      end do
      end do

*
      nflop=nflop+8
      end
************************************************************************
      function ydn (f)
*
*  f is centered on (i,j,k)
*  ydn(f) is centered on (i,j-.5,k), fifth order stagger
*
*-----------------------------------------------------------------------
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, ydn
!hpf$ align with TheCube:: f, ydn
*
      if (my .eq. 1) then
        ydn=f
        return
      endif

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

      if (my .eq. 1) then
	do k=1,mz
	do i=1,mx
	  ydn(i,1,k)=f(i,1,k)
	enddo
	enddo
	return
      endif

*
C$ do across local(i,j,k)
      do k=1,mz
        do i=1,mx
          ydn(i,my-1,k) = (
     &                   a*(f(i,my-2,k)+f(i,my-1,k)) +
     &                   b*(f(i,my-3,k)+f(i,my  ,k)) + 
     &                   c*(f(i,my-4,k)+f(i,1   ,k)))
          ydn(i,my  ,k) = (
     &                   a*(f(i,my-1,k)+f(i,my  ,k)) +
     &                   b*(f(i,my-2,k)+f(i,1   ,k)) + 
     &                   c*(f(i,my-3,k)+f(i,2   ,k)))
          ydn(i,1   ,k) = (
     &                   a*(f(i,my  ,k)+f(i,1   ,k)) +
     &                   b*(f(i,my-1,k)+f(i,2   ,k)) + 
     &                   c*(f(i,my-2,k)+f(i,3   ,k)))
          ydn(i,2   ,k) = (
     &                   a*(f(i,1   ,k)+f(i,2   ,k)) +
     &                   b*(f(i,my  ,k)+f(i,3   ,k)) + 
     &                   c*(f(i,my-1,k)+f(i,4   ,k)))
          ydn(i,3   ,k) = (
     &                   a*(f(i,2   ,k)+f(i,3   ,k)) +
     &                   b*(f(i,1   ,k)+f(i,4   ,k)) + 
     &                   c*(f(i,my  ,k)+f(i,5   ,k)))
        end do
        do j=3,my-3
         do i=1,mx
          ydn(i,j+1,k) = (
     &                   a*(f(i,j   ,k)+f(i,j+1 ,k)) +
     &                   b*(f(i,j-1 ,k)+f(i,j+2 ,k)) + 
     &                   c*(f(i,j-2 ,k)+f(i,j+3 ,k)))
         end do
        end do
      end do
*
      nflop=nflop+8
      end
************************************************************************
      function zdn (f)
*
*  f is centered on (i,j,k)
*  zdn(f) is centered on (i,j,k-.5), fifth order stagger
*
*-----------------------------------------------------------------------
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, zdn
!hpf$ align with TheCube:: f, zdn
*
      if (mz .eq. 1) then
        zdn=f
        return
      endif

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

C*$* ASSERT DO (SERIAL)
      do j=1,my
        do i=1,mx
          zdn(i,j,mz-1) = (
     &                   a*(f(i,j,mz-2)+f(i,j,mz-1)) +
     &                   b*(f(i,j,mz-3)+f(i,j,mz  )) +
     &                   c*(f(i,j,mz-4)+f(i,j,1   )))
          zdn(i,j,mz  ) = (
     &                   a*(f(i,j,mz-1)+f(i,j,mz  )) +
     &                   b*(f(i,j,mz-2)+f(i,j,1   )) +
     &                   c*(f(i,j,mz-3)+f(i,j,2   )))
          zdn(i,j,1   ) = (
     &                   a*(f(i,j,mz  )+f(i,j,1   )) +
     &                   b*(f(i,j,mz-1)+f(i,j,2   )) +
     &                   c*(f(i,j,mz-2)+f(i,j,3   )))
          zdn(i,j,2   ) = (
     &                   a*(f(i,j,1   )+f(i,j,2   )) +
     &                   b*(f(i,j,mz  )+f(i,j,3   )) +
     &                   c*(f(i,j,mz-1)+f(i,j,4   )))
          zdn(i,j,3   ) = (
     &                   a*(f(i,j,2   )+f(i,j,3   )) +
     &                   b*(f(i,j,1   )+f(i,j,4   )) +
     &                   c*(f(i,j,mz  )+f(i,j,5   )))
        end do
      end do

C$ do across local(i,j,k)
      do k=3,mz-3
        do j=1,my
          do i=1,mx
            zdn(i,j,k+1) = (
     &                   a*(f(i,j,k   )+f(i,j,k+1 )) +
     &                   b*(f(i,j,k-1 )+f(i,j,k+2 )) + 
     &                   c*(f(i,j,k-2 )+f(i,j,k+3 )))
          end do
        end do
      end do
*
      nflop=nflop+8
      end
************************************************************************
      function ddxup (f)
*
*  f is centered on (i-.5,j,k)
*  ddxup(f) is centered on (i,j,k), sixth order stagger
*
*-----------------------------------------------------------------------
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddxup
!hpf$ align with TheCube:: f, ddxup

      if (mx .eq. 1) then
        ddxup=0.
        return
      endif

      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

C$ do across local(i,j,k)
      do k=1,mz
      do j=1,my
        ddxup(mx-2,j,k) = (
     &                 a*(f(mx-1,j,k)-f(mx-2,j,k))
     &               + b*(f(mx  ,j,k)-f(mx-3,j,k))
     &               + c*(f(1   ,j,k)-f(mx-4,j,k)))
        ddxup(mx-1,j,k) = (
     &                 a*(f(mx  ,j,k)-f(mx-1,j,k))
     &               + b*(f(1   ,j,k)-f(mx-2,j,k))
     &               + c*(f(2   ,j,k)-f(mx-3,j,k)))
        ddxup(mx  ,j,k) = (
     &                 a*(f(1   ,j,k)-f(mx  ,j,k))
     &               + b*(f(2   ,j,k)-f(mx-1,j,k))
     &               + c*(f(3   ,j,k)-f(mx-2,j,k)))
        ddxup(1  ,j,k) = (
     &                 a*(f(2   ,j,k)-f(1   ,j,k))
     &               + b*(f(3   ,j,k)-f(mx  ,j,k))
     &               + c*(f(4   ,j,k)-f(mx-1,j,k)))
        ddxup(2  ,j,k) = (
     &                 a*(f(3   ,j,k)-f(2  ,j,k))
     &               + b*(f(4   ,j,k)-f(1  ,j,k))
     &               + c*(f(5   ,j,k)-f(mx ,j,k)))
      end do
      do j=1,my
        do i=3,mx-3
          ddxup(i  ,j,k) = (
     &                 a*(f(i+1 ,j,k)-f(i   ,j,k))
     &               + b*(f(i+2 ,j,k)-f(i-1 ,j,k))
     &               + c*(f(i+3 ,j,k)-f(i-2 ,j,k)))
        end do
      end do
      end do

*
      nflop=nflop+8
      end
************************************************************************
      function ddyup (f)
*
*  f is centered on (i,j-.5,k)
*  ddyup(f) is centered on (i,j,k), sixth order stagger
*
*-----------------------------------------------------------------------
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddyup
!hpf$ align with TheCube:: f, ddyup

      if (my .eq. 1) then
        ddyup=0.
        return
      endif

      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

C$ do across local(i,j,k)
      do k=1,mz
        do i=1,mx
          ddyup(i,my-2,k) = (
     &                   a*(f(i,my-1,k)-f(i,my-2,k)) +
     &                   b*(f(i,my  ,k)-f(i,my-3,k)) + 
     &                   c*(f(i,1   ,k)-f(i,my-4,k)))
          ddyup(i,my-1,k) = (
     &                   a*(f(i,my  ,k)-f(i,my-1,k)) +
     &                   b*(f(i,1   ,k)-f(i,my-2,k)) + 
     &                   c*(f(i,2   ,k)-f(i,my-3,k)))
          ddyup(i,my  ,k) = (
     &                   a*(f(i,1   ,k)-f(i,my  ,k)) +
     &                   b*(f(i,2   ,k)-f(i,my-1,k)) + 
     &                   c*(f(i,3   ,k)-f(i,my-2,k)))
          ddyup(i,1   ,k) = (
     &                   a*(f(i,2   ,k)-f(i,1   ,k)) +
     &                   b*(f(i,3   ,k)-f(i,my  ,k)) + 
     &                   c*(f(i,4   ,k)-f(i,my-1,k)))
          ddyup(i,2   ,k) = (
     &                   a*(f(i,3   ,k)-f(i,2   ,k)) +
     &                   b*(f(i,4   ,k)-f(i,1   ,k)) + 
     &                   c*(f(i,5   ,k)-f(i,my  ,k)))
        end do
        do j=3,my-3
         do i=1,mx
          ddyup(i,j  ,k) = (
     &                   a*(f(i,j+1 ,k)-f(i,j   ,k)) +
     &                   b*(f(i,j+2 ,k)-f(i,j-1 ,k)) + 
     &                   c*(f(i,j+3 ,k)-f(i,j-2 ,k)))
         end do
        end do
      end do

*
      nflop=nflop+8
      end
************************************************************************
      function ddzup (f)
*
*  f is centered on (i,j,k-.5)
*  ddzup(f) is centered on (i,j,k), sixth order stagger
*
*-----------------------------------------------------------------------
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddzup
!hpf$ align with TheCube:: f, ddzup

      if (mz .eq. 1) then
        ddzup=0.
        return
      endif

      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

C*$* ASSERT DO (SERIAL)
      do j=1,my
        do i=1,mx
          ddzup(i,j,mz-2) = (
     &                   a*(f(i,j,mz-1)-f(i,j,mz-2)) +
     &                   b*(f(i,j,mz  )-f(i,j,mz-3)) +
     &                   c*(f(i,j,1   )-f(i,j,mz-4)))
          ddzup(i,j,mz-1) = (
     &                   a*(f(i,j,mz  )-f(i,j,mz-1)) +
     &                   b*(f(i,j,1   )-f(i,j,mz-2)) +
     &                   c*(f(i,j,2   )-f(i,j,mz-3)))
          ddzup(i,j,mz  ) = (
     &                   a*(f(i,j,1   )-f(i,j,mz  )) +
     &                   b*(f(i,j,2   )-f(i,j,mz-1)) +
     &                   c*(f(i,j,3   )-f(i,j,mz-2)))
          ddzup(i,j,1   ) = (
     &                   a*(f(i,j,2   )-f(i,j,1   )) +
     &                   b*(f(i,j,3   )-f(i,j,mz  )) +
     &                   c*(f(i,j,4   )-f(i,j,mz-1)))
          ddzup(i,j,2   ) = (
     &                   a*(f(i,j,3   )-f(i,j,2   )) +
     &                   b*(f(i,j,4   )-f(i,j,1   )) +
     &                   c*(f(i,j,5  )-f(i,j,mz  )))
        end do
      end do

C$ do across local(i,j,k)
      do k=3,mz-3
        do j=1,my
          do i=1,mx
            ddzup(i,j,k) = (
     &                   a*(f(i,j,k+1 )-f(i,j,k   )) +
     &                   b*(f(i,j,k+2 )-f(i,j,k-1 )) + 
     &                   c*(f(i,j,k+3 )-f(i,j,k-2 )))
          end do
        end do
      end do
*
      nflop=nflop+8
      end
************************************************************************
      function ddxdn (f)
*
*  f is centered on (i,j,k)
*  ddxdn(f) is centered on (i-.5,j,k), sixth order stagger
*
*-----------------------------------------------------------------------
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddxdn
!hpf$ align with TheCube:: f, ddxdn

      if (mx .eq. 1) then
        ddxdn=0.
        return
      endif

      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

C$ do across local(i,j,k)
      do k=1,mz
       do j=1,my
        ddxdn(mx-1,j,k) = (
     &                 a*(f(mx-1,j,k)-f(mx-2,j,k))
     &               + b*(f(mx  ,j,k)-f(mx-3,j,k))
     &               + c*(f(1   ,j,k)-f(mx-4,j,k)))
        ddxdn(mx  ,j,k) = (
     &                 a*(f(mx  ,j,k)-f(mx-1,j,k))
     &               + b*(f(1   ,j,k)-f(mx-2,j,k))
     &               + c*(f(2   ,j,k)-f(mx-3,j,k)))
        ddxdn(1   ,j,k) = (
     &                 a*(f(1   ,j,k)-f(mx  ,j,k))
     &               + b*(f(2   ,j,k)-f(mx-1,j,k))
     &               + c*(f(3   ,j,k)-f(mx-2,j,k)))
        ddxdn(2  ,j,k) = (
     &                 a*(f(2   ,j,k)-f(1   ,j,k))
     &               + b*(f(3   ,j,k)-f(mx  ,j,k))
     &               + c*(f(4   ,j,k)-f(mx-1,j,k)))
        ddxdn(3  ,j,k) = (
     &                 a*(f(3   ,j,k)-f(2   ,j,k))
     &               + b*(f(4   ,j,k)-f(1   ,j,k))
     &               + c*(f(5   ,j,k)-f(mx  ,j,k)))
       end do
       do j=1,my
         do i=3,mx-3
          ddxdn(i+1,j,k) = (
     &                 a*(f(i+1 ,j,k)-f(i   ,j,k))
     &               + b*(f(i+2 ,j,k)-f(i-1 ,j,k))
     &               + c*(f(i+3 ,j,k)-f(i-2 ,j,k)))
        end do
       end do
      end do
*
      nflop=nflop+8
      end
************************************************************************
      function ddydn (f)
*
*  f is centered on (i,j,k)
*  ddydn(f) is centered on (i,j-.5,k), sixth order stagger
*
*-----------------------------------------------------------------------
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddydn
!hpf$ align with TheCube:: f, ddydn

      if (my .eq. 1) then
        ddydn=0.
        return
      endif

      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

C$ do across local(i,j,k)
      do k=1,mz
        do i=1,mx
          ddydn(i,my-1,k) = (
     &                   a*(f(i,my-1,k)-f(i,my-2,k)) +
     &                   b*(f(i,my  ,k)-f(i,my-3,k)) + 
     &                   c*(f(i,1   ,k)-f(i,my-4,k)))
          ddydn(i,my  ,k) = (
     &                   a*(f(i,my  ,k)-f(i,my-1,k)) +
     &                   b*(f(i,1   ,k)-f(i,my-2,k)) + 
     &                   c*(f(i,2   ,k)-f(i,my-3,k)))
          ddydn(i,1  ,k) = (
     &                   a*(f(i,1   ,k)-f(i,my  ,k)) +
     &                   b*(f(i,2   ,k)-f(i,my-1,k)) + 
     &                   c*(f(i,3   ,k)-f(i,my-2,k)))
          ddydn(i,2   ,k) = (
     &                   a*(f(i,2   ,k)-f(i,1   ,k)) +
     &                   b*(f(i,3   ,k)-f(i,my  ,k)) + 
     &                   c*(f(i,4   ,k)-f(i,my-1,k)))
          ddydn(i,3   ,k) = (
     &                   a*(f(i,3   ,k)-f(i,2   ,k)) +
     &                   b*(f(i,4   ,k)-f(i,1   ,k)) + 
     &                   c*(f(i,5   ,k)-f(i,my  ,k)))
        end do
        do j=3,my-3
         do i=1,mx
          ddydn(i,j+1 ,k) = (
     &                   a*(f(i,j+1 ,k)-f(i,j   ,k)) +
     &                   b*(f(i,j+2 ,k)-f(i,j-1 ,k)) + 
     &                   c*(f(i,j+3 ,k)-f(i,j-2 ,k)))
         end do
        end do
      end do

*
      nflop=nflop+8
      end
************************************************************************
      function ddzdn (f)
*
*  f is centered on (i,j,k)
*  ddzdn(f) is centered on (i,j,k-.5), sixth order stagger
*
*-----------------------------------------------------------------------
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddzdn
!hpf$ align with TheCube:: f, ddzdn

      if (mz .eq. 1) then
        ddzdn=0.
        return
      endif

      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

C*$* ASSERT DO (SERIAL)
      do j=1,my
        do i=1,mx
          ddzdn(i,j,mz-1) = (
     &                   a*(f(i,j,mz-1)-f(i,j,mz-2)) +
     &                   b*(f(i,j,mz  )-f(i,j,mz-3)) +
     &                   c*(f(i,j,1   )-f(i,j,mz-4)))
          ddzdn(i,j,mz  ) = (
     &                   a*(f(i,j,mz  )-f(i,j,mz-1)) +
     &                   b*(f(i,j,1   )-f(i,j,mz-2)) +
     &                   c*(f(i,j,2   )-f(i,j,mz-3)))
          ddzdn(i,j,1   ) = (
     &                   a*(f(i,j,1   )-f(i,j,mz  )) +
     &                   b*(f(i,j,2   )-f(i,j,mz-1)) +
     &                   c*(f(i,j,3   )-f(i,j,mz-2)))
          ddzdn(i,j,2   ) = (
     &                   a*(f(i,j,2   )-f(i,j,1   )) +
     &                   b*(f(i,j,3   )-f(i,j,mz  )) +
     &                   c*(f(i,j,4   )-f(i,j,mz-1)))
          ddzdn(i,j,3   ) = (
     &                   a*(f(i,j,3   )-f(i,j,2   )) +
     &                   b*(f(i,j,4   )-f(i,j,1   )) +
     &                   c*(f(i,j,5   )-f(i,j,mz  )))
        end do
      end do
C$ do across local(i,j,k)
      do k=3,mz-3
        do j=1,my
          do i=1,mx
            ddzdn(i,j,k+1) = (
     &                   a*(f(i,j,k+1 )-f(i,j,k   )) +
     &                   b*(f(i,j,k+2 )-f(i,j,k-1 )) + 
     &                   c*(f(i,j,k+3 )-f(i,j,k-2 )))
          end do
        end do
      end do
*
      nflop=nflop+8
      end

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

      if (mx .eq. 1) then
        xup1=f
        return
      endif

      a = .5
*
C$ do across local(i,j,k)
      do k=1,mz
       do j=1,my
        xup1(mx  ,j,k) = (
     &                 a*(f(mx  ,j,k)+f(1   ,j,k)))
       end do
       do j=1,my
        do i=1,mx-1
          xup1(i  ,j,k) = (
     &                  a*(f(i  ,j,k)+f(i+1 ,j,k)))
        end do
       end do
      end do
*
      nflop=nflop+2
      end
************************************************************************
      function yup1 (f)
*
*  f is centered on (i,j-.5,k)
*  yup1(f) is centered on (i,j,k), first order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, yup1
!hpf$ align with TheCube:: f, yup1
*
*-----------------------------------------------------------------------
      if (my .eq. 1) then
        yup1=f
        return
      endif

      a = .5

      if (my .eq. 1) then
	do k=1,mz
	do i=1,mx
	  yup1(i,1,k)=f(i,1,k)
	enddo
	enddo
	return
      endif
*
C$ do across local(i,j,k)
      do k=1,mz
        do i=1,mx
          yup1(i,my  ,k) = (
     &                   a*(f(i,my  ,k)+f(i,1   ,k)))
        end do
        do j=1,my-1
         do i=1,mx
          yup1(i,j  ,k) = (
     &                   a*(f(i,j   ,k)+f(i,j+1 ,k)))
         end do
        end do
      end do
*
      nflop=nflop+2
      end
************************************************************************
      function zup1 (f)
*
*  f is centered on (i,j,k-.5)
*  zup1(f) is centered on (i,j,k), first order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, zup1
!hpf$ align with TheCube:: f, zup1
*
*-----------------------------------------------------------------------
      if (mz .eq. 1) then
        zup1=f
        return
      endif

      a = .5
*
C*$* ASSERT DO (SERIAL)
      do j=1,my
        do i=1,mx
          zup1(i,j,mz  ) = (
     &                   a*(f(i,j,mz  )+f(i,j,1   )))
        end do
      end do
C$ do across local(i,j,k)
      do k=1,mz-1
        do j=1,my
          do i=1,mx
            zup1(i,j,k) = (
     &                   a*(f(i,j,k   )+f(i,j,k+1 )))
          end do
        end do
      end do

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

      if (mx .eq. 1) then
        xdn1=f
        return
      endif

      a = .5

C$ do across local(i,j,k)
      do k=1,mz
       do j=1,my
        xdn1(1  ,j,k) = (
     &                 a*(f(mx  ,j,k)+f(1   ,j,k)))
       end do
       do j=1,my
        do i=1,mx-1
          xdn1(i+1,j,k) = (
     &                  a*(f(i  ,j,k)+f(i+1 ,j,k)))
        end do
       end do
      end do
*
      nflop=nflop+2
      end
************************************************************************
      function ydn1 (f)
*
*  f is centered on (i,j,k)
*  ydn1(f) is centered on (i,j-.5,k), first order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, ydn1
!hpf$ align with TheCube:: f, ydn1
*
*-----------------------------------------------------------------------
      if (my .eq. 1) then
        ydn1=f
        return
      endif

      a = .5

      if (my .eq. 1) then
	do k=1,mz
	do i=1,mx
	  ydn1(i,1,k)=f(i,1,k)
	enddo
	enddo
	return
      endif
*
C$ do across local(i,j,k)
      do k=1,mz
        do i=1,mx
          ydn1(i,1   ,k) = (
     &                   a*(f(i,my  ,k)+f(i,1   ,k)))
        end do
        do j=1,my-1
         do i=1,mx
          ydn1(i,j+1,k) = (
     &                   a*(f(i,j   ,k)+f(i,j+1 ,k)))
         end do
        end do
      end do
*
      nflop=nflop+2
      end
************************************************************************
      function zdn1 (f)
*
*  f is centered on (i,j,k)
*  zdn1(f) is centered on (i,j,k-.5), first order stagger
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, zdn1
!hpf$ align with TheCube:: f, zdn1
*
*-----------------------------------------------------------------------
      if (mz .eq. 1) then
        zdn1=f
        return
      endif

      a = .5
*
C*$* ASSERT DO (SERIAL)
      do j=1,my
        do i=1,mx
          zdn1(i,j,1   ) = (
     &                   a*(f(i,j,mz  )+f(i,j,1   )))
        end do
      end do
C$ do across local(i,j,k)
      do k=1,mz-1
        do j=1,my
          do i=1,mx
            zdn1(i,j,k+1) = (
     &                   a*(f(i,j,k   )+f(i,j,k+1 )))
          end do
        end do
      end do
*
      nflop=nflop+2
      end
************************************************************************
      function ddxup1 (f)
*
*  f is centered on (i-.5,j,k)
*  ddxup1(f) is centered on (i,j,k), first order stagger
*
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddxup1
!hpf$ align with TheCube:: f, ddxup1

      if (mx .eq. 1) then
        ddxup1=0.
        return
      endif

      a = 1./dx

C$ do across local(i,j,k)
      do k=1,mz
      do j=1,my
        ddxup1(mx  ,j,k) = (
     &                 a*(f(1   ,j,k)-f(mx  ,j,k)))
      end do
      do j=1,my
        do i=1,mx-1
          ddxup1(i  ,j,k) = (
     &                 a*(f(i+1 ,j,k)-f(i   ,j,k)))
        end do
      end do
      end do
*
      nflop=nflop+2
      end
************************************************************************
      function ddyup1 (f)
*
*  f is centered on (i,j-.5,k)
*  ddyup1(f) is centered on (i,j,k), first order stagger
*
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddyup1
!hpf$ align with TheCube:: f, ddyup1

      if (my .eq. 1) then
        ddyup1=0.
        return
      endif

      a = 1./dy

C$ do across local(i,j,k)
      do k=1,mz
        do i=1,mx
          ddyup1(i,my  ,k) = (
     &                   a*(f(i,1   ,k)-f(i,my  ,k)))
        end do
        do j=1,my-1
         do i=1,mx
          ddyup1(i,j  ,k) = (
     &                   a*(f(i,j+1 ,k)-f(i,j   ,k)))
         end do
        end do
      end do
*
      nflop=nflop+2
      end
************************************************************************
      function ddzup1 (f)
*
*  f is centered on (i,j,k-.5)
*  ddzup1(f) is centered on (i,j,k), first order stagger
*
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddzup1
!hpf$ align with TheCube:: f, ddzup1

      if (mz .eq. 1) then
        ddzup1=0.
        return
      endif

      a = 1./dz

C*$* ASSERT DO (SERIAL)
      do j=1,my
        do i=1,mx
          ddzup1(i,j,mz  ) = (
     &                   a*(f(i,j,1   )-f(i,j,mz  )))
        end do
      end do
C$ do across local(i,j,k)
      do k=1,mz-1
        do j=1,my
          do i=1,mx
            ddzup1(i,j,k) = (
     &                   a*(f(i,j,k+1 )-f(i,j,k   )))
          end do
        end do
      end do
*
      nflop=nflop+2
      end
************************************************************************
      function ddxdn1 (f)
*
*  f is centered on (i,j,k)
*  ddxdn1(f) is centered on (i-.5,j,k), first order stagger
*
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddxdn1
!hpf$ align with TheCube:: f, ddxdn1

      if (mx .eq. 1) then
        ddxdn1=0.
        return
      endif

      a = 1./dx

C$ do across local(i,j,k)
      do k=1,mz
      do j=1,my
        ddxdn1(1   ,j,k) = (
     &                 a*(f(1   ,j,k)-f(mx  ,j,k)))
      end do
      do j=1,my
        do i=1,mx-1
          ddxdn1(i+1,j,k) = (
     &                 a*(f(i+1 ,j,k)-f(i   ,j,k)))
        end do
      end do
      end do
*
      nflop=nflop+2
      end
************************************************************************
      function ddydn1 (f)
*
*  f is centered on (i,j,k)
*  ddydn1(f) is centered on (i,j-.5,k), first order stagger
*
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddydn1
!hpf$ align with TheCube:: f, ddydn1

      if (my .eq. 1) then
        ddydn1=0.
        return
      endif

      a = 1./dy

C$ do across local(i,j,k)
      do k=1,mz
        do i=1,mx
          ddydn1(i,1   ,k) = (
     &                   a*(f(i,1   ,k)-f(i,my  ,k)))
        end do
        do j=1,my-1
         do i=1,mx
          ddydn1(i,j+1,k) = (
     &                   a*(f(i,j+1 ,k)-f(i,j   ,k)))
         end do
        end do
      end do
*
      nflop=nflop+2
      end
************************************************************************
      function ddzdn1 (f)
*
*  f is centered on (i,j,k)
*  ddzdn1(f) is centered on (i,j,k-.5), first order stagger
*
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz):: f, ddzdn1
!hpf$ align with TheCube:: f, ddzdn1

      if (mz .eq. 1) then
        ddzdn1=0.
        return
      endif

      a = 1./dz

C*$* ASSERT DO (SERIAL)
      do j=1,my
        do i=1,mx
          ddzdn1(i,j,1   ) = (
     &                   a*(f(i,j,1   )-f(i,j,mz  )))
        end do
      end do
C$ do across local(i,j,k)
      do k=1,mz-1
        do j=1,my
          do i=1,mx
            ddzdn1(i,j,k+1) = (
     &                   a*(f(i,j,k+1 )-f(i,j,k   )))
          end do
        end do
      end do
*
      nflop=nflop+2
      end
*-----------------------------------------------------------------------
      function difx1(f)
*
*  Finite difference one zone left to two zones right.
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, difx1
*-----------------------------------------------------------------------
      if (mx .eq. 1) then
        difx1 = 0.0
        return
      endif

      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
*-----------------------------------------------------------------------
      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
*-----------------------------------------------------------------------
      if (mz .eq. 1) then
	  difz1 = 0.0
	  return
      endif

      difz1 = f - cshift(f,dim=3,shift=1)
      nflop=nflop+1
      END
*-----------------------------------------------------------------------
      function difx2(f)
*
*  Finite difference one zone left to two zones right.
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, difx2
*-----------------------------------------------------------------------
      if (mx .eq. 1) then
        difx2 = 0.0
        return
      endif

      difx2 = cshift(f,dim=1,shift=-1) - cshift(f,dim=1,shift=1)
      nflop=nflop+1
      END
*-----------------------------------------------------------------------
      function dify2(f)
*
*  Finite difference one zone left to two zones right.
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, dify2
*-----------------------------------------------------------------------
      if (my .eq. 1) then
	  dify2 = 0.0
	  return
      endif

      dify2 = cshift(f,dim=2,shift=-1) - cshift(f,dim=2,shift=1)
      nflop=nflop+1
      END
*-----------------------------------------------------------------------
      function difz2(f)
*
*  Finite difference one zone left to two zones right.
*
      include 'cparam.inc'
      real, dimension(mx,my,mz):: f, difz2
*-----------------------------------------------------------------------
      if (mz .eq. 1) then
	  difz2 = 0.0
	  return
      endif

      difz2 = cshift(f,dim=3,shift=-1) - 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

      max3 = f
      if(mx .gt. 1) then
        max3 = amax1(max3,
     &        cshift(max3,dim=1,shift=+1),
     &        cshift(max3,dim=1,shift=-1))
      endif
      if (my .gt. 1) then
        max3 = amax1(max3,
     &        cshift(max3,dim=2,shift=+1),
     &        cshift(max3,dim=2,shift=-1))
      endif
      if (mz .gt. 1) then
        max3 = amax1(max3,
     &      cshift(max3,dim=3,shift=+1),
     &      cshift(max3,dim=3,shift=-1))
      endif
      nflop=nflop+9
      end
*----------------------------------------------------------------------
      function smooth3 (f)
      include 'cparam.inc'
      include 'cdata.inc'
*
      real, dimension (mx,my,mz) :: f, smooth3

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

      nflop=nflop+9
      end
************************************************************************
      subroutine init_derivs ()
*
      end
*-----------------------------------------------------------------------
      subroutine init_stagger ()
*
      end
