! WRF:MODEL_LAYER:PHYSICS
 
    MODULE module_diffusion_em

    USE module_configure
    USE module_bc
    USE module_state_description
    USE module_big_step_utilities_em
    USE module_model_constants    
    USE module_wrf_error

    CONTAINS

!=======================================================================
!=======================================================================

    SUBROUTINE cal_deform_and_div( config_flags, u, v, w, div,       &
                                   defor11, defor22, defor33,        &
                                   defor12, defor13, defor23,        &
                                   nba_rij, n_nba_rij,               & !JDM
                                   u_base, v_base, msfux, msfuy,     &
                                   msfvx, msfvy, msftx, msfty,       &
                                   rdx, rdy, dn, dnw, rdz, rdzw,     &
                                   fnm, fnp, cf1, cf2, cf3, zx, zy,  &
                                   ids, ide, jds, jde, kds, kde,     &
                                   ims, ime, jms, jme, kms, kme,     &
                                   its, ite, jts, jte, kts, kte      )

! History:     Sep 2003  Changes by Jason Knievel and George Bryan, NCAR
!              Oct 2001  Converted to mass core by Bill Skamarock, NCAR
!              ...        ...

! Purpose:     This routine calculates deformation and 3-d divergence.

! References:  Klemp and Wilhelmson (JAS 1978)
!              Chen and Dudhia (NCAR WRF physics report 2000)

!-----------------------------------------------------------------------
! Comments 10-MAR-05
! Equations 13a-f, Chen and Dudhia 2000, Appendix A:
! Eqn 13a: D11=defor11= 2m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
! Eqn 13b: D22=defor22= 2m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
! Eqn 13c: D33=defor33= 2 * partial dw/dz [SIMPLER FORM]
! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
!                              partial dpsi/dx * partial dv^/dpsi +
!                              partial dpsi/dy * partial du^/dpsi)
! Eqn 13e: D13=defor13= m^2 * (partial dw^/dX + partial dpsi/dx * partial dw^/dpsi)
!                           + partial du/dz [SIMPLER FORM]
! Eqn 13f: D23=defor23= m^2 * (partial dw^/dY + partial dpsi/dy * partial dw^/dpsi)
!                           + partial dv/dz [SIMPLER FORM]
!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: ids, ide, jds, jde, kds, kde, &
       ims, ime, jms, jme, kms, kme, &
       its, ite, jts, jte, kts, kte

    REAL, INTENT( IN )  &
    :: rdx, rdy, cf1, cf2, cf3

    REAL, DIMENSION( kms:kme ), INTENT( IN )  &
    :: fnm, fnp, dn, dnw, u_base, v_base

    REAL, DIMENSION( ims:ime , jms:jme ),  INTENT( IN )  &
    :: msfux, msfuy, msfvx, msfvy, msftx, msfty

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
    ::  u, v, w, zx, zy, rdz, rdzw

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &
    :: defor11, defor22, defor33, defor12, defor13, defor23, div 

   INTEGER, INTENT(  IN ) :: n_nba_rij !JDM

   REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_rij), INTENT(INOUT) & !JDM
   :: nba_rij


! Local variables.

    INTEGER  &
    :: i, j, k, ktf, ktes1, ktes2, i_start, i_end, j_start, j_end

    REAL  &
    :: tmp, tmpzx, tmpzy, tmpzeta_z, cft1, cft2

    REAL, DIMENSION( its:ite, jts:jte )  &
    :: mm, zzavg, zeta_zd12

    REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 )  &
    :: tmp1, hat, hatavg

! End declarations.
!-----------------------------------------------------------------------

! Comments 10-MAR-2005
! Treat all differentials as 'div-style' [or 'curl-style'],
! i.e., du/dX becomes (in map coordinate space) mx*my * d(u/my)/dx,
! NB - all equations referred to here are from Chen and Dudhia 2002, from the
! WRF physics documents web pages:
! http://www.mmm.ucar.edu/wrf/users/docs/wrf-doc-physics.pdf

!=======================================================================
! In the following section, calculate 3-d divergence and the first three
! (defor11, defor22, defor33) of six deformation terms.

    ktes1   = kte-1
    ktes2   = kte-2

    cft2    = - 0.5 * dnw(ktes1) / dn(ktes1)
    cft1    = 1.0 - cft2

    ktf     = MIN( kte, kde-1 )

    i_start = its
    i_end   = MIN( ite, ide-1 )
    j_start = jts
    j_end   = MIN( jte, jde-1 )

! Square the map scale factor.

    DO j = j_start, j_end
    DO i = i_start, i_end
      mm(i,j) = msftx(i,j) * msfty(i,j)
    END DO
    END DO

!-----------------------------------------------------------------------
! Calculate du/dx.

! Apply a coordinate transformation to zonal velocity, u.

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end+1
      hat(i,k,j) = u(i,k,j) / msfuy(i,j)
    END DO
    END DO
    END DO

! Average in x and z.

    DO j=j_start,j_end
    DO k=kts+1,ktf
    DO i=i_start,i_end
      hatavg(i,k,j) = 0.5 *  &
                    ( fnm(k) * ( hat(i,k  ,j) + hat(i+1,  k,j) ) +  &
                      fnp(k) * ( hat(i,k-1,j) + hat(i+1,k-1,j) ) )
    END DO
    END DO
    END DO

! Extrapolate to top and bottom of domain (to w levels).

    DO j = j_start, j_end
    DO i = i_start, i_end
      hatavg(i,1,j)   =  0.5 * (  &
                         cf1 * hat(i  ,1,j) +  &
                         cf2 * hat(i  ,2,j) +  &
                         cf3 * hat(i  ,3,j) +  &
                         cf1 * hat(i+1,1,j) +  &
                         cf2 * hat(i+1,2,j) +  &
                         cf3 * hat(i+1,3,j) )
      hatavg(i,kte,j) =  0.5 * (  &
                        cft1 * ( hat(i,ktes1,j) + hat(i+1,ktes1,j) )  +  &
                        cft2 * ( hat(i,ktes2,j) + hat(i+1,ktes2,j) ) )
    END DO
    END DO

    ! Comments 10-MAR-05
    ! Eqn 13a: D11=defor11= 2m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
    ! Below, D11 is set = 2*tmp1
    ! => tmp1 = m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
    ! tmpzx = averaged value of dpsi/dx (=zx)

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tmpzx       = 0.25 * (  &
                    zx(i,k  ,j) + zx(i+1,k  ,j) +  &
                    zx(i,k+1,j) + zx(i+1,k+1,j) )
      tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) *tmpzx * rdzw(i,k,j)
      ! tmp1 to here = partial dpsi/dx * partial du^/dpsi:
    END DO
    END DO
    END DO

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tmp1(i,k,j) = mm(i,j) * ( rdx * ( hat(i+1,k,j) - hat(i,k,j) ) -  &
                    tmp1(i,k,j))
    END DO
    END DO
    END DO

! End calculation of du/dx.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Calculate defor11 (2*du/dx).
! Comments 10-MAR-05
! Eqn 13a: D11=defor11= 2 m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
!                     = 2*tmp1

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      defor11(i,k,j) = 2.0 * tmp1(i,k,j)
    END DO
    END DO
    END DO

! End calculation of defor11.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Calculate zonal divergence (du/dx) and add it to the divergence array.

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      div(i,k,j) = tmp1(i,k,j)
    END DO
    END DO
    END DO

! End calculation of zonal divergence.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Calculate dv/dy.

! Apply a coordinate transformation to meridional velocity, v.

    DO j = j_start, j_end+1
    DO k = kts, ktf
    DO i = i_start, i_end
      ! Because msfvx at the poles will be undefined (1./0.), we will have
      ! trouble.  But we are OK since v at the poles is 0., and that takes
      ! precedence in this case.
      IF ((config_flags%polar) .AND. ((j == jds) .OR. (j == jde))) THEN
         hat(i,k,j) = 0.
      ELSE ! normal code
      hat(i,k,j) = v(i,k,j) / msfvx(i,j)
      ENDIF
    END DO
    END DO
    END DO

! Account for the slope in y of eta surfaces.

    DO j=j_start,j_end
    DO k=kts+1,ktf
    DO i=i_start,i_end
      hatavg(i,k,j) = 0.5 * (  &
                      fnm(k) * ( hat(i,k  ,j) + hat(i,k  ,j+1) ) +  &
                      fnp(k) * ( hat(i,k-1,j) + hat(i,k-1,j+1) ) )
    END DO
    END DO
    END DO

! Extrapolate to top and bottom of domain (to w levels).

    DO j = j_start, j_end
    DO i = i_start, i_end
      hatavg(i,1,j)   =  0.5 * (  &
                         cf1 * hat(i,1,j  ) +  &
                         cf2 * hat(i,2,j  ) +  &
                         cf3 * hat(i,3,j  ) +  &
                         cf1 * hat(i,1,j+1) +  &
                         cf2 * hat(i,2,j+1) +  &
                         cf3 * hat(i,3,j+1) )
      hatavg(i,kte,j) =  0.5 * (  &
                        cft1 * ( hat(i,ktes1,j) + hat(i,ktes1,j+1) ) +  &
                        cft2 * ( hat(i,ktes2,j) + hat(i,ktes2,j+1) ) )
    END DO
    END DO

    ! Comments 10-MAR-05
    ! Eqn 13b: D22=defor22= 2m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
    ! Below, D22 is set = 2*tmp1
    ! => tmp1 = m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
    ! tmpzy = averaged value of dpsi/dy (=zy)

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tmpzy       =  0.25 * (  &
                     zy(i,k  ,j) + zy(i,k  ,j+1) +  &
                     zy(i,k+1,j) + zy(i,k+1,j+1)  )
      tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * tmpzy * rdzw(i,k,j)
      ! tmp1 to here = partial dpsi/dy * partial dv^/dpsi:
    END DO
    END DO
    END DO

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tmp1(i,k,j) = mm(i,j) * (  &
                    rdy * ( hat(i,k,j+1) - hat(i,k,j) ) - tmp1(i,k,j) )
    END DO
    END DO
    END DO

! End calculation of dv/dy.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Calculate defor22 (2*dv/dy).
! Comments 10-MAR-05
! Eqn 13b: D22=defor22= 2 m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
!                     = 2*tmp1

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      defor22(i,k,j) = 2.0 * tmp1(i,k,j)
    END DO
    END DO
    END DO

! End calculation of defor22.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Calculate meridional divergence (dv/dy) and add it to the divergence
! array.

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      div(i,k,j) = div(i,k,j) + tmp1(i,k,j)
    END DO
    END DO
    END DO

! End calculation of meridional divergence.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Comments 10-MAR-05
! Eqn 13c: D33=defor33= 2 * partial dw/dz
! Below, D33 is set = 2*tmp1
! => tmp1 = partial dw/dz

! Calculate dw/dz.

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tmp1(i,k,j) = ( w(i,k+1,j) - w(i,k,j) ) * rdzw(i,k,j)
    END DO
    END DO
    END DO

! End calculation of dw/dz.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Calculate defor33 (2*dw/dz).

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      defor33(i,k,j) = 2.0 * tmp1(i,k,j)
    END DO
    END DO
    END DO

! End calculation of defor33.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Calculate vertical divergence (dw/dz) and add it to the divergence
! array.

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      div(i,k,j) = div(i,k,j) + tmp1(i,k,j)
    END DO
    END DO
    END DO

! End calculation of vertical divergence. 
!-----------------------------------------------------------------------

! Three-dimensional divergence is now finished and values are in array
! "div."  Also, the first three (defor11, defor22, defor33) of six
! deformation terms are now calculated at pressure points.
!=======================================================================

! Comments 10-MAR-2005
! Treat all differentials as 'div-style' [or 'curl-style'],
! i.e., du/dY becomes (in map coordinate space) mx*my * d(u/mx)/dy,
!       dv/dX becomes (in map coordinate space) mx*my * d(v/my)/dx,
! (see e.g. Haltiner and Williams p. 441)

!=======================================================================
! Calculate the final three deformations (defor12, defor13, defor23) at 
! vorticity points.

    i_start = its
    i_end   = ite
    j_start = jts
    j_end   = jte

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested) i_start = MAX( ids+1, its )
    IF ( config_flags%open_xe .OR. config_flags%specified .OR. & 
         config_flags%nested) i_end   = MIN( ide-1, ite )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested) j_start = MAX( jds+1, jts )
    IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
         config_flags%nested) j_end   = MIN( jde-1, jte )
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = ite


!-----------------------------------------------------------------------
! Calculate du/dy.

! First, calculate an average mapscale factor.

! Comments 10-MAR-05
! du/dy => need u map scale factor in x (which is defined at u points)
! averaged over j and j-1
! dv/dx => need v map scale factor in y (which is defined at v points)
! averaged over i and i-1

    DO j = j_start, j_end
    DO i = i_start, i_end
      mm(i,j) = 0.25 * ( msfux(i,j-1) + msfux(i,j) ) * ( msfvy(i-1,j) + msfvy(i,j) )
    END DO
    END DO

! Apply a coordinate transformation to zonal velocity, u.

    DO j =j_start-1, j_end
    DO k =kts, ktf
    DO i =i_start, i_end
      ! Fixes to set_physical_bc2/3d for polar boundary conditions 
      ! remove issues with loop over j
      hat(i,k,j) = u(i,k,j) / msfux(i,j)
    END DO
    END DO
    END DO

! Average in y and z.

    DO j=j_start,j_end
    DO k=kts+1,ktf
    DO i=i_start,i_end
      hatavg(i,k,j) = 0.5 * (  &
                      fnm(k) * ( hat(i,k  ,j-1) + hat(i,k  ,j) ) +  &
                      fnp(k) * ( hat(i,k-1,j-1) + hat(i,k-1,j) ) )
    END DO
    END DO
    END DO

! Extrapolate to top and bottom of domain (to w levels).

    DO j = j_start, j_end
    DO i = i_start, i_end
      hatavg(i,1,j)   =  0.5 * (  &
                         cf1 * hat(i,1,j-1) +  &
                         cf2 * hat(i,2,j-1) +  &
                         cf3 * hat(i,3,j-1) +  &
                         cf1 * hat(i,1,j  ) +  &
                         cf2 * hat(i,2,j  ) +  &
                         cf3 * hat(i,3,j  ) )
      hatavg(i,kte,j) =  0.5 * (  &
                        cft1 * ( hat(i,ktes1,j-1) + hat(i,ktes1,j) ) +  &
                        cft2 * ( hat(i,ktes2,j-1) + hat(i,ktes2,j) ) )
    END DO
    END DO

    ! tmpzy = averaged value of dpsi/dy (=zy) on vorticity grid
    ! tmp1  = partial dpsi/dy * partial du^/dpsi
    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tmpzy       = 0.25 * (  &
                    zy(i-1,k  ,j) + zy(i,k  ,j) +  &
                    zy(i-1,k+1,j) + zy(i,k+1,j) )
      tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) *  &
                    0.25 * tmpzy * ( rdzw(i,k,j) + rdzw(i-1,k,j) + &
                                     rdzw(i-1,k,j-1) + rdzw(i,k,j-1) )
    END DO
    END DO
    END DO

! End calculation of du/dy.
!---------------------------------------------------------------------- 

!-----------------------------------------------------------------------
! Add the first term to defor12 (du/dy+dv/dx) at vorticity points.

! Comments 10-MAR-05
! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
!                              partial dpsi/dx * partial dv^/dpsi +
!                              partial dpsi/dy * partial du^/dpsi)
! Here deal with m^2 * (partial du^/dY + partial dpsi/dy * partial du^/dpsi)
! Still need to add v^ terms: 
!   m^2 * (partial dv^/dX + partial dpsi/dx * partial dv^/dpsi)

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      defor12(i,k,j) = mm(i,j) * (  &
                       rdy * ( hat(i,k,j) - hat(i,k,j-1) ) - tmp1(i,k,j) )
    END DO
    END DO
    END DO

! End addition of the first term to defor12.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Calculate dv/dx.

! Apply a coordinate transformation to meridional velocity, v.

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start-1, i_end
       hat(i,k,j) = v(i,k,j) / msfvy(i,j)
    END DO
    END DO
    END DO

! Account for the slope in x of eta surfaces.

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      hatavg(i,k,j) = 0.5 * (  &
                      fnm(k) * ( hat(i-1,k  ,j) + hat(i,k  ,j) ) +  &
                      fnp(k) * ( hat(i-1,k-1,j) + hat(i,k-1,j) ) )
    END DO
    END DO
    END DO

! Extrapolate to top and bottom of domain (to w levels).

    DO j = j_start, j_end
    DO i = i_start, i_end
       hatavg(i,1,j)   =  0.5 * (  &
                          cf1 * hat(i-1,1,j) +  &
                          cf2 * hat(i-1,2,j) +  &
                          cf3 * hat(i-1,3,j) +  &
                          cf1 * hat(i  ,1,j) +  &
                          cf2 * hat(i  ,2,j) +  &
                          cf3 * hat(i  ,3,j) )
       hatavg(i,kte,j) =  0.5 * (  &
                         cft1 * ( hat(i,ktes1,j) + hat(i-1,ktes1,j) ) +  &
                         cft2 * ( hat(i,ktes2,j) + hat(i-1,ktes2,j) ) )
    END DO
    END DO

    ! Fixes to set_physical_bc2/3d have made any check for polar B.C.'s
    ! unnecessary in this place.  zx, rdzw, and hatavg are all defined
    ! in places they need to be and the values at the poles are replications
    ! of the values one grid point in, so the averaging over j and j-1 works
    ! to act as just using the value at j or j-1 (with out extra code).
    !
    ! tmpzx = averaged value of dpsi/dx (=zx) on vorticity grid
    ! tmp1  = partial dpsi/dx * partial dv^/dpsi
    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tmpzx       = 0.25 * (  &
                    zx(i,k  ,j-1) + zx(i,k  ,j) +  &
                    zx(i,k+1,j-1) + zx(i,k+1,j) )
      tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) *  &
                    0.25 * tmpzx * ( rdzw(i,k,j) + rdzw(i,k,j-1) + &
                                     rdzw(i-1,k,j-1) + rdzw(i-1,k,j) )
    END DO
    END DO
    END DO

! End calculation of dv/dx.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Add the second term to defor12 (du/dy+dv/dx) at vorticity points.

! Comments 10-MAR-05
! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
!                              partial dpsi/dx * partial dv^/dpsi +
!                              partial dpsi/dy * partial du^/dpsi)
! Here adding v^ terms:
!    m^2 * (partial dv^/dX + partial dpsi/dx * partial dv^/dpsi)

  IF ( config_flags%sfs_opt .GT. 0 ) THEN ! NBA--

!JDM____________________________________________________________________
!
!     s12 =  du/dy + dv/dx 
!         = (du/dy - dz/dy*du/dz) + (dv/dx - dz/dx*dv/dz)
!            ______defor12______             ___tmp1___
!
!     r12 =  du/dy - dv/dx 
!         = (du/dy - dz/dy*du/dz) - (dv/dx - dz/dx*dv/dz)
!            ______defor12______             ___tmp1___
!_______________________________________________________________________


    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end

      nba_rij(i,k,j,P_r12) = defor12(i,k,j) -  &    
                             mm(i,j) * (   &                            
                             rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) ) 

      defor12(i,k,j) = defor12(i,k,j) +  &
                       mm(i,j) * (  &
                       rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
    END DO
    END DO
    END DO

! End addition of the second term to defor12.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Update the boundary for defor12 (might need to change later).
 
    IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1 ) THEN
      DO j = jts, jte
      DO k = kts, kte
        defor12(ids,k,j) = defor12(ids+1,k,j)
        nba_rij(ids,k,j,P_r12) = nba_rij(ids+1,k,j,P_r12) 
      END DO
      END DO
    END IF
 
    IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
      DO k = kts, kte
      DO i = its, ite
        defor12(i,k,jds) = defor12(i,k,jds+1)
        nba_rij(i,k,jds,P_r12) = nba_rij(i,k,jds+1,P_r12) 
      END DO
      END DO
    END IF

    IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
      DO j = jts, jte
      DO k = kts, kte
        defor12(ide,k,j) = defor12(ide-1,k,j)
        nba_rij(ide,k,j,P_r12) = nba_rij(ide-1,k,j,P_r12) 
      END DO
      END DO
    END IF

    IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
      DO k = kts, kte
      DO i = its, ite
        defor12(i,k,jde) = defor12(i,k,jde-1)
        nba_rij(i,k,jde,P_r12) = nba_rij(i,k,jde-1,P_r12) 
      END DO
      END DO
    END IF

  ELSE ! NOT NBA--------------------------------------------------------

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      defor12(i,k,j) = defor12(i,k,j) +  &
                       mm(i,j) * (  &
                       rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
    END DO
    END DO
    END DO

! End addition of the second term to defor12.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Update the boundary for defor12 (might need to change later).
 
    IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1 ) THEN
      DO j = jts, jte
      DO k = kts, kte
        defor12(ids,k,j) = defor12(ids+1,k,j)
      END DO
      END DO
    END IF
 
    IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
      DO k = kts, kte
      DO i = its, ite
        defor12(i,k,jds) = defor12(i,k,jds+1)
      END DO
      END DO
    END IF

    IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
      DO j = jts, jte
      DO k = kts, kte
        defor12(ide,k,j) = defor12(ide-1,k,j)
      END DO
      END DO
    END IF

    IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
      DO k = kts, kte
      DO i = its, ite
        defor12(i,k,jde) = defor12(i,k,jde-1)
      END DO
      END DO
    END IF

  ENDIF ! NBA-----------------------------------------------------------

! End update of boundary for defor12.
!-----------------------------------------------------------------------

! Comments 10-MAR-05
! Further deformation terms not needed for 2-dimensional Smagorinsky diffusion,
! so those terms have not been dealt with yet.
! A "y" has simply been added to all map scale factors to allow the model to
! compile without errors.

!-----------------------------------------------------------------------
! Calculate dw/dx.

    i_start = its
    i_end   = MIN( ite, ide-1 )
    j_start = jts
    j_end   = MIN( jte, jde-1 )

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested) i_start = MAX( ids+1, its )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested) j_start = MAX( jds+1, jts )

    IF ( config_flags%periodic_x ) i_start = its
    IF ( config_flags%periodic_x ) i_end = MIN( ite, ide )
    IF ( config_flags%periodic_y ) j_end = MIN( jte, jde )

! Square the mapscale factor.

    DO j = jts, jte
    DO i = its, ite
      mm(i,j) = msfux(i,j) * msfuy(i,j)
    END DO
    END DO

! Apply a coordinate transformation to vertical velocity, w.  This is for both
! defor13 and defor23.

    DO j = j_start, j_end
    DO k = kts, kte
    DO i = i_start, i_end
      hat(i,k,j) = w(i,k,j) / msfty(i,j)
    END DO
    END DO
    END DO

    i = i_start-1
    DO j = j_start, MIN( jte, jde-1 )
    DO k = kts, kte
      hat(i,k,j) = w(i,k,j) / msfty(i,j)
    END DO
    END DO

    j = j_start-1
    DO k = kts, kte
    DO i = i_start, MIN( ite, ide-1 )
      hat(i,k,j) = w(i,k,j) / msfty(i,j)
    END DO
    END DO

! QUESTION: What is this for?

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      hatavg(i,k,j) = 0.25 * (  &
                      hat(i  ,k  ,j) +  &
                      hat(i  ,k+1,j) +  &
                      hat(i-1,k  ,j) +  &
                      hat(i-1,k+1,j) )
    END DO
    END DO
    END DO

! Calculate dw/dx.

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      tmp1(i,k,j) = ( hatavg(i,k,j) - hatavg(i,k-1,j) ) * zx(i,k,j) *  &
                    0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
    END DO
    END DO
    END DO

