c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BPDT                       ######
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 bpdt(nx,ny,nz,pprt,u,v,w,pdteb,pdtwb,pdtnb,pdtsb,
     :           pbar,rhobar,c2rho)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of pprt on the lateral boundaries
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/14/92.
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       pprt     Perturbation pressure at a given time level (Pascal)
c       u        x component of velocity at a given time level (m/s)
c       v        y component of velocity at a given time level (m/s)
c       w        Vertical component of velocity in Cartesian
c                coordinates at a given time level (m/s).
c
c     OUTPUT:
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     WORKING ARRAYS:
c
c       c2rho    rhobar*Cs**2, where Cs is sound wave speed . 
c                A local array.
c
c#######################################################################
c

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

      integer nx,ny,nz          ! The number of grid points in 3 directions

      real pprt  (nx,ny,nz)     ! Perturbation pressure (Pascal)
      real u     (nx,ny,nz)     ! Total u-velocity at tfuture (m/s).
      real v     (nx,ny,nz)     ! Total v-velocity at tfuture (m/s).
      real w     (nx,ny,nz)     ! Total w-velocity at tfuture (m/s).

      real pdteb (ny,nz)        ! Time tendency of pprt at east boundary
      real pdtwb (ny,nz)        ! Time tendency of pprt at west boundary
      real pdtnb (nx,nz)        ! Time tendency of pprt at north boundary
      real pdtsb (nx,nz)        ! Time tendency of pprt at south boundary

      real pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real c2rho (nx,ny,nz)     ! rhobar*Cs**2, where Cs is sound wave 
                                ! speed . A local array.

c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      real div                  ! Velocitry divergence defined as
                                ! div=du/dx+dv/dy+dw/dz.
 
      integer i,j,k
      real    dudx,dvdy,dwdz
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'phycst.inc'      ! Unchanging physical constants
      include 'globcst.inc'     ! Global constants that control model
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      IF(csopt.eq.1) THEN       ! Original definition of sound speed.

        do 1010 k=1,nz-1
        do 1010 j=1,ny-1
          c2rho( 1  ,j,k)= cpdcv*pbar( 1  ,j,k)
          c2rho(nx-1,j,k)= cpdcv*pbar(nx-1,j,k)
1010    continue

        do 1020 k=1,nz-1
        do 1020 i=1,nx-1
          c2rho(i, 1  ,k)= cpdcv*pbar(i, 1  ,k)
          c2rho(i,ny-1,k)= cpdcv*pbar(i,ny-1,k)
1020    continue

      ELSEIF(csopt.eq.2) THEN   ! Original speed multiplied by a factor

        do 1030 k=1,nz-1
        do 1030 j=1,ny-1
          c2rho( 1  ,j,k)= csfactr**2 *cpdcv*pbar( 1  ,j,k)
          c2rho(nx-1,j,k)= csfactr**2 *cpdcv*pbar(nx-1,j,k)
1030    continue

        do 1040 k=1,nz-1
        do 1040 i=1,nx-1
          c2rho(i, 1  ,k)= csfactr**2 *cpdcv*pbar(i, 1  ,k)
          c2rho(i,ny-1,k)= csfactr**2 *cpdcv*pbar(i,ny-1,k)
1040    continue

      ELSE                      ! Specified constant sound speed.

        do 1050 k=1,nz-1
        do 1050 j=1,ny-1
          c2rho( 1  ,j,k)= csound**2 * rhobar( 1  ,j,k)
          c2rho(nx-1,j,k)= csound**2 * rhobar(nx-1,j,k)
1050    continue

        do 1060 k=1,nz-1
        do 1060 i=1,nx-1
          c2rho(i, 1  ,k)= csound**2 * rhobar(i, 1  ,k)
          c2rho(i,ny-1,k)= csound**2 * rhobar(i,ny-1,k)
1060    continue

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of pprt on the west boundary (pdtwb)
c
c#######################################################################
c
      do 1110 k=1,nz-1
      do 1110 j=1,ny-1
        dudx=(u( 2 ,j,k)-u(1,j,k))*dxinv
        dvdy=(v(1,j+1,k)-v(1,j,k))*dyinv
        dwdz=(w(1,j,k+1)-w(1,j,k))*dzinv
        div=dudx+dvdy+dwdz

        pdtwb(j,k)=c2rho(1,j,k)*div

1110  continue

c
c#######################################################################
c
c     Calculate the time tendency of pprt on the east boundary (pdteb)
c
c#######################################################################
c
      do 1120 k=1,nz-1
      do 1120 j=1,ny-1
        dudx=(u( nx , j, k)-u(nx-1,j,k))*dxinv
        dvdy=(v(nx-1,j+1,k)-v(nx-1,j,k))*dyinv
        dwdz=(w(nx-1,j,k+1)-w(nx-1,j,k))*dzinv
        div=dudx+dvdy+dwdz

        pdteb(j,k)=c2rho(nx-1,j,k)*div

1120  continue

c
c#######################################################################
c
c     Calculate the time tendency of pprt on the south boundary (pdtsb)
c
c#######################################################################
c
      do 1310 k=1,nz-1
      do 1310 i=1,nx-1
        dudx=(u(i+1,1,k)-u(i,1,k))*dxinv
        dvdy=(v(i, 2 ,k)-v(i,1,k))*dyinv
        dwdz=(w(i,1,k+1)-w(i,1,k))*dzinv
        div=dudx+dvdy+dwdz

        pdtsb(i,k)=c2rho(i,1,k)*div

1310  continue

c
c#######################################################################
c
c     Calculate the time tendency of pprt on the north boundary (pdtnb)
c
c#######################################################################
c
      do 1410 k=1,nz-1
      do 1410 i=1,nx-1
        dudx=(u(i+1,ny-1,k)-u(i,ny-1,k))*dxinv
        dvdy=(v( i, ny , k)-v(i,ny-1,k))*dyinv
        dwdz=(w(i,ny-1,k+1)-w(i,ny-1,k))*dzinv
        div=dudx+dvdy+dwdz

        pdtnb(i,k)=c2rho(i,ny-1,k)*div

1410  continue

      RETURN
      END

