c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE AAMULT                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE AAMULT(a,b,nx,ny,nz,
     :           ibgn,iend,jbgn,jend,kbgn,kend, ab)

c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the element-wise product of arrays a and b.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
c
c     6/01/2 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       First dimension of arrays a, b and ab
c       ny       Second dimension of arrays a, b and ab
c       nz       Third dimension of arrays a, b and ab
c
c       a        First multiplier array
c       b        Second multiplier array
c
c       ibgn     Index in first dimension to begin multiplication
c       iend     Index in first dimension to end multiplication
c       jbgn     Index in second dimension to begin multiplication
c       jend     Index in second dimension to end multiplication
c       kbgn     Index in third dimension to begin multiplication
c       kend     Index in third dimension to end multiplication
c
c     OUTPUT:
c
c       ab       Element-wise product a*b over range specified by 
c                the starting and ending indices
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx, ny, nz        ! Number of grid points in 3 directions
c
      real a (nx,ny,nz)         ! Input array 1
      real b (nx,ny,nz)         ! Input array 2
c
      integer ibgn,iend,jbgn,jend,kbgn,kend     
c
      real ab(nx,ny,nz)         ! Product array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,j,k
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c     Beginning of executable code...
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      DO 100 k=kbgn,kend
      DO 100 j=jbgn,jend
      DO 100 i=ibgn,iend

        ab(i,j,k)=a(i,j,k)*b(i,j,k)

100   CONTINUE

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE AVGX                       ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE AVGX(a, onvf, nx,ny,nz,
     :           ibgn,iend,jbgn,jend,kbgn,kend, aavg)

c
c#######################################################################
c
c     PURPOSE:
c
c     Perform a spatial average operation on array a in the x direction. 
c     The averaged variable is defined at the midpoint between the two 
c     points for the two input values.
c     If the averaged variable aavg is on the grid volume face, onvf =1
c     otherwise, onvf = 0
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
c
c     6/01/2 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       First dimension of arrays a and aavg
c       ny       Second dimension of arrays a and aavg
c       nz       Third dimension of arrays a and aavg
c
c       a        Input array 
c
c       onvf     Integer grid point indicator
c                If the averaged variable aavg is on the grid volume face, 
c                onvf =1; otherwise, onvf = 0
c
c       ibgn     Index in first dimension to begin multiplication
c       iend     Index in first dimension to end multiplication
c       jbgn     Index in second dimension to begin multiplication
c       jend     Index in second dimension to end multiplication
c       kbgn     Index in third dimension to begin multiplication
c       kend     Index in third dimension to end multiplication
c
c     OUTPUT:
c
c       aavg     Result of average operation on array a in x direction
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx, ny, nz        ! Number of grid points in 3 directions
c
      real a   (nx,ny,nz)       ! Input array
      integer onvf
c
      integer ibgn,iend,jbgn,jend,kbgn,kend
c
      real aavg(nx,ny,nz)       ! Averaged array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,j,k,iright,ileft
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c     Beginning of executable code...
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      IF ( onvf.eq.1) THEN

        iright = 0
        ileft = -1

      ELSE

        iright = 1
        ileft = 0

      END IF

      DO 100 k=kbgn,kend
      DO 100 j=jbgn,jend
      DO 100 i=ibgn,iend

        aavg(i,j,k)=(a(i+iright,j,k)
     :              +a(i+ileft ,j,k))*0.5

100   CONTINUE

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE AVGY                       ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE AVGY(a, onvf, nx,ny,nz,
     :           ibgn,iend,jbgn,jend,kbgn,kend, aavg)

c
c#######################################################################
c
c     PURPOSE:
c
c     Perform a spatial average operation on array a in the y direction.
c     The averaged variable is defined at the midpoint between the two 
c     points for the two input values.
c     If the averaged variable aavg is on the grid volume face, onvf =1
c     otherwise, onvf = 0
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
c
c     6/01/2 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       First dimension of arrays a and aavg
c       ny       Second dimension of arrays a and aavg
c       nz       Third dimension of arrays a and aavg
c
c       a        Input array
c
c       onvf     Integer grid point indicator
c                If the averaged variable aavg is on the grid volume face, 
c                onvf =1; otherwise, onvf = 0
c
c       ibgn     Index in first dimension to begin operation
c       iend     Index in first dimension to end operation
c       jbgn     Index in second dimension to begin operation
c       jend     Index in second dimension to end operation
c       kbgn     Index in third dimension to begin operation
c       kend     Index in third dimension to end operation
c
c     OUTPUT:
c
c       aavg     Result of average operation on array a in y direction
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx, ny, nz        ! Number of grid points in 3 directions
c
      real a   (nx,ny,nz)       ! Input array
      integer onvf