! End calculation of dw/dx.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Add the first term (dw/dx) to defor13 (dw/dx+du/dz) at vorticity
! points.

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      defor13(i,k,j) = mm(i,j) * (  &
                       rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
    END DO
    END DO
    END DO

    DO j = j_start, j_end
    DO i = i_start, i_end
      defor13(i,kts,j  ) = 0.0
      defor13(i,ktf+1,j) = 0.0
    END DO
    END DO

! End addition of the first term to defor13.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Calculate du/dz.

    IF ( config_flags%mix_full_fields ) THEN

      DO j = j_start, j_end
      DO k = kts+1, ktf
      DO i = i_start, i_end
        tmp1(i,k,j) = ( u(i,k,j) - u(i,k-1,j) ) *  &
                      0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
      END DO
      END DO
      END DO

    ELSE

      DO j = j_start, j_end
      DO k = kts+1, ktf
      DO i = i_start, i_end
        tmp1(i,k,j) = ( u(i,k,j) - u_base(k) - u(i,k-1,j) + u_base(k-1) ) *  &
                      0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
      END DO
      END DO
      END DO

    END IF

!-----------------------------------------------------------------------
! Add the second term (du/dz) to defor13 (dw/dx+du/dz) at vorticity
! points.


  IF ( config_flags%sfs_opt .GT. 0 ) THEN ! NBA--

!JDM____________________________________________________________________
!
!     s13 = du/dz +  dw/dx  
!         = du/dz + (dw/dx - dz/dx*dw/dz)
!         = tmp1  +  ______defor13______
!
!     r13 = du/dz -  dw/dx 
!         = du/dz - (dw/dx - dz/dx*dw/dz) 
!         = tmp1  -  ______defor13______  
!_______________________________________________________________________

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      nba_rij(i,k,j,P_r13) =  tmp1(i,k,j) - defor13(i,k,j)   
      defor13(i,k,j) = defor13(i,k,j) + tmp1(i,k,j)
    END DO
    END DO
    END DO

    DO j = j_start, j_end !change for different surface B. C. 
    DO i = i_start, i_end
      nba_rij(i,kts  ,j,P_r13) = 0.0
      nba_rij(i,ktf+1,j,P_r13) = 0.0
    END DO
    END DO

  ELSE ! NOT NBA--------------------------------------------------------

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      defor13(i,k,j) = defor13(i,k,j) + tmp1(i,k,j)
    END DO
    END DO
    END DO

  ENDIF ! NBA-----------------------------------------------------------

! End addition of the second term to defor13.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Calculate dw/dy.

    i_start = its
    i_end   = MIN( ite, ide-1 )
    j_start = jts
    j_end   = MIN( jte, jde-1 )

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested) i_start = MAX( ids+1, its )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested) j_start = MAX( jds+1, jts )
    IF ( config_flags%periodic_y ) j_end = MIN( jte, jde )
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )

! Square mapscale factor.

    DO j = jts, jte
    DO i = its, ite
      mm(i,j) = msfvx(i,j) * msfvy(i,j)
    END DO
    END DO

! Apply a coordinate transformation to vertical velocity, w.  Added by CW 7/19/07

    DO j = j_start, j_end
    DO k = kts, kte
    DO i = i_start, i_end
      hat(i,k,j) = w(i,k,j) / msftx(i,j)
    END DO
    END DO
    END DO

    i = i_start-1
    DO j = j_start, MIN( jte, jde-1 )
    DO k = kts, kte
      hat(i,k,j) = w(i,k,j) / msftx(i,j)
    END DO
    END DO

    j = j_start-1
    DO k = kts, kte
    DO i = i_start, MIN( ite, ide-1 )
      hat(i,k,j) = w(i,k,j) / msftx(i,j)
    END DO
    END DO

! QUESTION: What is this for?

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      hatavg(i,k,j) = 0.25 * (  &
                      hat(i,k  ,j  ) +  &
                      hat(i,k+1,j  ) +  &
                      hat(i,k  ,j-1) +  &
                      hat(i,k+1,j-1) )
    END DO
    END DO
    END DO

! Calculate dw/dy and store in tmp1.

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      tmp1(i,k,j) = ( hatavg(i,k,j) - hatavg(i,k-1,j) ) * zy(i,k,j) *  &
                    0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
    END DO
    END DO
    END DO

! End calculation of dw/dy.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Add the first term (dw/dy) to defor23 (dw/dy+dv/dz) at vorticity
! points.

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      defor23(i,k,j) = mm(i,j) * (  &
                       rdy * ( hat(i,k,j) - hat(i,k,j-1) ) - tmp1(i,k,j) )
    END DO
    END DO
    END DO

    DO j = j_start, j_end
    DO i = i_start, i_end
      defor23(i,kts,j  ) = 0.0
      defor23(i,ktf+1,j) = 0.0
    END DO
    END DO

! End addition of the first term to defor23.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Calculate dv/dz.

    IF ( config_flags%mix_full_fields ) THEN

      DO j = j_start, j_end
      DO k = kts+1, ktf
      DO i = i_start, i_end
        tmp1(i,k,j) = ( v(i,k,j) - v(i,k-1,j) ) *  &
                      0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
      END DO
      END DO
      END DO

    ELSE

      DO j = j_start, j_end
      DO k = kts+1, ktf
      DO i = i_start, i_end
        tmp1(i,k,j) = ( v(i,k,j) - v_base(k) - v(i,k-1,j) + v_base(k-1) ) *  &
                      0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
      END DO
      END DO
      END DO

    END IF

! End calculation of dv/dz.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Add the second term (dv/dz) to defor23 (dw/dy+dv/dz) at vorticity
! points.

  IF ( config_flags%sfs_opt .GT. 0 ) THEN ! NBA--

!JDM___________________________________________________________________
!
!     s23 = dv/dz +  dw/dy  
!         = dv/dz + (dw/dy - dz/dy*dw/dz)
!           tmp1  +  ______defor23______
!
!     r23 = dv/dz -  dw/dy 
!         = dv/dz - (dw/dy - dz/dy*dw/dz) 
!         = tmp1  -  ______defor23______  

! Add tmp1 to defor23.

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      nba_rij(i,k,j,P_r23) = tmp1(i,k,j) - defor23(i,k,j)  
      defor23(i,k,j) = defor23(i,k,j) + tmp1(i,k,j)
    END DO
    END DO
    END DO

    DO j = j_start, j_end
      DO i = i_start, i_end
        nba_rij(i,kts  ,j,P_r23) = 0.0
        nba_rij(i,ktf+1,j,P_r23) = 0.0
      END DO
    END DO

! End addition of the second term to defor23.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Update the boundary for defor13 and defor23 (might need to change
! later).

    IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1) THEN
      DO j = jts, jte
      DO k = kts, kte
        defor13(ids,k,j) = defor13(ids+1,k,j)
        defor23(ids,k,j) = defor23(ids+1,k,j)
        nba_rij(ids,k,j,P_r13) = nba_rij(ids+1,k,j,P_r13) 
        nba_rij(ids,k,j,P_r23) = nba_rij(ids+1,k,j,P_r23) 
      END DO
      END DO
    END IF

    IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
      DO k = kts, kte
      DO i = its, ite
        defor13(i,k,jds) = defor13(i,k,jds+1)
        defor23(i,k,jds) = defor23(i,k,jds+1)
        nba_rij(i,k,jds,P_r13) = nba_rij(i,k,jds+1,P_r13) 
        nba_rij(i,k,jds,P_r23) = nba_rij(i,k,jds+1,P_r23) 
      END DO
      END DO
    END IF

    IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
      DO j = jts, jte
      DO k = kts, kte
        defor13(ide,k,j) = defor13(ide-1,k,j)
        defor23(ide,k,j) = defor23(ide-1,k,j)
        nba_rij(ide,k,j,P_r13) = nba_rij(ide-1,k,j,P_r13) 
        nba_rij(ide,k,j,P_r23) = nba_rij(ide-1,k,j,P_r23) 
      END DO
      END DO
    END IF

    IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
      DO k = kts, kte
      DO i = its, ite
        defor13(i,k,jde) = defor13(i,k,jde-1)
        defor23(i,k,jde) = defor23(i,k,jde-1)
        nba_rij(i,k,jde,P_r13) = nba_rij(i,k,jde-1,P_r13) 
        nba_rij(i,k,jde,P_r23) = nba_rij(i,k,jde-1,P_r23) 
      END DO
      END DO
    END IF

  ELSE ! NOT NBA--------------------------------------------------------

! Add tmp1 to defor23.

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      defor23(i,k,j) = defor23(i,k,j) + tmp1(i,k,j)
    END DO
    END DO
    END DO

! End addition of the second term to defor23.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Update the boundary for defor13 and defor23 (might need to change
! later).

    IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1) THEN
      DO j = jts, jte
      DO k = kts, kte
        defor13(ids,k,j) = defor13(ids+1,k,j)
        defor23(ids,k,j) = defor23(ids+1,k,j)
      END DO
      END DO
    END IF

    IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
      DO k = kts, kte
      DO i = its, ite
        defor13(i,k,jds) = defor13(i,k,jds+1)
        defor23(i,k,jds) = defor23(i,k,jds+1)
      END DO
      END DO
    END IF

    IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
      DO j = jts, jte
      DO k = kts, kte
        defor13(ide,k,j) = defor13(ide-1,k,j)
        defor23(ide,k,j) = defor23(ide-1,k,j)
      END DO
      END DO
    END IF

    IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
      DO k = kts, kte
      DO i = its, ite
        defor13(i,k,jde) = defor13(i,k,jde-1)
        defor23(i,k,jde) = defor23(i,k,jde-1)
      END DO
      END DO
    END IF

  ENDIF ! NBA-----------------------------------------------------------

! End update of boundary for defor13 and defor23.
!-----------------------------------------------------------------------

! The second three (defor12, defor13, defor23) of six deformation terms
! are now calculated at vorticity points.
!=======================================================================

    END SUBROUTINE cal_deform_and_div

!=======================================================================
!=======================================================================

    SUBROUTINE calculate_km_kh( config_flags, dt,                        &
                                dampcoef, zdamp, damp_opt,               &
                                xkmh, xkmv, xkhh, xkhv,                  &
                                BN2, khdif, kvdif, div,                  &
                                defor11, defor22, defor33,               &
                                defor12, defor13, defor23,               &
                                tke, p8w, t8w, theta, t, p, moist,       &
                                dn, dnw, dx, dy, rdz, rdzw, isotropic,   &
                                n_moist, cf1, cf2, cf3, warm_rain,       &
                                mix_upper_bound,                         &
                                msftx, msfty,                            &
                                ids, ide, jds, jde, kds, kde,            &
                                ims, ime, jms, jme, kms, kme,            &
                                its, ite, jts, jte, kts, kte             )

! History:     Sep 2003  Changes by George Bryan and Jason Knievel, NCAR
!              Oct 2001  Converted to mass core by Bill Skamarock, NCAR
!              ...       ...

! Purpose:     This routine calculates exchange coefficients for the TKE
!              scheme.

! References:  Klemp and Wilhelmson (JAS 1978)
!              Deardorff (B-L Meteor 1980)
!              Chen and Dudhia (NCAR WRF physics report 2000)

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags   

    INTEGER, INTENT( IN )  &
    :: n_moist, damp_opt, isotropic,  & 
       ids, ide, jds, jde, kds, kde,  &
       ims, ime, jms, jme, kms, kme,  &
       its, ite, jts, jte, kts, kte 

    LOGICAL, INTENT( IN )  &
    :: warm_rain

    REAL, INTENT( IN )  &
    :: dx, dy, zdamp, dt, dampcoef, cf1, cf2, cf3, khdif, kvdif

    REAL, DIMENSION( kms:kme ), INTENT( IN )  &
    :: dnw, dn

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, n_moist ), INTENT( INOUT )  &
    :: moist

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &
    :: xkmv, xkmh, xkhv, xkhh, BN2  

    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),  INTENT( IN )  &
    :: defor11, defor22, defor33, defor12, defor13, defor23,      &
       div, rdz, rdzw, p8w, t8w, theta, t, p

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &
    :: tke

    REAL, INTENT( IN )  &
    :: mix_upper_bound

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    :: msftx, msfty

! Local variables.

    INTEGER  &
    :: i_start, i_end, j_start, j_end, ktf, i, j, k

! End declarations.
!-----------------------------------------------------------------------

    ktf     = MIN( kte, kde-1 )
    i_start = its
    i_end   = MIN( ite, ide-1 )
    j_start = jts
    j_end   = MIN( jte, jde-1 )

    CALL calculate_N2( config_flags, BN2, moist,           &
                       theta, t, p, p8w, t8w,              &
                       dnw, dn, rdz, rdzw,                 &
                       n_moist, cf1, cf2, cf3, warm_rain,  &
                       ids, ide, jds, jde, kds, kde,       &
                       ims, ime, jms, jme, kms, kme,       &
                       its, ite, jts, jte, kts, kte        )

! Select a scheme for calculating diffusion coefficients.

    km_coef: SELECT CASE( config_flags%km_opt )

      CASE (1)
            CALL isotropic_km( config_flags, xkmh, xkmv,                &
                               xkhh, xkhv, khdif, kvdif,                &
                               ids, ide, jds, jde, kds, kde,            &
                               ims, ime, jms, jme, kms, kme,            &
                               its, ite, jts, jte, kts, kte             )
      CASE (2)  
            CALL tke_km(       config_flags, xkmh, xkmv,                &
                               xkhh, xkhv, BN2, tke, p8w, t8w, theta,   &
                               rdz, rdzw, dx, dy, dt, isotropic,        &
                               mix_upper_bound, msftx, msfty,           &
                               ids, ide, jds, jde, kds, kde,            &
                               ims, ime, jms, jme, kms, kme,            &
                               its, ite, jts, jte, kts, kte             )
      CASE (3)  
            CALL smag_km(      config_flags, xkmh, xkmv,                &
                               xkhh, xkhv, BN2, div,                    &
                               defor11, defor22, defor33,               &
                               defor12, defor13, defor23,               &
                               rdzw, dx, dy, dt, isotropic,             &
                               mix_upper_bound, msftx, msfty,           &
                               ids, ide, jds, jde, kds, kde,            &
                               ims, ime, jms, jme, kms, kme,            &
                               its, ite, jts, jte, kts, kte             )
      CASE (4)  
            CALL smag2d_km(    config_flags, xkmh, xkmv,                &
                               xkhh, xkhv, defor11, defor22, defor12,   &
                               rdzw, dx, dy, msftx, msfty,              &
                               ids, ide, jds, jde, kds, kde,            &
                               ims, ime, jms, jme, kms, kme,            &
                               its, ite, jts, jte, kts, kte             )
      CASE DEFAULT
            CALL wrf_error_fatal( 'Please choose diffusion coefficient scheme' )

    END SELECT km_coef

    IF ( damp_opt .eq. 1 ) THEN
      CALL cal_dampkm( config_flags, xkmh, xkhh, xkmv, xkhv,    &
                       dx, dy, dt, dampcoef, rdz, rdzw, zdamp,  &
                       msftx, msfty,                            &
                       ids, ide, jds, jde, kds, kde,            &
                       ims, ime, jms, jme, kms, kme,            &
                       its, ite, jts, jte, kts, kte             )
    END IF

    END SUBROUTINE calculate_km_kh

!=======================================================================

SUBROUTINE cal_dampkm( config_flags,xkmh,xkhh,xkmv,xkhv,                       &
                       dx,dy,dt,dampcoef,                                      &
                       rdz, rdzw ,zdamp,                                       &
                       msftx, msfty,                                           &
                       ids,ide, jds,jde, kds,kde,                              &
                       ims,ime, jms,jme, kms,kme,                              &
                       its,ite, jts,jte, kts,kte                              )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type) , INTENT(IN   ) :: config_flags

   INTEGER ,          INTENT(IN   )           :: ids, ide, jds, jde, kds, kde, &
                                                 ims, ime, jms, jme, kms, kme, &
                                                 its, ite, jts, jte, kts, kte

   REAL    ,          INTENT(IN   )           :: zdamp,dx,dy,dt,dampcoef


   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT)    ::     xkmh , &
                                                                         xkhh , &
                                                                         xkmv , &
                                                                         xkhv 

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   )    ::     rdz,   &
                                                                         rdzw

   REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN   )             ::     msftx, &
                                                                         msfty
! LOCAL VARS

   INTEGER :: i_start, i_end, j_start, j_end, ktf, ktfm1, i, j, k
   REAL    :: kmmax,kmmvmax,degrad90,dz,tmp
   REAL    :: ds
   REAL ,     DIMENSION( its:ite )                                ::   deltaz
   REAL , DIMENSION( its:ite, kts:kte, jts:jte)                   ::   dampk,dampkv

! End declarations.
!-----------------------------------------------------------------------

   ktf = min(kte,kde-1)
   ktfm1 = ktf-1

   i_start = its
   i_end   = MIN(ite,ide-1)
   j_start = jts
   j_end   = MIN(jte,jde-1)

! keep upper damping diffusion away from relaxation zones at boundaries if used
   IF(config_flags%specified .OR. config_flags%nested)THEN
     i_start = MAX(i_start,ids+config_flags%spec_bdy_width-1)
     i_end   = MIN(i_end,ide-config_flags%spec_bdy_width)
     j_start = MAX(j_start,jds+config_flags%spec_bdy_width-1)
     j_end   = MIN(j_end,jde-config_flags%spec_bdy_width)
   ENDIF

   kmmax=dx*dx/dt
   degrad90=DEGRAD*90.
   DO j = j_start, j_end

      k=ktf
      DO i = i_start, i_end
         ! Unmodified dx used above may produce very large diffusivities
         ! when msftx is very large.  And the above formula ignores the fact
         ! that dy may now be different from dx as well.  Let's fix that by
         ! defining a "ds" as the minimum of the "real-space" (physical
         ! distance) values of dx and dy, and then using that smallest value
         ! to calculate a point-by-point kmmax
         ds = MIN(dx/msftx(i,j),dy/msfty(i,j))
         kmmax=ds*ds/dt

!         deltaz(i)=0.5*dnw(k)/zeta_z(i,j)
!         dz=dnw(k)/zeta_z(i,j)
         dz = 1./rdzw(i,k,j)
         deltaz(i) = 0.5*dz

         kmmvmax=dz*dz/dt
         tmp=min(deltaz(i)/zdamp,1.)
         dampk(i,k,j)=cos(degrad90*tmp)*cos(degrad90*tmp)*kmmax*dampcoef
         dampkv(i,k,j)=cos(degrad90*tmp)*cos(degrad90*tmp)*kmmvmax*dampcoef
! set upper limit on vertical K (based on horizontal K)
         dampkv(i,k,j)=min(dampkv(i,k,j),dampk(i,k,j))

      ENDDO

      DO k = ktfm1,kts,-1
      DO i = i_start, i_end
         ! Unmodified dx used above may produce very large diffusivities
         ! when msftx is very large.  And the above formula ignores the fact
         ! that dy may now be different from dx as well.  Let's fix that by
         ! defining a "ds" as the minimum of the "real-space" (physical
         ! distance) values of dx and dy, and then using that smallest value
         ! to calculate a point-by-point kmmax
         ds = MIN(dx/msftx(i,j),dy/msfty(i,j))
         kmmax=ds*ds/dt

!         deltaz(i)=deltaz(i)+dn(k)/zeta_z(i,j)
!         dz=dnw(k)/zeta_z(i,j)
         dz = 1./rdz(i,k,j)
         deltaz(i) = deltaz(i) + dz
         dz = 1./rdzw(i,k,j)

         kmmvmax=dz*dz/dt
         tmp=min(deltaz(i)/zdamp,1.)
         dampk(i,k,j)=cos(degrad90*tmp)*cos(degrad90*tmp)*kmmax*dampcoef
         dampkv(i,k,j)=cos(degrad90*tmp)*cos(degrad90*tmp)*kmmvmax*dampcoef
! set upper limit on vertical K (based on horizontal K)
         dampkv(i,k,j)=min(dampkv(i,k,j),dampk(i,k,j))
      ENDDO
      ENDDO

   ENDDO

   DO j = j_start, j_end
   DO k = kts,ktf
   DO i = i_start, i_end
      xkmh(i,k,j)=max(xkmh(i,k,j),dampk(i,k,j))
      xkhh(i,k,j)=max(xkhh(i,k,j),dampk(i,k,j))
      xkmv(i,k,j)=max(xkmv(i,k,j),dampkv(i,k,j))
      xkhv(i,k,j)=max(xkhv(i,k,j),dampkv(i,k,j))
   ENDDO
   ENDDO
   ENDDO

END SUBROUTINE cal_dampkm

!=======================================================================
!=======================================================================

    SUBROUTINE calculate_N2( config_flags, BN2, moist,           &
                             theta, t, p, p8w, t8w,              &
                             dnw, dn, rdz, rdzw,                 &
                             n_moist, cf1, cf2, cf3, warm_rain,  &
                             ids, ide, jds, jde, kds, kde,       &
                             ims, ime, jms, jme, kms, kme,       &
                             its, ite, jts, jte, kts, kte        )

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: n_moist,  &
       ids, ide, jds, jde, kds, kde, &
       ims, ime, jms, jme, kms, kme, &
       its, ite, jts, jte, kts, kte

    LOGICAL, INTENT( IN )  &
    :: warm_rain

    REAL, INTENT( IN )  &
    :: cf1, cf2, cf3

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &
    :: BN2

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
    :: rdz, rdzw, theta, t, p, p8w, t8w 

    REAL, DIMENSION( kms:kme ), INTENT( IN )  &
    :: dnw, dn

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), INTENT( INOUT )  &
    :: moist

! Local variables.

    INTEGER  &
    :: i, j, k, ktf, ispe, ktes1, ktes2,  &
       i_start, i_end, j_start, j_end

    REAL  &
    :: coefa, thetaep1, thetaem1, qc_cr, es, tc, qlpqi, qsw, qsi,  &
       tmpdz, xlvqv, thetaesfc, thetasfc, qvtop, qvsfc, thetatop, thetaetop

    REAL, DIMENSION( its:ite, jts:jte )  &
    :: tmp1sfc, tmp1top

    REAL, DIMENSION( its:ite, kts:kte, jts:jte )  &
    :: tmp1, qvs, qctmp

! End declarations.
!-----------------------------------------------------------------------

    qc_cr   = 0.00001  ! in Kg/Kg

    ktf     = MIN( kte, kde-1 )
    ktes1   = kte-1
    ktes2   = kte-2

    i_start = its
    i_end   = MIN( ite, ide-1 )
    j_start = jts
    j_end   = MIN( jte, jde-1 )

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested) i_start = MAX( ids+1, its )
    IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
         config_flags%nested) i_end   = MIN( ide-2, ite )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested) j_start = MAX( jds+1, jts )
    IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
         config_flags%nested) j_end   = MIN( jde-2 ,jte )
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
 
    IF ( P_QC .GT. PARAM_FIRST_SCALAR) THEN
      DO j = j_start, j_end
      DO k = kts, ktf
      DO i = i_start, i_end
        qctmp(i,k,j) = moist(i,k,j,P_QC)
      END DO
      END DO
      END DO
    ELSE
      DO j = j_start, j_end
      DO k = kts, ktf
      DO i = i_start, i_end
         qctmp(i,k,j) = 0.0
      END DO
      END DO
      END DO
    END IF
 
    DO j = jts, jte
    DO k = kts, kte
    DO i = its, ite
      tmp1(i,k,j) = 0.0
    END DO
    END DO
    END DO
 
    DO j = jts,jte
    DO i = its,ite
      tmp1sfc(i,j) = 0.0
      tmp1top(i,j) = 0.0
    END DO
    END DO
 
    DO ispe = PARAM_FIRST_SCALAR, n_moist
      IF ( ispe .EQ. P_QV .OR. ispe .EQ. P_QC .OR. ispe .EQ. P_QI) THEN
        DO j = j_start, j_end
        DO k = kts, ktf
        DO i = i_start, i_end
          tmp1(i,k,j) = tmp1(i,k,j) + moist(i,k,j,ispe)
        END DO
        END DO
        END DO
 
        DO j = j_start, j_end
        DO i = i_start, i_end
          tmp1sfc(i,j) = tmp1sfc(i,j) +  &
                         cf1 * moist(i,1,j,ispe) +  &
                         cf2 * moist(i,2,j,ispe) +  &
                         cf3 * moist(i,3,j,ispe)
          tmp1top(i,j) = tmp1top(i,j) +  &
                         moist(i,ktes1,j,ispe) + &
                         ( moist(i,ktes1,j,ispe) - moist(i,ktes2,j,ispe) ) *  &
                         0.5 * dnw(ktes1) / dn(ktes1)
        END DO
        END DO
      END IF
    END DO

