*     Randall R. Correll
*     Mijan F. Huq
*     Scott A. Klasky
*     Richard A. Matzner
*
*     COPYRIGHT 1995-1996

      subroutine diff(nx,ny,nz,h,f,df,dir,m,d,nl)
      implicit none
      integer nx,ny,nz,nl
      real*8 h
      real*8 f(nx,ny,nz)
      real*8 df(nx,ny,nz)
      real*8 d(3*nl)
      integer m(3*nl)
      character dir

      call diff_int(nx,ny,nz,h,f,df,dir)
      call diff_ob(nx,ny,nz,h,f,df,dir)

      if (nl.gt.1) call diff_ib(nx,ny,nz,h,f,df,dir,m,d,nl)

      return
      end

      subroutine diff_int(nx,ny,nz,h,f,df,dir)
      implicit none
      integer nx,ny,nz
      real*8 h
      real*8 f(nx,ny,nz)
      real*8 df(nx,ny,nz)
      character dir

      if (dir.eq.'x') then
         call diffx_int(nx,ny,nz,h,f,df)
      else if (dir.eq.'y') then
         call diffy_int(nx,ny,nz,h,f,df)
      else if (dir.eq.'z') then
         call diffz_int(nx,ny,nz,h,f,df)
      else
         call err('diff_int','wrong derivative')   
      end if

      return
      end

      subroutine diff_ob(nx,ny,nz,h,f,df,dir)
      implicit none
      integer nx,ny,nz
      real*8 h
      real*8 f(nx,ny,nz)
      real*8 df(nx,ny,nz)
      character dir

      if (dir.eq.'x') then
         call diffx_ob(nx,ny,nz,h,f,df)
      else if (dir.eq.'y') then
         call diffy_ob(nx,ny,nz,h,f,df)
      else if (dir.eq.'z') then
         call diffz_ob(nx,ny,nz,h,f,df)
      else
         call err('diff_ob','wrong derivative')   
      end if

      return
      end


      subroutine diff_ib(nx,ny,nz,h,f,df,dir,m,d,nl)
      implicit none
      integer nx,ny,nz,nl
      real*8 h
      real*8 f(nx,ny,nz)
      real*8 df(nx,ny,nz)
      real*8 d(3*nl)
      integer m(3*nl)
      character dir

      if (dir.eq.'x') then
         call diffx_ib(nx,ny,nz,h,f,df,m,d,nl)
      else if (dir.eq.'y') then
         call diffy_ib(nx,ny,nz,h,f,df,m,d,nl)
      else if (dir.eq.'z') then
         call diffz_ib(nx,ny,nz,h,f,df,m,d,nl)
      else
         call err('diff_ib','wrong derivative')
      end if

      return
      end


      subroutine diffx_int(nx,ny,nz,h,f,df)
      implicit none
      integer nx,ny,nz
      real*8 h,f(nx,ny,nz),df(nx,ny,nz),h2
      integer i,j,k

      h2 = 0.5e0 / h

      do k = 2 , nz - 1
         do j = 2 , ny - 1
            do i = 2 , nx - 1

               df(i,j,k) = h2 * ( f(i+1,j,k) - f(i-1,j,k) )
            end do
         end do
      end do

      return
      end

      subroutine diffy_int(nx,ny,nz,h,f,df)
      implicit none
      integer nx,ny,nz
      real*8 h,f(nx,ny,nz),df(nx,ny,nz),h2
      integer i,j,k
      
      h2 = 0.5e0 / h

      do k = 2 , nz - 1
         do j = 2 , ny - 1
            do i = 2 , nx - 1

               df(i,j,k) = h2 * ( f(i,j+1,k) - f(i,j-1,k) )
            end do                                          
         end do
      end do

      return
      end

      subroutine diffz_int(nx,ny,nz,h,f,df)
      implicit none
      integer nx,ny,nz
      real*8 h,f(nx,ny,nz),df(nx,ny,nz),h2
      integer i,j,k
      
      h2 = 0.5e0 / h

      do k = 2 , nz - 1
         do j = 2 , ny - 1
            do i = 2 , nx - 1

               df(i,j,k) = h2 * ( f(i,j,k+1) - f(i,j,k-1) )
            end do                                          
         end do
      end do

      return
      end


      subroutine diffx_ob(nx,ny,nz,h,f,df)
      implicit none
      integer nx,ny,nz
      real*8 h,f(nx,ny,nz),df(nx,ny,nz),h2,hi
      integer i,j,k

      hi = 1.0e0/h
      h2 = 0.5e0 / h

      i = 1
      do k = 1 , nz
         do j = 1 , ny
            df(i,j,k) = h2 * ( -3.0e0*f(i,j,k) 
     $                + 4.0e0*f(i+1,j,k) - f(i+2,j,k) )
         end do
      end do

      i = nx 
      do k = 1 , nz
         do j = 1 , ny
            df(i,j,k) = h2 * (  3.0e0*f(i,j,k)
     $                - 4.0e0*f(i-1,j,k) + f(i-2,j,k) )
         end do
      end do

      do k = 1 , nz 
         do i = 2 , nx-1
            df(i,1,k) = h2 *(f(i+1,1,k)-f(i-1,1,k) )
            df(i,ny,k) = h2 *(f(i+1,ny,k)-f(i-1,ny,k) )
         end do
      end do

      do j = 1 , ny
         do i = 2 , nx-1
            df(i,j,1) = h2 *(f(i+1,j,1)-f(i-1,j,1) )
            df(i,j,nz) = h2 *(f(i+1,j,nz)-f(i-1,j,nz) )
         end do
      end do

      return
      end

      subroutine diffy_ob(nx,ny,nz,h,f,df)
      implicit none
      integer nx,ny,nz
      real*8 h,f(nx,ny,nz),df(nx,ny,nz),h2,hi
      integer i,j,k

      hi = 1.0e0/h
      h2 = 0.5e0 / h

      j = 1
      do k = 1 , nz
         do i = 1 , nx
            df(i,j,k) = h2 * ( -3.0e0*f(i,j,k) 
     $                + 4.0e0*f(i,j+1,k) - f(i,j+2,k) )
         end do
      end do

      j = nx 
      do k = 1 , nz
         do i = 1 , nx
            df(i,j,k) = h2 * (  3.0e0*f(i,j,k)
     $                - 4.0e0*f(i,j-1,k) + f(i,j-2,k) )
         end do
      end do

      do k = 1 , nz
         do j = 2 , ny-1
            df(1,j,k) = h2 *(f(1,j+1,k)-f(1,j-1,k) )
            df(nx,j,k) = h2 *(f(nx,j+1,k)-f(nx,j-1,k) )
         end do
      end do

      do i = 1 , ny
         do j = 2 , nx -1
            df(i,j,1) = h2 *(f(i,j+1,1)-f(i,j-1,1) )
            df(i,j,nz) = h2 *(f(i,j+1,nz)-f(i,j-1,nz) )
         end do
      end do


      return
      end

      subroutine diffz_ob(nx,ny,nz,h,f,df)
      implicit none
      integer nx,ny,nz
      real*8 h,f(nx,ny,nz),df(nx,ny,nz),h2,hi
      integer i,j,k

      hi = 1.0e0/h
      h2 = 0.5e0 / h

      k = 1
      do j = 1 , ny
         do i = 1 , nx
            df(i,j,k) = h2 * ( -3.0e0*f(i,j,k) 
     $                + 4.0e0*f(i,j,k+1) - f(i,j,k+2) )
         end do
      end do

      k = nz
      do j = 1 , ny
         do i = 1 , nx
            df(i,j,k) = h2 * (  3.0e0*f(i,j,k)
     $                - 4.0e0*f(i,j,k-1) + f(i,j,k-2) )
         end do
      end do

      do i = 1 , nz
         do k = 2 , nx-1
            df(i,1,k) = h2 *(f(i,1,k+1)-f(i,1,k-1) )
            df(i,ny,k) = h2 *(f(i,ny,k+1)-f(i,ny,k-1) )
         end do
      end do

      do j = 1 , nz
         do k = 2 , ny-1
            df(1,j,k) = h2 *(f(1,j,k+1)-f(1,j,k-1) )
            df(nx,j,k) = h2 *(f(nx,j,k+1)-f(nx,j,k-1) )
         end do
      end do

      return
      end
