C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE MIXUVW                     ######
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 MIXUVW(nx,ny,nz,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :           usflx,vsflx, x,y,z,zp, j1,j2,j3,
     :           umix,vmix,wmix, km,
     :           tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9)
c
c#######################################################################
c
c     PURPOSE:
c
C     Calculate the total mixing terms (turbulent mixing and the
c     externally imposed computational mixing) for the momentum equations. 
c     This subroutine also calculates the turbulent mixing coefficient km 
c     which will be used to calculate mixing terms for temperature and 
c     water quantities. The mixing coefficient km is based on Smagorinsky's
c     formulation. For the computational mixing, there are two options: 
c     second order and fourth order mixing. The mixing applies only to 
c     perturbations from the base state. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue)
c     Added full documentation.
c
c     6/1/92 (M. Xue and H. Jin) 
c     Further facelift. 
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity at a given time level (m/s)
c       v        y component of velocity at a given time level (m/s)
c       w        Vertical component of velocity in Cartesian
c                coordinates at a given time level (m/s).
c       ptprt    Perturbation potential temperature at a given time level (K)
c       pprt     Perturbation pressure at a given time level (Pascal)
c       qv       Water vapor specific humidity at a given time level (kg/kg)
c       qc       Cloud water mixing ratio at a given time level (kg/kg)
c       qr       Rainwater mixing ratio at a given time level (kg/kg)
c       qi       Cloud ice mixing ratio at a given time level (kg/kg)
c       qs       Snow mixing ratio at a given time level (kg/kg)
c       qh       Hail mixing ratio at a given time level (kg/kg)
c
c       ubar     Base state x velocity component (m/s)
c       vbar     Base state y velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state air density (kg/m**3)
c       qvbar    Base state water vapor specific humidity (kg/kg)
c
c       usflx    Surface flux of u-momentum
c       vsflx    Surface flux of v-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     OUTPUT:
c
c       umix     Total mixing in u equation (kg/(m*s)**2)
c       vmix     Total mixing in v equation (kg/(m*s)**2)
c       wmix     Total mixing in w equation (kg/(m*s)**2)
c       km       Turbulent mixing coefficient for momentum ( m**2/s ).
c
c     WORK ARRAYS:
c
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       tem7     Temporary work array.
c       tem8     Temporary work array.
c       tem9     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
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real u     (nx,ny,nz)     ! Total u-velocity (m/s)
      real v     (nx,ny,nz)     ! Total v-velocity (m/s)
      real w     (nx,ny,nz)     ! Total w-velocity (m/s)
      real ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
      real pprt  (nx,ny,nz)     ! Perturbation pressure (Pascal)
c
      real qv    (nx,ny,nz)     ! Water vapor specific humidity (kg/kg)
      real qc    (nx,ny,nz)     ! Cloud water mixing ratio (kg/kg)
      real qr    (nx,ny,nz)     ! Rain water mixing ratio (kg/kg)
      real qi    (nx,ny,nz)     ! Cloud ice mixing ratio (kg/kg)
      real qs    (nx,ny,nz)     ! Snow mixing ratio (kg/kg)
      real qh    (nx,ny,nz)     ! Hail mixing ratio (kg/kg)
c
      real ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
      real vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
      real ptbar (nx,ny,nz)     ! Base state potential temperature (K)
      real pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
      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)

      real usflx (nx,ny)        ! surface flux of u-momentum (kg/(m*s**2))
      real vsflx (nx,ny)        ! surface flux of v-momentum (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 umix  (nx,ny,nz)     ! Turbulent mixing on u-momentum.
      real vmix  (nx,ny,nz)     ! Turbulent mixing on v-momentum.
      real wmix  (nx,ny,nz)     ! Turbulent mixing on w-momentum.
      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )
c
      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
      real tem7  (nx,ny,nz)     ! Temporary work array
      real tem8  (nx,ny,nz)     ! Temporary work array
      real tem9  (nx,ny,nz)     ! Temporary work array
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer nxyz
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external flzero
      external tmixuvw
      external cmix2uvw
      external cmix4uvw
      external rdmpuvw
c
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Control parameters for subgrid scalar turbulent mixing are passed
c     in globcst.inc, these include:
c
c     tmixopt: =0,1,2 or 3 for zero turb. mixing, constant coeff. mixing,
c              Smagorinsky mixing and Smargorinsky+constant coeff. mixing.
c     tmixcst: the constant mixing coeff. for option 1 and 3.
c
c#######################################################################
c
      IF( tmixopt.eq.0) THEN 

        nxyz = nx*ny*nz
        CALL flzero(umix,nxyz)
        CALL flzero(vmix,nxyz)
        CALL flzero(wmix,nxyz)

      ELSE

        CALL tmixuvw(nx,ny,nz,
     :       u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :       ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :       usflx,vsflx, x,y,z,zp, j1,j2,j3,
     :       umix,vmix,wmix,km, 
     :       tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9)

      ENDIF

      IF( kcmix2 .eq. 1) THEN 
c
c#######################################################################
c
c     calculate second order computational mixing terms and add
c     them to the turbulent mixing terms umix, vmix, wmix.
c
c#######################################################################
c
        CALL cmix2uvw(nx,ny,nz,
     :                u,v,w, ubar,vbar,rhobar,
     :                umix,vmix,wmix, 
     :                tem1,tem2,tem3)

      ENDIF

      IF( kcmix4 .eq. 1) THEN 
c
c#######################################################################
c
c     Calculate fourth order computational mixing terms and add
c     them to the turbulent mixing terms umix, vmix, wmix.
c
c#######################################################################
c
        CALL cmix4uvw(nx,ny,nz,
     :                u,v,w, ubar,vbar,rhobar,
     :                umix,vmix,wmix, 
     :                tem1,tem2,tem3,tem4)

      ENDIF

      IF( raydmp.ne.0) THEN 
c
c#######################################################################
c
c    calculate upper level Rayleigh damping on u, v and w perturbations. 
c    The terms are accumulated in array umix, vmix and wmix.
c
c#######################################################################
c
        CALL rdmpuvw(nx,ny,nz,
     :               u,v,w, ubar,vbar,rhobar, zp,
     :               umix,vmix,wmix, tem1)

      ENDIF

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE MIXPT                      ######
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 MIXPT(nx,ny,nz,
     :           ptprt,ptbar,rhobar,km,ptsflx,x,y,z,zp,j1,j2,j3,
     :           ptmix,
     :           tem1,tem2,tem3,tem4,tem5,tem6)
