c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE SOLVPT                     ######
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 SOLVPT(nx,ny,nz, dtbig1, 
     :           ptprt,u,v,w, ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :           ptbar,rhobar,km,ptsflx,
     :           x,y,z,zp, j1,j2,j3,
     :           ptadv,ptmix,ptsrc,tem1,tem2,tem3,tem4,tem5,tem6)
c
C#######################################################################
c
c     PURPOSE:
c
c     Coordinate the time integration of the potential temperature 
c     equation.
c
C#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (M. Xue)
c     Added full documentation.
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       dtbig1   The big time step size (s)
c
c       ptprt    Perturbation potential temperature at times tpast and
c                tpresent (K)
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        Vertical component of Cartesian velocity at all t-levels (m/s)
c
c       ptdteb   Time tendency of the ptprt field at the east boundary
c       ptdtwb   Time tendency of the ptprt field at the west boundary
c       ptdtnb   Time tendency of the ptprt field at the north boundary
c       ptdtsb   Time tendency of the ptprt field at the south boundary
c
c       ptbar    Base state potential temperature (K)
c       rhobar   Base state density (kg/m**3)
c
c       km       Turbulent mixing coefficient for momentum.
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c
c       j1       Coordinate transform Jacobian -d(zp)/dx
c       j2       Coordinate transform Jacobian -d(zp)/dy
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c
c     OUTPUT:
c
c       ptprt    Perturbation potential temperature at time tfuture (K)
c
c     WORK ARRAYS:
c
c       ptadv    Advection term of potential temperature eq. (kg*K/(m**3*s)) 
c                A local array.
c       ptmix    Total mixing in potential temperature eq. (kg*K/(m**3*s))
c                A local array.
c       ptsrc    Other source/sink terms of potential temperature eq.  
c                (kg*K/(m**3*s)). A local array.
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c       tem3     Temporary work array.
c       tem4     Temporary work array.
c       tem5     Temporary work array.
c       tem6     Temporary work array.
c          
C#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none             ! Force explicit declarations
c
      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        ! Number of grid points in 3 directions

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

      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 ptdteb(ny,nz)        ! Time tendency of ptprt at east boundary
      real ptdtwb(ny,nz)        ! Time tendency of ptprt at west boundary
      real ptdtnb(nx,nz)        ! Time tendency of ptprt at north boundary
      real ptdtsb(nx,nz)        ! Time tendency of ptprt at south boundary

      real ptbar (nx,ny,nz)     ! Base state potential temperature (K)
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)

      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )
      real ptsflx(nx,ny)        ! Surface flux of heat (K*kg/(m*s**2))

      real x     (nx)           ! The x-coord. of the physical and
                                ! computational grid. Defined at u-point.
      real y     (ny)           ! The y-coord. of the physical and
                                ! computational grid. Defined at v-point.
      real z     (nz)           ! The z-coord. of the computational grid.
                                ! Defined at w-point on the staggered grid.
      real zp    (nx,ny,nz)     ! The physical height coordinate defined at
                                ! w-point of the staggered grid.

      real j1    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( x ).
      real j2    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( y ).
      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! d( zp )/d( z ).
c
      real ptadv (nx,ny,nz)     ! Advection term of potential temperature
                                ! (kg*K/(m**3*s)). A local array.
      real ptmix (nx,ny,nz)     ! Total mixing in potential temperature
                                ! (kg*K/(m**3*s)). A local array.
      real ptsrc (nx,ny,nz)     ! Other source/sink terms of potential 
                                ! temperature (kg*K/(m**3*s)). A local array.

      real tem1  (nx,ny,nz)     ! Temporary work array
      real tem2  (nx,ny,nz)     ! Temporary work array
      real tem3  (nx,ny,nz)     ! Temporary work array
      real tem4  (nx,ny,nz)     ! Temporary work array
      real tem5  (nx,ny,nz)     ! Temporary work array
      real tem6  (nx,ny,nz)     ! Temporary work array
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c

      real dtbig1               ! Local big time step
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external advpt
      external mixpt
      external srcpt
      external steppt
      external bcpt
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Compute the advection term of the potential temperature
c     equation and store it in array ptadv.
c
c#######################################################################
c
      CALL advpt(nx,ny,nz,ptprt,u,v,w, rhobar,ptbar,
     :           ptadv,
     :           tem1,tem2,tem3,tem4,tem5)
