!+
! NAME:
!	Deriv_dx
! PURPOSE:
!       Partial Derivative of a 3-d array with respect to the first
!       (x) index
! CALLING SEQUENCE:
!	call Deriv_dx(x, f, nx, ny, nz, use_saved, d, cxarr)
! INPUTS:
!	X:  Variable to differentiate with respect to
!	f:  Function to be differentiated, must be 3-d with the first
!	index having the same number of elements as X
!       nx, ny, nz are the dimensions of the array f
!                  The array x has dimensions of nx
!       use_saved: if set to 1, use the values of cxarr that are there
!       cxarr is an array of coefficieints used in the derivatives,
!                   an array of nx by 3
! OUTPUTS:
!	d = the derivative df/dx at each point
! HISTORY:
!	17-May-2005, jmm, jimm@ssl.berkeley.edu, from deriv_dx.pro
!       2-mar-2007, implicit none, and assignment statements added by Dave Bercik
!-
      subroutine deriv_dx(x, f, nx, ny, nz, use_saved, d, cxarr)
      implicit none

      integer :: nx, ny, nz, use_saved
      integer :: j, k, l, nx1
      real*8 :: x(nx), f(nx, ny, nz), d(nx, ny, nz), cxarr(nx, 3)
      real*8,allocatable:: x01(:), x02(:), x12(:)
      real*8, save :: sx1, sx2, sx3, rx1, rx2, rx3
      nx1 = nx-1
! dy/dx=y0*(2x-x1-x2)/(x01*x02)+y1*(2x-x0-x2)/(x10*x12)+y2*(2x-x0-x1)/(x20*x21)
! Where: x01 = x0-x1, x02 = x0-x2, x12 = x1-x2, etc.
! where x1 is x at the point you get the derivative, x0 is the previous point
! and x2 is the next point
      if(use_saved.eq.0) then
!initialize grids
         allocate(x01(nx), x02(nx), x12(nx))
         x12(1:nx) = 0.0
         x01(1:nx) = 0.0
         x02(1:nx) = 0.0
         cxarr(1:nx,1:3) = 0.0
         do 10 j = 2, nx-1
            x01(j) = x(j-1)-x(j)
            x02(j) = x(j-1)-x(j+1)
            x12(j) = x(j)-x(j+1)
!cx's are obtained from x01,x02,x12
            cxarr(j,1) = (x12(j)/(x01(j)*x02(j)))
            cxarr(j,2) = (1./x12(j) - 1./x01(j))
            cxarr(j,3) = (x01(j)/(x02(j)*x12(j)))
 10      continue
         sx1 = (x01(2)+x02(2))/(x01(2)*x02(2))
         sx2 = x02(2)/(x01(2)*x12(2))
         sx3 = x01(2)/(x02(2)*x12(2))
         rx1 = x12(nx1)/(x01(nx1)*x02(nx1))
         rx2 = x02(nx1)/(x01(nx1)*x12(nx1))
         rx3 = (x02(nx1)+x12(nx1))/(x02(nx1)*x12(nx1))
         deallocate(x01, x02, x12)
      end if
!Do the derivative
      forall(j=2:nx-1, k=1:ny, l=1:nz)
            d(j, k, l) = cxarr(j, 1)*f(j-1, k, l) + 
     &           cxarr(j,2)*f(j, k, l) - cxarr(j,3)*f(j+1, k, l)
      end forall
      forall(k=1:ny, l=1:nz)
!Formulae for the first and last points:
            d(1, k, l) = f(1, k, l)*sx1 -
     &           f(2, k, l)*sx2 +  f(3, k, l)*sx3
            d(nx, k, l) = -f(nx-2, k, l)*rx1 +
     &           f(nx-1, k, l)*rx2 - f(nx, k, l)*rx3
      end forall
      Return
      End
!+
! NAME:
!	Deriv_dy
! PURPOSE:
!       Partial Derivative of a 3-d array with respect to the second
!       (y) index
! CALLING SEQUENCE:
!	call Deriv_dy(y, f, nx, ny, nz, use_saved, d, cyarr)
! INPUTS:
!	y:  Variable to differentiate with respect to
!	f:  Subroutine to be differentiated, must be 3-d with the second
!	index having the same number of elements as y
!       nx, ny, nz are the dimensions of the array f
!                           The array y has dimensions of ny
!       use_saved: if set to 1, use current values of cyarr
!       cyarr are coffeicients for the derivative
! OUTPUTS:
!	d = the derivative df/dy at each point
! HISTORY:
!	17-May-2005, jmm, jimm@ssl.berkeley.edu, from deriv_dx.pro
!       2-mar-2007, implicit none, and assignment statements added by Dave Bercik
!-
      subroutine deriv_dy(y, f, nx, ny, nz, use_saved, d, cyarr)
      implicit none

      integer :: nx, ny, nz, use_saved
      integer :: j, k, l, ny1
      real*8 :: y(ny), f(nx, ny, nz), d(nx, ny, nz), cyarr(ny, 3)
      real*8, allocatable:: y01(:), y02(:), y12(:)
      real*8, save :: sy1, sy2, sy3, ry1, ry2, ry3
      ny1 = ny-1
! dy/dx=y0*(2x-x1-x2)/(x01*x02)+y1*(2x-x0-x2)/(x10*x12)+y2*(2x-x0-x1)/(x20*x21)
! Where: x01 = x0-x1, x02 = x0-x2, x12 = x1-x2, etc.
! where x1 is x at the point you get the derivative, x0 is the previous point
! and x2 is the next point
      if(use_saved.eq.0) then
