!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CMIX2UVW                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cmix2uvw(nx,ny,nz,                                           & 1
           u,v,w, ubar,vbar,rhostr,                                     &
           umix,vmix,wmix,                                              &
           tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the second order computational mixing terms for the momentum
!  equations. Computational mixing is applied to velocity perturbations
!  only. These terms are placed in the arrays umix, vmix and wmix.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/91
!
!  MODIFICATION HISTORY:
!
!  5/05/92 (M. Xue)
!  Added full documentation.
!
!  6/3/92 (M. Xue and H. Jin)
!  Further facelift.
!
!  4/2/92 (M. Xue and H. Jin)
!  Modify to include terrain
!
!  5/19/1998 (M. Xue)
!  Reformulated the computational mixing terms to move rhostr
!  outside the inner-most derivative.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    u        x component of velocity at a given time level (m/s)
!    v        y component of velocity at a given time level (m/s)
!    w        Verticle component of Cartesian velocity at a
!             given time level (m/s)
!
!    umix     Array containing turbulent mixing on u
!    vmix     Array containing turbulent mixing on v
!    wmix     Array containing turbulent mixing on w
!
!    ubar     Base state zonal velocity component (m/s)
!    vbar     Base state meridional velocity component (m/s)
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!
!  OUTPUT:
!
!    umix     Turbulent and computational mixing on u.
!    vmix     Turbulent and computational mixing on v.
!    wmix     Turbulent and computational mixing on w.
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  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 :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.

  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.

  REAL :: tem1  (nx,ny,nz)     ! Temporary work array
  REAL :: tem2  (nx,ny,nz)     ! Temporary work array
  REAL :: tem3  (nx,ny,nz)     ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j, k
  REAL :: temx,temy,temz
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'grid.inc'          ! Grid parameters
  INCLUDE 'globcst.inc'

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!  IF the coefficients of computational mixing in both horizontal
!  and vertical directions are zero, exit this subroutine.
!
!-----------------------------------------------------------------------
!
  IF( cfcmh2 == 0.0 .AND. cfcmv2 == 0.0 ) RETURN

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx
        tem3(i,j,k)=u(i,j,k)-ubar(i,j,k)
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the second order
!  horizontal computational mixing on u'
!
!-----------------------------------------------------------------------
!
  IF( cfcmh2 /= 0.0) THEN

    temx = cfcmh2/(dx*dx)
    temy = cfcmh2/(4.0*dy*dy)

    DO k=2,nz-2
      DO j=1,ny-1
        DO i=1,nx-1
          tem1(i,j,k)=(tem3(i+1,j,k)-tem3(i,j,k))*rhostr(i,j,k)*temx
        END DO
      END DO
    END DO

    DO k=2,nz-2
      DO j=2,ny-1
        DO i=2,nx-1
          tem2(i,j,k)=(tem3(i,j,k)-tem3(i,j-1,k))*                      &
              (rhostr(i,j  ,k)+rhostr(i-1,j  ,k)                        &
              +rhostr(i,j-1,k)+rhostr(i-1,j-1,k))*temy
        END DO
      END DO
    END DO
!
    DO k=2,nz-2
      DO i=2,nx-1
        tem2(i, 1,k) = 0.0
        tem2(i,ny,k) = 0.0
      END DO
    END DO

    DO k=2,nz-2
      DO j=1,ny-1
        DO i=2,nx-1
          umix(i,j,k)=umix(i,j,k)+(tem1(i,j,k)-tem1(i-1,j,k)            &
                                  +tem2(i,j+1,k)-tem2(i,j,k))
        END DO
      END DO
    END DO

  END IF
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the second order
!  vertical computational mixing on u'
!
!-----------------------------------------------------------------------
!
  IF( cfcmv2 /= 0.0) THEN
!
    temz = cfcmv2/(4.0*dz*dz)

    DO k=2,nz-1
      DO j=1,ny-1
        DO i=2,nx-1
          tem2(i,j,k)=(tem3(i,j,k)-tem3(i,j,k-1))*                      &
              (rhostr(i,j,k  )+rhostr(i-1,j,k  )                        &
              +rhostr(i,j,k-1)+rhostr(i-1,j,k-1))*temz
        END DO
      END DO
    END DO
!
    DO k=2,nz-2
      DO j=1,ny-1
        DO i=2,nx-1
          umix(i,j,k)=umix(i,j,k)+(tem2(i,j,k+1)-tem2(i,j,k))
        END DO
      END DO
    END DO

  END IF
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the second order
!  horizontal computational mixing on v'
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny
      DO i=1,nx-1
        tem3(i,j,k)=v(i,j,k)-vbar(i,j,k)
      END DO
    END DO
  END DO

  IF( cfcmh2 /= 0.0) THEN

    temx = cfcmh2/(4.0*dx*dx)
    temy = cfcmh2/(dy*dy)

    DO k=2,nz-2
      DO j=2,ny-1
        DO i=2,nx-1
          tem1(i,j,k)=(tem3(i,j,k)-tem3(i-1,j,k))*                      &
              (rhostr(i,j  ,k)+rhostr(i  ,j-1,k)                        &
              +rhostr(i-1,j,k)+rhostr(i-1,j-1,k))*temx
        END DO
      END DO
    END DO

    DO k=2,nz-2
      DO j=2,ny-1
        tem1(1, j,k) = 0.0
        tem1(nx,j,k) = 0.0
      END DO
    END DO

    DO k=2,nz-2
      DO j=1,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=(tem3(i,j+1,k)-tem3(i,j,k))*rhostr(i,j,k)*temy
        END DO
      END DO
    END DO
!
    DO k=2,nz-2
      DO j=2,ny-1
        DO i=1,nx-1
          vmix(i,j,k)=vmix(i,j,k)+(tem1(i+1,j,k)-tem1(i,j,k)            &
                                  +tem2(i,j,k)-tem2(i,j-1,k))
        END DO
      END DO
    END DO

  END IF
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the second order
!  vertical computational mixing on v'
!
!-----------------------------------------------------------------------
!
  IF( cfcmv2 /= 0.0) THEN

    temz = cfcmv2/(4.0*dz*dz)

    DO k=2,nz-1
      DO j=2,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=(tem3(i,j,k)-tem3(i,j,k-1))*                      &
              (rhostr(i,j,k  )+rhostr(i,j-1,k  )                        &
              +rhostr(i,j,k-1)+rhostr(i,j-1,k-1))*temz
        END DO
      END DO
    END DO
!
    DO k=2,nz-2
      DO j=2,ny-1
        DO i=1,nx-1
          vmix(i,j,k)=vmix(i,j,k)+(tem2(i,j,k+1)-tem2(i,j,k))
        END DO
      END DO
    END DO

  END IF
!
!-----------------------------------------------------------------------
!
!  Second order computational mixing on w
!
!-----------------------------------------------------------------------
!
  IF( cfcmh2 /= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the horizontal
!  computational mixing on w.
!
!-----------------------------------------------------------------------
!
    temx = cfcmh2/(4.0*dx*dx)
    temy = cfcmh2/(4.0*dy*dy)

    DO k=2,nz-1
      DO j=1,ny-1
        DO i=2,nx-1
          tem1(i,j,k)=(w(i,j,k)-w(i-1,j,k))                             &
                     *(rhostr(i,j  ,k)+rhostr(i-1,j  ,k)+               &
                       rhostr(i,j,k-1)+rhostr(i-1,j,k-1))*temx
        END DO
      END DO
    END DO

    DO k=2,nz-1
      DO j=1,ny-1
        tem1(1 ,j,k)=0.0
        tem1(nx,j,k)=0.0
      END DO
    END DO

    DO k=2,nz-1
      DO j=2,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=(w(i,j,k)-w(i,j-1,k))                             &
                     *(rhostr(i,j,  k)+rhostr(i,j-1,k)+                 &
                       rhostr(i,j,k-1)+rhostr(i,j-1,k-1))*temy
        END DO
      END DO
    END DO

    DO k=2,nz-1
      DO i=1,nx-1
        tem2(i,1 ,k)=0.0
        tem2(i,ny,k)=0.0
      END DO
    END DO

!-----------------------------------------------------------------------
!
!  Note that tem1 at i=1 and nx are zero, and tem2 at j=1,ny are zero,
!  equivalent to assuming zero normal gradient of s at the boundaries.
!
!-----------------------------------------------------------------------

    DO k=2,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          wmix(i,j,k)=wmix(i,j,k)+(tem1(i+1,j,k)-tem1(i,j,k)+           &
                                   tem2(i,j+1,k)-tem2(i,j,k))
        END DO
      END DO
    END DO

  END IF

  IF( cfcmv2 /= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the vertical
!  computational mixing on w.
!
!-----------------------------------------------------------------------
!
    temz = cfcmv2/(dz*dz)

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem1(i,j,k)=(w(i,j,k+1)-w(i,j,k))*rhostr(i,j,k)*temz
        END DO
      END DO
    END DO

    DO k=2,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          wmix(i,j,k)=wmix(i,j,k)+(tem1(i,j,k)-tem1(i,j,k-1))
        END DO
      END DO
    END DO

  END IF

  RETURN
END SUBROUTINE cmix2uvw

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CMIX2S                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cmix2s(nx,ny,nz, s ,rhostr, smix, tem1,tem2,tem3) 4
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Generic routine to calculate the second order computational mixing
!  for scalar s. The computational mixing is accumulated into array
!  smix on exit.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/91
!
!  MODIFICATION HISTORY:
!
!  5/05/92 (M. Xue)
!  Added full documentation.
!
!  6/3/92 (M. Xue and H. Jin)
!  Further facelift.
!
!  10/17/93 (M. Xue)
!  Bug fixes fix in line 3 of loop 111 and 112.
!
!  10/17/94 (M. Xue)
!  Subroutines CMIX2PT and CMIX2Q combined into a single generic
!  routine CMIX2S.
!
!  5/19/1998 (M. Xue)
!  Reformulated the computational mixing terms to move rhostr
!  outside the inner-most derivative.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    u        x component of velocity at a given time level (m/s)
!    v        y component of velocity at a given time level (m/s)
!    w        Verticle component of Cartesian velocity at a
!             given time level (m/s)
!
!    s        scalar field to which the compuational mixing is to
!             be applied.
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!
!    smix     Array containing turbulent mixing s on input.
!
!  OUTPUT:
!
!    smix     Turbulent mixing plus computational mixing on s.
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  REAL :: s     (nx,ny,nz)     ! Scalar to be mixed.

  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.

  REAL :: smix  (nx,ny,nz)     ! Total mixing on 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
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j, k
  REAL :: cfcmh2s,cfcmv2s, temx,temy,temz
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!

  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid parameters
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!
!-----------------------------------------------------------------------
!
!  If the coefficients of computational mixing in both the horizontal
!  and vertical are zero, exit this subroutine.
!
!-----------------------------------------------------------------------
!
  cfcmh2s = cfcmh2 * scmixfctr
  cfcmv2s = cfcmv2 * scmixfctr

  IF( cfcmh2s == 0.0 .AND. cfcmv2s == 0.0 ) RETURN
!
!-----------------------------------------------------------------------
!
!  Second order computational mixing on s.
!
!-----------------------------------------------------------------------
!
  IF( cfcmh2s /= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the horizontal
!  computational mixing.
!
!-----------------------------------------------------------------------
!
    temx = cfcmh2s/(2.0*dx*dx)
    temy = cfcmh2s/(2.0*dy*dy)


    DO k=1,nz-1
      DO j=1,ny-1
        DO i=2,nx-1
          tem1(i,j,k)=(s(i,j,k)-s(i-1,j,k))                             &
                     *(rhostr(i,j,k)+rhostr(i-1,j,k))*temx
        END DO
      END DO
    END DO

    DO k=1,nz-1
      DO j=1,ny-1
        tem1(1 ,j,k)=0.0
        tem1(nx,j,k)=0.0
      END DO
    END DO

    DO k=1,nz-1
      DO j=2,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=(s(i,j,k)-s(i,j-1,k))                             &
                     *(rhostr(i,j,k)+rhostr(i,j-1,k))*temy
        END DO
      END DO
    END DO

    DO k=1,nz-1
      DO i=1,nx-1
        tem2(i,1 ,k)=0.0
        tem2(i,ny,k)=0.0
      END DO
    END DO

!-----------------------------------------------------------------------
!
!  Note that tem1 at i=1 and nx are zero, and tem2 at j=1,ny are zero,
!  equivalent to assuming zero normal gradient of s at the boundaries.
!
!-----------------------------------------------------------------------

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          smix(i,j,k)=smix(i,j,k)+(tem1(i+1,j,k)-tem1(i,j,k)+           &
                                   tem2(i,j+1,k)-tem2(i,j,k))
        END DO
      END DO
    END DO

  END IF

  IF( cfcmv2s /= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the vertical
!  computational mixing.
!
!-----------------------------------------------------------------------
!
    temz = cfcmv2s/(2.0*dz*dz)

    DO k=2,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem1(i,j,k)=(s(i,j,k)-s(i,j,k-1))                             &
                     *(rhostr(i,j,k)+rhostr(i,j,k-1))*temz
        END DO
      END DO
    END DO

    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,1 )=0.0
        tem1(i,j,nz)=0.0
      END DO
    END DO

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          smix(i,j,k)=smix(i,j,k)+(tem1(i,j,k+1)-tem1(i,j,k))
        END DO
      END DO
    END DO

  END IF

  RETURN
END SUBROUTINE cmix2s
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CMIX4UVW                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cmix4uvw(nx,ny,nz,                                           & 1,60
           u,v,w, ubar,vbar,rhostr,                                     &
           umix,vmix,wmix,                                              &
           tem1,tem2,tem3,tem4,tem5)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the fourth order computational mixing terms for the momentum
!  equations. These terms are placed in the arrays umix, vmix and wmix
!  and are applied to momentum perturbations only.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/91
!
!  MODIFICATION HISTORY:
!
!  5/05/92 (M. Xue)
!  Added full documentation.
!
!  6/3/92 (M. Xue and H. Jin)
!  Further facelift.
!
!  4/2/92 (M. Xue and H. Jin)
!  Modify to include terrain
!
!  5/25/1998 (M. Xue)
!  Reformulated the computational mixing terms to move rhostr
!  outside the inner-most derivative.
!
!  11/9/2001 (M. Xue, D. Weber, and X. Jin)
!  Added monotonic computational mixing and 6-th order mixing option
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    u        x component of velocity at a given time level (m/s)
!    v        y component of velocity at a given time level (m/s)
!    w        Verticle component of Cartesian velocity at a
!             given time level (m/s)
!
!    ubar     Base state zonal velocity component (m/s)
!    vbar     Base state meridional velocity component (m/s)
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!
!    umix     Array containing turbulent mixing on u (kg/(m*s)**2)
!    vmix     Array containing turbulent mixing on v (kg/(m*s)**2)
!    wmix     Array containing turbulent mixing on w (kg/(m*s)**2)
!
!  OUTPUT:
!
!    umix     Turbulent and computational mixing on u (kg/(m*s)**2)
!    vmix     Turbulent and computational mixing on v (kg/(m*s)**2)
!    wmix     Turbulent and computational mixing on w (kg/(m*s)**2)
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem5     Temporary work array.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  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 :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.

  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)

  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
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j, k
  REAL :: temx, temy, temz
  INTEGER :: kount
  INTEGER :: mptag1,mptag2
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!

  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'        ! Grid parameters
  INCLUDE 'phycst.inc'
  INCLUDE 'bndry.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!-----------------------------------------------------------------------

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!  If the coefficients of computational mixing in both the horizontal
!  and vertical are zero, exit this subroutine.
!
!-----------------------------------------------------------------------
!
  IF( cfcmh4 == 0.0 .AND. cfcmv4 == 0.0 ) RETURN

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx
        tem5(i,j,k)=u(i,j,k)-ubar(i,j,k)
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Calculate ustr=rhostr*u', vstr=rhostr*v', wstr=rhostr*w'
!
!-----------------------------------------------------------------------
!


  IF( cfcmh4 /= 0.0 ) THEN

    temx = cfcmh4/(    dx**4)
    temy = cfcmh4/(4.0*dy**4)

    DO k=2,nz-2
      DO j=1,ny-1
        DO i=1,nx-1
          tem1(i,j,k)=(tem5(i+1,j,k)-tem5(i,j,k))*rhostr(i,j,k)*temx
        END DO
      END DO
    END DO

    DO k=2,nz-2
      DO j=2,ny-1
        DO i=2,nx-1
          tem2(i,j,k)=(tem5(i,j,k)-tem5(i,j-1,k))*                      &
              (rhostr(i,j  ,k)+rhostr(i-1,j  ,k)                        &
              +rhostr(i,j-1,k)+rhostr(i-1,j-1,k))*temy
        END DO
      END DO
    END DO
!
    DO k=2,nz-2
      DO i=2,nx-1
        tem2(i, 1,k) = 0.0
        tem2(i,ny,k) = 0.0
      END DO
    END DO

    DO k=2,nz-2
      DO j=1,ny-1
        DO i=2,nx-1
          tem3(i,j,k)=tem1(i,j,k)-tem1(i-1,j,k)
          tem4(i,j,k)=tem2(i,j+1,k)-tem2(i,j,k)
        END DO
      END DO
    END DO

!call test_dump (tem3,"XXXcmix_tem3",nx,ny,nz,1,0)
!call test_dump (tem4,"XXXcmix_tem4",nx,ny,nz,1,0)
    IF (mp_opt > 0) THEN
      CALL acct_interrupt(mp_acct)
      CALL mpsend2dew(tem3,nx,ny,nz,ebc,wbc,1,mptag1,tem1)
      CALL mpsend2dns(tem4,nx,ny,nz,nbc,sbc,1,mptag2,tem1)
      CALL mprecv2dew(tem3,nx,ny,nz,ebc,wbc,1,mptag1,tem1)
      CALL mprecv2dns(tem4,nx,ny,nz,nbc,sbc,1,mptag2,tem1)
    END IF
    CALL acct_interrupt(bc_acct)
    CALL budifxx(tem3, nx,ny,nz,1,ny-1,2,nz-2,ebc,wbc)
    CALL budifyy(tem4, nx,ny,nz,2,nx-1,2,nz-2,nbc,sbc)
    CALL acct_stop_inter
!call test_dump2(tem3,"XXXAcmix_tem3",nx,ny,nz,1,nx,2,ny-2,2,nz-2)
!call test_dump2(tem4,"XXXAcmix_tem4",nx,ny,nz,2,nx-1,1,ny-1,2,nz-2)

    IF(cmix_opt == 0) THEN     !  non-monotonic 4th order comp. mixing

      DO k=2,nz-2
        DO j=2,ny-2
          DO i=2,nx-1
            umix(i,j,k)=umix(i,j,k)-(                                   &
                (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k)) &
                +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 1) THEN  !  4th order monotonic computational mixing

      DO k=2,nz-2
        DO j=2,ny-2
          DO i=1,nx-1
            tem1(i,j,k)=(tem3(i+1,j,k)-tem3(i,j,k))
            IF(-tem1(i,j,k)*(tem5(i+1,j,k)-tem5(i,j,k)) < 0.0) tem1(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=2,nx-1
            tem2(i,j,k)=(tem4(i,j,k)-tem4(i,j-1,k))
            IF(-tem2(i,j,k)*(tem5(i,j,k)-tem5(i,j-1,k)) < 0.0) tem2(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=2,ny-2
          DO i=2,nx-1
            umix(i,j,k)=umix(i,j,k)-                                    &
                (tem1(i,j,k)-tem1(i-1,j,k)+tem2(i,j+1,k)-tem2(i,j,k))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 2 .OR.cmix_opt == 3 ) THEN  ! 6th order 
                                             ! = 2 6th order no mono...
                                             ! = 3 6th order mono...

      DO k=2,nz-2
        DO j=2,ny-2
          DO i=2,nx-1
            tem1(i,j,k)=(tem3(i+1,j,k)-tem3(i,j,k))                     &
                       -(tem3(i,j,k)-tem3(i-1,j,k))
            tem2(i,j,k)=(tem4(i,j+1,k)-tem4(i,j,k))                     &
                       -(tem4(i,j,k)-tem4(i,j-1,k))
          END DO
        END DO
      END DO

      IF (mp_opt > 0) THEN
        CALL acct_interrupt(mp_acct)
        CALL mpsend2dew(tem1,nx,ny,nz,ebc,wbc,1,mptag1,tem3)
        CALL mpsend2dns(tem2,nx,ny,nz,nbc,sbc,1,mptag2,tem3)
        CALL mprecv2dew(tem1,nx,ny,nz,ebc,wbc,1,mptag1,tem3)
        CALL mprecv2dns(tem2,nx,ny,nz,nbc,sbc,1,mptag2,tem3)
      END IF
      CALL acct_interrupt(bc_acct)
      CALL budifxx(tem1, nx,ny,nz,1,ny-1,2,nz-2,ebc,wbc)
      CALL budifyy(tem2, nx,ny,nz,2,nx-1,2,nz-2,nbc,sbc)
      CALL acct_stop_inter

      kount = 0
      DO k=2,nz-2
        DO j=2,ny-2
          DO i=1,nx-1
            tem3(i,j,k)=(tem1(i+1,j,k)-tem1(i,j,k))
            IF(cmix_opt == 3.AND.                                       &
                  tem3(i,j,k)*(tem5(i+1,j,k)-tem5(i,j,k)) < 0.0)THEN
              tem3(i,j,k)=0.0
              IF(j == 2) kount = kount + 1
            END IF
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=2,nx-1
            tem4(i,j,k)=(tem2(i,j,k)-tem2(i,j-1,k))
            IF(cmix_opt == 3.AND.                                       &
                tem4(i,j,k)*(tem5(i,j,k)-tem5(i,j-1,k)) < 0.0)          &
                tem4(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=2,ny-2
          DO i=2,nx-1
            umix(i,j,k)=umix(i,j,k)+                                    &
                (tem3(i,j,k)-tem3(i-1,j,k)+tem4(i,j+1,k)-tem4(i,j,k))
          END DO
        END DO
      END DO

    END IF

    !wdt update
    IF (sbc == 4) THEN ! do them for open BC only
      j = 1
      DO k=2,nz-2
        DO i=2,nx-1
          umix(i,j,k)=umix(i,j,k)-(                                     &
               (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))  &
              +(tem4(i,j+1,k)-tem4(i,j,k))- tem4(i,j,k))
        END DO
      END DO
    END IF
    IF (nbc == 4) THEN ! do them for open BC only
      j = ny-1
      DO k=2,nz-2
        DO i=2,nx-1
          umix(i,j,k)=umix(i,j,k)-(                                     &
               (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))  &
              +(             -tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)))
        END DO
      END DO
    END IF

  END IF

  IF( cfcmv4 /= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the vertical
!  computational mixing on u'
!
!-----------------------------------------------------------------------
!
    temz = cfcmv4/(4.0*dz**4)

    DO k=2,nz-1
      DO j=1,ny-1
        DO i=2,nx-1
          tem1(i,j,k)=(tem5(i,j,k)-tem5(i,j,k-1))*                      &
              (rhostr(i,j,k  )+rhostr(i-1,j,k  )                        &
              +rhostr(i,j,k-1)+rhostr(i-1,j,k-1))*temz
        END DO
      END DO
    END DO
!
    DO k=2,nz-2
      DO j=1,ny-1
        DO i=2,nx-1
          tem2(i,j,k)=tem1(i,j,k+1)-tem1(i,j,k)
        END DO
      END DO
    END DO

    CALL budifzz(tem2, nx,ny,nz,2,nx-1,1,ny-1, tbc, bbc)

    IF( cmix_opt == 0) THEN  ! no monotonic application.

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=2,nx-1
            umix(i,j,k)=umix(i,j,k)                                     &
                -((tem2(i,j,k+1)-tem2(i,j,k))-(tem2(i,j,k)-tem2(i,j,k-1)))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 1) THEN !  fourth order monotonic 

      DO k=2,nz-1
        DO j=1,ny-1
          DO i=2,nx-1
            tem1(i,j,k)=tem2(i,j,k)-tem2(i,j,k-1)
            IF(-tem1(i,j,k)*(tem5(i,j,k)-tem5(i,j,k-1)) < 0.0) tem1(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=2,nx-1
            umix(i,j,k)=umix(i,j,k)-(tem1(i,j,k+1)-tem1(i,j,k))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 2 .OR.cmix_opt == 3 ) THEN  ! 6th order 
                                             ! = 2 6th order no mono...
                                             ! = 3 6th order mono...

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=2,nx-1
            tem3(i,j,k)= (tem2(i,j,k+1)-tem2(i,j,k))                    &
                        -(tem2(i,j,k)-tem2(i,j,k-1))
          END DO
        END DO
      END DO
      CALL budifzz(tem3, nx,ny,nz,2,nx-1,1,ny-1, tbc, bbc)

      DO k=2,nz-1
        DO j=1,ny-1
          DO i=2,nx-1
            tem1(i,j,k)=tem3(i,j,k)-tem3(i,j,k-1)
            IF(cmix_opt == 3.AND.                                       &
                  tem1(i,j,k)*(tem5(i,j,k)-tem5(i,j,k-1)) < 0.0)THEN
              tem1(i,j,k)=0.0
              IF(j == 2)kount = kount+1
            END IF
          END DO
        END DO
      END DO

      PRINT*,'On-off switch active for u at ',                          &
              FLOAT(kount)/(2*(nx-2)*(nz-2))                            &
             ,' percent of times at t=', curtim

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=2,nx-1
            umix(i,j,k)=umix(i,j,k)+(tem1(i,j,k+1)-tem1(i,j,k))
          END DO
        END DO
      END DO

    END IF

  END IF

!
!-----------------------------------------------------------------------
!
!  Fourth-order computational mixing on v'.
!
!-----------------------------------------------------------------------

  DO k=1,nz-1
    DO j=1,ny
      DO i=1,nx-1
        tem5(i,j,k)=v(i,j,k)-vbar(i,j,k)
      END DO
    END DO
  END DO

  IF( cfcmh4 /= 0.0 ) THEN
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the horizontal
!  computational mixing on v'
!
!-----------------------------------------------------------------------
!

    temx = cfcmh4/(4.0*dx**4)
    temy = cfcmh4/(dy**4)

    DO k=2,nz-2
      DO j=2,ny-1
        DO i=2,nx-1
          tem1(i,j,k)=(tem5(i,j,k)-tem5(i-1,j,k))*                      &
              (rhostr(i,j  ,k)+rhostr(i  ,j-1,k)                        &
              +rhostr(i-1,j,k)+rhostr(i-1,j-1,k))*temx
        END DO
      END DO
    END DO

    DO k=2,nz-2
      DO j=2,ny-1
        tem1(1, j,k) = 0.0
        tem1(nx,j,k) = 0.0
      END DO
    END DO

    DO k=2,nz-2
      DO j=1,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=(tem5(i,j+1,k)-tem5(i,j,k))*rhostr(i,j,k)*temy
        END DO
      END DO
    END DO
!
    DO k=2,nz-2
      DO j=2,ny-1
        DO i=1,nx-1
          tem3(i,j,k)=tem1(i+1,j,k)-tem1(i,j,k)
          tem4(i,j,k)=tem2(i,j,k)-tem2(i,j-1,k)
        END DO
      END DO
    END DO


!call test_dump (tem3,"XXX1cmix_tem3",nx,ny,nz,2,0)
!call test_dump (tem4,"XXX1cmix_tem4",nx,ny,nz,2,0)
    IF (mp_opt > 0) THEN
      CALL acct_interrupt(mp_acct)
      CALL mpsend2dew(tem3,nx,ny,nz,ebc,wbc,2,mptag1,tem1)
      CALL mpsend2dns(tem4,nx,ny,nz,nbc,sbc,2,mptag2,tem1)
      CALL mprecv2dew(tem3,nx,ny,nz,ebc,wbc,2,mptag1,tem1)
      CALL mprecv2dns(tem4,nx,ny,nz,nbc,sbc,2,mptag2,tem1)
    END IF
    CALL acct_interrupt(bc_acct)
    CALL bvdifxx(tem3, nx,ny,nz,2,ny-1,2,nz-2,ebc, wbc)
    CALL bvdifyy(tem4, nx,ny,nz,1,nx-1,2,nz-2,nbc, sbc)
    CALL acct_stop_inter

    IF( cmix_opt == 0) THEN  !  4th order non-monotonic computational mixing

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=2,nx-2
            vmix(i,j,k)=vmix(i,j,k)-(                                   &
                (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k)) &
                +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 1) THEN  ! 4th order monotonic

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=2,nx-1
            tem1(i,j,k)=(tem3(i,j,k)-tem3(i-1,j,k))
            IF(-tem1(i,j,k)*(tem5(i,j,k)-tem5(i-1,j,k)) < 0.0) tem1(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=2,nx-2
            tem2(i,j,k)=(tem4(i,j+1,k)-tem4(i,j,k))
            IF(-tem2(i,j,k)*(tem5(i,j+1,k)-tem5(i,j,k)) < 0.0) tem2(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=2,nx-2
            vmix(i,j,k)=vmix(i,j,k)-                                    &
                 (tem1(i+1,j,k)-tem1(i,j,k)+tem2(i,j,k)-tem2(i,j-1,k))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 2 .OR.cmix_opt == 3 ) THEN  ! 6th order 
                                             ! = 2 6th order no mono...
                                             ! = 3 6th order mono...

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=2,nx-2
            tem1(i,j,k)=(tem3(i+1,j,k)-tem3(i,j,k))                     &
                       -(tem3(i,j,k)-tem3(i-1,j,k))
            tem2(i,j,k)=(tem4(i,j+1,k)-tem4(i,j,k))                     &
                       -(tem4(i,j,k)-tem4(i,j-1,k))
          END DO
        END DO
      END DO

      IF (mp_opt > 0) THEN
        CALL acct_interrupt(mp_acct)
        CALL mpsend2dew(tem1,nx,ny,nz,ebc,wbc,2,mptag1,tem3)
        CALL mpsend2dns(tem2,nx,ny,nz,nbc,sbc,2,mptag2,tem3)
        CALL mprecv2dew(tem1,nx,ny,nz,ebc,wbc,2,mptag1,tem3)
        CALL mprecv2dns(tem2,nx,ny,nz,nbc,sbc,2,mptag2,tem3)
      END IF
      CALL acct_interrupt(bc_acct)
      CALL bvdifxx(tem1, nx,ny,nz,2,ny-1,2,nz-2,ebc, wbc)
      CALL bvdifyy(tem2, nx,ny,nz,1,nx-1,2,nz-2,nbc, sbc)
      CALL acct_stop_inter

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=2,nx-1
            tem3(i,j,k)=(tem1(i,j,k)-tem1(i-1,j,k))
            IF(cmix_opt == 3.AND.                                       &
                tem3(i,j,k)*(tem5(i,j,k)-tem5(i-1,j,k)) < 0.0)          &
                tem3(i,j,k)=0.0
          END DO
        END DO
      END DO

!      kount = 0

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=2,nx-2
            tem4(i,j,k)=(tem2(i,j+1,k)-tem2(i,j,k))
            IF(cmix_opt == 3.AND.                                       &
                  tem4(i,j,k)*(tem5(i,j+1,k)-tem5(i,j,k)) < 0.0)THEN
              tem4(i,j,k)=0.0
!           if(j.eq.2) kount = kount + 1
            END IF
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=2,nx-2
            vmix(i,j,k)=vmix(i,j,k)+                                    &
                 (tem3(i+1,j,k)-tem3(i,j,k)+tem4(i,j,k)-tem4(i,j-1,k))
          END DO
        END DO
      END DO

    END IF

    !wdt update
    IF (ebc == 4) THEN ! Do them for open BC only
      i = nx-1
      DO k=2,nz-2
        DO j=2,ny-1
          vmix(i,j,k)=vmix(i,j,k)-(                                     &
              (             -tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))   &
              +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)))
        END DO
      END DO
    END IF
    IF (wbc == 4) THEN ! Do them for open BC only
      i = 1
      DO k=2,nz-2
        DO j=2,ny-1
          vmix(i,j,k)=vmix(i,j,k)-(                                     &
              (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k))                 &
              +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)))
        END DO
      END DO
    END IF

  END IF