! Calculate saturation mixing ratio.

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tc         = t(i,k,j) - SVPT0
      es         = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t(i,k,j) - SVP3 ) )
      qvs(i,k,j) = EP_2 * es / ( p(i,k,j) - es )
    END DO
    END DO
    END DO
 
    DO j = j_start, j_end
    DO k = kts+1, ktf-1
    DO i = i_start, i_end
      tmpdz = 1.0 / rdz(i,k,j) + 1.0 / rdz(i,k+1,j)
      IF ( moist(i,k,j,P_QV) .GE. qvs(i,k,j) .OR. qctmp(i,k,j) .GE. qc_cr) THEN
        xlvqv      = XLV * moist(i,k,j,P_QV)
        coefa      = ( 1.0 + xlvqv / R_d / t(i,k,j) ) / &
                     ( 1.0 + XLV * xlvqv / Cp / R_v / t(i,k,j) / t(i,k,j) ) /  &
                     theta(i,k,j)
        thetaep1   = theta(i,k+1,j) *  &
                     ( 1.0 + XLV * qvs(i,k+1,j) / Cp / t(i,k+1,j) )
        thetaem1   = theta(i,k-1,j) *  &
                     ( 1.0 + XLV * qvs(i,k-1,j) / Cp / t(i,k-1,j) )
        BN2(i,k,j) = g * ( coefa * ( thetaep1 - thetaem1 ) / tmpdz -  &
                     ( tmp1(i,k+1,j) - tmp1(i,k-1,j) ) / tmpdz )
      ELSE
        BN2(i,k,j) = g * ( (theta(i,k+1,j) - theta(i,k-1,j) ) /  &
                     theta(i,k,j) / tmpdz +  &
                     1.61 * ( moist(i,k+1,j,P_QV) - moist(i,k-1,j,P_QV) ) / &
                     tmpdz -   &
                     ( tmp1(i,k+1,j) - tmp1(i,k-1,j) ) / tmpdz )
      ENDIF
    END DO
    END DO
    END DO

    k = kts
    DO j = j_start, j_end
    DO i = i_start, i_end
      tmpdz     = 1.0 / rdz(i,k+1,j) + 0.5 / rdzw(i,k,j)
      thetasfc  = T8w(i,kts,j) / ( p8w(i,k,j) / p1000mb )**( R_d / Cp )
      IF ( moist(i,k,j,P_QV) .GE. qvs(i,k,j) .OR. qctmp(i,k,j) .GE. qc_cr) THEN
        qvsfc     = cf1 * qvs(i,1,j) +  &
                    cf2 * qvs(i,2,j) +  &
                    cf3 * qvs(i,3,j)
        xlvqv      = XLV * moist(i,k,j,P_QV)
        coefa      = ( 1.0 + xlvqv / R_d / t(i,k,j) ) /  &
                     ( 1.0 + XLV * xlvqv / Cp / R_v / t(i,k,j) / t(i,k,j) ) /  &
                     theta(i,k,j)
        thetaep1   = theta(i,k+1,j) *  &
                     ( 1.0 + XLV * qvs(i,k+1,j) / Cp / t(i,k+1,j) )
        thetaesfc  = thetasfc *  &
                     ( 1.0 + XLV * qvsfc / Cp / t8w(i,kts,j) )
        BN2(i,k,j) = g * ( coefa * ( thetaep1 - thetaesfc ) / tmpdz -  &
                     ( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz )
      ELSE
        qvsfc     = cf1 * moist(i,1,j,P_QV) +  &
                    cf2 * moist(i,2,j,P_QV) +  &
                    cf3 * moist(i,3,j,P_QV)
!        BN2(i,k,j) = g * ( ( theta(i,k+1,j) - thetasfc ) /  &
!                     theta(i,k,j) / tmpdz +  &
!                     1.61 * ( moist(i,k+1,j,P_QV) - qvsfc ) /  &
!                     tmpdz -  &
!                     ( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz  )
!...... MARTA: change in computation of BN2 at the surface, WCS 040331

        tmpdz= 1./rdzw(i,k,j) ! controlare come calcola rdzw
        BN2(i,k,j) = g * ( ( theta(i,k+1,j) - theta(i,k,j)) /  &
                     theta(i,k,j) / tmpdz +  &
                     1.61 * ( moist(i,k+1,j,P_QV) - qvsfc ) /  &
                     tmpdz -  &
                     ( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz  )
! end of MARTA/WCS change

      ENDIF
    END DO
    END DO
 

!...... MARTA: change in computation of BN2 at the top, WCS 040331
    DO j = j_start, j_end
    DO i = i_start, i_end
       BN2(i,ktf,j)=BN2(i,ktf-1,j)
    END DO
    END DO   
! end of MARTA/WCS change

    END SUBROUTINE calculate_N2

!=======================================================================
!=======================================================================

SUBROUTINE isotropic_km( config_flags,                                         &
                         xkmh,xkmv,xkhh,xkhv,khdif,kvdif,                      &
                         ids,ide, jds,jde, kds,kde,                            &
                         ims,ime, jms,jme, kms,kme,                            &
                         its,ite, jts,jte, kts,kte                            )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type) , INTENT(IN   ) :: config_flags

   INTEGER ,          INTENT(IN   )           :: ids, ide, jds, jde, kds, kde, &
                                                 ims, ime, jms, jme, kms, kme, &
                                                 its, ite, jts, jte, kts, kte

   REAL    ,          INTENT(IN   )           :: khdif,kvdif               

   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) ::     xkmh, &
                                                                         xkmv, &
                                                                         xkhh, &
                                                                         xkhv
! LOCAL VARS

   INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
   REAL    :: khdif3,kvdif3

! End declarations.
!-----------------------------------------------------------------------

   ktf = kte

   i_start = its
   i_end   = MIN(ite,ide-1)
   j_start = jts
   j_end   = MIN(jte,jde-1)

!   khdif3=khdif*3.
!   kvdif3=kvdif*3.
   khdif3=khdif/prandtl
   kvdif3=kvdif/prandtl

   DO j = j_start, j_end
   DO k = kts, ktf
   DO i = i_start, i_end
      xkmh(i,k,j)=khdif
      xkmv(i,k,j)=kvdif
      xkhh(i,k,j)=khdif3
      xkhv(i,k,j)=kvdif3
   ENDDO
   ENDDO
   ENDDO

END SUBROUTINE isotropic_km

!=======================================================================
!=======================================================================

SUBROUTINE smag_km( config_flags,xkmh,xkmv,xkhh,xkhv,BN2,                      &
                    div,defor11,defor22,defor33,defor12,                       &
                    defor13,defor23,                                           &
                    rdzw,dx,dy,dt,isotropic,                                   &
                    mix_upper_bound, msftx, msfty,                             &
                    ids,ide, jds,jde, kds,kde,                                 &
                    ims,ime, jms,jme, kms,kme,                                 &
                    its,ite, jts,jte, kts,kte                                  )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type) , INTENT(IN   ) :: config_flags

   INTEGER ,          INTENT(IN   )           :: ids, ide, jds, jde, kds, kde, &
                                                 ims, ime, jms, jme, kms, kme, &
                                                 its, ite, jts, jte, kts, kte

   INTEGER ,          INTENT(IN   )           :: isotropic
   REAL    ,          INTENT(IN   )           :: dx, dy, dt, mix_upper_bound


   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) ::      BN2, &
                                                                         rdzw

   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) ::     xkmh, &
                                                                         xkmv, &
                                                                         xkhh, &
                                                                         xkhv

   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ),  INTENT(IN   )      ::      &    
                                                                      defor11, &
                                                                      defor22, &
                                                                      defor33, &
                                                                      defor12, &
                                                                      defor13, &
                                                                      defor23, &
                                                                          div
   REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN   ) ::                msftx, &
                                                                        msfty
! LOCAL VARS

   INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
   REAL    :: deltas, tmp, pr, mlen_h, mlen_v, c_s

   REAL, DIMENSION( its:ite , kts:kte , jts:jte )                 ::     def2

! End declarations.
!-----------------------------------------------------------------------

   ktf = min(kte,kde-1)

   i_start = its
   i_end   = MIN(ite,ide-1)
   j_start = jts
   j_end   = MIN(jte,jde-1)

   IF ( config_flags%open_xs .or. config_flags%specified .or. &
        config_flags%nested) i_start = MAX(ids+1,its)
   IF ( config_flags%open_xe .or. config_flags%specified .or. &
        config_flags%nested) i_end   = MIN(ide-2,ite)
   IF ( config_flags%open_ys .or. config_flags%specified .or. &
        config_flags%nested) j_start = MAX(jds+1,jts)
   IF ( config_flags%open_ye .or. config_flags%specified .or. &
        config_flags%nested) j_end   = MIN(jde-2,jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )

   pr = prandtl
   c_s = config_flags%c_s

   do j=j_start,j_end
   do k=kts,ktf
   do i=i_start,i_end
      def2(i,k,j)=0.5*(defor11(i,k,j)*defor11(i,k,j) + &
                       defor22(i,k,j)*defor22(i,k,j) + &
                       defor33(i,k,j)*defor33(i,k,j))
   enddo
   enddo
   enddo

   do j=j_start,j_end
   do k=kts,ktf
   do i=i_start,i_end
      tmp=0.25*(defor12(i  ,k,j)+defor12(i  ,k,j+1)+ &
                defor12(i+1,k,j)+defor12(i+1,k,j+1))
      def2(i,k,j)=def2(i,k,j)+tmp*tmp
   enddo
   enddo
   enddo

   do j=j_start,j_end
   do k=kts,ktf
   do i=i_start,i_end
      tmp=0.25*(defor13(i  ,k+1,j)+defor13(i  ,k,j)+ &
                defor13(i+1,k+1,j)+defor13(i+1,k,j))
      def2(i,k,j)=def2(i,k,j)+tmp*tmp
   enddo
   enddo
   enddo

   do j=j_start,j_end
   do k=kts,ktf
   do i=i_start,i_end
      tmp=0.25*(defor23(i,k+1,j  )+defor23(i,k,j  )+ &
                defor23(i,k+1,j+1)+defor23(i,k,j+1))
      def2(i,k,j)=def2(i,k,j)+tmp*tmp
   enddo
   enddo
   enddo
!
   IF (isotropic .EQ. 0) THEN
      DO j = j_start, j_end
      DO k = kts, ktf
      DO i = i_start, i_end
         mlen_h=sqrt(dx/msftx(i,j) * dy/msfty(i,j))
         mlen_v= 1./rdzw(i,k,j)
         tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr)
         tmp=tmp**0.5
         xkmh(i,k,j)=max(c_s*c_s*mlen_h*mlen_h*tmp, 1.0E-6*mlen_h*mlen_h )
         xkmh(i,k,j)=min(xkmh(i,k,j), mix_upper_bound * mlen_h * mlen_h / dt )
         xkmv(i,k,j)=max(c_s*c_s*mlen_v*mlen_v*tmp, 1.0E-6*mlen_v*mlen_v )
         xkmv(i,k,j)=min(xkmv(i,k,j), mix_upper_bound * mlen_v * mlen_v / dt )
         xkhh(i,k,j)=xkmh(i,k,j)/pr
         xkhh(i,k,j)=min(xkhh(i,k,j), mix_upper_bound * mlen_h * mlen_h / dt )
         xkhv(i,k,j)=xkmv(i,k,j)/pr
         xkhv(i,k,j)=min(xkhv(i,k,j), mix_upper_bound * mlen_v * mlen_v / dt )
      ENDDO
      ENDDO
      ENDDO
   ELSE
      DO j = j_start, j_end
      DO k = kts, ktf
      DO i = i_start, i_end
         deltas=(dx/msftx(i,j) * dy/msfty(i,j)/rdzw(i,k,j))**0.33333333
         tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr)
         tmp=tmp**0.5
         xkmh(i,k,j)=max(c_s*c_s*deltas*deltas*tmp, 1.0E-6*deltas*deltas )
         xkmh(i,k,j)=min(xkmh(i,k,j), mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt )
         xkmv(i,k,j)=xkmh(i,k,j)
         xkmv(i,k,j)=min(xkmv(i,k,j), mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt )
         xkhh(i,k,j)=xkmh(i,k,j)/pr
         xkhh(i,k,j)=min(xkhh(i,k,j), mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt )
         xkhv(i,k,j)=xkmv(i,k,j)/pr
         xkhv(i,k,j)=min(xkhv(i,k,j), mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt )
      ENDDO
      ENDDO
      ENDDO
   ENDIF

END SUBROUTINE smag_km

!=======================================================================
!=======================================================================

SUBROUTINE smag2d_km( config_flags,xkmh,xkmv,xkhh,xkhv,                        &
                    defor11,defor22,defor12,                                   &
                    rdzw,dx,dy,msftx, msfty,                                   &
                    ids,ide, jds,jde, kds,kde,                                 &
                    ims,ime, jms,jme, kms,kme,                                 &
                    its,ite, jts,jte, kts,kte                                  )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type) , INTENT(IN   ) :: config_flags

   INTEGER ,          INTENT(IN   )           :: ids, ide, jds, jde, kds, kde, &
                                                 ims, ime, jms, jme, kms, kme, &
                                                 its, ite, jts, jte, kts, kte

   REAL    ,          INTENT(IN   )           :: dx, dy


   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) ::     rdzw

   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) ::     xkmh, &
                                                                         xkmv, &
                                                                         xkhh, &
                                                                         xkhv

   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ),  INTENT(IN   )      ::      &    
                                                                      defor11, &
                                                                      defor22, &
                                                                      defor12

   REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN   ) ::                msftx, &
                                                                        msfty
! LOCAL VARS

   INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
   REAL    :: deltas, tmp, pr, mlen_h, c_s

   REAL, DIMENSION( its:ite , kts:kte , jts:jte )                 ::     def2

! End declarations.
!-----------------------------------------------------------------------

   ktf = min(kte,kde-1)

   i_start = its
   i_end   = MIN(ite,ide-1)
   j_start = jts
   j_end   = MIN(jte,jde-1)

   IF ( config_flags%open_xs .or. config_flags%specified .or. &
        config_flags%nested) i_start = MAX(ids+1,its)
   IF ( config_flags%open_xe .or. config_flags%specified .or. &
        config_flags%nested) i_end   = MIN(ide-2,ite)
   IF ( config_flags%open_ys .or. config_flags%specified .or. &
        config_flags%nested) j_start = MAX(jds+1,jts)
   IF ( config_flags%open_ye .or. config_flags%specified .or. &
        config_flags%nested) j_end   = MIN(jde-2,jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )

   pr=prandtl
   c_s = config_flags%c_s

   do j=j_start,j_end
   do k=kts,ktf
   do i=i_start,i_end
      def2(i,k,j)=0.25*((defor11(i,k,j)-defor22(i,k,j))*(defor11(i,k,j)-defor22(i,k,j)))
      tmp=0.25*(defor12(i  ,k,j)+defor12(i  ,k,j+1)+ &
                defor12(i+1,k,j)+defor12(i+1,k,j+1))
      def2(i,k,j)=def2(i,k,j)+tmp*tmp
   enddo
   enddo
   enddo
!
      DO j = j_start, j_end
      DO k = kts, ktf
      DO i = i_start, i_end
         mlen_h=sqrt(dx/msftx(i,j) * dy/msfty(i,j))
         tmp=sqrt(def2(i,k,j))
!        xkmh(i,k,j)=max(c_s*c_s*mlen_h*mlen_h*tmp, 1.0E-6*mlen_h*mlen_h )
         xkmh(i,k,j)=c_s*c_s*mlen_h*mlen_h*tmp
         xkmh(i,k,j)=min(xkmh(i,k,j), 10.*mlen_h )
         xkmv(i,k,j)=0.
         xkhh(i,k,j)=xkmh(i,k,j)/pr
         xkhv(i,k,j)=0.
      ENDDO
      ENDDO
      ENDDO

END SUBROUTINE smag2d_km

!=======================================================================
!=======================================================================

    SUBROUTINE tke_km( config_flags, xkmh, xkmv, xkhh, xkhv,         &
                       bn2, tke, p8w, t8w, theta,                    &
                       rdz, rdzw, dx,dy, dt, isotropic,              &
                       mix_upper_bound, msftx, msfty,                &
                       ids, ide, jds, jde, kds, kde,                 &
                       ims, ime, jms, jme, kms, kme,                 &
                       its, ite, jts, jte, kts, kte                  )

! History:     Sep 2003   Changes by Jason Knievel and George Bryan, NCAR
!              Oct 2001   Converted to mass core by Bill Skamarock, NCAR
!              ...        ...

! Purpose:     This routine calculates the exchange coefficients for the
!              TKE turbulence parameterization.

! References:  Klemp and Wilhelmson (JAS 1978)
!              Chen and Dudhia (NCAR WRF physics report 2000)

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: ids, ide, jds, jde, kds, kde,  &
       ims, ime, jms, jme, kms, kme,  &
       its, ite, jts, jte, kts, kte

    INTEGER, INTENT( IN )  :: isotropic
    REAL, INTENT( IN )  &
    :: dx, dy, dt

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
    :: tke, p8w, t8w, theta, rdz, rdzw, bn2

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &
    :: xkmh, xkmv, xkhh, xkhv

    REAL, INTENT( IN )  &
    :: mix_upper_bound

   REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN   ) ::     msftx, &
                                                             msfty
! Local variables.

    REAL, DIMENSION( its:ite, kts:kte, jts:jte )  &
    :: l_scale

    REAL, DIMENSION( its:ite, kts:kte, jts:jte )  &
    :: dthrdn

    REAL  &
    :: deltas, tmp, mlen_s, mlen_h, mlen_v, tmpdz,  &
       thetasfc, thetatop, minkx, pr_inv, pr_inv_h, pr_inv_v, c_k

    INTEGER  &
    :: i_start, i_end, j_start, j_end, ktf, i, j, k

    REAL, PARAMETER :: tke_seed_value = 1.e-06
    REAL            :: tke_seed
    REAL, PARAMETER :: epsilon = 1.e-10

! End declarations.
!-----------------------------------------------------------------------

    ktf     = MIN( kte, kde-1 )
    i_start = its
    i_end   = MIN( ite, ide-1 )
    j_start = jts
    j_end   = MIN( jte, jde-1 )

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested) i_start = MAX( ids+1, its )
    IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
         config_flags%nested) i_end   = MIN( ide-2, ite )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested) j_start = MAX( jds+1, jts )
    IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
         config_flags%nested) j_end   = MIN( jde-2, jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )

! in the absence of surface drag or a surface heat flux, there
! is no way to generate tke without pre-existing tke.  Use
! tke_seed if the drag and flux are off.

    c_k = config_flags%c_k
    tke_seed = tke_seed_value
    if( (config_flags%tke_drag_coefficient .gt. epsilon) .or.  &
        (config_flags%tke_heat_flux .gt. epsilon)  ) tke_seed = 0.

    DO j = j_start, j_end
    DO k = kts+1, ktf-1
    DO i = i_start, i_end
      tmpdz         = 1.0 / ( rdz(i,k+1,j) + rdz(i,k,j) )
      dthrdn(i,k,j) = ( theta(i,k+1,j) - theta(i,k-1,j) ) / tmpdz
    END DO
    END DO
    END DO

    k = kts
    DO j = j_start, j_end
    DO i = i_start, i_end
      tmpdz         = 1.0 / ( rdzw(i,k+1,j) + rdzw(i,k,j) )
      thetasfc      = T8w(i,kts,j) / ( p8w(i,k,j) / p1000mb )**( R_d / Cp )
      dthrdn(i,k,j) = ( theta(i,k+1,j) - thetasfc ) / tmpdz
    END DO
    END DO

    k = ktf
    DO j = j_start, j_end
    DO i = i_start, i_end
      tmpdz         = 1.0 / rdz(i,k,j) + 0.5 / rdzw(i,k,j)
      thetatop      = T8w(i,kde,j) / ( p8w(i,kde,j) / p1000mb )**( R_d / Cp )
      dthrdn(i,k,j) = ( thetatop - theta(i,k-1,j) ) / tmpdz
    END DO
    END DO

    IF ( isotropic .EQ. 0 ) THEN
      DO j = j_start, j_end
      DO k = kts, ktf
      DO i = i_start, i_end
        mlen_h = SQRT( dx/msftx(i,j) * dy/msfty(i,j) )
        tmp    = SQRT( MAX( tke(i,k,j), tke_seed ) )
        deltas = 1.0 / rdzw(i,k,j)
        mlen_v = deltas
        IF ( dthrdn(i,k,j) .GT. 0.) THEN
          mlen_s = 0.76 * tmp / ( ABS( g / theta(i,k,j) * dthrdn(i,k,j) ) )**0.5
          mlen_v = MIN( mlen_v, mlen_s )
        END IF
        xkmh(i,k,j)  = MAX( c_k * tmp * mlen_h, 1.0E-6 * mlen_h * mlen_h )
        xkmh(i,k,j)  = MIN( xkmh(i,k,j), mix_upper_bound * mlen_h *mlen_h / dt )
        xkmv(i,k,j)  = MAX( c_k * tmp * mlen_v, 1.0E-6 * deltas * deltas )
        xkmv(i,k,j)  = MIN( xkmv(i,k,j), mix_upper_bound * deltas *deltas / dt )
        pr_inv_h     = 1./prandtl
        pr_inv_v     = 1.0 + 2.0 * mlen_v / deltas
        xkhh(i,k,j)  = xkmh(i,k,j) * pr_inv_h
        xkhv(i,k,j)  = xkmv(i,k,j) * pr_inv_v
      END DO
      END DO
      END DO
    ELSE
      CALL calc_l_scale( config_flags, tke, BN2, l_scale,      &
                         i_start, i_end, ktf, j_start, j_end,  &
                         dx, dy, rdzw, msftx, msfty,           &
                         ids, ide, jds, jde, kds, kde,         &
                         ims, ime, jms, jme, kms, kme,         &
                         its, ite, jts, jte, kts, kte          )
      DO j = j_start, j_end
      DO k = kts, ktf
      DO i = i_start, i_end
        tmp          = SQRT( MAX( tke(i,k,j), tke_seed ) )
        deltas       = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333
        xkmh(i,k,j)  = c_k * tmp * l_scale(i,k,j)
        xkmh(i,k,j)  = MIN( mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt,  xkmh(i,k,j) )
        xkmv(i,k,j)  = c_k * tmp * l_scale(i,k,j)
        xkmv(i,k,j)  = MIN( mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt ,  xkmv(i,k,j) )
        pr_inv       = 1.0 + 2.0 * l_scale(i,k,j) / deltas
        xkhh(i,k,j)  = MIN( mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt, xkmh(i,k,j) * pr_inv )
        xkhv(i,k,j)  = MIN( mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt, xkmv(i,k,j) * pr_inv )
      END DO
      END DO
      END DO
    END IF

    END SUBROUTINE tke_km

