c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BCU                        ######
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 BCU(nx,ny,nz,dtsml,
     :           u, udteb,udtwb,udtnb,udtsb,
     :           ebc,wbc,nbc,sbc,tbc,bbc)
c
C#######################################################################
c
c     PURPOSE:
c
c     Set the boundary conditions for the u-velocity component. Please note
c     that the values at the corner points may depend on the order that e-w,
c     n-s and t-b boundary conditions are applied.
c
C#######################################################################
c      
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/01/92 (M. Xue)
c     Added full documentation.
c
c     6/04/92 (M. Xue and H. Jin)
c     Further facelift. 
c
c     6/07/92 (M. Xue)
c     Modified to take in u at time tfuture only. 
c
c     6/15/92 (M. Xue and H. Jin)
c     Implemented open boundary condition 
c
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       dtsml    The small time step size (s)
c
c       u        Interior domain values of u-velocity at tfuture (m/s) 
c
c       udteb    Time tendency of the u field at the east boundary
c       udtwb    Time tendency of the u field at the west boundary
c       udtnb    Time tendency of the u field at the north boundary
c       udtsb    Time tendency of the u field at the south boundary
c
c       ebc      Parameter defining east   boundary condition type.
c       wbc      Parameter defining west   boundary condition type.
c       nbc      Parameter defining north  boundary condition type.
c       sbc      Parameter defining south  boundary condition type.
c       tbc      Parameter defining top    boundary condition type.
c       bbc      Parameter defining bottom boundary condition type.
c
c     OUTPUT:
c
c       u        The u-velocity over the entire domain at tfuture (m/s)
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 dtsml                ! The small time step size (s)
c
      real u  (nx,ny,nz)        ! Total u-velocity at tfuture (m/s)
c
      real udteb (ny,nz)        ! Time tendency of u at east boundary
      real udtwb (ny,nz)        ! Time tendency of u at west boundary
      real udtnb (nx,nz)        ! Time tendency of u at north boundary
      real udtsb (nx,nz)        ! Time tendency of u at south boundary
c
c#######################################################################
c
c     The following integer parameters define the type of condition
c     at each boundary.
c
c      1 for rigid wall (mirror) type boundary condition.
c      2 for periodic boundary condition.
c      3 for zero normal gradient boundary condition. 
c      4 for open (radiation) boundary condition.
c      5 for user (externally) specified boundary condition.
c
c#######################################################################
c
      integer ebc   ! Parameter defining east   boundary condition type.
      integer wbc   ! Parameter defining west   boundary condition type.
      integer nbc   ! Parameter defining north  boundary condition type.
      integer sbc   ! Parameter defining south  boundary condition type.
      integer tbc   ! Parameter defining top    boundary condition type.
      integer bbc   ! Parameter defining bottom boundary condition type.
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

c
C#######################################################################
c
c     Set the west boundary conditions
c
C#######################################################################
c      
      IF(wbc.eq.1) THEN             ! Wall boundary condition
      
        DO 1110 k=2,nz-2
        DO 1110 j=1,ny-1
          u(1,j,k)=-u(3,j,k)
1110    CONTINUE

      ELSEIF(wbc.eq.2) THEN         ! Periodic boundary condition.
        
        DO 1120 k=2,nz-2
        DO 1120 j=1,ny-1
          u(1,j,k)=u(nx-2,j,k)
1120    CONTINUE

      ELSEIF(wbc.eq.3) THEN         ! Zero normal gradient condition.
        
        DO 1130 k=2,nz-2
        DO 1130 j=1,ny-1
          u(1,j,k)=u(3,j,k)
1130    CONTINUE

      ELSEIF(wbc.eq.4 .or. wbc.eq.5) THEN 
                                    ! Radiation or user specified condition.
        
        DO 1140 k=2,nz-2
        DO 1140 j=2,ny-2
          u(1,j,k)=u(1,j,k)+udtwb(j,k)*dtsml
1140    CONTINUE

      ELSE

        write(6,900) 'BCU', wbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the east boundary conditions