c
      integer ibgn,iend,jbgn,jend,kbgn,kend
c
      real aavg(nx,ny,nz)       ! Average array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,j,k,jright,jleft
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c     Beginning of executable code...
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      IF ( onvf.eq.1) THEN

        jright = 0
        jleft = -1

      ELSE

        jright = 1
        jleft = 0

      END IF

      DO 100 k=kbgn,kend
      DO 100 i=ibgn,iend
      DO 100 j=jbgn,jend

        aavg(i,j,k)=(a(i,j+jright,k)+a(i,j+jleft,k))*0.5

100   CONTINUE

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE AVGZ                       ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE AVGZ(a, onvf, nx,ny,nz,
     :           ibgn,iend,jbgn,jend,kbgn,kend, aavg)

c
c#######################################################################
c
c     PURPOSE:
c
c     Perform a spatial average operation on array a in the z direction.
c     The averaged variable is defined at the midpoint between the two 
c     points for the two input values.
c     If the averaged variable aavg is on the grid volume face, onvf =1
c     otherwise, onvf = 0
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
c
c     6/01/2 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       First dimension of arrays a and aavg
c       ny       Second dimension of arrays a and aavg
c       nz       Third dimension of arrays a and aavg
c
c       a        Input array
c
c       onvf     Integer grid point indicator
c                If the averaged variable aavg is on the grid volume face, 
c                onvf =1; otherwise, onvf = 0
c
c       ibgn     Index in first dimension to begin operation
c       iend     Index in first dimension to end operation
c       jbgn     Index in second dimension to begin operation
c       jend     Index in second dimension to end operation
c       kbgn     Index in third dimension to begin operation
c       kend     Index in third dimension to end operation
c
c     OUTPUT:
c
c       aavg     Result of average operation on array a in z direction
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx, ny, nz        ! Number of grid points in 3 directions
c
      real a   (nx,ny,nz)       ! Input array
      integer onvf
c
      integer ibgn,iend,jbgn,jend,kbgn,kend
c
      real aavg(nx,ny,nz)       ! Average array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,j,k,kup,kdown
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c     Beginning of executable code...
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      IF ( onvf.eq.1) THEN

        kup = 0
        kdown = -1

      ELSE

        kup = 1
        kdown = 0

      END IF

      DO 100 k=kbgn,kend
      DO 100 j=jbgn,jend
      DO 100 i=ibgn,iend

        aavg(i,j,k)=(a(i,j,k+kup)+a(i,j,k+kdown ))*0.5

100   CONTINUE

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE DIFX                       ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE DIFX(a, onvf, nx,ny,nz,
     :           ibgn,iend,jbgn,jend,kbgn,kend, dx, adifx)

c
c#######################################################################
c
c     PURPOSE:
c
c     Perform finite difference operation on an array a in the x direction.
c     adifx = d( a ) / dx. The output variable is defined at the midpoint 
c     between the two points whose values are differenced.
c     If the output variable adifx is defined on the grid volume face when 
c     onvf =1; otherwise, onvf = 0.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
c
c     6/01/2 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       First dimension of arrays a and adifx
c       ny       Second dimension of arrays a and adifx
c       nz       Third dimension of arrays a and adfix
c
c       a        Input array
c
c       onvf     Integer grid point indicator
c                If the averaged variable aavg is on the grid volume face, 
c                onvf =1; otherwise, onvf = 0
c
c       ibgn     Index in first dimension to begin operation
c       iend     Index in first dimension to end operation
c       jbgn     Index in second dimension to begin operation
c       jend     Index in second dimension to end operation
c       kbgn     Index in third dimension to begin operation
c       kend     Index in third dimension to end operation
c
c       dx       Grid spacing in x direction (m)
c
c     OUTPUT:
c
c       adifx    Differenced array del(a)/delx
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx, ny, nz        ! Number of grid points in 3 directions
c
      real a    (nx,ny,nz)      ! Input array
      integer onvf
c
      integer ibgn,iend,jbgn,jend,kbgn,kend
c
      real dx
      real adifx(nx,ny,nz)      
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,j,k,iright,ileft
      real dxinv
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c     Beginning of executable code...
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      IF ( onvf.eq.1) THEN

        iright = 0
        ileft = -1

      ELSE

        iright = 1
        ileft = 0

      END IF

      dxinv = 1.0/dx

      DO 100 k=kbgn,kend
      DO 100 j=jbgn,jend
      DO 100 i=ibgn,iend

        adifx(i,j,k)=(a(i+iright,j,k)-a(i+ileft ,j,k))*dxinv

100   CONTINUE

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE DIFY                       ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE DIFY(a, onvf, nx,ny,nz,
     :           ibgn,iend,jbgn,jend,kbgn,kend, dy, adify)

