c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ADVUVW                     ######
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 ADVUVW(nx,ny,nz,u,v,w,rhobar,
     :                  uadv,vadv,wadv, 
     :                  urho,vrho,wrho,tem1,tem2)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the advection terms uadv, vadv and wadv of the u, v and w 
c     equations. These terms are cast in equivalent advection form.  
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/20/91.
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (M. Xue) 
c     Added full documentation.
c
c     5/29/92 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       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 velocity in Cartesian 
c                coordinates at all time levels (m/s).
c
c       rhobar   Base state air density (kg/m**3)
c
c     OUTPUT:
c
c       uadv     Advection terms of u-equation.  (kg/m**3)*(m/s)/s
c       vadv     Advection terms of v-equation.  (kg/m**3)*(m/s)/s
c       wadv     Advection terms of the w-equation.  (kg/m**3)*(m/s)/s
c
c     WORK ARRAYS:
c
c       urho     u * rhobar
c       vrho     v * rhobar
c       wrho     w * rhobar
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c
c      (These arrays are defined and used locally (i.e. inside this
c       subroutine), they may also be passed into routines called by
c       this one. Exiting the call to this subroutine, these temporary
c       work arrays may be used for other purposes, and therefore their
c       contents may be overwritten. Please examine the usage of work
c       arrays before you alter the code.)
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
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)
c
      integer nx,ny,nz          ! Number of grid points in 3 directions

      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 rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real uadv  (nx,ny,nz)     ! advection terms of u-eqn (kg/m**3)*(m/s)/s
      real vadv  (nx,ny,nz)     ! advection terms of v-eqn (kg/m**3)*(m/s)/s
      real wadv  (nx,ny,nz)     ! advection terms of the w-eq. 
                                ! (kg/m**3)*(m/s)/s
c
      real urho  (nx,ny,nz)     ! u * rhobar
      real vrho  (nx,ny,nz)     ! v * rhobar
      real wrho  (nx,ny,nz)     ! w * rhobar
c
      real tem1  (nx,ny,nz)     ! temporary work array
      real tem2  (nx,ny,nz)     ! temporary work array
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer tlevel
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external uvwrho
      external advu
      external advv
      external advw
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
       tlevel=tpresent
c
c#######################################################################
c
c     Calculate urho=rhobar*u, vrho=rhobar*v, wrho=rhobar*w
c
c#######################################################################
c
      CALL uvwrho(nx,ny,nz,
     :            u(1,1,1,tlevel),v(1,1,1,tlevel),w(1,1,1,tlevel),
     :            rhobar, urho,vrho,wrho, tem1)
c
C#######################################################################
c
c     Calculate uadv:
c
C#######################################################################
c
      CALL advu(nx,ny,nz,u(1,1,1,tlevel),urho,vrho,wrho, uadv,
     :          tem1,tem2)
c
C#######################################################################
c
c     Calculate vadv:
c
C#######################################################################
c
      CALL advv(nx,ny,nz,v(1,1,1,tlevel),urho,vrho,wrho, vadv,
     :          tem1,tem2)
c
C#######################################################################
c
c     Calculate wadv:
c
C#######################################################################
c
      CALL advw(nx,ny,nz,w(1,1,1,tlevel),urho,vrho,wrho, wadv,
     :          tem1,tem2)
      
      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ADVP                       ######
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 ADVP(nx,ny,nz,pprt,u,v,w,rhobar,
     :                padv, 
     :                tem1,tem2)
 
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the advection of total pressure.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/20/91.
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (M. Xue) 
c     Added full documentation.
c
c     5/29/92 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       pprt     Perturbation pressure at all time levels (Pascal)
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 velocity in Cartesian 
c                coordinates at al time levels (m/s).
c
c       rhobar   Base state air density (kg/m**3)
c
c     OUTPUT:
c
c       padv     Advection term of pressure equation (kg/(m*s**3))
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c
c      (These arrays are defined and used locally (i.e. inside this
c       subroutine), they may also be passed into routines called by
c       this one. Exiting the call to this subroutine, these temporary
c       work arrays may be used for other purposes, and therefore their
c       contents may be overwritten. Please examine the usage of work
c       arrays before you alter the code.)
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
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)
c
      integer nx,ny,nz          ! Number of grid points in 3 directions

      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 pprt  (nx,ny,nz,nt)  ! Perturbation pressure from that
                                ! of base state atmosphere (Pascal)
