c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE INITBDT                    ######
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 INITBDT(nx,ny,nz,
     :           udteb , udtwb , udtnb , udtsb ,
     :           vdteb , vdtwb , vdtnb , vdtsb ,
     :           wdteb , wdtwb , wdtnb , wdtsb ,
     :           ptdteb, ptdtwb, ptdtnb, ptdtsb,
     :           pdteb , pdtwb , pdtnb , pdtsb ,
     :           qvdteb, qvdtwb, qvdtnb, qvdtsb,
     :           qcdteb, qcdtwb, qcdtnb, qcdtsb,
     :           qrdteb, qrdtwb, qrdtnb, qrdtsb,
     :           qidteb, qidtwb, qidtnb, qidtsb,
     :           qsdteb, qsdtwb, qsdtnb, qsdtsb,
     :           qhdteb, qhdtwb, qhdtnb, qhdtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Set the time tendencies of all time dependent variables on the
c     lateral boundaries for the initial time step. These tendencies are 
c     set to zero and stored in their respective 2-D arrays. they are used 
c     to set the boundary conditions of the corresponding variables for 
c     the first time step integration.
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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     OUTPUT:
c
c       udteb    Time tendency of u at east  boundary (m/s**2)
c       udtwb    Time tendency of u at west  boundary (m/s**2)
c       udtnb    Time tendency of u at north boundary (m/s**2)
c       udtsb    Time tendency of u at south boundary (m/s**2)
c
c       vdteb    Time tendency of v at east  boundary (m/s**2)
c       vdtwb    Time tendency of v at west  boundary (m/s**2)
c       vdtnb    Time tendency of v at north boundary (m/s**2)
c       vdtsb    Time tendency of v at south boundary (m/s**2)
c
c       wdteb    Time tendency of w at east  boundary (m/s**2)
c       wdtwb    Time tendency of w at west  boundary (m/s**2)
c       wdtnb    Time tendency of w at north boundary (m/s**2)
c       wdtsb    Time tendency of w at south boundary (m/s**2)
c
c       pdteb    Time tendency of pprt at east  boundary (Pescal/s)
c       pdtwb    Time tendency of pprt at west  boundary (Pescal/s)
c       pdtnb    Time tendency of pprt at north boundary (Pescal/s)
c       pdtsb    Time tendency of pprt at south boundary (Pescal/s)
c
c      ptdteb    Time tendency of ptprt at east  boundary (K/s)
c      ptdtwb    Time tendency of ptprt at west  boundary (K/s)
c      ptdtnb    Time tendency of ptprt at north boundary (K/s)
c      ptdtsb    Time tendency of ptprt at south boundary (K/s)
c
c      qvdteb    Time tendency of qv at east  boundary (kg/(kg*s))
c      qvdtwb    Time tendency of qv at west  boundary (kg/(kg*s))
c      qvdtnb    Time tendency of qv at north boundary (kg/(kg*s))
c      qvdtsb    Time tendency of qv at south boundary (kg/(kg*s))
c
c      qcdteb    Time tendency of qc at east  boundary (kg/(kg*s))
c      qcdtwb    Time tendency of qc at west  boundary (kg/(kg*s))
c      qcdtnb    Time tendency of qc at north boundary (kg/(kg*s))
c      qcdtsb    Time tendency of qc at south boundary (kg/(kg*s))
c
c      qrdteb    Time tendency of qr at east  boundary (kg/(kg*s))
c      qrdtwb    Time tendency of qr at west  boundary (kg/(kg*s))
c      qrdtnb    Time tendency of qr at north boundary (kg/(kg*s))
c      qrdtsb    Time tendency of qr at south boundary (kg/(kg*s))
c
c      qidteb    Time tendency of qi at east  boundary (kg/(kg*s))
c      qidtwb    Time tendency of qi at west  boundary (kg/(kg*s))
c      qidtnb    Time tendency of qi at north boundary (kg/(kg*s))
c      qidtsb    Time tendency of qi at south boundary (kg/(kg*s))
c
c      qsdteb    Time tendency of qs at east  boundary (kg/(kg*s))
c      qsdtwb    Time tendency of qs at west  boundary (kg/(kg*s))
c      qsdtnb    Time tendency of qs at north boundary (kg/(kg*s))
c      qsdtsb    Time tendency of qs at south boundary (kg/(kg*s))
c
c      qhdteb    Time tendency of qh at east  boundary (kg/(kg*s))
c      qhdtwb    Time tendency of qh at west  boundary (kg/(kg*s))
c      qhdtnb    Time tendency of qh at north boundary (kg/(kg*s))
c      qhdtsb    Time tendency of qh at south boundary (kg/(kg*s))
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 udteb (ny,nz)  ! T-tendency of u at e-boundary (m/s**2)
      real udtwb (ny,nz)  ! T-tendency of u at w-boundary (m/s**2)
      real udtnb (nx,nz)  ! T-tendency of u at n-boundary (m/s**2)
      real udtsb (nx,nz)  ! T-tendency of u at s-boundary (m/s**2)

      real vdteb (ny,nz)  ! T-tendency of v at e-boundary (m/s**2)
      real vdtwb (ny,nz)  ! T-tendency of v at w-boundary (m/s**2)
      real vdtnb (nx,nz)  ! T-tendency of v at n-boundary (m/s**2)
      real vdtsb (nx,nz)  ! T-tendency of v at s-boundary (m/s**2)

      real wdteb (ny,nz)  ! T-tendency of w at e-boundary (m/s**2)
      real wdtwb (ny,nz)  ! T-tendency of w at w-boundary (m/s**2)
      real wdtnb (nx,nz)  ! T-tendency of w at n-boundary (m/s**2)
      real wdtsb (nx,nz)  ! T-tendency of w at s-boundary (m/s**2)

      real ptdteb(ny,nz)  ! T-tendency of ptprt at e-boundary (K/s) 
      real ptdtwb(ny,nz)  ! T-tendency of ptprt at w-boundary (K/s) 
      real ptdtnb(nx,nz)  ! T-tendency of ptprt at n-boundary (K/s) 
      real ptdtsb(nx,nz)  ! T-tendency of ptprt at s-boundary (K/s) 

      real pdteb (ny,nz)  ! T-tendency of pprt at e-boundary (Pescal/s)
      real pdtwb (ny,nz)  ! T-tendency of pprt at w-boundary (Pescal/s)
      real pdtnb (nx,nz)  ! T-tendency of pprt at n-boundary (Pescal/s)
      real pdtsb (nx,nz)  ! T-tendency of pprt at s-boundary (Pescal/s)

      real qvdteb(ny,nz)  ! T-tendency of qv at e-boundary (kg/(kg*s))
      real qvdtwb(ny,nz)  ! T-tendency of qv at w-boundary (kg/(kg*s))
      real qvdtnb(nx,nz)  ! T-tendency of qv at n-boundary (kg/(kg*s))
      real qvdtsb(nx,nz)  ! T-tendency of qv at s-boundary (kg/(kg*s))

      real qcdteb(ny,nz)  ! T-tendency of qc at e-boundary (kg/(kg*s))
      real qcdtwb(ny,nz)  ! T-tendency of qc at w-boundary (kg/(kg*s))
      real qcdtnb(nx,nz)  ! T-tendency of qc at n-boundary (kg/(kg*s))
      real qcdtsb(nx,nz)  ! T-tendency of qc at s-boundary (kg/(kg*s))

      real qrdteb(ny,nz)  ! T-tendency of qr at e-boundary (kg/(kg*s))
      real qrdtwb(ny,nz)  ! T-tendency of qr at w-boundary (kg/(kg*s))
      real qrdtnb(nx,nz)  ! T-tendency of qr at n-boundary (kg/(kg*s))
      real qrdtsb(nx,nz)  ! T-tendency of qr at s-boundary (kg/(kg*s))

      real qidteb(ny,nz)  ! T-tendency of qi at e-boundary (kg/(kg*s))
      real qidtwb(ny,nz)  ! T-tendency of qi at w-boundary (kg/(kg*s))
      real qidtnb(nx,nz)  ! T-tendency of qi at n-boundary (kg/(kg*s))
      real qidtsb(nx,nz)  ! T-tendency of qi at s-boundary (kg/(kg*s))

      real qsdteb(ny,nz)  ! T-tendency of qs at e-boundary (kg/(kg*s))
      real qsdtwb(ny,nz)  ! T-tendency of qs at w-boundary (kg/(kg*s))
      real qsdtnb(nx,nz)  ! T-tendency of qs at n-boundary (kg/(kg*s))
      real qsdtsb(nx,nz)  ! T-tendency of qs at s-boundary (kg/(kg*s))

      real qhdteb(ny,nz)  ! T-tendency of qh at e-boundary (kg/(kg*s))
      real qhdtwb(ny,nz)  ! T-tendency of qh at w-boundary (kg/(kg*s))
      real qhdtnb(nx,nz)  ! T-tendency of qh at n-boundary (kg/(kg*s))
      real qhdtsb(nx,nz)  ! T-tendency of qh at s-boundary (kg/(kg*s))

c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external zerobdt
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
C#######################################################################
c
c     Set the time tendencies of u on the lateral boundaries
c     at initial time step
c
C#######################################################################
c

      CALL zerobdt(nx,ny,nz,udteb,udtwb,udtnb,udtsb)

c
C#######################################################################
c
c     Set the time tendencies of v on the lateral boundaries
c     at initial time step
c
C#######################################################################
c

      CALL zerobdt(nx,ny,nz,vdteb,vdtwb,vdtnb,vdtsb)

c
C#######################################################################
c
c     Set the time tendencies of w on the lateral boundaries
c     at initial time step
c
C#######################################################################
c

      CALL zerobdt(nx,ny,nz,wdteb,wdtwb,wdtnb,wdtsb)

c
C#######################################################################
c
c     Set the time tendencies of ptprt on the lateral boundaries
c     at initial time step
c
C#######################################################################
c

      CALL zerobdt(nx,ny,nz,ptdteb,ptdtwb,ptdtnb,ptdtsb)
c
C#######################################################################
c
c     Set the time tendencies of pprt on the lateral boundaries
c     at initial time step
c
C#######################################################################
c

      CALL zerobdt(nx,ny,nz,pdteb,pdtwb,pdtnb,pdtsb)

c
C#######################################################################
c
c     Set the time tendencies of qv on the lateral boundaries
c     at initial time step
c
C#######################################################################
c

      CALL zerobdt(nx,ny,nz,qvdteb,qvdtwb,qvdtnb,qvdtsb)

c
C#######################################################################
c
c     Set the time tendencies of qc on the lateral boundaries
c     at initial time step
c
C#######################################################################
c

      CALL zerobdt(nx,ny,nz,qcdteb,qcdtwb,qcdtnb,qcdtsb)

c
C#######################################################################
c
c     Set the time tendencies of qr on the lateral boundaries
c     at initial time step
c
C#######################################################################
c

      CALL zerobdt(nx,ny,nz,qrdteb,qrdtwb,qrdtnb,qrdtsb)

c
C#######################################################################
c
c     Set the time tendencies of qi on the lateral boundaries
c     at initial time step
c
C#######################################################################
c

      CALL zerobdt(nx,ny,nz,qidteb,qidtwb,qidtnb,qidtsb)

c
C#######################################################################
c
c     Set the time tendencies of qs on the lateral boundaries
c     at initial time step
c
C#######################################################################
c

      CALL zerobdt(nx,ny,nz,qsdteb,qsdtwb,qsdtnb,qsdtsb)

c
C#######################################################################
c
c     Set the time tendencies of qh on the lateral boundaries
c     at initial time step
c
C#######################################################################
c

      CALL zerobdt(nx,ny,nz,qhdteb,qhdtwb,qhdtnb,qhdtsb)

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ZEROBDT                    ######
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 ZEROBDT(nx,ny,nz,phidteb,phidtwb,phidtnb,phidtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Set the time tendencies of a time dependent variable on the lateral 
c     boundaries to zero for the initial time step.
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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     OUTPUT:
c
c       phidteb    Time tendency of phi at east  boundary 
c       phidtwb    Time tendency of phi at west  boundary 
c       phidtnb    Time tendency of phi at north boundary 
c       phidtsb    Time tendency of phi at south boundary 
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 phidteb (ny,nz)   ! T-tendency of phi at e-boundary 
      real phidtwb (ny,nz)   ! T-tendency of phi at w-boundary 
      real phidtnb (nx,nz)   ! T-tendency of phi at n-boundary 
      real phidtsb (nx,nz)   ! T-tendency of phi at s-boundary 
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 time tendency of phi on the west and east boundary to zero
c
c#######################################################################
c
      DO 1110 k=1,nz
      DO 1110 j=1,ny

        phidtwb(j,k)=0.0
        phidteb(j,k)=0.0

1110  CONTINUE

c
c#######################################################################
c
c     Set the time tendency of phi on the south and north boundary to zero
c
c#######################################################################
c
      DO 1120 k=1,nz
      DO 1120 i=1,nx

        phidtsb(i,k)=0.0
        phidtnb(i,k)=0.0

1120  CONTINUE

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE RADBDT                     ######
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 RADBDT(nx,ny,nz,dtbig1,
     :           u, v, w, ptprt, pprt, qv, qc, qr, qi, qs, qh,
     :           udteb , udtwb , udtnb , udtsb ,
     :           vdteb , vdtwb , vdtnb , vdtsb ,
     :           wdteb , wdtwb , wdtnb , wdtsb ,
     :           ptdteb, ptdtwb, ptdtnb, ptdtsb,
     :           pdteb , pdtwb , pdtnb , pdtsb ,
     :           qvdteb, qvdtwb, qvdtnb, qvdtsb,
     :           qcdteb, qcdtwb, qcdtnb, qcdtsb,
     :           qrdteb, qrdtwb, qrdtnb, qrdtsb,
     :           qidteb, qidtwb, qidtnb, qidtsb,
     :           qsdteb, qsdtwb, qsdtnb, qsdtsb,
     :           qhdteb, qhdtwb, qhdtnb, qhdtsb)
c 
c####################################################################### 
c 
c     PURPOSE:
c
c     Set the time tendencies of all time dependent variables on the lateral
c     boundaries according to Orlanski type radiation boundary conditions. 
c     These tendencies are stored in their respective 2-D arrays and will be 
c     used to set the boundary conditions of corresponding variables.
c     These tendencies are calculated near the end of one big time step
c     before the arrays are swapped in time. They are valid at the time level
c     "tfuture" before the swap and are carried into the next time step to 
c     calculate the boundary values of the respective variables at the new 
c     time level "tfuture".
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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       u        x-component of velocity at all time levels (m/s).
c       v        y-component of velocity at all time levels (m/s).
c       w        z-component of velocity at all time levels (m/s).
c       ptprt    Perturbation potential temperature at all time levels (K)
c       qv       Water vapor specific humidity at all time levels (kg/kg)
c       qc       Cloud water mixing ratio at all time levels (kg/kg)
c       qr       Rainwater mixing ratio at all time levels (kg/kg)
c       qi       Cloud ice mixing ratio at all time levels (kg/kg)
c       qs       Snow mixing ratio at all time levels (kg/kg)
c       qh       Hail mixing ratio at all time levels (kg/kg)
c
c     OUTPUT:
c
c       udteb    Time tendency of u at east  boundary (m/s**2)
c       udtwb    Time tendency of u at west  boundary (m/s**2)
c       udtnb    Time tendency of u at north boundary (m/s**2)
c       udtsb    Time tendency of u at south boundary (m/s**2)
c
c       vdteb    Time tendency of v at east  boundary (m/s**2)
c       vdtwb    Time tendency of v at west  boundary (m/s**2)
c       vdtnb    Time tendency of v at north boundary (m/s**2)
c       vdtsb    Time tendency of v at south boundary (m/s**2)
c
c       wdteb    Time tendency of w at east  boundary (m/s**2)
c       wdtwb    Time tendency of w at west  boundary (m/s**2)
c       wdtnb    Time tendency of w at north boundary (m/s**2)
c       wdtsb    Time tendency of w at south boundary (m/s**2)
c
c       pdteb    Time tendency of pprt at east  boundary (Pescal/s)
c       pdtwb    Time tendency of pprt at west  boundary (Pescal/s)
c       pdtnb    Time tendency of pprt at north boundary (Pescal/s)
c       pdtsb    Time tendency of pprt at south boundary (Pescal/s)
c
c      ptdteb    Time tendency of ptprt at east  boundary (K/s)
c      ptdtwb    Time tendency of ptprt at west  boundary (K/s)
c      ptdtnb    Time tendency of ptprt at north boundary (K/s)
c      ptdtsb    Time tendency of ptprt at south boundary (K/s)
c
c      qvdteb    Time tendency of qv at east  boundary (kg/(kg*s))
c      qvdtwb    Time tendency of qv at west  boundary (kg/(kg*s))
c      qvdtnb    Time tendency of qv at north boundary (kg/(kg*s))
c      qvdtsb    Time tendency of qv at south boundary (kg/(kg*s))
c
c      qcdteb    Time tendency of qc at east  boundary (kg/(kg*s))
c      qcdtwb    Time tendency of qc at west  boundary (kg/(kg*s))
c      qcdtnb    Time tendency of qc at north boundary (kg/(kg*s))
c      qcdtsb    Time tendency of qc at south boundary (kg/(kg*s))
c
c      qrdteb    Time tendency of qr at east  boundary (kg/(kg*s))
c      qrdtwb    Time tendency of qr at west  boundary (kg/(kg*s))
c      qrdtnb    Time tendency of qr at north boundary (kg/(kg*s))
c      qrdtsb    Time tendency of qr at south boundary (kg/(kg*s))
c
c      qidteb    Time tendency of qi at east  boundary (kg/(kg*s))
c      qidtwb    Time tendency of qi at west  boundary (kg/(kg*s))
c      qidtnb    Time tendency of qi at north boundary (kg/(kg*s))
c      qidtsb    Time tendency of qi at south boundary (kg/(kg*s))
c
c      qsdteb    Time tendency of qs at east  boundary (kg/(kg*s))
c      qsdtwb    Time tendency of qs at west  boundary (kg/(kg*s))
c      qsdtnb    Time tendency of qs at north boundary (kg/(kg*s))
c      qsdtsb    Time tendency of qs at south boundary (kg/(kg*s))
c
c      qhdteb    Time tendency of qh at east  boundary (kg/(kg*s))
c      qhdtwb    Time tendency of qh at west  boundary (kg/(kg*s))
c      qhdtnb    Time tendency of qh at north boundary (kg/(kg*s))
c      qhdtsb    Time tendency of qh at south boundary (kg/(kg*s))
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1               ! The big time step size (s)

      real u     (nx,ny,nz,nt)  ! Total u-velocity (m/s).
      real v     (nx,ny,nz,nt)  ! Total v-velocity (m/s).
      real w     (nx,ny,nz,nt)  ! Total w-velocity (m/s).
      real ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K).
      real pprt  (nx,ny,nz,nt)  ! Perturbation pressure (Pascal)

      real qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg).
      real qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg).
      real qr    (nx,ny,nz,nt)  ! Rain water mixing ratio (kg/kg).
      real qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (kg/kg).
      real qs    (nx,ny,nz,nt)  ! Snow mixing ratio (kg/kg).
      real qh    (nx,ny,nz,nt)  ! Hail mixing ratio (kg/kg).

      real udteb (ny,nz)  ! T-tendency of u at e-boundary (m/s**2)
      real udtwb (ny,nz)  ! T-tendency of u at w-boundary (m/s**2)
      real udtnb (nx,nz)  ! T-tendency of u at n-boundary (m/s**2)
      real udtsb (nx,nz)  ! T-tendency of u at s-boundary (m/s**2)

      real vdteb (ny,nz)  ! T-tendency of v at e-boundary (m/s**2)
      real vdtwb (ny,nz)  ! T-tendency of v at w-boundary (m/s**2)
      real vdtnb (nx,nz)  ! T-tendency of v at n-boundary (m/s**2)
      real vdtsb (nx,nz)  ! T-tendency of v at s-boundary (m/s**2)

      real wdteb (ny,nz)  ! T-tendency of w at e-boundary (m/s**2)
      real wdtwb (ny,nz)  ! T-tendency of w at w-boundary (m/s**2)
      real wdtnb (nx,nz)  ! T-tendency of w at n-boundary (m/s**2)
      real wdtsb (nx,nz)  ! T-tendency of w at s-boundary (m/s**2)

      real ptdteb(ny,nz)  ! T-tendency of ptprt at e-boundary (K/s)
      real ptdtwb(ny,nz)  ! T-tendency of ptprt at w-boundary (K/s)
      real ptdtnb(nx,nz)  ! T-tendency of ptprt at n-boundary (K/s)
      real ptdtsb(nx,nz)  ! T-tendency of ptprt at s-boundary (K/s)

      real pdteb (ny,nz)  ! T-tendency of pprt at e-boundary (Pescal/s)
      real pdtwb (ny,nz)  ! T-tendency of pprt at w-boundary (Pescal/s)
      real pdtnb (nx,nz)  ! T-tendency of pprt at n-boundary (Pescal/s)
      real pdtsb (nx,nz)  ! T-tendency of pprt at s-boundary (Pescal/s)

      real qvdteb(ny,nz)  ! T-tendency of qv at e-boundary (kg/(kg*s))
      real qvdtwb(ny,nz)  ! T-tendency of qv at w-boundary (kg/(kg*s))
      real qvdtnb(nx,nz)  ! T-tendency of qv at n-boundary (kg/(kg*s))
      real qvdtsb(nx,nz)  ! T-tendency of qv at s-boundary (kg/(kg*s))

      real qcdteb(ny,nz)  ! T-tendency of qc at e-boundary (kg/(kg*s))
      real qcdtwb(ny,nz)  ! T-tendency of qc at w-boundary (kg/(kg*s))
      real qcdtnb(nx,nz)  ! T-tendency of qc at n-boundary (kg/(kg*s))
      real qcdtsb(nx,nz)  ! T-tendency of qc at s-boundary (kg/(kg*s))

      real qrdteb(ny,nz)  ! T-tendency of qr at e-boundary (kg/(kg*s))
      real qrdtwb(ny,nz)  ! T-tendency of qr at w-boundary (kg/(kg*s))
      real qrdtnb(nx,nz)  ! T-tendency of qr at n-boundary (kg/(kg*s))
      real qrdtsb(nx,nz)  ! T-tendency of qr at s-boundary (kg/(kg*s))

      real qidteb(ny,nz)  ! T-tendency of qi at e-boundary (kg/(kg*s))
      real qidtwb(ny,nz)  ! T-tendency of qi at w-boundary (kg/(kg*s))
      real qidtnb(nx,nz)  ! T-tendency of qi at n-boundary (kg/(kg*s))
      real qidtsb(nx,nz)  ! T-tendency of qi at s-boundary (kg/(kg*s))

      real qsdteb(ny,nz)  ! T-tendency of qs at e-boundary (kg/(kg*s))
      real qsdtwb(ny,nz)  ! T-tendency of qs at w-boundary (kg/(kg*s))
      real qsdtnb(nx,nz)  ! T-tendency of qs at n-boundary (kg/(kg*s))
      real qsdtsb(nx,nz)  ! T-tendency of qs at s-boundary (kg/(kg*s))

      real qhdteb(ny,nz)  ! T-tendency of qh at e-boundary (kg/(kg*s))
      real qhdtwb(ny,nz)  ! T-tendency of qh at w-boundary (kg/(kg*s))
      real qhdtnb(nx,nz)  ! T-tendency of qh at n-boundary (kg/(kg*s))
      real qhdtsb(nx,nz)  ! T-tendency of qh at s-boundary (kg/(kg*s))