c
c#######################################################################
c
c     PURPOSE:
c
c     Perform finite difference operation on an array a in the y direction.
c     adify = d( a ) / dy. The output variable is defined at the midpoint 
c     between the two points whose values are differenced.
c     If the output variable adify is defined on the grid volume face when 
c     onvf =1; otherwise, onvf = 0.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
c
c     6/01/2 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       First dimension of arrays a and adify
c       ny       Second dimension of arrays a and adify
c       nz       Third dimension of arrays a and adify
c
c       a        Input array 
c
c       onvf     Integer grid point indicator
c                If the averaged variable aavg is on the grid volume face, 
c                onvf =1; otherwise, onvf = 0
c
c       ibgn     Index in first dimension to begin operation
c       iend     Index in first dimension to end operation
c       jbgn     Index in second dimension to begin operation
c       jend     Index in second dimension to end operation
c       kbgn     Index in third dimension to begin operation
c       kend     Index in third dimension to end operation
c
c       dy       Grid spacing in y direction (m)
c
c     OUTPUT:
c
c       adify    Differenced array del(a)/dely
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx, ny, nz        ! Number of grid points in 3 directions
c
      real a   (nx,ny,nz)       ! Input array
      integer onvf
c
      integer ibgn,iend,jbgn,jend,kbgn,kend
c
      real dy
      real adify(nx,ny,nz)      
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,j,k,jright,jleft
      real dyinv
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c     Beginning of executable code...
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      IF ( onvf.eq.1) THEN

        jright = 0
        jleft = -1

      ELSE

        jright = 1
        jleft = 0

      END IF

      dyinv = 1.0/dy

      DO 100 k=kbgn,kend
      DO 100 i=ibgn,iend
      DO 100 j=jbgn,jend

        adify(i,j,k)=(a(i,j+jright,k)-a(i,j+jleft,k))*dyinv

100   CONTINUE

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE DIFZ                       ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE DIFZ(a, onvf, nx,ny,nz,
     :           ibgn,iend,jbgn,jend,kbgn,kend, dz, adifz)

c
c#######################################################################
c
c     PURPOSE:
c
c     Perform finite difference operation on an array a in the z direction.
c     adifz = d( a ) / dz. The output variable is defined at the midpoint 
c     between the two points whose values are differenced.
c     If the output variable adifz is defined on the grid volume face when 
c     onvf =1; otherwise, onvf = 0.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
c
c     6/01/2 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       First dimension of arrays a and adifz
c       ny       Second dimension of arrays a and adifz
c       nz       Third dimension of arrays a and adifz
c
c       a        Input array 
c
c       onvf     Integer grid point indicator
c                If the averaged variable aavg is on the grid volume face, 
c                onvf =1; otherwise, onvf = 0
c
c       ibgn     Index in first dimension to begin operation
c       iend     Index in first dimension to end operation
c       jbgn     Index in second dimension to begin operation
c       jend     Index in second dimension to end operation
c       kbgn     Index in third dimension to begin operation
c       kend     Index in third dimension to end operation
c
c     OUTPUT:
c
c       adifz    Differenced array del(a)/delz
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx, ny, nz        ! Number of grid points in 3 directions
c
      real a    (nx,ny,nz)      ! input array
      integer onvf
c
      integer ibgn,iend,jbgn,jend,kbgn,kend
c
      real dz
      real adifz(nx,ny,nz)      ! Differenced array del(a)/delz
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,j,k,kup,kdown
      real dzinv
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c     Beginning of executable code...
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      IF ( onvf.eq.1) THEN

        kup = 0
        kdown = -1

      else

        kup = 1
        kdown = 0

      END IF

      dzinv = 1.0/dz

      DO 100 k=kbgn,kend
      DO 100 j=jbgn,jend
      DO 100 i=ibgn,iend

        adifz(i,j,k)=(a(i ,j,k+kup)-a(i ,j,k+kdown ))*dzinv

100   CONTINUE

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE DIFXX                      ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE DIFXX(a, nx,ny,nz,
     :           ibgn,iend,jbgn,jend,kbgn,kend, dx, adifxx)

c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the second order the difference of an array a in the x direction.
c     The operator is defined as del**2(a)/delx**2.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
c
c     6/01/2 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       First dimension of arrays a and adifxx
c       ny       Second dimension of arrays a and adifxx
c       nz       Third dimension of arrays a and adifxx
c
c       a        Input array 
c
c       ibgn     Index in first dimension to begin operation
c       iend     Index in first dimension to end operation
c       jbgn     Index in second dimension to begin operation
c       jend     Index in second dimension to end operation
c       kbgn     Index in third dimension to begin operation
c       kend     Index in third dimension to end operation
c
c       dx       Grid spacing in x direction (m)
c
c     OUTPUT:
c
c       adifxx   Differenced array del**2(a)/delx**2
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
c     implicit none
c
      integer nx, ny, nz        ! Number of grid points in 3 directions