!initialize grids
         allocate(y01(ny), y02(ny), y12(ny))
         y12(1:ny) = 0.0
         y01(1:ny) = 0.0
         y02(1:ny) = 0.0
         cyarr(1:ny,1:3) = 0.0
         do 10 k = 2, ny-1
            y12(k) = y(k)-y(k+1)
            y01(k) = y(k-1)-y(k)
            y02(k) = y(k-1)-y(k+1)
!cy's are obtained from y01,y02,y12
            cyarr(k, 1) = (y12(k)/(y01(k)*y02(k)))
            cyarr(k, 2) = (1./y12(k) - 1./y01(k))
            cyarr(k, 3) = (y01(k)/(y02(k)*y12(k)))
 10      continue
         sy1 = (y01(2)+y02(2))/(y01(2)*y02(2))
         sy2 = y02(2)/(y01(2)*y12(2))
         sy3 = y01(2)/(y02(2)*y12(2))
         ry1 = y12(ny1)/(y01(ny1)*y02(ny1))
         ry2 = y02(ny1)/(y01(ny1)*y12(ny1))
         ry3 = (y02(ny1)+y12(ny1))/(y02(ny1)*y12(ny1))
         deallocate(y01, y02, y12)
      end if
!Do the derivative
      forall (j=1:nx, k=2:ny-1, l=1:nz)
            d(j, k, l) = cyarr(k, 1)*f(j, k-1, l) + 
     &           cyarr(k, 2)*f(j, k, l) - cyarr(k, 3)*f(j, k+1, l)
      end forall
! Formulae for the first and last points:
      forall (j=1:nx,  l=1:nz)
            d(j, 1, l) = f(j, 1, l)*sy1 -
     &           f(j, 2, l)*sy2 + f(j, 3, l)*sy3
            d(j, ny, l) = -f(j, ny-2, l)*ry1 +
     &           f(j, ny-1, l)*ry2 - f(j, ny, l)*ry3
      end forall
      Return
      End
!+
! NAME:
!	deriv_dz
! PURPOSE:
!       Partial Derivative of a 3-d array with respect to the third
!       (z) index
! CALLING SEQUENCE:
!	call deriv_dz(z, f, nx, ny, nz, use_saved, d, czarr)
! INPUTS:
!	z:  Variable to differentiate with respect to
!	f:  Function to be differentiated, must be 3-d with the third
!	index having the same number of elements as z
!       nx, ny, nz are the dimensions of the array f
!                           The array z has dimensions of nz
!       use_saved: if set to 1, use the current values of czarr 
!       czarr are coefficients for the derivative
! OUTPUTS:
!	d = the derivative df/dz at each point
! HISTORY:
!	17-May-2005, jmm, jimm@ssl.berkeley.edu, from deriv_dx.pro
!       2-mar-2007, implicit none, and assignment statements added by Dave Bercik
!-
      subroutine deriv_dz(z, f, nx, ny, nz, use_saved, d, czarr)
      implicit none

      integer :: nx, ny, nz, use_saved
      integer :: j, k, l, nz1
      real*8 :: z(nz), f(nx, ny, nz), d(nx, ny, nz), czarr(nz, 3)
      real*8, allocatable:: z01(:), z02(:), z12(:)
      real*8, save :: sz1, sz2, sz3, rz1, rz2, rz3
      nz1 = nz-1
! dy/dx=y0*(2x-x1-x2)/(x01*x02)+y1*(2x-x0-x2)/(x10*x12)+y2*(2x-x0-x1)/(x20*x21)
! Where: x01 = x0-x1, x02 = x0-x2, x12 = x1-x2, etc.
! where x1 is x at the point you get the derivative, x0 is the previous point
! and x2 is the next point
      if(use_saved.eq.0) then
         allocate(z01(nz), z02(nz), z12(nz))
         z12(1:nz) = 0.0
         z01(1:nz) = 0.0
         z02(1:nz) = 0.0
         czarr(1:nz, 1:3) = 0.0
         do 10 l = 2, nz-1
            z12(l) = z(l)-z(l+1)
            z01(l) = z(l-1)-z(l)
            z02(l) = z(l-1)-z(l+1)
!cz's are obtained from z01,z02,z12
            czarr(l,1) = (z12(l)/(z01(l)*z02(l)))
            czarr(l,2) = (1./z12(l) - 1./z01(l))
            czarr(l,3) = (z01(l)/(z02(l)*z12(l)))
 10      continue
         sz1 = (z01(2)+z02(2))/(z01(2)*z02(2))
         sz2 = z02(2)/(z01(2)*z12(2))
         sz3 = z01(2)/(z02(2)*z12(2))
         rz1 = z12(nz1)/(z01(nz1)*z02(nz1))
         rz2 = z02(nz1)/(z01(nz1)*z12(nz1))
         rz3 = (z02(nz1)+z12(nz1))/(z02(nz1)*z12(nz1))
         deallocate(z01, z02, z12)
      end if
!Do the derivative
      forall (j=1:nx, k=1:ny, l=2:nz-1)
            d(j, k, l) = czarr(l,1)*f(j, k, l-1) + 
     &           czarr(l,2)*f(j, k, l) - czarr(l,3)*f(j, k, l+1)
      end forall
! Formulae for the first and last points:
      forall (j=1:nx, k=1:ny)
            d(j, k, 1) = f(j, k, 1)*sz1 -
     &           f(j, k, 2)*sz2 + f(j, k, 3)*sz3
            d(j, k, nz) = -f(j, k, nz-2)*rz1 +
     &           f(j, k, nz-1)*rz2 - f(j, k, nz)*rz3
      end forall
      Return
      End