c
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real padv  (nx,ny,nz)     ! Advection term of the pressure equation
c                               ! (kg/(m*s**3)
c
      real tem1(nx,ny,nz)       ! Temporary work array
      real tem2(nx,ny,nz)       ! Temporary work array
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i,j,k,onvf,tlevel
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'phycst.inc'        ! Physical constants
 
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external advcts
      external avgz
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      tlevel=tpresent
c
C#######################################################################
c
c     Advection of the perturbation pressure 
c
C#######################################################################
c
      CALL advcts(nx,ny,nz,pprt(1,1,1,tlevel),u(1,1,1,tlevel),
     :            v(1,1,1,tlevel),w(1,1,1,tlevel), padv, tem1,tem2)
c
C#######################################################################
c
c     Calculate the advection of base state pressure, which is the added
c     to array padv to yield the advection of total pressure.
c
c     The formulation of advection of base state pressure is simplified 
c     by making use of the base state hydrostatic relation. 
c
C#######################################################################
c
      onvf = 0
      CALL avgz(w(1,1,1,tlevel), onvf,
     :     nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-2, tem1)
      
      DO 300 k=2,nz-2
      DO 300 j=2,ny-2
      DO 300 i=2,nx-2
        padv(i,j,k)=padv(i,j,k)
     :              -g*tem1(i,j,k)*rhobar(i,j,k)
300   CONTINUE

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                 SUBROUTINE ADVPT                     ######
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 ADVPT(nx,ny,nz,ptprt,u,v,w, rhobar,ptbar,
     :           ptadv,
     :           urho,vrho,wrho,tem1,tem2)

c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the advection of total potential temperature
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/20/91.
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (M. Xue) 
c     Added full documentation.
c
c     5/29/92 (K. Brewster)
c     Further facelift and order of argument list changed.
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       rhobar   Base state air density (kg/m**3)
c       ptbar    Base state potential temperature (K)
c
c     OUTPUT:
c
c       ptadv    Advection term of potential temperature eqn (kg/m**3)*K/s
c
c     WORK ARRAYS:
c
c       urho     u * rhobar
c       vrho     v * rhobar
c       wrho     w * rhobar
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c
c      (These arrays are defined and used locally (i.e. inside this
c       subroutine), they may also be passed into routines called by
c       this one. Exiting the call to this subroutine, these temporary
c       work arrays may be used for other purposes, and therefore their
c       contents may be overwritten. Please examine the usage of work
c       arrays before you alter the code.)
c
c#######################################################################
c
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
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)
c
      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)
c
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real ptbar (nx,ny,nz)     ! Base state potential temperature (K)
c
      real ptadv(nx,ny,nz)      ! Advection of potential temperature
c
      real urho  (nx,ny,nz)     ! u * rhobar
      real vrho  (nx,ny,nz)     ! v * rhobar
      real wrho  (nx,ny,nz)     ! w * rhobar
c
      real tem1  (nx,ny,nz)     ! temporary work array
      real tem2  (nx,ny,nz)     ! temporary work array

c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i,j,k,onvf,tlevel
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'globcst.inc'       ! Global model control parameters 
c
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external uvwrho
      external advcts
      external difz
      external aamult
      external avgz
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      tlevel=tpresent
c
c#######################################################################
c
c     Calculate urho=rhobar*u, vrho=rhobar*v, wrho=rhobar*w
c
c#######################################################################
c
      CALL uvwrho(nx,ny,nz,
     :            u(1,1,1,tlevel),v(1,1,1,tlevel),w(1,1,1,tlevel),
     :            rhobar, urho,vrho,wrho, tem1)
c
C#######################################################################
c
c     Advection of the perturbation potential temperature
c
C#######################################################################
c
      CALL advcts(nx,ny,nz,ptprt(1,1,1,tlevel),urho,vrho,wrho, 
     :            ptadv, tem1,tem2)

c
C#######################################################################
c
c     Advection of base state potential temperature which is added to
c     array ptadv to yield the advection of total potential temperature.
c
C#######################################################################
c
      onvf = 1
      CALL difz(ptbar, onvf,
     :     nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-1, dz, tem1)
c
c     At this point tem1 has d(ptbar)/dz
c
      CALL aamult(wrho, tem1,
     :     nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-1, tem2)