c
c#######################################################################
c
c     Compute the mixing term of the potential temperature equation.  
c     This includes both physical and computational mixing.  
c     Store it in array ptmix.
c
c#######################################################################
c
      CALL mixpt(nx,ny,nz,
     :           ptprt(1,1,1,1),ptbar,rhobar,km,ptsflx,
     :           x,y,z,zp,j1,j2,j3,
     :           ptmix,
     :           tem1,tem2,tem3,tem4,tem5,tem6)
c
c#######################################################################
c
c     Compute the source and sink terms of the potential temperature
c     equation, including forcings for phase change and radiation
c     Store it in array ptsrc.
c
c#######################################################################

      CALL srcpt(nx,ny,nz,ptprt,ptbar,ptsrc)
c
c#######################################################################
c
c     Integrate forward by one timestep the potential temperature
c     equation, yielding ptprt at tfuture.
c
c#######################################################################
c
      CALL steppt(nx,ny,nz, dtbig1,
     :            ptprt,rhobar,ptadv,ptmix,ptsrc)
c
c#######################################################################
c
c     Set the boundary conditions for the potential temperature
c
c#######################################################################
c
      CALL bcpt(nx,ny,nz, dtbig1, ptprt, ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :          ebc,wbc,nbc,sbc,tbc,bbc)

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE SOLVQV                     ######
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
c
      SUBROUTINE SOLVQV(nx,ny,nz, dtbig1,
     :           qv,u,v,w, qvdteb,qvdtwb,qvdtnb,qvdtsb,
     :           rhobar,qvbar, km, qvsflx,
     :           x,y,z,zp, j1,j2,j3,
     :           qadv,qmix,tem1,tem2,tem3,tem4,tem5,tem6)

c
c########################################################################
c
c     PURPOSE:
c
c     Coordinate the time integration of the equation for water vapor 
c     specific humidity qv, and the setting of boundary conditions.
c
c########################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91  
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (M. Xue)
c     Added full documentation.
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       dtbig1   The big time step size (s)
c
c       qv       Water vapor specific humidity at tpast and tpresent (kg/kg) 
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        Vertical component of Cartesian velocity at all time 
c                levels (m/s)
c
c       qvdteb   Time tendency of the qv field at the east boundary
c       qvdtwb   Time tendency of the qv field at the west boundary
c       qvdtnb   Time tendency of the qv field at the north boundary
c       qvdtsb   Time tendency of the qv field at the south boundary
c
c       rhobar   Base state density (kg/m**3)
c       qvbar    Base state water vapor specific humidity (kg/kg)
c
c       km       Turbulent mixing coefficient for momentum.
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c
c       j1       Coordinate transform Jacobian -d(zp)/dx
c       j2       Coordinate transform Jacobian -d(zp)/dy
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c
c     OUTPUT:
c
c       qv       Water vapor specific humidity at time tfuture (kg/kg) 
c
c     WORK ARRAYS:
c
c       qadv     Advection term of water vapor eq. (kg/(m**3*s)).
c                A local array.
c       qmix     Total mixing in water vapor eq. (kg/(m**3*s)). 
c                A local array.
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c       tem3     Temporary work array.
c       tem4     Temporary work array.
c       tem5     Temporary work array.
c       tem6     Temporary work array.
c
C#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none             ! Force explicit declarations
c
      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        ! Number of grid points in 3 directions

      real qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg)
      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)