!=======================================================================
!=======================================================================

    SUBROUTINE calc_l_scale( config_flags, tke, BN2, l_scale,      &
                             i_start, i_end, ktf, j_start, j_end,  &
                             dx, dy, rdzw, msftx, msfty,           &
                             ids, ide, jds, jde, kds, kde,         &
                             ims, ime, jms, jme, kms, kme,         &
                             its, ite, jts, jte, kts, kte          )

! History:     Sep 2003   Written by Bryan and Knievel, NCAR

! Purpose:     This routine calculates the length scale, based on stability,
!              for TKE parameterization of subgrid-scale turbulence.

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: i_start, i_end, ktf, j_start, j_end,  &
       ids, ide, jds, jde, kds, kde,         &
       ims, ime, jms, jme, kms, kme,         &
       its, ite, jts, jte, kts, kte

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
    :: BN2, tke, rdzw

    REAL, INTENT( IN )  &
    :: dx, dy

    REAL, DIMENSION( its:ite, kts:kte, jts:jte ), INTENT( OUT )  &
    :: l_scale

    REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN   ) ::     msftx, &
                                                              msfty
! Local variables.

    INTEGER  &
    :: i, j, k

    REAL  &
    :: deltas, tmp

! End declarations.
!-----------------------------------------------------------------------

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      deltas         = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333
      l_scale(i,k,j) = deltas

      IF ( BN2(i,k,j) .gt. 1.0e-6 ) THEN
        tmp            = SQRT( MAX( tke(i,k,j), 1.0e-6 ) )
        l_scale(i,k,j) = 0.76 * tmp / SQRT( BN2(i,k,j) )
        l_scale(i,k,j) = MIN( l_scale(i,k,j), deltas)
        l_scale(i,k,j) = MAX( l_scale(i,k,j), 0.001 * deltas )
      END IF

    END DO
    END DO
    END DO

    END SUBROUTINE calc_l_scale

!=======================================================================
!=======================================================================

SUBROUTINE horizontal_diffusion_2 ( rt_tendf, ru_tendf, rv_tendf, rw_tendf,    &
                                    tke_tendf,                                 &
                                    moist_tendf, n_moist,                      &
                                    chem_tendf, n_chem,                        &
                                    scalar_tendf, n_scalar,                    &
                                    tracer_tendf, n_tracer,                    &
                                    thp, theta, mu, tke, config_flags,         &
                                    defor11, defor22, defor12,                 &
                                    defor13, defor23,                          &
                                    nba_mij, n_nba_mij,                        & !JDM
                                    div,                                       &
                                    moist, chem, scalar,tracer,                &
                                    msfux, msfuy, msfvx, msfvy,                &
                                    msftx, msfty, xkmh, xkhh,km_opt,           &
                                    rdx, rdy, rdz, rdzw, fnm, fnp,             &
                                    cf1, cf2, cf3, zx, zy, dn, dnw,            &
                                    ids, ide, jds, jde, kds, kde,              &
                                    ims, ime, jms, jme, kms, kme,              &
                                    its, ite, jts, jte, kts, kte               )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

   INTEGER ,        INTENT(IN   ) ::        ids, ide, jds, jde, kds, kde, &
                                            ims, ime, jms, jme, kms, kme, &
                                            its, ite, jts, jte, kts, kte

   INTEGER ,        INTENT(IN   ) ::        n_moist,n_chem,n_scalar,n_tracer,km_opt

   REAL ,           INTENT(IN   ) ::        cf1, cf2, cf3

   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnm
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnp
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: dnw
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::  dn

   REAL , DIMENSION( ims:ime, jms:jme) ,         INTENT(IN   ) ::   msfux, &
                                                                    msfuy, &
                                                                    msfvx, &
                                                                    msfvy, &
                                                                    msftx, &
                                                                    msfty, &
                                                                      mu

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::rt_tendf,&
                                                                 ru_tendf,&
                                                                 rv_tendf,&
                                                                 rw_tendf,&
                                                                tke_tendf

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist),                 &
          INTENT(INOUT) ::                                    moist_tendf

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem),                  &
          INTENT(INOUT) ::                                     chem_tendf

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar),                &
          INTENT(INOUT) ::                                   scalar_tendf

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_tracer),                &
          INTENT(INOUT) ::                                   tracer_tendf

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist),                 &
          INTENT(IN   ) ::                                          moist

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem),                  &
          INTENT(IN   ) ::                                          chem 

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) ,               &
          INTENT(IN   ) ::                                         scalar 

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_tracer) ,               &
          INTENT(IN   ) ::                                         tracer 

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::defor11, &
                                                                 defor22, &
                                                                 defor12, &
                                                                 defor13, &
                                                                 defor23, &
                                                                     div, &
                                                                    xkmh, &
                                                                    xkhh, &
                                                                      zx, &
                                                                      zy, &
                                                                   theta, &
                                                                     thp, &
                                                                     tke, &
                                                                     rdz, &
                                                                    rdzw


   REAL ,                                        INTENT(IN   ) ::    rdx, &
                                                                     rdy
   INTEGER, INTENT(  IN ) :: n_nba_mij !JDM 

   REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) &   !JDM 
   :: nba_mij

! LOCAL VARS
   
   INTEGER :: im, ic, is

!  REAL , DIMENSION(its-1:ite+1, kts:kte, jts-1:jte+1)       ::     xkhh

! End declarations.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Call diffusion subroutines.

    CALL horizontal_diffusion_u_2( ru_tendf, mu, config_flags,             &
                                   defor11, defor12, div,                  &
                                   nba_mij, n_nba_mij,                     & !JDM
                                   tke(ims,kms,jms),                       &
                                   msfux, msfuy, xkmh, rdx, rdy, fnm, fnp, &
                                   zx, zy, rdzw,                           &
                                   ids, ide, jds, jde, kds, kde,           &
                                   ims, ime, jms, jme, kms, kme,           &
                                   its, ite, jts, jte, kts, kte           )

    CALL horizontal_diffusion_v_2( rv_tendf, mu, config_flags,             &
                                   defor12, defor22, div,                  &
                                   nba_mij, n_nba_mij,                     & !JDM
                                   tke(ims,kms,jms),                       &
                                   msfvx, msfvy, xkmh, rdx, rdy, fnm, fnp, &
                                   zx, zy, rdzw,                           &
                                   ids, ide, jds, jde, kds, kde,           &
                                   ims, ime, jms, jme, kms, kme,           &
                                   its, ite, jts, jte, kts, kte           )

    CALL horizontal_diffusion_w_2( rw_tendf, mu, config_flags,             &
                                   defor13, defor23, div,                  &
                                   nba_mij, n_nba_mij,                     & !JDM
                                   tke(ims,kms,jms),                       &
                                   msftx, msfty, xkmh, rdx, rdy, fnm, fnp, &
                                   zx, zy, rdz,                            &
                                   ids, ide, jds, jde, kds, kde,           &
                                   ims, ime, jms, jme, kms, kme,           &
                                   its, ite, jts, jte, kts, kte           )

    CALL horizontal_diffusion_s  ( rt_tendf, mu, config_flags, thp,        &
                                   msftx, msfty, msfux, msfuy,             &
                                   msfvx, msfvy, xkhh, rdx, rdy,           &
                                   fnm, fnp, cf1, cf2, cf3,                &
                                   zx, zy, rdz, rdzw, dnw, dn,             &
                                   .false.,                                &
                                   ids, ide, jds, jde, kds, kde,           &
                                   ims, ime, jms, jme, kms, kme,           &
                                   its, ite, jts, jte, kts, kte           )

    IF (km_opt .eq. 2)                                                     &
    CALL horizontal_diffusion_s  ( tke_tendf(ims,kms,jms),                 &
                                   mu, config_flags,                       &
                                   tke(ims,kms,jms),                       &
                                   msftx, msfty, msfux, msfuy,             &
                                   msfvx, msfvy, xkhh, rdx, rdy,           &
                                   fnm, fnp, cf1, cf2, cf3,                &
                                   zx, zy, rdz, rdzw, dnw, dn,             &
                                   .true.,                                 &
                                   ids, ide, jds, jde, kds, kde,           &
                                   ims, ime, jms, jme, kms, kme,           &
                                   its, ite, jts, jte, kts, kte           )

    IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN 

      moist_loop: do im = PARAM_FIRST_SCALAR, n_moist

          CALL horizontal_diffusion_s( moist_tendf(ims,kms,jms,im),       &
                                       mu, config_flags,                  &
                                       moist(ims,kms,jms,im),             &
                                       msftx, msfty, msfux, msfuy,        &
                                       msfvx, msfvy, xkhh, rdx, rdy,      &
                                       fnm, fnp, cf1, cf2, cf3,           &
                                       zx, zy, rdz, rdzw, dnw, dn,        &
                                       .false.,                           &
                                       ids, ide, jds, jde, kds, kde,      &
                                       ims, ime, jms, jme, kms, kme,      &
                                       its, ite, jts, jte, kts, kte      )

      ENDDO moist_loop

    ENDIF

    IF (n_chem .ge. PARAM_FIRST_SCALAR) THEN 

      chem_loop: do ic = PARAM_FIRST_SCALAR, n_chem

        CALL horizontal_diffusion_s( chem_tendf(ims,kms,jms,ic),     &
                                     mu, config_flags,                 &
                                     chem(ims,kms,jms,ic),           &
                                     msftx, msfty, msfux, msfuy,       &
                                     msfvx, msfvy, xkhh, rdx, rdy,     &
                                     fnm, fnp, cf1, cf2, cf3,          &
                                     zx, zy, rdz, rdzw, dnw, dn,       &
                                     .false.,                          &
                                     ids, ide, jds, jde, kds, kde,     &
                                     ims, ime, jms, jme, kms, kme,     &
                                     its, ite, jts, jte, kts, kte     )

      ENDDO chem_loop

    ENDIF

    IF (n_tracer .ge. PARAM_FIRST_SCALAR) THEN 

      tracer_loop: do ic = PARAM_FIRST_SCALAR, n_tracer

        CALL horizontal_diffusion_s( tracer_tendf(ims,kms,jms,ic),     &
                                     mu, config_flags,                 &
                                     tracer(ims,kms,jms,ic),           &
                                     msftx, msfty, msfux, msfuy,       &
                                     msfvx, msfvy, xkhh, rdx, rdy,     &
                                     fnm, fnp, cf1, cf2, cf3,          &
                                     zx, zy, rdz, rdzw, dnw, dn,       &
                                     .false.,                          &
                                     ids, ide, jds, jde, kds, kde,     &
                                     ims, ime, jms, jme, kms, kme,     &
                                     its, ite, jts, jte, kts, kte     )

      ENDDO tracer_loop

    ENDIF
    IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN 

      scalar_loop: do is = PARAM_FIRST_SCALAR, n_scalar

        CALL horizontal_diffusion_s( scalar_tendf(ims,kms,jms,is),     &
                                     mu, config_flags,                 &
                                     scalar(ims,kms,jms,is),           &
                                     msftx, msfty, msfux, msfuy,       &
                                     msfvx, msfvy, xkhh, rdx, rdy,     &
                                     fnm, fnp, cf1, cf2, cf3,          &
                                     zx, zy, rdz, rdzw, dnw, dn,       &
                                     .false.,                          &
                                     ids, ide, jds, jde, kds, kde,     &
                                     ims, ime, jms, jme, kms, kme,     &
                                     its, ite, jts, jte, kts, kte     )

      ENDDO scalar_loop

    ENDIF

    END SUBROUTINE horizontal_diffusion_2

!=======================================================================
!=======================================================================

SUBROUTINE horizontal_diffusion_u_2( tendency, mu, config_flags,          &
                                     defor11, defor12, div,               &
                                     nba_mij, n_nba_mij,                  & !JDM
                                     tke,                                 &
                                     msfux, msfuy,                        &
                                     xkmh, rdx, rdy, fnm, fnp,            &
                                     zx, zy, rdzw,                        &
                                     ids, ide, jds, jde, kds, kde,        &
                                     ims, ime, jms, jme, kms, kme,        &
                                     its, ite, jts, jte, kts, kte        )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

   INTEGER ,        INTENT(IN   ) ::        ids, ide, jds, jde, kds, kde, &
                                            ims, ime, jms, jme, kms, kme, &
                                            its, ite, jts, jte, kts, kte

   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnm
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnp

   REAL , DIMENSION( ims:ime, jms:jme) ,         INTENT(IN   ) ::  msfux, &
                                                                   msfuy, &
                                                                      mu

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::   rdzw  
                                                                    
 
   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::defor11, &
                                                                 defor12, &
                                                                     div, &   
                                                                     tke, &   
                                                                    xkmh, &
                                                                      zx, &
                                                                      zy

   INTEGER, INTENT(  IN ) :: n_nba_mij !JDM

   REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij),  INTENT(INOUT) &  !JDM     
   :: nba_mij

   REAL ,                                        INTENT(IN   ) ::    rdx, &
                                                                     rdy
! Local data
   
   INTEGER :: i, j, k, ktf

   INTEGER :: i_start, i_end, j_start, j_end
   INTEGER :: is_ext,ie_ext,js_ext,je_ext  

   REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1)    :: titau1avg, &
                                                              titau2avg, &
                                                                 titau1, & 
                                                                 titau2, & 
                                                                 xkxavg, & 
                                                                  rravg
! new
!                                                                 zxavg, & 
!                                                                 zyavg
   REAL :: mrdx, mrdy, rcoup

   REAL :: tmpzy, tmpzeta_z

   REAL :: term1, term2, term3

! End declarations.
!-----------------------------------------------------------------------

   ktf=MIN(kte,kde-1)
 
!-----------------------------------------------------------------------
! u :   p (.), u(|), w(-)
!       
!       p  u  p  u                                  u     u
!
! p  |  .  |  .  |  .  |   k+1                |  .  |  .  |  .  |   k+1
!           
! w     - 13  -     -      k+1                     13               k+1 
!
! p  |  11 O 11  |  .  |   k                  |  12 O 12  |  .  |   k      
!
! w     - 13  -     -      k                       13               k  
!
! p  |  .  |  .  |  .  |   k-1                |  .  |  .  |  .  |   k-1
!
!      i-1 i  i i+1                          j-1 j  j j+1 j+1         
!

   i_start = its
   i_end   = ite
   j_start = jts
   j_end   = MIN(jte,jde-1)

   IF ( config_flags%open_xs .or. config_flags%specified .or. &
        config_flags%nested) i_start = MAX(ids+1,its)
   IF ( config_flags%open_xe .or. config_flags%specified .or. &
        config_flags%nested) i_end   = MIN(ide-1,ite)
   IF ( config_flags%open_ys .or. config_flags%specified .or. &
        config_flags%nested) j_start = MAX(jds+1,jts)
   IF ( config_flags%open_ye .or. config_flags%specified .or. &
        config_flags%nested) j_end   = MIN(jde-2,jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = ite

! titau1 = titau11 
   is_ext=1
   ie_ext=0
   js_ext=0
   je_ext=0
   CALL cal_titau_11_22_33( config_flags, titau1,            &
                            mu, tke, xkmh, defor11,          &
                            nba_mij(ims,kms,jms,P_m11),      & !JDM
                            is_ext, ie_ext, js_ext, je_ext,  &
                            ids, ide, jds, jde, kds, kde,    &
                            ims, ime, jms, jme, kms, kme,    &
                            its, ite, jts, jte, kts, kte     )

! titau2 = titau12
   is_ext=0
   ie_ext=0
   js_ext=0
   je_ext=1
   CALL cal_titau_12_21( config_flags, titau2,            &
                         mu, xkmh, defor12,               &
                         nba_mij(ims,kms,jms,P_m12),      & !JDM
                         is_ext, ie_ext, js_ext, je_ext,  &
                         ids, ide, jds, jde, kds, kde,    &
                         ims, ime, jms, jme, kms, kme,    &
                         its, ite, jts, jte, kts, kte     )

! titau1avg = titau11avg
! titau2avg = titau12avg 

   DO j = j_start, j_end
   DO k = kts+1,ktf
   DO i = i_start, i_end
      titau1avg(i,k,j)=0.5*(fnm(k)*(titau1(i-1,k  ,j)+titau1(i,k  ,j))+ &
                            fnp(k)*(titau1(i-1,k-1,j)+titau1(i,k-1,j)))
      titau2avg(i,k,j)=0.5*(fnm(k)*(titau2(i,k  ,j+1)+titau2(i,k  ,j))+ &
                            fnp(k)*(titau2(i,k-1,j+1)+titau2(i,k-1,j)))
      tmpzy = 0.25*( zy(i-1,k,j  )+zy(i,k,j  )+ &
                     zy(i-1,k,j+1)+zy(i,k,j+1)  )
!      tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j))
!      titau1avg(i,k,j)=titau1avg(i,k,j)*zx(i,k,j)*tmpzeta_z
!      titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy    *tmpzeta_z

      titau1avg(i,k,j)=titau1avg(i,k,j)*zx(i,k,j)
      titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy    

   ENDDO
   ENDDO
   ENDDO
!
   DO j = j_start, j_end
   DO i = i_start, i_end
      titau1avg(i,kts,j)=0.
      titau1avg(i,ktf+1,j)=0.
      titau2avg(i,kts,j)=0.
      titau2avg(i,ktf+1,j)=0.
   ENDDO
   ENDDO
!
   DO j = j_start, j_end
   DO k = kts,ktf
   DO i = i_start, i_end

      mrdx=msfux(i,j)*rdx
      mrdy=msfuy(i,j)*rdy
      tendency(i,k,j)=tendency(i,k,j)-                                    &
           (mrdx*(titau1(i,k,j  )-titau1(i-1,k,j))+                       &
            mrdy*(titau2(i,k,j+1)-titau2(i,k,j  ))-                       &
            msfuy(i,j)*rdzw(i,k,j)*((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ &
                                   (titau2avg(i,k+1,j)-titau2avg(i,k,j))  &
                                  )                                      )
   ENDDO
   ENDDO
   ENDDO

END SUBROUTINE horizontal_diffusion_u_2

!=======================================================================
!=======================================================================

SUBROUTINE horizontal_diffusion_v_2( tendency, mu, config_flags,          &
                                     defor12, defor22, div,               &
                                     nba_mij, n_nba_mij,                  & !JDM
                                     tke,                                 &
                                     msfvx, msfvy,                        &
                                     xkmh, rdx, rdy, fnm, fnp,            &
                                     zx, zy, rdzw,                        &
                                     ids, ide, jds, jde, kds, kde,        &
                                     ims, ime, jms, jme, kms, kme,        &
                                     its, ite, jts, jte, kts, kte        )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

   INTEGER ,        INTENT(IN   ) ::        ids, ide, jds, jde, kds, kde, &
                                            ims, ime, jms, jme, kms, kme, &
                                            its, ite, jts, jte, kts, kte

   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnm
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnp

   REAL , DIMENSION( ims:ime, jms:jme) ,         INTENT(IN   ) ::  msfvx, &
                                                                   msfvy, &
                                                                      mu

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::defor12, &
                                                                 defor22, &
                                                                     div, &
                                                                     tke, &
                                                                    xkmh, &
                                                                      zx, &
                                                                      zy, &
                                                                    rdzw

   INTEGER,  INTENT(  IN ) :: n_nba_mij !JDM

   REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij),  INTENT(INOUT) &  !JDM     
   :: nba_mij

   REAL ,                                        INTENT(IN   ) ::    rdx, &
                                                                     rdy

! Local data

   INTEGER :: i, j, k, ktf

   INTEGER :: i_start, i_end, j_start, j_end
   INTEGER :: is_ext,ie_ext,js_ext,je_ext  

   REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1)    :: titau1avg, &
                                                              titau2avg, &
                                                                 titau1, &
                                                                 titau2, &
                                                                 xkxavg, &
                                                                  rravg
! new
!                                                                 zxavg, &
!                                                                 zyavg

   REAL :: mrdx, mrdy, rcoup

   REAL :: tmpzx, tmpzeta_z

! End declarations.
!-----------------------------------------------------------------------

   ktf=MIN(kte,kde-1)
 
!-----------------------------------------------------------------------
! v :   p (.), v(+), w(-)
!       
!       p  v  p  v                                  v     v
!
! p  +  .  +  .  +  .  +   k+1                +  .  +  .  +  .  +   k+1
!           
! w     - 23  -     -      k+1                     23               k+1 
!
! p  +  22 O 22  +  .  +   k                  +  21 O 21  +  .  +   k      
!
! w     - 23  -     -      k                       23               k  
!
! p  +  .  +  .  +  .  +   k-1                +  .  +  .  +  .  +   k-1
!
!      j-1 j  j j+1                          i-1 i  i i+1 i+1         
!

   i_start = its
   i_end   = MIN(ite,ide-1)
   j_start = jts
   j_end   = jte

   IF ( config_flags%open_xs .or. config_flags%specified .or. &
        config_flags%nested) i_start = MAX(ids+1,its)
   IF ( config_flags%open_xe .or. config_flags%specified .or. &
        config_flags%nested) i_end   = MIN(ide-2,ite)
   IF ( config_flags%open_ys .or. config_flags%specified .or. &
        config_flags%nested) j_start = MAX(jds+1,jts)
   IF ( config_flags%open_ye .or. config_flags%specified .or. &
        config_flags%nested) j_end   = MIN(jde-1,jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)

! titau1 = titau21
   is_ext=0
   ie_ext=1
   js_ext=0
   je_ext=0
   CALL cal_titau_12_21( config_flags, titau1,          &
                         mu, xkmh, defor12,             &
                         nba_mij(ims,kms,jms,P_m12),    & !JDM
                         is_ext,ie_ext,js_ext,je_ext,   &
                         ids, ide, jds, jde, kds, kde,  &
                         ims, ime, jms, jme, kms, kme,  &
                         its, ite, jts, jte, kts, kte   )

! titau2 = titau22
   is_ext=0
   ie_ext=0
   js_ext=1
   je_ext=0
   CALL cal_titau_11_22_33( config_flags, titau2,           &
                            mu, tke, xkmh, defor22,         &
                            nba_mij(ims,kms,jms,P_m22),     & !JDM
                            is_ext, ie_ext, js_ext, je_ext, &
                            ids, ide, jds, jde, kds, kde,   &
                            ims, ime, jms, jme, kms, kme,   &
                            its, ite, jts, jte, kts, kte    )

   DO j = j_start, j_end
   DO k = kts+1,ktf
   DO i = i_start, i_end
      titau1avg(i,k,j)=0.5*(fnm(k)*(titau1(i+1,k  ,j)+titau1(i,k  ,j))+ &
                            fnp(k)*(titau1(i+1,k-1,j)+titau1(i,k-1,j)))
      titau2avg(i,k,j)=0.5*(fnm(k)*(titau2(i,k  ,j-1)+titau2(i,k  ,j))+ &
                            fnp(k)*(titau2(i,k-1,j-1)+titau2(i,k-1,j)))

      tmpzx = 0.25*( zx(i,k,j  )+zx(i+1,k,j  )+ &
                     zx(i,k,j-1)+zx(i+1,k,j-1)  )


      titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx
      titau2avg(i,k,j)=titau2avg(i,k,j)*zy(i,k,j)


   ENDDO
   ENDDO
   ENDDO

   DO j = j_start, j_end
   DO i = i_start, i_end
      titau1avg(i,kts,j)=0.
      titau1avg(i,ktf+1,j)=0.
      titau2avg(i,kts,j)=0.
      titau2avg(i,ktf+1,j)=0.
   ENDDO
   ENDDO