c
c     At this point tem2 has w*rho*d(ptbar)/dx
c
      onvf = 0
      CALL avgz(tem2 , onvf,
     :     nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-2, tem1)
c
c     At this point tem1 has w*rho*d(ptbar)/dx averaged in z
c
      DO 300 k=2,nz-2
      DO 300 j=2,ny-2
      DO 300 i=2,nx-2
        ptadv(i,j,k)=ptadv(i,j,k)+ tem1(i,j,k)
300   CONTINUE

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                 SUBROUTINE ADVQ                      ######
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 ADVQ(nx,ny,nz,q,u,v,w,rhobar,
     :                qadv,
     :                urho,vrho,wrho,tem1,tem2)

c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the advection of a water/ice quantity. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     05/20/92.
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (M. Xue) 
c     Added full documentation.
c
c     5/29/92 (K. Brewster)
c     Further facelift and order of argument list changed.
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 of a water or ice quantity (kg/kg)
c
c       rhobar   Base state air density (kg/m**3)
c
c     OUTPUT:
c
c       qadv     Water or ice mixing ratio advection term in 
c                conservation equation (kg/m**3)*(kg/kg)/s
c
c     WORK ARRAYS:
c
c       urho     u * rhobar
c       vrho     v * rhobar
c       wrho     w * rhobar
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c
c      (These arrays are defined and used locally (i.e. inside this
c       subroutine), they may also be passed into routines called by
c       this one. Exiting the call to this subroutine, these temporary
c       work arrays may be used for other purposes, and therefore their
c       contents may be overwritten. Please examine the usage of work
c       arrays before you alter the code.)
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
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)
c
      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 rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real qadv(nx,ny,nz)       ! Advction of a water/ice variable 
c
      real urho  (nx,ny,nz)     ! u * rhobar
      real vrho  (nx,ny,nz)     ! v * rhobar
      real wrho  (nx,ny,nz)     ! w * rhobar
c
      real tem1  (nx,ny,nz)     ! temporary work array
      real tem2  (nx,ny,nz)     ! temporary work array
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer tlevel
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'globcst.inc'       ! Global model control parameters 
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external uvwrho
      external advcts
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      tlevel=tpresent
c
c#######################################################################
c
c     Calculate urho=rhobar*u, vrho=rhobar*v, wrho=rhobar*w
c
c#######################################################################
c
      CALL uvwrho(nx,ny,nz,
     :            u(1,1,1,tlevel),v(1,1,1,tlevel),w(1,1,1,tlevel),
     :            rhobar, urho,vrho,wrho, tem1)
c
C#######################################################################
c
c     Advection of the water/ice mixing ratio 
c
C#######################################################################
c
      CALL advcts(nx,ny,nz,q(1,1,1,tlevel),urho,vrho,wrho,
     :            qadv, tem1,tem2)

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                 SUBROUTINE ADVU                      ######
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 ADVU(nx,ny,nz,u,urho,vrho,wrho, 
     :                uadv, 
     :                tem1,tem2)

c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the advection term of the u-equation. This term is cast 
c     in equivalent advection form.  
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/20/91.
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (M. Xue) 
c     Added full documentation.
c
c     5/29/92 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity at a given time level (m/s)
c
c       urho     u * rhobar
c       vrho     v * rhobar
c       wrho     w * rhobar
c
c     OUTPUT:
c
c       uadv     Advection term of the u equation (kg/m**3)*(m/s)/s
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c
c      (These arrays are defined and used locally (i.e. inside this
c       subroutine), they may also be passed into routines called by
c       this one. Exiting the call to this subroutine, these temporary
c       work arrays may be used for other purposes, and therefore their
c       contents may be overwritten. Please examine the usage of work
c       arrays before you alter the code.)
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      real u     (nx,ny,nz)     ! Total u-velocity at a given time level(m/s)
      real urho  (nx,ny,nz)     ! u*rhobar        
      real vrho  (nx,ny,nz)     ! v*rhobar
      real wrho  (nx,ny,nz)     ! w*rhobar
c
      real uadv  (nx,ny,nz)     ! Advection term of the u-eq. 
                                ! (kg/m**3)*(m/s)/s