c
      real qvdteb(ny,nz)        ! Time tendency of qv at east boundary
      real qvdtwb(ny,nz)        ! Time tendency of qv at west boundary
      real qvdtnb(nx,nz)        ! Time tendency of qv at north boundary
      real qvdtsb(nx,nz)        ! Time tendency of qv at south boundary
c
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                                ! (kg/kg)
c
      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )
      real qvsflx(nx,ny)        ! Surface flux of moisture (kg/(m**2*s))
c
      real x     (nx)           ! The x-coord. of the physical and
                                ! computational grid. Defined at u-point.
      real y     (ny)           ! The y-coord. of the physical and
                                ! computational grid. Defined at v-point.
      real z     (nz)           ! The z-coord. of the computational grid.
                                ! Defined at w-point on the staggered grid.
      real zp    (nx,ny,nz)     ! The physical height coordinate defined at
                                ! w-point of the staggered grid.
c
      real j1    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( x ).
      real j2    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( y ).
      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! d( zp )/d( z ).
c
      real qadv(nx,ny,nz)       ! Advection of water vapor (kg/(m**3*s))
                                ! A local array.
      real qmix(nx,ny,nz)       ! Total mixing of water vapor (kg/(m**3*s))
                                ! A local array.
      real tem1  (nx,ny,nz)     ! Temporary work array
      real tem2  (nx,ny,nz)     ! Temporary work array
      real tem3  (nx,ny,nz)     ! Temporary work array
      real tem4  (nx,ny,nz)     ! Temporary work array
      real tem5  (nx,ny,nz)     ! Temporary work array
      real tem6  (nx,ny,nz)     ! Temporary work array
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      real dtbig1               ! Local big time step
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external advq
      external mixqv
      external stepq
      external bcq
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Compute the advection term of the water vapor specific humidity
c     equation and store the result in array qadv.
c
c#######################################################################
c
      CALL advq(nx,ny,nz,qv,u,v,w,rhobar,
     :          qadv, 
     :          tem1,tem2,tem3,tem4,tem5)
c
c#######################################################################
c
c     Compute the mixing terms for the water vapor specific humidity
c     equation. This includes both physical and computational
c     mixing.  Store the result in array qmix.
c
c#######################################################################
c
      CALL mixqv(nx,ny,nz,
     :           qv(1,1,1,1),rhobar,qvbar,km,qvsflx, 
     :           x,y,z,zp, j1,j2,j3,
     :           qmix,
     :           tem1,tem2,tem3,tem4,tem5,tem6)

c
c#######################################################################
c
c     Integrate forward by one timestep the water vapor specific humidity
c     equation, yielding qv at tfuture.
c
c#######################################################################
c
      CALL stepq(nx,ny,nz, dtbig1, qv, qadv,qmix, rhobar)
c
c#######################################################################
c
c     Set the boundary conditions for the water vapor specific humidity
c
c#######################################################################

      CALL bcq(nx,ny,nz, dtbig1, qv, qvdteb,qvdtwb,qvdtnb,qvdtsb,
     :         ebc,wbc,nbc,sbc,tbc,bbc)

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE SOLVQ                      ######
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 SOLVQ(nx,ny,nz, dtbig1,
     :           q,u,v,w, qdteb,qdtwb,qdtnb,qdtsb, rhobar,km,
     :           x,y,z,zp, j1,j2,j3,
     :           qadv,qmix,tem1,tem2,tem3,tem4,tem5,tem6)
