
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CMIX2UVW                   ######
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 CMIX2UVW(nx,ny,nz,
     :           u,v,w, ubar,vbar,rhobar,
     :           umix,vmix,wmix,
     :           tem1,tem2,tem3)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the second order computational mixing terms for the momentum 
c     equations. Computational mixing is applied to velocity perturbation only.
c     These terms are placed in the arrays umix, vmix and wmix. 
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/3/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        Verticle component of Cartesian velocity at a
c                given time level (m/s)
c
c       ubar     Base state zonal velocity component (m/s)
c       vbar     Base state meridional velocity component (m/s)
c       rhobar   Base state density (kg/m**3)
c
c       umix     Array containing turbulent mixing on u 
c       vmix     Array containing turbulent mixing on v 
c       wmix     Array containing turbulent mixing on w 
c
c     OUTPUT:
c
c       umix     Turbulent and computational mixing on u.
c       vmix     Turbulent and computational mixing on v.
c       wmix     Turbulent and computational mixing on w.
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c       tem3     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

      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 ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
      real vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real umix  (nx,ny,nz)     ! Total mixing in u-eq.
      real vmix  (nx,ny,nz)     ! Total mixing in v-eq.
      real wmix  (nx,ny,nz)     ! Total mixing in w-eq.
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
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      integer onvf
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'phycst.inc'
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external difxx
      external difyy
      external difzz
      external avgz
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     IF the coefficients of computational mixing in both horizontal
c     and vertical direction are zero, exit this subroutine. 
c
c#######################################################################
c

      IF( cfcmh2.eq.0.0 .and. cfcmv2.eq.0.0 ) RETURN
c
c#######################################################################
c
c     Second order computational mixing on u'
c
c#######################################################################
c
      DO 100 k=1,nz-1
      DO 100 j=1,ny-1
      DO 100 i=1,nx
        tem1(i,j,k)=(u(i,j,k)-ubar(i,j,k))
     :              *rhobar(min(nx-1,i) ,j,k)