!
   DO j = j_start, j_end
   DO k = kts,ktf
   DO i = i_start, i_end
       
      mrdx=msfvx(i,j)*rdx
      mrdy=msfvy(i,j)*rdy
      tendency(i,k,j)=tendency(i,k,j)-                                    &
           (mrdy*(titau2(i  ,k,j)-titau2(i,k,j-1))+                       &
            mrdx*(titau1(i+1,k,j)-titau1(i,k,j  ))-                       &
           msfvy(i,j)*rdzw(i,k,j)*((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ &
                                   (titau2avg(i,k+1,j)-titau2avg(i,k,j))  &
                                )			                  &
           )

   ENDDO
   ENDDO
   ENDDO

END SUBROUTINE horizontal_diffusion_v_2

!=======================================================================
!=======================================================================

SUBROUTINE horizontal_diffusion_w_2( tendency, mu, config_flags,          &
                                     defor13, defor23, div,               &
                                     nba_mij, n_nba_mij,                  & !JDM
                                     tke,                                 &
                                     msftx, msfty,                        &
                                     xkmh, rdx, rdy, fnm, fnp,            &
                                     zx, zy, rdz,                         &
                                     ids, ide, jds, jde, kds, kde,        &
                                     ims, ime, jms, jme, kms, kme,        &
                                     its, ite, jts, jte, kts, kte        )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

   INTEGER ,        INTENT(IN   ) ::        ids, ide, jds, jde, kds, kde, &
                                            ims, ime, jms, jme, kms, kme, &
                                            its, ite, jts, jte, kts, kte

   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnm
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnp

   REAL , DIMENSION( ims:ime, jms:jme) ,         INTENT(IN   ) ::  msftx, &
                                                                   msfty, &
                                                                      mu

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::defor13, &
                                                                 defor23, &
                                                                     div, &
                                                                     tke, &
                                                                    xkmh, &
                                                                      zx, &
                                                                      zy, &
                                                                     rdz

   INTEGER,  INTENT(  IN ) :: n_nba_mij !JDM

   REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij),  INTENT(INOUT) &   !JDM    
   :: nba_mij

   REAL ,                                        INTENT(IN   ) ::    rdx, &
                                                                     rdy

! Local data

   INTEGER :: i, j, k, ktf

   INTEGER :: i_start, i_end, j_start, j_end
   INTEGER :: is_ext,ie_ext,js_ext,je_ext  

   REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1)    :: titau1avg, &
                                                              titau2avg, &
                                                                 titau1, &
                                                                 titau2, &
                                                                 xkxavg, &
                                                                  rravg
! new
!                                                                 zxavg, &
!                                                                 zyavg

   REAL :: mrdx, mrdy, rcoup

   REAL :: tmpzx, tmpzy, tmpzeta_z

! End declarations.
!-----------------------------------------------------------------------

   ktf=MIN(kte,kde-1)
 
!-----------------------------------------------------------------------
! w :   p (.), u(|), v(+), w(-)
!       
!       p  u  p  u                               p  v  p  v 
!
! w     -     -     -      k+1             w     -     -     -      k+1 
!
! p     .  | 33  |  .      k               p     .  + 33  +  .      k      
!
! w     -  31 O 31  -      k               w     -  32 O 32  -      k   
!
! p     .  | 33  |  .      k-1             p     .  | 33  |  .      k-1 
!
! w     -     -     -      k-1             w     -     -     -      k-1 
!
!      i-1 i  i i+1                             j-1 j  j j+1         
!
   i_start = its
   i_end   = MIN(ite,ide-1)
   j_start = jts
   j_end   = MIN(jte,jde-1)

   IF ( config_flags%open_xs .or. config_flags%specified .or. &
        config_flags%nested) i_start = MAX(ids+1,its)
   IF ( config_flags%open_xe .or. config_flags%specified .or. &
        config_flags%nested) i_end   = MIN(ide-2,ite)
   IF ( config_flags%open_ys .or. config_flags%specified .or. &
        config_flags%nested) j_start = MAX(jds+1,jts)
   IF ( config_flags%open_ye .or. config_flags%specified .or. &
        config_flags%nested) j_end   = MIN(jde-2,jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)

! titau1 = titau31
   is_ext=0
   ie_ext=1
   js_ext=0
   je_ext=0
   CALL cal_titau_13_31( config_flags, titau1, defor13,   &
                         nba_mij(ims,kms,jms,P_m13),      & !JDM
                         mu, xkmh, fnm, fnp,              &
                         is_ext, ie_ext, js_ext, je_ext,  &
                         ids, ide, jds, jde, kds, kde,    &
                         ims, ime, jms, jme, kms, kme,    &
                         its, ite, jts, jte, kts, kte     )

! titau2 = titau32
   is_ext=0
   ie_ext=0
   js_ext=0
   je_ext=1
   CALL cal_titau_23_32( config_flags, titau2, defor23,   &
                         nba_mij(ims,kms,jms,P_m23),      & !JDM
                         mu, xkmh, fnm, fnp,              &
                         is_ext, ie_ext, js_ext, je_ext,  &
                         ids, ide, jds, jde, kds, kde,    &
                         ims, ime, jms, jme, kms, kme,    &
                         its, ite, jts, jte, kts, kte     )

! titau1avg = titau31avg * zx * zeta_z = titau13avg * zx * zeta_z
! titau2avg = titau32avg * zy * zeta_z = titau23avg * zy * zeta_z

   DO j = j_start, j_end
   DO k = kts,ktf
   DO i = i_start, i_end
      titau1avg(i,k,j)=0.25*(titau1(i+1,k+1,j)+titau1(i,k+1,j)+ &
                             titau1(i+1,k  ,j)+titau1(i,k  ,j))
      titau2avg(i,k,j)=0.25*(titau2(i,k+1,j+1)+titau2(i,k+1,j)+ &
                             titau2(i,k  ,j+1)+titau2(i,k  ,j))
! new
      tmpzx  =0.25*( zx(i,k  ,j)+zx(i+1,k  ,j)+ &
                     zx(i,k+1,j)+zx(i+1,k+1,j)  )
      tmpzy  =0.25*( zy(i,k  ,j)+zy(i,k  ,j+1)+ &
                     zy(i,k+1,j)+zy(i,k+1,j+1)  )

      titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx
      titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy
!      titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx*zeta_z(i,j)
!      titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy*zeta_z(i,j)
   ENDDO
   ENDDO
   ENDDO

   DO j = j_start, j_end
   DO i = i_start, i_end
      titau1avg(i,ktf+1,j)=0.
      titau2avg(i,ktf+1,j)=0.
   ENDDO
   ENDDO

   DO j = j_start, j_end
   DO k = kts+1,ktf
   DO i = i_start, i_end

      mrdx=msftx(i,j)*rdx
      mrdy=msfty(i,j)*rdy

      tendency(i,k,j)=tendency(i,k,j)-                                 &
           (mrdx*(titau1(i+1,k,j)-titau1(i,k,j))+                      &
            mrdy*(titau2(i,k,j+1)-titau2(i,k,j))-                      &
           msfty(i,j)*rdz(i,k,j)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
                                  titau2avg(i,k,j)-titau2avg(i,k-1,j)  &
                               )				       &
           )
!            msft(i,j)/dn(k)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
!                                titau2avg(i,k,j)-titau2avg(i,k-1,j)  &
!                               )				     &
!           )
   ENDDO
   ENDDO
   ENDDO

END SUBROUTINE horizontal_diffusion_w_2

!=======================================================================
!=======================================================================

SUBROUTINE horizontal_diffusion_s (tendency, mu, config_flags, var,       &
                                   msftx, msfty, msfux, msfuy,            &
                                   msfvx, msfvy, xkhh, rdx, rdy,          &
                                   fnm, fnp, cf1, cf2, cf3,               &
                                   zx, zy, rdz, rdzw, dnw, dn,            &
                                   doing_tke,                             &
                                   ids, ide, jds, jde, kds, kde,          &
                                   ims, ime, jms, jme, kms, kme,          &
                                   its, ite, jts, jte, kts, kte           )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

   INTEGER ,        INTENT(IN   ) ::        ids, ide, jds, jde, kds, kde, &
                                            ims, ime, jms, jme, kms, kme, &
                                            its, ite, jts, jte, kts, kte

   LOGICAL,         INTENT(IN   ) ::        doing_tke

   REAL , INTENT(IN   )           ::        cf1, cf2, cf3

   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnm
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnp
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::     dn
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    dnw

   REAL , DIMENSION( ims:ime, jms:jme) ,         INTENT(IN   ) ::   msfux
   REAL , DIMENSION( ims:ime, jms:jme) ,         INTENT(IN   ) ::   msfuy
   REAL , DIMENSION( ims:ime, jms:jme) ,         INTENT(IN   ) ::   msfvx
   REAL , DIMENSION( ims:ime, jms:jme) ,         INTENT(IN   ) ::   msfvy
   REAL , DIMENSION( ims:ime, jms:jme) ,         INTENT(IN   ) ::   msftx
   REAL , DIMENSION( ims:ime, jms:jme) ,         INTENT(IN   ) ::   msfty

   REAL , DIMENSION( ims:ime, jms:jme) ,         INTENT(IN   ) ::   mu

!  REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1),                 &
!         INTENT(IN   ) ::                                         xkhh

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::         &
                                                                    xkhh, &
                                                                     rdz, &
                                                                     rdzw

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::    var, &
                                                                      zx, &
                                                                      zy

   REAL ,                                        INTENT(IN   ) ::    rdx, &
                                                                     rdy

! Local data

   INTEGER :: i, j, k, ktf

   INTEGER :: i_start, i_end, j_start, j_end

   REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1)    ::     H1avg, &
                                                                  H2avg, &
                                                                     H1, &
                                                                     H2, &
                                                                 xkxavg
! new
!                                                                 zxavg, &
!                                                                 zyavg

   REAL , DIMENSION( its:ite, kts:kte, jts:jte)            ::  tmptendf

   REAL    :: mrdx, mrdy, rcoup
   REAL    :: tmpzx, tmpzy, tmpzeta_z, rdzu, rdzv
   INTEGER :: ktes1,ktes2

! End declarations.
!-----------------------------------------------------------------------

   ktf=MIN(kte,kde-1)
 
!-----------------------------------------------------------------------
! scalars:   t (.), u(|), v(+), w(-)
!       
!       t  u  t  u                               t  v  t  v 
!
! w     -     3     -      k+1             w     -     3     -      k+1 
!
! t     .  1  O  1  .      k               t     .  2  O  2  .      k      
!
! w     -     3     -      k               w     -     3     -      k   
!
! t     .  |  .  |  .      k-1             t     .  +  .  +  .      k-1 
!
! w     -     -     -      k-1             w     -     -     -      k-1 
!
! t    i-1 i  i i+1                             j-1 j  j j+1         
!

   ktes1=kte-1
   ktes2=kte-2

   i_start = its
   i_end   = MIN(ite,ide-1)
   j_start = jts
   j_end   = MIN(jte,jde-1)

   IF ( config_flags%open_xs .or. config_flags%specified .or. &
        config_flags%nested) i_start = MAX(ids+1,its)
   IF ( config_flags%open_xe .or. config_flags%specified .or. &
        config_flags%nested) i_end   = MIN(ide-2,ite)
   IF ( config_flags%open_ys .or. config_flags%specified .or. &
        config_flags%nested) j_start = MAX(jds+1,jts)
   IF ( config_flags%open_ye .or. config_flags%specified .or. &
        config_flags%nested) j_end   = MIN(jde-2,jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)

! diffusion of the TKE needs mutiple 2

   IF ( doing_tke ) THEN
      DO j = j_start, j_end
      DO k = kts,ktf
      DO i = i_start, i_end
         tmptendf(i,k,j)=tendency(i,k,j)
      ENDDO
      ENDDO
      ENDDO
   ENDIF

! H1 = partial var over partial x

   DO j = j_start, j_end
   DO k = kts, ktf
   DO i = i_start, i_end + 1
! new
!     zxavg(i,k,j) =0.5*( zx(i-1,k,j)+ zx(i,k,j))
      xkxavg(i,k,j)=0.5*(xkhh(i-1,k,j)+xkhh(i,k,j))
   ENDDO
   ENDDO
   ENDDO

   DO j = j_start, j_end
   DO k = kts+1, ktf
   DO i = i_start, i_end + 1
      H1avg(i,k,j)=0.5*(fnm(k)*(var(i-1,k  ,j)+var(i,k  ,j))+  &
                        fnp(k)*(var(i-1,k-1,j)+var(i,k-1,j)))
   ENDDO
   ENDDO
   ENDDO

   DO j = j_start, j_end
   DO i = i_start, i_end + 1
      H1avg(i,kts  ,j)=0.5*(cf1*var(i  ,1,j)+cf2*var(i  ,2,j)+ &
                            cf3*var(i  ,3,j)+cf1*var(i-1,1,j)+  &
                            cf2*var(i-1,2,j)+cf3*var(i-1,3,j))
      H1avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
                            var(i,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)+ &
                            var(i-1,ktes1,j)+(var(i-1,ktes1,j)- &
                            var(i-1,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1))
   ENDDO
   ENDDO

   DO j = j_start, j_end
   DO k = kts, ktf
   DO i = i_start, i_end + 1
! new
      tmpzx = 0.5*( zx(i,k,j)+ zx(i,k+1,j))
      rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
      H1(i,k,j)=-msfuy(i,j)*xkxavg(i,k,j)*(                      &
                 rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx*         &
                     (H1avg(i,k+1,j)-H1avg(i,k,j))*rdzu )

!      tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j))
!      H1(i,k,j)=-msfuy(i,j)*xkxavg(i,k,j)*(                         &
!                 rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx*tmpzeta_z*  &
!                     (H1avg(i,k+1,j)-H1avg(i,k,j))/dnw(k))
   ENDDO
   ENDDO
   ENDDO

! H2 = partial var over partial y

   DO j = j_start, j_end + 1
   DO k = kts, ktf
   DO i = i_start, i_end
! new
!     zyavg(i,k,j) =0.5*( zy(i,k,j-1)+ zy(i,k,j))
      xkxavg(i,k,j)=0.5*(xkhh(i,k,j-1)+xkhh(i,k,j))
   ENDDO
   ENDDO
   ENDDO

   DO j = j_start, j_end + 1
   DO k = kts+1,   ktf
   DO i = i_start, i_end
! new
      H2avg(i,k,j)=0.5*(fnm(k)*(var(i,k  ,j-1)+var(i,k  ,j))+  &
                        fnp(k)*(var(i,k-1,j-1)+var(i,k-1,j)))
   ENDDO
   ENDDO
   ENDDO

   DO j = j_start, j_end + 1
   DO i = i_start, i_end
      H2avg(i,kts  ,j)=0.5*(cf1*var(i,1,j  )+cf2*var(i  ,2,j)+ &
                            cf3*var(i,3,j  )+cf1*var(i,1,j-1)+  &
                            cf2*var(i,2,j-1)+cf3*var(i,3,j-1))
      H2avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
                            var(i,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)+ &
                            var(i,ktes1,j-1)+(var(i,ktes1,j-1)- &
                            var(i,ktes2,j-1))*0.5*dnw(ktes1)/dn(ktes1))
   ENDDO
   ENDDO

   DO j = j_start, j_end + 1
   DO k = kts, ktf
   DO i = i_start, i_end
! new
      tmpzy = 0.5*( zy(i,k,j)+ zy(i,k+1,j))
      rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
      H2(i,k,j)=-msfvy(i,j)*xkxavg(i,k,j)*(                       &
                 rdy*(var(i,k,j)-var(i,k,j-1)) - tmpzy*          &
                     (H2avg(i ,k+1,j)-H2avg(i,k,j))*rdzv)

!      tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i,j-1))
!      H2(i,k,j)=-msfvy(i,j)*xkxavg(i,k,j)*(                         &
!                 rdy*(var(i,k,j)-var(i,k,j-1)) - tmpzy*tmpzeta_z*  &
!                     (H2avg(i ,k+1,j)-H2avg(i,k,j))/dnw(k))
   ENDDO
   ENDDO
   ENDDO

   DO j = j_start, j_end
   DO k = kts+1, ktf
   DO i = i_start, i_end
      H1avg(i,k,j)=0.5*(fnm(k)*(H1(i+1,k  ,j)+H1(i,k  ,j))+  &
                        fnp(k)*(H1(i+1,k-1,j)+H1(i,k-1,j)))
      H2avg(i,k,j)=0.5*(fnm(k)*(H2(i,k  ,j+1)+H2(i,k  ,j))+  &
                        fnp(k)*(H2(i,k-1,j+1)+H2(i,k-1,j)))
! new
!     zxavg(i,k,j)=fnm(k)*zx(i,k,j)+fnp(k)*zx(i,k-1,j)
!     zyavg(i,k,j)=fnm(k)*zy(i,k,j)+fnp(k)*zy(i,k-1,j)

! H1avg(i,k,j)=zx*H1avg*zeta_z
! H2avg(i,k,j)=zy*H2avg*zeta_z

      tmpzx = 0.5*( zx(i,k,j)+ zx(i+1,k,j  ))
      tmpzy = 0.5*( zy(i,k,j)+ zy(i  ,k,j+1))

      H1avg(i,k,j)=H1avg(i,k,j)*tmpzx
      H2avg(i,k,j)=H2avg(i,k,j)*tmpzy

!      H1avg(i,k,j)=H1avg(i,k,j)*tmpzx*zeta_z(i,j)
!      H2avg(i,k,j)=H2avg(i,k,j)*tmpzy*zeta_z(i,j)
   ENDDO
   ENDDO
   ENDDO
 
   DO j = j_start, j_end
   DO i = i_start, i_end
      H1avg(i,kts  ,j)=0.
      H1avg(i,ktf+1,j)=0.
      H2avg(i,kts  ,j)=0.
      H2avg(i,ktf+1,j)=0.
   ENDDO
   ENDDO

   DO j = j_start, j_end
   DO k = kts,ktf
   DO i = i_start, i_end

      mrdx=msftx(i,j)*rdx
      mrdy=msfty(i,j)*rdy

      tendency(i,k,j)=tendency(i,k,j)-                      &
           (mrdx*0.5*((mu(i+1,j)+mu(i,j))*H1(i+1,k,j)-      &
                      (mu(i-1,j)+mu(i,j))*H1(i  ,k,j))+     &
            mrdy*0.5*((mu(i,j+1)+mu(i,j))*H2(i,k,j+1)-      &
                      (mu(i,j-1)+mu(i,j))*H2(i,k,j  ))-     &
           msfty(i,j)*mu(i,j)*(H1avg(i,k+1,j)-H1avg(i,k,j)+ &
                       H2avg(i,k+1,j)-H2avg(i,k,j)          &
                                )*rdzw(i,k,j)               &
                                                          )

   ENDDO
   ENDDO
   ENDDO
           
   IF ( doing_tke ) THEN
      DO j = j_start, j_end
      DO k = kts,ktf
      DO i = i_start, i_end
          tendency(i,k,j)=tmptendf(i,k,j)+2.* &
                          (tendency(i,k,j)-tmptendf(i,k,j))
      ENDDO
      ENDDO
      ENDDO
   ENDIF

END SUBROUTINE horizontal_diffusion_s

!=======================================================================
!=======================================================================

SUBROUTINE vertical_diffusion_2   ( ru_tendf, rv_tendf, rw_tendf, rt_tendf,   &
                                    tke_tendf, moist_tendf, n_moist,          &
                                    chem_tendf, n_chem,                       &
                                    scalar_tendf, n_scalar,                   &
                                    tracer_tendf, n_tracer,                   &
                                    u_2, v_2,                                 &
                                    thp,u_base,v_base,t_base,qv_base,mu,tke,  &
                                    config_flags,defor13,defor23,defor33,     &
                                    nba_mij, n_nba_mij,                       & !JDM
                                    div,                                      &
                                    moist,chem,scalar,tracer,xkmv,xkhv,km_opt,&
                                    fnm, fnp, dn, dnw, rdz, rdzw,             &
                                    hfx, qfx, ust, rho,                       &
                                    ids, ide, jds, jde, kds, kde,             &
                                    ims, ime, jms, jme, kms, kme,             &
                                    its, ite, jts, jte, kts, kte              )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

   INTEGER ,        INTENT(IN   ) ::        ids, ide, jds, jde, kds, kde, &
                                            ims, ime, jms, jme, kms, kme, &
                                            its, ite, jts, jte, kts, kte

   INTEGER ,        INTENT(IN   ) ::        n_moist,n_chem,n_scalar,n_tracer,km_opt

   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fnm
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fnp
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: dnw
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::  dn
   REAL , DIMENSION( ims:ime , jms:jme ) ,  INTENT(IN   )      ::  mu

   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: qv_base
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::  u_base
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::  v_base
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::  t_base

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::ru_tendf,&
                                                                 rv_tendf,&
                                                                 rw_tendf,&
                                                                tke_tendf,&
                                                                rt_tendf  

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist),                 &
          INTENT(INOUT) ::                                    moist_tendf

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem),                  &
          INTENT(INOUT) ::                                     chem_tendf

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) ,               &
          INTENT(INOUT) ::                                   scalar_tendf
   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_tracer) ,               &
          INTENT(INOUT) ::                                   tracer_tendf

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist),                 &
          INTENT(INOUT) ::                                          moist

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem),                  &
          INTENT(INOUT) ::                                           chem

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) ,               &
          INTENT(IN   ) ::                                         scalar
   REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_tracer) ,               &
          INTENT(IN   ) ::                                         tracer

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::defor13, &
                                                                 defor23, &
                                                                 defor33, &
                                                                     div, &
                                                                    xkmv, &
                                                                    xkhv, &
                                                                     tke, &
                                                                     rdz, &
                                                                     u_2, &
                                                                     v_2, &
                                                                    rdzw

   INTEGER, INTENT(  IN ) :: n_nba_mij !JDM

   REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) &  !JDM     
   :: nba_mij

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN   )   :: rho  
   REAL , DIMENSION( ims:ime, jms:jme), INTENT(INOUT)            :: hfx,  &
                                                                    qfx
   REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN   )            :: ust
   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::    thp

! LOCAL VAR

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme)  ::    var_mix

   INTEGER :: im, i,j,k
   INTEGER :: i_start, i_end, j_start, j_end

!  REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: xkhv

!***************************************************************************
!***************************************************************************
!MODIFICA VARIABILI PER I FLUSSI
!
    REAL :: V0_u,V0_v,tao_xz,tao_yz,ustar,cd0
    REAL :: xsfc,psi1,vk2,zrough,lnz
    REAL :: heat_flux, moist_flux, heat_flux0
    REAL :: cpm
!
!FINE MODIFICA VARIABILI PER I FLUSSI
!***************************************************************************
!

! End declarations.
!-----------------------------------------------------------------------

   i_start = its
   i_end   = MIN(ite,ide-1)
   j_start = jts
   j_end   = MIN(jte,jde-1)
!
!-----------------------------------------------------------------------

      CALL vertical_diffusion_u_2( ru_tendf, config_flags, mu,    &
                                   defor13, xkmv,                 &
                                   nba_mij, n_nba_mij,            & !JDM
                                   dnw, rdzw, fnm, fnp,           &
                                   ids, ide, jds, jde, kds, kde,  &
                                   ims, ime, jms, jme, kms, kme,  &
                                   its, ite, jts, jte, kts, kte  )


      CALL vertical_diffusion_v_2( rv_tendf, config_flags, mu,    &
                                   defor23, xkmv,                 &
                                   nba_mij, n_nba_mij,            & !JDM
                                   dnw, rdzw, fnm, fnp,           &
                                   ids, ide, jds, jde, kds, kde,  &
                                   ims, ime, jms, jme, kms, kme,  &
                                   its, ite, jts, jte, kts, kte  )

      CALL vertical_diffusion_w_2( rw_tendf, config_flags, mu,    &
                                   defor33, tke(ims,kms,jms),     &
                                   nba_mij, n_nba_mij,            & !JDM
                                   div, xkmv,                     &
                                   dn, rdz,                       &  
                                   ids, ide, jds, jde, kds, kde,  &
                                   ims, ime, jms, jme, kms, kme,  &
                                   its, ite, jts, jte, kts, kte  )

