/wrfv2_fire/dyn_em/module_diffusion_em.F

Large files files are truncated, but you can click here to view the full file

! WRF:MODEL_LAYER:PHYSICS
 
    MODULE module_diffusion_em

    USE module_bc, only: set_physical_bc3d
    USE module_state_description, only: p_m23, p_m13, p_m22, p_m33, p_r23, p_r13, p_r12, p_m12, p_m11
    USE module_big_step_utilities_em, only: grid_config_rec_type, param_first_scalar, p_qv, p_qi, p_qc
    USE module_model_constants    

    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.
!-------…