c
      real tem1  (nx,ny,nz)     ! temporary work array
      real tem2  (nx,ny,nz)     ! temporary work array
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i,j,k,onvf
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'globcst.inc'       ! Global model control parameters 
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external avgx
      external avgy
      external avgz
      external difx
      external dify
      external difz
      external aamult
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
C#######################################################################
c
c     Calculate rho*u* du/dx = avgx( avgx( urho )*difx( u ) )
c
C#######################################################################
c
      onvf = 0
      CALL avgx(urho, onvf,
     :      nx,ny,nz, 1,nx-1, 2,ny-2, 2,nz-2, tem1)
c
      onvf = 0
      CALL difx(u, onvf,
     :      nx,ny,nz, 1,nx-1, 2,ny-2, 2,nz-2, dx, tem2)
c
c      At this point tem1 contains avgx(u*rho)
c      At this point tem2 contains du/dx
c
      CALL aamult(tem1, tem2,
     :      nx,ny,nz, 1,nx-1, 2,ny-2, 2,nz-2, tem2)
c
c      At this point tem2 contains (avgx(u*rho))*du/dx
c
      onvf = 1
      CALL avgx(tem2, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-2, 2,nz-2, uadv)
c
C#######################################################################
c
c     Calculate rho*v* du/dy = avgy( avgx( vrho )* dify( u ) )
c
C#######################################################################
c
      onvf = 1
      CALL avgx(vrho, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-1, 2,nz-2, tem1)
c
      onvf = 1
      CALL dify(u, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-1, 2,nz-2, dy, tem2)
c
c      At this point tem1 contains avgx(v*rho)
c      At this point tem2 contains du/dy
c
      CALL aamult(tem1, tem2,
     :      nx,ny,nz, 2,nx-1, 2,ny-1, 2,nz-2, tem2)
c
c     At this point tem2 contains (avgx(v*rho)) * du/dy
c
      onvf = 0
      CALL avgy(tem2, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-2, 2,nz-2, tem1)
c
C#######################################################################
c
c     Sum the x and y contributions to u advection
c     At this point tem1 contains the y contribution to u advection
c
C#######################################################################
c
      DO 300 k=2,nz-2
      DO 300 j=2,ny-2
      DO 300 i=2,nx-1
        uadv(i,j,k)=uadv(i,j,k)+tem1(i,j,k)
300   CONTINUE
c
C#######################################################################
c
c     Calculate rho*w* du/dz = avgz( avgx( wrho ) * difz ( u ) )
c
C#######################################################################
c
      onvf = 1
      CALL avgx(wrho, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-2, 2,nz-1, tem1)
c
      onvf = 1
      CALL difz(u, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-2, 2,nz-1, dz, tem2)
c
c     At this point tem1 contains avgx(w*rho)
c     At this point tem2 contains du/dz
c
      CALL aamult(tem1, tem2,
     :      nx,ny,nz, 2,nx-1, 2,ny-2, 2,nz-1, tem2)
c
c     At this point tem2 contains avgx(w*rho) * du/dz
c
      onvf = 0
      CALL avgz(tem2, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-2, 2,nz-2, tem1)
c
C#######################################################################
c
c     Sum the x,y and z contributions to u advection
c     At this point tem1 contains the z contribution to u advection
c
C#######################################################################
c
      DO 400 k=2,nz-2
      DO 400 j=2,ny-2
      DO 400 i=2,nx-1
        uadv(i,j,k)=uadv(i,j,k)+tem1(i,j,k)
400   CONTINUE
c
      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ADVV                       ######
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 ADVV(nx,ny,nz,v,urho,vrho,wrho, 
     :                vadv, 
     :                tem1,tem2)

c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the advection term of the v-equation. This term is cast 
c     in equivalent advection form.  
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/20/91.
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (M. Xue) 
c     Added full documentation.
c
c     5/29/92 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       v        y component of velocity at a given time level (m/s)
c
c       urho     u * rhobar
c       vrho     v * rhobar
c       wrho     w * rhobar
c
c     OUTPUT:
c
c       vadv     Advection term in v equation (kg/m**3)*(m/s)/s
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c
c      (These arrays are defined and used locally (i.e. inside this
c       subroutine), they may also be passed into routines called by
c       this one. Exiting the call to this subroutine, these temporary
c       work arrays may be used for other purposes, and therefore their
c       contents may be overwritten. Please examine the usage of work
c       arrays before you alter the code.)
c
c#######################################################################
c
      real v     (nx,ny,nz)     ! Total v-velocity at a given time level(m/s)
      real urho  (nx,ny,nz)     ! u*rhobar        
      real vrho  (nx,ny,nz)     ! v*rhobar
      real wrho  (nx,ny,nz)     ! w*rhobar