c
c########################################################################
c
c     PURPOSE:
c
c     Coordinate the time integration of the equations  for the water substance 
c     quantities qc,qr,qi,qs and qh, i.e. all water variables except water 
c     vapor specific humidity qv.
c
c########################################################################
c
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (K. Droegemeier and M. Xue)
c     Added full documentation.
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       dtbig1   The big time step size (s)
c
c       q        One of the liquid or ice variables at tpast and 
c                tpresent (kg/kg) 
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        Vertical component of Cartesian velocity at all 
c                time levels (m/s)
c
c       qdteb    Time tendency of liquid/ice variable q at east boundary
c       qdtwb    Time tendency of liquid/ice variable q at west boundary
c       qdtnb    Time tendency of liquid/ice variable q at north boundary
c       qdtsb    Time tendency of liquid/ice variable q at south boundary
c
c       rhobar   Base state density (kg/m**3)
c
c       km       Turbulent mixing coefficient for momentum.
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c
c       j1       Coordinate transform Jacobian -d(zp)/dx
c       j2       Coordinate transform Jacobian -d(zp)/dy
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c
c     OUTPUT:
c
c       q        Liquid/ice variable q at tfuture (kg/kg) 
c
c     WORK ARRAYS:
c
c       qadv     Advection term of water variable eq. (kg/(m**3*s)).
c                A local array.
c       qmix     Total mixing in water variable eq. (kg/(m**3*s)).
c                A local array.
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c       tem3     Temporary work array.
c       tem4     Temporary work array.
c       tem5     Temporary work array.
c       tem6     Temporary work array.
c
C#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none             ! Force explicit declarations
c
      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        ! Number of grid points in 3 directions

      real q     (nx,ny,nz,nt)  ! One of the water/ice variables (kg/kg) 
      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)
c
      real qdteb(ny,nz)         ! Time tendency of liquid/ice at e-boundary
      real qdtwb(ny,nz)         ! Time tendency of liquid/ice at w-boundary
      real qdtnb(nx,nz)         ! Time tendency of liquid/ice at n-boundary
      real qdtsb(nx,nz)         ! Time tendency of liquid/ice at s-boundary

      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )
c
      real x     (nx)           ! The x-coord. of the physical and
                                ! computational grid. Defined at u-point.
      real y     (ny)           ! The y-coord. of the physical and
                                ! computational grid. Defined at v-point.
      real z     (nz)           ! The z-coord. of the computational grid.
                                ! Defined at w-point on the staggered grid.
      real zp    (nx,ny,nz)     ! The physical height coordinate defined at
                                ! w-point of the staggered grid.
c
      real j1    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( x ).
      real j2    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( y ).
      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! d( zp )/d( z ).
c
      real qadv(nx,ny,nz)       ! Advection of water/ice substance 
      				  ! (kg/(m**3*s)). A local array.
      real qmix(nx,ny,nz)       ! Total mixing of water/ice substance 
                                ! (kg/(m**3*s)). A local array.
      real tem1  (nx,ny,nz)     ! Temporary work array
      real tem2  (nx,ny,nz)     ! Temporary work array
      real tem3  (nx,ny,nz)     ! Temporary work array
      real tem4  (nx,ny,nz)     ! Temporary work array
      real tem5  (nx,ny,nz)     ! Temporary work array
      real tem6  (nx,ny,nz)     ! Temporary work array
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      real dtbig1               ! Local big time step
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external advcts
      external mixq
      external stepq
      external bcq
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Compute the advection term in the equation for a water substance 
c     q and store the result in qadv.
c
c#######################################################################
c
      CALL advq(nx,ny,nz,q,u,v,w,rhobar,
     :          qadv, 
     :          tem1,tem2,tem3,tem4,tem5)
c
c#######################################################################
c
c     Compute the mixing terms for the q equation, including both
c     physical and computational mixing.  Store the result in
c     array qmix.  
c
c#######################################################################
c
      CALL mixq(nx,ny,nz,
     :          q(1,1,1,1),rhobar,km,
     :          x,y,z,zp, j1,j2,j3,
     :          qmix,
     :          tem1,tem2,tem3,tem4,tem5,tem6)

c
c#######################################################################
c
c     Integrate forward by one timestep the q equation, yielding 
c     q at tfuture. 
c
c#######################################################################
c
      CALL stepq(nx,ny,nz, dtbig1, q, qadv,qmix, rhobar)