c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external bdtu1d
      external bdtv1d
      external bdtw1d
      external bdtpt1d
      external bdtp1d
      external bdtq1d
      external bdtu3d
      external bdtv3d
      external bdtw3d
      external bdtpt3d
      external bdtp3d
      external bdtq3d
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
C#######################################################################
c
c     Compute the time tendencies of u on the lateral boundaries 
c
C#######################################################################
c
      IF( rbcopt .eq.1 ) THEN         !

        CALL bdtu1d(nx,ny,nz,dtbig1,u,udteb,udtwb,udtnb,udtsb)

      ELSEIF ( rbcopt .eq.2 ) THEN

        CALL bdtu3d(nx,ny,nz,dtbig1,u,udteb,udtwb,udtnb,udtsb)

      ENDIF

c
C#######################################################################
c
c     Compute the time tendencies of v on the lateral boundaries 
c
C#######################################################################
c
      IF( rbcopt .eq.1 ) THEN         !

        CALL bdtv1d(nx,ny,nz,dtbig1,v,vdteb,vdtwb,vdtnb,vdtsb)

      ELSEIF ( rbcopt .eq.2 ) THEN

        CALL bdtv3d(nx,ny,nz,dtbig1,v,vdteb,vdtwb,vdtnb,vdtsb)

      ENDIF

c
C#######################################################################
c
c     Compute the time tendencies of w on the lateral boundaries 
c
C#######################################################################
c
      IF( rbcopt .eq.1 ) THEN         !

        CALL bdtw1d(nx,ny,nz,dtbig1,w,wdteb,wdtwb,wdtnb,wdtsb)

      ELSEIF ( rbcopt .eq.2 ) THEN

        CALL bdtw3d(nx,ny,nz,dtbig1,w,wdteb,wdtwb,wdtnb,wdtsb)

      ENDIF

c
C#######################################################################
c
c     Compute the time tendencies of ptprt on the lateral boundaries 
c
C#######################################################################
c
      IF( rbcopt .eq.1 ) THEN         !

        CALL bdtpt1d(nx,ny,nz,dtbig1,ptprt,ptdteb,ptdtwb,ptdtnb,ptdtsb)

      ELSEIF ( rbcopt .eq.2 ) THEN

        CALL bdtpt3d(nx,ny,nz,dtbig1,ptprt,ptdteb,ptdtwb,ptdtnb,ptdtsb)

      ENDIF

c
C#######################################################################
c
c     Compute the time tendencies of ptprt on the lateral boundaries 
c
C#######################################################################
c
      IF( rbcopt .eq.1 ) THEN         !

        CALL bdtp1d(nx,ny,nz,dtbig1,pprt,pdteb,pdtwb,pdtnb,pdtsb)

      ELSEIF ( rbcopt .eq.2 ) THEN

        CALL bdtp3d(nx,ny,nz,dtbig1,pprt,pdteb,pdtwb,pdtnb,pdtsb)

      ENDIF

c
C#######################################################################
c
c     Compute the time tendencies of qv on the lateral boundaries 
c
C#######################################################################
c
      IF( rbcopt .eq.1 ) THEN         !

        CALL bdtq1d(nx,ny,nz,dtbig1,qv,qvdteb,qvdtwb,qvdtnb,qvdtsb)

      ELSEIF ( rbcopt .eq.2 ) THEN

        CALL bdtq3d(nx,ny,nz,dtbig1,qv,qvdteb,qvdtwb,qvdtnb,qvdtsb)

      ENDIF

c
C#######################################################################
c
c     Compute the time tendencies of qc on the lateral boundaries 
c
C#######################################################################
c
      IF( rbcopt .eq.1 ) THEN         !

        CALL bdtq1d(nx,ny,nz,dtbig1,qc,qcdteb,qcdtwb,qcdtnb,qcdtsb)

      ELSEIF ( rbcopt .eq.2 ) THEN

        CALL bdtq3d(nx,ny,nz,dtbig1,qc,qcdteb,qcdtwb,qcdtnb,qcdtsb)

      ENDIF

c
C#######################################################################
c
c     Compute the time tendencies of qr on the lateral boundaries 
c
C#######################################################################
c
      IF( rbcopt .eq.1 ) THEN         !

        CALL bdtq1d(nx,ny,nz,dtbig1,qr,qrdteb,qrdtwb,qrdtnb,qrdtsb)

      ELSEIF ( rbcopt .eq.2 ) THEN

        CALL bdtq3d(nx,ny,nz,dtbig1,qr,qrdteb,qrdtwb,qrdtnb,qrdtsb)

      ENDIF

c
C#######################################################################
c
c     Compute the time tendencies of qi on the lateral boundaries 
c
C#######################################################################
c
      IF( rbcopt .eq.1 ) THEN         !

        CALL bdtq1d(nx,ny,nz,dtbig1,qi,qidteb,qidtwb,qidtnb,qidtsb)

      ELSEIF ( rbcopt .eq.2 ) THEN

        CALL bdtq3d(nx,ny,nz,dtbig1,qi,qidteb,qidtwb,qidtnb,qidtsb)

      ENDIF

c
C#######################################################################
c
c     Compute the time tendencies of qs on the lateral boundaries 
c
C#######################################################################
c
      IF( rbcopt .eq.1 ) THEN         !

        CALL bdtq1d(nx,ny,nz,dtbig1,qs,qsdteb,qsdtwb,qsdtnb,qsdtsb)

      ELSEIF ( rbcopt .eq.2 ) THEN

        CALL bdtq3d(nx,ny,nz,dtbig1,qs,qsdteb,qsdtwb,qsdtnb,qsdtsb)

      ENDIF