!*****************************************
!*****************************************
!  MODIFICA al flusso di momento alla parete
!
  vflux: SELECT CASE( config_flags%isfflx )
  CASE (0) ! Assume cd a constant, specified in namelist
    cd0 = config_flags%tke_drag_coefficient  ! constant drag coefficient
                                             ! set in namelist.input
!
!calcolo del modulo della velocita
    DO j = j_start, j_end
    DO i = i_start, ite
       V0_u=0.
       tao_xz=0.
       V0_u=    sqrt((u_2(i,kts,j)**2) +         &
                        (((v_2(i  ,kts,j  )+          &
                           v_2(i  ,kts,j+1)+          &
                           v_2(i-1,kts,j  )+          &
                           v_2(i-1,kts,j+1))/4)**2))+epsilon
       tao_xz=cd0*V0_u*u_2(i,kts,j)
       ru_tendf(i,kts,j)=ru_tendf(i,kts,j)            &
                         -0.25*(mu(i,j)+mu(i-1,j))*tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j))
    ENDDO
    ENDDO
 
!
    DO j = j_start, jte
    DO i = i_start, i_end
       V0_v=0.
       tao_yz=0.
       V0_v=    sqrt((v_2(i,kts,j)**2) +         &
                        (((u_2(i  ,kts,j  )+          &
                           u_2(i  ,kts,j-1)+          &
                           u_2(i+1,kts,j  )+          &
                           u_2(i+1,kts,j-1))/4)**2))+epsilon
       tao_yz=cd0*V0_v*v_2(i,kts,j)
       rv_tendf(i,kts,j)=rv_tendf(i,kts,j)            &
                         -0.25*(mu(i,j)+mu(i,j-1))*tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1))
    ENDDO
    ENDDO

  CASE (1,2) ! ustar computed from surface routine
    DO j = j_start, j_end
    DO i = i_start, ite
       V0_u=0.
       tao_xz=0.
       V0_u=    sqrt((u_2(i,kts,j)**2) +         &
                        (((v_2(i  ,kts,j  )+          &
                           v_2(i  ,kts,j+1)+          &
                           v_2(i-1,kts,j  )+          &
                           v_2(i-1,kts,j+1))/4)**2))+epsilon
       ustar=0.5*(ust(i,j)+ust(i-1,j))
       tao_xz=ustar*ustar*u_2(i,kts,j)/V0_u
       ru_tendf(i,kts,j)=ru_tendf(i,kts,j)            &
                         -0.25*(mu(i,j)+mu(i-1,j))*tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j))
    ENDDO
    ENDDO
 
    DO j = j_start, jte
    DO i = i_start, i_end
       V0_v=0.
       tao_yz=0.
       V0_v=    sqrt((v_2(i,kts,j)**2) +         &
                        (((u_2(i  ,kts,j  )+          &
                           u_2(i  ,kts,j-1)+          &
                           u_2(i+1,kts,j  )+          &
                           u_2(i+1,kts,j-1))/4)**2))+epsilon
       ustar=0.5*(ust(i,j)+ust(i,j-1))
       tao_yz=ustar*ustar*v_2(i,kts,j)/V0_v
       rv_tendf(i,kts,j)=rv_tendf(i,kts,j)            &
                         -0.25*(mu(i,j)+mu(i,j-1))*tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1))
    ENDDO
    ENDDO

  CASE DEFAULT
    CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
  END SELECT vflux

!
!  FINE MODIFICA al flusso di momento alla parete
!*****************************************
!*****************************************

   IF ( config_flags%mix_full_fields ) THEN

     DO j=jts,min(jte,jde-1)
     DO k=kts,kte-1
     DO i=its,min(ite,ide-1)
       var_mix(i,k,j) = thp(i,k,j)
     ENDDO
     ENDDO
     ENDDO

   ELSE

     DO j=jts,min(jte,jde-1)
     DO k=kts,kte-1
     DO i=its,min(ite,ide-1)
       var_mix(i,k,j) = thp(i,k,j) - t_base(k)
     ENDDO
     ENDDO
     ENDDO

   END IF

   CALL vertical_diffusion_s( rt_tendf, config_flags, var_mix, mu, xkhv, &
                              dn, dnw, rdz, rdzw, fnm, fnp,          &
                              .false.,                               &
                              ids, ide, jds, jde, kds, kde,          &
                              ims, ime, jms, jme, kms, kme,          &
                              its, ite, jts, jte, kts, kte          )


!*****************************************
!*****************************************
!MODIFICA al flusso di calore
!
!
  hflux: SELECT CASE( config_flags%isfflx )
  CASE (0,2) ! with fixed surface heat flux given in the namelist
    heat_flux = config_flags%tke_heat_flux  ! constant heat flux value
                                            ! set in namelist.input
    DO j = j_start, j_end
    DO i = i_start, i_end
       cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV)) 
       hfx(i,j)=heat_flux*cp*rho(i,1,j)         ! provided for output only
       rt_tendf(i,kts,j)=rt_tendf(i,kts,j)  &
            +mu(i,j)*heat_flux*rdzw(i,kts,j)
    ENDDO
    ENDDO

  CASE (1) ! use surface heat flux computed from surface routine
    DO j = j_start, j_end
    DO i = i_start, i_end

       cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV))
       heat_flux = hfx(i,j)/cpm/rho(i,1,j)
       rt_tendf(i,kts,j)=rt_tendf(i,kts,j)  &
            +mu(i,j)*heat_flux*rdzw(i,kts,j)

    ENDDO
    ENDDO

  CASE DEFAULT
    CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
  END SELECT hflux

!
! FINE MODIFICA al flusso di calore
!*****************************************
!*****************************************

   If (km_opt .eq. 2) then
   CALL vertical_diffusion_s( tke_tendf(ims,kms,jms),               &
                              config_flags, tke(ims,kms,jms),       &
                              mu, xkhv,                             &
                              dn, dnw, rdz, rdzw, fnm, fnp,         &
                              .true.,                               &
                              ids, ide, jds, jde, kds, kde,         &
                              ims, ime, jms, jme, kms, kme,         &
                              its, ite, jts, jte, kts, kte         )
   endif
 
   IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN 

     moist_loop: do im = PARAM_FIRST_SCALAR, n_moist

       IF ( (.not. config_flags%mix_full_fields) .and. (im == P_QV) ) THEN

         DO j=jts,min(jte,jde-1)
         DO k=kts,kte-1
         DO i=its,min(ite,ide-1)
          var_mix(i,k,j) = moist(i,k,j,im) - qv_base(k)
         ENDDO
         ENDDO
         ENDDO

       ELSE

         DO j=jts,min(jte,jde-1)
         DO k=kts,kte-1
         DO i=its,min(ite,ide-1)
          var_mix(i,k,j) = moist(i,k,j,im)
         ENDDO
         ENDDO
         ENDDO

       END IF


          CALL vertical_diffusion_s( moist_tendf(ims,kms,jms,im),         &
                                     config_flags, var_mix,               &
                                     mu, xkhv,                            &
                                     dn, dnw, rdz, rdzw, fnm, fnp,        &
                                     .false.,                             &
                                     ids, ide, jds, jde, kds, kde,        &
                                     ims, ime, jms, jme, kms, kme,        &
                                     its, ite, jts, jte, kts, kte        )

!*****************************************
!*****************************************
!MODIFICATIONS for water vapor flux
!
!
  qflux: SELECT CASE( config_flags%isfflx )
  CASE (0)
! do nothing
  CASE (1,2) ! with surface moisture flux 
    IF ( im == P_QV ) THEN
       DO j = j_start, j_end
       DO i = i_start, i_end
          moist_flux = qfx(i,j)/rho(i,1,j)/(1.+moist(i,kts,j,P_QV))
          moist_tendf(i,kts,j,im)=moist_tendf(i,kts,j,im)  &
               +mu(i,j)*moist_flux*rdzw(i,kts,j)
       ENDDO
       ENDDO
    ENDIF

  CASE DEFAULT
    CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
  END SELECT qflux
!
! END MODIFICATIONS for water vapor flux
!*****************************************
!*****************************************
     ENDDO moist_loop

   ENDIF

   IF (n_chem .ge. PARAM_FIRST_SCALAR) THEN 

     chem_loop: do im = PARAM_FIRST_SCALAR, n_chem

          CALL vertical_diffusion_s( chem_tendf(ims,kms,jms,im),         &
                                     config_flags, chem(ims,kms,jms,im), &
                                     mu, xkhv,                             &
                                     dn, dnw, rdz, rdzw, fnm, fnp,         &
                                     .false.,                              &
                                     ids, ide, jds, jde, kds, kde,         &
                                     ims, ime, jms, jme, kms, kme,         &
                                     its, ite, jts, jte, kts, kte         )
     ENDDO chem_loop

   ENDIF

   IF (n_tracer .ge. PARAM_FIRST_SCALAR) THEN 

     tracer_loop: do im = PARAM_FIRST_SCALAR, n_tracer

          CALL vertical_diffusion_s( tracer_tendf(ims,kms,jms,im),         &
                                     config_flags, tracer(ims,kms,jms,im), &
                                     mu, xkhv,                             &
                                     dn, dnw, rdz, rdzw, fnm, fnp,         &
                                     .false.,                              &
                                     ids, ide, jds, jde, kds, kde,         &
                                     ims, ime, jms, jme, kms, kme,         &
                                     its, ite, jts, jte, kts, kte         )
     ENDDO tracer_loop

   ENDIF


   IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN 

     scalar_loop: do im = PARAM_FIRST_SCALAR, n_scalar

          CALL vertical_diffusion_s( scalar_tendf(ims,kms,jms,im),         &
                                     config_flags, scalar(ims,kms,jms,im), &
                                     mu, xkhv,                             &
                                     dn, dnw, rdz, rdzw, fnm, fnp,         &
                                     .false.,                              &
                                     ids, ide, jds, jde, kds, kde,         &
                                     ims, ime, jms, jme, kms, kme,         &
                                     its, ite, jts, jte, kts, kte         )
     ENDDO scalar_loop

   ENDIF

END SUBROUTINE vertical_diffusion_2

!=======================================================================
!=======================================================================

SUBROUTINE vertical_diffusion_u_2( tendency, config_flags, mu,            &
                                   defor13, xkmv,                         &
                                   nba_mij, n_nba_mij,                    & !JDM
                                   dnw, rdzw, fnm, fnp,                   &
                                   ids, ide, jds, jde, kds, kde,          &
                                   ims, ime, jms, jme, kms, kme,          &
                                   its, ite, jts, jte, kts, kte          )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

   INTEGER ,         INTENT(IN   ) ::       ids, ide, jds, jde, kds, kde, &
                                            ims, ime, jms, jme, kms, kme, &
                                            its, ite, jts, jte, kts, kte

   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnm
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnp
   REAL , DIMENSION( kms:kme ) ,            INTENT(IN   )      :: dnw
!   REAL , DIMENSION( ims:ime , jms:jme ) ,  INTENT(IN   )      :: zeta_z

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency

   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) ,                       &
                                            INTENT(IN   )      ::defor13, &
                                                                    xkmv, &
                                                                      rdzw

   INTEGER, INTENT(  IN ) :: n_nba_mij !JDM

   REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) &   !JDM    
   :: nba_mij

   REAL , DIMENSION( ims:ime , jms:jme ) ,  INTENT(IN   )      :: mu

! LOCAL VARS

   INTEGER :: i, j, k, ktf

   INTEGER :: i_start, i_end, j_start, j_end
   INTEGER :: is_ext,ie_ext,js_ext,je_ext  

   REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1)        :: titau3

   REAL , DIMENSION( its:ite, jts:jte)                         ::  zzavg

   REAL :: rdzu

! End declarations.
!-----------------------------------------------------------------------

   ktf=MIN(kte,kde-1)
  
   i_start = its
   i_end   = ite
   j_start = jts
   j_end   = MIN(jte,jde-1)

   IF ( config_flags%open_xs .or. config_flags%specified .or. &
        config_flags%nested) i_start = MAX(ids+1,its)
   IF ( config_flags%open_xe .or. config_flags%specified .or. &
        config_flags%nested) i_end   = MIN(ide-1,ite)
   IF ( config_flags%open_ys .or. config_flags%specified .or. &
        config_flags%nested) j_start = MAX(jds+1,jts)
   IF ( config_flags%open_ye .or. config_flags%specified .or. &
        config_flags%nested) j_end   = MIN(jde-2,jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = ite

! titau3 = titau13
   is_ext=0
   ie_ext=0
   js_ext=0
   je_ext=0
   CALL cal_titau_13_31( config_flags, titau3, defor13,   &
                         nba_mij(ims,kms,jms,P_m13),      & !JDM
                         mu, xkmv, fnm, fnp,              &
                         is_ext, ie_ext, js_ext, je_ext,  &
                         ids, ide, jds, jde, kds, kde,    &
                         ims, ime, jms, jme, kms, kme,    &
                         its, ite, jts, jte, kts, kte     )
!
      DO j = j_start, j_end
      DO k=kts+1,ktf
      DO i = i_start, i_end

         rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
         tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j)-titau3(i,k,j))

      ENDDO
      ENDDO
      ENDDO

! ******** MODIF...
!  we will pick up the surface drag (titau3(i,kts,j)) later
!
       DO j = j_start, j_end
       k=kts
       DO i = i_start, i_end
 
          rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
          tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j))
       ENDDO
       ENDDO
! ******** MODIF...

END SUBROUTINE vertical_diffusion_u_2

!=======================================================================
!=======================================================================

SUBROUTINE vertical_diffusion_v_2( tendency, config_flags, mu,            &
                                   defor23, xkmv,                         &
                                   nba_mij, n_nba_mij,                    & !JDM
                                   dnw, rdzw, fnm, fnp,                   &
                                   ids, ide, jds, jde, kds, kde,          &
                                   ims, ime, jms, jme, kms, kme,          &
                                   its, ite, jts, jte, kts, kte          )
!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

   INTEGER ,         INTENT(IN   ) ::       ids, ide, jds, jde, kds, kde, &
                                            ims, ime, jms, jme, kms, kme, &
                                            its, ite, jts, jte, kts, kte
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnm
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnp
   REAL , DIMENSION( kms:kme ) ,            INTENT(IN   )      :: dnw
!   REAL , DIMENSION( ims:ime , jms:jme ) ,  INTENT(IN   )      :: zeta_z

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency

   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) ,                       &
                                            INTENT(IN   )      ::defor23, &
                                                                    xkmv, &
                                                                    rdzw

   INTEGER, INTENT(  IN ) :: n_nba_mij !JDM

   REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij),  INTENT(INOUT) &     !JDM  
   :: nba_mij

   REAL , DIMENSION( ims:ime , jms:jme ) ,  INTENT(IN   )      :: mu

! LOCAL VARS

   INTEGER :: i, j, k, ktf

   INTEGER :: i_start, i_end, j_start, j_end
   INTEGER :: is_ext,ie_ext,js_ext,je_ext  

   REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1)        :: titau3

   REAL , DIMENSION( its:ite, jts:jte)                         ::  zzavg

   REAL  :: rdzv

! End declarations.
!-----------------------------------------------------------------------

   ktf=MIN(kte,kde-1)
  
   i_start = its
   i_end   = MIN(ite,ide-1)
   j_start = jts
   j_end   = jte

   IF ( config_flags%open_xs .or. config_flags%specified .or. &
        config_flags%nested) i_start = MAX(ids+1,its)
   IF ( config_flags%open_xe .or. config_flags%specified .or. &
        config_flags%nested) i_end   = MIN(ide-2,ite)
   IF ( config_flags%open_ys .or. config_flags%specified .or. &
        config_flags%nested) j_start = MAX(jds+1,jts)
   IF ( config_flags%open_ye .or. config_flags%specified .or. &
        config_flags%nested) j_end   = MIN(jde-1,jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)

! titau3 = titau23
   is_ext=0
   ie_ext=0
   js_ext=0
   je_ext=0
   CALL cal_titau_23_32( config_flags, titau3, defor23,   &
                         nba_mij(ims,kms,jms,P_m23),      & !JDM
                         mu, xkmv, fnm, fnp,              &
                         is_ext, ie_ext, js_ext, je_ext,  &
                         ids, ide, jds, jde, kds, kde,    &
                         ims, ime, jms, jme, kms, kme,    &
                         its, ite, jts, jte, kts, kte     )

   DO j = j_start, j_end
   DO k = kts+1,ktf
   DO i = i_start, i_end

      rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
      tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j)-titau3(i,k,j))

   ENDDO
   ENDDO
   ENDDO

! ******** MODIF...
!  we will pick up the surface drag (titau3(i,kts,j)) later
!
       DO j = j_start, j_end
       k=kts
       DO i = i_start, i_end
 
          rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
          tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j))
 
       ENDDO
       ENDDO
! ******** MODIF...

END SUBROUTINE vertical_diffusion_v_2

!=======================================================================
!=======================================================================

SUBROUTINE vertical_diffusion_w_2(tendency, config_flags, mu,             &
                                defor33, tke,                             &
                                nba_mij, n_nba_mij,                       & !JDM
                                div, xkmv,                                &
                                dn, rdz,                                  &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                its, ite, jts, jte, kts, kte              )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

   INTEGER ,         INTENT(IN   ) ::       ids, ide, jds, jde, kds, kde, &
                                            ims, ime, jms, jme, kms, kme, &
                                            its, ite, jts, jte, kts, kte

   REAL , DIMENSION( kms:kme ) ,            INTENT(IN   )      ::  dn

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency

   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) ,                       &
                                            INTENT(IN   )      ::defor33, &
                                                                     tke, &
                                                                     div, &
                                                                    xkmv, &
                                                                     rdz

   INTEGER, INTENT(  IN ) :: n_nba_mij !JDM  

   REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij),  INTENT(INOUT) &   !JDM      
   :: nba_mij

   REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN   ) :: mu

! LOCAL VARS

   INTEGER :: i, j, k, ktf

   INTEGER :: i_start, i_end, j_start, j_end
   INTEGER :: is_ext,ie_ext,js_ext,je_ext  

   REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1)        :: titau3

! End declarations.
!-----------------------------------------------------------------------

   ktf=MIN(kte,kde-1)
  
   i_start = its
   i_end   = MIN(ite,ide-1)
   j_start = jts
   j_end   = MIN(jte,jde-1)

   IF ( config_flags%open_xs .or. config_flags%specified .or. &
        config_flags%nested) i_start = MAX(ids+1,its)
   IF ( config_flags%open_xe .or. config_flags%specified .or. &
        config_flags%nested) i_end   = MIN(ide-2,ite)
   IF ( config_flags%open_ys .or. config_flags%specified .or. &
        config_flags%nested) j_start = MAX(jds+1,jts)
   IF ( config_flags%open_ye .or. config_flags%specified .or. &
        config_flags%nested) j_end   = MIN(jde-2,jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)

! titau3 = titau33
   is_ext=0
   ie_ext=0
   js_ext=0
   je_ext=0
   CALL cal_titau_11_22_33( config_flags, titau3,            &
                            mu, tke, xkmv, defor33,          &
                            nba_mij(ims,kms,jms,P_m33),      & !JDM
                            is_ext, ie_ext, js_ext, je_ext,  &
                            ids, ide, jds, jde, kds, kde,    &
                            ims, ime, jms, jme, kms, kme,    &
                            its, ite, jts, jte, kts, kte     )

!   DO j = j_start, j_end
!   DO k = kts+1, ktf
!   DO i = i_start, i_end
!      titau3(i,k,j)=titau3(i,k,j)*zeta_z(i,j)
!   ENDDO
!   ENDDO
!   ENDDO

   DO j = j_start, j_end
   DO k = kts+1, ktf
   DO i = i_start, i_end
      tendency(i,k,j)=tendency(i,k,j)-rdz(i,k,j)*(titau3(i,k,j)-titau3(i,k-1,j))
   ENDDO
   ENDDO
   ENDDO

END SUBROUTINE vertical_diffusion_w_2

!=======================================================================
!=======================================================================

SUBROUTINE vertical_diffusion_s( tendency, config_flags, var, mu, xkhv,   &
                                 dn, dnw, rdz, rdzw, fnm, fnp,            &
                                 doing_tke,                               &
                                 ids, ide, jds, jde, kds, kde,            &
                                 ims, ime, jms, jme, kms, kme,            &
                                 its, ite, jts, jte, kts, kte            )

!-----------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

   INTEGER ,         INTENT(IN   ) ::       ids, ide, jds, jde, kds, kde, &
                                            ims, ime, jms, jme, kms, kme, &
                                            its, ite, jts, jte, kts, kte

   LOGICAL,         INTENT(IN   ) ::        doing_tke

   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnm
   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) ::    fnp
   REAL , DIMENSION( kms:kme ) ,            INTENT(IN   )      ::  dn
   REAL , DIMENSION( kms:kme ) ,            INTENT(IN   )      :: dnw

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency

   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(IN) ::   xkhv

   REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) ::   mu

   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) ,                       &
                                            INTENT(IN   )      ::    var, &
                                                                     rdz, &
                                                                    rdzw
! LOCAL VARS

   INTEGER :: i, j, k, ktf

   INTEGER :: i_start, i_end, j_start, j_end

   REAL , DIMENSION( its:ite, kts:kte, jts:jte)            ::        H3, &
                                                                 xkxavg, &
                                                                  rravg

   REAL , DIMENSION( its:ite, kts:kte, jts:jte)            ::  tmptendf

! End declarations.
!-----------------------------------------------------------------------

   ktf=MIN(kte,kde-1)
  
   i_start = its
   i_end   = MIN(ite,ide-1)
   j_start = jts
   j_end   = MIN(jte,jde-1)

   IF ( config_flags%open_xs .or. config_flags%specified .or. &
        config_flags%nested) i_start = MAX(ids+1,its)
   IF ( config_flags%open_xe .or. config_flags%specified .or. &
        config_flags%nested) i_end   = MIN(ide-2,ite)
   IF ( config_flags%open_ys .or. config_flags%specified .or. &
        config_flags%nested) j_start = MAX(jds+1,jts)
   IF ( config_flags%open_ye .or. config_flags%specified .or. &
        config_flags%nested) j_end   = MIN(jde-2,jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)

   IF (doing_tke) THEN
      DO j = j_start, j_end
      DO k = kts,ktf
      DO i = i_start, i_end
         tmptendf(i,k,j)=tendency(i,k,j)
      ENDDO
      ENDDO
      ENDDO
   ENDIF

! H3

   xkxavg = 0.

   DO j = j_start, j_end
   DO k = kts+1,ktf
   DO i = i_start, i_end
      xkxavg(i,k,j)=fnm(k)*xkhv(i,k,j)+fnp(k)*xkhv(i,k-1,j)
      H3(i,k,j)=-xkxavg(i,k,j)*(var(i,k,j)-var(i,k-1,j))*rdz(i,k,j)
!      H3(i,k,j)=-xkxavg(i,k,j)*zeta_z(i,j)* &
!                 (var(i,k,j)-var(i,k-1,j))/dn(k)
   ENDDO
   ENDDO
   ENDDO

   DO j = j_start, j_end
   DO i = i_start, i_end
      H3(i,kts,j)=0.
      H3(i,ktf+1,j)=0.