c
      real vadv  (nx,ny,nz)     ! Advection term of the v-eq. 
                                ! (kg/m**3)*(m/s)/s
c
      real tem1  (nx,ny,nz)     ! temporary work array
      real tem2  (nx,ny,nz)     ! temporary work array
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i,j,k,onvf
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'globcst.inc'       ! Global model control parameters 
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external avgx
      external avgy
      external avgz
      external difx
      external dify
      external difz
      external aamult
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
C#######################################################################
c
c     Calculate rho*u* dv/dx  = avgx( avgy( urho ) * difx( v ) ) 
c
C#######################################################################
c
      onvf = 1
      CALL avgy(urho, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-1, 2,nz-2, tem1)
c
      onvf = 1
      CALL difx(v, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-1, 2,nz-2, dx, tem2)
c
c     At this point tem1 contains avgy(u*rho)
c     At this point tem2 contains dv/dx
c
      CALL aamult(tem1, tem2,
     :      nx,ny,nz, 2,nx-1, 2,ny-1, 2,nz-2, tem2)
c
C#######################################################################
c
c     Store the x contribution to the v advection in the vadv array.
c
c     At this point tem2 contains avgy(u*rho) * dv/dx
c
C#######################################################################
c
      onvf = 0
      CALL avgx(tem2, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-1, 2,nz-2, vadv)

c
C#######################################################################
c
c     Calculate rho*v* dv/dy = avgy( avgy( vrho ) * dify( v ) )  
c
C#######################################################################
c
      onvf = 0
      CALL avgy(vrho, onvf,
     :      nx,ny,nz, 2,nx-2, 1,ny-1, 2,nz-2, tem1)
c
      onvf = 0
      CALL dify(v, onvf,
     :      nx,ny,nz, 2,nx-2, 1,ny-1, 2,nz-2, dy, tem2)
c
c     At this point tem1 contains avgy(v*rho)
c     At this point tem2 contains dv/dy
c
      CALL aamult(tem1, tem2,
     :      nx,ny,nz, 2,nx-2, 1,ny-1, 2,nz-2, tem2)
c
c     At this point tem2 contains avgy(v*rho) * dv/dy
c
      onvf = 1
      CALL avgy(tem2, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-1, 2,nz-2, tem1)
c
C#######################################################################
c
c     Sum the x and y contributions to v advection
c     At this point tem1 contains the y contribution to the v advection
c
C#######################################################################
c
      DO 300 k=2,nz-2
      DO 300 j=2,ny-1
      DO 300 i=2,nx-2
        vadv(i,j,k)=vadv(i,j,k)+tem1(i,j,k)
300   CONTINUE
c
C#######################################################################
c
c     Calculate rho*w* dv/dz = avgz( avgy( wrho ) * difz( v ) ) 
c
C#######################################################################
c

      onvf = 1
      CALL avgy(wrho, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-1, 2,nz-1, tem1)
c
      onvf = 1
      CALL difz(v, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-1, 2,nz-1, dz, tem2)
c
c     At this point tem1 contains avgy(w*rho)
c     At this point tem2 contains dv/dz
c
      CALL aamult(tem1, tem2,
     :      nx,ny,nz, 2,nx-2, 2,ny-1, 2,nz-1, tem2)
c
c     At this point tem2 contains avgy(w*rho)*dv/dz
c
      onvf = 0
      CALL avgz(tem2, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-1, 2,nz-2, tem1)
c
c#######################################################################
c
c     Sum the x,y and z contributions to v advection
c     At this point tem1 contains the z contribution to v advection
c
c#######################################################################
c
      DO 400 k=2,nz-2
      DO 400 j=2,ny-1
      DO 400 i=2,nx-2
        vadv(i,j,k)=vadv(i,j,k)+tem1(i,j,k)
400   CONTINUE

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ADVW                       ######
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 ADVW(nx,ny,nz,w,urho,vrho,wrho, 
     :                wadv, 
     :                tem1,tem2)