c
C#######################################################################
c
c     Compute the time tendencies of qh on the lateral boundaries 
c
C#######################################################################
c
      IF( rbcopt .eq.1 ) THEN         !

        CALL bdtq1d(nx,ny,nz,dtbig1,qh,qhdteb,qhdtwb,qhdtnb,qhdtsb)

      ELSEIF ( rbcopt .eq.2 ) THEN

        CALL bdtq3d(nx,ny,nz,dtbig1,qh,qhdteb,qhdtwb,qhdtnb,qhdtsb)

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTU1D                     ######
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 BDTU1D(nx,ny,nz,dtbig1, u,udteb,udtwb,udtnb,udtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of u on the lateral boundaries
c     using the radiation boundary condition of Orlanski.
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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       u        x-component of velocity at all time levels (m/s).
c
c     OUTPUT:
c
c       udteb    Time tendency of u at east  boundary (m/s**2)
c       udtwb    Time tendency of u at west  boundary (m/s**2)
c       udtnb    Time tendency of u at north boundary (m/s**2)
c       udtsb    Time tendency of u at south boundary (m/s**2)
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real u (nx,ny,nz,nt) ! Total u-velocity (m/s).

      real udteb (ny,nz)  ! T-tendency of u at e-boundary (m/s**2)
      real udtwb (ny,nz)  ! T-tendency of u at w-boundary (m/s**2)
      real udtnb (nx,nz)  ! T-tendency of u at n-boundary (m/s**2)
      real udtsb (nx,nz)  ! T-tendency of u at s-boundary (m/s**2)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of u on the west boundary (udtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=1,nz-1
        DO 1110 j=1,ny-1

          tem  =-u(2,j,k,3)-u(2,j,k,1)+2.*u(3,j,k,2)
          cdtdx=(u(2,j,k,1)-u(2,j,k,3))/(sign(eps,tem)+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)

          udtwb(j,k)=cdtdx*(-u(2,j,k,3)+u(1,j,k,2))/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of u on the east boundary (udteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdx must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=1,nz-1
        DO 1210 j=1,ny-1

          tem  =u(nx-1,j,k,3)+u(nx-1,j,k,1)-2.*u(nx-2,j,k,2)
          cdtdx=(u(nx-1,j,k,1)-u(nx-1,j,k,3))/(sign(eps,tem)+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)

          udteb(j,k)=cdtdx*(u(nx-1,j,k,3)-u(nx,j,k,2))/
     :               (dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of u on the south boundary (udtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1310 k=1,nz-1
        DO 1310 i=1,nx

          tem  =-u(i,2,k,3)-u(i,2,k,1)+2.*u(i,3,k,2)
          cdtdy=(u(i,2,k,1)-u(i,2,k,3))/(sign(eps,tem)+tem)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          udtsb(i,k)=cdtdy*(-u(i,2,k,3)+u(i,1,k,2))/
     :               (dtbig1*(1.-cdtdy))

1310    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of u on the north boundary (udtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=1,nz-1
        DO 1410 i=1,nx

          tem  =u(i,ny-2,k,3)+u(i,ny-2,k,1)-2.*u(i,ny-3,k,2)
          cdtdy=(u(i,ny-2,k,1)-u(i,ny-2,k,3))/(sign(eps,tem)+tem)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          udtnb(i,k)=cdtdy*(u(i,ny-2,k,3)-u(i,ny-1,k,2))/
     :               (dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c#######################################################################
c
c     Calculate the time tendency on the southwest corner of u boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=1,nz-1

          temx =(-u(2,2,k,3)-u(2,2,k,1)+2.*u(3,2,k,2))/(2.*dx)
          temy =(-u(2,2,k,3)-u(2,2,k,1)+2.*u(2,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(u(2,2,k,1)-u(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(u(2,2,k,1)-u(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          udtsb(1,k)=(cdtdx*(-u(2,2,k,3)+u(1,2,k,2))
     :               +cdtdy*(-u(2,2,k,3)+u(2,1,k,2)))/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          udtwb(1,k)=udtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of u boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=1,nz-1

          temx =(u(nx-1,2,k,3)+u(nx-1,2,k,1)-2.*u(nx-2,2,k,2))/(2.*dx)
          temy =(-u(nx-1,2,k,3)-u(nx-1,2,k,1)+2.*u(nx-1,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(u(nx-1,2,k,1)-u(nx-1,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(u(nx-1,2,k,1)-u(nx-1,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          udtsb(nx,k)=(cdtdx*(u(nx-1,2,k,3)-u(nx,2,k,2))
     :                +cdtdy*(-u(nx-1,2,k,3)+u(nx-1,1,k,2)))/
     :                (dtbig1*(1.+cdtdx-cdtdy))
          udteb( 1,k)=udtsb(nx,k)

1520    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of u boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=1,nz-1

          temx =(-u(2,ny-2,k,3)-u(2,ny-2,k,1)+2.*u(3,ny-2,k,2))/(2.*dx)
          temy =(u(2,ny-2,k,3)+u(2,ny-2,k,1)-2.*u(2,ny-3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(u(2,ny-2,k,1)-u(2,ny-2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(u(2,ny-2,k,1)-u(2,ny-2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          udtnb( 1  ,k)=(cdtdx*(-u(2,ny-2,k,3)+u(1,ny-2,k,2))
     :                  +cdtdy*(u(2,ny-2,k,3)-u(2,ny-1,k,2)))/
     :                  (dtbig1*(1.-cdtdx+cdtdy))
          udtwb(ny-1,k)=udtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of u boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=1,nz-1

          temx =(u(nx-1,ny-2,k,3)+u(nx-1,ny-2,k,1)
     :          -2.*u(nx-2,ny-2,k,2))/(2.*dx)
          temy =(u(nx-1,ny-2,k,3)+u(nx-1,ny-2,k,1)
     :          -2.*u(nx-1,ny-3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(u(nx-1,ny-2,k,1)-u(nx-1,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(u(nx-1,ny-2,k,1)-u(nx-1,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          udtnb(nx  ,k)=(cdtdx*(u(nx-1,ny-2,k,3)-u(nx,ny-2,k,2))
     :                  +cdtdy*(u(nx-1,ny-2,k,3)-u(nx-1,ny-1,k,2)))/
     :                  (dtbig1*(1.+cdtdx+cdtdy))
          udteb(ny-1,k)=udtnb(nx,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTV1D                     ######
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 BDTV1D(nx,ny,nz,dtbig1, v,vdteb,vdtwb,vdtnb,vdtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of v on the lateral boundaries
c     using the radiation boundary condition of Orlanski
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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       v        y-component of velocity at all time levels (m/s).
c
c     OUTPUT:
c
c       vdteb    Time tendency of v at east  boundary (m/s**2)
c       vdtwb    Time tendency of v at west  boundary (m/s**2)
c       vdtnb    Time tendency of v at north boundary (m/s**2)
c       vdtsb    Time tendency of v at south boundary (m/s**2)
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real v (nx,ny,nz,nt) ! Total v-velocity (m/s).

      real vdteb (ny,nz)  ! T-tendency of v at e-boundary (m/s**2)
      real vdtwb (ny,nz)  ! T-tendency of v at w-boundary (m/s**2)
      real vdtnb (nx,nz)  ! T-tendency of v at n-boundary (m/s**2)
      real vdtsb (nx,nz)  ! T-tendency of v at s-boundary (m/s**2)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of v on the west boundary (vdtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=1,nz-1
        DO 1110 j=1,ny

          tem  =-v(2,j,k,3)-v(2,j,k,1)+2.*v(3,j,k,2)
          cdtdx=(v(2,j,k,1)-v(2,j,k,3))/(sign(eps,tem)+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)

          vdtwb(j,k)=cdtdx*(-v(2,j,k,3)+v(1,j,k,2))/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of v on the east boundary (vdteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdx must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=1,nz-1
        DO 1210 j=1,ny

          tem  =v(nx-2,j,k,3)+v(nx-2,j,k,1)-2.*v(nx-3,j,k,2)
          cdtdx=(v(nx-2,j,k,1)-v(nx-2,j,k,3))/(sign(eps,tem)+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)

          vdteb(j,k)=cdtdx*(v(nx-2,j,k,3)-v(nx-1,j,k,2))/
     :               (dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of v on the south boundary (vdtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

      DO 1310 k=1,nz-1
      DO 1310 i=1,nx-1

        tem  =-v(i,2,k,3)-v(i,2,k,1)+2.*v(i,3,k,2)
        cdtdy=(v(i,2,k,1)-v(i,2,k,3))/(sign(eps,tem)+tem)
        cdtdy=min(cdtdy,0.0)
        cdtdy=max(cdtdy,-1.0)

        vdtsb(i,k)=cdtdy*(-v(i,2,k,3)+v(i,1,k,2))/
     :             (dtbig1*(1.-cdtdy))

1310  CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of v on the north boundary (vdtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=1,nz-1
        DO 1410 i=1,nx-1
          tem  =v(i,ny-1,k,3)+v(i,ny-1,k,1)-2.*v(i,ny-2,k,2)
          cdtdy=(v(i,ny-1,k,1)-v(i,ny-1,k,3))/(sign(eps,tem)+tem)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          vdtnb(i,k)=cdtdy*(v(i,ny-1,k,3)-v(i,ny,k,2))/
     :               (dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southwest corner of v boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=1,nz-1

          temx =(-v(2,2,k,3)-v(2,2,k,1)+2.*v(3,2,k,2))/(2.*dx)
          temy =(-v(2,2,k,3)-v(2,2,k,1)+2.*v(2,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(v(2,2,k,1)-v(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(v(2,2,k,1)-v(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          vdtsb(1,k)=(cdtdx*(-v(2,2,k,3)+v(1,2,k,2))
     :               +cdtdy*(-v(2,2,k,3)+v(2,1,k,2)))/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          vdtwb(1,k)=vdtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of v boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=1,nz-1

          temx =(v(nx-2,2,k,3)+v(nx-2,2,k,1)-2.*v(nx-3,2,k,2))/(2.*dx)
          temy =(-v(nx-2,2,k,3)-v(nx-2,2,k,1)+2.*v(nx-2,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(v(nx-2,2,k,1)-v(nx-2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(v(nx-2,2,k,1)-v(nx-2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          vdtsb(nx-1,k)=(cdtdx*(v(nx-2,2,k,3)-v(nx-1,2,k,2))
     :                  +cdtdy*(-v(nx-2,2,k,3)+v(nx-2,1,k,2)))/
     :                  (dtbig1*(1.+cdtdx-cdtdy))
          vdteb( 1  ,k)=vdtsb(nx-1,k)

1520    CONTINUE

      ENDIF
c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of v boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=1,nz-1

          temx =(-v(2,ny-1,k,3)-v(2,ny-1,k,1)+2.*v(3,ny-1,k,2))/(2.*dx)
          temy =(v(2,ny-1,k,3)+v(2,ny-1,k,1)-2.*v(2,ny-2,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(v(2,ny-1,k,1)-v(2,ny-1,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(v(2,ny-1,k,1)-v(2,ny-1,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          vdtnb( 1,k)=(cdtdx*(-v(2,ny-1,k,3)+v(1,ny-1,k,2))
     :                  +cdtdy*(v(2,ny-1,k,3)-v(2,ny,k,2)))/
     :                  (dtbig1*(1.-cdtdx+cdtdy))
          vdtwb(ny,k)=vdtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of v boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=1,nz-1

          temx =(v(nx-2,ny-1,k,3)+v(nx-2,ny-1,k,1)
     :          -2.*v(nx-3,ny-1,k,2))/(2.*dx)
          temy =(v(nx-2,ny-1,k,3)+v(nx-2,ny-1,k,1)
     :          -2.*v(nx-2,ny-2,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(v(nx-2,ny-1,k,1)-v(nx-2,ny-1,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(v(nx-2,ny-1,k,1)-v(nx-2,ny-1,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          vdtnb(nx-1,k)=(cdtdx*(v(nx-2,ny-1,k,3)-v(nx-1,ny-1,k,2))
     :                  +cdtdy*(v(nx-2,ny-1,k,3)-v(nx-2,ny,k,2)))/
     :                  (dtbig1*(1.+cdtdx+cdtdy))
          vdteb(ny  ,k)=vdtnb(nx-1,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTW1D                     ######
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 BDTW1D(nx,ny,nz,dtbig1, w,wdteb,wdtwb,wdtnb,wdtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of w on the lateral boundaries
c     using the radiation boundary condition of Orlanski.
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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       w        z-component of velocity at all time levels (m/s).
c
c     OUTPUT:
c
c       wdteb    Time tendency of w at east  boundary (m/s**2)
c       wdtwb    Time tendency of w at west  boundary (m/s**2)
c       wdtnb    Time tendency of w at north boundary (m/s**2)
c       wdtsb    Time tendency of w at south boundary (m/s**2)
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real w (nx,ny,nz,nt) ! Total w-velocity (m/s).

      real wdteb (ny,nz)  ! T-tendency of w at e-boundary (m/s**2)
      real wdtwb (ny,nz)  ! T-tendency of w at w-boundary (m/s**2)
      real wdtnb (nx,nz)  ! T-tendency of w at n-boundary (m/s**2)
      real wdtsb (nx,nz)  ! T-tendency of w at s-boundary (m/s**2)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of w on the west boundary (wdtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=1,nz
        DO 1110 j=1,ny-1

          tem  =-w(2,j,k,3)-w(2,j,k,1)+2.*w(3,j,k,2)
          cdtdx=(w(2,j,k,1)-w(2,j,k,3))/(sign(eps,tem)+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)

          wdtwb(j,k)=cdtdx*(-w(2,j,k,3)+w(1,j,k,2))/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of w on the east boundary (wdteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=1,nz
        DO 1210 j=1,ny-1

          tem  =w(nx-2,j,k,3)+w(nx-2,j,k,1)-2.*w(nx-3,j,k,2)
          cdtdx=(w(nx-2,j,k,1)-w(nx-2,j,k,3))/(sign(eps,tem)+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)

          wdteb(j,k)=cdtdx*(w(nx-2,j,k,3)-w(nx-1,j,k,2))/
     :               (dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of w on the south boundary (wdtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1310 k=1,nz
        DO 1310 i=1,nx-1

          tem  =-w(i,2,k,3)-w(i,2,k,1)+2.*w(i,3,k,2)
          cdtdy=(w(i,2,k,1)-w(i,2,k,3))/(sign(eps,tem)+tem)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          wdtsb(i,k)=cdtdy*(-w(i,2,k,3)+w(i,1,k,2))/
     :               (dtbig1*(1.-cdtdy))

1310    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of w on the north boundary (wdtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=1,nz
        DO 1410 i=1,nx-1

          tem  =w(i,ny-2,k,3)+w(i,ny-2,k,1)-2.*w(i,ny-3,k,2)
          cdtdy=(w(i,ny-2,k,1)-w(i,ny-2,k,3))/(sign(eps,tem)+tem)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          wdtnb(i,k)=cdtdy*(w(i,ny-2,k,3)-w(i,ny-1,k,2))/
     :               (dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southwest corner of w boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=1,nz

          temx =(-w(2,2,k,3)-w(2,2,k,1)+2.*w(3,2,k,2))/(2.*dx)
          temy =(-w(2,2,k,3)-w(2,2,k,1)+2.*w(2,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(w(2,2,k,1)-w(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(w(2,2,k,1)-w(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          wdtsb(1,k)=(cdtdx*(-w(2,2,k,3)+w(1,2,k,2))
     :               +cdtdy*(-w(2,2,k,3)+w(2,1,k,2)))/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          wdtwb(1,k)=wdtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of w boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=1,nz

          temx =(w(nx-2,2,k,3)+w(nx-2,2,k,1)-2.*w(nx-3,2,k,2))/(2.*dx)
          temy =(-w(nx-2,2,k,3)-w(nx-2,2,k,1)+2.*w(nx-2,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(w(nx-2,2,k,1)-w(nx-2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(w(nx-2,2,k,1)-w(nx-2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          wdtsb(nx-1,k)=(cdtdx*(w(nx-2,2,k,3)-w(nx-1,2,k,2))
     :                  +cdtdy*(-w(nx-2,2,k,3)+w(nx-2,1,k,2)))/
     :                  (dtbig1*(1.+cdtdx-cdtdy))
          wdteb( 1  ,k)=wdtsb(nx-1,k)

1520    CONTINUE

      ENDIF
c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of w boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=1,nz

          temx =(-w(2,ny-2,k,3)-w(2,ny-2,k,1)+2.*w(3,ny-2,k,2))/(2.*dx)
          temy =(w(2,ny-2,k,3)+w(2,ny-2,k,1)-2.*w(2,ny-3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(w(2,ny-2,k,1)-w(2,ny-2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(w(2,ny-2,k,1)-w(2,ny-2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          wdtnb( 1  ,k)=(cdtdx*(-w(2,ny-2,k,3)+w(1,ny-2,k,2))
     :                  +cdtdy*(w(2,ny-2,k,3)-w(2,ny-1,k,2)))/
     :                  (dtbig1*(1.-cdtdx+cdtdy))
          wdtwb(ny-1,k)=wdtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of w boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=1,nz

          temx =(w(nx-2,ny-2,k,3)+w(nx-2,ny-2,k,1)
     :          -2.*w(nx-3,ny-2,k,2))/(2.*dx)
          temy =(w(nx-2,ny-2,k,3)+w(nx-2,ny-2,k,1)
     :          -2.*w(nx-2,ny-3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(w(nx-2,ny-2,k,1)-w(nx-2,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(w(nx-2,ny-2,k,1)-w(nx-2,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          wdtnb(nx-1,k)=(cdtdx*(w(nx-2,ny-2,k,3)-w(nx-1,ny-2,k,2))
     :                  +cdtdy*(w(nx-2,ny-2,k,3)-w(nx-2,ny-1,k,2)))/
     :                  (dtbig1*(1.+cdtdx+cdtdy))
          wdteb(ny-1,k)=wdtnb(nx-1,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTP1D                     ######
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 BDTP1D(nx,ny,nz,dtbig1,pprt,pdteb,pdtwb,pdtnb,pdtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of the perturbation pressure on the lateral 
c     boundaries using the radiation boundary condition of Orlanski.
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 all time levels (Pascal)
c
c     OUTPUT:
c
c       pdteb    Time tendency of perturbation pressre at east  
c                boundary (Pescal/s)
c       pdtwb    Time tendency of perturbation pressre at west  
c                boundary (Pescal/s)
c       pdtnb    Time tendency of perturbation pressre at north 
c                boundary (Pescal/s)
c       pdtsb    Time tendency of perturbation pressre at south 
c                boundary (Pescal/s)
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real pprt(nx,ny,nz,nt)    ! Perturbation pressure (Pascal)

      real pdteb (ny,nz)  ! T-tendency of pprt at e-boundary (Pescal/s)
      real pdtwb (ny,nz)  ! T-tendency of pprt at w-boundary (Pescal/s)
      real pdtnb (nx,nz)  ! T-tendency of pprt at n-boundary (Pescal/s)
      real pdtsb (nx,nz)  ! T-tendency of pprt at s-boundary (Pescal/s)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of pprt on the west boundary (pdtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=1,nz-1
        DO 1110 j=1,ny-1

          tem  =-pprt(2,j,k,3)-pprt(2,j,k,1)+2.*pprt(3,j,k,2)
          cdtdx=(pprt(2,j,k,1)-pprt(2,j,k,3))/(sign(eps,tem)+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)

          pdtwb(j,k)=cdtdx*(-pprt(2,j,k,3)+pprt(1,j,k,2))/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of pprt on the east boundary (pdteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdx must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=1,nz-1
        DO 1210 j=1,ny-1

          tem  =pprt(nx-2,j,k,3)+pprt(nx-2,j,k,1)-2.*pprt(nx-3,j,k,2)
          cdtdx=(pprt(nx-2,j,k,1)-pprt(nx-2,j,k,3))/(sign(eps,tem)+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)

          pdteb(j,k)=cdtdx*(pprt(nx-2,j,k,3)-pprt(nx-1,j,k,2))/
     :                (dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of pprt on the south boundary (pdtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1310 k=1,nz-1
        DO 1310 i=1,nx-1

          tem  =-pprt(i,2,k,3)-pprt(i,2,k,1)+2.*pprt(i,3,k,2)
          cdtdy=(pprt(i,2,k,1)-pprt(i,2,k,3))/
     :          (sign(eps,tem)+tem)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          pdtsb(i,k)=cdtdy*(-pprt(i,2,k,3)+pprt(i,1,k,2))/
     :               (dtbig1*(1.-cdtdy))

1310    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of pprt on the north boundary (pdtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=1,nz-1
        DO 1410 i=1,nx-1

          tem  =pprt(i,ny-2,k,3)+pprt(i,ny-2,k,1)-2.*pprt(i,ny-3,k,2)
          cdtdy=(pprt(i,ny-2,k,1)-pprt(i,ny-2,k,3))/
     :          (sign(eps,tem)+tem)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          pdtnb(i,k)=cdtdy*(pprt(i,ny-2,k,3)-pprt(i,ny-1,k,2))/
     :                (dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southwest corner of p boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=1,nz-1

          temx =(-pprt(2,2,k,3)-pprt(2,2,k,1)
     :          +2.*pprt(3,2,k,2))/(2.*dx)
          temy =(-pprt(2,2,k,3)-pprt(2,2,k,1)
     :          +2.*pprt(2,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(pprt(2,2,k,1)-pprt(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(pprt(2,2,k,1)-pprt(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          pdtsb(1,k)=(cdtdx*(-pprt(2,2,k,3)+pprt(1,2,k,2))
     :               +cdtdy*(-pprt(2,2,k,3)+pprt(2,1,k,2)))/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          pdtwb(1,k)=pdtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of p boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=1,nz-1

          temx =(pprt(nx-2,2,k,3)+pprt(nx-2,2,k,1)
     :          -2.*pprt(nx-3,2,k,2))/(2.*dx)
          temy =(-pprt(nx-2,2,k,3)-pprt(nx-2,2,k,1)
     :          +2.*pprt(nx-2,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(pprt(nx-2,2,k,1)-pprt(nx-2,2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(pprt(nx-2,2,k,1)-pprt(nx-2,2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          pdtsb(nx-1,k)=(cdtdx*(pprt(nx-2,2,k,3)-pprt(nx-1,2,k,2))
     :                  +cdtdy*(-pprt(nx-2,2,k,3)+pprt(nx-2,1,k,2)))/
     :                  (dtbig1*(1.+cdtdx-cdtdy))
          pdteb( 1  ,k)=pdtsb(nx-1,k)

1520    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of p boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=1,nz-1

          temx =(-pprt(2,ny-2,k,3)-pprt(2,ny-2,k,1)
     :          +2.*pprt(3,ny-2,k,2))/(2.*dx)
          temy =(pprt(2,ny-2,k,3)+pprt(2,ny-2,k,1)
     :          -2.*pprt(2,ny-3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(pprt(2,ny-2,k,1)-pprt(2,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(pprt(2,ny-2,k,1)-pprt(2,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          pdtnb( 1  ,k)=(cdtdx*(-pprt(2,ny-2,k,3)+pprt(1,ny-2,k,2))
     :                  +cdtdy*(pprt(2,ny-2,k,3)-pprt(2,ny-1,k,2)))/
     :                  (dtbig1*(1.-cdtdx+cdtdy))
          pdtwb(ny-1,k)=pdtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of p boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=1,nz-1

          temx =(pprt(nx-2,ny-2,k,3)+pprt(nx-2,ny-2,k,1)
     :          -2.*pprt(nx-3,ny-2,k,2))/(2.*dx)
          temy =(pprt(nx-2,ny-2,k,3)+pprt(nx-2,ny-2,k,1)
     :          -2.*pprt(nx-2,ny-3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(pprt(nx-2,ny-2,k,1)-pprt(nx-2,ny-2,k,3))/(2.*dx)
     :          *temx/(eps+tem)
          cdtdy=(pprt(nx-2,ny-2,k,1)-pprt(nx-2,ny-2,k,3))/(2.*dy)
     :          *temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          pdtnb(nx-1,k)=(cdtdx*(pprt(nx-2,ny-2,k,3)-pprt(nx-1,ny-2,k,2))
     :                  +cdtdy*(pprt(nx-2,ny-2,k,3)-pprt(nx-2,ny-1,k,2))
     :                  )/(dtbig1*(1.+cdtdx+cdtdy))
          pdteb(ny-1,k)=pdtnb(nx-1,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTPT1D                    ######
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 BDTPT1D(nx,ny,nz,dtbig1,
     :           ptprt,ptdteb,ptdtwb,ptdtnb,ptdtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of the potential temperature perturbation
c     on the lateral boundaries using the radiation boundary condition of 
c     Orlanski.
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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       ptprt    Perturbation potential temperature at all time levels (K)
c
c     OUTPUT:
c
c       ptdteb    Time tendency of pt at east  boundary (K/s)   
c       ptdtwb    Time tendency of pt at west  boundary (K/s)   
c       ptdtnb    Time tendency of pt at north boundary (K/s)   
c       ptdtsb    Time tendency of pt at south boundary (K/s)   
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real ptprt(nx,ny,nz,nt) ! Perturbation potential temperature (K).

      real ptdteb(ny,nz)  ! T-tendency of ptprt at e-boundary (K/s)
      real ptdtwb(ny,nz)  ! T-tendency of ptprt at w-boundary (K/s)
      real ptdtnb(nx,nz)  ! T-tendency of ptprt at n-boundary (K/s)
      real ptdtsb(nx,nz)  ! T-tendency of ptprt at s-boundary (K/s)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of ptprt on the west boundary (ptdtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=1,nz-1
        DO 1110 j=1,ny-1

          tem  =-ptprt(2,j,k,3)-ptprt(2,j,k,1)+2.*ptprt(3,j,k,2)
          cdtdx=(ptprt(2,j,k,1)-ptprt(2,j,k,3))/(sign(eps,tem)+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)

          ptdtwb(j,k)=cdtdx*(-ptprt(2,j,k,3)+ptprt(1,j,k,2))/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of ptprt on the east boundary (ptdteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdx must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=1,nz-1
        DO 1210 j=1,ny-1

          tem  =ptprt(nx-2,j,k,3)+ptprt(nx-2,j,k,1)-2.*ptprt(nx-3,j,k,2)
          cdtdx=(ptprt(nx-2,j,k,1)-ptprt(nx-2,j,k,3))/
     :          (sign(eps,tem)+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)

          ptdteb(j,k)=cdtdx*(ptprt(nx-2,j,k,3)-ptprt(nx-1,j,k,2))/
     :                (dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of ptprt on the south boundary (ptdtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1310 k=1,nz-1
        DO 1310 i=1,nx-1

          tem  =-ptprt(i,2,k,3)-ptprt(i,2,k,1)+2.*ptprt(i,3,k,2)
          cdtdy=(ptprt(i,2,k,1)-ptprt(i,2,k,3))/
     :          (sign(eps,tem)+tem)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          ptdtsb(i,k)=cdtdy*(-ptprt(i,2,k,3)+ptprt(i,1,k,2))/
     :                (dtbig1*(1.-cdtdy))

1310    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of ptprt on the north boundary (ptdtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=1,nz-1
        DO 1410 i=1,nx-1

          tem  =ptprt(i,ny-2,k,3)+ptprt(i,ny-2,k,1)-2.*ptprt(i,ny-3,k,2)
          cdtdy=(ptprt(i,ny-2,k,1)-ptprt(i,ny-2,k,3))/
     :          (sign(eps,tem)+tem)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          ptdtnb(i,k)=cdtdy*(ptprt(i,ny-2,k,3)-ptprt(i,ny-1,k,2))/
     :                (dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southwest corner of pt boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=1,nz-1

          temx =(-ptprt(2,2,k,3)-ptprt(2,2,k,1)
     :          +2.*ptprt(3,2,k,2))/(2.*dx)
          temy =(-ptprt(2,2,k,3)-ptprt(2,2,k,1)
     :          +2.*ptprt(2,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(ptprt(2,2,k,1)-ptprt(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(ptprt(2,2,k,1)-ptprt(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          ptdtsb(1,k)=(cdtdx*(-ptprt(2,2,k,3)+ptprt(1,2,k,2))
     :               +cdtdy*(-ptprt(2,2,k,3)+ptprt(2,1,k,2)))/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          ptdtwb(1,k)=ptdtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of pt boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=1,nz-1

          temx =(ptprt(nx-2,2,k,3)+ptprt(nx-2,2,k,1)
     :          -2.*ptprt(nx-3,2,k,2))/(2.*dx)
          temy =(-ptprt(nx-2,2,k,3)-ptprt(nx-2,2,k,1)
     :          +2.*ptprt(nx-2,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(ptprt(nx-2,2,k,1)-ptprt(nx-2,2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(ptprt(nx-2,2,k,1)-ptprt(nx-2,2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          ptdtsb(nx-1,k)=(cdtdx*(ptprt(nx-2,2,k,3)-ptprt(nx-1,2,k,2))
     :                  +cdtdy*(-ptprt(nx-2,2,k,3)+ptprt(nx-2,1,k,2))
     :                  )/(dtbig1*(1.+cdtdx-cdtdy))
          ptdteb( 1  ,k)=ptdtsb(nx-1,k)

1520    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of pt boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=1,nz-1

          temx =(-ptprt(2,ny-2,k,3)-ptprt(2,ny-2,k,1)
     :          +2.*ptprt(3,ny-2,k,2))/(2.*dx)
          temy =(ptprt(2,ny-2,k,3)+ptprt(2,ny-2,k,1)
     :          -2.*ptprt(2,ny-3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(ptprt(2,ny-2,k,1)-ptprt(2,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(ptprt(2,ny-2,k,1)-ptprt(2,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          ptdtnb( 1  ,k)=(cdtdx*(-ptprt(2,ny-2,k,3)+ptprt(1,ny-2,k,2))
     :                  +cdtdy*(ptprt(2,ny-2,k,3)-ptprt(2,ny-1,k,2)))/
     :                  (dtbig1*(1.-cdtdx+cdtdy))
          ptdtwb(ny-1,k)=ptdtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of pt boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=1,nz-1

          temx =(ptprt(nx-2,ny-2,k,3)+ptprt(nx-2,ny-2,k,1)
     :          -2.*ptprt(nx-3,ny-2,k,2))/(2.*dx)
          temy =(ptprt(nx-2,ny-2,k,3)+ptprt(nx-2,ny-2,k,1)
     :          -2.*ptprt(nx-2,ny-3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(ptprt(nx-2,ny-2,k,1)-ptprt(nx-2,ny-2,k,3))/(2.*dx)
     :          *temx/(eps+tem)
          cdtdy=(ptprt(nx-2,ny-2,k,1)-ptprt(nx-2,ny-2,k,3))/(2.*dy)
     :          *temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          ptdtnb(nx-1,k)=(cdtdx*(ptprt(nx-2,ny-2,k,3)
     :                  -ptprt(nx-1,ny-2,k,2))
     :                  +cdtdy*(ptprt(nx-2,ny-2,k,3)
     :                  -ptprt(nx-2,ny-1,k,2))
     :                  )/(dtbig1*(1.+cdtdx+cdtdy))
          ptdteb(ny-1,k)=ptdtnb(nx-1,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTQ1D                     ######
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 BDTQ1D(nx,ny,nz,dtbig1, q,qdteb,qdtwb,qdtnb,qdtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of q on the lateral boundaries
c     using the radiation boundary condition of Orlanski.
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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       q        Mixing ratio for one of the water variables (kg/kg)
c
c     OUTPUT:
c
c       qdteb    Time tendency of q at east  boundary (kg/(kg*s))
c       qdtwb    Time tendency of q at west  boundary (kg/(kg*s))
c       qdtnb    Time tendency of q at north boundary (kg/(kg*s))
c       qdtsb    Time tendency of q at south boundary (kg/(kg*s))
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real q(nx,ny,nz,nt)  ! Mixing ratio for one of the water
                           ! variables (kg/kg)

      real qdteb (ny,nz)   ! T-tendency of q at e-boundary (kg/(kg*s))
      real qdtwb (ny,nz)   ! T-tendency of q at w-boundary (kg/(kg*s))
      real qdtnb (nx,nz)   ! T-tendency of q at n-boundary (kg/(kg*s))
      real qdtsb (nx,nz)   ! T-tendency of q at s-boundary (kg/(kg*s))
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of q on the west boundary (qdtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=1,nz-1
        DO 1110 j=1,ny-1

          tem  =-q(2,j,k,3)-q(2,j,k,1)+2.*q(3,j,k,2)
          cdtdx=(q(2,j,k,1)-q(2,j,k,3))/(sign(eps,tem)+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)

          qdtwb(j,k)=cdtdx*(-q(2,j,k,3)+q(1,j,k,2))/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of q on the east boundary (qdteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdx must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=1,nz-1
        DO 1210 j=1,ny-1

          tem  =q(nx-2,j,k,3)+q(nx-2,j,k,1)-2.*q(nx-3,j,k,2)
          cdtdx=(q(nx-2,j,k,1)-q(nx-2,j,k,3))/(sign(eps,tem)+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)

          qdteb(j,k)=cdtdx*(q(nx-2,j,k,3)-q(nx-1,j,k,2))/
     :               (dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of q on the south boundary (qdtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1310 k=1,nz-1
        DO 1310 i=1,nx-1

          tem  =-q(i,2,k,3)-q(i,2,k,1)+2.*q(i,3,k,2) 
          cdtdy=(q(i,2,k,1)-q(i,2,k,3))/(sign(eps,tem)+tem)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          qdtsb(i,k)=cdtdy*(-q(i,2,k,3)+q(i,1,k,2))/
     :               (dtbig1*(1.-cdtdy))

1310    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of q on the north boundary (qdtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=1,nz-1
        DO 1410 i=1,nx-1

          tem  =q(i,ny-2,k,3)+q(i,ny-2,k,1)-2.*q(i,ny-3,k,2)
          cdtdy=(q(i,ny-2,k,1)-q(i,ny-2,k,3))/(sign(eps,tem)+tem)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          qdtnb(i,k)=cdtdy*(q(i,ny-2,k,3)-q(i,ny-1,k,2))/
     :               (dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southwest corner of q boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=2,nz-1

          temx =(-q(2,2,k,3)-q(2,2,k,1)+2.*q(3,2,k,2))/(2.*dx)
          temy =(-q(2,2,k,3)-q(2,2,k,1)+2.*q(2,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(q(2,2,k,1)-q(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(q(2,2,k,1)-q(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          qdtsb(1,k)=(cdtdx*(-q(2,2,k,3)+q(1,2,k,2))
     :               +cdtdy*(-q(2,2,k,3)+q(2,1,k,2)))/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          qdtwb(1,k)=qdtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of q boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=2,nz-1

          temx =(q(nx-2,2,k,3)+q(nx-2,2,k,1)-2.*q(nx-3,2,k,2))/(2.*dx)
          temy =(-q(nx-2,2,k,3)-q(nx-2,2,k,1)+2.*q(nx-2,3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(q(nx-2,2,k,1)-q(nx-2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(q(nx-2,2,k,1)-q(nx-2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)

          qdtsb(nx-1,k)=(cdtdx*(q(nx-2,2,k,3)-q(nx-1,2,k,2))
     :                  +cdtdy*(-q(nx-2,2,k,3)+q(nx-2,1,k,2)))/
     :                  (dtbig1*(1.+cdtdx-cdtdy))
          qdteb( 1  ,k)=qdtsb(nx-1,k)

1520    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of q boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=2,nz-1

          temx =(-q(2,ny-2,k,3)-q(2,ny-2,k,1)+2.*q(3,ny-2,k,2))/(2.*dx)
          temy =(q(2,ny-2,k,3)+q(2,ny-2,k,1)-2.*q(2,ny-3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(q(2,ny-2,k,1)-q(2,ny-2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(q(2,ny-2,k,1)-q(2,ny-2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          qdtnb( 1  ,k)=(cdtdx*(-q(2,ny-2,k,3)+q(1,ny-2,k,2))
     :                  +cdtdy*(q(2,ny-2,k,3)-q(2,ny-1,k,2)))/
     :                  (dtbig1*(1.-cdtdx+cdtdy))
          qdtwb(ny-1,k)=qdtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of q boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=2,nz-1

          temx =(q(nx-2,ny-2,k,3)+q(nx-2,ny-2,k,1)
     :          -2.*q(nx-3,ny-2,k,2))/(2.*dx)
          temy =(q(nx-2,ny-2,k,3)+q(nx-2,ny-2,k,1)
     :          -2.*q(nx-2,ny-3,k,2))/(2.*dy)
          tem  =temx**2+temy**2
          cdtdx=(q(nx-2,ny-2,k,1)-q(nx-2,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(q(nx-2,ny-2,k,1)-q(nx-2,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)

          qdtnb(nx-1,k)=(cdtdx*(q(nx-2,ny-2,k,3)-q(nx-1,ny-2,k,2))
     :                  +cdtdy*(q(nx-2,ny-2,k,3)-q(nx-2,ny-1,k,2)))/
     :                  (dtbig1*(1.+cdtdx+cdtdy))
          qdteb(ny-1,k)=qdtnb(nx-1,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTU3D                     ######
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 BDTU3D(nx,ny,nz,dtbig1, u,udteb,udtwb,udtnb,udtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of u on the lateral boundaries,
c     using the radiation boundary condition of Raymond and Kuo.
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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       u        x-component of velocity at all time levels (m/s).
c
c     OUTPUT:
c
c       udteb    Time tendency of u at east  boundary (m/s**2)
c       udtwb    Time tendency of u at west  boundary (m/s**2)
c       udtnb    Time tendency of u at north boundary (m/s**2)
c       udtsb    Time tendency of u at south boundary (m/s**2)
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real u (nx,ny,nz,nt) ! Total u-velocity (m/s).

      real udteb (ny,nz)  ! T-tendency of u at e-boundary (m/s**2)
      real udtwb (ny,nz)  ! T-tendency of u at w-boundary (m/s**2)
      real udtnb (nx,nz)  ! T-tendency of u at n-boundary (m/s**2)
      real udtsb (nx,nz)  ! T-tendency of u at s-boundary (m/s**2)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy, temz

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
      real cdtdz           ! Estimated phase speed times dt/dz .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of u on the west boundary (udtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=1,nz-1
        DO 1110 j=2,ny-2

          temx =(-u(2,j,k,3)-u(2,j,k,1)+2.*u(3,j,k,2))/(2.*dx)
          temy =(u(2,j+1,k,2)-u(2,j-1,k,2))/(2.*dy)
          temz =(u(2,j,k+1,2)-u(2,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(u(2,j,k,1)-u(2,j,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(u(2,j,k,1)-u(2,j,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(u(2,j,k,1)-u(2,j,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          udtwb(j,k)=(cdtdx*(-u(2,j,k,3)+u(1,j,k,2))
     :               -cdtdy*(u(2,j+1,k,3)-u(2,j-1,k,3))/2.
     :               -cdtdz*(u(2,j,k+1,3)-u(2,j,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of u on the east boundary (udteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdx must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=1,nz-1
        DO 1210 j=2,ny-2

          temx =(u(nx-1,j,k,3)+u(nx-1,j,k,1)-2.*u(nx-2,j,k,2))/(2.*dx)
          temy =(u(nx-1,j+1,k,2)-u(nx-1,j-1,k,2))/(2.*dy)
          temz =(u(nx-1,j,k+1,2)-u(nx-1,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(u(nx-1,j,k,1)-u(nx-1,j,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(u(nx-1,j,k,1)-u(nx-1,j,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(u(nx-1,j,k,1)-u(nx-1,j,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          udteb(j,k)=(cdtdx*(u(nx-1,j,k,3)-u(nx,j,k,2))
     :               -cdtdy*(u(nx-1,j+1,k,3)-u(nx-1,j-1,k,3))/2.
     :               -cdtdz*(u(nx-1,j,k+1,3)-u(nx-1,j,k-1,3))/2.)/
     :               (dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of u on the south boundary (udtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1310 k=1,nz-1
        DO 1310 i=2,nx-1

          temx =(u(i+1,2,k,2)-u(i-1,2,k,2))/(2.*dx)
          temy =(-u(i,2,k,3)-u(i,2,k,1)+2.*u(i,3,k,2))/(2.*dy)
          temz =(u(i,2,k+1,2)-u(i,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(u(i,2,k,1)-u(i,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(u(i,2,k,1)-u(i,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(u(i,2,k,1)-u(i,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          udtsb(i,k)=(-cdtdx*(u(i+1,2,k,3)-u(i-1,2,k,3))/2.
     :               +cdtdy*(-u(i,2,k,3)+u(i,1,k,2))
     :               -cdtdz*(u(i,2,k+1,3)-u(i,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdy))

1310    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of u on the north boundary (udtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=1,nz-1
        DO 1410 i=2,nx-1

          temx =(u(i+1,ny-2,k,2)-u(i-1,ny-2,k,2))/(2.*dx)
          temy =(u(i,ny-2,k,3)+u(i,ny-2,k,1)-2.*u(i,ny-3,k,2))/(2.*dy)
          temz =(u(i,ny-2,k+1,2)-u(i,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(u(i,ny-2,k,1)-u(i,ny-2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(u(i,ny-2,k,1)-u(i,ny-2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(u(i,ny-2,k,1)-u(i,ny-2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          udtnb(i,k)=(-cdtdx*(u(i+1,ny-2,k,3)-u(i-1,ny-2,k,3))/2.
     :               +cdtdy*(u(i,ny-2,k,3)-u(i,ny-1,k,2))
     :               -cdtdz*(u(i,ny-2,k+1,3)-u(i,ny-2,k-1,3))/2.)/
     :               (dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c#######################################################################
c
c     Calculate the time tendency on the southwest corner of u boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=1,nz-1

          temx =(-u(2,2,k,3)-u(2,2,k,1)+2.*u(3,2,k,2))/(2.*dx)
          temy =(-u(2,2,k,3)-u(2,2,k,1)+2.*u(2,3,k,2))/(2.*dy)
          temz =(u(2,2,k+1,2)-u(2,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(u(2,2,k,1)-u(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(u(2,2,k,1)-u(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(u(2,2,k,1)-u(2,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          udtsb(1,k)=(cdtdx*(-u(2,2,k,3)+u(1,2,k,2))
     :               +cdtdy*(-u(2,2,k,3)+u(2,1,k,2))
     :               -cdtdz*(u(2,2,k+1,3)-u(2,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          udtwb(1,k)=udtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of u boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=1,nz-1

          temx =(u(nx-1,2,k,3)+u(nx-1,2,k,1)-2.*u(nx-2,2,k,2))/(2.*dx)
          temy =(-u(nx-1,2,k,3)-u(nx-1,2,k,1)+2.*u(nx-1,3,k,2))/(2.*dy)
          temz =(u(nx-1,2,k+1,2)-u(nx-1,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(u(nx-1,2,k,1)-u(nx-1,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(u(nx-1,2,k,1)-u(nx-1,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(u(nx-1,2,k,1)-u(nx-1,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          udtsb(nx,k)=(cdtdx*(u(nx-1,2,k,3)-u(nx,2,k,2))
     :                +cdtdy*(-u(nx-1,2,k,3)+u(nx-1,1,k,2))
     :                -cdtdz*(u(nx-1,2,k+1,3)-u(nx-1,2,k-1,3))/2.)/
     :                (dtbig1*(1.+cdtdx-cdtdy))
          udteb( 1,k)=udtsb(nx,k)

1520    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of u boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=1,nz-1

          temx =(-u(2,ny-2,k,3)-u(2,ny-2,k,1)+2.*u(3,ny-2,k,2))/(2.*dx)
          temy =(u(2,ny-2,k,3)+u(2,ny-2,k,1)-2.*u(2,ny-3,k,2))/(2.*dy)
          temz =(u(2,ny-2,k+1,2)-u(2,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(u(2,ny-2,k,1)-u(2,ny-2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(u(2,ny-2,k,1)-u(2,ny-2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(u(2,ny-2,k,1)-u(2,ny-2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          udtnb( 1  ,k)=(cdtdx*(-u(2,ny-2,k,3)+u(1,ny-2,k,2))
     :                  +cdtdy*(u(2,ny-2,k,3)-u(2,ny-1,k,2))
     :                  -cdtdz*(u(2,ny-2,k+1,3)-u(2,ny-2,k-1,3))/2.)/
     :                  (dtbig1*(1.-cdtdx+cdtdy))
          udtwb(ny-1,k)=udtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of u boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=1,nz-1

          temx =(u(nx-1,ny-2,k,3)+u(nx-1,ny-2,k,1)
     :          -2.*u(nx-2,ny-2,k,2))/(2.*dx)
          temy =(u(nx-1,ny-2,k,3)+u(nx-1,ny-2,k,1)
     :          -2.*u(nx-1,ny-3,k,2))/(2.*dy)
          temz =(u(nx-1,ny-2,k+1,2)-u(nx-1,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(u(nx-1,ny-2,k,1)-u(nx-1,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(u(nx-1,ny-2,k,1)-u(nx-1,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(u(nx-1,ny-2,k,1)-u(nx-1,ny-2,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          udtnb(nx  ,k)=(cdtdx*(u(nx-1,ny-2,k,3)-u(nx,ny-2,k,2))
     :                  +cdtdy*(u(nx-1,ny-2,k,3)-u(nx-1,ny-1,k,2))
     :                  -cdtdz*(u(nx-1,ny-2,k+1,3)-u(nx-1,ny-2,k-1,3))
     :                  /2.)/(dtbig1*(1.+cdtdx+cdtdy))
          udteb(ny-1,k)=udtnb(nx,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTV3D                     ######
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 BDTV3D(nx,ny,nz,dtbig1, v,vdteb,vdtwb,vdtnb,vdtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of v on the lateral boundaries,
c     using the radiation boundary condition of Raymond and Kuo. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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       v        y-component of velocity at all time levels (m/s).
c
c     OUTPUT:
c
c       vdteb    Time tendency of v at east  boundary (m/s**2)
c       vdtwb    Time tendency of v at west  boundary (m/s**2)
c       vdtnb    Time tendency of v at north boundary (m/s**2)
c       vdtsb    Time tendency of v at south boundary (m/s**2)
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real v (nx,ny,nz,nt) ! Total v-velocity (m/s).

      real vdteb (ny,nz)  ! T-tendency of v at e-boundary (m/s**2)
      real vdtwb (ny,nz)  ! T-tendency of v at w-boundary (m/s**2)
      real vdtnb (nx,nz)  ! T-tendency of v at n-boundary (m/s**2)
      real vdtsb (nx,nz)  ! T-tendency of v at s-boundary (m/s**2)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy, temz

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
      real cdtdz           ! Estimated phase speed times dt/dz .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of v on the west boundary (vdtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=1,nz-1
        DO 1110 j=2,ny-1

          temx =(-v(2,j,k,3)-v(2,j,k,1)+2.*v(3,j,k,2))/(2.*dx)
          temy =(v(2,j+1,k,2)-v(2,j-1,k,2))/(2.*dy)
          temz =(v(2,j,k+1,2)-v(2,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(v(2,j,k,1)-v(2,j,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(v(2,j,k,1)-v(2,j,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(v(2,j,k,1)-v(2,j,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          vdtwb(j,k)=(cdtdx*(-v(2,j,k,3)+v(1,j,k,2))
     :               -cdtdy*(v(2,j+1,k,3)-v(2,j-1,k,3))/2.
     :               -cdtdz*(v(2,j,k+1,3)-v(2,j,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of v on the east boundary (vdteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdx must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=1,nz-1
        DO 1210 j=2,ny-1

          temx =(v(nx-2,j,k,3)+v(nx-2,j,k,1)-2.*v(nx-3,j,k,2))/(2.*dx)
          temy =(v(nx-2,j+1,k,2)-v(nx-2,j-1,k,2))/(2.*dy)
          temz =(v(nx-2,j,k+1,2)-v(nx-2,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(v(nx-2,j,k,1)-v(nx-2,j,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(v(nx-2,j,k,1)-v(nx-2,j,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(v(nx-2,j,k,1)-v(nx-2,j,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          vdteb(j,k)=(cdtdx*(v(nx-2,j,k,3)-v(nx-1,j,k,2))
     :               -cdtdy*(v(nx-2,j+1,k,3)-v(nx-2,j-1,k,3))/2.
     :               -cdtdz*(v(nx-2,j,k+1,3)-v(nx-2,j,k-1,3))/2.)/
     :               (dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of v on the south boundary (vdtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1310 k=1,nz-1
        DO 1310 i=2,nx-2

          temx =(v(i+1,2,k,2)-v(i-1,2,k,2))/(2.*dx)
          temy =(-v(i,2,k,3)-v(i,2,k,1)+2.*v(i,3,k,2))/(2.*dy)
          temz =(v(i,2,k+1,2)-v(i,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(v(i,2,k,1)-v(i,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(v(i,2,k,1)-v(i,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(v(i,2,k,1)-v(i,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          vdtsb(i,k)=(-cdtdx*(v(i+1,2,k,3)-v(i-1,2,k,3))/2.
     :               +cdtdy*(-v(i,2,k,3)+v(i,1,k,2))
     :               -cdtdz*(v(i,2,k+1,3)-v(i,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdy))

1310  CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of v on the north boundary (vdtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=1,nz-1
        DO 1410 i=2,nx-2

          temx =(v(i+1,ny-1,k,2)-v(i-1,ny-1,k,2))/(2.*dx)
          temy =(v(i,ny-1,k,3)+v(i,ny-1,k,1)-2.*v(i,ny-2,k,2))/(2.*dy)
          temz =(v(i,ny-1,k+1,2)-v(i,ny-1,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(v(i,ny-1,k,1)-v(i,ny-1,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(v(i,ny-1,k,1)-v(i,ny-1,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(v(i,ny-1,k,1)-v(i,ny-1,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          vdtnb(i,k)=(-cdtdx*(v(i+1,ny-1,k,3)-v(i-1,ny-1,k,3))/2.
     :               +cdtdy*(v(i,ny-1,k,3)-v(i,ny,k,2))
     :               -cdtdz*(v(i,ny-1,k+1,3)-v(i,ny-1,k-1,3))/2.)/
     :               (dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southwest corner of v boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=1,nz-1

          temx =(-v(2,2,k,3)-v(2,2,k,1)+2.*v(3,2,k,2))/(2.*dx)
          temy =(-v(2,2,k,3)-v(2,2,k,1)+2.*v(2,3,k,2))/(2.*dy)
          temz =(v(2,2,k+1,2)-v(2,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(v(2,2,k,1)-v(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(v(2,2,k,1)-v(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(v(2,2,k,1)-v(2,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          vdtsb(1,k)=(cdtdx*(-v(2,2,k,3)+v(1,2,k,2))
     :               +cdtdy*(-v(2,2,k,3)+v(2,1,k,2))
     :               -cdtdz*(v(2,2,k+1,3)-v(2,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          vdtwb(1,k)=vdtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of v boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=1,nz-1

          temx =(v(nx-2,2,k,3)+v(nx-2,2,k,1)-2.*v(nx-3,2,k,2))/(2.*dx)
          temy =(-v(nx-2,2,k,3)-v(nx-2,2,k,1)+2.*v(nx-2,3,k,2))/(2.*dy)
          temz =(v(nx-2,2,k+1,2)-v(nx-2,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(v(nx-2,2,k,1)-v(nx-2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(v(nx-2,2,k,1)-v(nx-2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(v(nx-2,2,k,1)-v(nx-2,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          vdtsb(nx-1,k)=(cdtdx*(v(nx-2,2,k,3)-v(nx-1,2,k,2))
     :                  +cdtdy*(-v(nx-2,2,k,3)+v(nx-2,1,k,2))
     :                  -cdtdz*(v(nx-2,2,k+1,3)-v(nx-2,2,k-1,3))/2.)/
     :                  (dtbig1*(1.+cdtdx-cdtdy))
          vdteb( 1  ,k)=vdtsb(nx-1,k)

1520    CONTINUE

      ENDIF
c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of v boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=1,nz-1

          temx =(-v(2,ny-1,k,3)-v(2,ny-1,k,1)+2.*v(3,ny-1,k,2))/(2.*dx)
          temy =(v(2,ny-1,k,3)+v(2,ny-1,k,1)-2.*v(2,ny-2,k,2))/(2.*dy)
          temz =(v(2,ny-1,k+1,2)-v(2,ny-1,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(v(2,ny-1,k,1)-v(2,ny-1,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(v(2,ny-1,k,1)-v(2,ny-1,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(v(2,ny-1,k,1)-v(2,ny-1,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          vdtnb( 1,k)=(cdtdx*(-v(2,ny-1,k,3)+v(1,ny-1,k,2))
     :                +cdtdy*(v(2,ny-1,k,3)-v(2,ny,k,2))
     :                -cdtdz*(v(2,ny-1,k+1,3)-v(2,ny-1,k-1,3))/2.)/
     :                (dtbig1*(1.-cdtdx+cdtdy))
          vdtwb(ny,k)=vdtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of v boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=1,nz-1

          temx =(v(nx-2,ny-1,k,3)+v(nx-2,ny-1,k,1)
     :          -2.*v(nx-3,ny-1,k,2))/(2.*dx)
          temy =(v(nx-2,ny-1,k,3)+v(nx-2,ny-1,k,1)
     :          -2.*v(nx-2,ny-2,k,2))/(2.*dy)
          temz =(v(nx-2,ny-1,k+1,2)-v(nx-2,ny-1,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(v(nx-2,ny-1,k,1)-v(nx-2,ny-1,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(v(nx-2,ny-1,k,1)-v(nx-2,ny-1,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(v(nx-2,ny-1,k,1)-v(nx-2,ny-1,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          vdtnb(nx-1,k)=(cdtdx*(v(nx-2,ny-1,k,3)-v(nx-1,ny-1,k,2))
     :                  +cdtdy*(v(nx-2,ny-1,k,3)-v(nx-2,ny,k,2))
     :                  -cdtdz*(v(nx-2,ny-1,k+1,3)-v(nx-2,ny-1,k-1,3))
     :                  /2.)/(dtbig1*(1.+cdtdx+cdtdy))
          vdteb(ny  ,k)=vdtnb(nx-1,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTW3D                     ######
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 BDTW3D(nx,ny,nz,dtbig1, w,wdteb,wdtwb,wdtnb,wdtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of w on the lateral boundaries,
c     using the radiation boundary condition of Raymond and Kuo.
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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       w        z-component of velocity at all time levels (m/s).
c
c     OUTPUT:
c
c       wdteb    Time tendency of w at east  boundary (m/s**2)
c       wdtwb    Time tendency of w at west  boundary (m/s**2)
c       wdtnb    Time tendency of w at north boundary (m/s**2)
c       wdtsb    Time tendency of w at south boundary (m/s**2)
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real w (nx,ny,nz,nt) ! Total w-velocity (m/s).

      real wdteb (ny,nz)  ! T-tendency of w at e-boundary (m/s**2)
      real wdtwb (ny,nz)  ! T-tendency of w at w-boundary (m/s**2)
      real wdtnb (nx,nz)  ! T-tendency of w at n-boundary (m/s**2)
      real wdtsb (nx,nz)  ! T-tendency of w at s-boundary (m/s**2)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy, temz

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
      real cdtdz           ! Estimated phase speed times dt/dz .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of w on the west boundary (wdtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=1,nz
        DO 1110 j=2,ny-2

          temx =(-w(2,j,k,3)-w(2,j,k,1)+2.*w(3,j,k,2))/(2.*dx)
          temy =(w(2,j+1,k,2)-w(2,j-1,k,2))/(2.*dy)
          temz =(w(2,j,k+1,2)-w(2,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(w(2,j,k,1)-w(2,j,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(w(2,j,k,1)-w(2,j,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(w(2,j,k,1)-w(2,j,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          wdtwb(j,k)=(cdtdx*(-w(2,j,k,3)+w(1,j,k,2))
     :               -cdtdy*(w(2,j+1,k,3)-w(2,j-1,k,3))/2.
     :               -cdtdz*(w(2,j,k+1,3)-w(2,j,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of w on the east boundary (wdteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=1,nz
        DO 1210 j=2,ny-2

          temx =(w(nx-2,j,k,3)+w(nx-2,j,k,1)-2.*w(nx-3,j,k,2))/(2.*dx)
          temy =(w(nx-2,j+1,k,2)-w(nx-2,j-1,k,2))/(2.*dy)
          temz =(w(nx-2,j,k+1,2)-w(nx-2,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(w(nx-2,j,k,1)-w(nx-2,j,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(w(nx-2,j,k,1)-w(nx-2,j,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(w(nx-2,j,k,1)-w(nx-2,j,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          wdteb(j,k)=(cdtdx*(w(nx-2,j,k,3)-w(nx-1,j,k,2))
     :               -cdtdy*(w(nx-2,j+1,k,3)-w(nx-2,j-1,k,3))/2.
     :               -cdtdz*(w(nx-2,j,k+1,3)-w(nx-2,j,k-1,3))/2.)/
     :               (dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of w on the south boundary (wdtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1310 k=1,nz
        DO 1310 i=2,nx-2

          temx =(w(i+1,2,k,2)-w(i-1,2,k,2))/(2.*dx)
          temy =(-w(i,2,k,3)-w(i,2,k,1)+2.*w(i,3,k,2))/(2.*dy)
          temz =(w(i,2,k+1,2)-w(i,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(w(i,2,k,1)-w(i,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(w(i,2,k,1)-w(i,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(w(i,2,k,1)-w(i,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          wdtsb(i,k)=(-cdtdx*(w(i+1,2,k,3)-w(i-1,2,k,3))/2.
     :               +cdtdy*(-w(i,2,k,3)+w(i,1,k,2))
     :               -cdtdz*(w(i,2,k+1,3)-w(i,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdy))

1310    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of w on the north boundary (wdtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=1,nz
        DO 1410 i=2,nx-2

          temx =(w(i+1,ny-2,k,2)-w(i-1,ny-2,k,2))/(2.*dx)
          temy =(w(i,ny-2,k,3)+w(i,ny-2,k,1)-2.*w(i,ny-3,k,2))/(2.*dy)
          temz =(w(i,ny-2,k+1,2)-w(i,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(w(i,ny-2,k,1)-w(i,ny-2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(w(i,ny-2,k,1)-w(i,ny-2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(w(i,ny-2,k,1)-w(i,ny-2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          wdtnb(i,k)=(-cdtdx*(w(i+1,ny-2,k,3)-w(i-1,ny-2,k,3))/2.
     :               +cdtdy*(w(i,ny-2,k,3)-w(i,ny-1,k,2))
     :               -cdtdz*(w(i,ny-2,k+1,3)-w(i,ny-2,k-1,3))/2.)/
     :               (dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southwest corner of w boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=1,nz

          temx =(-w(2,2,k,3)-w(2,2,k,1)+2.*w(3,2,k,2))/(2.*dx)
          temy =(-w(2,2,k,3)-w(2,2,k,1)+2.*w(2,3,k,2))/(2.*dy)
          temz =(w(2,2,k+1,2)-w(2,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(w(2,2,k,1)-w(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(w(2,2,k,1)-w(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(w(2,2,k,1)-w(2,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          wdtsb(1,k)=(cdtdx*(-w(2,2,k,3)+w(1,2,k,2))
     :               +cdtdy*(-w(2,2,k,3)+w(2,1,k,2))
     :               -cdtdz*(w(2,2,k+1,3)-w(2,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          wdtwb(1,k)=wdtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of w boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=1,nz

          temx =(w(nx-2,2,k,3)+w(nx-2,2,k,1)-2.*w(nx-3,2,k,2))/(2.*dx)
          temy =(-w(nx-2,2,k,3)-w(nx-2,2,k,1)+2.*w(nx-2,3,k,2))/(2.*dy)
          temz =(w(nx-2,2,k+1,2)-w(nx-2,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(w(nx-2,2,k,1)-w(nx-2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(w(nx-2,2,k,1)-w(nx-2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(w(nx-2,2,k,1)-w(nx-2,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          wdtsb(nx-1,k)=(cdtdx*(w(nx-2,2,k,3)-w(nx-1,2,k,2))
     :                  +cdtdy*(-w(nx-2,2,k,3)+w(nx-2,1,k,2))
     :                  -cdtdz*(w(nx-2,2,k+1,3)-w(nx-2,2,k-1,3))/2.)/
     :                  (dtbig1*(1.+cdtdx-cdtdy))
          wdteb( 1  ,k)=wdtsb(nx-1,k)

1520    CONTINUE

      ENDIF
c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of w boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=1,nz

          temx =(-w(2,ny-2,k,3)-w(2,ny-2,k,1)+2.*w(3,ny-2,k,2))/(2.*dx)
          temy =(w(2,ny-2,k,3)+w(2,ny-2,k,1)-2.*w(2,ny-3,k,2))/(2.*dy)
          temz =(w(2,ny-2,k+1,2)-w(2,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(w(2,ny-2,k,1)-w(2,ny-2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(w(2,ny-2,k,1)-w(2,ny-2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(w(2,ny-2,k,1)-w(2,ny-2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          wdtnb( 1  ,k)=(cdtdx*(-w(2,ny-2,k,3)+w(1,ny-2,k,2))
     :                  +cdtdy*(w(2,ny-2,k,3)-w(2,ny-1,k,2))
     :                  -cdtdz*(w(2,ny-2,k+1,3)-w(2,ny-2,k-1,3))/2.)/
     :                  (dtbig1*(1.-cdtdx+cdtdy))
          wdtwb(ny-1,k)=wdtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of w boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=1,nz

          temx =(w(nx-2,ny-2,k,3)+w(nx-2,ny-2,k,1)
     :          -2.*w(nx-3,ny-2,k,2))/(2.*dx)
          temy =(w(nx-2,ny-2,k,3)+w(nx-2,ny-2,k,1)
     :          -2.*w(nx-2,ny-3,k,2))/(2.*dy)
          temz =(w(nx-2,ny-2,k+1,2)-w(nx-2,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(w(nx-2,ny-2,k,1)-w(nx-2,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(w(nx-2,ny-2,k,1)-w(nx-2,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(w(nx-2,ny-2,k,1)-w(nx-2,ny-2,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          wdtnb(nx-1,k)=(cdtdx*(w(nx-2,ny-2,k,3)-w(nx-1,ny-2,k,2))
     :                  +cdtdy*(w(nx-2,ny-2,k,3)-w(nx-2,ny-1,k,2))
     :                  -cdtdz*(w(nx-2,ny-2,k+1,3)-w(nx-2,ny-2,k-1,3))
     :                  /2.)/(dtbig1*(1.+cdtdx+cdtdy))
          wdteb(ny-1,k)=wdtnb(nx-1,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTP3D                     ######
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 BDTP3D(nx,ny,nz,dtbig1,pprt,pdteb,pdtwb,pdtnb,pdtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of the perturbation pressure on the lateral 
c     boundaries, using radiation boundary condition of Raymond and Kuo.
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 all time levels (Pascal)
c
c     OUTPUT:
c
c       pdteb    Time tendency of perturbation pressre at east  
c                boundary (Pescal/s)
c       pdtwb    Time tendency of perturbation pressre at west  
c                boundary (Pescal/s)
c       pdtnb    Time tendency of perturbation pressre at north 
c                boundary (Pescal/s)
c       pdtsb    Time tendency of perturbation pressre at south 
c                boundary (Pescal/s)
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real pprt(nx,ny,nz,nt)    ! Perturbation pressure (Pascal)

      real pdteb (ny,nz)  ! T-tendency of pprt at e-boundary (Pescal/s)
      real pdtwb (ny,nz)  ! T-tendency of pprt at w-boundary (Pescal/s)
      real pdtnb (nx,nz)  ! T-tendency of pprt at n-boundary (Pescal/s)
      real pdtsb (nx,nz)  ! T-tendency of pprt at s-boundary (Pescal/s)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy, temz

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
      real cdtdz           ! Estimated phase speed times dt/dz .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of pprt on the west boundary (pdtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=1,nz-1
        DO 1110 j=2,ny-2

          temx =(-pprt(2,j,k,3)-pprt(2,j,k,1)+2.*pprt(3,j,k,2))/
     :          (2.*dx)
          temy =(pprt(2,j+1,k,2)-pprt(2,j-1,k,2))/(2.*dy)
          temz =(pprt(2,j,k+1,2)-pprt(2,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(pprt(2,j,k,1)-pprt(2,j,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(pprt(2,j,k,1)-pprt(2,j,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(pprt(2,j,k,1)-pprt(2,j,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          pdtwb(j,k)=(cdtdx*(-pprt(2,j,k,3)+pprt(1,j,k,2))
     :               -cdtdy*(pprt(2,j+1,k,3)-pprt(2,j-1,k,3))/2.
     :               -cdtdz*(pprt(2,j,k+1,3)-pprt(2,j,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of pprt on the east boundary (pdteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdx must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=1,nz-1
        DO 1210 j=2,ny-2

          temx =(pprt(nx-2,j,k,3)+pprt(nx-2,j,k,1)
     :          -2.*pprt(nx-3,j,k,2))/(2.*dx)
          temy =(pprt(nx-2,j+1,k,2)-pprt(nx-2,j-1,k,2))/(2.*dy)
          temz =(pprt(nx-2,j,k+1,2)-pprt(nx-2,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(pprt(nx-2,j,k,1)-pprt(nx-2,j,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(pprt(nx-2,j,k,1)-pprt(nx-2,j,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(pprt(nx-2,j,k,1)-pprt(nx-2,j,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          pdteb(j,k)=(cdtdx*(pprt(nx-2,j,k,3)-pprt(nx-1,j,k,2))
     :               -cdtdy*(pprt(nx-2,j+1,k,3)-pprt(nx-2,j-1,k,3))/2.
     :               -cdtdz*(pprt(nx-2,j,k+1,3)-pprt(nx-2,j,k-1,3))/2.
     :               )/(dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of pprt on the south boundary (pdtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1310 k=1,nz-1
        DO 1310 i=2,nx-2

          temx =(pprt(i+1,2,k,2)-pprt(i-1,2,k,2))/(2.*dx)
          temy =(-pprt(i,2,k,3)-pprt(i,2,k,1)+2.*pprt(i,3,k,2))/
     :          (2.*dy)
          temz =(pprt(i,2,k+1,2)-pprt(i,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(pprt(i,2,k,1)-pprt(i,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(pprt(i,2,k,1)-pprt(i,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(pprt(i,2,k,1)-pprt(i,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          pdtsb(i,k)=(-cdtdx*(pprt(i+1,2,k,3)-pprt(i-1,2,k,3))/2.
     :               +cdtdy*(-pprt(i,2,k,3)+pprt(i,1,k,2))
     :               -cdtdz*(pprt(i,2,k+1,3)-pprt(i,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdy))

1310    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of pprt on the north boundary (pdtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=1,nz-1
        DO 1410 i=2,nx-2

          temx =(pprt(i+1,ny-2,k,2)-pprt(i-1,ny-2,k,2))/(2.*dx)
          temy =(pprt(i,ny-2,k,3)+pprt(i,ny-2,k,1)
     :          -2.*pprt(i,ny-3,k,2))/(2.*dy)
          temz =(pprt(i,ny-2,k+1,2)-pprt(i,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(pprt(i,ny-2,k,1)-pprt(i,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(pprt(i,ny-2,k,1)-pprt(i,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(pprt(i,ny-2,k,1)-pprt(i,ny-2,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          pdtnb(i,k)=(-cdtdx*(pprt(i+1,ny-2,k,3)
     :                -pprt(i-1,ny-2,k,3))/2.
     :                +cdtdy*(pprt(i,ny-2,k,3)-pprt(i,ny-1,k,2))
     :                -cdtdz*(pprt(i,ny-2,k+1,3)-pprt(i,ny-2,k-1,3)
     :                )/2.)/(dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southwest corner of p boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=1,nz-1

          temx =(-pprt(2,2,k,3)-pprt(2,2,k,1)+2.*pprt(3,2,k,2))/
     :          (2.*dx)
          temy =(-pprt(2,2,k,3)-pprt(2,2,k,1)+2.*pprt(2,3,k,2))/
     :          (2.*dy)
          temz =(pprt(2,2,k+1,2)-pprt(2,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(pprt(2,2,k,1)-pprt(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(pprt(2,2,k,1)-pprt(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(pprt(2,2,k,1)-pprt(2,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          pdtsb(1,k)=(cdtdx*(-pprt(2,2,k,3)+pprt(1,2,k,2))
     :               +cdtdy*(-pprt(2,2,k,3)+pprt(2,1,k,2))
     :               -cdtdz*(pprt(2,2,k+1,3)-pprt(2,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          pdtwb(1,k)=pdtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of p boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=1,nz-1

          temx =(pprt(nx-2,2,k,3)+pprt(nx-2,2,k,1)
     :          -2.*pprt(nx-3,2,k,2))/(2.*dx)
          temy =(-pprt(nx-2,2,k,3)-pprt(nx-2,2,k,1)
     :          +2.*pprt(nx-2,3,k,2))/(2.*dy)
          temz =(pprt(nx-2,2,k+1,2)-pprt(nx-2,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(pprt(nx-2,2,k,1)-pprt(nx-2,2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(pprt(nx-2,2,k,1)-pprt(nx-2,2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(pprt(nx-2,2,k,1)-pprt(nx-2,2,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          pdtsb(nx-1,k)=(cdtdx*(pprt(nx-2,2,k,3)-pprt(nx-1,2,k,2))
     :                  +cdtdy*(-pprt(nx-2,2,k,3)+pprt(nx-2,1,k,2))
     :                  -cdtdz*(pprt(nx-2,2,k+1,3)-pprt(nx-2,2,k-1,3)
     :                  )/2.)/(dtbig1*(1.+cdtdx-cdtdy))
          pdteb( 1  ,k)=pdtsb(nx-1,k)

1520    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of p boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=1,nz-1

          temx =(-pprt(2,ny-2,k,3)-pprt(2,ny-2,k,1)
     :          +2.*pprt(3,ny-2,k,2))/(2.*dx)
          temy =(pprt(2,ny-2,k,3)+pprt(2,ny-2,k,1)
     :          -2.*pprt(2,ny-3,k,2))/(2.*dy)
          temz =(pprt(2,ny-2,k+1,2)-pprt(2,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(pprt(2,ny-2,k,1)-pprt(2,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(pprt(2,ny-2,k,1)-pprt(2,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(pprt(2,ny-2,k,1)-pprt(2,ny-2,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          pdtnb( 1  ,k)=(cdtdx*(-pprt(2,ny-2,k,3)+pprt(1,ny-2,k,2))
     :                  +cdtdy*(pprt(2,ny-2,k,3)-pprt(2,ny-1,k,2))
     :                  -cdtdz*(pprt(2,ny-2,k+1,3)-pprt(2,ny-2,k-1,3)
     :                  )/2.)/(dtbig1*(1.-cdtdx+cdtdy))
          pdtwb(ny-1,k)=pdtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of p boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=1,nz-1

          temx =(pprt(nx-2,ny-2,k,3)+pprt(nx-2,ny-2,k,1)
     :          -2.*pprt(nx-3,ny-2,k,2))/(2.*dx)
          temy =(pprt(nx-2,ny-2,k,3)+pprt(nx-2,ny-2,k,1)
     :          -2.*pprt(nx-2,ny-3,k,2))/(2.*dy)
          temz =(pprt(nx-2,ny-2,k+1,2)-pprt(nx-2,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(pprt(nx-2,ny-2,k,1)-pprt(nx-2,ny-2,k,3))/(2.*dx)
     :          *temx/(eps+tem)
          cdtdy=(pprt(nx-2,ny-2,k,1)-pprt(nx-2,ny-2,k,3))/(2.*dy)
     :          *temy/(eps+tem)
          cdtdz=(pprt(nx-2,ny-2,k,1)-pprt(nx-2,ny-2,k,3))/(2.*dz)
     :          *temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          pdtnb(nx-1,k)=(cdtdx*(pprt(nx-2,ny-2,k,3)
     :                  -pprt(nx-1,ny-2,k,2))
     :                  +cdtdy*(pprt(nx-2,ny-2,k,3)
     :                  -pprt(nx-2,ny-1,k,2))
     :                  -cdtdz*(pprt(nx-2,ny-2,k+1,3)
     :                  -pprt(nx-2,ny-2,k-1,3))
     :                  /2.)/(dtbig1*(1.+cdtdx+cdtdy))
          pdteb(ny-1,k)=pdtnb(nx-1,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTPT3D                    ######
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 BDTPT3D(nx,ny,nz,dtbig1,
     :           ptprt,ptdteb,ptdtwb,ptdtnb,ptdtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of the potential temperature perturbation
c     on the lateral boundaries using the radiation boundary condition of 
c     Raymond and Kuo.
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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       ptprt    Perturbation potential temperature at all time levels (K)
c
c     OUTPUT:
c
c       ptdteb    Time tendency of pt at east  boundary (K/s)   
c       ptdtwb    Time tendency of pt at west  boundary (K/s)   
c       ptdtnb    Time tendency of pt at north boundary (K/s)   
c       ptdtsb    Time tendency of pt at south boundary (K/s)   
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real ptprt(nx,ny,nz,nt) ! Perturbation potential temperature (K).

      real ptdteb(ny,nz)  ! T-tendency of ptprt at e-boundary (K/s)
      real ptdtwb(ny,nz)  ! T-tendency of ptprt at w-boundary (K/s)
      real ptdtnb(nx,nz)  ! T-tendency of ptprt at n-boundary (K/s)
      real ptdtsb(nx,nz)  ! T-tendency of ptprt at s-boundary (K/s)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy, temz

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
      real cdtdz           ! Estimated phase speed times dt/dz .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of ptprt on the west boundary (ptdtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=1,nz-1
        DO 1110 j=2,ny-2

          temx =(-ptprt(2,j,k,3)-ptprt(2,j,k,1)+2.*ptprt(3,j,k,2))/
     :          (2.*dx)
          temy =(ptprt(2,j+1,k,2)-ptprt(2,j-1,k,2))/(2.*dy)
          temz =(ptprt(2,j,k+1,2)-ptprt(2,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(ptprt(2,j,k,1)-ptprt(2,j,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(ptprt(2,j,k,1)-ptprt(2,j,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(ptprt(2,j,k,1)-ptprt(2,j,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          ptdtwb(j,k)=(cdtdx*(-ptprt(2,j,k,3)+ptprt(1,j,k,2))
     :               -cdtdy*(ptprt(2,j+1,k,3)-ptprt(2,j-1,k,3))/2.
     :               -cdtdz*(ptprt(2,j,k+1,3)-ptprt(2,j,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of ptprt on the east boundary (ptdteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdx must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=1,nz-1
        DO 1210 j=2,ny-2

          temx =(ptprt(nx-2,j,k,3)+ptprt(nx-2,j,k,1)
     :          -2.*ptprt(nx-3,j,k,2))/(2.*dx)
          temy =(ptprt(nx-2,j+1,k,2)-ptprt(nx-2,j-1,k,2))/(2.*dy)
          temz =(ptprt(nx-2,j,k+1,2)-ptprt(nx-2,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(ptprt(nx-2,j,k,1)-ptprt(nx-2,j,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(ptprt(nx-2,j,k,1)-ptprt(nx-2,j,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(ptprt(nx-2,j,k,1)-ptprt(nx-2,j,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          ptdteb(j,k)=(cdtdx*(ptprt(nx-2,j,k,3)-ptprt(nx-1,j,k,2))
     :               -cdtdy*(ptprt(nx-2,j+1,k,3)-ptprt(nx-2,j-1,k,3))/2.
     :               -cdtdz*(ptprt(nx-2,j,k+1,3)-ptprt(nx-2,j,k-1,3))/2.
     :               )/(dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of ptprt on the south boundary (ptdtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1310 k=1,nz-1
        DO 1310 i=2,nx-2

          temx =(ptprt(i+1,2,k,2)-ptprt(i-1,2,k,2))/(2.*dx)
          temy =(-ptprt(i,2,k,3)-ptprt(i,2,k,1)+2.*ptprt(i,3,k,2))/
     :          (2.*dy)
          temz =(ptprt(i,2,k+1,2)-ptprt(i,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(ptprt(i,2,k,1)-ptprt(i,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(ptprt(i,2,k,1)-ptprt(i,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(ptprt(i,2,k,1)-ptprt(i,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          ptdtsb(i,k)=(-cdtdx*(ptprt(i+1,2,k,3)-ptprt(i-1,2,k,3))/2.
     :               +cdtdy*(-ptprt(i,2,k,3)+ptprt(i,1,k,2))
     :               -cdtdz*(ptprt(i,2,k+1,3)-ptprt(i,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdy))

1310    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of ptprt on the north boundary (ptdtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=1,nz-1
        DO 1410 i=2,nx-2

          temx =(ptprt(i+1,ny-2,k,2)-ptprt(i-1,ny-2,k,2))/(2.*dx)
          temy =(ptprt(i,ny-2,k,3)+ptprt(i,ny-2,k,1)
     :          -2.*ptprt(i,ny-3,k,2))/(2.*dy)
          temz =(ptprt(i,ny-2,k+1,2)-ptprt(i,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(ptprt(i,ny-2,k,1)-ptprt(i,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(ptprt(i,ny-2,k,1)-ptprt(i,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(ptprt(i,ny-2,k,1)-ptprt(i,ny-2,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          ptdtnb(i,k)=(-cdtdx*(ptprt(i+1,ny-2,k,3)
     :                -ptprt(i-1,ny-2,k,3))/2.
     :                +cdtdy*(ptprt(i,ny-2,k,3)-ptprt(i,ny-1,k,2))
     :                -cdtdz*(ptprt(i,ny-2,k+1,3)-ptprt(i,ny-2,k-1,3)
     :                )/2.)/(dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southwest corner of pt boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=1,nz-1

          temx =(-ptprt(2,2,k,3)-ptprt(2,2,k,1)+2.*ptprt(3,2,k,2))/
     :          (2.*dx)
          temy =(-ptprt(2,2,k,3)-ptprt(2,2,k,1)+2.*ptprt(2,3,k,2))/
     :          (2.*dy)
          temz =(ptprt(2,2,k+1,2)-ptprt(2,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(ptprt(2,2,k,1)-ptprt(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(ptprt(2,2,k,1)-ptprt(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(ptprt(2,2,k,1)-ptprt(2,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          ptdtsb(1,k)=(cdtdx*(-ptprt(2,2,k,3)+ptprt(1,2,k,2))
     :               +cdtdy*(-ptprt(2,2,k,3)+ptprt(2,1,k,2))
     :               -cdtdz*(ptprt(2,2,k+1,3)-ptprt(2,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          ptdtwb(1,k)=ptdtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of pt boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=1,nz-1

          temx =(ptprt(nx-2,2,k,3)+ptprt(nx-2,2,k,1)
     :          -2.*ptprt(nx-3,2,k,2))/(2.*dx)
          temy =(-ptprt(nx-2,2,k,3)-ptprt(nx-2,2,k,1)
     :          +2.*ptprt(nx-2,3,k,2))/(2.*dy)
          temz =(ptprt(nx-2,2,k+1,2)-ptprt(nx-2,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(ptprt(nx-2,2,k,1)-ptprt(nx-2,2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(ptprt(nx-2,2,k,1)-ptprt(nx-2,2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(ptprt(nx-2,2,k,1)-ptprt(nx-2,2,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          ptdtsb(nx-1,k)=(cdtdx*(ptprt(nx-2,2,k,3)-ptprt(nx-1,2,k,2))
     :                  +cdtdy*(-ptprt(nx-2,2,k,3)+ptprt(nx-2,1,k,2))
     :                  -cdtdz*(ptprt(nx-2,2,k+1,3)-ptprt(nx-2,2,k-1,3)
     :                  )/2.)/(dtbig1*(1.+cdtdx-cdtdy))
          ptdteb( 1  ,k)=ptdtsb(nx-1,k)

1520    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of pt boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=1,nz-1

          temx =(-ptprt(2,ny-2,k,3)-ptprt(2,ny-2,k,1)
     :          +2.*ptprt(3,ny-2,k,2))/(2.*dx)
          temy =(ptprt(2,ny-2,k,3)+ptprt(2,ny-2,k,1)
     :          -2.*ptprt(2,ny-3,k,2))/(2.*dy)
          temz =(ptprt(2,ny-2,k+1,2)-ptprt(2,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(ptprt(2,ny-2,k,1)-ptprt(2,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(ptprt(2,ny-2,k,1)-ptprt(2,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(ptprt(2,ny-2,k,1)-ptprt(2,ny-2,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          ptdtnb( 1  ,k)=(cdtdx*(-ptprt(2,ny-2,k,3)+ptprt(1,ny-2,k,2))
     :                  +cdtdy*(ptprt(2,ny-2,k,3)-ptprt(2,ny-1,k,2))
     :                  -cdtdz*(ptprt(2,ny-2,k+1,3)-ptprt(2,ny-2,k-1,3)
     :                  )/2.)/(dtbig1*(1.-cdtdx+cdtdy))
          ptdtwb(ny-1,k)=ptdtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of pt boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=1,nz-1

          temx =(ptprt(nx-2,ny-2,k,3)+ptprt(nx-2,ny-2,k,1)
     :          -2.*ptprt(nx-3,ny-2,k,2))/(2.*dx)
          temy =(ptprt(nx-2,ny-2,k,3)+ptprt(nx-2,ny-2,k,1)
     :          -2.*ptprt(nx-2,ny-3,k,2))/(2.*dy)
          temz =(ptprt(nx-2,ny-2,k+1,2)-ptprt(nx-2,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(ptprt(nx-2,ny-2,k,1)-ptprt(nx-2,ny-2,k,3))/(2.*dx)
     :          *temx/(eps+tem)
          cdtdy=(ptprt(nx-2,ny-2,k,1)-ptprt(nx-2,ny-2,k,3))/(2.*dy)
     :          *temy/(eps+tem)
          cdtdz=(ptprt(nx-2,ny-2,k,1)-ptprt(nx-2,ny-2,k,3))/(2.*dz)
     :          *temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          ptdtnb(nx-1,k)=(cdtdx*(ptprt(nx-2,ny-2,k,3)
     :                  -ptprt(nx-1,ny-2,k,2))
     :                  +cdtdy*(ptprt(nx-2,ny-2,k,3)
     :                  -ptprt(nx-2,ny-1,k,2))
     :                  -cdtdz*(ptprt(nx-2,ny-2,k+1,3)
     :                  -ptprt(nx-2,ny-2,k-1,3))
     :                  /2.)/(dtbig1*(1.+cdtdx+cdtdy))
          ptdteb(ny-1,k)=ptdtnb(nx-1,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BDTQ3D                     ######
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 BDTQ3D(nx,ny,nz,dtbig1, q,qdteb,qdtwb,qdtnb,qdtsb)
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the time tendencies of q on the lateral boundaries
c     using the radiation boundary condition of Raymond and Kuo.
c
c#######################################################################
c
c     AUTHOR: Ming Xue and Hao Jin
c     6/10/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       q        Mixing ratio for one of the water variables (kg/kg)
c
c     OUTPUT:
c
c       qdteb    Time tendency of q at east  boundary (kg/(kg*s))
c       qdtwb    Time tendency of q at west  boundary (kg/(kg*s))
c       qdtnb    Time tendency of q at north boundary (kg/(kg*s))
c       qdtsb    Time tendency of q at south boundary (kg/(kg*s))
c
c#######################################################################
c

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

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

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

      real dtbig1          ! The big time step size (s)

      real q(nx,ny,nz,nt)  ! Mixing ratio for one of the water
                           ! variables (kg/kg)

      real qdteb (ny,nz)   ! T-tendency of q at e-boundary (kg/(kg*s))
      real qdtwb (ny,nz)   ! T-tendency of q at w-boundary (kg/(kg*s))
      real qdtnb (nx,nz)   ! T-tendency of q at n-boundary (kg/(kg*s))
      real qdtsb (nx,nz)   ! T-tendency of q at s-boundary (kg/(kg*s))
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      real    eps, tem, temx, temy, temz

      real cdtdx           ! Estimated phase speed times dt/dx .
      real cdtdy           ! Estimated phase speed times dt/dy .
      real cdtdz           ! Estimated phase speed times dt/dz .
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'globcst.inc'     ! Global model control parameters
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      eps=1.0e-30
c
c#######################################################################
c
c     Calculate the time tendency of q on the west boundary (qdtwb)
c     eps is added to the denominator to avoid division by zero.
c     For the west boundary, cdtdx must be negative but not less than -1
c
c#######################################################################
c
      IF(wbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1110 k=2,nz-1
        DO 1110 j=2,ny-2

          temx =(-q(2,j,k,3)-q(2,j,k,1)+2.*q(3,j,k,2))/(2.*dx)
          temy =(q(2,j+1,k,2)-q(2,j-1,k,2))/(2.*dy)
          temz =(q(2,j,k+1,2)-q(2,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(q(2,j,k,1)-q(2,j,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(q(2,j,k,1)-q(2,j,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(q(2,j,k,1)-q(2,j,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          qdtwb(j,k)=(cdtdx*(-q(2,j,k,3)+q(1,j,k,2))
     :               -cdtdy*(q(2,j+1,k,3)-q(2,j-1,k,3))/2.
     :               -cdtdz*(q(2,j,k+1,3)-q(2,j,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx))

1110    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of q on the east boundary (qdteb)
c     eps is added to the denominator to avoid division by zero.
c     For the east boundary, cdtdx must be positive but not larger than 1
c
c#######################################################################
c
      IF(ebc.eq.4) THEN             ! Radiation boundary condition.

        DO 1210 k=2,nz-1
        DO 1210 j=2,ny-2

          temx =(q(nx-2,j,k,3)+q(nx-2,j,k,1)-2.*q(nx-3,j,k,2))/(2.*dx)
          temy =(q(nx-2,j+1,k,2)-q(nx-2,j-1,k,2))/(2.*dy)
          temz =(q(nx-2,j,k+1,2)-q(nx-2,j,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(q(nx-2,j,k,1)-q(nx-2,j,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(q(nx-2,j,k,1)-q(nx-2,j,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(q(nx-2,j,k,1)-q(nx-2,j,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          qdteb(j,k)=(cdtdx*(q(nx-2,j,k,3)-q(nx-1,j,k,2))
     :               -cdtdy*(q(nx-2,j+1,k,3)-q(nx-2,j-1,k,3))/2.
     :               -cdtdz*(q(nx-2,j,k+1,3)-q(nx-2,j,k-1,3))/2.)/
     :               (dtbig1*(1.+cdtdx))

1210    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of q on the south boundary (qdtsb)
c     eps is added to the denominator to avoid division by zero.
c     For the south boundary, cdtdy must be negative but not less than -1
c
c#######################################################################
c
      IF(sbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1310 k=2,nz-1
        DO 1310 i=2,nx-2

          temx =(q(i+1,2,k,2)-q(i-1,2,k,2))/(2.*dx)
          temy =(-q(i,2,k,3)-q(i,2,k,1)+2.*q(i,3,k,2))/(2.*dy)
          temz =(q(i,2,k+1,2)-q(i,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(q(i,2,k,1)-q(i,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(q(i,2,k,1)-q(i,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(q(i,2,k,1)-q(i,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          qdtsb(i,k)=(-cdtdx*(q(i+1,2,k,3)-q(i-1,2,k,3))/2.
     :               +cdtdy*(-q(i,2,k,3)+q(i,1,k,2))
     :               -cdtdz*(q(i,2,k+1,3)-q(i,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdy))

1310    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency of q on the north boundary (qdtnb)
c     eps is added to the denominator to avoid division by zero.
c     For the north boundary, cdtdy must be positive but not larger than 1
c
c#######################################################################
c
      IF(nbc.eq.4) THEN             ! Radiation boundary condition.

        DO 1410 k=2,nz-1
        DO 1410 i=2,nx-2

          temx =(q(i+1,ny-2,k,2)-q(i-1,ny-2,k,2))/(2.*dx)
          temy =(q(i,ny-2,k,3)+q(i,ny-2,k,1)-2.*q(i,ny-3,k,2))/(2.*dy)
          temz =(q(i,ny-2,k+1,2)-q(i,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(q(i,ny-2,k,1)-q(i,ny-2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(q(i,ny-2,k,1)-q(i,ny-2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(q(i,ny-2,k,1)-q(i,ny-2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,-1.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          qdtnb(i,k)=(-cdtdx*(q(i+1,ny-2,k,3)-q(i-1,ny-2,k,3))/2.
     :               +cdtdy*(q(i,ny-2,k,3)-q(i,ny-1,k,2))
     :               -cdtdz*(q(i,ny-2,k+1,3)-q(i,ny-2,k-1,3))/2.)/
     :               (dtbig1*(1.+cdtdy))

1410    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southwest corner of q boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.wbc.eq.4) THEN

        DO 1510 k=2,nz-1

          temx =(-q(2,2,k,3)-q(2,2,k,1)+2.*q(3,2,k,2))/(2.*dx)
          temy =(-q(2,2,k,3)-q(2,2,k,1)+2.*q(2,3,k,2))/(2.*dy)
          temz =(q(2,2,k+1,2)-q(2,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(q(2,2,k,1)-q(2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(q(2,2,k,1)-q(2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(q(2,2,k,1)-q(2,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          qdtsb(1,k)=(cdtdx*(-q(2,2,k,3)+q(1,2,k,2))
     :               +cdtdy*(-q(2,2,k,3)+q(2,1,k,2))
     :               -cdtdz*(q(2,2,k+1,3)-q(2,2,k-1,3))/2.)/
     :               (dtbig1*(1.-cdtdx-cdtdy))
          qdtwb(1,k)=qdtsb(1,k)

1510    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the southeast corner of q boundary
c
c#######################################################################
c
      IF(sbc.eq.4.and.ebc.eq.4) THEN

        DO 1520 k=2,nz-1

          temx =(q(nx-2,2,k,3)+q(nx-2,2,k,1)-2.*q(nx-3,2,k,2))/(2.*dx)
          temy =(-q(nx-2,2,k,3)-q(nx-2,2,k,1)+2.*q(nx-2,3,k,2))/(2.*dy)
          temz =(q(nx-2,2,k+1,2)-q(nx-2,2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(q(nx-2,2,k,1)-q(nx-2,2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(q(nx-2,2,k,1)-q(nx-2,2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(q(nx-2,2,k,1)-q(nx-2,2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=min(cdtdy,0.0)
          cdtdy=max(cdtdy,-1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          qdtsb(nx-1,k)=(cdtdx*(q(nx-2,2,k,3)-q(nx-1,2,k,2))
     :                  +cdtdy*(-q(nx-2,2,k,3)+q(nx-2,1,k,2))
     :                  -cdtdz*(q(nx-2,2,k+1,3)-q(nx-2,2,k-1,3))/2.)/
     :                  (dtbig1*(1.+cdtdx-cdtdy))
          qdteb( 1  ,k)=qdtsb(nx-1,k)

1520    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northwest corner of q boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.wbc.eq.4) THEN

        DO 1530 k=2,nz-1

          temx =(-q(2,ny-2,k,3)-q(2,ny-2,k,1)+2.*q(3,ny-2,k,2))/(2.*dx)
          temy =(q(2,ny-2,k,3)+q(2,ny-2,k,1)-2.*q(2,ny-3,k,2))/(2.*dy)
          temz =(q(2,ny-2,k+1,2)-q(2,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(q(2,ny-2,k,1)-q(2,ny-2,k,3))/(2.*dx)*temx/(eps+tem)
          cdtdy=(q(2,ny-2,k,1)-q(2,ny-2,k,3))/(2.*dy)*temy/(eps+tem)
          cdtdz=(q(2,ny-2,k,1)-q(2,ny-2,k,3))/(2.*dz)*temz/(eps+tem)
          cdtdx=min(cdtdx,0.0)
          cdtdx=max(cdtdx,-1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          qdtnb( 1  ,k)=(cdtdx*(-q(2,ny-2,k,3)+q(1,ny-2,k,2))
     :                  +cdtdy*(q(2,ny-2,k,3)-q(2,ny-1,k,2))
     :                  -cdtdz*(q(2,ny-2,k+1,3)-q(2,ny-2,k-1,3))/2.)/
     :                  (dtbig1*(1.-cdtdx+cdtdy))
          qdtwb(ny-1,k)=qdtnb(1,k)

1530    CONTINUE

      ENDIF

c
c#######################################################################
c
c     Calculate the time tendency on the northeast corner of q boundary
c
c#######################################################################
c
      IF(nbc.eq.4.and.ebc.eq.4) THEN

        DO 1540 k=2,nz-1

          temx =(q(nx-2,ny-2,k,3)+q(nx-2,ny-2,k,1)
     :          -2.*q(nx-3,ny-2,k,2))/(2.*dx)
          temy =(q(nx-2,ny-2,k,3)+q(nx-2,ny-2,k,1)
     :          -2.*q(nx-2,ny-3,k,2))/(2.*dy)
          temz =(q(nx-2,ny-2,k+1,2)-q(nx-2,ny-2,k-1,2))/(2.*dz)
          tem  =temx**2+temy**2+temz**2
          cdtdx=(q(nx-2,ny-2,k,1)-q(nx-2,ny-2,k,3))/
     :          (2.*dx)*temx/(eps+tem)
          cdtdy=(q(nx-2,ny-2,k,1)-q(nx-2,ny-2,k,3))/
     :          (2.*dy)*temy/(eps+tem)
          cdtdz=(q(nx-2,ny-2,k,1)-q(nx-2,ny-2,k,3))/
     :          (2.*dz)*temz/(eps+tem)
          cdtdx=max(cdtdx,0.0)
          cdtdx=min(cdtdx,1.0)
          cdtdy=max(cdtdy,0.0)
          cdtdy=min(cdtdy,1.0)
          cdtdz=max(cdtdz,-1.0)
          cdtdz=min(cdtdz,1.0)

          qdtnb(nx-1,k)=(cdtdx*(q(nx-2,ny-2,k,3)-q(nx-1,ny-2,k,2))
     :                  +cdtdy*(q(nx-2,ny-2,k,3)-q(nx-2,ny-1,k,2))
     :                  -cdtdz*(q(nx-2,ny-2,k+1,3)-q(nx-2,ny-2,k-1,3))
     :                  /2.)/(dtbig1*(1.+cdtdx+cdtdy))
          qdteb(ny-1,k)=qdtnb(nx-1,k)

1540    CONTINUE

      ENDIF

      RETURN
      END