c
c#######################################################################
c

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

      ELSEIF(ebc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1220 k=2,nz-2
        DO 1220 j=1,ny-1
          u(nx,j,k)=u(3,j,k)
1220    CONTINUE

      ELSEIF(ebc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1230 k=2,nz-2
        DO 1230 j=1,ny-1
          u(nx,j,k)=u(nx-2,j,k)
1230    CONTINUE

      ELSEIF(ebc.eq.4 .or. ebc.eq.5) THEN 
                                    ! Radiation or user specified condition.

        DO 1240 k=2,nz-2
        DO 1240 j=2,ny-2
          u(nx,j,k)=u(nx,j,k)+udteb(j,k)*dtsml
1240    CONTINUE

      ELSE

        write(6,900) 'BCU', ebc 
        STOP

      ENDIF
c
C#######################################################################
c
c     Set the north boundary conditions
c
C#######################################################################
c
      IF(nbc.eq.1) THEN             ! Wall boundary condition
      
        DO 1310 k=2,nz-2
        DO 1310 i=1,nx
          u(i,ny-1,k)=u(i,ny-2,k)
1310    CONTINUE

      ELSEIF(nbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1320 k=2,nz-2
        DO 1320 i=1,nx
          u(i,ny-1,k)=u(i,2,k)
1320    CONTINUE

      ELSEIF(nbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1330 k=2,nz-2
        DO 1330 i=1,nx
          u(i,ny-1,k)=u(i,ny-2,k)
1330    CONTINUE

      ELSEIF(nbc.eq.4 .or. nbc.eq.5) THEN 
                                    ! Radiation or user specified condition.

        DO 1340 k=2,nz-2
        DO 1340 i=2,nx-1
          u(i,ny-1,k)=u(i,ny-1,k)+udtnb(i,k)*dtsml
1340    CONTINUE

      ELSE

        write(6,900) 'BCU', nbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the south boundary conditions
c
c#######################################################################
c
      
      IF(sbc.eq.1) THEN             ! Wall boundary condition
      
        DO 1410 k=2,nz-2
        DO 1410 i=1,nx
          u(i,1,k)=u(i,2,k)
1410    CONTINUE

      ELSEIF(sbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1420 k=2,nz-2
        DO 1420 i=1,nx
          u(i,1,k)=u(i,ny-2,k)
1420    CONTINUE

      ELSEIF(sbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1430 k=2,nz-2
        DO 1430 i=1,nx
          u(i,1,k)=u(i,2,k)
1430    CONTINUE

      ELSEIF(sbc.eq.4 .or. sbc.eq.5) THEN
                                    ! Radiation or user specified condition.

        DO 1440 k=2,nz-2
        DO 1440 i=2,nx-1
          u(i,1,k)=u(i,1,k)+udtsb(i,k)*dtsml
1440    CONTINUE

      ELSE

        write(6,900) 'BCU', sbc 
        STOP

      ENDIF
c
C#######################################################################
c
c     Set the top boundary conditions
c
C#######################################################################
c
      IF(tbc.eq.1) THEN             ! Rigid lid boundary condition
      
        DO 1510 j=1,ny-1
        DO 1510 i=1,nx
          u(i,j,nz-1)=u(i,j,nz-2)
1510    CONTINUE

      ELSEIF(tbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1520 j=1,ny-1
        DO 1520 i=1,nx
          u(i,j,nz-1)=u(i,j,2)
1520    CONTINUE

      ELSEIF(tbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1530 j=1,ny-1
        DO 1530 i=1,nx
          u(i,j,nz-1)=u(i,j,nz-2)
1530    CONTINUE

      ELSE

        write(6,900) 'BCU', tbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the bottom boundary conditions
c
c#######################################################################
c

      IF(bbc.eq.1) THEN             ! Non-penetrative ground condition
      
        DO 1610 j=1,ny-1
        DO 1610 i=1,nx
          u(i,j,1)=u(i,j,2)
1610    CONTINUE

      ELSEIF(bbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1620 j=1,ny-1
        DO 1620 i=1,nx
          u(i,j,1)=u(i,j,nz-2)
1620    CONTINUE

      ELSEIF(bbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1630 j=1,ny-1
        DO 1630 i=1,nx
          u(i,j,1)=u(i,j,2)
1630    CONTINUE

      ELSE

        write(6,900) 'BCU', bbc 
        STOP

      ENDIF

c
c#######################################################################
c
c     IF the south and west boundaries are radiation conditions, the
c     southwest corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1710 k=2,nz-2
          u(1,1,k)=u(1,1,k)+udtsb(1,k)*dtsml
1710    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the south and east boundaries are radiation conditions, the
c     southeast corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1720 k=2,nz-2
          u(nx,1,k)=u(nx,1,k)+udtsb(nx,k)*dtsml
1720    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the north and west boundaries are radiation conditions, the 
c     northwest corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1730 k=2,nz-2
          u(1,ny-1,k)=u(1,ny-1,k)+udtnb(1,k)*dtsml
1730    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the north and east boundaries are radiation conditions, the
c     northeast corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1740 k=2,nz-2
          u(nx,ny-1,k)=u(nx,ny-1,k)+udtnb(nx,k)*dtsml
1740    CONTINUE

      ENDIF

      RETURN

900   format(1x,'Invalid boundary condition option found in ',a,
     :       /1x,'The option was ',i3,' Job stopped.')
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BCV                        ######
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 BCV(nx,ny,nz,dtsml,
     :           v, vdteb,vdtwb,vdtnb,vdtsb,
     :           ebc,wbc,nbc,sbc,tbc,bbc)
c
C#######################################################################
c
c     PURPOSE:
c
c     Set the boundary conditions for the v-velocity component. Please note 
c     that the values at the corner points may depend on the order that e-w, 
c     n-s and t-b boundary conditions are applied.
c
C#######################################################################
c      
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/01/92 (M. Xue)
c     Added full documentation.
c
c     6/04/92 (M. Xue and H. Jin)
c     Further facelift. 
c
c     6/07/92 (M. Xue)
c     Modified to take in v at time tfuture only. 
c
c     6/15/92 (M. Xue and H. Jin)
c     Implemented open boundary condition 
c
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       dtsml    The small time step size (s)
c
c       v        Interior domain values of v-velocity at tfuture (m/s) 
c
c       vdteb    Time tendency of the v field at the east boundary
c       vdtwb    Time tendency of the v field at the west boundary
c       vdtnb    Time tendency of the v field at the north boundary
c       vdtsb    Time tendency of the v field at the south boundary
c
c       ebc      Parameter defining east   boundary condition type.
c       wbc      Parameter defining west   boundary condition type.
c       nbc      Parameter defining north  boundary condition type.
c       sbc      Parameter defining south  boundary condition type.
c       tbc      Parameter defining top    boundary condition type.
c       bbc      Parameter defining bottom boundary condition type.
c
c     OUTPUT:
c
c       v        The v-velocity over the entire domain at tfuture (m/s)
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
c
      real dtsml                ! The small time step size (s)
c
      real v     (nx,ny,nz)     ! Total v-velocity at tfuture (m/s)
c
      real vdteb (ny,nz)        ! Time tendency of v at east boundary
      real vdtwb (ny,nz)        ! Time tendency of v at west boundary
      real vdtnb (nx,nz)        ! Time tendency of v at north boundary
      real vdtsb (nx,nz)        ! Time tendency of v at south boundary
c
c#######################################################################
c
c     The following integer parameters define the type of condition
c     at each boundary.
c
c      1 for rigid wall (mirror) type boundary condition.
c      2 for periodic boundary condition.
c      3 for zero normal gradient boundary condition. 
c      4 for open (radiation) boundary condition.
c      5 for user (externally) specified boundary condition.
c
c#######################################################################
c
      integer ebc   ! Parameter defining east   boundary condition type.
      integer wbc   ! Parameter defining west   boundary condition type.
      integer nbc   ! Parameter defining north  boundary condition type.
      integer sbc   ! Parameter defining south  boundary condition type.
      integer tbc   ! Parameter defining top    boundary condition type.
      integer bbc   ! Parameter defining bottom boundary condition type.
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

c
C#######################################################################
c
c     Set the west boundary conditions
c
C#######################################################################
c
      IF(wbc.eq.1) THEN             ! Wall boundary condition

        DO 1110 k=2,nz-2
        DO 1110 j=1,ny
          v(1,j,k)=v(2,j,k)
1110    CONTINUE

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

        DO 1120 k=2,nz-2
        DO 1120 j=1,ny
          v(1,j,k)=v(nx-2,j,k)
1120    CONTINUE

      ELSEIF(wbc.eq.3) THEN         ! Zero normal gradient condition.

        DO 1130 k=2,nz-2
        DO 1130 j=1,ny
          v(1,j,k)=v(2,j,k)
1130    CONTINUE

      ELSEIF(wbc.eq.4 .or. wbc.eq.5) THEN
                                    ! Radiation or user specified condition.

        DO 1140 k=2,nz-2
        DO 1140 j=2,ny-1
          v(1,j,k)=v(1,j,k)+vdtwb(j,k)*dtsml
1140    CONTINUE

      ELSE

        write(6,900) 'BCV', wbc
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the east boundary conditions
c
c#######################################################################
c

      IF(ebc.eq.1) THEN             ! Wall boundary condition

        DO 1210 k=2,nz-2
        DO 1210 j=1,ny
          v(nx-1,j,k)=v(nx-2,j,k)
1210    CONTINUE

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

        DO 1220 k=2,nz-2
        DO 1220 j=1,ny
          v(nx-1,j,k)=v(2,j,k)
1220    CONTINUE

      ELSEIF(ebc.eq.3) THEN         ! Zero normal gradient condition.

        DO 1230 k=2,nz-2
        DO 1230 j=1,ny
          v(nx-1,j,k)=v(nx-2,j,k)
1230    CONTINUE

      ELSEIF(ebc.eq.4 .or. ebc.eq.5) THEN
                                    ! Radiation or user specified condition.

        DO 1240 k=2,nz-2
        DO 1240 j=2,ny-1
          v(nx-1,j,k)=v(nx-1,j,k)+vdteb(j,k)*dtsml
1240    CONTINUE

      ELSE

        write(6,900) 'BCV', ebc
        STOP

      ENDIF
c
C#######################################################################
c
c     Set the north boundary conditions
c
C#######################################################################
c
      IF(nbc.eq.1) THEN             ! Wall boundary condition
      
        DO 1310 k=2,nz-2
        DO 1310 i=1,nx-1
          v(i,ny,k)=-v(i,ny-2,k)
1310    CONTINUE

      ELSEIF(nbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1320 k=2,nz-2
        DO 1320 i=1,nx-1
          v(i,ny,k)=v(i,3,k)
1320    CONTINUE

      ELSEIF(nbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1330 k=2,nz-2
        DO 1330 i=1,nx-1
          v(i,ny,k)=v(i,ny-2,k)
1330    CONTINUE

      ELSEIF(nbc.eq.4 .or. nbc.eq.5) THEN
                                    ! Radiation or user specified condition.

        DO 1340 k=2,nz-2
        DO 1340 i=2,nx-2
          v(i,ny,k)=v(i,ny,k)+vdtnb(i,k)*dtsml
1340    CONTINUE

      ELSE

        write(6,900) 'BCV', nbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the south boundary conditions
c
c#######################################################################
c
      
      IF(sbc.eq.1) THEN             ! Wall boundary condition
      
        DO 1410 k=2,nz-2
        DO 1410 i=1,nx-1
          v(i,1,k)=-v(i,3,k)
1410    CONTINUE

      ELSEIF(sbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1420 k=2,nz-2
        DO 1420 i=1,nx-1
          v(i,1,k)=v(i,ny-2,k)
1420    CONTINUE

      ELSEIF(sbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1430 k=2,nz-2
        DO 1430 i=1,nx-1
          v(i,1,k)=v(i,3,k)
1430    CONTINUE

      ELSEIF(sbc.eq.4 .or. sbc.eq.5) THEN
                                    ! Radiation or user specified condition.
        
        DO 1440 k=2,nz-2
        DO 1440 i=2,nx-2
          v(i,1,k)=v(i,1,k)+vdtsb(i,k)*dtsml
1440    CONTINUE

      ELSE

        write(6,900) 'BCV', sbc 
        STOP

      ENDIF
c
C#######################################################################
c
c     Set the top boundary conditions
c
C#######################################################################
c
      IF(tbc.eq.1) THEN             ! Rigid lid boundary condition
      
        DO 1510 j=1,ny
        DO 1510 i=1,nx-1
          v(i,j,nz-1)=v(i,j,nz-2)
1510    CONTINUE

      ELSEIF(tbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1520 j=1,ny
        DO 1520 i=1,nx-1
          v(i,j,nz-1)=v(i,j,2)
1520    CONTINUE

      ELSEIF(tbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1530 j=1,ny
        DO 1530 i=1,nx-1
          v(i,j,nz-1)=v(i,j,nz-2)
1530    CONTINUE

      ELSE

        write(6,900) 'BCV', tbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the bottom boundary conditions
c
c#######################################################################
c

      IF(bbc.eq.1) THEN             ! Non-penetrative ground condition
      
        DO 1610 j=1,ny
        DO 1610 i=1,nx-1
          v(i,j,1)=v(i,j,2)
1610    CONTINUE

      ELSEIF(bbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1620 j=1,ny
        DO 1620 i=1,nx-1
          v(i,j,1)=v(i,j,nz-2)
1620    CONTINUE

      ELSEIF(bbc.eq.3) THEN         ! Zero normal gradient condition.
        
        DO 1630 j=1,ny
        DO 1630 i=1,nx-1
          v(i,j,1)=v(i,j,2)
1630    CONTINUE

      ELSE

        write(6,900) 'BCV', bbc 
        STOP

      ENDIF

c
c#######################################################################
c
c     IF the south and west boundaries are radiation conditions, the
c     southwest corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1710 k=2,nz-2
          v(1,1,k)=v(1,1,k)+vdtsb(1,k)*dtsml
1710    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the south and east boundaries are radiation conditions, the
c     southeast corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1720 k=2,nz-2
          v(nx-1,1,k)=v(nx-1,1,k)+vdtsb(nx-1,k)*dtsml
1720    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the north and west boundaries are radiation conditions, the
c     northwest corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1730 k=2,nz-2
          v(1,ny,k)=v(1,ny,k)+vdtnb(1,k)*dtsml
1730    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the north and east boundaries are radiation conditions, the
c     northeast corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1740 k=2,nz-2
          v(nx-1,ny,k)=v(nx-1,ny,k)+vdtnb(nx-1,k)*dtsml
1740    CONTINUE

      ENDIF

      RETURN

900   format(1x,'Invalid boundary condition option found in ',a,
     :       /1x,'The option was ',i3,' Job stopped.')
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BCW                        ######
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 BCW(nx,ny,nz,dtsml,
     :           w, wdteb,wdtwb,wdtnb,wdtsb,
     :           ebc,wbc,nbc,sbc,tbc,bbc)
c
C#######################################################################
c
c     PURPOSE:
c
c     Set the boundary conditions for the w-velocity component. Please note 
c     that the values at the corner points may depend on the order that e-w, 
c     n-s and t-b boundary conditions are applied.
c
C#######################################################################
c      
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/01/92 (M. Xue)
c     Added full documentation.
c
c     6/04/92 (M. Xue and H. Jin)
c     Further facelift. 
c
c     6/07/92 (M. Xue)
c     Modified to take in w at time tfuture only. 
c
c     6/15/92 (M. Xue and H. Jin)
c     Implemented open boundary condition 
c
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       dtsml    The small time step size (s)
c
c       w        Interior domain values of w-velocity at tfuture (m/s) 
c
c       wdteb    Time tendency of the w field at the east boundary
c       wdtwb    Time tendency of the w field at the west boundary
c       wdtnb    Time tendency of the w field at the north boundary
c       wdtsb    Time tendency of the w field at the south boundary
c
c       ebc      Parameter defining east   boundary condition type.
c       wbc      Parameter defining west   boundary condition type.
c       nbc      Parameter defining north  boundary condition type.
c       sbc      Parameter defining south  boundary condition type.
c       tbc      Parameter defining top    boundary condition type.
c       bbc      Parameter defining bottom boundary condition type.
c
c     OUTPUT:
c
c       w        The w-velocity over the entire domain at tfuture (m/s)
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
c
      real dtsml                ! The small time step size (s)
c
      real w     (nx,ny,nz)     ! Total w-velocity at tfuture (m/s)
c
      real wdteb (ny,nz)        ! Time tendency of w at east boundary
      real wdtwb (ny,nz)        ! Time tendency of w at west boundary
      real wdtnb (nx,nz)        ! Time tendency of w at north boundary
      real wdtsb (nx,nz)        ! Time tendency of w at south boundary
c
c#######################################################################
c
c     The following integer parameters define the type of condition
c     at each boundary.
c
c      1 for rigid wall (mirror) type boundary condition.
c      2 for periodic boundary condition.
c      3 for zero normal gradient boundary condition. 
c      4 for open (radiation) boundary condition.
c      5 for user (externally) specified boundary condition.
c
c#######################################################################
c
      integer ebc   ! Parameter defining east   boundary condition type.
      integer wbc   ! Parameter defining west   boundary condition type.
      integer nbc   ! Parameter defining north  boundary condition type.
      integer sbc   ! Parameter defining south  boundary condition type.
      integer tbc   ! Parameter defining top    boundary condition type.
      integer bbc   ! Parameter defining bottom boundary condition type.
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

c
C#######################################################################
c
c     Set the west boundary conditions
c
C#######################################################################
c      
      IF(wbc.eq.1) THEN            ! Wall boundary condition
      
        DO 1110 k=1,nz 
        DO 1110 j=1,ny-1
          w(1,j,k)=w(2,j,k)
1110    CONTINUE

      ELSEIF(wbc.eq.2) THEN        ! Periodic boundary condition.
        
        DO 1120 k=1,nz 
        DO 1120 j=1,ny-1
          w(1,j,k)=w(nx-2,j,k)
1120    CONTINUE

      ELSEIF(wbc.eq.3) THEN        ! Zero normal gradient condition.
        
        DO 1130 k=1,nz 
        DO 1130 j=1,ny-1
          w(1,j,k)=w(2,j,k)
1130    CONTINUE

      ELSEIF(wbc.eq.4 .or. wbc.eq.5) THEN
                                    ! Radiation or user specified condition.
        
        DO 1140 k=1,nz 
        DO 1140 j=2,ny-2
          w(1,j,k)=w(1,j,k)+wdtwb(j,k)*dtsml
1140    CONTINUE

      ELSE

        write(6,900) 'BCW', wbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the east boundary conditions
c
c#######################################################################
c

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

      ELSEIF(ebc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1220 k=1,nz 
        DO 1220 j=1,ny-1
          w(nx-1,j,k)=w(2,j,k)
1220    CONTINUE

      ELSEIF(ebc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1230 k=1,nz 
        DO 1230 j=1,ny-1
          w(nx-1,j,k)=w(nx-2,j,k)
1230    CONTINUE

      ELSEIF(ebc.eq.4 .or. ebc.eq.5) THEN 
                                    ! Radiation or user specified condition.
        
        DO 1240 k=1,nz 
        DO 1240 j=2,ny-2
          w(nx-1,j,k)=w(nx-1,j,k)+wdteb(j,k)*dtsml
1240    CONTINUE

      ELSE

        write(6,900) 'BCW', ebc 
        STOP

      ENDIF
c
C#######################################################################
c
c     Set the north boundary conditions
c
C#######################################################################
c
      IF(nbc.eq.1) THEN             ! Wall boundary condition
      
        DO 1310 k=1,nz 
        DO 1310 i=1,nx-1
          w(i,ny-1,k)=w(i,ny-2,k)
1310    CONTINUE

      ELSEIF(nbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1320 k=1,nz 
        DO 1320 i=1,nx-1
          w(i,ny-1,k)=w(i,2,k)
1320    CONTINUE

      ELSEIF(nbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1330 k=1,nz 
        DO 1330 i=1,nx-1
          w(i,ny-1,k)=w(i,ny-2,k)
1330    CONTINUE

      ELSEIF(nbc.eq.4 .or. nbc.eq.5) THEN
                                    ! Radiation or user specified condition.
        
        DO 1340 k=1,nz 
        DO 1340 i=2,nx-2
          w(i,ny-1,k)=w(i,ny-1,k)+wdtnb(i,k)*dtsml
1340    CONTINUE

      ELSE

        write(6,900) 'BCW', nbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the south boundary conditions
c
c#######################################################################
c
      
      IF(sbc.eq.1) THEN             ! Wall boundary condition
      
        DO 1410 k=1,nz 
        DO 1410 i=1,nx-1
          w(i,1,k)=w(i,2,k)
1410    CONTINUE

      ELSEIF(sbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1420 k=1,nz 
        DO 1420 i=1,nx-1
          w(i,1,k)=w(i,ny-2,k)
1420    CONTINUE

      ELSEIF(sbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1430 k=1,nz 
        DO 1430 i=1,nx-1
          w(i,1,k)=w(i,2,k)
1430    CONTINUE

      ELSEIF(sbc.eq.4 .or. sbc.eq.5) THEN
                                    ! Radiation or user specified condition.
        
        DO 1440 k=1,nz 
        DO 1440 i=2,nx-2
          w(i,1,k)=w(i,1,k)+wdtsb(i,k)*dtsml
1440    CONTINUE

      ELSE

        write(6,900) 'BCW', sbc 
        STOP

      ENDif

c
c#######################################################################
c
c     Set the top boundary conditions
c
c#######################################################################
c
      IF(tbc.eq.1) THEN             ! Rigid lid boundary condition

        DO 1510 j=1,ny-1
        DO 1510 i=1,nx-1
          w(i,j,nz)=-w(i,j,nz-2)
          w(i,j,nz-1)=0.0
1510    CONTINUE

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

        DO 1520 j=1,ny-1
        DO 1520 i=1,nx-1
          w(i,j,nz)=w(i,j,3)
1520    CONTINUE

      ELSEIF(tbc.eq.3) THEN         ! Zero normal gradient condition.

        DO 1530 j=1,ny-1
        DO 1530 i=1,nx-1
          w(i,j,nz)=w(i,j,nz-1)
1530    CONTINUE

      ELSE

        write(6,900) 'BCW', tbc
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the bottom boundary conditions
c
c#######################################################################
c

      IF(bbc.eq.1) THEN             ! Non-penetrative ground condition

        DO 1610 j=1,ny-1
        DO 1610 i=1,nx-1
          w(i,j,1)=-w(i,j,3)
          w(i,j,2)=0.0
1610    CONTINUE

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

        DO 1620 j=1,ny-1
        DO 1620 i=1,nx-1
          w(i,j,1)=w(i,j,nz-2)
1620    CONTINUE

      ELSEIF(bbc.eq.3) THEN         ! Zero normal gradient condition.

        DO 1630 j=1,ny-1
        DO 1630 i=1,nx-1
          w(i,j,1)=w(i,j,2)
1630    CONTINUE

      ELSE

        write(6,900) 'BCW', bbc
        STOP

      ENDIF

c
c#######################################################################
c
c     IF the south and west boundaries are radiation conditions, the
c     southwest corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1710 k=1,nz
          w(1,1,k)=w(1,1,k)+wdtsb(1,k)*dtsml
1710    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the south and east boundaries are radiation conditions, the
c     southeast corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1720 k=1,nz
          w(nx-1,1,k)=w(nx-1,1,k)+wdtsb(nx-1,k)*dtsml
1720    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the north and west boundaries are radiation conditions, the
c     northwest corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1730 k=1,nz
          w(1,ny-1,k)=w(1,ny-1,k)+wdtnb(1,k)*dtsml
1730    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the north and east boundaries are radiation conditions, the
c     northeast corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1740 k=1,nz
          w(nx-1,ny-1,k)=w(nx-1,ny-1,k)+wdtnb(nx-1,k)*dtsml
1740    CONTINUE

      ENDIF

      RETURN

900   format(1x,'Invalid boundary condition option found in ',a,
     :       /1x,'The option was ',i3,' Job stopped.')
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BCP                        ######
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 BCP(nx,ny,nz,dtsml,
     :           pprt, pdteb,pdtwb,pdtnb,pdtsb,
     :           ebc,wbc,nbc,sbc,tbc,bbc)
c
c#######################################################################
c
c     PURPOSE:
c
c     Set the boundary conditions for the perturbation pressure, pprt.
c
c#######################################################################
c
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/01/92 (M. Xue)
c     Added full documentation.
c
c     6/07/92 (M. Xue)
c     Modified to take in pprt at time tfuture only.
c
c     6/15/92 (M. Xue and H. Jin)
c     Implemented open boundary condition 
c
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       dtsml    The small time step size (s)
c
c       pprt     Interior domain perturbation pressure at tfuture (Pascal)
c
c       pdteb    Time tendency of the pprt field at the east boundary
c       pdtwb    Time tendency of the pprt field at the west boundary
c       pdtnb    Time tendency of the pprt field at the north boundary
c       pdtsb    Time tendency of the pprt field at the south boundary
c
c       ebc      Parameter defining east   boundary condition type.
c       wbc      Parameter defining west   boundary condition type.
c       nbc      Parameter defining north  boundary condition type.
c       sbc      Parameter defining south  boundary condition type.
c       tbc      Parameter defining top    boundary condition type.
c       bbc      Parameter defining bottom boundary condition type.
c
c     OUTPUT:
c
c       pprt     Purturbation pressure over the entire domain at tfuture
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
c
      real dtsml                ! The small time step size (s)
c
      real pprt  (nx,ny,nz)     ! Perturbation presure at tfuture (Pascal) 
c
      real pdteb (ny,nz)        ! Time tendency of pprt field at east boundary
      real pdtwb (ny,nz)        ! Time tendency of pprt field at west boundary
      real pdtnb (nx,nz)        ! Time tendency of pprt field at north boundary
      real pdtsb (nx,nz)        ! Time tendency of pprt field at south boundary
c
c#######################################################################
c
c     The following integer parameters define the type of condition
c     at each boundary.
c
c      1 for rigid wall (mirror) type boundary condition.
c      2 for periodic boundary condition.
c      3 for zero normal gradient boundary condition. 
c      4 for open (radiation) boundary condition.
c      5 for user (externally) specified boundary condition.
c
c#######################################################################
c
      integer ebc   ! Parameter defining east   boundary condition type.
      integer wbc   ! Parameter defining west   boundary condition type.
      integer nbc   ! Parameter defining north  boundary condition type.
      integer sbc   ! Parameter defining south  boundary condition type.
      integer tbc   ! Parameter defining top    boundary condition type.
      integer bbc   ! Parameter defining bottom boundary condition type.
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

c
c#######################################################################
c
c     Set the west boundary conditions
c
c#######################################################################
c      
      IF(wbc.eq.1) THEN            ! Wall boundary condition
         
        DO 1110 k=2,nz-2
        DO 1110 j=1,ny-1
          pprt(1,j,k)=pprt(2,j,k)
1110    CONTINUE

      ELSEIF(wbc.eq.2) THEN        ! Periodic boundary condition.
      
        DO 1120 k=2,nz-2
        DO 1120 j=1,ny-1
          pprt(1,j,k)=pprt(nx-2,j,k)
1120    CONTINUE

      ELSEIF(wbc.eq.3) THEN        ! Zero normal gradient condition.
      
        DO 1130 k=2,nz-2
        DO 1130 j=1,ny-1
          pprt(1,j,k)=pprt(2,j,k)
1130    CONTINUE

      ELSEIF(wbc.eq.4 .or. wbc.eq.5) THEN
                                   ! Radiation or user specified condition.
      
        DO 1140 k=2,nz-2
        DO 1140 j=2,ny-2
          pprt(1,j,k)=pprt(1,j,k)+pdtwb(j,k)*dtsml
1140    CONTINUE

      ELSE

        write(6,900) 'BCP', wbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the east boundary conditions
c
c#######################################################################
c      
      IF(ebc.eq.1) THEN             ! Wall boundary condition
      
        DO 1210 k=2,nz-2
        DO 1210 j=1,ny-1
          pprt(nx-1,j,k)=pprt(nx-2,j,k)
1210    CONTINUE

      ELSEIF(ebc.eq.2) THEN         ! Periodic boundary condition.
      
        DO 1220 k=2,nz-2
        DO 1220 j=1,ny-1
          pprt(nx-1,j,k)=pprt(2,j,k)
1220    CONTINUE

      ELSEIF(ebc.eq.3) THEN         ! Zero normal gradient condition.
      
        DO 1230 k=2,nz-2
        DO 1230 j=1,ny-1
          pprt(nx-1,j,k)=pprt(nx-2,j,k)
1230    CONTINUE

      ELSEIF(ebc.eq.4 .or. ebc.eq.5) THEN
                                    ! Radiation or user specified condition.
      
        DO 1240 k=2,nz-2
        DO 1240 j=2,ny-2
          pprt(nx-1,j,k)=pprt(nx-1,j,k)+pdteb(j,k)*dtsml
1240    CONTINUE

      ELSE

        write(6,900) 'BCP', ebc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the north boundary conditions
c
c#######################################################################
c
      IF(nbc.eq.1) THEN             ! Wall boundary condition
      
        DO 1310 k=2,nz-2
        DO 1310 i=1,nx-1
          pprt(i,ny-1,k)=pprt(i,ny-2,k)
1310    CONTINUE

      ELSEIF(nbc.eq.2) THEN         ! Periodic boundary condition. 
      
        DO 1320 k=2,nz-2
        DO 1320 i=1,nx-1
          pprt(i,ny-1,k)=pprt(i,2,k)
1320    CONTINUE

      ELSEIF(nbc.eq.3) THEN         ! Zero normal gradient condition. 
      
        DO 1330 k=2,nz-2
        DO 1330 i=1,nx-1
          pprt(i,ny-1,k)=pprt(i,ny-2,k)
1330    CONTINUE

      ELSEIF(nbc.eq.4 .or. nbc.eq.5) THEN
                                    ! Radiation or user specified condition.

      
        DO 1340 k=2,nz-2
        DO 1340 i=2,nx-2
          pprt(i,ny-1,k)=pprt(i,ny-1,k)+pdtnb(i,k)*dtsml
1340    CONTINUE

      ELSE

        write(6,900) 'BCP', nbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the south boundary conditions
c
c#######################################################################
c
      IF(sbc.eq.1) THEN             ! Wall boundary condition
      
        DO 1410 k=2,nz-2
        DO 1410 i=1,nx-1
          pprt(i,1,k)=pprt(i,2,k)
1410    CONTINUE

      ELSEIF(sbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1420 k=2,nz-2
        DO 1420 i=1,nx-1
          pprt(i,1,k)=pprt(i,ny-2,k)
1420    CONTINUE

      ELSEIF(sbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1430 k=2,nz-2
        DO 1430 i=1,nx-1
          pprt(i,1,k)=pprt(i,2,k)
1430    CONTINUE

      ELSEIF(sbc.eq.4 .or. sbc.eq.5) THEN 
                                    ! Radiation or user specified condition.
        
        DO 1440 k=2,nz-2
        DO 1440 i=2,nx-2
          pprt(i,1,k)=pprt(i,1,k)+pdtsb(i,k)*dtsml
1440    CONTINUE

      ELSE

        write(6,900) 'BCP', sbc 
        STOP

      ENDIF
c
C#######################################################################
c
c     Set the top boundary conditions
c
C#######################################################################
c
      IF(tbc.eq.1) THEN             ! Rigid lid boundary condition
      
        DO 1510 j=1,ny-1
        DO 1510 i=1,nx-1
          pprt(i,j,nz-1)=pprt(i,j,nz-2)
1510    CONTINUE

      ELSEIF(tbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1520 j=1,ny-1
        DO 1520 i=1,nx-1
          pprt(i,j,nz-1)=pprt(i,j,2)
1520    CONTINUE

      ELSEIF(tbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1530 j=1,ny-1
        DO 1530 i=1,nx-1
          pprt(i,j,nz-1)=pprt(i,j,nz-2)
1530    CONTINUE

      ELSE

        write(6,900) 'BCP', tbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the bottom boundary conditions
c
c#######################################################################
c

      IF(bbc.eq.1) THEN             ! Non-penetrative ground condition
      
        DO 1610 j=1,ny-1
        DO 1610 i=1,nx-1
          pprt(i,j,1)=pprt(i,j,2)
1610    CONTINUE

      ELSEIF(bbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1620 j=1,ny-1
        DO 1620 i=1,nx-1
          pprt(i,j,1)=pprt(i,j,nz-2)
1620    CONTINUE

      ELSEIF(bbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1630 j=1,ny-1
        DO 1630 i=1,nx-1
          pprt(i,j,1)=pprt(i,j,2)
1630    CONTINUE

      ELSE

        write(6,900) 'BCP', bbc 
        STOP

      ENDIF
            
c
c#######################################################################
c
c     IF the south and west boundaries are radiation conditions, the
c     southwest corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1710 k=2,nz-2
          pprt(1,1,k)=pprt(1,1,k)+pdtsb(1,k)*dtsml
1710    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the south and east boundaries are radiation conditions, the
c     southeast corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1720 k=2,nz-2
          pprt(nx-1,1,k)=pprt(nx-1,1,k)+pdtsb(nx-1,k)*dtsml
1720    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the north and west boundaries are radiation conditions, the
c     northwest corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1730 k=2,nz-2
          pprt(1,ny-1,k)=pprt(1,ny-1,k)+pdtnb(1,k)*dtsml
1730    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the north and east boundaries are radiation conditions, the
c     northeast corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1740 k=2,nz-2
          pprt(nx-1,ny-1,k)=pprt(nx-1,ny-1,k)+pdtnb(nx-1,k)*dtsml
1740    CONTINUE

      ENDIF

      RETURN

900   format(1x,'Invalid boundary condition option found in ',a,
     :       /1x,'The option was ',i3,' Job stopped.')
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BCPT                       ######
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 BCPT(nx,ny,nz,dtbig,
     :           ptprt, ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :           ebc,wbc,nbc,sbc,tbc,bbc)
c
c#######################################################################
c
c     PURPOSE:
c
c     Set the boundary conditions for the potential temperature, ptprt.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/01/92 (M. Xue)
c     Added full documentation.
c
c     6/04/92 (M. Xue and H. Jin)
c     Further facelift. 
c
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       dtbig    The large time step size (s)
c
c       ptprt    Perturbation potential temperature at all time levels(K)
c
c       ptdteb   Time tendency of the ptprt field at the east boundary
c       ptdtwb   Time tendency of the ptprt field at the west boundary
c       ptdtnb   Time tendency of the ptprt field at the north boundary
c       ptdtsb   Time tendency of the ptprt field at the south boundary
c
c       ebc      Parameter defining east   boundary condition type.
c       wbc      Parameter defining west   boundary condition type.
c       nbc      Parameter defining north  boundary condition type.
c       sbc      Parameter defining south  boundary condition type.
c       tbc      Parameter defining top    boundary condition type.
c       bbc      Parameter defining bottom boundary condition type.
c
c     OUTPUT:
c
c       ptprt    ptprt over the entire domain at tfuture (K)
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
c
      real dtbig                ! The big time step size (s)
c
      integer nt                ! The no. of t-levels of t-dependent arrays.
      integer tpast             ! Index of time level for past time.
      integer tpresent          ! Index of time level for present time.
      integer tfuture           ! Index of time level for furture time.
c
      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)
c
      real ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
c
      real ptdteb(ny,nz)        ! Time tendency of ptprt at east boundary
      real ptdtwb(ny,nz)        ! Time tendency of ptprt at west boundary
      real ptdtnb(nx,nz)        ! Time tendency of ptprt at north boundary
      real ptdtsb(nx,nz)        ! Time tendency of ptprt at south boundary
c
c#######################################################################
c
c     The following integer parameters define the type of condition
c     at each boundary.
c
c      1 for rigid wall (mirror) type boundary condition.
c      2 for periodic boundary condition.
c      3 for zero normal gradient boundary condition. 
c      4 for open (radiation) boundary condition.
c      5 for user (externally) specified boundary condition.
c
c#######################################################################
c
      integer ebc   ! Parameter defining east   boundary condition type.
      integer wbc   ! Parameter defining west   boundary condition type.
      integer nbc   ! Parameter defining north  boundary condition type.
      integer sbc   ! Parameter defining south  boundary condition type.
      integer tbc   ! Parameter defining top    boundary condition type.
      integer bbc   ! Parameter defining bottom boundary condition type.
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external bcsclr
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
C#######################################################################
c
c     Set the boundary conditions for ptprt.
c
C#######################################################################
c      
      CALL bcsclr(nx,ny,nz,dtbig,
     :            ptprt(1,1,1,tpast),ptprt(1,1,1,tpresent),
     :            ptprt(1,1,1,tfuture),ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :            ebc,wbc,nbc,sbc,tbc,bbc)
            
      RETURN

900   format(1x,'Invalid boundary condition option found in ',a,
     :       /1x,'The option was ',i3,' Job stopped.')
      END
      
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BCQ                        ######
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 BCQ(nx,ny,nz,dtbig,
     :           q, qdteb,qdtwb,qdtnb,qdtsb,
     :           ebc,wbc,nbc,sbc,tbc,bbc)
c
c#######################################################################
c
c     PURPOSE:
c
c     Set the boundary conditions for all of the water quantities
c
c#######################################################################
c
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/01/92 (M. Xue)
c     Added full documentation.
c
c     6/04/92 (M. Xue and H. Jin)
c     Further facelift. 
c
C#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       dtbig    The large time step size (s)
c
c       q        mixing ratio for one of the water variables at all
c                time levels (kg/kg)
c
c       qdteb    Time tendency of the q field at the east boundary
c       qdtwb    Time tendency of the q field at the west boundary
c       qdtnb    Time tendency of the q field at the north boundary
c       qdtsb    Time tendency of the q field at the south boundary
c
c       ebc      Parameter defining east   boundary condition type.
c       wbc      Parameter defining west   boundary condition type.
c       nbc      Parameter defining north  boundary condition type.
c       sbc      Parameter defining south  boundary condition type.
c       tbc      Parameter defining top    boundary condition type.
c       bbc      Parameter defining bottom boundary condition type.
c
c     OUTPUT:
c
c       q        Array q over the entire domain at tfuture (m/s)
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
c
      real dtbig                ! The big time step size (s)
c
      integer nt                ! The no. of t-levels of t-dependent arrays.
      integer tpast             ! Index of time level for past time.
      integer tpresent          ! Index of time level for present time.
      integer tfuture           ! Index of time level for furture time.
c
      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)
c
      real q     (nx,ny,nz,nt)  ! mixing ratio for one of the water/ice
                                ! variables (kg/kg)
c
      real qdteb(ny,nz)         ! Time tendency of q at east boundary
      real qdtwb(ny,nz)         ! Time tendency of q at west boundary
      real qdtnb(nx,nz)         ! Time tendency of q at north boundary
      real qdtsb(nx,nz)         ! Time tendency of q at south boundary
c
c#######################################################################
c
c     The following integer parameters define the type of condition
c     at each boundary.
c
c      1 for rigid wall (mirror) type boundary condition.
c      2 for periodic boundary condition.
c      3 for zero normal gradient boundary condition. 
c      4 for open (radiation) boundary condition.
c      5 for user (externally) specified boundary condition.
c
c#######################################################################
c
      integer ebc   ! Parameter defining east   boundary condition type.
      integer wbc   ! Parameter defining west   boundary condition type.
      integer nbc   ! Parameter defining north  boundary condition type.
      integer sbc   ! Parameter defining south  boundary condition type.
      integer tbc   ! Parameter defining top    boundary condition type.
      integer bbc   ! Parameter defining bottom boundary condition type.
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external bcsclr
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
C#######################################################################
c
c     Set the boundary conditions for q.
c
C#######################################################################
c
      CALL bcsclr(nx,ny,nz,dtbig,
     :            q(1,1,1,tpast),q(1,1,1,tpresent),q(1,1,1,tfuture),
     :            qdteb,qdtwb,qdtnb,qdtsb,
     :            ebc,wbc,nbc,sbc,tbc,bbc)

      RETURN
      END

      
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BCSCLR                     ######
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 BCSCLR(nx,ny,nz,dtbig,
     :           s1,s2,s3,sdteb,sdtwb,sdtnb,sdtsb,
     :           ebc,wbc,nbc,sbc,tbc,bbc)
c
C#######################################################################
c
c     PURPOSE:
c
c     Set the boundary conditions for a scalar s. Please note that the values 
c     at the cornor points may depend on the order that e-w, n-s and t-b 
c     boundary conditions are applied.
c
C#######################################################################
c      
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/01/92 (M. Xue)
c     Added full documentation.
c
c     6/04/92 (M. Xue and H. Jin)
c     Further facelift. 
c
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       dtbig    The large time step size (s)
c
c       s1       A scalar variable at time tpast
c       s2       A scalar variable at time tpresent
c       s3       A scalar variable at time tfuture
c
c       sdteb    Time tendency of the s field at the east boundary
c       sdtwb    Time tendency of the s field at the west boundary
c       sdtnb    Time tendency of the s field at the north boundary
c       sdtsb    Time tendency of the s field at the south boundary
c
c       ebc      Parameter defining east   boundary condition type.
c       wbc      Parameter defining west   boundary condition type.
c       nbc      Parameter defining north  boundary condition type.
c       sbc      Parameter defining south  boundary condition type.
c       tbc      Parameter defining top    boundary condition type.
c       bbc      Parameter defining bottom boundary condition type.
c
c     OUTPUT:
c
c       s3       Scalar array s3 over the entire domain tfuture 
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
c
      real dtbig                ! The big time step size (s)
c
      real s1    (nx,ny,nz)     ! A scalar variable at time tpast.
      real s2    (nx,ny,nz)     ! A scalar variable at time tpresent.
      real s3    (nx,ny,nz)     ! A scalar variable at time tfuture.
c
      real sdteb(ny,nz)  ! Time tendency of s field at east boundary
      real sdtwb(ny,nz)  ! Time tendency of s field at west boundary
      real sdtnb(nx,nz)  ! Time tendency of s field at north boundary
      real sdtsb(nx,nz)  ! Time tendency of s field at south boundary
c
c#######################################################################
c
c     The following integer parameters define the type of condition
c     at each boundary.
c
c      1 for rigid wall (mirror) type boundary condition.
c      2 for periodic boundary condition.
c      3 for zero normal gradient boundary condition. 
c      4 for open (radiation) boundary condition.
c      5 for user (externally) specified boundary condition.
c
c#######################################################################
c
      integer ebc   ! Parameter defining east   boundary condition type.
      integer wbc   ! Parameter defining west   boundary condition type.
      integer nbc   ! Parameter defining north  boundary condition type.
      integer sbc   ! Parameter defining south  boundary condition type.
      integer tbc   ! Parameter defining top    boundary condition type.
      integer bbc   ! Parameter defining bottom boundary condition type.
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

c
c#######################################################################
c
c     Set the west boundary conditions
c
c#######################################################################
c      
      IF(wbc.eq.1) THEN            ! Wall boundary condition
         
        DO 1110 k=2,nz-2
        DO 1110 j=1,ny-1
          s3(1,j,k)=s3(2,j,k)
1110    CONTINUE

      ELSEIF(wbc.eq.2) THEN        ! Periodic boundary condition.
      
        DO 1120 k=2,nz-2
        DO 1120 j=1,ny-1
          s3(1,j,k)=s3(nx-2,j,k)
1120    CONTINUE

      ELSEIF(wbc.eq.3) THEN        ! Zero normal gradient condition.
      
        DO 1130 k=2,nz-2
        DO 1130 j=1,ny-1
          s3(1,j,k)=s3(2,j,k)
1130    CONTINUE

      ELSEIF(wbc.eq.4 .or. wbc.eq.5) THEN
                                    ! Radiation or user specified condition.
      
        DO 1140 k=2,nz-2
        DO 1140 j=2,ny-2
          s3(1,j,k)=s1(1,j,k)+sdtwb(j,k)*2.*dtbig
1140    CONTINUE

      ELSE

        write(6,900) 'BCSCLR', wbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the east boundary conditions
c
c#######################################################################
c      
      IF(ebc.eq.1) THEN             ! Wall boundary condition
      
        DO 1210 k=2,nz-2
        DO 1210 j=1,ny-1
          s3(nx-1,j,k)=s3(nx-2,j,k)
1210    CONTINUE

      ELSEIF(ebc.eq.2) THEN         ! Periodic boundary condition.
      
        DO 1220 k=2,nz-2
        DO 1220 j=1,ny-1
          s3(nx-1,j,k)=s3(2,j,k)
1220    CONTINUE

      ELSEIF(ebc.eq.3) THEN         ! Zero normal gradient condition.
      
        DO 1230 k=2,nz-2
        DO 1230 j=1,ny-1
          s3(nx-1,j,k)=s3(nx-2,j,k)
1230    CONTINUE

      ELSEIF(ebc.eq.4 .or. ebc.eq.5) THEN
                                    ! Radiation or user specified condition.
      
        DO 1240 k=2,nz-2
        DO 1240 j=2,ny-2
          s3(nx-1,j,k)=s1(nx-1,j,k)+sdteb(j,k)*2.*dtbig
1240    CONTINUE

      ELSE

        write(6,900) 'BCSCLR', ebc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the north boundary conditions
c
c#######################################################################
c
      IF(nbc.eq.1) THEN             ! Wall boundary condition
      
        DO 1310 k=2,nz-2
        DO 1310 i=1,nx-1
          s3(i,ny-1,k)=s3(i,ny-2,k)
1310    CONTINUE

      ELSEIF(nbc.eq.2) THEN         ! Periodic boundary condition. 
      
        DO 1320 k=2,nz-2
        DO 1320 i=1,nx-1
          s3(i,ny-1,k)=s3(i,2,k)
1320    CONTINUE

      ELSEIF(nbc.eq.3) THEN         ! Zero normal gradient condition. 
      
        DO 1330 k=2,nz-2
        DO 1330 i=1,nx-1
          s3(i,ny-1,k)=s3(i,ny-2,k)
1330    CONTINUE

      ELSEIF(nbc.eq.4 .or. nbc.eq.5) THEN
                                    ! Radiation or user specified condition.
      
        DO 1340 k=2,nz-2
        DO 1340 i=2,nx-2
          s3(i,ny-1,k)=s1(i,ny-1,k)+sdtnb(i,k)*2.*dtbig
1340    CONTINUE

      ELSE

        write(6,900) 'BCSCLR', nbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the south boundary conditions
c
c#######################################################################
c
      IF(sbc.eq.1) THEN             ! Wall boundary condition
      
        DO 1410 k=2,nz-2
        DO 1410 i=1,nx-1
          s3(i,1,k)=s3(i,2,k)
1410    CONTINUE

      ELSEIF(sbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1420 k=2,nz-2
        DO 1420 i=1,nx-1
          s3(i,1,k)=s3(i,ny-2,k)
1420    CONTINUE

      ELSEIF(sbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1430 k=2,nz-2
        DO 1430 i=1,nx-1
          s3(i,1,k)=s3(i,2,k)
1430    CONTINUE

      ELSEIF(sbc.eq.4 .or. sbc.eq.5) THEN
                                    ! Radiation or user specified condition.
        
        DO 1440 k=2,nz-2
        DO 1440 i=2,nx-2
          s3(i,1,k)=s1(i,1,k)+sdtsb(i,k)*2.*dtbig
1440    CONTINUE

      ELSE

        write(6,900) 'BCSCLR', sbc 
        STOP

      ENDIF
c
C#######################################################################
c
c     Set the top boundary conditions
c
C#######################################################################
c
      IF(tbc.eq.1) THEN             ! Rigid lid boundary condition
      
        DO 1510 j=1,ny-1
        DO 1510 i=1,nx-1
          s3(i,j,nz-1)=s3(i,j,nz-2)
1510    CONTINUE

      ELSEIF(tbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1520 j=1,ny-1
        DO 1520 i=1,nx-1
          s3(i,j,nz-1)=s3(i,j,2)
1520    CONTINUE

      ELSEIF(tbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1530 j=1,ny-1
        DO 1530 i=1,nx-1
          s3(i,j,nz-1)=s3(i,j,nz-2)
1530    CONTINUE

      ELSE

        write(6,900) 'BCSCLR', tbc 
        STOP

      ENDIF
c
c#######################################################################
c
c     Set the bottom boundary conditions
c
c#######################################################################
c

      IF(bbc.eq.1) THEN             ! Non-penetrative ground condition
      
        DO 1610 j=1,ny-1
        DO 1610 i=1,nx-1
          s3(i,j,1)=s3(i,j,2)
1610    CONTINUE

      ELSEIF(bbc.eq.2) THEN         ! Periodic boundary condition. 
        
        DO 1620 j=1,ny-1
        DO 1620 i=1,nx-1
          s3(i,j,1)=s3(i,j,nz-2)
1620    CONTINUE

      ELSEIF(bbc.eq.3) THEN         ! Zero normal gradient condition. 
        
        DO 1630 j=1,ny-1
        DO 1630 i=1,nx-1
          s3(i,j,1)=s3(i,j,2)
1630    CONTINUE

      ELSE

        write(6,900) 'BCSCLR', bbc 
        STOP

      ENDIF

c
c#######################################################################
c
c     IF the south and west boundaries are radiation conditions, the
c     southwest corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1710 k=2,nz-2
          s3(1,1,k)=s1(1,1,k)+sdtsb(1,k)*2.*dtbig
1710    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the south and east boundaries are radiation conditions, the
c     southeast corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1720 k=2,nz-2
          s3(nx-1,1,k)=s1(nx-1,1,k)+sdtsb(nx-1,k)*2.*dtbig
1720    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the north and west boundaries are radiation conditions, the
c     northwest corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1730 k=2,nz-2
          s3(1,ny-1,k)=s1(1,ny-1,k)+sdtnb(1,k)*2.*dtbig
1730    CONTINUE

      ENDIF

c
c#######################################################################
c
c     IF the north and east boundaries are radiation conditions, the
c     northeast corner should be determinded by radiation conditions.
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1740 k=2,nz-2
          s3(nx-1,ny-1,k)=s1(nx-1,ny-1,k)+sdtnb(nx-1,k)*2.*dtbig
1740    CONTINUE

      ENDIF

      RETURN

900   format(1x,'Invalid boundary condition option found in ',a,
     :       /1x,'The option was ',i3,' Job stopped.')
      END