!      H3(i,kts,j)=H3(i,kts+1,j)
!      H3(i,ktf+1,j)=H3(i,ktf,j)
   ENDDO
   ENDDO

   DO j = j_start, j_end
   DO k = kts,ktf
   DO i = i_start, i_end
      tendency(i,k,j)=tendency(i,k,j)  &
                       -mu(i,j)*(H3(i,k+1,j)-H3(i,k,j))*rdzw(i,k,j)
   ENDDO
   ENDDO
   ENDDO

   IF (doing_tke) THEN
      DO j = j_start, j_end
      DO k = kts,ktf
      DO i = i_start, i_end
          tendency(i,k,j)=tmptendf(i,k,j)+2.* &
                          (tendency(i,k,j)-tmptendf(i,k,j))
      ENDDO
      ENDDO
      ENDDO
   ENDIF

END SUBROUTINE vertical_diffusion_s

!=======================================================================
!=======================================================================

    SUBROUTINE cal_titau_11_22_33( config_flags, titau,              &
                                   mu, tke, xkx, defor,              &
                                   mtau,                             & !JDM
                                   is_ext, ie_ext, js_ext, je_ext,   &
                                   ids, ide, jds, jde, kds, kde,     &
                                   ims, ime, jms, jme, kms, kme,     &
                                   its, ite, jts, jte, kts, kte      )

! History:     Sep 2003  Changes by George Bryan and Jason Knievel, NCAR
!              Oct 2001  Converted to mass core by Bill Skamarock, NCAR
!              Aug 2000  Original code by Shu-Hua Chen, UC-Davis

! Purpose:     This routine calculates stress terms (taus) for use in
!              the calculation of production of TKE by sheared wind

! References:  Klemp and Wilhelmson (JAS 1978)
!              Deardorff (B-L Meteor 1980)
!              Chen and Dudhia (NCAR WRF physics report 2000)

! Key:

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: ids, ide, jds, jde, kds, kde,  &
       ims, ime, jms, jme, kms, kme,  &
       its, ite, jts, jte, kts, kte

    INTEGER, INTENT( IN )  &
    :: is_ext, ie_ext, js_ext, je_ext  

    REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT )  &
    :: titau 

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
    :: defor, xkx, tke

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &  !JDM
    :: mtau                            

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    :: mu

! Local variables.

    INTEGER  &
    :: i, j, k, ktf, i_start, i_end, j_start, j_end

! End declarations.
!-----------------------------------------------------------------------

    ktf = MIN( kte, kde-1 )

    i_start = its
    i_end   = ite
    j_start = jts
    j_end   = jte

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested) i_start = MAX( ids+1, its )
    IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
         config_flags%nested) i_end   = MIN( ide-1, ite )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested) j_start = MAX( jds+1, jts )
    IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
         config_flags%nested) j_end   = MIN( jde-1, jte )
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = ite

    i_start = i_start - is_ext
    i_end   = i_end   + ie_ext   
    j_start = j_start - js_ext
    j_end   = j_end   + je_ext   

    IF ( config_flags%sfs_opt .GT. 0 ) THEN ! USE NBA MODEL SFS STRESSES

      DO j = j_start, j_end
      DO k = kts, ktf
      DO i = i_start, i_end

        titau(i,k,j) = mu(i,j) * mtau(i,k,j)

      END DO
      END DO
      END DO  

    ELSE !NOT NBA 

      IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT

        DO j = j_start, j_end
        DO k = kts, ktf
        DO i = i_start, i_end

          titau(i,k,j) = - mu(i,j) * xkx(i,k,j) * defor(i,k,j)
          mtau(i,k,j) = - xkx(i,k,j) * defor(i,k,j) 
 
        END DO
        END DO
        END DO

      ELSE !NO STRESS OUTPUT 

        DO j = j_start, j_end
        DO k = kts, ktf
        DO i = i_start, i_end

          titau(i,k,j) = - mu(i,j) * xkx(i,k,j) * defor(i,k,j)

        END DO
        END DO
        END DO

      ENDIF 

    ENDIF 

    END SUBROUTINE cal_titau_11_22_33

!=======================================================================
!=======================================================================

    SUBROUTINE cal_titau_12_21( config_flags, titau,             &
                                mu, xkx, defor,                  &
                                mtau,                            & !JDM
                                is_ext, ie_ext, js_ext, je_ext,  &
                                ids, ide, jds, jde, kds, kde,    &
                                ims, ime, jms, jme, kms, kme,    &
                                its, ite, jts, jte, kts, kte     )

! History:     Sep 2003   Modifications by George Bryan and Jason Knievel, NCAR
!              Oct 2001   Converted to mass core by Bill Skamarock, NCAR
!              Aug 2000   Original code by Shu-Hua Chen, UC-Davis

! Pusrpose     This routine calculates the stress terms (taus) for use in
!              the calculation of production of TKE by sheared wind

! References:  Klemp and Wilhelmson (JAS 1978)
!              Deardorff (B-L Meteor 1980)
!              Chen and Dudhia (NCAR WRF physics report 2000)

! Key:

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: ids, ide, jds, jde, kds, kde,  &
       ims, ime, jms, jme, kms, kme,  &
       its, ite, jts, jte, kts, kte

    INTEGER, INTENT( IN )  &
    :: is_ext, ie_ext, js_ext, je_ext  

    REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT )  &
    :: titau 
 
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
    :: defor, xkx

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  & !JDM
    :: mtau                              

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    :: mu

! Local variables.

    INTEGER  &
    :: i, j, k, ktf, i_start, i_end, j_start, j_end

    REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 )  &
    :: xkxavg  

    REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 )  &
    :: muavg

! End declarations.
!-----------------------------------------------------------------------

    ktf = MIN( kte, kde-1 )

! Needs one more point in the x and y directions.

    i_start = its
    i_end   = ite
    j_start = jts
    j_end   = jte

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested ) i_start = MAX( ids+1, its )
    IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
         config_flags%nested ) i_end   = MIN( ide-1, ite )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested ) j_start = MAX( jds+1, jts )
    IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
         config_flags%nested ) j_end   = MIN( jde-1, jte )
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = ite

    i_start = i_start - is_ext
    i_end   = i_end   + ie_ext   
    j_start = j_start - js_ext
    j_end   = j_end   + je_ext   

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      xkxavg(i,k,j) = 0.25 * ( xkx(i-1,k,j  ) + xkx(i,k,j  ) +  &
                               xkx(i-1,k,j-1) + xkx(i,k,j-1) )
    END DO
    END DO
    END DO

    DO j = j_start, j_end
    DO i = i_start, i_end
      muavg(i,j) = 0.25 * ( mu(i-1,j  ) + mu(i,j  ) +  &
                            mu(i-1,j-1) + mu(i,j-1) )
    END DO
    END DO

! titau12 or titau21

    IF ( config_flags%sfs_opt .GT. 0 ) THEN ! USE NBA MODEL SFS STRESSES
 
      DO j = j_start, j_end
      DO k = kts, ktf
      DO i = i_start, i_end

        titau(i,k,j) = muavg(i,j)  * mtau(i,k,j) 

      END DO
      END DO
      END DO

    ELSE ! NOT NBA
  
      IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT

        DO j = j_start, j_end
        DO k = kts, ktf
        DO i = i_start, i_end

          titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
          mtau(i,k,j) = - xkxavg(i,k,j) * defor(i,k,j) 

        END DO
        END DO
        END DO

      ELSE ! NO STRESS OUTPUT

        DO j = j_start, j_end
        DO k = kts, ktf
        DO i = i_start, i_end

          titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j) 

        END DO
        END DO
        END DO

      ENDIF

    ENDIF 

    END SUBROUTINE cal_titau_12_21

!=======================================================================

    SUBROUTINE cal_titau_13_31( config_flags, titau,             &
                                defor,                           & 
                                mtau,                            & !JDM
                                mu, xkx, fnm, fnp,               &
                                is_ext, ie_ext, js_ext, je_ext,  &
                                ids, ide, jds, jde, kds, kde,    &
                                ims, ime, jms, jme, kms, kme,    &
                                its, ite, jts, jte, kts, kte     )

! History:     Sep 2003   Modifications by George Bryan and Jason Knievel, NCAR
!              Oct 2001   Converted to mass core by Bill Skamarock, NCAR
!              Aug 2000   Original code by Shu-Hua Chen, UC-Davis

! Purpose:     This routine calculates the stress terms (taus) for use in
!              the calculation of production of TKE by sheared wind

! References:  Klemp and Wilhelmson (JAS 1978)
!              Deardorff (B-L Meteor 1980)
!              Chen and Dudhia (NCAR WRF physics report 2000)

! Key:

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: ids, ide, jds, jde, kds, kde,  &
       ims, ime, jms, jme, kms, kme,  &
       its, ite, jts, jte, kts, kte

    INTEGER, INTENT( IN )  &
    :: is_ext, ie_ext, js_ext, je_ext  

    REAL, DIMENSION( kms:kme ), INTENT( IN )  &
    :: fnm, fnp

    REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT )  &
    :: titau 
 
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme), INTENT( IN )  &
    :: defor, xkx

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  & !JDM
    :: mtau                                      

    REAL, DIMENSION( ims:ime, jms:jme), INTENT( IN )  &
    :: mu

! Local variables.

    INTEGER  &
    :: i, j, k, ktf, i_start, i_end, j_start, j_end

    REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 )  &
    :: xkxavg 

    REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 )  &
    :: muavg

! End declarations.
!-----------------------------------------------------------------------

    ktf = MIN( kte, kde-1 )

! Find ide-1 and jde-1 for averaging to p point.

    i_start = its
    i_end   = ite
    j_start = jts
    j_end   = MIN( jte, jde-1 )

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested) i_start = MAX( ids+1, its )
    IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
         config_flags%nested) i_end   = MIN( ide-1, ite )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested) j_start = MAX( jds+1, jts )
    IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
         config_flags%nested) j_end   = MIN( jde-2, jte )
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = ite

    i_start = i_start - is_ext
    i_end   = i_end   + ie_ext   
    j_start = j_start - js_ext
    j_end   = j_end   + je_ext   

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k  ,j) + xkx(i-1,k  ,j) ) +  &
                              fnp(k) * ( xkx(i,k-1,j) + xkx(i-1,k-1,j) ) )
    END DO
    END DO
    END DO

    DO j = j_start, j_end
    DO i = i_start, i_end
      muavg(i,j) = 0.5 * ( mu(i,j) + mu(i-1,j) )
    END DO
    END DO

    IF ( config_flags%sfs_opt .GT. 0 ) THEN ! USE NBA MODEL SFS STRESSES
 
      DO j = j_start, j_end
      DO k = kts+1, ktf
      DO i = i_start, i_end
         titau(i,k,j) = muavg(i,j) * mtau(i,k,j) 
      ENDDO
      ENDDO
      ENDDO

    ELSE ! NOT NBA
 
      IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT

        DO j = j_start, j_end
        DO k = kts+1, ktf
        DO i = i_start, i_end

          titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
          mtau(i,k,j) = - xkxavg(i,k,j) * defor(i,k,j)
 
        ENDDO
        ENDDO
        ENDDO

      ELSE ! NO STRESS OUTPUT

        DO j = j_start, j_end
        DO k = kts+1, ktf
        DO i = i_start, i_end

          titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j) 

        ENDDO
        ENDDO
        ENDDO

      ENDIF  

    ENDIF 

    DO j = j_start, j_end
    DO i = i_start, i_end
      titau(i,kts  ,j) = 0.0
      titau(i,ktf+1,j) = 0.0
    ENDDO
    ENDDO

    END SUBROUTINE cal_titau_13_31

!=======================================================================
!=======================================================================

    SUBROUTINE cal_titau_23_32( config_flags, titau, defor,      &
                                mtau,                            & !JDM 
                                mu, xkx, fnm, fnp,               &
                                is_ext, ie_ext, js_ext, je_ext,  &
                                ids, ide, jds, jde, kds, kde,    &
                                ims, ime, jms, jme, kms, kme,    &
                                its, ite, jts, jte, kts, kte     )

! History:     Sep 2003  Changes by George Bryan and Jason Knievel, NCAR
!              Oct 2001  Converted to mass core by Bill Skamarock, NCAR
!              Aug 2000  Original code by Shu-Hua Chen, UC-Davis

! Purpose:     This routine calculates stress terms (taus) for use in
!              the calculation of production of TKE by sheared wind

! References:  Klemp and Wilhelmson (JAS 1978)
!              Deardorff (B-L Meteor 1980)
!              Chen and Dudhia (NCAR WRF physics report 2000)

! Key:

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: ids, ide, jds, jde, kds, kde,  &
       ims, ime, jms, jme, kms, kme,  &
       its, ite, jts, jte, kts, kte

    INTEGER, INTENT( IN )  &
    :: is_ext,ie_ext,js_ext,je_ext  

    REAL, DIMENSION( kms:kme ), INTENT( IN )  &
    :: fnm, fnp

    REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT )  &  
    :: titau 
 
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
    :: defor, xkx
  
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  & !JDM
    :: mtau                                             

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    ::  mu

! Local variables.

    INTEGER  &
    :: i, j, k, ktf, i_start, i_end, j_start, j_end

    REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 )  &
    :: xkxavg 
                                                                   
    REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 )  &
    :: muavg

! End declarations.
!-----------------------------------------------------------------------

     ktf = MIN( kte, kde-1 )

! Find ide-1 and jde-1 for averaging to p point.

    i_start = its
    i_end   = MIN( ite, ide-1 )
    j_start = jts
    j_end   = jte

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested) i_start = MAX( ids+1, its )
    IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
         config_flags%nested) i_end   = MIN( ide-2, ite )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested) j_start = MAX( jds+1, jts )
    IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
         config_flags%nested) j_end   = MIN( jde-1, jte )
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )

    i_start = i_start - is_ext
    i_end   = i_end   + ie_ext   
    j_start = j_start - js_ext
    j_end   = j_end   + je_ext   

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k  ,j) + xkx(i,k  ,j-1) ) +  &
                              fnp(k) * ( xkx(i,k-1,j) + xkx(i,k-1,j-1) ) )
    END DO
    END DO
    END DO
 
    DO j = j_start, j_end
    DO i = i_start, i_end
      muavg(i,j) = 0.5 * ( mu(i,j) + mu(i,j-1) )
    END DO
    END DO
 
    IF ( config_flags%sfs_opt .EQ. 1 ) THEN ! USE NBA MODEL SFS STRESSES

      DO j = j_start, j_end
      DO k = kts+1, ktf
      DO i = i_start, i_end

        titau(i,k,j) = muavg(i,j) * mtau(i,k,j)

      END DO
      END DO
      END DO

    ELSE ! NOT NBA

      IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT

        DO j = j_start, j_end
        DO k = kts+1, ktf
        DO i = i_start, i_end

          titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
          mtau(i,k,j) = - xkxavg(i,k,j) * defor(i,k,j) 

        END DO
        END DO
        END DO

      ELSE ! NO STRESS OUTPUT

        DO j = j_start, j_end
        DO k = kts+1, ktf
        DO i = i_start, i_end

          titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j) 

        END DO
        END DO
        END DO

      ENDIF 

    ENDIF 

    DO j = j_start, j_end
    DO i = i_start, i_end
      titau(i,kts  ,j) = 0.0
      titau(i,ktf+1,j) = 0.0
    END DO
    END DO

    END SUBROUTINE cal_titau_23_32

!=======================================================================
!=======================================================================

SUBROUTINE phy_bc ( config_flags,div,defor11,defor22,defor33,              &
                    defor12,defor13,defor23,xkmh,xkmv,xkhh,xkhv,tke,       &
                    RUBLTEN, RVBLTEN,                                      &
                    RUCUTEN, RVCUTEN,                                      &
                    RUSHTEN, RVSHTEN,                                      &
                    ids, ide, jds, jde, kds, kde,                          &
                    ims, ime, jms, jme, kms, kme,                          &
                    ips, ipe, jps, jpe, kps, kpe,                          &
                    its, ite, jts, jte, kts, kte                           )

!------------------------------------------------------------------------------
! Begin declarations.

   IMPLICIT NONE

   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

   INTEGER ,        INTENT(IN   ) ::        ids, ide, jds, jde, kds, kde, &
                                            ims, ime, jms, jme, kms, kme, &
                                            ips, ipe, jps, jpe, kps, kpe, &
                                            its, ite, jts, jte, kts, kte

   REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::RUBLTEN, &
                                                                 RVBLTEN, &
                                                                 RUCUTEN, &
                                                                 RVCUTEN, &
                                                                 RUSHTEN, &
                                                                 RVSHTEN, &
                                                                 defor11, &
                                                                 defor22, &
                                                                 defor33, &
                                                                 defor12, &
                                                                 defor13, &
                                                                 defor23, &
                                                                    xkmh, &
                                                                    xkmv, &
                                                                    xkhh, &
                                                                    xkhv, &
                                                                     tke, &
                                                                     div