!
  IF( cfcmv4 /= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the vertical
!  computational mixing on v'
!
!-----------------------------------------------------------------------
!
    temz = cfcmv4/(4.0*dz**4)

    DO k=2,nz-1
      DO j=2,ny-1
        DO i=1,nx-1
          tem1(i,j,k)=(tem5(i,j,k)-tem5(i,j,k-1))*                      &
              (rhostr(i,j,k  )+rhostr(i,j-1,k  )                        &
              +rhostr(i,j,k-1)+rhostr(i,j-1,k-1))*temz
        END DO
      END DO
    END DO
!
    DO k=2,nz-2
      DO j=2,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=tem1(i,j,k+1)-tem1(i,j,k)
        END DO
      END DO
    END DO

    CALL bvdifzz(tem2, nx,ny,nz,1,nx-1,2,ny-1,tbc, bbc)

    IF( cmix_opt == 0) THEN  !  no monotonic application

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=1,nx-1
            vmix(i,j,k)=vmix(i,j,k)                                     &
                -((tem2(i,j,k+1)-tem2(i,j,k))-(tem2(i,j,k)-tem2(i,j,k-1)))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 1) THEN  !  4th order mono...

      DO k=2,nz-1
        DO j=2,ny-1
          DO i=1,nx-1
            tem1(i,j,k)=tem2(i,j,k)-tem2(i,j,k-1)
            IF(-tem1(i,j,k)*(tem5(i,j,k)-tem5(i,j,k-1)) < 0.0) tem1(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=1,nx-1
            vmix(i,j,k)=vmix(i,j,k)-(tem1(i,j,k+1)-tem1(i,j,k))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 2 .OR.cmix_opt == 3 ) THEN  ! 6th order 
                                             ! = 2 6th order no mono...
                                             ! = 3 6th order mono...

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=1,nx-1
            tem3(i,j,k)=(tem2(i,j,k+1)-tem2(i,j,k))                     &
                       -(tem2(i,j,k)-tem2(i,j,k-1))
          END DO
        END DO
      END DO
      CALL bvdifzz(tem3, nx,ny,nz,1,nx-1,2,ny-1,tbc, bbc)

      DO k=2,nz-1
        DO j=2,ny-1
          DO i=1,nx-1
            tem1(i,j,k)=tem3(i,j,k)-tem3(i,j,k-1)
            IF(cmix_opt == 3.AND.                                       &
                  tem1(i,j,k)*(tem5(i,j,k)-tem5(i,j,k-1)) < 0.0)THEN
              tem1(i,j,k)=0.0