C------------------------------------------------------------------------------
      subroutine diffy_ib(nx,ny,nz,h,f,df,m,d,nl)
      implicit none
      integer nx,ny,nz,nl
      real*8 h,f(nx,ny,nz),df(nx,ny,nz),h2
      real*8 d(3*nl)
      integer m(3*nl)
      integer i,j,k,l
		real*8 d1,d2,d3,da

      h2 = 0.5e0 / h

		do l = 1 , nl

			i  = m(l)
			j  = m(nl+l)
			k  = m(2*nl+l)
			d1 = d(l)
			d2 = d(nl+l)
			d3 = d(2*nl+l)
			da = abs(d2)
         df(i,j,k) = h2*(d2*(-3.0e0*f(i,j,k)+4.0e0*f(i,j+int(d2),k)
     $                       -      f(i,j+2*int(d2),k)) +
     $               (1.0e0-da)*(f(i,j+1,k)-f(i,j-1,k) ) )

		end do

		return
		end
C------------------------------------------------------------------------------
      subroutine diffx_ib(nx,ny,nz,h,f,df,m,d,nl)
      implicit none
      integer nx,ny,nz,nl
      real*8 h,f(nx,ny,nz),df(nx,ny,nz),h2
      real*8 d(3*nl)
      integer m(3*nl)
      integer i,j,k,l
		real*8 d1,d2,d3,da

      h2 = 0.5e0 / h

		do l = 1 , nl

			i  = m(l)
			j  = m(nl+l)
			k  = m(2*nl+l)
			d1 = d(l)
			d2 = d(nl+l)
			d3 = d(2*nl+l)
			da = abs(d1)
         df(i,j,k) = h2*(d1*(-3.0e0*f(i,j,k)+4.0e0*f(i+int(d1),j,k)
     $                       -      f(i+2*int(d1),j,k)) +
     $               (1.0e0-da)*(f(i+1,j,k)-f(i-1,j,k) ) )

		end do

		return
		end
C------------------------------------------------------------------------------
      subroutine diffz_ib(nx,ny,nz,h,f,df,m,d,nl)
      implicit none
      integer nx,ny,nz,nl
      real*8 h,f(nx,ny,nz),df(nx,ny,nz),h2
      real*8 d(3*nl)
      integer m(3*nl)
      integer i,j,k,l
		real*8 d1,d2,d3,da

      h2 = 0.5e0 / h

		do l = 1 , nl

			i  = m(l)
			j  = m(nl+l)
			k  = m(2*nl+l)
			d1 = d(l)
			d2 = d(nl+l)
			d3 = d(2*nl+l)
			da = abs(d3)
         df(i,j,k) = h2*(d3*(-3.0e0*f(i,j,k)+4.0e0*f(i,j,k+int(d3))
     $                       -      f(i,j,k+2*int(d3))) +
     $               (1.0e0-da)*(f(i,j,k+1)-f(i,j,k-1) ) )

		end do

		return
		end