c
c#######################################################################
c
c     PURPOSE:
c
c
c     Calculate the advection term of the w-equation. This term is cast 
c     in equivalent advection form.  
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/20/91.
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (M. Xue) 
c     Added full documentation.
c
c     5/29/92 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       w        z component of velocity at a given time level (m/s)
c
c       urho     u * rhobar
c       vrho     v * rhobar
c       wrho     w * rhobar
c
c     OUTPUT:
c
c       wadv     Advection term in w equation (kg/m**3)*(m/s)/s
c
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c
c      (These arrays are defined and used locally (i.e. inside this
c       subroutine), they may also be passed into routines called by
c       this one. Exiting the call to this subroutine, these temporary
c       work arrays may be used for other purposes, and therefore their
c       contents may be overwritten. Please examine the usage of work
c       arrays before you alter the code.)
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
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)
c
      integer nx,ny,nz          ! Number of grid points in 3 directions

      real w     (nx,ny,nz,nt)  ! Total w-velocity (m/s)
      real urho  (nx,ny,nz)     ! u*rhobar        
      real vrho  (nx,ny,nz)     ! v*rhobar
      real wrho  (nx,ny,nz)     ! w*rhobar
c
      real wadv  (nx,ny,nz)     ! Advection term of w-eqn (kg/m**3)*(m/s)/s
c
      real tem1  (nx,ny,nz)     ! temporary work array
      real tem2  (nx,ny,nz)     ! temporary work array
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i,j,k,onvf
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'globcst.inc'       ! Global model control parameters 
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external avgx
      external avgy
      external avgz
      external difx
      external dify
      external difz
      external aamult
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
C#######################################################################
c
c     Calculate rho*u* dw/dx = avgx( avgz( urho ) * difx( w ) ) 
c
C#######################################################################
c
      onvf = 1
      CALL avgz(urho, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-2, 2,nz-1, tem1)
c
      onvf = 1
      CALL difx(w, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-2, 2,nz-1, dx, tem2)
c
c     At this point tem1 contains avgz(u*rho)
c     At this point tem2 contains dw/dx
c
      CALL aamult(tem1, tem2,
     :      nx,ny,nz, 2,nx-1, 2,ny-2, 2,nz-1, tem2)
c
c     At this point tem2 contains avgz(u*rho)*dw/dx
c
      onvf = 0
      CALL avgx(tem2, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-1, wadv)
c
C#######################################################################
c
c     Calculate rho*v* dw/dy = avgy( avgz( vrho ) * dify ( w ) ) 
c
C#######################################################################
c
      onvf = 1
      CALL avgz(vrho, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-1, 2,nz-1, tem1)
c
      onvf = 1
      CALL dify(w, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-1, 2,nz-1, dy, tem2)
c
c     At this point tem1 contains avgz(v*rho)
c     At this point tem2 contains dw/dy
c
      CALL aamult(tem1, tem2,
     :      nx,ny,nz, 2,nx-2, 2,ny-1, 2,nz-1, tem2)
c
c     At this point tem2 contains avgz(v*rho)*dw/dy
c
      onvf = 0
      CALL avgy(tem2, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-1, tem1)
c
c#######################################################################
c
c     Sum the x and y contributions to w advection
c     At this point tem1 contains the y contribution to w advection
c
c#######################################################################
c
      DO 300 k=2,nz-1
      DO 300 j=2,ny-2
      DO 300 i=2,nx-2
        wadv(i,j,k)=wadv(i,j,k)+tem1(i,j,k)
300   CONTINUE
c
C#######################################################################
c
c     Calculate rho*w* dw/dz = avgz( avgz( wrho ) * difz( w ) ) 
c
C#######################################################################
c
      onvf = 0
      CALL avgz(wrho, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 1,nz-1, tem1)
c
      onvf = 0
      CALL difz(w, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 1,nz-1, dz, tem2)
c
c     At this point tem1 contains avgz(w*rho)
c     At this point tem2 contains dw/dz
c
      CALL aamult(tem1, tem2,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 1,nz-1, tem2)
c
c     At this point tem2 contains avgz(w*rho)*dw/dz
c
      onvf = 1
      CALL avgz(tem2, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-1, tem1)
c
c#######################################################################
c
c     Sum the x,y and z contributions to w advection
c     At this point tem1 contains the z contribution to w advection
c
c#######################################################################
c
      DO 400 k=2,nz-1
      DO 400 j=2,ny-2
      DO 400 i=2,nx-2
        wadv(i,j,k)=wadv(i,j,k)+tem1(i,j,k)