c
c#######################################################################
c
c     Set the boundary conditions for q 
c
c#######################################################################

      CALL bcq(nx,ny,nz, dtbig1, q, qdteb,qdtwb,qdtnb,qdtsb,
     :         ebc,wbc,nbc,sbc,tbc,bbc)

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE SOLVUVW                    ######
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 SOLVUVW(nx,ny,nz, dtsml1,
     :           u,v,w,pprt, 
     :           udteb,udtwb,udtnb,udtsb,
     :           vdteb,vdtwb,vdtnb,vdtsb,wdteb,wdtwb,wdtnb,wdtsb,
     :           rhobar, 
     :           uforce,vforce,wforce, 
     :           x,y,z,zp, j1,j2,j3,
     :           upgrad,vpgrad,wpgrad, 
     :           urho,vrho,wrho,div,tem1)
c
C#######################################################################
c
c     PURPOSE:
c
c     Coordinate the time stepping of the three momentum equations and the 
c     setting of boundary conditions. Divergence damping is also handled 
c     by this routine. The acoustically inactive forcing terms in these 
c     equations have been evaluated prior to this subroutine and are passed 
c     into this routine.  
c
C#######################################################################
c
c
c     AUTHOR: Ming Xue
c     10/10/92
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (K. Droegemeier and M. Xue)
c     Added full documentation.
c
c     6/07/92 ( M. Xue)
c     Now only tfuture fields of u, v, w and pprt are passed in.
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 tfuture (m/s)
c       v        y component of velocity at tfuture (m/s)
c       w        Vertical component of Cartesian velocity at tfuture (m/s)
c
c       pprt     Perturbation pressure at tfuture (Pascal)
c
c       udteb    Time tendency of the u field at the east boundary
c       udtwb    Time tendency of the u field at the west boundary
c       udtnb    Time tendency of the u field at the north boundary
c       udtsb    Time tendency of the u field at the south boundary
c
c       vdteb    Time tendency of the v field at the east boundary
c       vdtwb    Time tendency of the v field at the west boundary
c       vdtnb    Time tendency of the v field at the north boundary
c       vdtsb    Time tendency of the v field at the south boundary
c
c       wdteb    Time tendency of the w field at the east boundary
c       wdtwb    Time tendency of the w field at the west boundary
c       wdtnb    Time tendency of the w field at the north boundary
c       wdtsb    Time tendency of the w field at the south boundary
c
c       rhobar   Base state density (kg/m**3)
c
c       uforce   Acoustically inactive forcing terms in u-eq.(kg/(m*s)**2)
c       vforce   Acoustically inactive forcing terms in v-eq.(kg/(m*s)**2)
c       wforce   Acoustically inactive forcing terms in w-eq.(kg/(m*s)**2)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c
c       j1       Coordinate transform Jacobian -d(zp)/dx
c       j2       Coordinate transform Jacobian -d(zp)/dy
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c
c     OUTPUT:
c
c       u        x component of velocity at time tfuture updated 
c                in time by dtsml1 (m/s)
c       v        y component of velocity at time tfuture updated 
c                in time by dtsml1 (m/s)
c       w        Vertical component of Cartesian velocity at tfuture 
c                updated in time by dtsml1 (m/s)
c
c     WORK ARRAYS:
c
c       upgrad   Pressure gradient term in u-eq. (kg/(m*s)**2)
c       vpgrad   Pressure gradient term in v-eq. (kg/(m*s)**2)
c       wpgrad   Pressure gradient term in w-eq. (kg/(m*s)**2)
c       urho     u * rhobar
c       vrho     v * rhobar
c       wrho     w * rhobar
c       div      Temporary array
c       tem1     Temporary array
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none             ! Force explicit declarations
c
      integer nx, ny, nz        ! Number of grid points in 3 directions

      real u     (nx,ny,nz)     ! u-velocity at tfuture (m/s)
      real v     (nx,ny,nz)     ! v-velocity at tfuture (m/s)
      real w     (nx,ny,nz)     ! w-velocity at tfuture (m/s)