! End declarations.
!-----------------------------------------------------------------------

   IF(config_flags%bl_pbl_physics .GT. 0) THEN

        CALL set_physical_bc3d( RUBLTEN , 't', config_flags,              &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

        CALL set_physical_bc3d( RVBLTEN , 't', config_flags,              &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   ENDIF

   IF(config_flags%cu_physics .GT. 0) THEN

        CALL set_physical_bc3d( RUCUTEN , 't', config_flags,              &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

        CALL set_physical_bc3d( RVCUTEN , 't', config_flags,              &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   ENDIF

   IF(config_flags%shcu_physics .GT. 0) THEN

        CALL set_physical_bc3d( RUSHTEN , 't', config_flags,              &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

        CALL set_physical_bc3d( RVSHTEN , 't', config_flags,              &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   ENDIF

   ! move out of the conditional, below; horiz coeffs needed for 
   ! all diff_opt cases.  JM

   CALL set_physical_bc3d( xkmh    , 't', config_flags,                   &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   CALL set_physical_bc3d( xkhh    , 't', config_flags,                   &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   IF(config_flags%diff_opt .eq. 2) THEN

   CALL set_physical_bc3d( xkmv    , 't', config_flags,                   &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   CALL set_physical_bc3d( xkhv    , 't', config_flags,                   &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   CALL set_physical_bc3d( div     , 't', config_flags,                   &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   CALL set_physical_bc3d( defor11 , 't', config_flags,                   &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   CALL set_physical_bc3d( defor22 , 't', config_flags,                   &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   CALL set_physical_bc3d( defor33 , 't', config_flags,                   &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   CALL set_physical_bc3d( defor12 , 'd', config_flags,                   &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   CALL set_physical_bc3d( defor13 , 'e', config_flags,                   &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   CALL set_physical_bc3d( defor23 , 'f', config_flags,                   &
                                ids, ide, jds, jde, kds, kde,             &
                                ims, ime, jms, jme, kms, kme,             &
                                ips, ipe, jps, jpe, kps, kpe,             &
                                its, ite, jts, jte, kts, kte              )

   ENDIF

END SUBROUTINE phy_bc 

!=======================================================================
!=======================================================================

    SUBROUTINE tke_rhs( tendency, BN2, config_flags,            &
                        defor11, defor22, defor33,              &
                        defor12, defor13, defor23,              &
                        u, v, w, div, tke, mu,                  &
                        theta, p, p8w, t8w, z, fnm, fnp,        &
                        cf1, cf2, cf3, msftx, msfty,            &
                        xkmh, xkmv, xkhv,                       &
                        rdx, rdy, dx, dy, dt, zx, zy,           &
                        rdz, rdzw, dn, dnw, isotropic,          &
                        hfx, qfx, qv, ust, rho,                 &
                        ids, ide, jds, jde, kds, kde,           &
                        ims, ime, jms, jme, kms, kme,           &
                        its, ite, jts, jte, kts, kte            )

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: ids, ide, jds, jde, kds, kde,  &
       ims, ime, jms, jme, kms, kme,  &
       its, ite, jts, jte, kts, kte

    INTEGER, INTENT( IN )  :: isotropic
    REAL, INTENT( IN )  &
    :: cf1, cf2, cf3, dt, rdx, rdy, dx, dy

    REAL, DIMENSION( kms:kme ), INTENT( IN )  &
    :: fnm, fnp, dnw, dn

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    :: msftx, msfty

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &
    :: tendency

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
    :: defor11, defor22, defor33, defor12, defor13, defor23,  &
       div, BN2, tke, xkmh, xkmv, xkhv, zx, zy, u, v, w, theta,  &
       p, p8w, t8w, z, rdz, rdzw

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    :: mu

    REAL, DIMENSION ( ims:ime, jms:jme ), INTENT( IN )   &
    :: hfx, ust, qfx
    REAL, DIMENSION ( ims:ime, kms:kme, jms:jme ), INTENT ( IN ) &
    :: qv, rho

! Local variables.

    INTEGER  &
    :: i, j, k, ktf, i_start, i_end, j_start, j_end

! End declarations.
!-----------------------------------------------------------------------

    CALL tke_shear(    tendency, config_flags,                &
                       defor11, defor22, defor33,             &
                       defor12, defor13, defor23,             &
                       u, v, w, tke, ust, mu, fnm, fnp,       &
                       cf1, cf2, cf3, msftx, msfty,           &
                       xkmh, xkmv,                            &
                       rdx, rdy, zx, zy, rdz, rdzw, dnw, dn,  &
                       ids, ide, jds, jde, kds, kde,          &
                       ims, ime, jms, jme, kms, kme,          &
                       its, ite, jts, jte, kts, kte           )

    CALL tke_buoyancy( tendency, config_flags, mu,            &
                       tke, xkhv, BN2, theta, dt,             &
                       hfx, qfx, qv,  rho,                    &
                       ids, ide, jds, jde, kds, kde,          &
                       ims, ime, jms, jme, kms, kme,          &
                       its, ite, jts, jte, kts, kte           ) 

    CALL tke_dissip(   tendency, config_flags,                &
                       mu, tke, bn2, theta, p8w, t8w, z,      &
                       dx, dy,rdz, rdzw, isotropic,           &
                       msftx, msfty,                          &
                       ids, ide, jds, jde, kds, kde,          &
                       ims, ime, jms, jme, kms, kme,          &
                       its, ite, jts, jte, kts, kte           )

! Set a lower limit on TKE.

    ktf     = MIN( kte, kde-1 )
    i_start = its
    i_end   = MIN( ite, ide-1 )
    j_start = jts
    j_end   = MIN( jte, jde-1 )

    IF ( config_flags%open_xs .or. config_flags%specified .or. &
         config_flags%nested) i_start = MAX(ids+1,its)
    IF ( config_flags%open_xe .or. config_flags%specified .or. &
         config_flags%nested) i_end   = MIN(ide-2,ite)
    IF ( config_flags%open_ys .or. config_flags%specified .or. &
         config_flags%nested) j_start = MAX(jds+1,jts)
    IF ( config_flags%open_ye .or. config_flags%specified .or. &
         config_flags%nested) j_end   = MIN(jde-2,jte)
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
 
    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tendency(i,k,j) = max( tendency(i,k,j), -mu(i,j) * max( 0.0 , tke(i,k,j) ) / dt )
    END DO
    END DO
    END DO

    END SUBROUTINE tke_rhs

!=======================================================================
!=======================================================================

    SUBROUTINE tke_buoyancy( tendency, config_flags, mu,    &
                             tke, xkhv, BN2, theta, dt,     &
                             hfx, qfx, qv,  rho,            &
                             ids, ide, jds, jde, kds, kde,  &
                             ims, ime, jms, jme, kms, kme,  &
                             its, ite, jts, jte, kts, kte   )

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: ids, ide, jds, jde, kds, kde,  &
       ims, ime, jms, jme, kms, kme,  &
       its, ite, jts, jte, kts, kte

    REAL, INTENT( IN )  &
    :: dt

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &
    :: tendency

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
    :: xkhv, tke, BN2, theta 

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    :: mu

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT ( IN ) &
    :: qv, rho

    REAL, DIMENSION(ims:ime, jms:jme ), INTENT ( IN ) :: hfx, qfx
 
! Local variables.

    INTEGER  &
    :: i, j, k, ktf

    INTEGER  &
    :: i_start, i_end, j_start, j_end

    REAL :: heat_flux, heat_flux0

    REAL :: cpm

! End declarations.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Add to the TKE tendency the term that accounts for production of TKE
! due to buoyant motions.

    ktf     = MIN( kte, kde-1 )
    i_start = its
    i_end   = MIN( ite, ide-1 )
    j_start = jts
    j_end   = MIN( jte, jde-1 )

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested ) i_start = MAX( ids+1, its )
    IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
         config_flags%nested ) i_end   = MIN( ide-2, ite )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested ) j_start = MAX( jds+1, jts )
    IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
         config_flags%nested ) j_end   = MIN( jde-2, jte )
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
 
    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end
      tendency(i,k,j) = tendency(i,k,j) - mu(i,j) * xkhv(i,k,j) * BN2(i,k,j)
    END DO
    END DO
    END DO

! MARTA: change in the computation of the tke's tendency  at the surface.
!  the buoyancy flux is the average of the surface heat flux (0.06) and the
!   flux at the first w level
!
! WCS 040331

  hflux: SELECT CASE( config_flags%isfflx )
  CASE (0,2) ! with fixed surface heat flux given in the namelist
   heat_flux0 = config_flags%tke_heat_flux  ! constant heat flux value
                                            ! set in namelist.input
! LES mods
   K=KTS
   DO j = j_start, j_end
   DO i = i_start, i_end 
      heat_flux = heat_flux0 
      tendency(i,k,j)= tendency(i,k,j) - &
                   mu(i,j)*((xkhv(i,k,j)*BN2(i,k,j))- (g/theta(i,k,j))*heat_flux)/2.

   ENDDO
   ENDDO   

  CASE (1) ! use surface heat flux computed from surface routine
   K=KTS
   DO j = j_start, j_end
   DO i = i_start, i_end 
      cpm = cp * (1. + 0.8*qv(i,k,j))
      heat_flux = (hfx(i,j)/cpm)/rho(i,k,j)

      tendency(i,k,j)= tendency(i,k,j) - &
                   mu(i,j)*((xkhv(i,k,j)*BN2(i,k,j))- (g/theta(i,k,j))*heat_flux)/2.

   ENDDO
   ENDDO   

  CASE DEFAULT
    CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
  END SELECT hflux

! end of MARTA/WCS change

! The tendency array now includes production of TKE from buoyant
! motions.
!-----------------------------------------------------------------------

    END SUBROUTINE tke_buoyancy

!=======================================================================
!=======================================================================

    SUBROUTINE tke_dissip( tendency, config_flags,            &
                           mu, tke, bn2, theta, p8w, t8w, z,  &
                           dx, dy, rdz, rdzw, isotropic,      &
                           msftx, msfty,                      &
                           ids, ide, jds, jde, kds, kde,      &
                           ims, ime, jms, jme, kms, kme,      &
                           its, ite, jts, jte, kts, kte       )

! History:     Sep 2003  Changes by George Bryan and Jason Knievel, NCAR
!              Oct 2001  Converted to mass core by Bill Skamarock, NCAR
!              Aug 2000  Original code by Shu-Hua Chen, UC-Davis

! Purpose:     This routine calculates dissipation of turbulent kinetic
!              energy.

! References:  Klemp and Wilhelmson (JAS 1978)
!              Deardorff (B-L Meteor 1980)
!              Chen and Dudhia (NCAR WRF physics report 2000)

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    ::  ids, ide, jds, jde, kds, kde,  &
        ims, ime, jms, jme, kms, kme,  &
        its, ite, jts, jte, kts, kte

    INTEGER, INTENT( IN )  :: isotropic
    REAL, INTENT( IN )  &
    :: dx, dy
 
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &
    :: tendency
 
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
    :: tke, bn2, theta, p8w, t8w, z, rdz, rdzw

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    :: mu

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    :: msftx, msfty
! Local variables.

    REAL, DIMENSION( its:ite, kts:kte, jts:jte )  &
    :: dthrdn

    REAL, DIMENSION( its:ite, kts:kte, jts:jte )  &
    :: l_scale

    REAL, DIMENSION( its:ite )  & 
    :: sumtke,  sumtkez

    INTEGER  &
    :: i, j, k, ktf, i_start, i_end, j_start, j_end

    REAL  &
    :: disp_len, deltas, coefc, tmpdz, len_s, thetasfc,  &
       thetatop, len_0, tketmp, tmp, ce1, ce2, c_k

! End declarations.
!-----------------------------------------------------------------------
    c_k = config_flags%c_k

    ce1 = ( c_k / 0.10 ) * 0.19
    ce2 = max( 0.0 , 0.93 - ce1 )

    ktf     = MIN( kte, kde-1 )
    i_start = its
    i_end   = MIN(ite,ide-1)
    j_start = jts
    j_end   = MIN(jte,jde-1)

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested) i_start = MAX( ids+1, its )
    IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
         config_flags%nested) i_end   = MIN( ide-2, ite )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested) j_start = MAX( jds+1, jts )
    IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
         config_flags%nested) j_end   = MIN( jde-2, jte )
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )

      CALL calc_l_scale( config_flags, tke, BN2, l_scale,      &
                         i_start, i_end, ktf, j_start, j_end,  &
                         dx, dy, rdzw, msftx, msfty,           &
                         ids, ide, jds, jde, kds, kde,         &
                         ims, ime, jms, jme, kms, kme,         &
                         its, ite, jts, jte, kts, kte          )
      DO j = j_start, j_end
      DO k = kts, ktf
      DO i = i_start, i_end
        deltas  = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333
        tketmp  = MAX( tke(i,k,j), 1.0e-6 )

! Apply Deardorff's (1980) "wall effect" at the bottom of the domain. 
! For LES with fine grid, no need for this wall effect!

        IF ( k .eq. kts .or. k .eq. ktf ) then
          coefc = 3.9
        ELSE
          coefc = ce1 + ce2 * l_scale(i,k,j) / deltas
        END IF

        tendency(i,k,j) = tendency(i,k,j) - &
                          mu(i,j) * coefc * tketmp**1.5 / l_scale(i,k,j)
      END DO
      END DO
      END DO

    END SUBROUTINE tke_dissip

!=======================================================================
!=======================================================================

    SUBROUTINE tke_shear( tendency, config_flags,                &
                          defor11, defor22, defor33,             &
                          defor12, defor13, defor23,             &
                          u, v, w, tke, ust, mu, fnm, fnp,       &
                          cf1, cf2, cf3, msftx, msfty,           &
                          xkmh, xkmv,                            &
                          rdx, rdy, zx, zy, rdz, rdzw, dn, dnw,  &
                          ids, ide, jds, jde, kds, kde,          &
                          ims, ime, jms, jme, kms, kme,          &
                          its, ite, jts, jte, kts, kte           )

! History:     Sep 2003   Rewritten by George Bryan and Jason Knievel,
!                         NCAR
!              Oct 2001   Converted to mass core by Bill Skamarock, NCAR
!              Aug 2000   Original code by Shu-Hua Chen, UC-Davis

! Purpose:     This routine calculates the production of turbulent
!              kinetic energy by stresses due to sheared wind.

! References:  Klemp and Wilhelmson (JAS 1978)
!              Deardorff (B-L Meteor 1980) 
!              Chen and Dudhia (NCAR WRF physics report 2000)

! Key:

! avg          temporary working array
! cf1          
! cf2         
! cf3
! defor11      deformation term ( du/dx + du/dx )
! defor12      deformation term ( dv/dx + du/dy ); same as defor21
! defor13      deformation term ( dw/dx + du/dz ); same as defor31
! defor22      deformation term ( dv/dy + dv/dy )
! defor23      deformation term ( dw/dy + dv/dz ); same as defor32
! defor33      deformation term ( dw/dz + dw/dz )
! div          3-d divergence
! dn
! dnw
! fnm
! fnp
! msftx
! msfty
! rdx
! rdy
! tendency
! titau        tau (stress tensor) with a tilde, indicating division by
!              a map-scale factor and the fraction of the total modeled
!              atmosphere beneath a given altitude (titau = tau/m/zeta)
! tke          turbulent kinetic energy

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: ids, ide, jds, jde, kds, kde,  &
       ims, ime, jms, jme, kms, kme,  &
       its, ite, jts, jte, kts, kte

    REAL, INTENT( IN )  &
    :: cf1, cf2, cf3, rdx, rdy

    REAL, DIMENSION( kms:kme ), INTENT( IN )  &
    :: fnm, fnp, dn, dnw

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    :: msftx, msfty

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &
    :: tendency

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &  
    :: defor11, defor22, defor33, defor12, defor13, defor23,    &
       tke, xkmh, xkmv, zx, zy, u, v, w, rdz, rdzw

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    :: mu

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN )  &
    :: ust

! Local variables.

    INTEGER  &
    :: i, j, k, ktf, ktes1, ktes2,      &
       i_start, i_end, j_start, j_end,  &
       is_ext, ie_ext, js_ext, je_ext   

    REAL  &
    :: mtau

    REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 )  &
    :: avg, titau, tmp2

    REAL, DIMENSION( its:ite, kts:kte, jts:jte )  &
    :: titau12, tmp1, zxavg, zyavg

    REAL :: absU, cd0, Cd

! End declarations.
!-----------------------------------------------------------------------

    ktf    = MIN( kte, kde-1 )
    ktes1  = kte-1
    ktes2  = kte-2
   
    i_start = its
    i_end   = MIN( ite, ide-1 )
    j_start = jts
    j_end   = MIN( jte, jde-1 )

    IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
         config_flags%nested ) i_start = MAX( ids+1, its )
    IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
         config_flags%nested ) i_end   = MIN( ide-2, ite )
    IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
         config_flags%nested ) j_start = MAX( jds+1, jts )
    IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
         config_flags%nested ) j_end   = MIN( jde-2, jte )
      IF ( config_flags%periodic_x ) i_start = its
      IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      zxavg(i,k,j) = 0.25 * ( zx(i,k  ,j) + zx(i+1,k  ,j) + &
                              zx(i,k+1,j) + zx(i+1,k+1,j)  )
      zyavg(i,k,j) = 0.25 * ( zy(i,k  ,j) + zy(i,k  ,j+1) + &
                              zy(i,k+1,j) + zy(i,k+1,j+1)  )
    END DO
    END DO
    END DO

! Begin calculating production of turbulence due to shear.  The approach
! is to add together contributions from six terms, each of which is the
! square of a deformation that is then multiplied by an exchange
! coefficiant.  The same exchange coefficient is assumed for horizontal
! and vertical coefficients for some of the terms (the vertical value is
! the one used). 

! For defor11.

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tendency(i,k,j) = tendency(i,k,j) + 0.5 *  &
                        mu(i,j) * xkmh(i,k,j) * ( ( defor11(i,k,j) )**2 )
    END DO
    END DO
    END DO

! For defor22.

    DO j = j_start, j_end 
    DO k = kts, ktf
    DO i = i_start, i_end
      tendency(i,k,j) = tendency(i,k,j) + 0.5 *  &
                        mu(i,j) * xkmh(i,k,j) * ( ( defor22(i,k,j) )**2 )
    END DO
    END DO
    END DO

! For defor33.

    DO j = j_start, j_end 
    DO k = kts, ktf
    DO i = i_start, i_end
      tendency(i,k,j) = tendency(i,k,j) + 0.5 *  &
                        mu(i,j) * xkmv(i,k,j) * ( ( defor33(i,k,j) )**2 )
    END DO
    END DO
    END DO

! For defor12.

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      avg(i,k,j) = 0.25 *  &
                   ( ( defor12(i  ,k,j)**2 ) + ( defor12(i  ,k,j+1)**2 ) +  &
                     ( defor12(i+1,k,j)**2 ) + ( defor12(i+1,k,j+1)**2 ) )
    END DO
    END DO
    END DO

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmh(i,k,j) * avg(i,k,j)
    END DO
    END DO
    END DO

! For defor13.

    DO j = j_start, j_end
    DO k = kts+1, ktf
    DO i = i_start, i_end+1
      tmp2(i,k,j) = defor13(i,k,j)
    END DO
    END DO
    END DO

    DO j = j_start, j_end
    DO i = i_start, i_end+1
      tmp2(i,kts  ,j) = 0.0
      tmp2(i,ktf+1,j) = 0.0
    END DO
    END DO

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      avg(i,k,j) = 0.25 *  &
                   ( ( tmp2(i  ,k+1,j)**2 ) + ( tmp2(i  ,k,j)**2 ) +  &
                     ( tmp2(i+1,k+1,j)**2 ) + ( tmp2(i+1,k,j)**2 ) )
    END DO
    END DO
    END DO

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j)
    END DO
    END DO
    END DO

!MARTA: add the drag at the surface; WCS 040331
    K=KTS

  uflux: SELECT CASE( config_flags%isfflx )
  CASE (0) ! Assume cd a constant, specified in namelist

    cd0 = config_flags%tke_drag_coefficient  ! drag coefficient set 
                                             ! in namelist.input
    DO j = j_start, j_end   
    DO i = i_start, i_end

      absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)
      Cd = cd0
      tendency(i,k,j) = tendency(i,k,j) +       &
           mu(i,j)*( (u(i,k,j)+u(i+1,k,j))*0.5* &
                     Cd*absU*(defor13(i,kts+1,j)+defor13(i+1,kts+1,j))*0.5 )

    END DO
    END DO

  CASE (1,2) ! ustar computed from surface routine

    DO j = j_start, j_end
    DO i = i_start, i_end

      absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)+epsilon
      Cd = (ust(i,j)**2)/(absU**2)
      tendency(i,k,j) = tendency(i,k,j) +       &
           mu(i,j)*( (u(i,k,j)+u(i+1,k,j))*0.5* &
                     Cd*absU*(defor13(i,kts+1,j)+defor13(i+1,kts+1,j))*0.5 )

    END DO
    END DO

  CASE DEFAULT
    CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
  END SELECT uflux
! end of MARTA/WCS change

! For defor23.

    DO j = j_start, j_end+1
    DO k = kts+1, ktf
    DO i = i_start, i_end
      tmp2(i,k,j) = defor23(i,k,j)
    END DO
    END DO
    END DO

    DO j = j_start, j_end+1
    DO i = i_start, i_end
      tmp2(i,kts,  j) = 0.0
      tmp2(i,ktf+1,j) = 0.0
    END DO
    END DO

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      avg(i,k,j) = 0.25 *  &
                   ( ( tmp2(i,k+1,j  )**2 ) + ( tmp2(i,k,j  )**2) +  &
                     ( tmp2(i,k+1,j+1)**2 ) + ( tmp2(i,k,j+1)**2) )
    END DO
    END DO
    END DO

    DO j = j_start, j_end
    DO k = kts, ktf
    DO i = i_start, i_end
      tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j)
    END DO
    END DO
    END DO

!MARTA: add the drag at the surface; WCS 040331
    K=KTS

  vflux: SELECT CASE( config_flags%isfflx )
  CASE (0) ! Assume cd a constant, specified in namelist

    cd0 = config_flags%tke_drag_coefficient   ! constant drag coefficient 
                                              ! set in namelist.input
    DO j = j_start, j_end   
    DO i = i_start, i_end

      absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)
      Cd = cd0
      tendency(i,k,j) = tendency(i,k,j) +       &
           mu(i,j)*( (v(i,k,j)+v(i,k,j+1))*0.5* &
                     Cd*absU*(defor23(i,kts+1,j)+defor23(i,kts+1,j+1))*0.5 )

    END DO
    END DO

  CASE (1,2) ! ustar computed from surface routine

    DO j = j_start, j_end   
    DO i = i_start, i_end

      absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)+epsilon
      Cd = (ust(i,j)**2)/(absU**2)
      tendency(i,k,j) = tendency(i,k,j) +       &
           mu(i,j)*( (v(i,k,j)+v(i,k,j+1))*0.5* &
                     Cd*absU*(defor23(i,kts+1,j)+defor23(i,kts+1,j+1))*0.5 )

    END DO
    END DO

  CASE DEFAULT
    CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
  END SELECT vflux
! end of MARTA/WCS change

    END SUBROUTINE tke_shear

!=======================================================================
!=======================================================================

    SUBROUTINE compute_diff_metrics( config_flags, ph, phb, z, rdz, rdzw,  &
                                     zx, zy, rdx, rdy,                     &
                                     ids, ide, jds, jde, kds, kde,         &
                                     ims, ime, jms, jme, kms, kme,         &
                                     its, ite, jts, jte, kts, kte         )

!-----------------------------------------------------------------------
! Begin declarations.

    IMPLICIT NONE

    TYPE( grid_config_rec_type ), INTENT( IN )  &
    :: config_flags

    INTEGER, INTENT( IN )  &
    :: ids, ide, jds, jde, kds, kde,  &
       ims, ime, jms, jme, kms, kme,  &
       its, ite, jts, jte, kts, kte

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
    :: ph, phb

    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT )  &
    :: rdz, rdzw, zx, zy, z

    REAL, INTENT( IN )  &
    :: rdx, rdy

! Local variables.

    INTEGER  &
    :: i, j, k, i_start, i_end, j_start, j_end, ktf

! End declarations.
!-----------------------------------------------------------------------

    ktf = MIN( kte, kde-1 )

! Bug fix, WCS, 22 april 2002.

! We need rdzw in halo for average to u and v points.

    j_start = jts-1
    j_end   = jte

! Begin with dz computations.

    DO j = j_start, j_end

      IF ( ( j_start >= jts ) .AND. ( j_end <= MIN( jte, jde-1 ) ) ) THEN
        i_start = its-1
        i_end   = ite
      ELSE
        i_start = its
        i_end   = MIN( ite, ide-1 )
      END IF

! Compute z at w points for rdz and rdzw computations.  We'll switch z
! to z at p points before returning

      DO k = 1, kte

! Bug fix, WCS, 22 april 2002

      DO i = i_start, i_end
        z(i,k,j) = ( ph(i,k,j) + phb(i,k,j) ) / g
      END DO
      END DO

      DO k = 1, ktf
      DO i = i_start, i_end
        rdzw(i,k,j) = 1.0 / ( z(i,k+1,j) - z(i,k,j) )
      END DO
      END DO

      DO k = 2, ktf
      DO i = i_start, i_end
        rdz(i,k,j) = 2.0 / ( z(i,k+1,j) - z(i,k-1,j) )
      END DO
      END DO

! Bug fix, WCS, 22 april 2002; added the following code

      DO i = i_start, i_end
        rdz(i,1,j) = 2./(z(i,2,j)-z(i,1,j))
      END DO

    END DO

! End bug fix.

! Now compute zx and zy; we'll assume that the halo for ph and phb is
! properly filled.

    i_start = its
    i_end   = MIN( ite, ide-1 )
    j_start = jts
    j_end   = MIN( jte, jde-1 )

    DO j = j_start, j_end
    DO k = 1, kte
    DO i = MAX( ids+1, its ), i_end
      zx(i,k,j) = rdx * ( phb(i,k,j) - phb(i-1,k,j) ) / g
    END DO
    END DO
    END DO

    DO j = j_start, j_end
    DO k = 1, kte
    DO i = MAX( ids+1, its ), i_end
      zx(i,k,j) = zx(i,k,j) + rdx * ( ph(i,k,j) - ph(i-1,k,j) ) / g
    END DO
    END DO
    END DO

    DO j = MAX( jds+1, jts ), j_end
    DO k = 1, kte
    DO i = i_start, i_end
      zy(i,k,j) = rdy * ( phb(i,k,j) - phb(i,k,j-1) ) / g
    END DO
    END DO
    END DO

    DO j = MAX( jds+1, jts ), j_end
    DO k = 1, kte
    DO i = i_start, i_end
      zy(i,k,j) = zy(i,k,j) + rdy * ( ph(i,k,j) - ph(i,k,j-1) ) / g
    END DO
    END DO
    END DO

! Some b.c. on zx and zy.

    IF ( .NOT. config_flags%periodic_x ) THEN

      IF ( ite == ide ) THEN
        DO j = j_start, j_end
        DO k = 1, ktf
          zx(ide,k,j) = 0.0
        END DO
        END DO
      END IF

      IF ( its == ids ) THEN
        DO j = j_start, j_end
        DO k = 1, ktf
          zx(ids,k,j) = 0.0
        END DO
        END DO
      END IF

    ELSE

      IF ( ite == ide ) THEN
        DO j=j_start,j_end
        DO k=1,ktf
         zx(ide,k,j) = rdx * ( phb(ide,k,j) - phb(ide-1,k,j) ) / g
        END DO
        END DO

        DO j = j_start, j_end
        DO k = 1, ktf
          zx(ide,k,j) = zx(ide,k,j) + rdx * ( ph(ide,k,j) - ph(ide-1,k,j) ) / g
        END DO
        END DO
      END IF

      IF ( its == ids ) THEN
        DO j = j_start, j_end
        DO k = 1, ktf
          zx(ids,k,j) = rdx * ( phb(ids,k,j) - phb(ids-1,k,j) ) / g
        END DO
        END DO

        DO j =j_start,j_end
        DO k =1,ktf
          zx(ids,k,j) = zx(ids,k,j) + rdx * ( ph(ids,k,j) - ph(ids-1,k,j) ) / g
        END DO
        END DO
      END IF

    END IF

    IF ( .NOT. config_flags%periodic_y ) THEN

      IF ( jte == jde ) THEN
        DO k =1, ktf
        DO i =i_start, i_end
          zy(i,k,jde) = 0.0
        END DO
        END DO
      END IF

      IF ( jts == jds ) THEN
        DO k =1, ktf
        DO i =i_start, i_end
          zy(i,k,jds) = 0.0
        END DO
        END DO
      END IF

    ELSE

      IF ( jte == jde ) THEN
        DO j=j_start, j_end
        DO k=1, ktf
          zy(i,k,jde) = rdy * ( phb(i,k,jde) - phb(i,k,jde-1) ) / g
        END DO
        END DO

        DO j = j_start, j_end
        DO k = 1, ktf
          zy(i,k,jde) = zy(i,k,jde) + rdy * ( ph(i,k,jde) - ph(i,k,jde-1) ) / g
        END DO
        END DO
      END IF

      IF ( jts == jds ) THEN
        DO j = j_start, j_end
        DO k = 1, ktf
          zy(i,k,jds) = rdy * ( phb(i,k,jds) - phb(i,k,jds-1) ) / g
        END DO
        END DO

        DO j = j_start, j_end
        DO k = 1, ktf
          zy(i,k,jds) = zy(i,k,jds) + rdy * ( ph(i,k,jds) - ph(i,k,jds-1) ) / g
        END DO
        END DO
      END IF

    END IF
      
! Calculate z at p points.

    DO j = j_start, j_end
      DO k = 1, ktf
      DO i = i_start, i_end
        z(i,k,j) = 0.5 *  &
                   ( ph(i,k,j) + phb(i,k,j) + ph(i,k+1,j) + phb(i,k+1,j) ) / g
      END DO
      END DO
    END DO

    END SUBROUTINE compute_diff_metrics

!=======================================================================
!=======================================================================

    END MODULE module_diffusion_em

!=======================================================================
!=======================================================================