100   CONTINUE


      IF( cfcmh2 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the second order
c     horizontal computational mixing on u'
c
c#######################################################################
c
        CALL difxx(tem1, nx,ny,nz, 2,nx-1,2,ny-2,2,nz-2, dx, tem2)
c
        CALL difyy(tem1, nx,ny,nz, 2,nx-1,2,ny-2,2,nz-2, dy, tem3)

        DO 110 k=2,nz-2
        DO 110 j=2,ny-2
        DO 110 i=2,nx-1
          umix(i,j,k)=umix(i,j,k)+cfcmh2*(tem2(i,j,k)+tem3(i,j,k))
110     CONTINUE

      ENDIF
c
      IF( cfcmv2 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the second order 
c     vertical computational mixing on u'
c
c#######################################################################
c
        CALL difzz(tem1, nx,ny,nz, 2,nx-1,2,ny-2,2,nz-2, dz, tem2)

        DO 120 k=2,nz-2
        DO 120 j=2,ny-2
        DO 120 i=2,nx-1
          umix(i,j,k)=umix(i,j,k)+cfcmv2*tem2(i,j,k)
120     CONTINUE

      ENDIF
c
c#######################################################################
c
c     Second order computational mixing on v'
c
c#######################################################################
c
      DO 200 k=1,nz-1
      DO 200 j=1,ny
      DO 200 i=1,nx-1
        tem1(i,j,k)=(v(i,j,k)-vbar(i,j,k))
     :              *rhobar(i ,min(ny-1,j) ,k)
200   CONTINUE


      IF( cfcmh2 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the horizontal 
c     computational mixing on v'
c
c#######################################################################
c
        CALL difxx(tem1, nx,ny,nz, 2,nx-2,2,ny-1,2,nz-2, dx, tem2)
c
        CALL difyy(tem1, nx,ny,nz, 2,nx-2,2,ny-1,2,nz-2, dy, tem3)

        DO 210 k=2,nz-2
        DO 210 j=2,ny-1
        DO 210 i=2,nx-2
          vmix(i,j,k)=vmix(i,j,k)+cfcmh2*(tem2(i,j,k)+tem3(i,j,k))
210     CONTINUE

      ENDIF
c
      IF( cfcmv2 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the vertical 
c     computational mixing on v'
c
c#######################################################################

        CALL difzz(tem1, nx,ny,nz, 2,nx-2,2,ny-1,2,nz-2, dz, tem2)

        DO 220 k=2,nz-2
        DO 220 j=2,ny-1
        DO 220 i=2,nx-2
          vmix(i,j,k)=vmix(i,j,k)+cfcmv2*tem2(i,j,k)
220     CONTINUE

      ENDIF

c
c#######################################################################
c
c     Second order computational mixing on w
c
c#######################################################################
c
      onvf = 1
      CALL avgz(rhobar, onvf,
     :      nx,ny,nz, 1,nx-1, 1,ny-1, 2,nz-1, tem2)

      DO 305 j=1,ny-1
      DO 305 i=1,nx-1
        tem2(i,j,1) = tem2(i,j,3)
        tem2(i,j,nz)= tem2(i,j,nz-2)
305   CONTINUE

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


      IF( cfcmh2 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the horizontal 
c     computational mixing on w
c
c#######################################################################
c
        CALL difxx(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-1, dx, tem2)
c
        CALL difyy(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-1, dy, tem3)

        DO 310 k=2,nz-1
        DO 310 j=2,ny-2
        DO 310 i=2,nx-2
          wmix(i,j,k)=wmix(i,j,k)+cfcmh2*(tem2(i,j,k)+tem3(i,j,k))
310     CONTINUE

      ENDIF
c
      IF( cfcmv2 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the vertical 
c     computational mixing on w
c
c#######################################################################
c
        CALL difzz(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-1, dz, tem2)

        DO 320 k=2,nz-1
        DO 320 j=2,ny-2
        DO 320 i=2,nx-2
          wmix(i,j,k)=wmix(i,j,k)+cfcmv2*tem2(i,j,k)
320     CONTINUE

      ENDIF

   
      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CMIX2PT                    ######
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 CMIX2PT(nx,ny,nz, ptprt,rhobar, ptmix, tem1,tem2,tem3)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the second order computational mixing term in the potential
c     temperature equation. The mixing is applied to the potential temperature 
c     perturbation only. This term is placed in the array ptmix. 
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/3/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        Verticle component of Cartesian velocity at a
c                given time level (m/s)

c       ptprt    Perturbation potential temperature at a given
c                time level (K)
c       rhobar   Base state density (kg/m**3)
c
c       ptmix    Array containing turbulent mixing on potential 
c                temperature (kg*K/(m**3 * s).
c 
c     OUTPUT:
c
c       ptmix    Turbulent and computational mixing on potential
c                temperature (kg*K/(m**3 * s).
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c       tem3     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 ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
c
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real ptmix (nx,ny,nz)     ! Total mixing in potential temperature
                                ! equation (kg*K/(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
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'
      include 'phycst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external aamult
      external difxx
      external difyy
      external difzz
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c
c#######################################################################
c
c     If the coefficients of computatioanl mixing in both horizontal
c     and vertical are zero, exit this subroutine. 
c
c#######################################################################
c

      IF( cfcmh2.eq.0.0 .and. cfcmv2.eq.0.0 ) RETURN
c
c
c#######################################################################
c
c     Second order computational mixing on ptprt
c
c#######################################################################
c
      CALL aamult(rhobar,ptprt, nx,ny,nz,
     :            1,nx-1,1,ny-1,1,nz-1, tem3)
c

      IF( cfcmh2 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the horizontal 
c     computational mixing on pt
c
c#######################################################################
c
        CALL difxx(tem3, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dx, tem1)

c
        CALL difyy(tem3, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dy, tem2)

        DO 110 k=2,nz-2
        DO 110 j=2,ny-2
        DO 110 i=2,nx-2
          ptmix(i,j,k)=ptmix(i,j,k)+cfcmh2*(tem1(i,j,k)+tem2(i,j,k))
110     CONTINUE

      ENDIF
c
      IF( cfcmv2 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the vertical 
c     computational mixing on pt
c
c#######################################################################
c
        CALL difzz(tem3, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dz, tem1)

        DO 120 k=2,nz-2
        DO 120 j=2,ny-2
        DO 120 i=2,nx-2
          ptmix(i,j,k)=ptmix(i,j,k)+cfcmv2*tem1(i,j,k)
120     CONTINUE

      ENDIF

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CMIX2Q                     ######
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 CMIX2Q(nx,ny,nz, q ,rhobar, qmix, tem1,tem2,tem3)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the second order computational mixing term for each of the 
c     water/ice variables. This term is placed in the array qmix. 
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/3/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 for one of the water/ice variables (kg/kg)
c
c       rhobar   Base state density (kg/m**3)
c
c       qmix     Array containing turbulent mixing on water/ice variable.
c
c     OUTPUT:
c
c       qmix     Turbulent and computational mixing on water/ice variable.
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c       tem3     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 q     (nx,ny,nz)     ! mixing ratio for one of the water/ice
                                ! variables (kg/kg)
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real qmix  (nx,ny,nz)     ! Total mixing in water/ice variable eq.
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
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'
      include 'phycst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external aamult
      external difxx
      external difyy
      external difzz
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     If the coefficients of computatioanl mixing in both horizontal
c     and vertical are zero, exit this subroutine. 
c
c#######################################################################
c
      IF( cfcmh2.eq.0.0 .and. cfcmv2.eq.0.0 ) RETURN

      CALL aamult(rhobar,q,nx,ny,nz, 1,nx-1,1,ny-1,1,nz-1, tem3)
c
c#######################################################################
c
c     Second order computational mixing on q 
c
c#######################################################################

      IF( cfcmh2 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the horizontal 
c     computational mixing on q
c
c#######################################################################
c
        CALL difxx(tem3, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dx, tem1)
c
        CALL difyy(tem3, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dy, tem2)

        DO 110 k=2,nz-2
        DO 110 j=2,ny-2
        DO 110 i=2,nx-2
          qmix(i,j,k)=qmix(i,j,k)+cfcmh2*(tem1(i,j,k)+tem2(i,j,k))
110     CONTINUE

      ENDIF
c
      IF( cfcmv2 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the vertical 
c     computational mixing on q 
c
c#######################################################################
c
        CALL difzz(tem3, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dz, tem1)

        DO 120 k=2,nz-2
        DO 120 j=2,ny-2
        DO 120 i=2,nx-2
          qmix(i,j,k)=qmix(i,j,k)+cfcmv2*tem1(i,j,k)
120     CONTINUE

      ENDIF

      RETURN
      END
C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CMIX4UVW                   ######
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 CMIX4UVW(nx,ny,nz,
     :           u,v,w, ubar,vbar,rhobar,
     :           umix,vmix,wmix,
     :           tem1,tem2,tem3,tem4)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the fourth order computational mixing terms for the momentum 
c     equations. These terms are placed in the arrays umix, vmix and wmix 
c     and are applied to momentum perturbations only.
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/3/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        Verticle component of Cartesian velocity at a
c                given time level (m/s)
c
c       ubar     Base state zonal velocity component (m/s)
c       vbar     Base state meridional velocity component (m/s)
c       rhobar   Base state density (kg/m**3)
c
c       umix     Array containing turbulent mixing on u (kg/(m*s)**2)
c       vmix     Array containing turbulent mixing on v (kg/(m*s)**2)
c       wmix     Array containing turbulent mixing on w (kg/(m*s)**2)
c
c     OUTPUT:
c
c       umix     Turbulent and computational mixing on u (kg/(m*s)**2)
c       vmix     Turbulent and computational mixing on v (kg/(m*s)**2)
c       wmix     Turbulent and computational mixing on w (kg/(m*s)**2)
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
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)
c
      real ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
      real vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real umix  (nx,ny,nz)     ! Total mixing in u-eq. (kg/(m*s)**2)
      real vmix  (nx,ny,nz)     ! Total mixing in v-eq. (kg/(m*s)**2)
      real wmix  (nx,ny,nz)     ! Total mixing in w-eq. (kg/(m*s)**2)
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
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer onvf
      integer i, j, k
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 difxx
      external difyy
      external difzz
      external budifxx
      external budifyy
      external budifzz
      external bvdifxx
      external bvdifyy
      external bvdifzz
      external bwdifxx
      external bwdifyy
      external bwdifzz
      external avgz
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     If the coefficients of computational mixing in both horizontal
c     and vertical are zero, exit this subroutine. 
c
c#######################################################################
c
      IF( cfcmh4.eq.0.0 .and. cfcmv4.eq.0.0 ) RETURN
c
c#######################################################################
c
c     Fourth order computational mixing on u'
c
c#######################################################################
c
      DO 100 k=1,nz-1
      DO 100 j=1,ny-1
      DO 100 i=1,nx
        tem1(i,j,k)=(u(i,j,k)-ubar(i,j,k))
     :              *rhobar(min(nx-1,i) ,j,k)
100   CONTINUE


      IF( cfcmh4 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the horizontal 
c     computational mixing on u'
c
c#######################################################################
c
        CALL difxx(tem1, nx,ny,nz, 2,nx-1,2,ny-2,2,nz-2, dx, tem2)

        CALL budifxx(tem2, nx,ny,nz, ebc, wbc)

        CALL difxx(tem2, nx,ny,nz, 2,nx-1,2,ny-2,2,nz-2, dx, tem3)
c
        CALL difyy(tem1, nx,ny,nz, 2,nx-1,2,ny-2,2,nz-2, dy, tem2)

        CALL budifyy(tem2, nx,ny,nz, nbc, sbc)

        CALL difyy(tem2, nx,ny,nz, 2,nx-1,2,ny-2,2,nz-2, dy, tem4)


        DO 110 k=2,nz-2
        DO 110 j=2,ny-2
        DO 110 i=2,nx-1
          umix(i,j,k)=umix(i,j,k)-cfcmh4*(tem3(i,j,k)+tem4(i,j,k))
110     CONTINUE

      ENDIF
c
      IF( cfcmv4 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the vertical 
c     computational mixing on u'
c
c#######################################################################
c
        CALL difzz(tem1, nx,ny,nz, 2,nx-1,2,ny-2,2,nz-2, dz, tem2)

        CALL budifzz(tem2, nx,ny,nz, tbc, bbc)

        CALL difzz(tem2, nx,ny,nz, 2,nx-1,2,ny-2,2,nz-2, dz, tem3)

        DO 120 k=2,nz-2
        DO 120 j=2,ny-2
        DO 120 i=2,nx-1
          umix(i,j,k)=umix(i,j,k)-cfcmv4*tem3(i,j,k)
120     CONTINUE

      ENDIF
c
c#######################################################################
c
c     Fourth order computational mixing on v'
c
c#######################################################################
c
      DO 200 k=1,nz-1
      DO 200 j=1,ny
      DO 200 i=1,nx-1
        tem1(i,j,k)=(v(i,j,k)-vbar(i,j,k))
     :              *rhobar(i ,min(ny-1,j) ,k)
200   CONTINUE


      IF( cfcmh4 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the horizontal 
c     computational mixing on v'
c
c#######################################################################
c
        CALL difxx(tem1, nx,ny,nz, 2,nx-2,2,ny-1,2,nz-2, dx, tem2)

        CALL bvdifxx(tem2, nx,ny,nz, ebc, wbc)

        CALL difxx(tem2, nx,ny,nz, 2,nx-2,2,ny-1,2,nz-2, dx, tem3)

c
        CALL difyy(tem1, nx,ny,nz, 2,nx-2,2,ny-1,2,nz-2, dy, tem2)

        CALL bvdifyy(tem2, nx,ny,nz, nbc, sbc)

        CALL difyy(tem2, nx,ny,nz, 2,nx-2,2,ny-1,2,nz-2, dy, tem4)


        DO 210 k=2,nz-2
        DO 210 j=2,ny-1
        DO 210 i=2,nx-2
          vmix(i,j,k)=vmix(i,j,k)-cfcmh4*(tem3(i,j,k)+tem4(i,j,k))
210     CONTINUE

      ENDIF
c
      IF( cfcmv4 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the vertical 
c     computational mixing on v'
c
c#######################################################################
c
        CALL difzz(tem1, nx,ny,nz, 2,nx-2,2,ny-1,2,nz-2, dz, tem2)

        CALL bvdifzz(tem2, nx,ny,nz, tbc, bbc)

        CALL difzz(tem2, nx,ny,nz, 2,nx-2,2,ny-1,2,nz-2, dz, tem3)

        DO 220 k=2,nz-2
        DO 220 j=2,ny-1
        DO 220 i=2,nx-2
          vmix(i,j,k)=vmix(i,j,k)-cfcmv4*tem3(i,j,k)
220     CONTINUE

      ENDIF

c
c#######################################################################
c
c     Fourth order computational mixing on w
c
c#######################################################################
c
      onvf = 1
      CALL avgz(rhobar, onvf,
     :      nx,ny,nz, 1,nx-1, 1,ny-1, 2,nz-1, tem1)

      DO 305 j=1,ny-1
      DO 305 i=1,nx-1
        tem1(i,j,1) = tem1(i,j,3)
        tem1(i,j,nz)= tem1(i,j,nz-2)
305   CONTINUE

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


      IF( cfcmh4 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the horizontal 
c     computational mixing on w
c
c#######################################################################
c
        CALL difxx(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-1, dx, tem2)

        CALL bwdifxx(tem2, nx,ny,nz, ebc, wbc)

        CALL difxx(tem2, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-1, dx, tem3)
c
        CALL difyy(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-1, dy, tem2)

        CALL bwdifyy(tem2, nx,ny,nz, nbc, sbc)

        CALL difyy(tem2, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-1, dy, tem4)

        DO 310 k=2,nz-1
        DO 310 j=2,ny-2
        DO 310 i=2,nx-2
          wmix(i,j,k)=wmix(i,j,k)-cfcmh4*(tem3(i,j,k)+tem4(i,j,k))
310     CONTINUE

      ENDIF
c
      IF( cfcmv4 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the vertical 
c     computational mixing on w
c
c#######################################################################
c
        CALL difzz(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-1, dz, tem2)

        CALL bwdifzz(tem2, nx,ny,nz, tbc, bbc)

        CALL difzz(tem2, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-1, dz, tem3)

        DO 320 k=2,nz-1
        DO 320 j=2,ny-2
        DO 320 i=2,nx-2
          wmix(i,j,k)=wmix(i,j,k)-cfcmv4*tem3(i,j,k)
320     CONTINUE

      ENDIF

   
      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CMIX4PT                    ######
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 CMIX4PT(nx,ny,nz, ptprt,rhobar, ptmix,
     :           tem1,tem2,tem3,tem4)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the fourth order computational mixing term for the potential
c     temperature equation. Computational mixing is applied to the potential
c     temperature perturbation only and is placed in the array ptmix. 
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/3/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
c                time level (K)
c       rhobar   Base state density (kg/m**3)
c
c       ptmix    Array containing turbulent mixing on potential
c                temperature.
c
c     OUTPUT:
c
c       ptmix    Turbulent and computational mixing on potential
c                temperature.
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
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)
c
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real ptmix (nx,ny,nz)     ! Total mixing in potential temperature
                                ! equation.
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
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'
      include 'phycst.inc'
      include 'bndry.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external aamult
      external difxx
      external difyy
      external difzz
      external bsdifxx
      external bsdifyy
      external bsdifzz
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     If the coefficients of computatioanl mixing in both horizontal
c     and vertical are zero, exit this subroutine. 
c
c#######################################################################
c

      IF( cfcmh4.eq.0.0 .and. cfcmv4.eq.0.0 ) RETURN

      CALL aamult(rhobar,ptprt, nx,ny,nz,
     :            1,nx-1,1,ny-1,1,nz-1, tem1)
c
c#######################################################################
c
c     Fourth order computational mixing on pt
c
c#######################################################################

      IF( cfcmh4 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the horizontal 
c     computational mixing on pt
c
c#######################################################################
c
        CALL difxx(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dx, tem2)

        CALL bsdifxx(tem2, nx,ny,nz, ebc, wbc)

        CALL difxx(tem2, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dx, tem3)

c
        CALL difyy(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dy, tem2)

        CALL bsdifyy(tem2, nx,ny,nz, nbc, sbc)

        CALL difyy(tem2, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dy, tem4)


        DO 110 k=2,nz-2
        DO 110 j=2,ny-2
        DO 110 i=2,nx-2
          ptmix(i,j,k)=ptmix(i,j,k)-cfcmh4*(tem3(i,j,k)+tem4(i,j,k))
110     CONTINUE

      ENDIF
c
      IF( cfcmv4 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the vertical 
c     computational mixing on pt
c
c#######################################################################
c
        CALL difzz(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dz, tem2)

        CALL bsdifzz(tem2, nx,ny,nz, tbc, bbc)

        CALL difzz(tem2, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dz, tem3)

        DO 120 k=2,nz-2
        DO 120 j=2,ny-2
        DO 120 i=2,nx-2
          ptmix(i,j,k)=ptmix(i,j,k)-cfcmv4*tem3(i,j,k)
120     CONTINUE

      ENDIF

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CMIX4Q                     ######
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 CMIX4Q(nx,ny,nz, q,rhobar, qmix,
     :           tem1,tem2,tem3,tem4)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the fourth order computational mixing term for each of the 
c     water/ice variable equation. Computational mixing is placed in the
c     array qmix. 
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/3/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 for one of the water/ice variables (kg/kg) 
c
c       rhobar   Base state density (kg/m**3)
c
c       qmix     Array containing turbulent mixing on water/ice variable.
c
c     OUTPUT:
c
c       qmix     Turbulent and computational mixing on water/ice variable.
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
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 q     (nx,ny,nz)     ! mixing ratio for one of the water/ice
                                ! variables (kg/kg)
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real qmix  (nx,ny,nz)     ! Total mixing in water/ice variable eq.
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
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'
      include 'phycst.inc'
      include 'bndry.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external aamult
      external difxx
      external difyy
      external difzz
      external bsdifxx
      external bsdifyy
      external bsdifzz
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     If the coefficients of computatioanl mixing in both horizontal
c     and vertical are zero, exit this subroutine. 
c
c#######################################################################
c

      IF( cfcmh4.eq.0.0 .and. cfcmv4.eq.0.0 ) RETURN

      CALL aamult(rhobar,q, nx,ny,nz,
     :            1,nx-1,1,ny-1,1,nz-1, tem1)
c
c#######################################################################
c
c     Fourth order computational mixing on q 
c
c#######################################################################

      IF( cfcmh4 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the horizontal 
c     computational mixing on q 
c
c#######################################################################
c
        CALL difxx(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dx, tem2)

        CALL bsdifxx(tem2, nx,ny,nz, ebc, wbc)

        CALL difxx(tem2, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dx, tem3)

c
        CALL difyy(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dy, tem2)

        CALL bsdifyy(tem2, nx,ny,nz, nbc, sbc)

        CALL difyy(tem2, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dy, tem4)


        DO 110 k=2,nz-2
        DO 110 j=2,ny-2
        DO 110 i=2,nx-2
          qmix(i,j,k)=qmix(i,j,k)-cfcmh4*(tem3(i,j,k)+tem4(i,j,k))
110     CONTINUE

      ENDIF
c
      IF( cfcmv4 .ne. 0.0) THEN 
c
c#######################################################################
c
c     When the coefficient is not zero, calculate the vertical 
c     computational mixing on q 
c
c#######################################################################
c
        CALL difzz(tem1, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dz, tem2)

        CALL bsdifzz(tem2, nx,ny,nz, tbc, bbc)

        CALL difzz(tem2, nx,ny,nz, 2,nx-2,2,ny-2,2,nz-2, dz, tem3)

        DO 120 k=2,nz-2
        DO 120 j=2,ny-2
        DO 120 i=2,nx-2
          qmix(i,j,k)=qmix(i,j,k)-cfcmv4*tem3(i,j,k)
120     CONTINUE

      ENDIF

      RETURN
      END