c
      real pprt  (nx,ny,nz)     ! Perturbation pressure at tfuture (Pascal)
c
      real udteb (ny,nz)        ! Time tendency of u at east boundary
      real udtwb (ny,nz)        ! Time tendency of u at west boundary
      real udtnb (nx,nz)        ! Time tendency of u at north boundary
      real udtsb (nx,nz)        ! Time tendency of u at south boundary

      real vdteb (ny,nz)        ! Time tendency of v at east boundary
      real vdtwb (ny,nz)        ! Time tendency of v at west boundary
      real vdtnb (nx,nz)        ! Time tendency of v at north boundary
      real vdtsb (nx,nz)        ! Time tendency of v at south boundary

      real wdteb (ny,nz)        ! Time tendency of w at east boundary
      real wdtwb (ny,nz)        ! Time tendency of w at west boundary
      real wdtnb (nx,nz)        ! Time tendency of w at north boundary
      real wdtsb (nx,nz)        ! Time tendency of w at south boundary
c
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real uforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                                ! in u-momentum equation (kg/(m*s)**2)
      real vforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                                ! in v-momentum equation (kg/(m*s)**2)
      real wforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                                ! in w-momentum equation (kg/(m*s)**2)
c
      real x     (nx)           ! The x-coord. of the physical and
                                ! computational grid. Defined at u-point.
      real y     (ny)           ! The y-coord. of the physical and
                                ! computational grid. Defined at v-point.
      real z     (nz)           ! The z-coord. of the computational grid.
                                ! Defined at w-point on the staggered grid.
      real zp    (nx,ny,nz)     ! The physical height coordinate defined at
                                ! w-point of the staggered grid.

      real j1    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( x ).
      real j2    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( y ).
      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! d( zp )/d( z ).
c
c#######################################################################
c
c     Temporary WORK ARRAYS: 
c
c#######################################################################
c
      real upgrad(nx,ny,nz)     ! Pressure gradient term in u-eq. 
      real vpgrad(nx,ny,nz)     ! Pressure gradient term in v-eq.
      real wpgrad(nx,ny,nz)     ! Pressure gradient term in w-eq.
      real urho  (nx,ny,nz)     ! u * rhobar
      real vrho  (nx,ny,nz)     ! v * rhobar
      real wrho  (nx,ny,nz)     ! w * rhobar
      real div   (nx,ny,nz)     ! Temporary array
      real tem1  (nx,ny,nz)     ! Temporary array

c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      real    dtsml1            ! Local small time step
      integer i, j, k
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'bndry.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external uvwrho
      external pgrad
      external stepu
      external stepv
      external stepw
      external bcu
      external bcv
      external bcw
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

C#######################################################################
c
c     If divergence damping is activated, divdmp = 1.  
c     Otherwise, this is skipped to save computations.
c
C#######################################################################


      IF( divdmp.eq.1) THEN 
c
C#######################################################################
c
c     The divergence damping terms are defined as 
c
c       d(cdvdmp*div)/dx for u eq., 
c       d(cdvdmp*div)/dy for v eq., 
c       d(cdvdmp*div)/dz for w eq., 
c
c     where 
c
c       div = difx(u*rhobar) + dify(v*rhobar) + difz(w*rhobar)
c
c     For convenience, the div. damping terms are grouped together 
c     with the pressure gradient terms, i.e. 
c
c       upgrad = d( pprt - cdvdmp*div )/dz
c       vpgrad = d( pprt - cdvdmp*div )/dy
c       wpgrad = d( pprt - cdvdmp*div )/dz
c
C#######################################################################
c

c
C#######################################################################
c
c     First compute arrays rhobar*u, rhobar*v and rhobar*w 
c     and store them in urho, vrho and wrho.
c
C#######################################################################

        CALL uvwrho(nx,ny,nz,u,v,w, rhobar, urho,vrho,wrho, tem1)