!           if(j.eq.2) kount = kount + 1
            END IF
          END DO
        END DO
      END DO

!      print*,'On-off switch active for v at ',
!    :            float(kount)/(2*(nx-1)*(nz-2))
!    :           ,' percent of times at t=', curtim


      DO k=2,nz-2
        DO j=2,ny-1
          DO i=1,nx-1
            vmix(i,j,k)=vmix(i,j,k)+(tem1(i,j,k+1)-tem1(i,j,k))
          END DO
        END DO
      END DO

    END IF

  END IF

  IF( cfcmh4 /= 0.0 ) THEN
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the horizontal
!  computational mixing on w
!
!-----------------------------------------------------------------------
!
    temx = cfcmh4/(4.0*dx**4)
    temy = cfcmh4/(4.0*dy**4)

    DO k=2,nz-1
      DO j=1,ny-1
        DO i=2,nx-1
          tem1(i,j,k)=(w(i,j,k)-w(i-1,j,k))                             &
                     *(rhostr(i,j  ,k)+rhostr(i-1,j  ,k)+               &
                       rhostr(i,j,k-1)+rhostr(i-1,j,k-1))*temx
        END DO
      END DO
    END DO

    DO k=2,nz-1
      DO j=1,ny-1
        tem1(1 ,j,k)=0.0
        tem1(nx,j,k)=0.0
      END DO
    END DO

    DO k=2,nz-1
      DO j=2,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=(w(i,j,k)-w(i,j-1,k))                             &
                     *(rhostr(i,j,  k)+rhostr(i,j-1,k)+                 &
                       rhostr(i,j,k-1)+rhostr(i,j-1,k-1))*temy
        END DO
      END DO
    END DO

    DO k=2,nz-1
      DO i=1,nx-1
        tem2(i,1 ,k)=0.0
        tem2(i,ny,k)=0.0
      END DO
    END DO