400   CONTINUE

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ADVCTS                     ######
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 ADVCTS(nx,ny,nz,s, urho,vrho,wrho,
     :                  sadv,
     :                  tem1,tem2)
 
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the advection term of the scalar equations. This term is 
c     cast in equivalent advection form.  
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/20/91.
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (M. Xue) 
c     Added full documentation.
c
c     5/29/92 (K. Brewster)
c     Further facelift.
c
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       s        scalar field at a given time level (scalar units)
c
c       urho     u * rhobar
c       vrho     v * rhobar
c       wrho     w * rhobar
c
c     OUTPUT:
c
c       sadv     Advection term in scalar equation (kg/m**3)*scalar unit/s
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c
c      (These arrays are defined and used locally (i.e. inside this
c       subroutine), they may also be passed into routines called by
c       this one. Exiting the call to this subroutine, these temporary
c       work arrays may be used for other purposes, and therefore their
c       contents may be overwritten. Please examine the usage of work
c       arrays before you alter the code.)
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      real s     (nx,ny,nz)     ! a scalar field to be advected.
      real urho  (nx,ny,nz)     ! u * rhobar
      real vrho  (nx,ny,nz)     ! v * rhobar
      real wrho  (nx,ny,nz)     ! w * rhobar
c
      real sadv  (nx,ny,nz)     ! advection term in scalar equation
c                               ! (kg/m**3)*(scalar units)/s
c
      real tem1  (nx,ny,nz)     ! temporary work array
      real tem2  (nx,ny,nz)     ! temporary work array
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i,j,k,onvf
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'globcst.inc'        ! Physical constants
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external avgx
      external avgy
      external avgz
      external difx
      external dify
      external difz
      external aamult
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
C#######################################################################
c
c     Calculate rho*u* ds/dx = avgx( urho * difx ( s ) )  
c
C#######################################################################
c
      onvf = 1
      CALL difx(s, onvf,
     :      nx,ny,nz, 2,nx-1, 2,ny-2, 2,nz-2, dx ,tem1)
c
c     At this point tem1 contains ds/dx
c
      CALL aamult(urho, tem1,
     :      nx,ny,nz, 2,nx-1, 2,ny-2, 2,nz-2, tem1)
c
c     At this point tem1 contains (u*rho)*ds/dx
c
      onvf = 0
      CALL avgx(tem1, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-2, sadv)
c
C#######################################################################
c
c     Calculate rho*v* ds/dy = avgy( vrho * dify( s ) )  
c
C#######################################################################
c
      onvf = 1
      CALL dify(s, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-1, 2,nz-2, dy, tem1)
c
c     At this point tem1 contains ds/dy
c
      CALL aamult(vrho, tem1,
     :      nx,ny,nz, 2,nx-2, 2,ny-1, 2,nz-2, tem1)
c
c     At this point tem1 contains (v*rho)*ds/dy
c
      onvf = 0
      CALL avgy(tem1, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-2, tem2)
c
c#######################################################################
c
c     Sum the x and y contributions to s advection
c     At this point tem2 contains y contribution to s advection
c
c#######################################################################
c
      DO 300 k=2,nz-2
      DO 300 j=2,ny-2
      DO 300 i=2,nx-2
        sadv(i,j,k)=sadv(i,j,k)+tem2(i,j,k)
300   CONTINUE
c
C#######################################################################
c
c     Calculate rho*w* ds/dz = avgz( wrho * difz( s ) )  
c
C#######################################################################
c
      onvf = 1
      CALL difz(s, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-1, dz, tem1)
c
c     At this point tem1 contains ds/dz
c
      CALL aamult(wrho, tem1,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-1, tem1)
c
c     At this point tem1 contains (w*rho)*ds/dz
c
      onvf = 0
      CALL avgz(tem1, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-1, tem2)
c
C#######################################################################
c
c     Sum the x, y and z contributions to s advection
c     At this point tem2 contains z contribution to s advection
c
C#######################################################################
c
      DO 400 k=2,nz-2
      DO 400 j=2,ny-2
      DO 400 i=2,nx-2
        sadv(i,j,k)=sadv(i,j,k)+tem2(i,j,k)
400   CONTINUE

      RETURN
      END