c
C#######################################################################
c
c     Compute the momentum divergence term, defined as
c
c     div = d(u*rhobar)/dx + d(v*rhobar)/dy + d(w*rhobar)/dz.
c
c     and combine the pressure and the diverence dmaping part into one
c     array so that pressure gradient force and the diverence damping 
c     can be calcualted in a single step.
c
C#######################################################################
c
        DO 100 k=1,nz-1
        DO 100 j=1,ny-1
        DO 100 i=1,nx-1
          div(i,j,k)=(urho(i+1,j,k)-urho(i,j,k))*dxinv
     :              +(vrho(i,j+1,k)-vrho(i,j,k))*dyinv
     :              +(wrho(i,j,k+1)-wrho(i,j,k))*dzinv
100     CONTINUE

        DO 200 k=1,nz-1
        DO 200 j=1,ny-1
        DO 200 i=1,nx-1
          tem1(i,j,k)= pprt(i,j,k) - cdvdmp*  div(i,j,k)
200     CONTINUE
 
      ELSE

        DO 300 k=1,nz-1
        DO 300 j=1,ny-1
        DO 300 i=1,nx-1
          tem1(i,j,k)= pprt(i,j,k)
300     CONTINUE

      ENDIF
c
C#######################################################################
c
c     Compute the pressure gradient force terms for the three
c     momentum equations and store the values in upgrad, vpgrad,
c     and wpgrad. 
c
c     When divergence damping is active, these arrays also contain
c     the divergence damping.
c
C#######################################################################
c
      CALL pgrad(nx,ny,nz, tem1, upgrad,vpgrad,wpgrad)
c
c#######################################################################
c
c     Integrate the u, v, and w equations forward by one small
c     timestep dtsml1 using a forward-in-time integration scheme
c     (that is, forward relative to the pressure gradient terms).
c
c#######################################################################
c
      CALL stepu(nx,ny,nz, dtsml1, u, uforce,upgrad,rhobar)

      CALL stepv(nx,ny,nz, dtsml1, v, vforce,vpgrad,rhobar)

      CALL stepw(nx,ny,nz, dtsml1, w, wforce,wpgrad,rhobar, tem1)
c
c#######################################################################
c
c     Set the boundary conditions for the three momentum equations.
c
c#######################################################################
c
      CALL bcu(nx,ny,nz, dtsml1, u, udteb,udtwb,udtnb,udtsb,
     :         ebc,wbc,nbc,sbc,tbc,bbc)

      CALL bcv(nx,ny,nz, dtsml1, v, vdteb,vdtwb,vdtnb,vdtsb,
     :         ebc,wbc,nbc,sbc,tbc,bbc)

      CALL bcw(nx,ny,nz, dtsml1, w, wdteb,wdtwb,wdtnb,wdtsb,
     :         ebc,wbc,nbc,sbc,tbc,bbc)

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE SOLVP                      ######
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 SOLVP(nx,ny,nz,dtsml1, 
     :           pprt,u,v,w, pdteb, pdtwb, pdtnb, pdtsb,
     :           rhobar,pbar, pforce,
     :           x,y,z,zp, j1,j2,j3,
     :           div,pdiv,tem1,tem2)