c
      real a    (nx,ny,nz)      ! Input array
c
      integer ibgn,iend,jbgn,jend,kbgn,kend
c
      real dx
      real adifxx(nx,ny,nz)     ! Differenced array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,j,k
      real dxinv2
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c     Beginning of executable code...
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      dxinv2 = 1.0/(dx*dx)

      DO 100 k=kbgn,kend
      DO 100 j=jbgn,jend
      DO 100 i=ibgn,iend
c
c#######################################################################
c
c     NOTE:
c   
c     The order of calucation in the following formula should not
c     be changed. The order affects the machine trunction error
c     in the result.
c 
c#######################################################################
c
        adifxx(i,j,k)=((a(i+1,j,k)-a(i,j,k))-(a(i,j,k)-a(i-1,j,k)))
     :               *dxinv2

100   CONTINUE

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE DIFYY                      ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE DIFYY(a, nx,ny,nz,
     :           ibgn,iend,jbgn,jend,kbgn,kend, dy, adifyy)

c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the second order finite difference of an array a in the y direction.
c     The operator is defined as del**2(a)/dely**2.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
c
c     6/01/2 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       First dimension of arrays a and adifyy
c       ny       Second dimension of arrays a and adifyy
c       nz       Third dimension of arrays a and adifyy
c
c       a        input array 
c
c       ibgn     index in first dimension to begin operation
c       iend     index in first dimension to end operation
c       jbgn     index in second dimension to begin operation
c       jend     index in second dimension to end operation
c       kbgn     index in third dimension to begin operation
c       kend     index in third dimension to end operation
c
c       dy       grid spacing in y direction (m)
c
c     OUTPUT:
c
c       adifyy   difference array del**2(a)/dely**2
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx, ny, nz         ! Number of grid points in 3 directions
c
      real a     (nx,ny,nz)      ! Input array
c
      integer ibgn,iend,jbgn,jend,kbgn,kend
c
      real dy
      real adifyy(nx,ny,nz)      ! Differenced array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,j,k
      real dyinv2
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c     Beginning of executable code...
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      dyinv2 = 1.0/(dy*dy)

      DO 100 k=kbgn,kend
      DO 100 i=ibgn,iend
      DO 100 j=jbgn,jend
c
c#######################################################################
c
c     NOTE:
c   
c     The order of calucation in the following formula should not
c     be changed. The order affects the machine trunction error
c     in the result.
c 
c#######################################################################
c
        adifyy(i,j,k)=((a(i,j+1,k)-a(i,j,k))-(a(i,j,k)-a(i,j-1,k)))
     :               *dyinv2

100   CONTINUE

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE DIFZZ                      ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE DIFZZ(a, nx,ny,nz,
     :           ibgn,iend,jbgn,jend,kbgn,kend, dz, adifzz)

c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the second order finite difference of an array a in the z direction.
c     The operator is defined as del**2(a)/delz**2.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
c
c     6/01/2 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       First dimension of arrays a and adifzz
c       ny       Second dimension of arrays a and adifzz
c       nz       Third dimension of arrays a and adifzz
c
c       a        Input array 
c
c       ibgn     Index in first dimension to begin operation
c       iend     Index in first dimension to end operation
c       jbgn     Index in second dimension to begin operation
c       jend     Index in second dimension to end operation
c       kbgn     Index in third dimension to begin operation
c       kend     Index in third dimension to end operation
c
c       dz       Grid spacing in the vertical direction (m)
c
c     OUTPUT:
c
c       adifzz   Differenced array del**2(a)/delz**2
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx, ny, nz         ! Number of grid points in 3 directions

      real a     (nx,ny,nz)      ! Input array
c
      integer ibgn,iend,jbgn,jend,kbgn,kend
c
      real dz
      real adifzz(nx,ny,nz)      ! Output difference array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,j,k
      real dzinv2
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c     Beginning of executable code...
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      dzinv2 = 1.0/(dz*dz)

      DO 100 k=kbgn,kend
      DO 100 j=jbgn,jend
      DO 100 i=ibgn,iend
c
c#######################################################################
c
c     NOTE:
c   
c     The order of calucation in the following formula should not
c     be changed. The order affects the machine trunction error
c     in the result.
c 
c#######################################################################
c
        adifzz(i,j,k)=((a(i,j,k+1)-a(i,j,k))-(a(i,j,k)-a(i,j,k-1)))
     :               *dzinv2

100   CONTINUE

      RETURN
      END
