c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BUDIFXX                    ######
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 BUDIFXX(d2dxu, nx,ny,nz, ebc, wbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the east and west boundary values of del**2(u-ubar)/delx**2
c     which appear in the calculation of the 4th order computational 
c     mixing term of the u-equation.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c
c       d2dxu        A 3-D array whose east and west boundary values 
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       ebc, wbc     Parameter defining the type of boundary conditions 
c                    for array a.
c  
c                    The values 1,2,3,4,5 of ebc and wbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c  
c     OUTPUT:
c  
c       d2dxu        Array a with updated east and west boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dxu(nx,ny,nz)      ! Input array
      integer ebc, wbc          ! Control parameter of east and west b.c.
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer j, k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     West boundary
c
c#######################################################################
c 
      IF( wbc.eq.1) THEN         ! Wall boundary condition

        DO 100 k=2,nz-2
        DO 100 j=2,ny-2
          d2dxu(1,j,k) = -d2dxu(3,j,k)
100     CONTINUE

      ELSEIF( wbc.eq.2) THEN     ! Periodic boundary condition

        DO 110 k=2,nz-2
        DO 110 j=2,ny-2
          d2dxu(1,j,k) = d2dxu(nx-2,j,k)
110     CONTINUE

      ELSE                       ! All the other 3 types of conditions

        DO 120 k=2,nz-2
        DO 120 j=2,ny-2
          d2dxu(1,j,k) = d2dxu(2,j,k)
120     CONTINUE

      ENDIF

c
c#######################################################################
c
c     East boundary
c
c#######################################################################
c 

      IF( ebc.eq.1) THEN         ! Wall boundary condition
     
        DO 200 k=2,nz-2
        DO 200 j=2,ny-2
          d2dxu(nx,j,k) = -d2dxu(nx-2,j,k)
200     CONTINUE

      ELSEIF(ebc.eq.2) THEN      ! Periodic boundary 

        DO 210 k=2,nz-2
        DO 210 j=2,ny-2
          d2dxu(nx,j,k) = d2dxu(3,j,k)
210     CONTINUE

      ELSE                       ! All the other 3 types of conditions

        DO 220 k=2,nz-2
        DO 220 j=2,ny-2
          d2dxu(nx,j,k) = d2dxu(nx-1,j,k)
220     CONTINUE

      ENDIF

      RETURN
      END

c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BUDIFYY                    ######
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 BUDIFYY(d2dyu, nx,ny,nz, nbc, sbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the north and south boundary values of del**2(u-ubar)/dely**2
c     which appear in the calculation of the 4th order computational
c     mixing term of the u-equation.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c
c       d2dyu        A 3-D array whose north and south boundary values
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       nbc, sbc     Parameter defining the type of boundary conditions 
c                    for array a.
c  
c                    The values 1,2,3,4,5 of nbc and sbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c  
c     OUTPUT:
c  
c       d2dyu        Array a with updated north and south boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dyu(nx,ny,nz)      ! Input aray
      integer nbc, sbc          ! Control parameter for north/south b.c.    
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     South boundary
c
c#######################################################################
c 
      IF( sbc.eq.2) THEN         ! Periodic boundary

        DO 110 k=2,nz-2
        DO 110 i=2,nx-1
          d2dyu(i,1,k) = d2dyu(i,ny-2,k)
110     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 120 k=2,nz-2
        DO 120 i=2,nx-1
          d2dyu(i,1,k) = d2dyu(i,2,k)
120     CONTINUE

      ENDIF

c
c#######################################################################
c
c     North boundary
c
c#######################################################################
c 
      IF(nbc.eq.2) THEN          ! Periodic boundary

        DO 210 k=2,nz-2
        DO 210 i=2,nx-1
          d2dyu(i,ny-1,k) = d2dyu(i,2,k)
210     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 220 k=2,nz-2
        DO 220 i=2,nx-1
          d2dyu(i,ny-1,k) = d2dyu(i,ny-2,k)
220     CONTINUE

      ENDIF

      RETURN
      END 

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BUDIFZZ                    ######
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 BUDIFZZ(d2dzu, nx,ny,nz, tbc, bbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the top and bottom boundary values of del**2(u-ubar)/delz**2
c     which appear in the calculation of the 4th order computational 
c     mixing term of the u-equation.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c
c       d2dzu        A 3-D array whose top and bottom boundary values 
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       tbc, bbc     Parameter defining the type of boundary conditions 
c                    for array a.
c  
c                    The values 1,2,3,4,5 of tbc and bbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c  
c     OUTPUT:
c  
c       d2dzu        Array a with updated top and bottom boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dzu(nx,ny,nz)      ! Input aray
      integer tbc, bbc          ! Control parameter for top/bottom b.c.    
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, j
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Top boundary: 
c
c#######################################################################
c

      IF(tbc.eq.2) THEN          ! Periodic top and bottom condition

        DO 110 j=2,ny-2
        DO 110 i=2,nx-1
          d2dzu(i,j,nz-1) = d2dzu(i,j,2)
110     CONTINUE

      ELSE                       ! Rigid lid and zero normal gradient 

        DO 120 j=2,ny-2
        DO 120 i=2,nx-1
          d2dzu(i,j,nz-1) = d2dzu(i,j,nz-2)
120     CONTINUE

      ENDIF
c
c#######################################################################
c
c     Bottom boundary: 
c
c#######################################################################
c
      IF(bbc.eq.2) THEN          ! Periodic top and bottom condition

        DO 210 j=2,ny-2
        DO 210 i=2,nx-1
          d2dzu(i,j,1) = d2dzu(i,j,nz-2)
210     CONTINUE

      ELSE                       ! Solid ground and zero normal gradient 

        DO 220 j=2,ny-2
        DO 220 i=2,nx-1
          d2dzu(i,j,1) = d2dzu(i,j,2)
220     CONTINUE
 
      ENDIF

      RETURN
      END 

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BVDIFXX                    ######
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 BVDIFXX(d2dvx,nx,ny,nz, ebc,wbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the east and west boundary values of del**2(v-vbar)/delx**2
c     which appear in the calculation of the 4th order computational 
c     mixing term of the v-equation.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c  
c       d2dvx        A 3-D array whose east and west boundary values 
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       ebc, wbc     Parameter defining the type of boundary conditions 
c                    for array a.
c  
c                    The values 1,2,3,4,5 of ebc and wbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c  
c     OUTPUT:
c  
c       d2dvx        Array a with updated east and west boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dvx(nx,ny,nz)      ! Input aray
      integer ebc, wbc          ! Control parameter for east/west b.c.    
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer j, k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     West boundary
c
c#######################################################################
c 
      IF(wbc.eq.2) THEN          ! Periodic boundary 

        DO 110 k=2,nz-2
        DO 110 j=2,ny-1
          d2dvx(1,j,k) = d2dvx(nx-2,j,k)
110     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 120 k=2,nz-2
        DO 120 j=2,ny-1
          d2dvx(1,j,k) = d2dvx(2,j,k)
120     CONTINUE

      ENDIF

c
c#######################################################################
c
c     East boundary
c
c#######################################################################
c 
      IF(ebc.eq.2) THEN          ! Periodic boundary 

        DO 210 k=2,nz-2
        DO 210 j=2,ny-1
          d2dvx(nx-1,j,k) = d2dvx(2,j,k)
210     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 220 k=2,nz-2
        DO 220 j=2,ny-1
          d2dvx(nx-1,j,k) = d2dvx(nx-2,j,k)
220     CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BVDIFYY                    ######
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 BVDIFYY(d2dyv,nx,ny,nz, nbc,sbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the north and south boundary values of del**2(v-vbar)/dely**2
c     which appear in the calculation of the 4th order computational 
c     mixing term of the v-equation.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c
c       d2dyv        A 3-D array whose north and south boundary values 
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       nbc, sbc     Parameter defining the type of boundary conditions 
c                    for array a.
c  
c                    The values 1,2,3,4,5 of nbc and sbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c  
c     OUTPUT:
c  
c       d2dyv        Array a with updated north and south boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dyv(nx,ny,nz)      ! Input aray
      integer nbc, sbc          ! Control parameter for north/south b.c.    
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     South boundary
c
c#######################################################################
c 

      IF(sbc.eq.1) THEN 

        DO 100 k=2,nz-2
        DO 100 i=2,nx-2
          d2dyv(i,1,k) = -d2dyv(i,3,k)
100     CONTINUE

      ELSEIF(sbc.eq.2) THEN      ! Periodic boundary 

        DO 110 k=2,nz-2
        DO 110 i=2,nx-2
          d2dyv(i,1,k) = d2dyv(i,ny-2,k)
110     CONTINUE

      ELSE                       ! All the other 3 types of conditions

        DO 120 k=2,nz-2
        DO 120 i=2,nx-2
          d2dyv(i,1,k) = d2dyv(i,2,k)
120     CONTINUE

      ENDIF

c
c#######################################################################
c
c     North boundary
c
c#######################################################################
c 
      IF(nbc.eq.1) THEN 

        DO 200 k=2,nz-2
        DO 200 i=2,nx-2
          d2dyv(i,ny,k) = -d2dyv(i,ny-2,k)
200     CONTINUE

      ELSEIF(nbc.eq.2) THEN      ! Periodic boundary 

        DO 210 k=2,nz-2
        DO 210 i=2,nx-2
          d2dyv(i,ny,k) = d2dyv(i,3,k)
210     CONTINUE

      ELSE                       ! All the other 3 types of conditions

        DO 220 k=2,nz-2
        DO 220 i=2,nx-2
          d2dyv(i,ny,k) = d2dyv(i,ny-1,k)
220     CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BVDIFZZ                    ######
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 BVDIFZZ(d2dzv,nx,ny,nz, tbc,bbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the top and bottom boundary values of del**2(v-vbar)/delz**2
c     which appear in the calculation of the 4th order computational 
c     mixing term of the v-equation.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c
c       d2dzv        A 3-D array whose top and bottom boundary values 
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       tbc, bbc     Parameter defining the type of boundary conditions 
c                    for array a.
c
c                    The values 1,2,3,4,5 of tbc and bbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c
c     OUTPUT:
c
c       d2dzv        Array a with updated top and bottom boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dzv(nx,ny,nz)      ! Input aray
      integer tbc, bbc          ! Control parameter for top/bottom b.c.    
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Top boundary:
c
c#######################################################################
c

      IF(tbc.eq.2) THEN          ! Periodic top and bottom condition

        DO 110 j=2,ny-1
        DO 110 i=2,nx-2
          d2dzv(i,j,nz-1) = d2dzv(i,j,2)
110     CONTINUE

      ELSE                       ! Rigid lid and zero normal gradient

        DO 120 j=2,ny-1
        DO 120 i=2,nx-2
          d2dzv(i,j,nz-1) = d2dzv(i,j,nz-2)
120     CONTINUE

      ENDIF
c
c#######################################################################
c
c     Bottom boundary:
c
c#######################################################################
c
      IF(bbc.eq.2) THEN          ! Periodic top and bottom condition

        DO 210 j=2,ny-1
        DO 210 i=2,nx-2
          d2dzv(i,j,1) = d2dzv(i,j,nz-2)
210     CONTINUE

      ELSE                       ! Solid ground and zero normal gradient

        DO 220 j=2,ny-1
        DO 220 i=2,nx-2
          d2dzv(i,j,1) = d2dzv(i,j,2)
220     CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BWDIFXX                    ######
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 BWDIFXX(d2dxw,nx,ny,nz, ebc,wbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the east and west boundary values of del**2(w-wbar)/delx**2
c     which appear in the calculation of the 4th order computational 
c     mixing term of the w-equation.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c
c       d2dxw        A 3-D array whose east and west boundary values 
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       ebc, wbc     Parameter defining the type of boundary conditions 
c                    for array a.
c  
c                    The values 1,2,3,4,5 of ebc and wbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c  
c     OUTPUT:
c  
c       d2dxw        Array a with updated east and west boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dxw(nx,ny,nz)      ! Input aray
      integer ebc, wbc          ! Control parameter for east/west b.c.    
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer j, k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     West boundary
c
c#######################################################################
c 
      IF(wbc.eq.2) THEN          ! Periodic boundary 

        DO 110 k=2,nz-1
        DO 110 j=2,ny-2
          d2dxw(1,j,k) = d2dxw(nx-2,j,k)
110     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 120 k=2,nz-1
        DO 120 j=2,ny-2
          d2dxw(1,j,k) = d2dxw(2,j,k)
120     CONTINUE

      ENDIF

c
c#######################################################################
c
c     East boundary
c
c#######################################################################
c 
      IF(ebc.eq.2) THEN          ! Periodic boundary  

        DO 210 k=2,nz-1
        DO 210 j=2,ny-2
          d2dxw(nx-1,j,k) = d2dxw(2,j,k)
210     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 220 k=2,nz-1
        DO 220 j=2,ny-2
          d2dxw(nx-1,j,k) = d2dxw(nx-2,j,k)
220     CONTINUE

      ENDIF
      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BWDIFYY                    ######
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 BWDIFYY(d2dyw,nx,ny,nz, nbc,sbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the north and south boundary values of del**2(w-wbar)/dely**2
c     which appear in the calculation of the 4th order computational 
c     mixing term of the w-equation.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c
c       d2dyw        A 3-D array whose north and south boundary values 
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       nbc, sbc     Parameter defining the type of boundary conditions 
c                    for array a.
c  
c                    The values 1,2,3,4,5 of nbc and sbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c  
c     OUTPUT:
c  
c       d2dyw        Array a with updated north and south boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dyw(nx,ny,nz)      ! Input aray
      integer nbc, sbc          ! Control parameter for north/south b.c.    
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
c     South boundary
c
c#######################################################################
c 

      IF(sbc.eq.2) THEN          ! Periodic boundary

        DO 110 k=2,nz-1
        DO 110 i=2,nx-2
          d2dyw(i,1,k) = d2dyw(i,ny-2,k)
110     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 120 k=2,nz-1
        DO 120 i=2,nx-2
          d2dyw(i,1,k) = d2dyw(i,2,k)
120     CONTINUE

      ENDIF

c
c#######################################################################
c
c     North boundary
c
c#######################################################################
c 
      IF(nbc.eq.2) THEN          ! Periodic boundary

        DO 210 k=2,nz-1
        DO 210 i=2,nx-2
          d2dyw(i,ny-1,k) = d2dyw(i,2,k)
210     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 220 k=2,nz-1
        DO 220 i=2,nx-2
          d2dyw(i,ny-1,k) = d2dyw(i,ny-2,k)
220     CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BWDIFZZ                    ######
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 BWDIFZZ(d2dzw,nx,ny,nz, tbc,bbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the top and bottom boundary values of del**2(w-wbar)/delz**2
c     which appear in the calculation of the 4th order computational 
c     mixing term of the w-equation.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c  
c       d2dzw        A 3-D array whose top and bottom boundary values 
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       tbc, bbc     Parameter defining the type of boundary conditions 
c                    for array a.
c  
c                    The values 1,2,3,4,5 of tbc and bbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c  
c     OUTPUT:
c  
c       d2dzw        Array a with updated top and bottom boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dzw(nx,ny,nz)      ! Input aray
      integer tbc, bbc          ! Control parameter for top/bottom b.c.    
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Top boundary:
c
c#######################################################################
c

      IF(tbc.eq.2) THEN          ! Rigid lid at the top. mirror condition 

        DO 100 j=2,ny-2
        DO 100 i=2,nx-2
          d2dzw(i,j,nz) = -d2dzw(i,j,nz-2)
100     CONTINUE

      ELSEIF(tbc.eq.2) THEN      ! Periodic top and bottom condition

        DO 110 j=2,ny-2
        DO 110 i=2,nx-2
          d2dzw(i,j,nz) = d2dzw(i,j,3)
110     CONTINUE

      ELSE                       ! Zero normal gradient

        DO 120 j=2,ny-2
        DO 120 i=2,nx-2
          d2dzw(i,j,nz) = d2dzw(i,j,nz-2)
120     CONTINUE

      ENDIF
c
c#######################################################################
c
c     Bottom boundary:
c
c#######################################################################
c
      IF(bbc.eq.2) THEN          ! Solid ground, mirror reflection 

        DO 200 j=2,ny-2
        DO 200 i=2,nx-2
          d2dzw(i,j,1) = -d2dzw(i,j,3)
200     CONTINUE

      ELSEIF(bbc.eq.2) THEN      ! Periodic top and bottom condition

        DO 210 j=2,ny-2
        DO 210 i=2,nx-2
          d2dzw(i,j,1) = d2dzw(i,j,nz-2)
210     CONTINUE

      ELSE                       ! Zero normal gradient

        DO 220 j=2,ny-2
        DO 220 i=2,nx-2
          d2dzw(i,j,1) = d2dzw(i,j,3)
220     CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BSDIFXX                    ######
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 BSDIFXX(d2dxs,nx,ny,nz, ebc,wbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the east and west boundary values of del**2(scalar)/delx**2
c     which appear in the calculation of the 4th order computational 
c     mixing term of the scalar equations.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c  
c       d2dzs        A 3-D array whose east and west boundary values 
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       ebc, wbc     Parameter defining the type of boundary conditions 
c                    for array a.
c  
c                    The values 1,2,3,4,5 of ebc and wbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c  
c     OUTPUT:
c  
c       d2dzs        Array a with updated east and west boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dxs(nx,ny,nz)      ! Input aray
      integer ebc, wbc          ! Control parameter for east/west b.c.    
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer j, k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     West boundary
c
c#######################################################################
c 
      IF(wbc.eq.2) THEN          ! Periodic boundary

        DO 110 k=2,nz-2
        DO 110 j=2,ny-2
          d2dxs(1,j,k) = d2dxs(nx-2,j,k)
110     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 120 k=2,nz-2
        DO 120 j=2,ny-2
          d2dxs(1,j,k) = d2dxs(2,j,k)
120     CONTINUE

      ENDIF

c
c#######################################################################
c
c     East boundary
c
c#######################################################################
c 
      IF(ebc.eq.2) THEN          ! Periodic boundary

        DO 210 k=2,nz-2
        DO 210 j=2,ny-2
          d2dxs(nx-1,j,k) = d2dxs(2,j,k)
210     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 220 k=2,nz-2
        DO 220 j=2,ny-2
          d2dxs(nx-1,j,k) = d2dxs(nx-2,j,k)
220     CONTINUE

      ENDIF
      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BSDIFYY                    ######
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 BSDIFYY(d2dys,nx,ny,nz, nbc,sbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the north and south boundary values of del**2(scalar)/dely**2
c     which appear in the calculation of the 4th order computational 
c     mixing term of the scalar equations.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c
c       d2dys        A 3-D array whose north and south boundary values 
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       nbc, sbc     Parameter defining the type of boundary conditions 
c                    for array a.
c  
c                    The values 1,2,3,4,5 of nbc and sbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c  
c     OUTPUT:
c  
c       d2dys        Array a with updated north and south boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dys(nx,ny,nz)      ! Input aray
      integer nbc, sbc          ! Control parameter for north/south b.c.    
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     South boundary
c
c#######################################################################
c 
      IF(sbc.eq.2) THEN          ! Periodic boundary

        DO 110 k=2,nz-2
        DO 110 i=2,nx-2
          d2dys(i,1,k) = d2dys(i,ny-2,k)
110     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 120 k=2,nz-2
        DO 120 i=2,nx-2
          d2dys(i,1,k) = d2dys(i,2,k)
120     CONTINUE

      ENDIF

c
c#######################################################################
c
c     North boundary
c
c#######################################################################
c 
      IF(nbc.eq.2) THEN          ! Periodic boundary

        DO 210 k=2,nz-2
        DO 210 i=2,nx-2
          d2dys(i,ny-1,k) = d2dys(i,2,k)
210     CONTINUE

      ELSE                       ! All the other 4 types of conditions

        DO 220 k=2,nz-2
        DO 220 i=2,nx-2
          d2dys(i,ny-1,k) = d2dys(i,ny-2,k)
220     CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BSDIFZZ                    ######
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 BSDIFZZ(d2dzs,nx,ny,nz, tbc,bbc)

c
c#######################################################################
c
c     PURPOSE:
c
c     Set the top and bottom boundary values of del**2(scalar)/delz**2
c     which appear in the calculation of the 4th order computational 
c     mixing term of the scalar equations.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     1/10/1992.
c
c     MODIFICATION HISTORY:
c
c     6/10/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c
c       d2dzs        A 3-D array whose top and bottom boundary values 
c                    are set in this routine.
c       nx, ny, nz   Dimensions of array a
c       tbc, bbc     Parameter defining the type of boundary conditions 
c                    for array a.
c  
c                    The values 1,2,3,4,5 of tbc and bbc stands for 
c                    rigid wall, periodic, open, user specified and zero 
c                    normal gradient boundary condition respectively.
c  
c     OUTPUT:
c
c       d2dzs        Array a with updated top and bottom boundary values. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx,ny,nz          ! Number of grid points in x, y and z
                                ! directions
      real d2dzs(nx,ny,nz)      ! Input aray
      integer tbc,bbc           ! Control parameter for top/bottom b.c.    
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Top boundary:
c
c#######################################################################
c

      IF(tbc.eq.2) THEN          ! Periodic top and bottom condition

        DO 110 j=2,ny-2
        DO 110 i=2,nx-2
          d2dzs(i,j,nz-1) = d2dzs(i,j,2)
110     CONTINUE

      ELSE                       ! Rigid lid and zero normal gradient

        DO 120 j=2,ny-2
        DO 120 i=2,nx-2
          d2dzs(i,j,nz-1) = d2dzs(i,j,nz-2)
120     CONTINUE

      ENDIF
c
c#######################################################################
c
c     Bottom boundary:
c
c#######################################################################
c
      IF(bbc.eq.2) THEN          ! Periodic top and bottom condition

        DO 210 j=2,ny-2
        DO 210 i=2,nx-2
          d2dzs(i,j,1) = d2dzs(i,j,nz-2)
210     CONTINUE

      ELSE                       ! Solid ground and zero normal gradient

        DO 220 j=2,ny-2
        DO 220 i=2,nx-2
          d2dzs(i,j,1) = d2dzs(i,j,2)
220     CONTINUE

      ENDIF

      RETURN
      END