c
c#######################################################################
c
c     PURPOSE:
c
c     Coordinates the time stepping of the pressure equation, 
c     and the application of the boundary conditions.
c 
c     The acoustically inactive terms in this equation have been 
c     previously evaluated and are passed into this routine. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (K. Droegemeier and M. Xue)
c     Added full documentation.
c
c     6/07/92 ( M. Xue)
c     Now only tfuture fields of pprt, u, v and w are passed in.
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       dtsml1   The local value of small time step size
c
c       pprt     Perturbation pressure at tfuture (Pascal)
c
c       u        x component of velocity at tfuture (m/s)
c       v        y component of velocity at tfuture (m/s)
c       w        z component of velocity at tfuture (m/s)
c
c       pdteb    Time tendency of the pressure field at the east boundary
c       pdtwb    Time tendency of the pressure field at the west boundary
c       pdtnb    Time tendency of the pressure field at the north boundary
c       pdtsb    Time tendency of the pressure field at the south boundary
c
c       rhobar   Base state density (kg/m**3)
c       pbar     Base state pressure (Pascal)
c
c       pforce   Acoustically inactive terms in pressure eq. (Pascal/s)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c
c       j1       Coordinate transform Jacobian -d(zp)/dx
c       j2       Coordinate transform Jacobian -d(zp)/dy
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c
c     OUTPUT:
c
c       pprt     Perturbation pressure at tfuture updated in time 
c                by stsml1 (Pascal)
c
c     WORK ARRAYS:
c
c       div      Velocity devergence (1/s), a local array
c       pdiv     Divergence term in pressure eq. (Pascal/s), a local array
c       tem1     Temporary work array
c       tem2     Temporary work array
c
C#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none             ! Force explicit declarations
c
      integer nx, ny, nz        ! Number of grid points in 3 directions

      real dtsml1               ! Local small time step size (s)

      real pprt  (nx,ny,nz)     ! Perturbation pressure at tfuture (Pascal)
c
      real u     (nx,ny,nz)     ! u-velocity at tfuture (m/s)
      real v     (nx,ny,nz)     ! v-velocity at tfuture (m/s)
      real w     (nx,ny,nz)     ! w-velocity at tfuture (m/s)
c
      real pdteb (ny,nz)        ! Time tendency of pressure at east boundary
      real pdtwb (ny,nz)        ! Time tendency of pressure at west boundary
      real pdtnb (nx,nz)        ! Time tendency of pressure at north boundary
      real pdtsb (nx,nz)        ! Time tendency of pressure at south boundary
c
      real pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real pforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                                ! in the pressure equation (Pascal/s)
c
      real x     (nx)           ! The x-coord. of the physical and
                                ! computational grid. Defined at u-point.
      real y     (ny)           ! The y-coord. of the physical and
                                ! computational grid. Defined at v-point.
      real z     (nz)           ! The z-coord. of the computational grid.
                                ! Defined at w-point on the staggered grid.
      real zp    (nx,ny,nz)     ! The physical height coordinate defined at
                                ! w-point of the staggered grid.

      real j1    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( x ).
      real j2    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( y ).
      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! d( zp )/d( z ).
c
c#######################################################################
c
c     Temporary work array: 
c
c#######################################################################
c
      real div   (nx,ny,nz)     ! Velocity devergence (1/s)
      real pdiv  (nx,ny,nz)     ! Divergence term in pressure eq. (Pascal/s)
      real tem1  (nx,ny,nz)     ! Temporary work array
      real tem2  (nx,ny,nz)     ! Temporary work array
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'bndry.inc'

c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external divgs
      external pdivrg
      external stepp
      external bcp
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Calculate the velocity divergence using newly updated velocity.
c
c#######################################################################
c

      CALL divgs(nx,ny,nz,u,v,w, div, tem1, tem2)

c
c#######################################################################
c
c     Compute the divergence term in the pressure equation.
c
c#######################################################################
c

      CALL pdivrg(nx,ny,nz, pbar,rhobar,div, pdiv)

c
c#######################################################################
c
c     Integrate forward the pressure equation by one small timestep, 
c     using a backard (relative to the divergence term) integration 
c     scheme.
c
c#######################################################################
c

      CALL stepp(nx,ny,nz, dtsml1, pprt, rhobar,pforce,pdiv)

c
c#######################################################################
c
c     Set boundary conditions for the pressure equation.
c
c#######################################################################
c

      CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb,
     :         ebc,wbc,nbc,sbc,tbc,bbc)

      RETURN
      END