c
c#######################################################################
c
c     PURPOSE:
c
C     Calculate the mixing term (turbulent mixing and an externally imposed 
c     computational mixing) for potential temperature equation. Turbulent 
c     mixing coefficient for potential temperature Kh is Km/Prandtl where 
c     KM is the mixing coefficient for momentum obtained in subroutine 
c     mixuvw, based on modified Smagorinsky formulation. For the computational 
c     mixing, there are two options: second order and fourth order mixing. 
c     The mixing applies only to perturbations from the base state. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
C
c     6/1/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       ptprt    Perturbation potential temperature at a given time level (K)
c       ptbar    Base state potential temperature (K)
c       rhobar   Base state air density (kg/m**3)
c
c       km       Turbulent mixing coefficient for momentum ( m**2/s ).
c
c       ptsflx   Surface flux of heat (K*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     OUTPUT:
c
c       ptmix    Total mixing in potential temperature equation 
c                (K*kg/(m**3 * s))
c
c     WORK ARRAYS:
c
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      (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
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
      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 ptmix(nx,ny,nz)      ! Mixing on potential 
                                ! temperature (K*kg/(m**3 * s))
c
      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
      integer nxyz
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c

      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external flzero
      external tmixpt
      external cmix2pt
      external cmix4pt
      external rdmppt
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      IF( tmixopt.eq.0) THEN 

        nxyz = nx*ny*nz
        CALL flzero(ptmix,nxyz)

      ELSE

        CALL tmixpt(nx,ny,nz, ptprt,ptbar,rhobar,km,ptsflx, 
     :             x,y,z,zp, j1,j2,j3,
     :             ptmix, 
     :             tem1,tem2,tem3,tem4,tem5)

      ENDIF

      IF( kcmix2 .eq. 1) THEN 
c
c#######################################################################
c
c     Compute the second order computational mixing term ptmix for
c     the potential temperature equation. 
c
c#######################################################################
c

        CALL cmix2pt(nx,ny,nz, ptprt,rhobar, ptmix, 
     :               tem1,tem2,tem3)


      ENDIF

      IF( kcmix4 .eq. 1) THEN 
c
c#######################################################################
c
c     calculate fourth order computational mixing term ptmix for 
c     the potential temperature equation.
c
c#######################################################################
c

        CALL cmix4pt(nx,ny,nz, ptprt,rhobar, ptmix, 
     :               tem1,tem2,tem3,tem4)

      ENDIF

      IF( raydmp.ne.0) THEN 
c
c#######################################################################
c
c    calculate upper level Rayleigh damping on u, v and w perturbations. 
c    The term is accumulated in array ptmix.
c
c#######################################################################
c
        CALL rdmppt(nx,ny,nz,
     :              ptprt ,rhobar, zp,
     :              ptmix, tem1)

      ENDIF

      RETURN
      END
C
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE MIXQV                      ######
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 MIXQV(nx,ny,nz,
     :           qv,rhobar,qvbar,km,qvsflx, x,y,z,zp,j1,j2,j3,
     :           qvmix,
     :           tem1,tem2,tem3,tem4,tem5,tem6)
c
c#######################################################################
c
c     PURPOSE:
c
C     Calculate the total mixing term (the turbulent mixing and an 
c     externally imposed computational mixing) for water vapor mixing 
c     ratio. The mixing term is different from those of other water
c     quantities because QV has a base state. Computational mixing 
c     applies to perturbation quantities only.
c
c     Turbulent mixing coefficient for water/ice variables Kh is
c     Km/Prandtl where KM is the mixing coefficient for momentum 
c     obtained in subroutine mixuvw, based on modified Smagorinsky 
c     formulation.
c
c     For the computational mixing, there are two options: second order
c     and fourth order mixing. This type of mixing only applies to 
c     perturbations from the base state. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
C
c     6/1/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       qv       Water vapor specific humidity at a given time level (kg/kg)
c       qvbar    Base state water vapor specific humidity (kg/kg)
c       rhobar   Base state air density (kg/m**3)
c
c       km       Turbulent mixing coefficient for momentum ( m**2/s ).
c
c       qvsflx   Surface flux of moisture (K*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     OUTPUT:
c
c       qvmix    Total mixing in water vapor equation (kg/(m**3 * s))
c
c     WORK ARRAYS:
c
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      (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     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real qv    (nx,ny,nz)     ! Water vapor specific humidity (kg/kg)
      real qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                                ! (kg/kg)
      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 qvsflx(nx,ny)        ! surface flux of moisture (kg/(m**2*s))

      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 qvmix(nx,ny,nz)      ! Mixing on water vapor specific humidity
                                ! (kg/(m**3 *s))

      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
      integer i, j, k
      integer nxyz
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external flzero
      external tmixqv
      external cmix2q 
      external cmix4q 
c
c#######################################################################
c

c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      IF( tmixopt.eq.0) THEN 

        nxyz = nx*ny*nz
        CALL flzero(qvmix,nxyz)

      ELSE 

        CALL tmixqv(nx,ny,nz, qv,rhobar,km,qvsflx, 
     :             x,y,z,zp, j1,j2,j3,
     :             qvmix,
     :             tem1,tem2,tem3,tem4,tem5)

      ENDIF

      IF( kcmix2 .eq. 1) THEN 
c
c#######################################################################
c
c     calculate second order computational mixing qvmix for qv eq. 
c
c#######################################################################
c
        DO 500 k=1,nz
        DO 500 j=1,ny
        DO 500 i=1,nx
          tem5(i,j,k)=qv(i,j,k)-qvbar(i,j,k)
500     CONTINUE

        CALL cmix2q(nx,ny,nz, tem5 ,rhobar, qvmix, 
     :              tem1,tem2,tem3)

      ENDIF

      IF( kcmix4 .eq. 1) THEN 
c
c#######################################################################
c
c     calculate fourth order computational mixing qvmix for qv eq. 
c
c#######################################################################
c
        DO 510 k=1,nz
        DO 510 j=1,ny
        DO 510 i=1,nx
          tem5(i,j,k)=qv(i,j,k)-qvbar(i,j,k)
510     CONTINUE

        CALL cmix4q(nx,ny,nz, tem5 ,rhobar, qvmix, 
     :              tem1,tem2,tem3,tem4)

      ENDIF

      RETURN
      END
C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE MIXQ                       ######
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 MIXQ(nx,ny,nz,
     :           q,rhobar,km,  x,y,z,zp,j1,j2,j3,
     :           qmix,
     :           tem1,tem2,tem3,tem4,tem5,tem6)
c
c#######################################################################
c
c     PURPOSE:
c
C     Calculate the total mixing term (which includes the turbulent mixing
c     as well as externally imposed computational mixing) for the water
c     quantities except water vapor.
c
c     Turbulent mixing coefficient for water/ice variables Kh is 
c     Km/Prandtl where KM is the mixing coefficient for momentum 
c     obtained in subroutine mixuvw, based on modified Smagorinsky 
c     formulation.
c
c     For computational mixing, there are two options: second order
c     and fourth order mixing. This type of mixing only applies to 
c     perturbations from the base state. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
C
c     6/1/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       q        mixing ratio of one of the water/ice quantities (kg/kg)
c       km       Turbulent mixing coefficient for momentum ( m**2/s ).
c       rhobar   Base state air density (kg/m**3)
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     OUTPUT:
c
c       qmix     Total mixing in water/ice substance 
c                equation (kg/(m**3 * s))
c
c     WORK ARRAYS:
c
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      (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     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real q     (nx,ny,nz)     ! One of the water/ice quantities (kg/kg)
      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)

      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 qmix(nx,ny,nz)       ! Mixing on water/ice substances
c
      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
      integer nxyz
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c

      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external flzero
      external tmixq
      external cmix2q
      external cmix4q
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      IF( tmixopt.eq.0) THEN 

        nxyz = nx*ny*nz
        CALL flzero(qmix,nxyz)

      ELSE 

        CALL tmixq(nx,ny,nz, q,rhobar,km, 
     :             x,y,z,zp, j1,j2,j3,
     :             qmix,
     :             tem1,tem2,tem3,tem4,tem5)

      ENDIF

      IF( kcmix2 .eq. 1) THEN 
c
c#######################################################################
c
c     Calculate second order computational mixing qmix for q eq. 
c
c#######################################################################
c
        CALL cmix2q(nx,ny,nz, q ,rhobar, qmix, 
     :              tem1,tem2,tem3)

      ENDIF

      IF( kcmix4 .eq. 1) THEN 
c
c#######################################################################
c
c     Calculate fourth order computational mixing qmix for q eq. 
c
c#######################################################################
c
        CALL cmix4q(nx,ny,nz, q ,rhobar, qmix, 
     :              tem1,tem2,tem3,tem4)

      ENDIF

      RETURN
      END
C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE TMIXUVW                    ######
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 TMIXUVW(nx,ny,nz,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :           usflx,vsflx, x,y,z,zp, j1,j2,j3,
     :           umix,vmix,wmix, km,
     :           nsqed,tau11,tau12,tau13,tau22,tau23,tau33,tem1,tem2)
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the turbulent mixing terms for the momentum equations. These 
c     terms are expressed in terms of the stress tensor.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
C
c     6/1/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity at a given time level (m/s)
c       v        y component of velocity at a given time level (m/s)
c       w        Vertical component of velocity in Cartesian
c                coordinates at a given time level (m/s).
c       ptprt    Perturbation potential temperature at a given time level (K)
c       pprt     Perturbation pressure at a given time level (Pascal)
c       qv       Water vapor specific humidity at a given time level (kg/kg)
c       qc       Cloud water mixing ratio at a given time level (kg/kg)
c       qr       Rainwater mixing ratio at a given time level (kg/kg)
c       qi       Cloud ice mixing ratio at a given time level (kg/kg)
c       qs       Snow mixing ratio at a given time level (kg/kg)
c       qh       Hail mixing ratio at a given time level (kg/kg)
c
c       ubar     Base state x velocity component (m/s)
c       vbar     Base state y velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state air density (kg/m**3)
c       qvbar    Base state water vapor specific humidity (kg/kg)
c
c       usflx    Surface flux of u-momentum
c       vsflx    Surface flux of v-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     OUTPUT:
c
c       umix     Total mixing in u equation (kg/(m*s)**2)
c       vmix     Total mixing in v equation (kg/(m*s)**2)
c       wmix     Total mixing in w equation (kg/(m*s)**2)
c       km       Turbulent mixing coefficient for momentum ( m**2/s ).
c
c     WORK ARRAYS:
c
c       real nsqed  Temporary array containing the Brunt Vaisala frquency
c       real tau11  Temporary work array
c       real tau12  Temporary work array
c       real tau13  Temporary work array
c       real tau22  Temporary work array
c       real tau23  Temporary work array
c       real tau33  Temporary work array
c       real tem1   Temporary work array
c       real 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
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions

      real u     (nx,ny,nz)     ! Total u-velocity (m/s)
      real v     (nx,ny,nz)     ! Total v-velocity (m/s)
      real w     (nx,ny,nz)     ! Total w-velocity (m/s)
      real ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
      real pprt  (nx,ny,nz)     ! Perturbation pressure (Pascal)
c
      real qv    (nx,ny,nz)     ! Water vapor specific humidity (kg/kg)
      real qc    (nx,ny,nz)     ! Cloud water mixing ratio (kg/kg)
      real qr    (nx,ny,nz)     ! Rain water mixing ratio (kg/kg)
      real qi    (nx,ny,nz)     ! Cloud ice mixing ratio (kg/kg)
      real qs    (nx,ny,nz)     ! Snow mixing ratio (kg/kg)
      real qh    (nx,ny,nz)     ! Hail mixing ratio (kg/kg)
c
      real ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
      real vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
      real ptbar (nx,ny,nz)     ! Base state potential temperature (K)
      real pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
      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)

      real usflx (nx,ny)        ! surface flux of u-momentum (kg/(m*s**2))
      real vsflx (nx,ny)        ! surface flux of v-momentum (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 umix  (nx,ny,nz)     ! Turbulent mixing on u-momentum.
      real vmix  (nx,ny,nz)     ! Turbulent mixing on v-momentum.
      real wmix  (nx,ny,nz)     ! Turbulent mixing on w-momentum.
      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )
c
      real nsqed (nx,ny,nz)     ! Temporary work array
      real tau11 (nx,ny,nz)     ! Temporary work array
      real tau12 (nx,ny,nz)     ! Temporary work array
      real tau13 (nx,ny,nz)     ! Temporary work array
      real tau22 (nx,ny,nz)     ! Temporary work array
      real tau23 (nx,ny,nz)     ! Temporary work array
      real tau33 (nx,ny,nz)     ! Temporary work array
      real tem1  (nx,ny,nz)     ! Temporary work array
      real tem2  (nx,ny,nz)     ! Temporary work array
c
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, j
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'phycst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external deform
      external cftmix
      external stress
      external divgu
      external divgv
      external divgw
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
C     Calculate the static stability parameter N squared (Brunt-Vaisala
c     frequency). tem1,tem2,tau11,tau12 are used as work arrays for 
c     this call.
c
c#######################################################################
c

      CALL stabnsq(nx,ny,nz,
     :             ptprt,pprt,qv,qc,qr,qi,qs,qh,ptbar,qvbar,pbar,j3,
     :             nsqed, tem1,tem2,tau11,tau12) 
c
c#######################################################################
c
c     Calculate the deformation tensor components Dij, which are stored
c     in arrays Tauij.
c
c#######################################################################
c
 
      CALL deform(nx,ny,nz,u,v,w,
     :            tau11,tau12,tau13,tau22,tau23,tau33)
 
c
c#######################################################################
c
c     Calculate the coefficient of turbulent mixing km using 
c     Smagorinsky's formulation
c
c#######################################################################
c

      CALL cftmix(nx,ny,nz,tau11,tau12,tau13,tau22,tau23,tau33,
     :            nsqed, km, tem1)

c
c#######################################################################
c
c     Calculate the stress tensor Tauij = Km * rhobar * Dij.
c
c#######################################################################
c

      CALL stress(nx,ny,nz,
     :            km,rhobar,tau11,tau12,tau13,tau22,tau23,tau33, tem1)

c
c
c#######################################################################
c
c
c     Set the surface momentum fluxes:
c
c     If sfcphy = 1, the surface momentum fluxes are set to the values in
c     usflx and vsflx, which are calcualted in the surface physics 
c     subroutine sfcflx. 
c
c     If sfcphy = 0, the fluxes are calculated in the same way as in the 
c     interior domain, i.e., using subgrid scale turbulence 
c     parameterization.
c
c#######################################################################
c

      IF( sfcphy.ne.0 ) THEN 
c
c#######################################################################
c
c     tau13 at the surface = usflx.
c
c#######################################################################
c
        DO 150 j=2,ny-2
        DO 150 i=2,nx-1
          tau13(i,j,2) = usflx(i,j)
150     CONTINUE
c  
c#######################################################################
c
c     tau23 at the surface = vsflx.
c
c#######################################################################
c
        DO 250 j=2,ny-1
        DO 250 i=2,nx-2
          tau23(i,j,2) = vsflx(i,j)
250     CONTINUE

      ENDIF
c
c#######################################################################
c
c     Calculate the divergence of stresses, which is the turbulent 
c     mixing on momentum u, v and w - umix, vmix and wmix.
c
c#######################################################################
c
      CALL divgu(nx,ny,nz, tau11,tau12,tau13, umix, tem1, tem2)

      CALL divgv(nx,ny,nz, tau12,tau22,tau23, vmix, tem1, tem2)

      CALL divgw(nx,ny,nz, tau13,tau23,tau33, wmix, tem1, tem2)

      RETURN
      END

C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE STABNSQ                    ######
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 STABNSQ(nx,ny,nz,
     :           ptprt,pprt,qv,qc,qr,qi,qs,qh,ptbar,qvbar,pbar,j3,nsqed, 
     :           tmprtr,qvs,qw,pt) 
c
c#######################################################################
c
c     PURPOSE:
c
C     Calculate the static stability parameter N squared (Brunt-Vaisala
c     frequency squared).  
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     9/15/1992 
c
c     MODIFICATION HISTORY:
c
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 (K)
c       pprt     Perturbation pressure (Pascal)
c       qv       Water vapor specific humidity (kg/kg)
c       qc       Cloud water mixing ratio (kg/kg)
c       qr       Rainwater mixing ratio (kg/kg)
c       qi       Cloud ice mixing ratio (kg/kg)
c       qs       Snow mixing ratio (kg/kg)
c       qh       Hail mixing ratio (kg/kg)
c       qvbar    Base state water vapor specific humidity (kg/kg)
c       pbar     Base state pressure (Pascal) 
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c     OUTPUT:
c
c       nsqed    Brunt-Vaisala frequency squared (1/s**2), a local array
c
c     WORK ARRAYS:
c
c       tmprtr   Work array for temperature
c       qvs      work array for saturation specific humidity
c       qw       Work array for total water mixing ratio
c       pt       Work arrays for total potential temperature 
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
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
      real pprt  (nx,ny,nz)     ! Perturbation pressure (Pascal)
      real ptbar (nx,ny,nz)     ! Base state potential temperature (K)
      real qv    (nx,ny,nz)     ! Water vapor specific humidity (kg/kg)
      real qc    (nx,ny,nz)     ! Cloud water mixing ratio (kg/kg)
      real qr    (nx,ny,nz)     ! Rain water mixing ratio (kg/kg)
      real qi    (nx,ny,nz)     ! Cloud ice mixing ratio (kg/kg)
      real qs    (nx,ny,nz)     ! Snow mixing ratio (kg/kg)
      real qh    (nx,ny,nz)     ! Hail mixing ratio (kg/kg)
      real qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
      real pbar  (nx,ny,nz)     ! Base state pressure (Pascal) 
      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! d( zp )/d( z ).

      real nsqed (nx,ny,nz)     ! Brunt-Vaisala frequency squared (1/s**2)

      real tmprtr(nx,ny,nz)     ! Temp. work array
      real qvs   (nx,ny,nz)     ! Temp. work array
      real qw    (nx,ny,nz)     ! Temp. work array
      real pt    (nx,ny,nz)     ! Temp. work array
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, j, k
c
      real ppi, dzinvd2
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'phycst.inc'
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Calculate the static stability n**2
c     The moist static stability follows that of Durran and Klemp (1982)
c
c#######################################################################
c
      DO 3000 k=1,nz-1
      DO 3000 j=1,ny-1
      DO 3000 i=1,nx-1

        ppi   = ((pbar(i,j,k)+pprt(i,j,k))/p0) ** (rd/cp)
        tmprtr(i,j,k) = (ptbar(i,j,k)+ptprt(i,j,k)) * ppi
        pt(i,j,k) =  ptbar(i,j,k)+ptprt(i,j,k)

3000  CONTINUE

      CALL satmr(nx,ny,nz,pbar, tmprtr, qvs)

      DO 3010 k=1,nz-1
      DO 3010 j=1,ny-1
      DO 3010 i=1,nx-1

        qw(i,j,k) = qvs(i,j,k)+qc(i,j,k)+qr(i,j,k)

3010  CONTINUE

      dzinvd2=1.0/(2*dz)

      DO 3100 k=2,nz-2
      DO 3100 j=1,ny-1
      DO 3100 i=1,nx-1

        IF( moist.eq.0. or. qc(i,j,k).le.0.0 ) THEN

          nsqed(i,j,k)=g*(pt(i,j,k+1)-pt(i,j,k-1))*dzinvd2
     :                 /(pt(i,j,k)*j3(i,j,k))

        ELSE

          nsqed(i,j,k)= g/j3(i,j,k) *
     :    ( (1+lathv*qvs(i,j,k)/(rd*tmprtr(i,j,k)))/
     :      (1+rd/rv*lathv*lathv*qvs(i,j,k)/(cp*rd*tmprtr(i,j,k)**2))
     :    * ((pt(i,j,k+1)-pt(i,j,k-1))/pt(i,j,k)
     :    +  lathv/(cp*tmprtr(i,j,k))*(qvs(i,j,k+1)-qvs(i,j,k-1)) )
     :    - (qw(i,j,k+1)-qw(i,j,k-1)) ) *dzinvd2

        ENDIF

3100  CONTINUE     
   
      DO 3200 j=1,ny-1
      DO 3200 i=1,nx-1
        nsqed(i,j,  1 )=nsqed(i,j,2)
        nsqed(i,j,nz-1)=nsqed(i,j,nz-2)
3200  CONTINUE        

      RETURN
      END
C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CFTMIX                     ######
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 CFTMIX(nx,ny,nz,
     :           tau11,tau12,tau13,tau22,tau23,tau33,nsqed, km, defsq)
c
c#######################################################################
c
c     PURPOSE:
c
C     Calculate the turbulent mixing coefficient km using modified 
c     Smagorinsky formulation.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
C
c     6/1/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       tau11    Deformation tensor component (1/s) 
c       tau12    Deformation tensor component (1/s)
c       tau13    Deformation tensor component (1/s)
c       tau22    Deformation tensor component (1/s)
c       tau23    Deformation tensor component (1/s)
c       tau33    Deformation tensor component (1/s)
c       nsqed    Brunt-Vaisala frequency (1/s**2), a local array
c
c     OUTPUT:
c
c       km       Turbulent mixing coefficient for momentum ( m**2/s )
c
c     WORK ARRAYS:
c
c       defsq    Deformation squared (1/s**2), a local 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
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real tau11(nx,ny,nz)      ! Deformation tensor component (1/s)
      real tau12(nx,ny,nz)      ! Deformation tensor component (1/s)
      real tau13(nx,ny,nz)      ! Deformation tensor component (1/s)
      real tau22(nx,ny,nz)      ! Deformation tensor component (1/s)
      real tau23(nx,ny,nz)      ! Deformation tensor component (1/s)
      real tau33(nx,ny,nz)      ! Deformation tensor component (1/s)

      real nsqed (nx,ny,nz)     ! Brunt-Vaisala frequency (1/s**2)
c
      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )
c
      real defsq (nx,ny,nz)     ! Deformation squared (1/s**2) 
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, j, k
c
      real tema, temb, temc
      real delta, coeff, prinv
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'phycst.inc'
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
c     Set the mixing coefficient km to a constant value TMIXCST.
c
c#######################################################################
c
      IF(tmixopt.eq.1) THEN 

      DO 50 k=1,nz-1
      DO 50 i=1,nx-1
      DO 50 j=1,ny-1
        km(i,j,k)=tmixcst
50    CONTINUE

      RETURN

      ENDIF
c
c#######################################################################
c
c     Calculate the magnitude of deformation Def**2.
c
c#######################################################################
c
      DO 2100 k=1,nz-1
      DO 2100 j=1,ny-1
      DO 2100 i=1,nx-1
        tema = (tau12(i,j,k)+tau12(i,j+1,k))
     :        +(tau12(i+1,j,k)+tau12(i+1,j+1,k))
        temb = (tau13(i,j,k)+tau13(i,j,k+1))
     :        +(tau13(i+1,j,k)+tau13(i+1,j,k+1))
        temc = (tau23(i,j,k)+tau23(i,j+1,k))
     :        +(tau23(i,j,k+1)+tau23(i,j+1,k+1))

        defsq(i,j,k) =(tau11(i,j,k)*tau11(i,j,k)
     :                +tau22(i,j,k)*tau22(i,j,k)
     :                +tau33(i,j,k)*tau33(i,j,k))*0.5
     :                +(tema*tema+temb*temb+temc*temc)*0.0625
2100  CONTINUE

c
c#######################################################################
c
c     Calculate the turbulent mixing coefficient km.
c     using modified Smagorinsky formulation.
c
c#######################################################################
c
      delta = grdscl 
      Coeff = ( 0.21 * delta )**2
      prinv = 1.0/prantl

      DO 3400 k=1,nz-1
      DO 3400 j=1,ny-1
      DO 3400 i=1,nx-1
         km(i,j,k)=coeff*
     :      sqrt( max(defsq(i,j,k)-nsqed(i,j,k)*prinv, 0.0) )
3400  CONTINUE        

      IF(tmixopt.eq.3) THEN 
c
c#######################################################################
c
c     Add to KM from Smagorinsky's formulation a constant.
c
c#######################################################################
c
      DO 60 k=1,nz-1
      DO 60 i=1,nx-1
      DO 60 j=1,ny-1
        km(i,j,k)=km(i,j,k)+tmixcst
60    CONTINUE

      ENDIF

      RETURN
      END

C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE DEFORM                     ######
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 DEFORM(nx,ny,nz,
     :           u,v,w,d11,d12,d13,d22,d23,d33)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the deformation tensor components Dij.
c
c     Dij = dv(i)/dx(j) + dv(j)/dx(i) - (2/3)*delta(i,j)dv(k)/dx(k),
c     where i,j,k represent summation subscripts and d indicates a partial
c     derivative.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
C
c     6/1/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity at a given time level (m/s)
c       v        y component of velocity at a given time level (m/s)
c       w        Vertical component of velocity in Cartesian
c                coordinates at a given time level (m/s).
c
c     OUTPUT:
c
c       d11      Deformation tensor component (1/s)
c       d12      Deformation tensor component (1/s)
c       d13      Deformation tensor component (1/s)
c       d21      Deformation tensor component (1/s)
c       d22      Deformation tensor component (1/s)
c       d23      Deformation tensor component (1/s)
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real u     (nx,ny,nz)     ! Total u-velocity (m/s)
      real v     (nx,ny,nz)     ! Total v-velocity (m/s)
      real w     (nx,ny,nz)     ! Total w-velocity (m/s)
c
      real d11   (nx,ny,nz)     ! Deformation tensor component (1/s)
      real d12   (nx,ny,nz)     ! Deformation tensor component (1/s)
      real d13   (nx,ny,nz)     ! Deformation tensor component (1/s)
      real d22   (nx,ny,nz)     ! Deformation tensor component (1/s)
      real d23   (nx,ny,nz)     ! Deformation tensor component (1/s)
      real d33   (nx,ny,nz)     ! Deformation tensor component (1/s)
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i,j,k
      real    dxinv2, dyinv2, dzinv2 
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'phycst.inc'
      include 'bndry.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external bcdefm
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
c     Calculate the deformation tensor components Dij, which are stored
c     in arrays dij.
c
c#######################################################################
c
      dxinv2 = 2/dx
      DO 1110 k=1,nz-1
      DO 1110 j=1,ny-1
      DO 1110 i=1,nx-1
        d11(i,j,k)=(u(i+1,j,k)-u(i,j,k))*dxinv2
1110  CONTINUE

      dyinv2 = 2/dy
      DO 1120 k=1,nz-1
      DO 1120 j=1,ny-1
      DO 1120 i=1,nx-1
        d22(i,j,k)=(v(i,j+1,k)-v(i,j,k))*dyinv2
1120  CONTINUE

      dzinv2 = 2/dz
      DO 1130 k=1,nz-1
      DO 1130 j=1,ny-1
      DO 1130 i=1,nx-1
        d33(i,j,k)=(w(i,j,k+1)-w(i,j,k))*dzinv2
1130  CONTINUE

      DO 1210 k=1,nz-1
      DO 1210 j=2,ny-1
      DO 1210 i=2,nx-1
        d12(i,j,k)=(u(i,j,k)-u(i,j-1,k))*dyinv
     :            +(v(i,j,k)-v(i-1,j,k))*dxinv
1210  CONTINUE


      DO 1220 k=2,nz-1
      DO 1220 j=1,ny-1
      DO 1220 i=2,nx-1
        d13(i,j,k)=(u(i,j,k)-u(i,j,k-1))*dzinv
     :            +(w(i,j,k)-w(i-1,j,k))*dxinv
1220  CONTINUE

      DO 1230 k=2,nz-1
      DO 1230 j=2,ny-1
      DO 1230 i=1,nx-1
        d23(i,j,k)=(v(i,j,k)-v(i,j,k-1))*dzinv
     :            +(w(i,j,k)-w(i,j-1,k))*dyinv
1230  CONTINUE
c
c#######################################################################
c
c     Boundary conditions for deformation tensor Dij, 
c
c#######################################################################
c 
      CALL bcdefm(nx,ny,nz,u,v,w, d12,d13,d23,
     :           ebc,wbc,nbc,sbc,tbc,bbc)

      RETURN
      END
C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE STRESS                     ######
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 STRESS(nx,ny,nz,
     :           km,rhobar,tau11,tau12,tau13,tau22,tau23,tau33, tem1)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the stress tensor tauij from the deformation tensor Dij 
c     (input as tauij) and the mixing coefficient KM.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
C
c     6/1/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       rhobar   Base state air density (kg/m**3)
c
c       km       Turbulent mixing coefficient for momentum ( m**2/s ).
c
c       tau11    Deformation tensor component (1/s) 
c       tau12    Deformation tensor component (1/s)
c       tau13    Deformation tensor component (1/s)
c       tau22    Deformation tensor component (1/s)
c       tau23    Deformation tensor component (1/s)
c       tau33    Deformation tensor component (1/s)
c
c
c     OUTPUT:
c
c       tau11    Stress tensor component (kg/(m*s**2)) 
c       tau12    Stress tensor component (kg/(m*s**2)) 
c       tau13    Stress tensor component (kg/(m*s**2)) 
c       tau22    Stress tensor component (kg/(m*s**2)) 
c       tau23    Stress tensor component (kg/(m*s**2)) 
c       tau33    Stress tensor component (kg/(m*s**2)) 
c
c     WORK ARRAYS:
c
c       tem1     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     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real tau11 (nx,ny,nz)     ! Deformation and stress tensor 
                                ! component (kg/(m*s**2))
      real tau12 (nx,ny,nz)     ! Deformation and stress tensor 
                                ! component (kg/(m*s**2))
      real tau13 (nx,ny,nz)     ! Deformation and stress tensor 
                                ! component (kg/(m*s**2))
      real tau22 (nx,ny,nz)     ! Deformation and stress tensor 
                                ! component (kg/(m*s**2))
      real tau23 (nx,ny,nz)     ! Deformation and stress tensor 
                                ! component (kg/(m*s**2))
      real tau33 (nx,ny,nz)     ! Deformation and stress tensor 
                                ! component (kg/(m*s**2))
c
      real tem1  (nx,ny,nz)     ! Temporary work array
c
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i,j,k
      real    tem 
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'phycst.inc'
c
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Calculate the stress tensor Tauij.
c
c#######################################################################
c
      DO 4100 k=1,nz-1
      DO 4100 j=1,ny-1
      DO 4100 i=1,nx-1
        tem1(i,j,k)=rhobar(i,j,k)*km(i,j,k)
4100  CONTINUE

      DO 4110 k=1,nz-1
      DO 4110 j=1,ny-1
      DO 4110 i=1,nx-1
        tau11(i,j,k)=tem1(i,j,k)*tau11(i,j,k)
4110  CONTINUE

      DO 4120 k=1,nz-1
      DO 4120 j=1,ny-1
      DO 4120 i=1,nx-1
        tau22(i,j,k)=tem1(i,j,k)*tau22(i,j,k)
4120  CONTINUE

      DO 4130 k=1,nz-1
      DO 4130 j=1,ny-1
      DO 4130 i=1,nx-1
        tau33(i,j,k)=tem1(i,j,k)*tau33(i,j,k)
4130  CONTINUE

      DO 4210 k=1,nz-1
      DO 4210 j=2,ny-1
      DO 4210 i=2,nx-1
        tem         =((tem1(i,j,k)+tem1(i-1,j,k))
     :              +(tem1(i,j-1,k)+tem1(i-1,j-1,k)))*0.25
        tau12(i,j,k)=tem*tau12(i,j,k)
4210  CONTINUE


      DO 4220 k=2,nz-1
      DO 4220 j=1,ny-1
      DO 4220 i=2,nx-1
        tem         =((tem1(i,j,k)+tem1(i-1,j,k))
     :               +(tem1(i,j,k-1)+tem1(i-1,j,k-1)))*0.25
        tau13(i,j,k)=tem*tau13(i,j,k)
4220  CONTINUE

      DO 4230 k=2,nz-1
      DO 4230 j=2,ny-1
      DO 4230 i=1,nx-1
        tem         =((tem1(i,j,k)+tem1(i,j-1,k))
     :              +(tem1(i,j,k-1)+tem1(i,j-1,k-1)))*0.25
        tau23(i,j,k)=tem*tau23(i,j,k)
4230  CONTINUE

      RETURN
      END
C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BCDEFM                     ######
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 BCDEFM(nx,ny,nz,u,v,w, d12,d13,d23,
     :           ebc,wbc,nbc,sbc,tbc,bbc)
c
c#######################################################################
c
c     PURPOSE:
c
c     Set the lateral boundary values of the deformation tensor Dij if i.ne.j.
c     For boundary option 2, periodic conditions are set, for all other
c     cases, zero gradient conditions are specified for Dij.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
C
c     6/1/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity at a given time level (m/s)
c       v        y component of velocity at a given time level (m/s)
c       w        Vertical component of velocity in Cartesian
c                coordinates at a given time level (m/s).
c
c       d12      Deformation tensor component (1/s) 
c       d13      Deformation tensor component (1/s)
c       d23      Deformation tensor component (1/s)
c
c     OUTPUT:
c
c       d12      Deformation tensor component including boundary values 
c       d13      Deformation tensor component including boundary values 
c       d23      Deformation tensor component including boundary values 
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      integer ebc,wbc,nbc,sbc,tbc,bbc
c
      real u     (nx,ny,nz)     ! Total u-velocity (m/s)
      real v     (nx,ny,nz)     ! Total v-velocity (m/s)
      real w     (nx,ny,nz)     ! Total w-velocity (m/s)

      real d12   (nx,ny,nz)     ! Deformation tensor component (1/s) 
      real d13   (nx,ny,nz)     ! Deformation tensor component (1/s)
      real d23   (nx,ny,nz)     ! Deformation tensor component (1/s)
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i,j,k
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
C     Boundary conditions for d12.
c
c#######################################################################
c 

c#######################################################################
c
c     X-boundary for d12.
c
c#######################################################################
c
      IF(ebc.eq.2) THEN 

        DO 110 k=1,nz-1
        DO 110 j=2,ny-1
          d12(1,j,k )=d12(nx-2,j,k)
110     CONTINUE

      ELSE

        DO 115 k=1,nz-1
        DO 115 j=2,ny-1
          d12(1,j,k )= - d12(3,j,k)
115     CONTINUE

      ENDIF

      IF(wbc.eq.2) THEN 

        DO 120 k=1,nz-1
        DO 120 j=2,ny-1
          d12(nx,j,k)=d12( 3 ,j,k)
120     CONTINUE

      ELSE

        DO 125 k=1,nz-1
        DO 125 j=2,ny-1
          d12(nx,j,k)=-d12(nx-2,j,k)
125     CONTINUE

      ENDIF

c
c#######################################################################
c
c     Y-boundary for d12.
c
c#######################################################################
c
      IF(sbc.eq.2) THEN 

        DO 130 k=1,nz-1
        DO 130 i=1,nx
          d12(i,1,k) =d12(i,ny-2,k)
130     CONTINUE

      ELSE

        DO 135 k=1,nz-1
        DO 135 i=1,nx
          d12(i,1,k) =-d12(i,3,k)
135     CONTINUE

      ENDIF

      IF(nbc.eq.2) THEN 

        DO 140 k=1,nz-1
        DO 140 i=1,nx
          d12(i,ny,k)=d12(i, 3 ,k)
140     CONTINUE

      ELSE

        DO 145 k=1,nz-1
        DO 145 i=1,nx
          d12(i,ny,k)=-d12(i,ny-2,k)
145     CONTINUE

      ENDIF
c
c
c#######################################################################
c
C     East-West and top-bottom boundary conditions for d13.
c
c#######################################################################
c

c#######################################################################
c
c     X b.c. for d13. 
c
c#######################################################################
c
      IF(ebc.eq.2) THEN 

        DO 210 j=1,ny-1
        DO 210 k=2,nz-1
          d13(1, j,k)=d13(nx-2,j,k)
210     CONTINUE

      ELSE

        DO 215 j=1,ny-1
        DO 215 k=2,nz-1
          d13(1, j,k)=-d13(3,j,k)
215     CONTINUE

      ENDIF

      IF(wbc.eq.2) THEN 

        DO 220 j=1,ny-1
        DO 220 k=2,nz-1
          d13(nx,j,k)=d13( 3 ,j,k)
220     CONTINUE

      ELSE

        DO 225 j=1,ny-1
        DO 225 k=2,nz-1
          d13(nx,j,k)=-d13(nx-2,j,k)
225     CONTINUE

      ENDIF
c
c#######################################################################
c
c     Top b.c. for d13
c
c#######################################################################
c
      IF(tbc.eq.2) THEN 

        DO 231 j=1,ny-1
        DO 231 i=1,nx
          d13(i,j,nz)=d13(i,j,3) 
231     CONTINUE

      ELSE

        DO 230 j=1,ny-1
        DO 230 i=1,nx
          d13(i,j,nz)=-d13(i,j,nz-2)
230     CONTINUE

      ENDIF
c 
c#######################################################################
c
c     Bottom b.c. for d13
c
c#######################################################################
c 
      IF(bbc.eq.2) THEN 

        DO 241 j=1,ny-1
        DO 241 i=1,nx
          d13(i,j,1 )=d13(i,j,nz-2)
241     CONTINUE

      ELSE

        DO 240 j=1,ny-1
        DO 240 i=1,nx
          d13(i,j,1 )=-d13(i,j,3)
240     CONTINUE

      ENDIF
c
c#######################################################################
c
C     Boundary conditions for d23.
c
c#######################################################################
c
c#######################################################################
c
c     Y b.c. for d23.
c
c#######################################################################
c
      IF(sbc.eq.2) THEN 

        DO 311 i=1,nx-1
        DO 311 k=2,nz-1
          d23(i,1 ,k)=d23(i,ny-2,k)
311     CONTINUE

      ELSE

        DO 312 i=1,nx-1
        DO 312 k=2,nz-1
          d23(i,1 ,k)=-d23(i,3,k)
312     CONTINUE

      ENDIF

      IF(nbc.eq.2) THEN 

        DO 321 i=1,nx-1
        DO 321 k=2,nz-1
          d23(i,ny,k)=d23(i, 3  ,k)
321     CONTINUE

      ELSE

        DO 322 i=1,nx-1
        DO 322 k=2,nz-1
          d23(i,ny,k)=-d23(i,ny-2,k)
322     CONTINUE

      ENDIF
c
c#######################################################################
c
c     Bottom b.c. for d23.
c
c#######################################################################
c
      IF(bbc.eq.2) THEN 

        DO 331 j=1,ny
        DO 331 i=1,nx-1
          d23(i,j,1 )=d23(i,j,nz-2)
331     CONTINUE

      ELSE

        DO 332 j=1,ny
        DO 332 i=1,nx-1
          d23(i,j,1 )=-d23(i,j,3)
332     CONTINUE

      ENDIF
c
c#######################################################################
c
c     Top b.c. for d23.
c
c#######################################################################
c
      IF(tbc.eq.2) THEN 

        DO 341 j=1,ny
        DO 341 i=1,nx-1
          d23(i,j,nz)=d23(i,j,3)
341     CONTINUE

      ELSE

        DO 342 j=1,ny
        DO 342 i=1,nx-1
          d23(i,j,nz)=-d23(i,j,nz-2)
342     CONTINUE

      ENDIF

      RETURN
      END
C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE TMIXPT                     ######
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 TMIXPT(nx,ny,nz, ptprt,ptbar,rhobar,km,ptsflx,
     :           x,y,z,zp, j1,j2,j3,
     :           ptmix,
     :           h1,h2,h3, tem1,tem2)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the turbulent mixing term for the potential temperature 
c     equation. The term is expressed in terms of turbulent heat fluxes. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
C
c     6/1/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       ptprt    Perturbation potential temperature at a given time level (K)
c       ptbar    Base state potential temperature (K)
c       rhobar   Base state air density (kg/m**3)
c
c       usflx    Surface flux of u-momentum
c       vsflx    Surface flux of v-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       km       Turbulent mixing coefficient for momentum ( m**2/s ).
c
c       ptsflx   Surface flux of heat (K*kg/(m*s**2))
c
c     OUTPUT:
c
c       ptmix    Total mixing in potential temperature
c                equation (K*kg/(m**3*s)).
c
c     WORK ARRAYS:
c
c       h1       Turbulent heat flux, a local array. 
c       h2       Turbulent heat flux, a local array.
c       h3       Turbulent heat flux, a local array.
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     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
      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 ptmix(nx,ny,nz)      ! Turbulent mixing on potential 
                                ! temperature (K*kg/(m**3 *s))
c
      real h1    (nx,ny,nz)     ! Turbulent heat flux 
      real h2    (nx,ny,nz)     ! Turbulent heat flux
      real h3    (nx,ny,nz)     ! Turbulent heat flux
      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
      real    tem
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'phycst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external divgs
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
c     Tem1 = ptbar+ptprt
c
c#######################################################################
c
      DO 1010 k=1,nz-1
      DO 1010 j=1,ny-1
      DO 1010 i=1,nx-1
        tem1(i,j,k)=ptbar(i,j,k)+ptprt(i,j,k)
1010  CONTINUE
c
c#######################################################################
c
c     Tem2 = rhobar*kh = rhobar*km/prantl
c
c#######################################################################
c
      tem = 1.0/prantl
      DO 1020 k=1,nz-1
      DO 1020 j=1,ny-1
      DO 1020 i=1,nx-1
        tem2(i,j,k)=rhobar(i,j,k)*km(i,j,k)*tem
1020  CONTINUE
c
c#######################################################################
c
c     Turbulent heat flux in x direction
c
c#######################################################################
c
      tem = 1.0/(2*dx)
      DO 1100 k=2,nz-2
      DO 1100 j=2,ny-2
      DO 1100 i=2,nx-1
        h1(i,j,k)=(tem2(i,j,k)+tem2(i-1,j,k))
     :           *(tem1(i,j,k)-tem1(i-1,j,k))*tem
1100  CONTINUE
c
c
c#######################################################################
c
c     Turbulent heat flux in y direction
c
c#######################################################################
c
      tem = 1.0/(2*dy)
      DO 1200 k=2,nz-2
      DO 1200 j=2,ny-1
      DO 1200 i=2,nx-2
        h2(i,j,k)=(tem2(i,j,k)+tem2(i,j-1,k))
     :           *(tem1(i,j,k)-tem1(i,j-1,k))*tem
1200  CONTINUE
c
c
c#######################################################################
c
c     Turbulent heat flux in z direction
c
c#######################################################################
c
      tem = 1.0/(2*dz)
      DO 1300 k=2,nz-1
      DO 1300 j=2,ny-2
      DO 1300 i=2,nx-2
        h3(i,j,k)=(tem2(i,j,k)+tem2(i,j,k-1))
     :           *(tem1(i,j,k)-tem1(i,j,k-1))*tem
1300  CONTINUE
c
c
c#######################################################################
c
c     The surface heat flux depends on the option chosen:
c
c     If sfcphy = 1, the surface heat flux is set to ptsflx, 
c     which is calculated in the surface physics subroutine sfcflx.
c
c     If sfcphy = 0, the flux is calculated in the same way as in the
c     interior, i.e., using subgrid scale turbulence parameterization.
c
c#######################################################################
c

      IF( sfcphy.ne.0 ) THEN
c
c#######################################################################
c
c     Heat flux h3 at the surface = ptsflx.
c
c#######################################################################
c
        DO 250 j=2,ny-2
        DO 250 i=2,nx-2
          h3(i,j,2) = ptsflx(i,j)
250     CONTINUE

      ENDIF

c
c#######################################################################
c
c     The divergence of heat flux, i.e. the turbulent mxing of heat.
c
c#######################################################################
c
      CALL divgs(nx,ny,nz, h1,h2,h3, ptmix, tem1,tem2)

      RETURN
      END
C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE TMIXQV                     ######
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 TMIXQV(nx,ny,nz, qv,rhobar,km,qvsflx,
     :           x,y,z,zp, j1,j2,j3,
     :           qvmix,
     :           h1,h2,h3, tem1,tem2)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the turbulent mixing term for any water substance equation. 
c     The term is expressed in terms of the turbulent flux of this water 
c     quantity.  
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue) 
c     Added full documentation.
C
c     6/1/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       qv       Water vapor specific humidity at a given time level (kg/kg)
c       rhobar   Base state air density (kg/m**3)
c
c       km       Turbulent mixing coefficient for momentum ( m**2/s ).
c
c       qvsflx  Surface flux of moisture (K*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     OUTPUT:
c
c       qvmix    Total mixing in water vapor equation (kg/(m**3 * s))
c
c     WORK ARRAYS:
c
c       h1       Turbulent flux of moisture. A local array. 
c       h2       Turbulent flux of moisture. A local array.
c       h3       Turbulent flux of moisture. A local array.
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     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real qv    (nx,ny,nz)     ! Water vapor specific humidity (kg/kg)
      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 qvsflx(nx,ny)        ! surface flux of moisture (kg/(m**2*s))

      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 qvmix(nx,ny,nz)      ! Turbulent mixing on water substance
                                ! kg/(m**3 *s)
c
      real h1    (nx,ny,nz)     ! Turbulent moisture flux.
      real h2    (nx,ny,nz)     ! Turbulent moisture flux.
      real h3    (nx,ny,nz)     ! Turbulent moisture flux.
      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
      real    tem
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c

      include 'globcst.inc'
      include 'phycst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external divgs
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
c     Tem1 = rhobar*kh = rhobar*km/prantl
c
c#######################################################################
c
      tem = 1.0/prantl
      DO 1020 k=1,nz-1
      DO 1020 j=1,ny-1
      DO 1020 i=1,nx-1
        tem1(i,j,k)=rhobar(i,j,k)*km(i,j,k)*tem
1020  CONTINUE
c
c#######################################################################
c
c     Turbulent heat flux in x direction
c
c#######################################################################
c
      tem = 1.0/(2*dx)
      DO 1100 k=2,nz-2
      DO 1100 j=2,ny-2
      DO 1100 i=2,nx-1
        h1(i,j,k)=(tem1(i,j,k)+tem1(i-1,j,k))
     :           *(qv(i,j,k)-qv(i-1,j,k))*tem
1100  CONTINUE
c
c#######################################################################
c
c     Turbulent heat flux in y direction
c
c#######################################################################
c
      tem = 1.0/(2*dy)
      DO 1200 k=2,nz-2
      DO 1200 j=2,ny-1
      DO 1200 i=2,nx-2
        h2(i,j,k)=(tem1(i,j,k)+tem1(i,j-1,k))
     :           *(qv(i,j,k)-qv(i,j-1,k))*tem
1200  CONTINUE
c
c#######################################################################
c
c     Turbulent heat flux in z direction
c
c#######################################################################
c
      tem = 1.0/(2*dz)
      DO 1300 k=2,nz-1
      DO 1300 j=2,ny-2
      DO 1300 i=2,nx-2
        h3(i,j,k)=(tem1(i,j,k)+tem1(i,j,k-1))
     :           *(qv(i,j,k)-qv(i,j,k-1))*tem
1300  CONTINUE
c
c
c#######################################################################
c
c     Set the surface moisture fluxes:
c
c     If sfcphy = 1, the surface moisture fluxes are set to ptsflx,
c     which is calcualted in the surface physics subroutine sfcflx.
c
c     If sfcphy = 0, the fluxes are calculated in the same way as in the
c     interior, i.e., using subgrid scale turbulence parameterization.
c
c#######################################################################
c

      IF( sfcphy.ne.0 ) THEN
c
c#######################################################################
c
c     Moisture flux h3 at the surface = ptsflx.
c
c#######################################################################
c
        DO 250 j=2,ny-2
        DO 250 i=2,nx-2
          h3(i,j,2) = qvsflx(i,j)
250     CONTINUE

      ENDIF

c
c#######################################################################
c
c     The divergence of heat fluxes, i.e. the turbulent mxing of heat.
c
c#######################################################################
c
      CALL divgs(nx,ny,nz, h1,h2,h3, qvmix, tem1,tem2)

      RETURN
      END
C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE TMIXQ                      ######
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 TMIXQ(nx,ny,nz, q,rhobar,km,
     :           x,y,z,zp, j1,j2,j3,
     :           qmix,
     :           h1,h2,h3, tem1,tem2)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the turbulent mixing term for any water substance equation. 
c     The term is expressed in terms of the turbulent fluxes of this water 
c     quantity.  
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue)
c     Added full documentation.
C
c     6/1/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       q        mixing ratio of one of the water/ice quantities at
c                a given time level (kg/kg)
c
c       rhobar   Base state air density (kg/m**3)
c
c       km       Turbulent mixing coefficient for momentum ( m**2/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     OUTPUT:
c
c       qmix     Total mixing in water/ice substance
c                equation (kg/(m**3 * s))
c
c     WORK ARRAYS:
c
c       h1       Turbulent flux a water quantity, a local array.
c       h2       Turbulent flux a water quantity, a local array.
c       h3       Turbulent flux a water quantity, a local array.
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 nx,ny,nz          ! Number of grid points in 3 directions
c
      real q     (nx,ny,nz)     ! Water/ice mixing ratio (kg/kg)
      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 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 qmix  (nx,ny,nz)     ! Turbulent mixing on water/ice 
                                ! substance (kg/(m**2 *s))
c
      real h1    (nx,ny,nz)     ! Turbulent flux of a water quantity
      real h2    (nx,ny,nz)     ! Turbulent flux of a water quantity
      real h3    (nx,ny,nz)     ! Turbulent flux of a water quantity
      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
      real    tem
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'phycst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external divgs
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Tem1 = rhobar*kh = rhobar*km/prantl
c
c#######################################################################
c
      tem = 1.0/prantl
      DO 1020 k=1,nz-1
      DO 1020 j=1,ny-1
      DO 1020 i=1,nx-1
        tem1(i,j,k)=rhobar(i,j,k)*km(i,j,k)*tem
1020  CONTINUE
c
c#######################################################################
c
c     Turbulent heat flux in x direction
c
c#######################################################################
c
      tem = 1.0/(2*dx)
      DO 1100 k=2,nz-2
      DO 1100 j=2,ny-2
      DO 1100 i=2,nx-1
        h1(i,j,k)=(tem1(i,j,k)+tem1(i-1,j,k))
     :           *(q(i,j,k)-q(i-1,j,k))*tem
1100  CONTINUE
c
c#######################################################################
c
c     Turbulent heat flux in y direction
c
c#######################################################################
c
      tem = 1.0/(2*dy)
      DO 1200 k=2,nz-2
      DO 1200 j=2,ny-1
      DO 1200 i=2,nx-2
        h2(i,j,k)=(tem1(i,j,k)+tem1(i,j-1,k))
     :           *(q(i,j,k)-q(i,j-1,k))*tem
1200  CONTINUE
c
c#######################################################################
c
c     Turbulent heat flux in z direction
c
c#######################################################################
c
      tem = 1.0/(2*dz)
      DO 1300 k=2,nz-1
      DO 1300 j=2,ny-2
      DO 1300 i=2,nx-2
        h3(i,j,k)=(tem1(i,j,k)+tem1(i,j,k-1))
     :           *(q(i,j,k)-q(i,j,k-1))*tem
1300  CONTINUE
c
c#######################################################################
c
c     The divergence of heat fluxes, i.e. the turbulent mxing of heat.
c
c#######################################################################
c
      CALL divgs(nx,ny,nz, h1,h2,h3, qmix, tem1,tem2)

      RETURN
      END