!-----------------------------------------------------------------------
!
!  Note that tem1 at i=1 and nx are zero, and tem2 at j=1,ny are zero,
!  equivalent to assuming zero normal gradient of s at the boundaries.
!
!-----------------------------------------------------------------------

    DO k=2,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem3(i,j,k)=tem1(i+1,j,k)-tem1(i,j,k)
          tem4(i,j,k)=tem2(i,j+1,k)-tem2(i,j,k)
        END DO
      END DO
    END DO

    IF (mp_opt > 0) THEN
      CALL acct_interrupt(mp_acct)
      CALL mpsend2dew(tem3,nx,ny,nz,ebc,wbc,0,mptag1,tem1)
      CALL mpsend2dns(tem4,nx,ny,nz,nbc,sbc,0,mptag2,tem1)
      CALL mprecv2dew(tem3,nx,ny,nz,ebc,wbc,0,mptag1,tem1)
      CALL mprecv2dns(tem4,nx,ny,nz,nbc,sbc,0,mptag2,tem1)
    END IF
    CALL acct_interrupt(bc_acct)
    CALL bwdifxx(tem3, nx,ny,nz,1,ny-1,2,nz-1,ebc, wbc)
    CALL bwdifyy(tem4, nx,ny,nz,1,nx-1,2,nz-1,nbc, sbc)
    CALL acct_stop_inter

    IF( cmix_opt == 0) THEN  !  no monotonic application

      DO k=2,nz-1
        DO j=2,ny-2
          DO i=2,nx-2
            wmix(i,j,k)=wmix(i,j,k)-(                                   &
                (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k)) &
                +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 1) THEN  !  4th order monotonic comp. mixing

      DO k=2,nz-1
        DO j=2,ny-2
          DO i=2,nx-1
            tem1(i,j,k)= tem3(i,j,k)-tem3(i-1,j,k)
            IF(-tem1(i,j,k)*(w(i,j,k)-w(i-1,j,k)) < 0.0) tem1(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-1
        DO j=2,ny-1
          DO i=2,nx-2
            tem2(i,j,k)= tem4(i,j,k)-tem4(i,j-1,k)
            IF(-tem2(i,j,k)*(w(i,j,k)-w(i,j-1,k)) < 0.0) tem2(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-1
        DO j=2,ny-2
          DO i=2,nx-2
            wmix(i,j,k)=wmix(i,j,k)-                                    &
                (tem1(i+1,j,k)-tem1(i,j,k)+tem2(i,j+1,k)-tem2(i,j,k))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 2 .OR.cmix_opt == 3 ) THEN  ! 6th order 
                                             ! = 2 6th order no mono...
                                             ! = 3 6th order mono...

      DO k=2,nz-1
        DO j=2,ny-2
          DO i=2,nx-2
            tem1(i,j,k)=(tem3(i+1,j,k)-tem3(i,j,k))                     &
                       -(tem3(i,j,k)-tem3(i-1,j,k))
            tem2(i,j,k)=(tem4(i,j+1,k)-tem4(i,j,k))                     &
                       -(tem4(i,j,k)-tem4(i,j-1,k))
          END DO
        END DO
      END DO
      IF (mp_opt > 0) THEN
        CALL acct_interrupt(mp_acct)
        CALL mpsend2dew(tem1,nx,ny,nz,ebc,wbc,0,mptag1,tem3)
        CALL mpsend2dns(tem2,nx,ny,nz,nbc,sbc,0,mptag2,tem3)
        CALL mprecv2dew(tem1,nx,ny,nz,ebc,wbc,0,mptag1,tem3)
        CALL mprecv2dns(tem2,nx,ny,nz,nbc,sbc,0,mptag2,tem3)
      END IF
      CALL acct_interrupt(bc_acct)
      CALL bwdifxx(tem1, nx,ny,nz,1,ny-1,2,nz-1,ebc, wbc)
      CALL bwdifyy(tem2, nx,ny,nz,1,nx-1,2,nz-1,nbc, sbc)
      CALL acct_stop_inter

      kount =0
      DO k=2,nz-1
        DO j=2,ny-2
          DO i=2,nx-1
            tem3(i,j,k)= tem1(i,j,k)-tem1(i-1,j,k)
            IF(cmix_opt == 3.AND. tem3(i,j,k)*(w(i,j,k)-w(i-1,j,k)) < 0.0) THEN
              tem3(i,j,k)=0.0
              IF(j == 2) kount = kount+1
            END IF
          END DO
        END DO
      END DO

      DO k=2,nz-1
        DO j=2,ny-1
          DO i=2,nx-2
            tem4(i,j,k)= tem2(i,j,k)-tem2(i,j-1,k)
            IF(cmix_opt == 3.AND. tem4(i,j,k)*(w(i,j,k)-w(i,j-1,k)) < 0.0) &
                tem4(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-1
        DO j=2,ny-2
          DO i=2,nx-2
            wmix(i,j,k)=wmix(i,j,k)+                                    &
                (tem3(i+1,j,k)-tem3(i,j,k)+tem4(i,j+1,k)-tem4(i,j,k))
          END DO
        END DO
      END DO

    END IF

    !wdt update
    IF (ebc == 4) THEN
      i = nx-1
      DO k=2,nz-1
        DO j=2,ny-2
          wmix(i,j,k)=wmix(i,j,k)-(                                     &
              (             -tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))   &
              +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
        END DO
      END DO
    END IF
    IF (wbc == 4) THEN
      i = 1
      DO k=2,nz-1
        DO j=2,ny-2
          wmix(i,j,k)=wmix(i,j,k)-(                                     &
              (tem3(i+1,j,k)-tem3(i,j,k))-tem3(i,j,k)                   &
              +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
        END DO
      END DO
    END IF
    IF (nbc == 4) THEN
      j = ny - 1
      DO k=2,nz-1
        DO i=2,nx-2
          wmix(i,j,k)=wmix(i,j,k)-(                                     &
              (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))   &
              +(             -tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
        END DO
      END DO
    END IF
    IF (sbc == 4) THEN
      j = 1
      DO k=2,nz-1
        DO i=2,nx-2
          wmix(i,j,k)=wmix(i,j,k)-(                                     &
              (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))   &
              +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)              ) )
        END DO
      END DO
    END IF
    IF (wbc == 4.AND.sbc == 4) THEN
      i = 1
      j = 1
      DO k=2,nz-1
        wmix(i,j,k)=wmix(i,j,k)-(                                       &
            (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k))                   &
            +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)) )
      END DO
    END IF
    IF (ebc == 4.AND.sbc == 4) THEN
      i = nx-1
      j = 1
      DO k=2,nz-1
        wmix(i,j,k)=wmix(i,j,k)-(                                       &
            (             -tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))     &
            +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)) )
      END DO
    END IF

    IF (wbc == 4.AND.nbc == 4) THEN
      i = 1
      j = ny-1
      DO k=2,nz-1
        wmix(i,j,k)=wmix(i,j,k)-(                                       &
            (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k))                   &
            +(             -tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
      END DO
    END IF
    IF (ebc == 4.AND.nbc == 4) THEN
      i = nx-1
      j = ny-1
      DO k=2,nz-1
        wmix(i,j,k)=wmix(i,j,k)-(                                       &
            (             -tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))     &
            +(             -tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
      END DO
    END IF

  END IF

  IF( cfcmv4 /= 0.0 ) THEN
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the vertical
!  computational mixing on w
!
!-----------------------------------------------------------------------
!
    temz = cfcmv4/(dz**4)

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem1(i,j,k)=(w(i,j,k+1)-w(i,j,k))*rhostr(i,j,k)*temz
        END DO
      END DO
    END DO

    DO k=2,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=tem1(i,j,k)-tem1(i,j,k-1)
        END DO
      END DO
    END DO

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

    IF(cmix_opt == 0) THEN  ! no monotonic application

      DO k=2,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            wmix(i,j,k)=wmix(i,j,k)-(                                   &
                (tem2(i,j,k+1)-tem2(i,j,k))-(tem2(i,j,k)-tem2(i,j,k-1)))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 1) THEN  !  4th order monotonic 

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem1(i,j,k)=(tem2(i,j,k+1)-tem2(i,j,k))
            IF(-tem1(i,j,k)*(w(i,j,k+1)-w(i,j,k)) < 0.0) tem1(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            wmix(i,j,k)=wmix(i,j,k)-(tem1(i,j,k)-tem1(i,j,k-1))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 2 .OR.cmix_opt == 3 ) THEN  ! 6th order 
                                             ! = 2 6th order no mono...
                                             ! = 3 6th order mono...

      DO k=2,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem3(i,j,k)=(tem2(i,j,k+1)-tem2(i,j,k))                     &
                       -(tem2(i,j,k)-tem2(i,j,k-1))
          END DO
        END DO
      END DO
      CALL bwdifzz(tem3, nx,ny,nz,1,nx-1,1,ny-1,tbc, bbc)

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem1(i,j,k)=(tem3(i,j,k+1)-tem3(i,j,k))
            IF(cmix_opt == 3.AND. tem1(i,j,k)*(w(i,j,k+1)-w(i,j,k)) < 0.0)THEN
              tem1(i,j,k)=0.0
              IF(j == 2) kount=kount+1
            END IF
          END DO
        END DO
      END DO

      PRINT*,'On-off switch active for w at ',                          &
              FLOAT(kount)/(2*(nx-1)*(nz-2))                            &
             ,' percent of times at t=', curtim


      DO k=2,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            wmix(i,j,k)=wmix(i,j,k)+(tem1(i,j,k)-tem1(i,j,k-1))
          END DO
        END DO
      END DO

    END IF

  END IF

  RETURN
END SUBROUTINE cmix4uvw

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CMIX4S                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cmix4s(nx,ny,nz, s ,rhostr, smix, tem1,tem2,tem3,tem4) 4,20
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Generic routine to calculate the fourth order computational mixing
!  for scalar s. The computational mixing is accumulated into array
!  smix on exit.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/91
!
!  MODIFICATION HISTORY:
!
!  5/05/92 (M. Xue)
!  Added full documentation.
!
!  6/3/92 (M. Xue and H. Jin)
!  Further facelift.
!
!  10/17/94 (M. Xue)
!  Subroutines CMIX4PT and CMIX4Q combined into a single generic
!  routine CMIX2S.
!
!  5/25/1998 (M. Xue)
!  Reformulated the computational mixing terms to move rhostr
!  outside the inner-most derivative.
!
!  11/9/2001 (M. Xue, D. Weber, and X. Jin)
!  Added monotonic computational mixing and 6-th order mixing option
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    s        scalar field to which the compuational mixing is to
!             be applied.
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!
!    smix     Array containing turbulent mixing on s.
!
!  OUTPUT:
!
!    smix     Turbulent mixing plus computational mixing on s.
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  REAL :: s     (nx,ny,nz)     ! Scalar to be mixed.

  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.

  REAL :: smix  (nx,ny,nz)     ! Total mixing on 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
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j,k
  REAL :: dxinv4, dyinv4
  REAL :: cfcmh4s,cfcmv4s, temx,temy,temz
  INTEGER :: kount
  INTEGER :: mptag1,mptag2
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!

  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid parameters
  INCLUDE 'phycst.inc'
  INCLUDE 'bndry.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!  If the coefficients of computational mixing in both the horizontal
!  and vertical are zero, exit this subroutine.
!
!-----------------------------------------------------------------------
!
  cfcmh4s = cfcmh4 * scmixfctr
  cfcmv4s = cfcmv4 * scmixfctr

  IF( cfcmh4s == 0.0 .AND. cfcmv4s == 0.0 ) RETURN

  dxinv4 = dxinv**4
  dyinv4 = dyinv**4

  IF( cfcmh4s /= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the 4th order horizontal
!  computational mixing on s.
!
!-----------------------------------------------------------------------
!
    temx = cfcmh4s/(2.0*dx**4)
    temy = cfcmh4s/(2.0*dy**4)


    DO k=1,nz-1
      DO j=1,ny-1
        DO i=2,nx-1
          tem1(i,j,k)=(s(i,j,k)-s(i-1,j,k))                             &
                     *(rhostr(i,j,k)+rhostr(i-1,j,k))*temx
        END DO
      END DO
    END DO

    DO k=1,nz-1
      DO j=1,ny-1
        tem1(1 ,j,k)=0.0
        tem1(nx,j,k)=0.0
      END DO
    END DO

    DO k=1,nz-1
      DO j=2,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=(s(i,j,k)-s(i,j-1,k))                             &
                     *(rhostr(i,j,k)+rhostr(i,j-1,k))*temy
        END DO
      END DO
    END DO

    DO k=1,nz-1
      DO i=1,nx-1
        tem2(i,1 ,k)=0.0
        tem2(i,ny,k)=0.0
      END DO
    END DO

!-----------------------------------------------------------------------
!
!  Note that tem1 at i=1 and nx are zero, and tem2 at j=1,ny are zero,
!  equivalent to assuming zero normal gradient of s at the boundaries.
!
!-----------------------------------------------------------------------

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem3(i,j,k)=tem1(i+1,j,k)-tem1(i,j,k)
          tem4(i,j,k)=tem2(i,j+1,k)-tem2(i,j,k)
        END DO
      END DO
    END DO

!call test_dump (tem3,"XXX3cmix_tem3",nx,ny,nz,0,0)
!call test_dump (tem4,"XXX3cmix_tem4",nx,ny,nz,0,0)
    IF (mp_opt > 0) THEN
      CALL acct_interrupt(mp_acct)
      CALL mpsend2dew(tem3,nx,ny,nz,ebc,wbc,0,mptag1,tem1)
      CALL mpsend2dns(tem4,nx,ny,nz,nbc,sbc,0,mptag2,tem1)
      CALL mprecv2dew(tem3,nx,ny,nz,ebc,wbc,0,mptag1,tem1)
      CALL mprecv2dns(tem4,nx,ny,nz,nbc,sbc,0,mptag2,tem1)
    END IF
    CALL acct_interrupt(bc_acct)
    CALL bsdifxx(tem3, nx,ny,nz,1,ny-1,2,nz-2,ebc, wbc)
    CALL bsdifyy(tem4, nx,ny,nz,1,nx-1,2,nz-2,nbc, sbc)
    CALL acct_stop_inter
!call test_dump2(tem3,"XXXA3cmix_tem3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(tem4,"XXXA3cmix_tem4",nx,ny,nz,2,nx-2,1,ny-1,2,nz-2)

    IF( cmix_opt == 0) THEN  !  no monotonic application

      DO k=2,nz-2
        DO j=2,ny-2
          DO i=2,nx-2
            smix(i,j,k)=smix(i,j,k)-(                                   &
                (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k)) &
                +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 1) THEN   ! 4th order monotonic.

      DO k=2,nz-2
        DO j=2,ny-2
          DO i=2,nx-1
            tem1(i,j,k)= tem3(i,j,k)-tem3(i-1,j,k)
            IF(-tem1(i,j,k)*(s(i,j,k)-s(i-1,j,k)) < 0.0) tem1(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=2,nx-2
            tem2(i,j,k)= tem4(i,j,k)-tem4(i,j-1,k)
            IF(-tem2(i,j,k)*(s(i,j,k)-s(i,j-1,k)) < 0.0) tem2(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=2,ny-2
          DO i=2,nx-2
            smix(i,j,k)=smix(i,j,k)-                                    &
                (tem1(i+1,j,k)-tem1(i,j,k)+tem2(i,j+1,k)-tem2(i,j,k))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 2 .OR.cmix_opt == 3 ) THEN  ! 6th order 
                                             ! = 2 6th order no mono...
                                             ! = 3 6th order mono...

      DO k=2,nz-2
        DO j=2,ny-2
          DO i=2,nx-2
            tem1(i,j,k)=(tem3(i+1,j,k)-tem3(i,j,k))                     &
                       -(tem3(i,j,k)-tem3(i-1,j,k))
            tem2(i,j,k)=(tem4(i,j+1,k)-tem4(i,j,k))                     &
                       -(tem4(i,j,k)-tem4(i,j-1,k))
          END DO
        END DO
      END DO

      IF (mp_opt > 0) THEN
        CALL acct_interrupt(mp_acct)
        CALL mpsend2dew(tem1,nx,ny,nz,ebc,wbc,0,mptag1,tem3)
        CALL mpsend2dns(tem2,nx,ny,nz,nbc,sbc,0,mptag2,tem3)
        CALL mprecv2dew(tem1,nx,ny,nz,ebc,wbc,0,mptag1,tem3)
        CALL mprecv2dns(tem2,nx,ny,nz,nbc,sbc,0,mptag2,tem3)
      END IF
      CALL acct_interrupt(bc_acct)
      CALL bsdifxx(tem1, nx,ny,nz,1,ny-1,2,nz-2,ebc, wbc)
      CALL bsdifyy(tem2, nx,ny,nz,1,nx-1,2,nz-2,nbc, sbc)
      CALL acct_stop_inter

      kount= 0
      DO k=2,nz-2
        DO j=2,ny-2
          DO i=2,nx-1
            tem3(i,j,k)= tem1(i,j,k)-tem1(i-1,j,k)
!           IF(scmix_opt == 1.AND.cmix_opt == 3.AND.                    &
            IF(cmix_opt == 3.AND.                    &
                  tem3(i,j,k)*(s(i,j,k)-s(i-1,j,k)) < 0.0)THEN
              tem3(i,j,k)=0.0
              IF(j == 2) kount=kount+1
            END IF
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=2,ny-1
          DO i=2,nx-2
            tem4(i,j,k)= tem2(i,j,k)-tem2(i,j-1,k)
!           IF(scmix_opt == 1.AND.cmix_opt == 3.AND.                    &
            IF(cmix_opt == 3.AND.                    &
                tem4(i,j,k)*(s(i,j,k)-s(i,j-1,k)) < 0.0)                &
                tem4(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=2,ny-2
          DO i=2,nx-2
            smix(i,j,k)=smix(i,j,k)+                                    &
                (tem3(i+1,j,k)-tem3(i,j,k)+tem4(i,j+1,k)-tem4(i,j,k))
          END DO
        END DO
      END DO

    END IF

!
!-----------------------------------------------------------------------
!
!  Calculate the mixing term on the boundary assuming that mixed
!  field has zero gradient outside the boundary.
!  The boundary values of smix are needed only for open boundary option.
!
!-----------------------------------------------------------------------
!
    !wdt update
    IF (wbc == 4) THEN
      i = 1
      DO k=2,nz-2
        DO j=2,ny-2
          smix(i,j,k)=smix(i,j,k)-(                                     &
              (tem3(i+1,j,k)-tem3(i,j,k))-tem3(i,j,k)                   &
              +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
        END DO
      END DO
    END IF
    IF (ebc == 4) THEN
      i = nx-1
      DO k=2,nz-2
        DO j=2,ny-2
          smix(i,j,k)=smix(i,j,k)-(                                     &
              (             -tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))   &
              +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
        END DO
      END DO
    END IF
    IF (sbc == 4) THEN
      j = 1
      DO k=2,nz-2
        DO i=2,nx-2
          smix(i,j,k)=smix(i,j,k)-(                                     &
              (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))   &
              +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)              ) )
        END DO
      END DO
    END IF
    IF (nbc == 4) THEN
      j = ny - 1
      DO k=2,nz-2
        DO i=2,nx-2
          smix(i,j,k)=smix(i,j,k)-(                                     &
              (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))   &
              +(             -tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
        END DO
      END DO
    END IF
    IF (wbc == 4.AND.sbc == 4) THEN
      i = 1
      j = 1
      DO k=2,nz-2
        smix(i,j,k)=smix(i,j,k)-(                                       &
            (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k))                   &
            +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)) )
      END DO
    END IF
    IF (ebc == 4.AND.sbc == 4) THEN
      i = nx-1
      j = 1
      DO k=2,nz-2
        smix(i,j,k)=smix(i,j,k)-(                                       &
            (             -tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))     &
            +(tem4(i,j+1,k)-tem4(i,j,k))-(tem4(i,j,k)) )
      END DO
    END IF
    IF (wbc == 4.AND.nbc == 4) THEN
      i = 1
      j = ny-1
      DO k=2,nz-2
        smix(i,j,k)=smix(i,j,k)-(                                       &
            (tem3(i+1,j,k)-tem3(i,j,k))-(tem3(i,j,k))                   &
            +(             -tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
      END DO
    END IF
    IF (ebc == 4.AND.nbc == 4) THEN
      i = nx-1
      j = ny-1
      DO k=2,nz-2
        smix(i,j,k)=smix(i,j,k)-(                                       &
            (             -tem3(i,j,k))-(tem3(i,j,k)-tem3(i-1,j,k))     &
            +(             -tem4(i,j,k))-(tem4(i,j,k)-tem4(i,j-1,k)) )
      END DO
    END IF

  END IF

  IF( cfcmv4s /= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  If the coefficient is not zero, calculate the vertical
!  computational mixing on s.
!
!-----------------------------------------------------------------------

    temz = cfcmv4s/(2.0*dz**4)

    DO k=2,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem1(i,j,k)=(s(i,j,k)-s(i,j,k-1))                             &
                     *(rhostr(i,j,k)+rhostr(i,j,k-1))*temz
        END DO
      END DO
    END DO

    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,1 )=0.0
        tem1(i,j,nz)=0.0
      END DO
    END DO

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=tem1(i,j,k+1)-tem1(i,j,k)
        END DO
      END DO
    END DO

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

    IF( cmix_opt == 0) THEN  !  no mono applied

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=1,nx-1
            smix(i,j,k)=smix(i,j,k)                                     &
                -((tem2(i,j,k+1)-tem2(i,j,k))-(tem2(i,j,k)-tem2(i,j,k-1)))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 1) THEN  !  4th order mono

      DO k=2,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem1(i,j,k)=tem2(i,j,k)-tem2(i,j,k-1)
            IF(-tem1(i,j,k)*(s(i,j,k)-s(i,j,k-1)) < 0.)tem1(i,j,k)=0.0
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=1,nx-1
            smix(i,j,k)=smix(i,j,k)-(tem1(i,j,k+1)-tem1(i,j,k))
          END DO
        END DO
      END DO

    ELSE IF( cmix_opt == 2 .OR.cmix_opt == 3 ) THEN  ! 6th order 
                                             ! = 2 6th order no mono...
                                             ! = 3 6th order mono...

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=1,nx-1
            tem3(i,j,k)=(tem2(i,j,k+1)-tem2(i,j,k))                     &
                       -(tem2(i,j,k)-tem2(i,j,k-1))
          END DO
        END DO
      END DO
      CALL bsdifzz(tem3, nx,ny,nz,1,nx-1,1,ny-1,tbc, bbc)

      DO k=2,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem1(i,j,k)=tem3(i,j,k)-tem3(i,j,k-1)
!           IF(scmix_opt == 1.AND.cmix_opt == 3.AND.                    &
            IF(cmix_opt == 3.AND.                    &
                  tem1(i,j,k)*(s(i,j,k)-s(i,j,k-1)) < 0.)THEN
              tem1(i,j,k)=0.0
              IF(j == 2) kount=kount+1
            END IF
          END DO
        END DO
      END DO

      PRINT*,'On-off switch active for pt at ',                         &
              FLOAT(kount)/(2*(nx-1)*(nz-2))                            &
             ,' percent of times at t=', curtim


      DO k=2,nz-2
        DO j=1,ny-1
          DO i=1,nx-1
            smix(i,j,k)=smix(i,j,k)+(tem1(i,j,k+1)-tem1(i,j,k))
          END DO
        END DO
      END DO

    END IF

  END IF

  RETURN
END SUBROUTINE cmix4s