eddy_diff.F90

  module eddy_diff

  !--------------------------------------------------------------------------------- !
  !                                                                                  !
  ! The University of Washington Moist Turbulence Scheme to compute eddy diffusion   ! 
  ! coefficients associated with dry and moist turbulences in the whole              !
  ! atmospheric layers.                                                              !
  !                                                                                  !
  ! For detailed description of the code and its performances, see                   !
  !                                                                                  !
  ! 1.'A new moist turbulence parametrization in the Community Atmosphere Model'     !
  !    by Christopher S. Bretherton and Sungsu Park. J. Climate. 2009. 22. 3422-3448 !
  ! 2.'The University of Washington shallow convection and moist turbulence schemes  !
  !    and their impact on climate simulations with the Community Atmosphere Model'  !
  !    by Sungsu Park and Christopher S. Bretherton. J. Climate. 2009. 22. 3449-3469 !
  !                                                                                  !
  ! For questions on the scheme and code, send an email to                           !
  !     Sungsu Park      at sungsup@ucar.edu (tel: 303-497-1375)                     !
  !     Chris Bretherton at breth@washington.edu                                     !
  !                                                                                  ! 
  ! Developed by Chris Bretherton at the University of Washington, Seattle, WA.      !
  !              Sungsu Park      at the CGD/NCAR, Boulder, CO.                      !
  ! Last coded on May.2006, Dec.2009 by Sungsu Park.                                 !
  !                                                                                  !  
  !--------------------------------------------------------------------------------- !

  use wv_saturation,    only: qsat

  implicit none
  private
  save

  public :: init_eddy_diff
  public :: trbintd
  public :: caleddy
  public :: ncvmax

  integer,          parameter :: r8 = selected_real_kind(12)    ! 8 byte real
  integer,          parameter :: i4 = selected_int_kind( 6)     ! 4 byte integer
  ! --------------------------------- !
  ! PBL Parameters used in the UW PBL !
  ! --------------------------------- !

  character,        parameter :: sftype         = 'l'           ! Method for calculating saturation fraction

  character(len=4), parameter :: choice_evhc    = 'maxi'        ! 'orig', 'ramp', 'maxi' : recommended to be used with choice_radf 
  character(len=6), parameter :: choice_radf    = 'maxi'        ! 'orig', 'ramp', 'maxi' : recommended to be used with choice_evhc 
  character(len=6), parameter :: choice_SRCL    = 'nonamb'      ! 'origin', 'remove', 'nonamb'
 
  character(len=6), parameter :: choice_tunl    = 'rampcl'      ! 'origin', 'rampsl'(Sungsu), 'rampcl'(Chris)
  real(r8),         parameter :: ctunl          =  2._r8        !  Maximum asympt leng = ctunl*tunl when choice_tunl = 'rampsl(cl)'
                                                                ! [ no unit ]
  character(len=6), parameter :: choice_leng    = 'origin'      ! 'origin', 'takemn'
  real(r8),         parameter :: cleng          =  3._r8        !  Order of 'leng' when choice_leng = 'origin' [ no unit ]
  character(len=6), parameter :: choice_tkes    = 'ibprod'      ! 'ibprod' (include tkes in computing bprod), 'ebprod'(exclude)

  real(r8)                    :: lbulk_max      =  40.e3_r8     ! Maximum master length scale designed to address issues in the
                                                                ! upper atmosphere where vertical model resolution is coarse [ m ].
                                                                ! In order not to disturb turbulence characteristics in the lower
                                                                ! troposphere, this should be set at least larger than ~ a few km.  
  real(r8), allocatable       :: leng_max(:)                    ! Maximum length scale designed to address issues in the upper
                                                                ! atmosphere.

  ! Parameters for 'sedimentation-entrainment feedback' for liquid stratus
  ! If .false.,  no sedimentation entrainment feedback ( i.e., use default evhc )

  logical,          parameter :: id_sedfact     = .false.
  real(r8),         parameter :: ased           =  9._r8        !  Valid only when id_sedfact = .true.

  ! --------------------------------------------------------------------------------------------------- !
  ! Parameters governing entrainment efficiency A = a1l(i)*evhc, evhc = 1 + a2l * a3l * L * ql / jt2slv !
  ! Here, 'ql' is cloud-top LWC and 'jt2slv' is the jump in 'slv' across                                !
  ! the cloud-top entrainment zone ( across two grid layers to consider full mixture )                  !
  ! --------------------------------------------------------------------------------------------------- !

  real(r8),         parameter :: a1l            =   0.10_r8     ! Dry entrainment efficiency for TKE closure
                                                                ! a1l = 0.2*tunl*erat^-1.5,
                                                                ! where erat = <e>/wstar^2 for dry CBL =  0.3.

  real(r8),         parameter :: a1i            =   0.2_r8      ! Dry entrainment efficiency for wstar closure
  real(r8),         parameter :: ccrit          =   0.5_r8      ! Minimum allowable sqrt(tke)/wstar.
                                                                ! Used in solving cubic equation for 'ebrk'
  real(r8),         parameter :: wstar3factcrit =   0.5_r8      ! 1/wstar3factcrit is the maximally allowed enhancement of
                                                                ! 'wstar3' due to entrainment.

  real(r8)                    :: a2l                            ! Moist entrainment enhancement param (recommended range : 10~30 )
  real(r8),         parameter :: a3l            =   0.8_r8      ! Approximation to a complicated thermodynamic parameters

  real(r8),         parameter :: jbumin         =   .001_r8     ! Minimum buoyancy jump at an entrainment jump, [m/s2]
  real(r8),         parameter :: evhcmax        =   10._r8      ! Upper limit of evaporative enhancement factor

  real(r8),         parameter :: onet           =   1._r8/3._r8 ! 1/3 power in wind gradient expression [ no unit ]
  integer                     :: ncvmax                         ! Max numbers of CLs (good to set to 'pver')
  real(r8),         parameter :: qmin           =   1.e-5_r8    ! Minimum grid-mean LWC counted as clouds [kg/kg]
  real(r8),         parameter :: ntzero         =   1.e-12_r8   ! Not zero (small positive number used in 's2')
  real(r8),         parameter :: b1             =   5.8_r8      ! TKE dissipation D = e^3/(b1*leng), e = b1*W.
  real(r8)                    :: b123                           ! b1**(2/3)
  real(r8),         parameter :: tunl           =   0.085_r8    ! Asympt leng = tunl*(turb lay depth)
  real(r8),         parameter :: alph1          =   0.5562_r8   ! alph1~alph5 : Galperin instability function parameters
  real(r8),         parameter :: alph2          =  -4.3640_r8   !               These coefficients are used to calculate 
  real(r8),         parameter :: alph3          = -34.6764_r8   !               'sh' and 'sm' from 'gh'.
  real(r8),         parameter :: alph4          =  -6.1272_r8   !
  real(r8),         parameter :: alph5          =   0.6986_r8   !
  real(r8),         parameter :: ricrit         =   0.19_r8     ! Critical Richardson number for turbulence.
                                                                ! Can be any value >= 0.19.
  real(r8),         parameter :: ae             =   1._r8       ! TKE transport efficiency [no unit]
  real(r8),         parameter :: rinc           =  -0.04_r8     ! Minimum W/<W> used for CL merging test 
  real(r8),         parameter :: wpertmin       =   1.e-6_r8    ! Minimum PBL eddy vertical velocity perturbation
  real(r8),         parameter :: wfac           =   1._r8       ! Ratio of 'wpert' to sqrt(tke) for CL.
  real(r8),         parameter :: tfac           =   1._r8       ! Ratio of 'tpert' to (w't')/wpert for CL.
                                                                ! Same ratio also used for q
  real(r8),         parameter :: fak            =   8.5_r8      ! Constant in surface temperature excess for stable STL.
                                                                ! [ no unit ]
  real(r8),         parameter :: rcapmin        =   0.1_r8      ! Minimum allowable e/<e> in a CL
  real(r8),         parameter :: rcapmax        =   2.0_r8      ! Maximum allowable e/<e> in a CL
  real(r8),         parameter :: tkemax         =  20._r8       ! TKE is capped at tkemax [m2/s2]

  logical,          parameter :: use_dw_surf    =  .true.       ! Used in 'zisocl'. Default is 'true'
                                                                ! If 'true', surface interfacial energy does not contribute
                                                                ! to the CL mean stability functions after finishing merging.
                                                                ! For this case, 'dl2n2_surf' is only used for a merging test
                                                                ! based on 'l2n2'
                                                                ! If 'false',surface interfacial enery explicitly contribute to
                                                                ! CL mean stability functions after finishing merging.
                                                                ! For this case, 'dl2n2_surf' and 'dl2s2_surf' are directly used
                                                                ! for calculating surface interfacial layer energetics

  logical,          parameter :: set_qrlzero    =  .false.      ! .true. ( .false.) : turning-off ( on) radiative-turbulence
                                                                ! interaction by setting qrl = 0.

  ! ------------------------------------------------------- !
  ! PBL constants set using values from other parts of code !
  ! ------------------------------------------------------- !

  real(r8)                    :: cpair                          ! Specific heat of dry air
  real(r8)                    :: rair                           ! Gas const for dry air
  real(r8)                    :: zvir                           ! rh2o/rair - 1
  real(r8)                    :: latvap                         ! Latent heat of vaporization
  real(r8)                    :: latice                         ! Latent heat of fusion
  real(r8)                    :: latsub                         ! Latent heat of sublimation
  real(r8)                    :: g                              ! Gravitational acceleration
  real(r8)                    :: vk                             ! Von Karman's constant

  integer                     :: ntop_turb                      ! Top interface level to which turbulent vertical diffusion
                                                                ! is applied ( = 1 )
  integer                     :: nbot_turb                      ! Bottom interface level to which turbulent vertical diff
                                                                ! is applied ( = pver )

  CONTAINS

  !============================================================================ !
  !                                                                             !
  !============================================================================ !
  
  subroutine init_eddy_diff( pver, gravx, cpairx, rairx, zvirx, & 
                             latvapx, laticex, ntop_eddy, nbot_eddy, vkx, &
                             eddy_lbulk_max, leng_max_in, &
                             eddy_moist_entrain_a2l, errstring)
    !---------------------------------------------------------------- ! 
    ! Purpose:                                                        !
    ! Initialize time independent constants/variables of PBL package. !
    !---------------------------------------------------------------- !
    
    ! --------- !
    ! Arguments !
    ! --------- !
    integer,  intent(in) :: pver       ! Number of vertical layers
    integer,  intent(in) :: ntop_eddy  ! Top interface level to which eddy vertical diffusivity is applied ( = 1 )
    integer,  intent(in) :: nbot_eddy  ! Bottom interface level to which eddy vertical diffusivity is applied ( = pver )
    real(r8), intent(in) :: gravx      ! Acceleration of gravity
    real(r8), intent(in) :: cpairx     ! Specific heat of dry air
    real(r8), intent(in) :: rairx      ! Gas constant for dry air
    real(r8), intent(in) :: zvirx      ! rh2o/rair - 1
    real(r8), intent(in) :: latvapx    ! Latent heat of vaporization
    real(r8), intent(in) :: laticex    ! Latent heat of fusion
    real(r8), intent(in) :: vkx        ! Von Karman's constant
    real(r8), intent(in) :: eddy_lbulk_max ! Maximum master length scale
    real(r8), intent(in) :: leng_max_in(pver) ! Maximum length scale for upper atmosphere
    real(r8), intent(in) :: eddy_moist_entrain_a2l ! Moist entrainment enhancement param

    character(len=*), intent(out) :: errstring

    integer              :: k          ! Vertical loop index

    errstring = ""

    ! --------------- !
    ! Basic constants !
    ! --------------- !

    ncvmax    = pver

    cpair     = cpairx
    rair      = rairx
    g         = gravx
    zvir      = zvirx
    latvap    = latvapx
    latice    = laticex
    latsub    = latvap + latice
    vk        = vkx
    ntop_turb = ntop_eddy
    nbot_turb = nbot_eddy
    b123      = b1**(2._r8/3._r8)
    a2l       = eddy_moist_entrain_a2l
    
    lbulk_max = eddy_lbulk_max

    allocate(leng_max(pver))
    leng_max = leng_max_in

  end subroutine init_eddy_diff

  !=============================================================================== !
  !                                                                                !
  !=============================================================================== !
  
  subroutine sfdiag( pcols   , pver    , ncol    , qt      , ql      , sl      , &
                     pi      , pm      , zi      , cld     , sfi     , sfuh    , &
                     sflh    , slslope , qtslope )
    !----------------------------------------------------------------------- ! 
    !                                                                        !
    ! Purpose: Interface for calculating saturation fractions  at upper and  ! 
    !          lower-half layers, & interfaces for use by turbulence scheme  !
    !                                                                        !
    ! Method : Various but 'l' should be chosen for consistency.             !
    !                                                                        ! 
    ! Author : B. Stevens and C. Bretherton (August 2000)                    !
    !          Sungsu Park. August 2006.                                     !
    !                       May.   2008.                                     ! 
    !                                                                        !  
    ! S.Park : The computed saturation fractions are repeatedly              !
    !          used to compute buoyancy coefficients in'trbintd' & 'caleddy'.!  
    !----------------------------------------------------------------------- !

    implicit none       

    ! --------------- !
    ! Input arguments !
    ! --------------- !

    integer,  intent(in)  :: pcols               ! Number of atmospheric columns   
    integer,  intent(in)  :: pver                ! Number of atmospheric layers   
    integer,  intent(in)  :: ncol                ! Number of atmospheric columns   

    real(r8), intent(in)  :: sl(pcols,pver)      ! Liquid water static energy [ J/kg ]
    real(r8), intent(in)  :: qt(pcols,pver)      ! Total water specific humidity [ kg/kg ]
    real(r8), intent(in)  :: ql(pcols,pver)      ! Liquid water specific humidity [ kg/kg ]
    real(r8), intent(in)  :: pi(pcols,pver+1)    ! Interface pressures [ Pa ]
    real(r8), intent(in)  :: pm(pcols,pver)      ! Layer mid-point pressures [ Pa ]
    real(r8), intent(in)  :: zi(pcols,pver+1)    ! Interface heights [ m ]
    real(r8), intent(in)  :: cld(pcols,pver)     ! Stratiform cloud fraction [ fraction ]
    real(r8), intent(in)  :: slslope(pcols,pver) ! Slope of 'sl' in each layer
    real(r8), intent(in)  :: qtslope(pcols,pver) ! Slope of 'qt' in each layer

    ! ---------------- !
    ! Output arguments !
    ! ---------------- !

    real(r8), intent(out) :: sfi(pcols,pver+1)   ! Interfacial layer saturation fraction [ fraction ]
    real(r8), intent(out) :: sfuh(pcols,pver)    ! Saturation fraction in upper half-layer [ fraction ]
    real(r8), intent(out) :: sflh(pcols,pver)    ! Saturation fraction in lower half-layer [ fraction ]

    ! --------------- !
    ! Local Variables !
    ! --------------- !

    integer               :: i                   ! Longitude index
    integer               :: k                   ! Vertical index
    integer               :: km1                 ! k-1
    integer               :: status              ! Status returned by function calls
    real(r8)              :: sltop, slbot        ! sl at top/bot of grid layer
    real(r8)              :: qttop, qtbot        ! qt at top/bot of grid layer
    real(r8)              :: tltop, tlbot  ! Liquid water temperature at top/bot of grid layer
    real(r8)              :: qxtop, qxbot        ! Sat excess at top/bot of grid layer
    real(r8)              :: qxm                 ! Sat excess at midpoint
    real(r8)              :: es               ! Saturation vapor pressure
    real(r8)              :: qs               ! Saturation spec. humidity
    real(r8)              :: cldeff(pcols,pver)  ! Effective Cloud Fraction [ fraction ]

    ! ----------------------- !
    ! Main Computation Begins ! 
    ! ----------------------- !

    sfi(1:ncol,:)    = 0._r8
    sfuh(1:ncol,:)   = 0._r8
    sflh(1:ncol,:)   = 0._r8
    cldeff(1:ncol,:) = 0._r8

    select case (sftype)
    case ('d')
       ! ----------------------------------------------------------------------- !
       ! Simply use the given stratus fraction ('horizontal' cloud partitioning) !
       ! ----------------------------------------------------------------------- !
       do k = ntop_turb + 1, nbot_turb
          km1 = k - 1
          do i = 1, ncol
             sfuh(i,k) = cld(i,k)
             sflh(i,k) = cld(i,k)
             sfi(i,k)  = 0.5_r8 * ( sflh(i,km1) + min( sflh(i,km1), sfuh(i,k) ) )
          end do
       end do
       do i = 1, ncol
          sfi(i,pver+1) = sflh(i,pver) 
       end do
    case ('l')
       ! ------------------------------------------ !
       ! Use modified stratus fraction partitioning !
       ! ------------------------------------------ !
       do k = ntop_turb + 1, nbot_turb
          km1 = k - 1
          do i = 1, ncol
             cldeff(i,k) = cld(i,k)
             sfuh(i,k)   = cld(i,k)
             sflh(i,k)   = cld(i,k)
             if( ql(i,k) .lt. qmin ) then
                 sfuh(i,k) = 0._r8
                 sflh(i,k) = 0._r8
             end if
           ! Modification : The contribution of ice should be carefully considered.
             if( choice_evhc .eq. 'ramp' .or. choice_radf .eq. 'ramp' ) then 
                 cldeff(i,k) = cld(i,k) * min( ql(i,k) / qmin, 1._r8 )
                 sfuh(i,k)   = cldeff(i,k)
                 sflh(i,k)   = cldeff(i,k)
             elseif( choice_evhc .eq. 'maxi' .or. choice_radf .eq. 'maxi' ) then 
                 cldeff(i,k) = cld(i,k)
                 sfuh(i,k)   = cldeff(i,k)
                 sflh(i,k)   = cldeff(i,k)
             endif
           ! At the stratus top, take the minimum interfacial saturation fraction
             sfi(i,k) = 0.5_r8 * ( sflh(i,km1) + min( sfuh(i,k), sflh(i,km1) ) )
           ! Modification : Currently sfi at the top and surface interfaces are set to be zero.
           !                Also, sfuh and sflh in the top model layer is set to be zero.
           !                However, I may need to set 
           !                         do i = 1, ncol
           !                            sfi(i,pver+1) = sflh(i,pver) 
           !                         end do
           !                for treating surface-based fog. 
           ! OK. I added below block similar to the other cases.
          end do
       end do
       do i = 1, ncol
          sfi(i,pver+1) = sflh(i,pver)
       end do
    case ('u')
       ! ------------------------------------------------------------------------- !
       ! Use unsaturated buoyancy - since sfi, sfuh, sflh have already been zeroed !
       ! nothing more need be done for this case.                                  !
       ! ------------------------------------------------------------------------- !
    case ('z')
       ! ------------------------------------------------------------------------- !
       ! Calculate saturation fraction based on whether the air just above or just !
       ! below the interface is saturated, i.e. with vertical cloud partitioning.  !
       ! The saturation fraction of the interfacial layer between mid-points k and !
       ! k+1 is computed by averaging the saturation fraction   of the half-layers !
       ! above and below the interface,  with a special provision   for cloud tops !
       ! (more cloud in the half-layer below than in the half-layer above).In each !
       ! half-layer, vertical partitioning of  cloud based on the slopes diagnosed !
       ! above is used.     Loop down through the layers, computing the saturation !
       ! fraction in each half-layer (sfuh for upper half, sflh for lower half).   !
       ! Once sfuh(i,k) is computed, use with sflh(i,k-1) to determine  saturation !
       ! fraction sfi(i,k) for interfacial layer k-0.5.                            !
       ! This is 'not' chosen for full consistent treatment of stratus fraction in !
       ! all physics schemes.                                                      !
       ! ------------------------------------------------------------------------- !
       do k = ntop_turb + 1, nbot_turb
          km1 = k - 1
          do i = 1, ncol
           ! Compute saturation excess at the mid-point of layer k
             sltop    = sl(i,k) + slslope(i,k) * ( pi(i,k) - pm(i,k) )      
             qttop    = qt(i,k) + qtslope(i,k) * ( pi(i,k) - pm(i,k) )
             tltop = ( sltop - g * zi(i,k) ) / cpair 
             call qsat( tltop, pi(i,k), es, qs)
             qxtop    = qttop - qs
             slbot    = sl(i,k) + slslope(i,k) * ( pi(i,k+1) - pm(i,k) )      
             qtbot    = qt(i,k) + qtslope(i,k) * ( pi(i,k+1) - pm(i,k) )
             tlbot = ( slbot - g * zi(i,k+1) ) / cpair 
             call qsat( tlbot, pi(i,k+1), es, qs)
             qxbot    = qtbot - qs
             qxm      = qxtop + ( qxbot - qxtop ) * ( pm(i,k) - pi(i,k) ) / ( pi(i,k+1) - pi(i,k) )
           ! Find the saturation fraction sfuh(i,k) of the upper half of layer k.
             if( ( qxtop .lt. 0._r8 ) .and. ( qxm .lt. 0._r8 ) ) then
                   sfuh(i,k) = 0._r8 
             else if( ( qxtop .gt. 0._r8 ) .and. ( qxm .gt. 0._r8 ) ) then
                   sfuh(i,k) = 1._r8  
             else ! Either qxm < 0 and qxtop > 0 or vice versa
                   sfuh(i,k) = max( qxtop, qxm ) / abs( qxtop - qxm )
             end if
           ! Combine with sflh(i) (still for layer k-1) to get interfac layer saturation fraction
             sfi(i,k) = 0.5_r8 * ( sflh(i,k-1) + min( sflh(i,k-1), sfuh(i,k) ) )
           ! Update sflh to be for the lower half of layer k.             
             if( ( qxbot .lt. 0._r8 ) .and. ( qxm .lt. 0._r8 ) ) then
                   sflh(i,k) = 0._r8 
             else if( ( qxbot .gt. 0._r8 ) .and. ( qxm .gt. 0._r8 ) ) then
                   sflh(i,k) = 1._r8 
             else ! Either qxm < 0 and qxbot > 0 or vice versa
                   sflh(i,k) = max( qxbot, qxm ) / abs( qxbot - qxm )
             end if
          end do  ! i
       end do ! k
       do i = 1, ncol
          sfi(i,pver+1) = sflh(i,pver)  ! Saturation fraction in the lowest half-layer. 
       end do
    end select

  return
  end subroutine sfdiag
  
  !=============================================================================== !
  !                                                                                !
  !=============================================================================== !
 
  subroutine trbintd( pcols   , pver    , ncol    ,                               &
                      z       , u       , v       ,                               &
                      t       , pmid    ,                                         &
                      s2      , n2      , ri      ,                               &
                      zi      , pi      , cld     ,                               &
                      qt      , qv      , ql      , qi      , sfi     , sfuh    , &
                      sflh    , sl      , slv     , slslope , qtslope ,           &
                      chs     , chu     , cms     , cmu     )
    !----------------------------------------------------------------------- !
    ! Purpose: Calculate buoyancy coefficients at all interfaces including   !
    !          surface. Also, computes the profiles of ( sl,qt,n2,s2,ri ).   !
    !          Note that (n2,s2,ri) are defined at each interfaces except    !
    !          surface.                                                      !
    !                                                                        !
    ! Author: B. Stevens  ( Extracted from pbldiff, August, 2000 )           !
    !         Sungsu Park ( August 2006, May. 2008 )                         !
    !----------------------------------------------------------------------- !

    implicit none

    ! --------------- !
    ! Input arguments !
    ! --------------- !

    integer,  intent(in)  :: pcols                            ! Number of atmospheric columns   
    integer,  intent(in)  :: pver                             ! Number of atmospheric layers   
    integer,  intent(in)  :: ncol                             ! Number of atmospheric columns
    real(r8), intent(in)  :: z(pcols,pver)                    ! Layer mid-point height above surface [ m ]
    real(r8), intent(in)  :: u(pcols,pver)                    ! Layer mid-point u [ m/s ]
    real(r8), intent(in)  :: v(pcols,pver)                    ! Layer mid-point v [ m/s ]
    real(r8), intent(in)  :: t(pcols,pver)                    ! Layer mid-point temperature [ K ]
    real(r8), intent(in)  :: pmid(pcols,pver)                 ! Layer mid-point pressure [ Pa ]
    real(r8), intent(in)  :: zi(pcols,pver+1)                 ! Interface height [ m ]
    real(r8), intent(in)  :: pi(pcols,pver+1)                 ! Interface pressure [ Pa ]
    real(r8), intent(in)  :: cld(pcols,pver)                  ! Stratus fraction
    real(r8), intent(in)  :: qv(pcols,pver)                   ! Water vapor specific humidity [ kg/kg ]
    real(r8), intent(in)  :: ql(pcols,pver)                   ! Liquid water specific humidity [ kg/kg ]
    real(r8), intent(in)  :: qi(pcols,pver)                   ! Ice water specific humidity [ kg/kg ]

    ! ---------------- !
    ! Output arguments !
    ! ---------------- !

    real(r8), intent(out) :: s2(pcols,pver)                   ! Interfacial ( except surface ) shear squared [ s-2 ]
    real(r8), intent(out) :: n2(pcols,pver)                   ! Interfacial ( except surface ) buoyancy frequency [ s-2 ]
    real(r8), intent(out) :: ri(pcols,pver)                   ! Interfacial ( except surface ) Richardson number, 'n2/s2'
 
    real(r8), intent(out) :: qt(pcols,pver)                   ! Total specific humidity [ kg/kg ]
    real(r8), intent(out) :: sfi(pcols,pver+1)                ! Interfacial layer saturation fraction [ fraction ]
    real(r8), intent(out) :: sfuh(pcols,pver)                 ! Saturation fraction in upper half-layer [ fraction ]
    real(r8), intent(out) :: sflh(pcols,pver)                 ! Saturation fraction in lower half-layer [ fraction ]
    real(r8), intent(out) :: sl(pcols,pver)                   ! Liquid water static energy [ J/kg ] 
    real(r8), intent(out) :: slv(pcols,pver)                  ! Liquid water virtual static energy [ J/kg ]
   
    real(r8), intent(out) :: chu(pcols,pver+1)                ! Heat buoyancy coef for dry states at all interfaces, finally.
                                                              ! [ unit ? ]
    real(r8), intent(out) :: chs(pcols,pver+1)                ! heat buoyancy coef for sat states at all interfaces, finally.
                                                              ! [ unit ? ]
    real(r8), intent(out) :: cmu(pcols,pver+1)                ! Moisture buoyancy coef for dry states at all interfaces, finally.
                                                              ! [ unit ? ]
    real(r8), intent(out) :: cms(pcols,pver+1)                ! Moisture buoyancy coef for sat states at all interfaces, finally.
                                                              ! [ unit ? ]
    real(r8), intent(out) :: slslope(pcols,pver)              ! Slope of 'sl' in each layer
    real(r8), intent(out) :: qtslope(pcols,pver)              ! Slope of 'qt' in each layer
 
    ! --------------- !
    ! Local Variables !
    ! --------------- ! 

    integer               :: i                                ! Longitude index
    integer               :: k, km1                           ! Level index
    integer               :: status                           ! Status returned by function calls

    real(r8)              :: qs(pcols,pver)                   ! Saturation specific humidity
    real(r8)              :: es(pcols,pver)                   ! Saturation vapor pressure
    real(r8)              :: gam(pcols,pver)                  ! (l/cp)*(d(qs)/dT)
    real(r8)              :: rdz                              ! 1 / (delta z) between midpoints
    real(r8)              :: dsldz                            ! 'delta sl / delta z' at interface
    real(r8)              :: dqtdz                            ! 'delta qt / delta z' at interface
    real(r8)              :: ch                               ! 'sfi' weighted ch at the interface
    real(r8)              :: cm                               ! 'sfi' weighted cm at the interface
    real(r8)              :: bfact                            ! Buoyancy factor in n2 calculations
    real(r8)              :: product                          ! Intermediate vars used to find slopes
    real(r8)              :: dsldp_a, dqtdp_a                 ! Slopes across interface above 
    real(r8)              :: dsldp_b(pcols), dqtdp_b(pcols)   ! Slopes across interface below

    ! ----------------------- !
    ! Main Computation Begins !
    ! ----------------------- !

    ! Calculate conservative scalars (qt,sl,slv) and buoyancy coefficients at the layer mid-points.
    ! Note that 'ntop_turb = 1', 'nbot_turb = pver'
    do k = ntop_turb, nbot_turb
       call qsat( t(1:ncol,k), pmid(1:ncol,k), es(1:ncol,k), qs(1:ncol,k), ncol, gam=gam(1:ncol,k))
       do i = 1, ncol
          qt(i,k)  = qv(i,k) + ql(i,k) + qi(i,k) 
          sl(i,k)  = cpair * t(i,k) + g * z(i,k) - latvap * ql(i,k) - latsub * qi(i,k)
          slv(i,k) = sl(i,k) * ( 1._r8 + zvir * qt(i,k) )
        ! Thermodynamic coefficients for buoyancy flux - in this loop these are
        ! calculated at mid-points; later,  they will be averaged to interfaces,
        ! where they will ultimately be used.  At the surface, the coefficients
        ! are taken from the lowest mid point.
          bfact    = g / ( t(i,k) * ( 1._r8 + zvir * qv(i,k) - ql(i,k) - qi(i,k) ) )
          chu(i,k) = ( 1._r8 + zvir * qt(i,k) ) * bfact / cpair
          chs(i,k) = ( ( 1._r8 + ( 1._r8 + zvir ) * gam(i,k) * cpair * t(i,k) / latvap ) / ( 1._r8 + gam(i,k) ) ) * bfact / cpair
          cmu(i,k) = zvir * bfact * t(i,k)
          cms(i,k) = latvap * chs(i,k)  -  bfact * t(i,k)
       end do
    end do

    do i = 1, ncol
       chu(i,pver+1) = chu(i,pver)
       chs(i,pver+1) = chs(i,pver)
       cmu(i,pver+1) = cmu(i,pver)
       cms(i,pver+1) = cms(i,pver)
    end do

    ! Compute slopes of conserved variables sl, qt within each layer k. 
    ! 'a' indicates the 'above' gradient from layer k-1 to layer k and 
    ! 'b' indicates the 'below' gradient from layer k   to layer k+1.
    ! We take a smaller (in absolute value)  of these gradients as the
    ! slope within layer k. If they have opposite signs,   gradient in 
    ! layer k is taken to be zero. I should re-consider whether   this
    ! profile reconstruction is the best or not.
    ! This is similar to the profile reconstruction used in the UWShCu. 

    do i = 1, ncol
     ! Slopes at endpoints determined by extrapolation
       slslope(i,pver) = ( sl(i,pver) - sl(i,pver-1) ) / ( pmid(i,pver) - pmid(i,pver-1) )
       qtslope(i,pver) = ( qt(i,pver) - qt(i,pver-1) ) / ( pmid(i,pver) - pmid(i,pver-1) )
       slslope(i,1)    = ( sl(i,2) - sl(i,1) ) / ( pmid(i,2) - pmid(i,1) )
       qtslope(i,1)    = ( qt(i,2) - qt(i,1) ) / ( pmid(i,2) - pmid(i,1) )
       dsldp_b(i)      = slslope(i,1)
       dqtdp_b(i)      = qtslope(i,1)
    end do

    do k = 2, pver - 1
       do i = 1, ncol
          dsldp_a    = dsldp_b(i)
          dqtdp_a    = dqtdp_b(i)
          dsldp_b(i) = ( sl(i,k+1) - sl(i,k) ) / ( pmid(i,k+1) - pmid(i,k) )
          dqtdp_b(i) = ( qt(i,k+1) - qt(i,k) ) / ( pmid(i,k+1) - pmid(i,k) )
          product    = dsldp_a * dsldp_b(i)
          if( product .le. 0._r8 ) then 
              slslope(i,k) = 0._r8
          else if( product .gt. 0._r8 .and. dsldp_a .lt. 0._r8 ) then 
              slslope(i,k) = max( dsldp_a, dsldp_b(i) )
          else if( product .gt. 0._r8 .and. dsldp_a .gt. 0._r8 ) then 
              slslope(i,k) = min( dsldp_a, dsldp_b(i) )
          end if
          product = dqtdp_a*dqtdp_b(i)
          if( product .le. 0._r8 ) then 
              qtslope(i,k) = 0._r8
          else if( product .gt. 0._r8 .and. dqtdp_a .lt. 0._r8 ) then 
              qtslope(i,k) = max( dqtdp_a, dqtdp_b(i) )
          else if( product .gt. 0._r8 .and. dqtdp_a .gt. 0._r8 ) then 
              qtslope(i,k) = min( dqtdp_a, dqtdp_b(i) )
          end if
       end do ! i
    end do ! k

    !  Compute saturation fraction at the interfacial layers for use in buoyancy
    !  flux computation.

    call sfdiag( pcols  , pver    , ncol    , qt      , ql      , sl      , & 
                 pi     , pmid    , zi      , cld     , sfi     , sfuh    , &
                 sflh   , slslope , qtslope )

    ! Calculate buoyancy coefficients at all interfaces (1:pver+1) and (n2,s2,ri) 
    ! at all interfaces except surface. Note 'nbot_turb = pver', 'ntop_turb = 1'.
    ! With the previous definition of buoyancy coefficients at the surface, the 
    ! resulting buoyancy coefficients at the top and surface interfaces becomes 
    ! identical to the buoyancy coefficients at the top and bottom layers. Note 
    ! that even though the dimension of (s2,n2,ri) is 'pver',  they are defined
    ! at interfaces ( not at the layer mid-points ) except the surface. 

    do k = nbot_turb, ntop_turb + 1, -1
       km1 = k - 1
       do i = 1, ncol
          rdz      = 1._r8 / ( z(i,km1) - z(i,k) )
          dsldz    = ( sl(i,km1) - sl(i,k) ) * rdz
          dqtdz    = ( qt(i,km1) - qt(i,k) ) * rdz 
          chu(i,k) = ( chu(i,km1) + chu(i,k) ) * 0.5_r8
          chs(i,k) = ( chs(i,km1) + chs(i,k) ) * 0.5_r8
          cmu(i,k) = ( cmu(i,km1) + cmu(i,k) ) * 0.5_r8
          cms(i,k) = ( cms(i,km1) + cms(i,k) ) * 0.5_r8
          ch       = chu(i,k) * ( 1._r8 - sfi(i,k) ) + chs(i,k) * sfi(i,k)
          cm       = cmu(i,k) * ( 1._r8 - sfi(i,k) ) + cms(i,k) * sfi(i,k)
          n2(i,k)  = ch * dsldz +  cm * dqtdz
          s2(i,k)  = ( ( u(i,km1) - u(i,k) )**2 + ( v(i,km1) - v(i,k) )**2) * rdz**2
          s2(i,k)  = max( ntzero, s2(i,k) )
          ri(i,k)  = n2(i,k) / s2(i,k)
       end do
    end do 
    do i = 1, ncol
       n2(i,1) = n2(i,2)
       s2(i,1) = s2(i,2)
       ri(i,1) = ri(i,2)
    end do

  return

  end subroutine trbintd

    ! ---------------------------------------------------------------------------- !
    !                                                                              !
    ! The University of Washington Moist Turbulence Scheme                         !
    !                                                                              !
    ! Authors : Chris Bretherton at the University of Washington, Seattle, WA      ! 
    !           Sungsu Park at the CGD/NCAR, Boulder, CO                           !
    !                                                                              !
    ! ---------------------------------------------------------------------------- !

    subroutine caleddy( pcols        , pver         , ncol        ,                             &
                        sl           , qt           , ql          , slv        , u            , &
                        v            , pi           , z           , zi         ,                &
                        qflx         , shflx        , slslope     , qtslope    ,                &
                        chu          , chs          , cmu         , cms        , sfuh         , &
                        sflh         , n2           , s2          , ri         , rrho         , &
                        pblh         , ustar        ,                                           &
                        kvh_in       , kvm_in       , kvh         , kvm        ,                &
                        tpert        , qpert        , qrlin       , kvf        , tke          , & 
                        wstarent     , bprod        , sprod       , minpblh    , wpert        , &
                        tkes         , went         , turbtype    ,                             &
                        kbase_o      , ktop_o       , ncvfin_o    ,                             & 
                        kbase_mg     , ktop_mg      , ncvfin_mg   ,                             & 
                        kbase_f      , ktop_f       , ncvfin_f    ,                             & 
                        wet_CL       , web_CL       , jtbu_CL     , jbbu_CL    ,                &
                        evhc_CL      , jt2slv_CL    , n2ht_CL     , n2hb_CL    , lwp_CL       , &
                        opt_depth_CL , radinvfrac_CL, radf_CL     , wstar_CL   , wstar3fact_CL, &
                        ebrk         , wbrk         , lbrk        , ricl       , ghcl         , & 
                        shcl         , smcl         ,                                           &
                        gh_a         , sh_a         , sm_a        , ri_a       , leng         , & 
                        wcap         , pblhp        , cld         , ipbl       , kpblh        , &
                        wsedl        , wsed_CL      , warnstring  , errstring)

    !--------------------------------------------------------------------------------- !
    !                                                                                  !
    ! Purpose : This is a driver routine to compute eddy diffusion coefficients        !
    !           for heat (sl), momentum (u, v), moisture (qt), and other  trace        !
    !           constituents.   This scheme uses first order closure for stable        !
    !           turbulent layers (STL). For convective layers (CL), entrainment        !
    !           closure is used at the CL external interfaces, which is coupled        !
    !           to the diagnosis of a CL regime mean TKE from the instantaneous        !
    !           thermodynamic and velocity profiles.   The CLs are diagnosed by        !
    !           extending original CL layers of moist static instability   into        !
    !           adjacent weakly stably stratified interfaces,   stopping if the        !
    !           stability is too strong.   This allows a realistic depiction of        !
    !           dry convective boundary layers with a downgradient approach.           !
    !                                                                                  !   
    ! NOTE:     This routine currently assumes ntop_turb = 1, nbot_turb = pver         !
    !           ( turbulent diffusivities computed at all interior interfaces )        !
    !           and will require modification to handle a different ntop_turb.         ! 
    !                                                                                  !
    ! Authors:  Sungsu Park and Chris Bretherton. 08/2006, 05/2008.                    !
    !                                                                                  ! 
    ! For details, see                                                                 !
    !                                                                                  !
    ! 1. 'A new moist turbulence parametrization in the Community Atmosphere Model'    !
    !     by Christopher S. Bretherton & Sungsu Park. J. Climate. 22. 3422-3448. 2009. !
    !                                                                                  !
    ! 2. 'The University of Washington shallow convection and moist turbulence schemes !
    !     and their impact on climate simulations with the Community Atmosphere Model' !
    !     by Sungsu Park & Christopher S. Bretherton. J. Climate. 22. 3449-3469. 2009. !
    !                                                                                  !
    ! For questions on the scheme and code, send an email to                           !
    !     sungsup@ucar.edu or breth@washington.edu                                     !
    !                                                                                  !
    !--------------------------------------------------------------------------------- !

    use pbl_utils, only: &
      compute_radf      ! Subroutine for computing radf    

    ! ---------------- !
    ! Inputs variables !
    ! ---------------- !

    implicit none
    integer,  intent(in) :: pcols                     ! Number of atmospheric columns   
    integer,  intent(in) :: pver                      ! Number of atmospheric layers   
    integer,  intent(in) :: ncol                      ! Number of atmospheric columns   
    real(r8), intent(in) :: u(pcols,pver)             ! U wind [ m/s ]
    real(r8), intent(in) :: v(pcols,pver)             ! V wind [ m/s ]
    real(r8), intent(in) :: sl(pcols,pver)            ! Liquid water static energy, cp * T + g * z - Lv * ql - Ls * qi [ J/kg ]
    real(r8), intent(in) :: slv(pcols,pver)           ! Liquid water virtual static energy, sl * ( 1 + 0.608 * qt ) [ J/kg ]
    real(r8), intent(in) :: qt(pcols,pver)            ! Total speccific humidity  qv + ql + qi [ kg/kg ] 
    real(r8), intent(in) :: ql(pcols,pver)            ! Liquid water specific humidity [ kg/kg ]
    real(r8), intent(in) :: pi(pcols,pver+1)          ! Interface pressures [ Pa ]
    real(r8), intent(in) :: z(pcols,pver)             ! Layer midpoint height above surface [ m ]
    real(r8), intent(in) :: zi(pcols,pver+1)          ! Interface height above surface, i.e., zi(pver+1) = 0 all over the globe
                                                      ! [ m ]
    real(r8), intent(in) :: chu(pcols,pver+1)         ! Buoyancy coeffi. unsaturated sl (heat) coef. at all interfaces.
                                                      ! [ unit ? ]
    real(r8), intent(in) :: chs(pcols,pver+1)         ! Buoyancy coeffi. saturated sl (heat) coef. at all interfaces.
                                                      ! [ unit ? ]
    real(r8), intent(in) :: cmu(pcols,pver+1)         ! Buoyancy coeffi. unsaturated qt (moisture) coef. at all interfaces
                                                      ! [ unit ? ]
    real(r8), intent(in) :: cms(pcols,pver+1)         ! Buoyancy coeffi. saturated qt (moisture) coef. at all interfaces
                                                      ! [ unit ? ]
    real(r8), intent(in) :: sfuh(pcols,pver)          ! Saturation fraction in upper half-layer [ fraction ]
    real(r8), intent(in) :: sflh(pcols,pver)          ! Saturation fraction in lower half-layer [ fraction ]
    real(r8), intent(in) :: n2(pcols,pver)            ! Interfacial (except surface) moist buoyancy frequency [ s-2 ]
    real(r8), intent(in) :: s2(pcols,pver)            ! Interfacial (except surface) shear frequency [ s-2 ]
    real(r8), intent(in) :: ri(pcols,pver)            ! Interfacial (except surface) Richardson number
    real(r8), intent(in) :: qflx(pcols)               ! Kinematic surface constituent ( water vapor ) flux [ kg/m2/s ]
    real(r8), intent(in) :: shflx(pcols)              ! Kinematic surface heat flux [ unit ? ] 
    real(r8), intent(in) :: slslope(pcols,pver)       ! Slope of 'sl' in each layer [ J/kg/Pa ]
    real(r8), intent(in) :: qtslope(pcols,pver)       ! Slope of 'qt' in each layer [ kg/kg/Pa ]
    real(r8), intent(in) :: qrlin(pcols,pver)         ! Input grid-mean LW heating rate : [ K/s ] * cpair * dp = [ W/kg*Pa ]
    real(r8), intent(in) :: wsedl(pcols,pver)         ! Sedimentation velocity of liquid stratus cloud droplet [ m/s ]
    real(r8), intent(in) :: ustar(pcols)              ! Surface friction velocity [ m/s ]
    real(r8), intent(in) :: rrho(pcols)               ! 1./bottom mid-point density. Specific volume [ m3/kg ]
    real(r8), intent(in) :: kvf(pcols,pver+1)         ! Free atmosphere eddy diffusivity [ m2/s ]
    logical,  intent(in) :: wstarent                  ! Switch for choosing wstar3 entrainment parameterization
    real(r8), intent(in) :: minpblh(pcols)            ! Minimum PBL height based on surface stress [ m ]
    real(r8), intent(in) :: kvh_in(pcols,pver+1)      ! kvh saved from last timestep or last iterative step [ m2/s ] 
    real(r8), intent(in) :: kvm_in(pcols,pver+1)      ! kvm saved from last timestep or last iterative step [ m2/s ]
    real(r8), intent(in) :: cld(pcols,pver)           ! Stratus Cloud Fraction [ fraction ]

    ! ---------------- !
    ! Output variables !
    ! ---------------- !

    real(r8), intent(out) :: kvh(pcols,pver+1)        ! Eddy diffusivity for heat, moisture, and tracers [ m2/s ]
    real(r8), intent(out) :: kvm(pcols,pver+1)        ! Eddy diffusivity for momentum [ m2/s ]
    real(r8), intent(out) :: pblh(pcols)              ! PBL top height [ m ]
    real(r8), intent(out) :: pblhp(pcols)             ! PBL top height pressure [ Pa ]
    real(r8), intent(out) :: tpert(pcols)             ! Convective temperature excess [ K ]
    real(r8), intent(out) :: qpert(pcols)             ! Convective humidity excess [ kg/kg ]
    real(r8), intent(out) :: wpert(pcols)             ! Turbulent velocity excess [ m/s ]
    real(r8), intent(out) :: tkes(pcols)              ! TKE at surface [ m2/s2 ] 
    real(r8), intent(out) :: went(pcols)              ! Entrainment rate at the PBL top interface [ m/s ] 
    real(r8), intent(out) :: tke(pcols,pver+1)        ! Turbulent kinetic energy [ m2/s2 ], 'tkes' at surface, pver+1.
    real(r8), intent(out) :: bprod(pcols,pver+1)      ! Buoyancy production [ m2/s3 ],     'bflxs' at surface, pver+1.
    real(r8), intent(out) :: sprod(pcols,pver+1)      ! Shear production [ m2/s3 ], (ustar(i)**3)/(vk*z(i,pver))
                                                      ! at surface, pver+1.
    integer(i4), intent(out) :: turbtype(pcols,pver+1) ! Turbulence type at each interface:
                                                      ! 0. = Non turbulence interface
                                                      ! 1. = Stable turbulence interface
                                                      ! 2. = CL interior interface ( if bflxs > 0, surface is this )
                                                      ! 3. = Bottom external interface of CL
                                                      ! 4. = Top external interface of CL.
                                                      ! 5. = Double entraining CL external interface 
    integer(i4), intent(out) :: ipbl(pcols)           ! If 1, PBL is CL, while if 0, PBL is STL.
    integer(i4), intent(out) :: kpblh(pcols)          ! Layer index containing PBL within or at the base interface
    real(r8), intent(out) :: wsed_CL(pcols,ncvmax)    ! Sedimentation velocity at the top of each CL [ m/s ]

    character(len=*), intent(out) :: warnstring
    character(len=*), intent(out) :: errstring

    ! --------------------------- !
    ! Diagnostic output variables !
    ! --------------------------- !

    real(r8) :: kbase_o(pcols,ncvmax)                 ! Original external base interface index of CL just after 'exacol'
    real(r8) :: ktop_o(pcols,ncvmax)                  ! Original external top  interface index of CL just after 'exacol'
    real(r8) :: ncvfin_o(pcols)                       ! Original number of CLs just after 'exacol'
    real(r8) :: kbase_mg(pcols,ncvmax)                ! kbase  just after extending-merging (after 'zisocl') but without SRCL
    real(r8) :: ktop_mg(pcols,ncvmax)                 ! ktop   just after extending-merging (after 'zisocl') but without SRCL
    real(r8) :: ncvfin_mg(pcols)                      ! ncvfin just after extending-merging (after 'zisocl') but without SRCL
    real(r8) :: kbase_f(pcols,ncvmax)                 ! Final kbase  after adding SRCL
    real(r8) :: ktop_f(pcols,ncvmax)                  ! Final ktop   after adding SRCL
    real(r8) :: ncvfin_f(pcols)                       ! Final ncvfin after adding SRCL
    real(r8) :: wet_CL(pcols,ncvmax)                  ! Entrainment rate at the CL top [ m/s ] 
    real(r8) :: web_CL(pcols,ncvmax)                  ! Entrainment rate at the CL base [ m/s ]
    real(r8) :: jtbu_CL(pcols,ncvmax)                 ! Buoyancy jump across the CL top [ m/s2 ]  
    real(r8) :: jbbu_CL(pcols,ncvmax)                 ! Buoyancy jump across the CL base [ m/s2 ]  
    real(r8) :: evhc_CL(pcols,ncvmax)                 ! Evaporative enhancement factor at the CL top
    real(r8) :: jt2slv_CL(pcols,ncvmax)               ! Jump of slv ( across two layers ) at CL top for use only in evhc [ J/kg ]
    real(r8) :: n2ht_CL(pcols,ncvmax)                 ! n2 defined at the CL top  interface
                                                      ! but using sfuh(kt)   instead of sfi(kt) [ s-2 ]
    real(r8) :: n2hb_CL(pcols,ncvmax)                 ! n2 defined at the CL base interface
                                                      ! but using sflh(kb-1) instead of sfi(kb) [ s-2 ]
    real(r8) :: lwp_CL(pcols,ncvmax)                  ! LWP in the CL top layer [ kg/m2 ]
    real(r8) :: opt_depth_CL(pcols,ncvmax)            ! Optical depth of the CL top layer
    real(r8) :: radinvfrac_CL(pcols,ncvmax)           ! Fraction of LW radiative cooling confined in the top portion of CL
    real(r8) :: radf_CL(pcols,ncvmax)                 ! Buoyancy production at the CL top due to radiative cooling [ m2/s3 ]
    real(r8) :: wstar_CL(pcols,ncvmax)                ! Convective velocity of CL including entrainment contribution finally [ m/s ]
    real(r8) :: wstar3fact_CL(pcols,ncvmax)           ! "wstar3fact" of CL. Entrainment enhancement of wstar3 (inverse)

    real(r8) :: gh_a(pcols,pver+1)                    ! Half of normalized buoyancy production, -l2n2/2e. [ no unit ]
    real(r8) :: sh_a(pcols,pver+1)                    ! Galperin instability function of heat-moisture at all interfaces [ no unit ]
    real(r8) :: sm_a(pcols,pver+1)                    ! Galperin instability function of momentum      at all interfaces [ no unit ]
    real(r8) :: ri_a(pcols,pver+1)                    ! Interfacial Richardson number                  at all interfaces [ no unit ]

    real(r8) :: ebrk(pcols,ncvmax)                    ! Net CL mean TKE [ m2/s2 ]
    real(r8) :: wbrk(pcols,ncvmax)                    ! Net CL mean normalized TKE [ m2/s2 ]
    real(r8) :: lbrk(pcols,ncvmax)                    ! Net energetic integral thickness of CL [ m ]
    real(r8) :: ricl(pcols,ncvmax)                    ! Mean Richardson number of CL ( l2n2/l2s2 )
    real(r8) :: ghcl(pcols,ncvmax)                    ! Half of normalized buoyancy production of CL                 
    real(r8) :: shcl(pcols,ncvmax)                    ! Instability function of heat and moisture of CL
    real(r8) :: smcl(pcols,ncvmax)                    ! Instability function of momentum of CL

    real(r8) :: leng(pcols,pver+1)                    ! Turbulent length scale [ m ], 0 at the surface.
    real(r8) :: wcap(pcols,pver+1)                    ! Normalized TKE [m2/s2], 'tkes/b1' at the surface and 'tke/b1' at
                                                      ! the top/bottom entrainment interfaces of CL assuming no transport.
    ! ------------------------ !
    ! Local Internal Variables !
    ! ------------------------ !

    logical :: belongcv(pcols,pver+1)                 ! True for interfaces in a CL (both interior and exterior are included)
    logical :: belongst(pcols,pver+1)                 ! True for stable turbulent layer interfaces (STL)
    logical :: in_CL                                  ! True if interfaces k,k+1 both in same CL.
    logical :: extend                                 ! True when CL is extended in zisocl
    logical :: extend_up                              ! True when CL is extended upward in zisocl
    logical :: extend_dn                              ! True when CL is extended downward in zisocl

    integer :: i                                      ! Longitude index
    integer :: k                                      ! Vertical index
    integer :: ks                                     ! Vertical index
    integer :: ncvfin(pcols)                          ! Total number of CL in column
    integer :: ncvf                                   ! Total number of CL in column prior to adding SRCL
    integer :: ncv                                    ! Index of current CL
    integer :: ncvnew                                 ! Index of added SRCL appended after regular CLs from 'zisocl'
    integer :: ncvsurf                                ! If nonzero, CL index based on surface
                                                      ! (usually 1, but can be > 1 when SRCL is based at sfc)
    integer :: kbase(pcols,ncvmax)                    ! Vertical index of CL base interface
    integer :: ktop(pcols,ncvmax)                     ! Vertical index of CL top interface
    integer :: kb, kt                                 ! kbase and ktop for current CL
    integer :: ktblw                                  ! ktop of the CL located at just below the current CL

    integer  :: ktopbl(pcols)                         ! PBL top height or interface index 
    real(r8) :: bflxs(pcols)                          ! Surface buoyancy flux [ m2/s3 ]
    real(r8) :: rcap                                  ! 'tke/ebrk' at all interfaces of CL.
                                                      ! Set to 1 at the CL entrainment interfaces
    real(r8) :: jtzm                                  ! Interface layer thickness of CL top interface [ m ]
    real(r8) :: jtsl                                  ! Jump of s_l across CL top interface [ J/kg ]
    real(r8) :: jtqt                                  ! Jump of q_t across CL top interface [ kg/kg ]
    real(r8) :: jtbu                                  ! Jump of buoyancy across CL top interface [ m/s2 ]
    real(r8) :: jtu                                   ! Jump of u across CL top interface [ m/s ]
    real(r8) :: jtv                                   ! Jump of v across CL top interface [ m/s ]
    real(r8) :: jt2slv                                ! Jump of slv ( across two layers ) at CL top for use only in evhc [ J/kg ]
    real(r8) :: radf                                  ! Buoyancy production at the CL top due to radiative cooling [ m2/s3 ]
    real(r8) :: jbzm                                  ! Interface layer thickness of CL base interface [ m ]
    real(r8) :: jbsl                                  ! Jump of s_l across CL base interface [ J/kg ]
    real(r8) :: jbqt                                  ! Jump of q_t across CL top interface [ kg/kg ]
    real(r8) :: jbbu                                  ! Jump of buoyancy across CL base interface [ m/s2 ]
    real(r8) :: jbu                                   ! Jump of u across CL base interface [ m/s ]
    real(r8) :: jbv                                   ! Jump of v across CL base interface [ m/s ]
    real(r8) :: ch                                    ! Buoyancy coefficients defined at the CL top and base interfaces
                                                      ! using CL internal
    real(r8) :: cm                                    ! sfuh(kt) and sflh(kb-1) instead of sfi(kt) and sfi(kb), respectively.
                                                      ! These are used for entrainment calculation at CL external interfaces
                                                      ! and SRCL identification.
    real(r8) :: n2ht                                  ! n2 defined at the CL top  interface
                                                      ! but using sfuh(kt)   instead of sfi(kt) [ s-2 ]
    real(r8) :: n2hb                                  ! n2 defined at the CL base interface
                                                      ! but using sflh(kb-1) instead of sfi(kb) [ s-2 ]
    real(r8) :: n2htSRCL                              ! n2 defined at the upper-half layer of SRCL.
                                                      ! This is used only for identifying SRCL.
                                                      ! n2htSRCL use SRCL internal slope sl and qt
                                                      ! as well as sfuh(kt) instead of sfi(kt) [ s-2 ]
    real(r8) :: gh                                    ! Half of normalized buoyancy production ( -l2n2/2e ) [ no unit ]
    real(r8) :: sh                                    ! Galperin instability function for heat and moisture
    real(r8) :: sm                                    ! Galperin instability function for momentum
    real(r8) :: lbulk                                 ! Depth of turbulent layer, Master length scale (not energetic length)
    real(r8) :: dzht                                  ! Thickness of top    half-layer [ m ]
    real(r8) :: dzhb                                  ! Thickness of bottom half-layer [ m ]
    real(r8) :: rootp                                 ! Sqrt(net CL-mean TKE including entrainment contribution) [ m/s ]     
    real(r8) :: evhc                                  ! Evaporative enhancement factor: (1+E)
                                                      ! with E = evap. cool. efficiency [ no unit ]
    real(r8) :: kentr                                 ! Effective entrainment diffusivity 'wet*dz', 'web*dz' [ m2/s ]
    real(r8) :: lwp                                   ! Liquid water path in the layer kt [ kg/m2 ]
    real(r8) :: opt_depth                             ! Optical depth of the layer kt [ no unit ]
    real(r8) :: radinvfrac                            ! Fraction of LW cooling in the layer kt
                                                      ! concentrated at the CL top [ no unit ]
    real(r8) :: wet                                   ! CL top entrainment rate [ m/s ]
    real(r8) :: web                                   ! CL bot entrainment rate [ m/s ]. Set to zero if CL is based at surface.
    real(r8) :: vyt                                   ! n2ht/n2 at the CL top  interface
    real(r8) :: vyb                                   ! n2hb/n2 at the CL base interface
    real(r8) :: vut                                   ! Inverse Ri (=s2/n2) at the CL top  interface
    real(r8) :: vub                                   ! Inverse Ri (=s2/n2) at the CL base interface
    real(r8) :: fact                                  ! Factor relating TKE generation to entrainment [ no unit ]
    real(r8) :: trma                                  ! Intermediate variables used for solving quadratic ( for gh from ri )
    real(r8) :: trmb                                  ! and cubic equations ( for ebrk: the net CL mean TKE )
    real(r8) :: trmc                                  !
    real(r8) :: trmp                                  !
    real(r8) :: trmq                                  !
    real(r8) :: qq                                    ! 
    real(r8) :: det                                   !
    real(r8) :: gg                                    ! Intermediate variable used for calculating stability functions of
                                                      ! SRCL or SBCL based at the surface with bflxs > 0.
    real(r8) :: dzhb5                                 ! Half thickness of the bottom-most layer of current CL regime
    real(r8) :: dzht5                                 ! Half thickness of the top-most layer of adjacent CL regime
                                                      ! just below current CL
    real(r8) :: qrlw(pcols,pver)                      ! Local grid-mean LW heating rate : [K/s] * cpair * dp = [ W/kg*Pa ]

    real(r8) :: cldeff(pcols,pver)                    ! Effective stratus fraction
    real(r8) :: qleff                                 ! Used for computing evhc
    real(r8) :: tunlramp                              ! Ramping tunl
    real(r8) :: leng_imsi                             ! For Kv = max(Kv_STL, Kv_entrain)
    real(r8) :: tke_imsi                              !
    real(r8) :: kvh_imsi                              !
    real(r8) :: kvm_imsi                              !
    real(r8) :: alph4exs                              ! For extended stability function in the stable regime
    real(r8) :: ghmin                                 !   

    real(r8) :: sedfact                               ! For 'sedimentation-entrainment feedback' 

    ! Local variables specific for 'wstar' entrainment closure

    real(r8) :: cet                                   ! Proportionality coefficient between wet and wstar3
    real(r8) :: ceb                                   ! Proportionality coefficient between web and wstar3
    real(r8) :: wstar                                 ! Convective velocity for CL [ m/s ]
    real(r8) :: wstar3                                ! Cubed convective velocity for CL [ m3/s3 ]
    real(r8) :: wstar3fact                            ! 1/(relative change of wstar^3 by entrainment)
    real(r8) :: rmin                                  ! sqrt(p)
    real(r8) :: fmin                                  ! f(rmin), where f(r) = r^3 - 3*p*r - 2q
    real(r8) :: rcrit                                 ! ccrit*wstar
    real(r8) :: fcrit                                 ! f(rcrit)
    logical     noroot                                ! True if f(r) has no root r > rcrit

    character(128) :: temp_string

    !-----------------------!
    ! Start of Main Program !
    !-----------------------!

    warnstring = ""
    errstring = ""
    
    ! Option: Turn-off LW radiative-turbulence interaction in PBL scheme
    !         by setting qrlw = 0.  Logical parameter 'set_qrlzero'  was
    !         defined in the first part of 'eddy_diff.F90' module. 

    if( set_qrlzero ) then
        qrlw(:,:) = 0._r8
    else
        qrlw(:ncol,:pver) = qrlin(:ncol,:pver)
    endif

    ! Define effective stratus fraction using the grid-mean ql.
    ! Modification : The contribution of ice should be carefully considered.
    !                This should be done in combination with the 'qrlw' and
    !                overlapping assumption of liquid and ice stratus. 

    do k = 1, pver
       do i = 1, ncol
          if( choice_evhc .eq. 'ramp' .or. choice_radf .eq. 'ramp' ) then 
              cldeff(i,k) = cld(i,k) * min( ql(i,k) / qmin, 1._r8 )
          else
              cldeff(i,k) = cld(i,k)
          endif
       end do
    end do

    ! For an extended stability function in the stable regime, re-define
    ! alph4exe and ghmin. This is for future work.

    if( ricrit .eq. 0.19_r8 ) then
        alph4exs = alph4
        ghmin    = -3.5334_r8
    elseif( ricrit .gt. 0.19_r8 ) then
        alph4exs = -2._r8 * b1 * alph2 / ( alph3 - 2._r8 * b1 * alph5 ) / ricrit
        ghmin    = -1.e10_r8
    else
       errstring = 'ricrit should be larger than 0.19 in UW PBL'
       return
    endif

    !
    ! Initialization of Diagnostic Output
    !

    do i = 1, ncol
       went(i)                  = 0._r8
       wet_CL(i,:ncvmax)        = 0._r8
       web_CL(i,:ncvmax)        = 0._r8
       jtbu_CL(i,:ncvmax)       = 0._r8
       jbbu_CL(i,:ncvmax)       = 0._r8
       evhc_CL(i,:ncvmax)       = 0._r8
       jt2slv_CL(i,:ncvmax)     = 0._r8
       n2ht_CL(i,:ncvmax)       = 0._r8
       n2hb_CL(i,:ncvmax)       = 0._r8                    
       lwp_CL(i,:ncvmax)        = 0._r8
       opt_depth_CL(i,:ncvmax)  = 0._r8
       radinvfrac_CL(i,:ncvmax) = 0._r8
       radf_CL(i,:ncvmax)       = 0._r8
       wstar_CL(i,:ncvmax)      = 0._r8          
       wstar3fact_CL(i,:ncvmax) = 0._r8
       ricl(i,:ncvmax)          = 0._r8
       ghcl(i,:ncvmax)          = 0._r8
       shcl(i,:ncvmax)          = 0._r8
       smcl(i,:ncvmax)          = 0._r8
       ebrk(i,:ncvmax)          = 0._r8
       wbrk(i,:ncvmax)          = 0._r8
       lbrk(i,:ncvmax)          = 0._r8
       gh_a(i,:pver+1)          = 0._r8
       sh_a(i,:pver+1)          = 0._r8
       sm_a(i,:pver+1)          = 0._r8
       ri_a(i,:pver+1)          = 0._r8
       ipbl(i)                  = 0
       kpblh(i)                 = pver
       wsed_CL(i,:ncvmax)       = 0._r8          
    end do  

    ! kvh and kvm are stored over timesteps in 'vertical_diffusion.F90' and 
    ! passed in as kvh_in and kvm_in.  However,  at the first timestep they
    ! need to be computed and these are done just before calling 'caleddy'.   
    ! kvm and kvh are also stored over iterative time step in the first part
    ! of 'eddy_diff.F90'

    ! Initialize kvh and kvm to kvf
    kvh(:,:) = kvf(:,:)
    kvm(:,:) = kvf(:,:)
    ! Zero diagnostic quantities for the new diffusion step.
    wcap(:,:) = 0._r8
    leng(:,:) = 0._r8
    tke(:,:)  = 0._r8
    turbtype(:,:) = 0


    ! Initialize 'bprod' [ m2/s3 ] and 'sprod' [ m2/s3 ] at all interfaces.
    ! Note this initialization is a hybrid initialization since 'n2' [s-2] and 's2' [s-2]
    ! are calculated from the given current initial profile, while 'kvh_in' [m2/s] and 
    ! 'kvm_in' [m2/s] are from the previous iteration or previous time step.
    ! This initially guessed 'bprod' and 'sprod' will be updated at the end of this 
    ! 'caleddy' subroutine for diagnostic output.
    ! This computation of 'brpod,sprod' below is necessary for wstar-based entrainment closure.

    do k = 2, pver
       do i = 1, ncol
            bprod(i,k) = -kvh_in(i,k) * n2(i,k)
            sprod(i,k) =  kvm_in(i,k) * s2(i,k)
       end do
    end do

    ! Set 'bprod' and 'sprod' at top and bottom interface.
    ! In calculating 'surface' (actually lowest half-layer) buoyancy flux,
    ! 'chu' at surface is defined to be the same as 'chu' at the mid-point
    ! of lowest model layer (pver) at the end of 'trbind'. The same is for
    ! the other buoyancy coefficients.  'sprod(i,pver+1)'  is defined in a
    ! consistent way as the definition of 'tkes' in the original code.
    ! ( Important Option ) If I want to isolate surface buoyancy flux from
    ! the other parts of CL regimes energetically even though bflxs > 0,
    ! all I should do is to re-define 'bprod(i,pver+1)=0' in the below 'do'
    ! block. Additionally for merging test of extending SBCL based on 'l2n2'
    ! in 'zisocl', I should use 'l2n2 = - wint / sh'  for similar treatment
    ! as previous code. All other parts of the code  are fully consistently
    ! treated by these change only.
    ! My future general convection scheme will use bflxs(i).

    do i = 1, ncol
       bprod(i,1) = 0._r8 ! Top interface
       sprod(i,1) = 0._r8 ! Top interface
       ch = chu(i,pver+1) * ( 1._r8 - sflh(i,pver) ) + chs(i,pver+1) * sflh(i,pver)   
       cm = cmu(i,pver+1) * ( 1._r8 - sflh(i,pver) ) + cms(i,pver+1) * sflh(i,pver)   
       bflxs(i) = ch * shflx(i) * rrho(i) + cm * qflx(i) * rrho(i)
       if( choice_tkes .eq. 'ibprod' ) then
           bprod(i,pver+1) = bflxs(i)
       else
           bprod(i,pver+1) = 0._r8
       endif
       sprod(i,pver+1) = (ustar(i)**3)/(vk*z(i,pver))
    end do

    ! Initially identify CL regimes in 'exacol'
    !    ktop  : Interface index of the CL top  external interface
    !    kbase : Interface index of the CL base external interface
    !    ncvfin: Number of total CLs
    ! Note that if surface buoyancy flux is positive ( bflxs = bprod(i,pver+1) > 0 ),
    ! surface interface is identified as an internal interface of CL. However, even
    ! though bflxs <= 0, if 'pver' interface is a CL internal interface (ri(pver)<0),
    ! surface interface is identified as an external interface of CL. If bflxs =< 0 
    ! and ri(pver) >= 0, then surface interface is identified as a stable turbulent
    ! intereface (STL) as shown at the end of 'caleddy'. Even though a 'minpblh' is
    ! passed into 'exacol', it is not used in the 'exacol'.

    call exacol( pcols, pver, ncol, ri, bflxs, minpblh, zi, ktop, kbase, ncvfin )

    ! Diagnostic output of CL interface indices before performing 'extending-merging'
    ! of CL regimes in 'zisocl'
    do i = 1, ncol
    do k = 1, ncvmax
       kbase_o(i,k) = real(kbase(i,k),r8)
       ktop_o(i,k)  = real(ktop(i,k),r8) 
       ncvfin_o(i)  = real(ncvfin(i),r8)
    end do
    end do 

    ! ----------------------------------- !
    ! Perform calculation for each column !
    ! ----------------------------------- !

    do i = 1, ncol

       ! Define Surface Interfacial Layer TKE, 'tkes'.
       ! In the current code, 'tkes' is used as representing TKE of surface interfacial
       ! layer (low half-layer of surface-based grid layer). In the code, when bflxs>0,
       ! surface interfacial layer is assumed to be energetically  coupled to the other
       ! parts of the CL regime based at the surface. In this sense, it is conceptually
       ! more reasonable to include both 'bprod' and 'sprod' in the definition of 'tkes'.
       ! Since 'tkes' cannot be negative, it is lower bounded by small positive number. 
       ! Note that inclusion of 'bprod' in the definition of 'tkes' may increase 'ebrk'
       ! and 'wstar3', and eventually, 'wet' at the CL top, especially when 'bflxs>0'.
       ! This might help to solve the problem of too shallow PBLH over the overcast Sc
       ! regime. If I want to exclude 'bprod(i,pver+1)' in calculating 'tkes' even when
       ! bflxs > 0, all I should to do is to set 'bprod(i,pver+1) = 0' in the above 
       ! initialization 'do' loop (explained above), NOT changing the formulation of
       ! tkes(i) in the below block. This is because for consistent treatment in the 
       ! other parts of the code also.
  
     ! tkes(i) = (b1*vk*z(i,pver)*sprod(i,pver+1))**(2._r8/3._r8)
       tkes(i) = max(b1*vk*z(i,pver)*(bprod(i,pver+1)+sprod(i,pver+1)), 1.e-7_r8)**(2._r8/3._r8)
       tkes(i) = min(tkes(i), tkemax)
       tke(i,pver+1)  = tkes(i)
       wcap(i,pver+1) = tkes(i)/b1

       ! Extend and merge the initially identified CLs, relabel the CLs, and calculate
       ! CL internal mean energetics and stability functions in 'zisocl'. 
       ! The CL nearest to the surface is CL(1) and the CL index, ncv, increases 
       ! with height. The following outputs are from 'zisocl'. Here, the dimension
       ! of below outputs are (pcols,ncvmax) (except the 'ncvfin(pcols)' and 
       ! 'belongcv(pcols,pver+1)) and 'ncv' goes from 1 to 'ncvfin'. 
       ! For 'ncv = ncvfin+1, ncvmax', below output are already initialized to be zero. 
       !      ncvfin       : Total number of CLs
       !      kbase(ncv)   : Base external interface index of CL
       !      ktop         : Top  external interface index of CL
       !      belongcv     : True if the interface (either internal or external) is CL  
       !      ricl         : Mean Richardson number of internal CL
       !      ghcl         : Normalized buoyancy production '-l2n2/2e' [no unit] of internal CL
       !      shcl         : Galperin instability function of heat-moisture of internal CL
       !      smcl         : Galperin instability function of momentum of internal CL
       !      lbrk, <l>int : Thickness of (energetically) internal CL (lint, [m])
       !      wbrk, <W>int : Mean normalized TKE of internal CL  ([m2/s2])
       !      ebrk, <e>int : Mean TKE of internal CL (b1*wbrk,[m2/s2])
       ! The ncvsurf is an identifier saying which CL regime is based at the surface.
       ! If 'ncvsurf=1', then the first CL regime is based at the surface. If surface
       ! interface is not a part of CL (neither internal nor external), 'ncvsurf = 0'.
       ! After identifying and including SRCLs into the normal CL regimes (where newly
       ! identified SRCLs are simply appended to the normal CL regimes using regime 
       ! indices of 'ncvfin+1','ncvfin+2' (as will be shown in the below SRCL part),..
       ! where 'ncvfin' is the final CL regime index produced after extending-merging 
       ! in 'zisocl' but before adding SRCLs), if any newly identified SRCL (e.g., 
       ! 'ncvfin+1') is based at surface, then 'ncvsurf = ncvfin+1'. Thus 'ncvsurf' can
       ! be 0, 1, or >1. 'ncvsurf' can be a useful diagnostic output.   

       ncvsurf = 0
       if( ncvfin(i) .gt. 0 ) then 
           call zisocl( pcols  , pver     , i        ,           &
                        z      , zi       , n2       , s2      , & 
                        bprod  , sprod    , bflxs    , tkes    , &
                        ncvfin , kbase    , ktop     , belongcv, &
                        ricl   , ghcl     , shcl     , smcl    , & 
                        lbrk   , wbrk     , ebrk     ,           & 
                        extend , extend_up, extend_dn, &
                        errstring)
           if (trim(errstring) /= "") return
           if( kbase(i,1) .eq. pver + 1 ) ncvsurf = 1
       else
           belongcv(i,:) = .false.
       endif

       ! Diagnostic output after finishing extending-merging process in 'zisocl'
       ! Since we are adding SRCL additionally, we need to print out these here.

       do k = 1, ncvmax
          kbase_mg(i,k) = real(kbase(i,k))
          ktop_mg(i,k)  = real(ktop(i,k)) 
          ncvfin_mg(i)  = real(ncvfin(i))
       end do 

       ! ----------------------- !
       ! Identification of SRCLs !
       ! ----------------------- !

     ! Modification : This cannot identify the 'cirrus' layer due to the condition of
     !                ql(i,k) .gt. qmin. This should be modified in future to identify
     !                a single thin cirrus layer.  
     !                Instead of ql, we may use cldn in future, including ice 
     !                contribution.

       ! ------------------------------------------------------------------------------ !
       ! Find single-layer radiatively-driven cloud-topped convective layers (SRCLs).   !
       ! SRCLs extend through a single model layer k, with entrainment at the top and   !
       ! bottom interfaces, unless bottom interface is the surface.                     !
       ! The conditions for an SRCL is identified are:                                  ! 
       !                                                                                !
       !   1. Cloud in the layer, k : ql(i,k) .gt. qmin = 1.e-5 [ kg/kg ]               !
       !   2. No cloud in the above layer (else assuming that some fraction of the LW   !
       !      flux divergence in layer k is concentrated at just below top interface    !
       !      of layer k is invalid). Then, this condition might be sensitive to the    !
       !      vertical resolution of grid.                                              !
       !   3. LW radiative cooling (SW heating is assumed uniformly distributed through !
       !      layer k, so not relevant to buoyancy production) in the layer k. However, !
       !      SW production might also contribute, which may be considered in a future. !
       !   4. Internal stratification 'n2ht' of upper-half layer should be unstable.    !
       !      The 'n2ht' is pure internal stratification of upper half layer, obtained  !
       !      using internal slopes of sl, qt in layer k (in contrast to conventional   !
       !      interfacial slope) and saturation fraction in the upper-half layer,       !
       !      sfuh(k) (in contrast to sfi(k)).                                          !
       !   5. Top and bottom interfaces not both in the same existing convective layer. !
       !      If SRCL is within the previouisly identified CL regimes, we don't define  !
       !      a new SRCL.                                                               !
       !   6. k >= ntop_turb + 1 = 2                                                    !
       !   7. Ri at the top interface > ricrit = 0.19 (otherwise turbulent mixing will  !
       !      broadly distribute the cloud top in the vertical, preventing localized    !
       !      radiative destabilization at the top interface).                          !
       !                                                                                !
       ! Note if 'k = pver', it identifies a surface-based single fog layer, possibly,  !
       ! warm advection fog. Note also the CL regime index of SRCLs itself increases    !
       ! with height similar to the regular CLs indices identified from 'zisocl'.       !
       ! ------------------------------------------------------------------------------ !

       ncv  = 1
       ncvf = ncvfin(i)

       if( choice_SRCL .eq. 'remove' ) goto 222 

       do k = nbot_turb, ntop_turb + 1, -1 ! 'k = pver, 2, -1' is a layer index.

          if( ql(i,k) .gt. qmin .and. ql(i,k-1) .lt. qmin .and. qrlw(i,k) .lt. 0._r8 &
                                .and. ri(i,k) .ge. ricrit ) then

              ! In order to avoid any confliction with the treatment of ambiguous layer,
              ! I need to impose an additional constraint that ambiguous layer cannot be
              ! SRCL. So, I added constraint that 'k+1' interface (base interface of k
              ! layer) should not be a part of previously identified CL. Since 'belongcv'
              ! is even true for external entrainment interfaces, below constraint is
              ! fully sufficient.
 
              if( choice_SRCL .eq. 'nonamb' .and. belongcv(i,k+1) ) then
                  go to 220 
              endif

              ch = ( 1._r8 - sfuh(i,k) ) * chu(i,k) + sfuh(i,k) * chs(i,k)
              cm = ( 1._r8 - sfuh(i,k) ) * cmu(i,k) + sfuh(i,k) * cms(i,k)

              n2htSRCL = ch * slslope(i,k) + cm * qtslope(i,k)

              if( n2htSRCL .le. 0._r8 ) then

                  ! Test if bottom and top interfaces are part of the pre-existing CL. 
                  ! If not, find appropriate index for the new SRCL. Note that this
                  ! calculation makes use of 'ncv set' obtained from 'zisocl'. The 
                  ! 'in_CL' is a parameter testing whether the new SRCL is already 
                  ! within the pre-existing CLs (.true.) or not (.false.). 

                  in_CL = .false.

                  do while ( ncv .le. ncvf )
                     if( ktop(i,ncv) .le. k ) then
                        if( kbase(i,ncv) .gt. k ) then 
                            in_CL = .true.
                        endif
                        exit             ! Exit from 'do while' loop if SRCL is within the CLs.
                     else
                        ncv = ncv + 1    ! Go up one CL
                     end if
                  end do ! ncv

                  if( .not. in_CL ) then ! SRCL is not within the pre-existing CLs.

                     ! Identify a new SRCL and add it to the pre-existing CL regime group.

                     ncvfin(i)       =  ncvfin(i) + 1
                     ncvnew          =  ncvfin(i)
                     ktop(i,ncvnew)  =  k
                     kbase(i,ncvnew) =  k+1
                     belongcv(i,k)   = .true.
                     belongcv(i,k+1) = .true.

                     ! Calculate internal energy of SRCL. There is no internal energy if
                     ! SRCL is elevated from the surface. Also, we simply assume neutral 
                     ! stability function. Note that this assumption of neutral stability
                     ! does not influence numerical calculation- stability functions here
                     ! are just for diagnostic output. In general SRCLs other than a SRCL 
                     ! based at surface with bflxs <= 0, there is no other way but to use
                     ! neutral stability function.  However, in case of SRCL based at the
                     ! surface,  we can explicitly calculate non-zero stability functions            
                     ! in a consistent way.   Even though stability functions of SRCL are
                     ! just diagnostic outputs not influencing numerical calculations, it
                     ! would be informative to write out correct reasonable values rather
                     ! than simply assuming neutral stability. I am doing this right now.
                     ! Similar calculations were done for the SBCL and when surface inter
                     ! facial layer was merged by overlying CL in 'ziscol'.

                     if( k .lt. pver ) then

                         wbrk(i,ncvnew) = 0._r8
                         ebrk(i,ncvnew) = 0._r8
                         lbrk(i,ncvnew) = 0._r8
                         ghcl(i,ncvnew) = 0._r8
                         shcl(i,ncvnew) = 0._r8
                         smcl(i,ncvnew) = 0._r8
                         ricl(i,ncvnew) = 0._r8

                     else ! Surface-based fog

                         if( bflxs(i) .gt. 0._r8 ) then    ! Incorporate surface TKE into CL interior energy
                                                           ! It is likely that this case cannot exist  since
                                                           ! if surface buoyancy flux is positive,  it would
                                                           ! have been identified as SBCL in 'zisocl' ahead. 
                             ebrk(i,ncvnew) = tkes(i)
                             lbrk(i,ncvnew) = z(i,pver)
                             wbrk(i,ncvnew) = tkes(i) / b1    
        
                             write(temp_string,*) 'Major mistake in SRCL: bflxs > 0 for surface-based SRCL'
                             warnstring = trim(warnstring)//temp_string
                             write(temp_string,*) 'bflxs = ', bflxs(i), &
                                  'ncvfin_o = ', ncvfin_o(i), &
                                  'ncvfin_mg = ', ncvfin_mg(i)
                             warnstring = trim(warnstring)//temp_string
                             do ks = 1, ncvmax
                                write(temp_string,*) 'ncv =', ks, ' ', kbase_o(i,ks), &
                                     ktop_o(i,ks), kbase_mg(i,ks), ktop_mg(i,ks)
                                warnstring = trim(warnstring)//temp_string
                             end do
                             errstring = 'CALEDDY: Major mistake in SRCL: bflxs > 0 for surface-based SRCL'
                             return
                         else                              ! Don't incorporate surface interfacial TKE into CL interior energy

                             ebrk(i,ncvnew) = 0._r8
                             lbrk(i,ncvnew) = 0._r8
                             wbrk(i,ncvnew) = 0._r8

                         endif

                         ! Calculate stability functions (ghcl, shcl, smcl, ricl) explicitly
                         ! using an reverse procedure starting from tkes(i). Note that it is
                         ! possible to calculate stability functions even when bflxs < 0.
                         ! Previous code just assumed neutral stability functions. Note that
                         ! since alph5 = 0.7 > 0, alph3 = -35 < 0, the denominator of gh  is
                         ! always positive if bflxs > 0. However, if bflxs < 0,  denominator
                         ! can be zero. For this case, we provide a possible maximum negative
                         ! value (the most stable state) to gh. Note also tkes(i) is always a
                         ! positive value by a limiter. Also, sprod(i,pver+1) > 0 by limiter.
                         
                         gg = 0.5_r8 * vk * z(i,pver) * bprod(i,pver+1) / ( tkes(i)**(3._r8/2._r8) )
                         if( abs(alph5-gg*alph3) .le. 1.e-7_r8 ) then
                           ! gh = -0.28_r8
                           ! gh = -3.5334_r8
                             gh = ghmin
                         else    
                             gh = gg / ( alph5 - gg * alph3 )
                         end if 
                       ! gh = min(max(gh,-0.28_r8),0.0233_r8)
                       ! gh = min(max(gh,-3.5334_r8),0.0233_r8)
                         gh = min(max(gh,ghmin),0.0233_r8)
                         ghcl(i,ncvnew) =  gh
                         shcl(i,ncvnew) =  max(0._r8,alph5/(1._r8+alph3*gh))
                         smcl(i,ncvnew) =  max(0._r8,(alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
                         ricl(i,ncvnew) = -(smcl(i,ncvnew)/shcl(i,ncvnew))*(bprod(i,pver+1)/sprod(i,pver+1))

                       ! 'ncvsurf' is CL regime index based at the surface. If there is no
                       ! such regime, then 'ncvsurf = 0'.
    
                         ncvsurf = ncvnew

                      end if

                  end if

              end if

          end if

   220 continue    

       end do ! End of 'k' loop where 'k' is a grid layer index running from 'pver' to 2

   222 continue

       ! -------------------------------------------------------------------------- !
       ! Up to this point, we identified all kinds of CL regimes :                  !
       !   1. A SBCL. By construction, 'bflxs > 0' for SBCL.                        !
       !   2. Surface-based CL with multiple layers and 'bflxs =< 0'                !
       !   3. Surface-based CL with multiple layers and 'bflxs > 0'                 !
       !   4. Regular elevated CL with two entraining interfaces                    ! 
       !   5. SRCLs. If SRCL is based at surface, it will be bflxs < 0.             !
       ! '1-4' were identified from 'zisocl' while '5' were identified separately   !
       ! after performing 'zisocl'. CL regime index of '1-4' increases with height  !
       ! ( e.g., CL = 1 is the CL regime nearest to the surface ) while CL regime   !
       ! index of SRCL is simply appended after the final index of CL regimes from  !
       ! 'zisocl'. However, CL regime indices of SRCLs itself increases with height !
       ! when there are multiple SRCLs, similar to the regular CLs from 'zisocl'.   !
       ! -------------------------------------------------------------------------- !

       ! Diagnostic output of final CL regimes indices
       
       do k = 1, ncvmax
          kbase_f(i,k) = real(kbase(i,k))
          ktop_f(i,k)  = real(ktop(i,k)) 
          ncvfin_f(i)  = real(ncvfin(i))
       end do 

       ! --------------------------------------------------------------------- !
       ! Compute radf for each CL in column by calling subroutine compute_radf !
       ! --------------------------------------------------------------------- !
       call compute_radf( choice_radf, i, pcols, pver, ncvmax,  ncvfin, ktop, qmin,         &
                          ql, pi, qrlw, g, cldeff, zi, chs, lwp_CL(i,:), opt_depth_CL(i,:), &
                          radinvfrac_CL(i,:), radf_CL(i,:) )

       ! ---------------------------------------- !
       ! Perform do loop for individual CL regime !
       ! ---------------------------------------- ! -------------------------------- !
       ! For individual CLs, compute                                                 !
       !   1. Entrainment rates at the CL top and (if any) base interfaces using     !
       !      appropriate entrainment closure (current code use 'wstar' closure).    !
       !   2. Net CL mean (i.e., including entrainment contribution) TKE (ebrk)      !
       !      and normalized TKE (wbrk).                                             ! 
       !   3. TKE (tke) and normalized TKE (wcap) profiles at all CL interfaces.     !
       !   4. ( kvm, kvh ) profiles at all CL interfaces.                            !
       !   5. ( bprod, sprod ) profiles at all CL interfaces.                        !
       ! Also calculate                                                              !
       !   1. PBL height as the top external interface of surface-based CL, if any.  !
       !   2. Characteristic excesses of convective 'updraft velocity (wpert)',      !
       !      'temperature (tpert)', and 'moisture (qpert)' in the surface-based CL, !
       !      if any, for use in the separate convection scheme.                     ! 
       ! If there is no surface-based CL, 'PBL height' and 'convective excesses' are !
       ! calculated later from surface-based STL (Stable Turbulent Layer) properties.!
       ! --------------------------------------------------------------------------- !

       ktblw = 0
       do ncv = 1, ncvfin(i)

          kt = ktop(i,ncv)
          kb = kbase(i,ncv)

          lwp        = lwp_CL(i,ncv)
          opt_depth  = opt_depth_CL(i,ncv)
          radinvfrac = radinvfrac_CL(i,ncv)
          radf       = radf_CL(i, ncv)

          ! Check whether surface interface is energetically interior or not.
          if( kb .eq. (pver+1) .and. bflxs(i) .le. 0._r8 ) then
              lbulk = zi(i,kt) - z(i,pver)
          else
              lbulk = zi(i,kt) - zi(i,kb)
          end if
          lbulk = min( lbulk, lbulk_max )

          ! Calculate 'turbulent length scale (leng)' and 'normalized TKE (wcap)'
          ! at all CL interfaces except the surface.  Note that below 'wcap' at 
          ! external interfaces are not correct. However, it does not influence 
          ! numerical calculation and correct normalized TKE at the entraining 
          ! interfaces will be re-calculated at the end of this 'do ncv' loop. 

          do k = min(kb,pver), kt, -1 
             if( choice_tunl .eq. 'rampcl' ) then
               ! In order to treat the case of 'ricl(i,ncv) >> 0' of surface-based SRCL
               ! with 'bflxs(i) < 0._r8', I changed ricl(i,ncv) -> min(0._r8,ricl(i,ncv))
               ! in the below exponential. This is necessary to prevent the model crash
               ! by too large values (e.g., 700) of ricl(i,ncv)   
                 tunlramp = ctunl*tunl*(1._r8-(1._r8-1._r8/ctunl)*exp(min(0._r8,ricl(i,ncv))))
                 tunlramp = min(max(tunlramp,tunl),ctunl*tunl)
             elseif( choice_tunl .eq. 'rampsl' ) then
                 tunlramp = ctunl*tunl
               ! tunlramp = 0.765_r8
             else
                 tunlramp = tunl
             endif
             if( choice_leng .eq. 'origin' ) then
                 leng(i,k) = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
               ! leng(i,k) = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
             else
                 leng(i,k) = min( vk*zi(i,k), tunlramp*lbulk )              
             endif
             leng(i,k) = min(leng_max(k), leng(i,k))
             wcap(i,k) = (leng(i,k)**2) * (-shcl(i,ncv)*n2(i,k)+smcl(i,ncv)*s2(i,k))
          end do ! k

          ! Calculate basic cross-interface variables ( jump condition ) across the 
          ! base external interface of CL.

          if( kb .lt. pver+1 ) then 

              jbzm = z(i,kb-1) - z(i,kb)                                      ! Interfacial layer thickness [m]
              jbsl = sl(i,kb-1) - sl(i,kb)                                    ! Interfacial jump of 'sl' [J/kg]
              jbqt = qt(i,kb-1) - qt(i,kb)                                    ! Interfacial jump of 'qt' [kg/kg]
              jbbu = n2(i,kb) * jbzm                                          ! Interfacial buoyancy jump [m/s2]
                                                                              ! considering saturation ( > 0 )
              jbbu = max(jbbu,jbumin)                                         ! Set minimum buoyancy jump, jbumin = 1.e-3
              jbu  = u(i,kb-1) - u(i,kb)                                      ! Interfacial jump of 'u' [m/s]
              jbv  = v(i,kb-1) - v(i,kb)                                      ! Interfacial jump of 'v' [m/s]
              ch   = (1._r8 -sflh(i,kb-1))*chu(i,kb) + sflh(i,kb-1)*chs(i,kb) ! Buoyancy coefficient just above the base interface
              cm   = (1._r8 -sflh(i,kb-1))*cmu(i,kb) + sflh(i,kb-1)*cms(i,kb) ! Buoyancy coefficient just above the base interface
              n2hb = (ch*jbsl + cm*jbqt)/jbzm                                 ! Buoyancy frequency [s-2]
                                                                              ! just above the base interface
              vyb  = n2hb*jbzm/jbbu                                           ! Ratio of 'n2hb/n2' at 'kb' interface
              vub  = min(1._r8,(jbu**2+jbv**2)/(jbbu*jbzm) )                  ! Ratio of 's2/n2 = 1/Ri' at 'kb' interface

          else 

            ! Below setting is necessary for consistent treatment when 'kb' is at the surface.
              jbbu = 0._r8
              n2hb = 0._r8
              vyb  = 0._r8
              vub  = 0._r8
              web  = 0._r8

          end if

          ! Calculate basic cross-interface variables ( jump condition ) across the 
          ! top external interface of CL. The meanings of variables are similar to
          ! the ones at the base interface.

          jtzm = z(i,kt-1) - z(i,kt)
          jtsl = sl(i,kt-1) - sl(i,kt)
          jtqt = qt(i,kt-1) - qt(i,kt)
          jtbu = n2(i,kt)*jtzm                                      ! Note : 'jtbu' is guaranteed positive by definition of CL top.
          jtbu = max(jtbu,jbumin)                                   ! But threshold it anyway to be sure.
          jtu  = u(i,kt-1) - u(i,kt)
          jtv  = v(i,kt-1) - v(i,kt)
          ch   = (1._r8 -sfuh(i,kt))*chu(i,kt) + sfuh(i,kt)*chs(i,kt) 
          cm   = (1._r8 -sfuh(i,kt))*cmu(i,kt) + sfuh(i,kt)*cms(i,kt) 
          n2ht = (ch*jtsl + cm*jtqt)/jtzm                       
          vyt  = n2ht*jtzm/jtbu                                  
          vut  = min(1._r8,(jtu**2+jtv**2)/(jtbu*jtzm))             

          ! Evaporative enhancement factor of entrainment rate at the CL top interface, evhc. 
          ! We take the full inversion strength to be 'jt2slv = slv(i,kt-2)-slv(i,kt)' 
          ! where 'kt-1' is in the ambiguous layer. However, for a cloud-topped CL overlain
          ! by another CL, it is possible that 'slv(i,kt-2) < slv(i,kt)'. To avoid negative
          ! or excessive evhc, we lower-bound jt2slv and upper-bound evhc.  Note 'jtslv' is
          ! used only for calculating 'evhc' : when calculating entrainment rate,   we will
          ! use normal interfacial buoyancy jump across CL top interface.

          evhc   = 1._r8
          jt2slv = 0._r8

        ! Modification : I should check whether below 'jbumin' produces reasonable limiting value.   
        !                In addition, our current formulation does not consider ice contribution. 

          if( choice_evhc .eq. 'orig' ) then

              if( ql(i,kt) .gt. qmin .and. ql(i,kt-1) .lt. qmin ) then 
                  jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
                  jt2slv = max( jt2slv, jbumin*slv(i,kt-1)/g )
                  evhc   = 1._r8 + a2l * a3l * latvap * ql(i,kt) / jt2slv
                  evhc   = min( evhc, evhcmax )
              end if

          elseif( choice_evhc .eq. 'ramp' ) then

              jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
              jt2slv = max( jt2slv, jbumin*slv(i,kt-1)/g )
              evhc   = 1._r8 + max(cldeff(i,kt)-cldeff(i,kt-1),0._r8) * a2l * a3l * latvap * ql(i,kt) / jt2slv
              evhc   = min( evhc, evhcmax )

          elseif( choice_evhc .eq. 'maxi' ) then

              qleff  = max( ql(i,kt-1), ql(i,kt) ) 
              jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
              jt2slv = max( jt2slv, jbumin*slv(i,kt-1)/g )
              evhc   = 1._r8 + a2l * a3l * latvap * qleff / jt2slv
              evhc   = min( evhc, evhcmax )

          endif

          ! ------------------------------------------------------------------- !
          ! Calculate 'wstar3' by summing buoyancy productions within CL from   !
          !   1. Interior buoyancy production ( bprod: fcn of TKE )             !
          !   2. Cloud-top radiative cooling                                    !
          !   3. Surface buoyancy flux contribution only when bflxs > 0.        !
          !      Note that master length scale, lbulk, has already been         !
          !      corrctly defined at the first part of this 'do ncv' loop       !
          !      considering the sign of bflxs.                                 !
          ! This 'wstar3' is used for calculation of entrainment rate.          !
          ! Note that this 'wstar3' formula does not include shear production   !
          ! and the effect of drizzle, which should be included later.          !
          ! Q : Strictly speaking, in calculating interior buoyancy production, ! 
          !     the use of 'bprod' is not correct, since 'bprod' is not correct !
          !     value but initially guessed value.   More reasonably, we should ! 
          !     use '-leng(i,k)*sqrt(b1*wcap(i,k))*shcl(i,ncv)*n2(i,k)' instead !
          !     of 'bprod(i,k)', although this is still an  approximation since !
          !     tke(i,k) is not exactly 'b1*wcap(i,k)'  due to a transport term.! 
          !     However since iterative calculation will be performed after all,! 
          !     below might also be OK. But I should test this alternative.     !
          ! ------------------------------------------------------------------- !      

          dzht   = zi(i,kt)  - z(i,kt)     ! Thickness of CL top half-layer
          dzhb   = z(i,kb-1) - zi(i,kb)    ! Thickness of CL bot half-layer
          wstar3 = radf * dzht
          do k = kt + 1, kb - 1 ! If 'kt = kb - 1', this loop will not be performed. 
               wstar3 =  wstar3 + bprod(i,k) * ( z(i,k-1) - z(i,k) )
             ! Below is an alternative which may speed up convergence.
             ! However, for interfaces merged into original CL, it can
             ! be 'wcap(i,k)<0' since 'n2(i,k)>0'.  Thus, I should use
             ! the above original one.
             ! wstar3 =  wstar3 - leng(i,k)*sqrt(b1*wcap(i,k))*shcl(i,ncv)*n2(i,k)* &
             !                    (z(i,k-1) - z(i,k))
          end do      
          if( kb .eq. (pver+1) .and. bflxs(i) .gt. 0._r8 ) then
             wstar3 = wstar3 + bflxs(i) * dzhb
           ! wstar3 = wstar3 + bprod(i,pver+1) * dzhb
          end if   
          wstar3 = max( 2.5_r8 * wstar3, 0._r8 )
   
          ! -------------------------------------------------------------- !
          ! Below single block is for 'sedimentation-entrainment feedback' !
          ! -------------------------------------------------------------- !          

          if( id_sedfact ) then
            ! wsed    = 7.8e5_r8*(ql(i,kt)/ncliq(i,kt))**(2._r8/3._r8)
              sedfact = exp(-ased*wsedl(i,kt)/(wstar3**(1._r8/3._r8)+1.e-6_r8))
              wsed_CL(i,ncv) = wsedl(i,kt)
              if( choice_evhc .eq. 'orig' ) then
                  if (ql(i,kt).gt.qmin .and. ql(i,kt-1).lt.qmin) then
                      jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
                      jt2slv = max(jt2slv, jbumin*slv(i,kt-1)/g)
                      evhc = 1._r8+sedfact*a2l*a3l*latvap*ql(i,kt) / jt2slv
                      evhc = min(evhc,evhcmax)
                  end if
              elseif( choice_evhc .eq. 'ramp' ) then
                  jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
                  jt2slv = max(jt2slv, jbumin*slv(i,kt-1)/g)
                  evhc = 1._r8+max(cldeff(i,kt)-cldeff(i,kt-1),0._r8)*sedfact*a2l*a3l*latvap*ql(i,kt) / jt2slv
                  evhc = min(evhc,evhcmax)
              elseif( choice_evhc .eq. 'maxi' ) then
                  qleff  = max(ql(i,kt-1),ql(i,kt))
                  jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
                  jt2slv = max(jt2slv, jbumin*slv(i,kt-1)/g)
                  evhc = 1._r8+sedfact*a2l*a3l*latvap*qleff / jt2slv
                  evhc = min(evhc,evhcmax)
              endif
          endif

          ! -------------------------------------------------------------------------- !
          ! Now diagnose CL top and bottom entrainment rates (and the contribution of  !
          ! top/bottom entrainments to wstar3) using entrainment closures of the form  !
          !                                                                            !        
          !                   wet = cet*wstar3, web = ceb*wstar3                       !
          !                                                                            !
          ! where cet and ceb depend on the entrainment interface jumps, ql, etc.      !
          ! No entrainment is diagnosed unless the wstar3 > 0. Note '1/wstar3fact' is  !
          ! a factor indicating the enhancement of wstar3 due to entrainment process.  !
          ! Q : Below setting of 'wstar3fact = max(..,0.5)'might prevent the possible  !
          !     case when buoyancy consumption by entrainment is  stronger than cloud  !
          !     top radiative cooling production. Is that OK ? No.  According to bulk  !
          !     modeling study, entrainment buoyancy consumption was always a certain  !
          !     fraction of other net productions, rather than a separate sum.  Thus,  !
          !     below max limit of wstar3fact is correct.   'wstar3fact = max(.,0.5)'  !
          !     prevents unreasonable enhancement of CL entrainment rate by cloud-top  !
          !     entrainment instability, CTEI.                                         !
          ! Q : Use of the same dry entrainment coefficient, 'a1i' both at the CL  top !
          !     and base interfaces may result in too small 'wstar3' and 'ebrk' below, !
          !     as was seen in my generalized bulk modeling study. This should be re-  !
          !     considered later                                                       !
          ! -------------------------------------------------------------------------- !
          
          if( wstar3 .gt. 0._r8 ) then
              cet = a1i * evhc / ( jtbu * lbulk )
              if( kb .eq. pver + 1 ) then 
                  wstar3fact = max( 1._r8 + 2.5_r8 * cet * n2ht * jtzm * dzht, wstar3factcrit )
              else    
                  ceb = a1i / ( jbbu * lbulk )
                  wstar3fact = max( 1._r8 + 2.5_r8 * cet * n2ht * jtzm * dzht &
                                          + 2.5_r8 * ceb * n2hb * jbzm * dzhb, wstar3factcrit )
              end if
              wstar3 = wstar3 / wstar3fact       
          else ! wstar3 == 0
              wstar3fact = 0._r8 ! This is just for dianostic output
              cet        = 0._r8
              ceb        = 0._r8
          end if 

          ! ---------------------------------------------------------------------------- !
          ! Calculate net CL mean TKE including entrainment contribution by solving a    !
          ! canonical cubic equation. The solution of cubic equ. is 'rootp**2 = ebrk'    !
          ! where 'ebrk' originally (before solving cubic eq.) was interior CL mean TKE, !
          ! but after solving cubic equation,  it is replaced by net CL mean TKE in the  !
          ! same variable 'ebrk'.                                                        !
          ! ---------------------------------------------------------------------------- !
          ! Solve cubic equation (canonical form for analytic solution)                  !
          !   r^3 - 3*trmp*r - 2*trmq = 0,   r = sqrt<e>                                 ! 
          ! to estimate <e> for CL, derived from layer-mean TKE balance:                 !
          !                                                                              !
          !   <e>^(3/2)/(b_1*<l>) \approx <B + S>   (*)                                  !
          !   <B+S> = (<B+S>_int * l_int + <B+S>_et * dzt + <B+S>_eb * dzb)/lbulk        !
          !   <B+S>_int = <e>^(1/2)/(b_1*<l>)*<e>_int                                    !
          !   <B+S>_et  = (-vyt+vut)*wet*jtbu + radf                                     !
          !   <B+S>_eb  = (-vyb+vub)*web*jbbu                                            !
          !                                                                              !
          ! where:                                                                       !
          !   <> denotes a vertical avg (over the whole CL unless indicated)             !
          !   l_int (called lbrk below) is aggregate thickness of interior CL layers     !
          !   dzt = zi(i,kt)-z(i,kt)   is thickness of top entrainment layer             !
          !   dzb = z(i,kb-1)-zi(i,kb) is thickness of bot entrainment layer             !
          !   <e>_int (called ebrk below) is the CL-mean TKE if only interior            !
          !                               interfaces contributed.                        !
          !   wet, web                  are top. bottom entrainment rates                !
          !                                                                              !
          ! For a single-level radiatively-driven convective layer, there are no         ! 
          ! interior interfaces so 'ebrk' = 'lbrk' = 0. If the CL goes to the            !
          ! surface, 'vyb' and 'vub' are set to zero before and 'ebrk' and 'lbrk'        !
          ! have already incorporated the surface interfacial layer contribution,        !
          ! so the same formulas still apply.                                            !
          !                                                                              !
          ! In the original formulation based on TKE,                                    !
          !    wet*jtbu = a1l*evhc*<e>^3/2/leng(i,kt)                                    ! 
          !    web*jbbu = a1l*<e>^3/2/leng(i,kt)                                         !
          !                                                                              !
          ! In the wstar formulation                                                     !
          !    wet*jtbu = a1i*evhc*wstar3/lbulk                                          !
          !    web*jbbu = a1i*wstar3/lbulk,                                              !
          ! ---------------------------------------------------------------------------- !

          fact = ( evhc * ( -vyt + vut ) * dzht + ( -vyb + vub ) * dzhb * leng(i,kb) / leng(i,kt) ) / lbulk

          if( wstarent ) then

              ! (Option 1) 'wstar' entrainment formulation 
              ! Here trmq can have either sign, and will usually be nonzero even for non-
              ! cloud topped CLs.  If trmq > 0, there will be two positive roots r; we take 
              ! the larger one. Why ? If necessary, we limit entrainment and wstar to prevent
              ! a solution with r < ccrit*wstar ( Why ? ) where we take ccrit = 0.5. 

              trma = 1._r8          
              trmp = ebrk(i,ncv) * ( lbrk(i,ncv) / lbulk ) / 3._r8 + ntzero
              trmq = 0.5_r8 * b1 * ( leng(i,kt)  / lbulk ) * ( radf * dzht + a1i * fact * wstar3 )

              ! Check if there is an acceptable root with r > rcrit = ccrit*wstar. 
              ! To do this, first find local minimum fmin of the cubic f(r) at sqrt(p), 
              ! and value fcrit = f(rcrit).

              rmin  = sqrt(trmp)
              fmin  = rmin * ( rmin * rmin - 3._r8 * trmp ) - 2._r8 * trmq
              wstar = wstar3**onet
              rcrit = ccrit * wstar
              fcrit = rcrit * ( rcrit * rcrit - 3._r8 * trmp ) - 2._r8 * trmq

              ! No acceptable root exists (noroot = .true.) if either:
              !    1) rmin < rcrit (in which case cubic is monotone increasing for r > rcrit)
              !       and f(rcrit) > 0.
              ! or 2) rmin > rcrit (in which case min of f(r) in r > rcrit is at rmin)
              !       and f(rmin) > 0.  
              ! In this case, we reduce entrainment and wstar3 such that r/wstar = ccrit;
              ! this changes the coefficients of the cubic.   It might be informative to
              ! check when and how many 'noroot' cases occur,  since when 'noroot',   we
              ! will impose arbitrary limit on 'wstar3, wet, web, and ebrk' using ccrit.

              noroot = ( ( rmin .lt. rcrit ) .and. ( fcrit .gt. 0._r8 ) ) &
                  .or. ( ( rmin .ge. rcrit ) .and. ( fmin  .gt. 0._r8 ) )
              if( noroot ) then ! Solve cubic for r
                  trma = 1._r8 - b1 * ( leng(i,kt) / lbulk ) * a1i * fact / ccrit**3
                  trma = max( trma, 0.5_r8 )  ! Limit entrainment enhancement of ebrk
                  trmp = trmp / trma 
                  trmq = 0.5_r8 * b1 * ( leng(i,kt) / lbulk ) * radf * dzht / trma
              end if   ! noroot

              ! Solve the cubic equation

              qq = trmq**2 - trmp**3
              if( qq .ge. 0._r8 ) then 
                  rootp = ( trmq + sqrt(qq) )**(1._r8/3._r8) + ( max( trmq - sqrt(qq), 0._r8 ) )**(1._r8/3._r8)
              else
                  rootp = 2._r8 * sqrt(trmp) * cos( acos( trmq / sqrt(trmp**3) ) / 3._r8 )
              end if
 
              ! Adjust 'wstar3' only if there is 'noroot'. 
              ! And calculate entrainment rates at the top and base interfaces.

              if( noroot )  wstar3 = ( rootp / ccrit )**3     ! Adjust wstar3 
              wet = cet * wstar3                              ! Find entrainment rates
              if( kb .lt. pver + 1 ) web = ceb * wstar3       ! When 'kb.eq.pver+1', it was set to web=0. 

          else !

              ! (Option.2) wstarentr = .false. Use original entrainment formulation.
              ! trmp > 0 if there are interior interfaces in CL, trmp = 0 otherwise.
              ! trmq > 0 if there is cloudtop radiative cooling, trmq = 0 otherwise.
             
              trma = 1._r8 - b1 * a1l * fact
              trma = max( trma, 0.5_r8 )  ! Prevents runaway entrainment instability
              trmp = ebrk(i,ncv) * ( lbrk(i,ncv) / lbulk ) / ( 3._r8 * trma )
              trmq = 0.5_r8 * b1 * ( leng(i,kt)  / lbulk ) * radf * dzht / trma

              qq = trmq**2 - trmp**3
              if( qq .ge. 0._r8 ) then 
                  rootp = ( trmq + sqrt(qq) )**(1._r8/3._r8) + ( max( trmq - sqrt(qq), 0._r8 ) )**(1._r8/3._r8)
              else ! Also part of case 3
                  rootp = 2._r8 * sqrt(trmp) * cos( acos( trmq / sqrt(trmp**3) ) / 3._r8 )
              end if   ! qq

             ! Find entrainment rates and limit them by free-entrainment values a1l*sqrt(e)

              wet = a1l * rootp * min( evhc * rootp**2 / ( leng(i,kt) * jtbu ), 1._r8 )   
              if( kb .lt. pver + 1 ) web = a1l * rootp * min( evhc * rootp**2 / ( leng(i,kb) * jbbu ), 1._r8 )

          end if ! wstarentr

          ! ---------------------------------------------------- !
          ! Finally, get the net CL mean TKE and normalized TKE  ! 
          ! ---------------------------------------------------- !

          ebrk(i,ncv) = rootp**2
          ebrk(i,ncv) = min(ebrk(i,ncv),tkemax) ! Limit CL-avg TKE used for entrainment
          wbrk(i,ncv) = ebrk(i,ncv)/b1  
        
          ! The only way ebrk = 0 is for SRCL which are actually radiatively cooled 
          ! at top interface. In this case, we remove 'convective' label from the 
          ! interfaces around this layer. This case should now be impossible, so 
          ! we flag it. Q: I can't understand why this case is impossible now. Maybe,
          ! due to various limiting procedures used in solving cubic equation ? 
          ! In case of SRCL, 'ebrk' should be positive due to cloud top LW radiative
          ! cooling contribution, although 'ebrk(internal)' of SRCL before including
          ! entrainment contribution (which include LW cooling contribution also) is
          ! zero. 

          if( ebrk(i,ncv) .le. 0._r8 ) then
              write(temp_string,*) 'CALEDDY: Warning, CL with zero TKE, i, kt, kb ', i, kt, kb
              warnstring = trim(warnstring)//temp_string
              belongcv(i,kt) = .false.
              belongcv(i,kb) = .false. 
          end if
          
          ! ----------------------------------------------------------------------- !
          ! Calculate complete TKE profiles at all CL interfaces, capped by tkemax. !
          ! We approximate TKE = <e> at entrainment interfaces. However when CL is  !
          ! based at surface, correct 'tkes' will be inserted to tke(i,pver+1).     !
          ! Note that this approximation at CL external interfaces do not influence !
          ! numerical calculation since 'e' at external interfaces are not used  in !
          ! actual numerical calculation afterward. In addition in order to extract !
          ! correct TKE averaged over the PBL in the cumulus scheme,it is necessary !
          ! to set e = <e> at the top entrainment interface.  Since net CL mean TKE !
          ! 'ebrk' obtained by solving cubic equation already includes tkes  ( tkes !
          ! is included when bflxs > 0 but not when bflxs <= 0 into internal ebrk ),!
          ! 'tkes' should be written to tke(i,pver+1)                               !
          ! ----------------------------------------------------------------------- !

          ! 1. At internal interfaces          
          do k = kb - 1, kt + 1, -1
             rcap = ( b1 * ae + wcap(i,k) / wbrk(i,ncv) ) / ( b1 * ae + 1._r8 )
             rcap = min( max(rcap,rcapmin), rcapmax )
             tke(i,k) = ebrk(i,ncv) * rcap
             tke(i,k) = min( tke(i,k), tkemax )
             kvh(i,k) = leng(i,k) * sqrt(tke(i,k)) * shcl(i,ncv)
             kvm(i,k) = leng(i,k) * sqrt(tke(i,k)) * smcl(i,ncv)
             bprod(i,k) = -kvh(i,k) * n2(i,k)
             sprod(i,k) =  kvm(i,k) * s2(i,k)
             turbtype(i,k) = 2                     ! CL interior interfaces.
          end do

          ! 2. At CL top entrainment interface
          kentr = wet * jtzm
          kvh(i,kt) = kentr
          kvm(i,kt) = kentr
          bprod(i,kt) = -kentr * n2ht + radf       ! I must use 'n2ht' not 'n2'
          sprod(i,kt) =  kentr * s2(i,kt)
          turbtype(i,kt) = 4                       ! CL top entrainment interface
          trmp = -b1 * ae / ( 1._r8 + b1 * ae )
          trmq = -(bprod(i,kt)+sprod(i,kt))*b1*leng(i,kt)/(1._r8+b1*ae)/(ebrk(i,ncv)**(3._r8/2._r8))
          rcap = compute_cubic(0._r8,trmp,trmq)**2._r8
          rcap = min( max(rcap,rcapmin), rcapmax )
          tke(i,kt)  = ebrk(i,ncv) * rcap
          tke(i,kt)  = min( tke(i,kt), tkemax )

          ! 3. At CL base entrainment interface and double entraining interfaces
          ! When current CL base is also the top interface of CL regime below,
          ! simply add the two contributions for calculating eddy diffusivity
          ! and buoyancy/shear production. Below code correctly works because
          ! we (CL regime index) always go from surface upward.

          if( kb .lt. pver + 1 ) then 

              kentr = web * jbzm

              if( kb .ne. ktblw ) then

                  kvh(i,kb) = kentr
                  kvm(i,kb) = kentr
                  bprod(i,kb) = -kvh(i,kb)*n2hb     ! I must use 'n2hb' not 'n2'
                  sprod(i,kb) =  kvm(i,kb)*s2(i,kb)
                  turbtype(i,kb) = 3                ! CL base entrainment interface
                  trmp = -b1*ae/(1._r8+b1*ae)
                  trmq = -(bprod(i,kb)+sprod(i,kb))*b1*leng(i,kb)/(1._r8+b1*ae)/(ebrk(i,ncv)**(3._r8/2._r8))
                  rcap = compute_cubic(0._r8,trmp,trmq)**2._r8
                  rcap = min( max(rcap,rcapmin), rcapmax )
                  tke(i,kb)  = ebrk(i,ncv) * rcap
                  tke(i,kb)  = min( tke(i,kb),tkemax )

              else
                  
                  kvh(i,kb) = kvh(i,kb) + kentr 
                  kvm(i,kb) = kvm(i,kb) + kentr
                ! dzhb5 : Half thickness of the lowest  layer of  current CL regime
                ! dzht5 : Half thickness of the highest layer of adjacent CL regime just below current CL. 
                  dzhb5 = z(i,kb-1) - zi(i,kb)
                  dzht5 = zi(i,kb) - z(i,kb)
                  bprod(i,kb) = ( dzht5*bprod(i,kb) - dzhb5*kentr*n2hb )     / ( dzhb5 + dzht5 )
                  sprod(i,kb) = ( dzht5*sprod(i,kb) + dzhb5*kentr*s2(i,kb) ) / ( dzhb5 + dzht5 )
                  trmp = -b1*ae/(1._r8+b1*ae)
                  trmq = -kentr*(s2(i,kb)-n2hb)*b1*leng(i,kb)/(1._r8+b1*ae)/(ebrk(i,ncv)**(3._r8/2._r8))
                  rcap = compute_cubic(0._r8,trmp,trmq)**2._r8
                  rcap = min( max(rcap,rcapmin), rcapmax )
                  tke_imsi = ebrk(i,ncv) * rcap
                  tke_imsi = min( tke_imsi, tkemax )
                  tke(i,kb)  = ( dzht5*tke(i,kb) + dzhb5*tke_imsi ) / ( dzhb5 + dzht5 )               
                  tke(i,kb)  = min(tke(i,kb),tkemax)
                  turbtype(i,kb) = 5                ! CL double entraining interface      
                 
              end if

           else

             ! If CL base interface is surface, compute similarly using wcap(i,kb)=tkes/b1    
             ! Even when bflx < 0, use the same formula in order to impose consistency of
             ! tke(i,kb) at bflx = 0._r8
 
             rcap = (b1*ae + wcap(i,kb)/wbrk(i,ncv))/(b1*ae + 1._r8)
             rcap = min( max(rcap,rcapmin), rcapmax )
             tke(i,kb) = ebrk(i,ncv) * rcap
             tke(i,kb) = min( tke(i,kb),tkemax )

          end if

          ! Calculate wcap at all interfaces of CL. Put a  minimum threshold on TKE
          ! to prevent possible division by zero.  'wcap' at CL internal interfaces
          ! are already calculated in the first part of 'do ncv' loop correctly.
          ! When 'kb.eq.pver+1', below formula produces the identical result to the
          ! 'tkes(i)/b1' if leng(i,kb) is set to vk*z(i,pver). Note  wcap(i,pver+1)
          ! is already defined as 'tkes(i)/b1' at the first part of caleddy.
          
          wcap(i,kt) = (bprod(i,kt)+sprod(i,kt))*leng(i,kt)/sqrt(max(tke(i,kt),1.e-6_r8))
          if( kb .lt. pver + 1 ) then
              wcap(i,kb) = (bprod(i,kb)+sprod(i,kb))*leng(i,kb)/sqrt(max(tke(i,kb),1.e-6_r8))
          end if

          ! Save the index of upper external interface of current CL-regime in order to
          ! handle the case when this interface is also the lower external interface of 
          ! CL-regime located just above. 

          ktblw = kt 

          ! Diagnostic Output

          wet_CL(i,ncv)        = wet
          web_CL(i,ncv)        = web
          jtbu_CL(i,ncv)       = jtbu
          jbbu_CL(i,ncv)       = jbbu
          evhc_CL(i,ncv)       = evhc
          jt2slv_CL(i,ncv)     = jt2slv
          n2ht_CL(i,ncv)       = n2ht
          n2hb_CL(i,ncv)       = n2hb          
          wstar_CL(i,ncv)      = wstar          
          wstar3fact_CL(i,ncv) = wstar3fact          

       end do        ! ncv
 
       ! Calculate PBL height and characteristic cumulus excess for use in the
       ! cumulus convection shceme. Also define turbulence type at the surface
       ! when the lowest CL is based at the surface. These are just diagnostic
       ! outputs, not influencing numerical calculation of current PBL scheme.
       ! If the lowest CL is based at the surface, define the PBL depth as the
       ! CL top interface. The same rule is applied for all CLs including SRCL.

       if( ncvsurf .gt. 0 ) then

           ktopbl(i) = ktop(i,ncvsurf)
           pblh(i)   = zi(i, ktopbl(i))
           pblhp(i)  = pi(i, ktopbl(i))
           wpert(i)  = max(wfac*sqrt(ebrk(i,ncvsurf)),wpertmin)
           tpert(i)  = max(abs(shflx(i)*rrho(i)/cpair)*tfac/wpert(i),0._r8)
           qpert(i)  = max(abs(qflx(i)*rrho(i))*tfac/wpert(i),0._r8)

           if( bflxs(i) .gt. 0._r8 ) then
               turbtype(i,pver+1) = 2 ! CL interior interface
           else
               turbtype(i,pver+1) = 3 ! CL external base interface
           endif

           ipbl(i)  = 1
           kpblh(i) = max(ktopbl(i)-1, 1)
           went(i)  = wet_CL(i,ncvsurf)
       end if ! End of the calculationf of te properties of surface-based CL.

       ! -------------------------------------------- !
       ! Treatment of Stable Turbulent Regime ( STL ) !
       ! -------------------------------------------- !

       ! Identify top and bottom most (internal) interfaces of STL except surface.
       ! Also, calculate 'turbulent length scale (leng)' at each STL interfaces.     

       belongst(i,1) = .false.   ! k = 1 (top interface) is assumed non-turbulent
       do k = 2, pver            ! k is an interface index
          belongst(i,k) = ( ri(i,k) .lt. ricrit ) .and. ( .not. belongcv(i,k) )
          if( belongst(i,k) .and. ( .not. belongst(i,k-1) ) ) then
              kt = k             ! Top interface index of STL
          elseif( .not. belongst(i,k) .and. belongst(i,k-1) ) then
              kb = k - 1         ! Base interface index of STL
              lbulk = z(i,kt-1) - z(i,kb)
              lbulk = min( lbulk, lbulk_max )
              do ks = kt, kb
                 if( choice_tunl .eq. 'rampcl' ) then
                     tunlramp = tunl
                 elseif( choice_tunl .eq. 'rampsl' ) then
                    tunlramp = max( 1.e-3_r8, ctunl * tunl * exp(-log(ctunl)*ri(i,ks)/ricrit) )
                  ! tunlramp = 0.065_r8 + 0.7_r8 * exp(-20._r8*ri(i,ks))
                 else
                    tunlramp = tunl
                 endif
                 if( choice_leng .eq. 'origin' ) then
                     leng(i,ks) = ( (vk*zi(i,ks))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
                   ! leng(i,ks) = vk*zi(i,ks) / (1._r8+vk*zi(i,ks)/(tunlramp*lbulk))
                 else
                     leng(i,ks) = min( vk*zi(i,ks), tunlramp*lbulk )              
                 endif
                 leng(i,ks) = min(leng_max(ks), leng(i,ks))
              end do
          end if
       end do ! k

       ! Now look whether STL extends to ground.  If STL extends to surface,
       ! re-define master length scale,'lbulk' including surface interfacial
       ! layer thickness, and re-calculate turbulent length scale, 'leng' at
       ! all STL interfaces again. Note that surface interface is assumed to
       ! always be STL if it is not CL.   
       
       belongst(i,pver+1) = .not. belongcv(i,pver+1)

       if( belongst(i,pver+1) ) then     ! kb = pver+1 (surface  STL)

           turbtype(i,pver+1) = 1        ! Surface is STL interface
          
           if( belongst(i,pver) ) then   ! STL includes interior
             ! 'kt' already defined above as the top interface of STL
               lbulk = z(i,kt-1)          
           else                          ! STL with no interior turbulence
               kt = pver+1
               lbulk = z(i,kt-1)
           end if
           lbulk = min( lbulk, lbulk_max )

           ! PBL height : Layer mid-point just above the highest STL interface
           ! Note in contrast to the surface based CL regime where  PBL height
           ! was defined at the top external interface, PBL height of  surface
           ! based STL is defined as the layer mid-point.

           ktopbl(i) = kt - 1
           pblh(i)   = z(i,ktopbl(i))
           pblhp(i)  = 0.5_r8 * ( pi(i,ktopbl(i)) + pi(i,ktopbl(i)+1) )          

           ! Re-calculate turbulent length scale including surface interfacial
           ! layer contribution to lbulk.

           do ks = kt, pver
              if( choice_tunl .eq. 'rampcl' ) then
                  tunlramp = tunl
              elseif( choice_tunl .eq. 'rampsl' ) then
                  tunlramp = max(1.e-3_r8,ctunl*tunl*exp(-log(ctunl)*ri(i,ks)/ricrit))
                ! tunlramp = 0.065_r8 + 0.7_r8 * exp(-20._r8*ri(i,ks))
              else
                  tunlramp = tunl
              endif
              if( choice_leng .eq. 'origin' ) then
                  leng(i,ks) = ( (vk*zi(i,ks))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
                ! leng(i,ks) = vk*zi(i,ks) / (1._r8+vk*zi(i,ks)/(tunlramp*lbulk))
              else
                  leng(i,ks) = min( vk*zi(i,ks), tunlramp*lbulk )              
              endif
             leng(i,ks) = min(leng_max(ks), leng(i,ks))
           end do ! ks

           ! Characteristic cumulus excess of surface-based STL.
           ! We may be able to use ustar for wpert.

           wpert(i) = 0._r8 
           tpert(i) = max(shflx(i)*rrho(i)/cpair*fak/ustar(i),0._r8) ! CCM stable-layer forms
           qpert(i) = max(qflx(i)*rrho(i)*fak/ustar(i),0._r8)

           ipbl(i)  = 0
           kpblh(i) = ktopbl(i)

       end if

       ! Calculate stability functions and energetics at the STL interfaces
       ! except the surface. Note that tke(i,pver+1) and wcap(i,pver+1) are
       ! already calculated in the first part of 'caleddy', kvm(i,pver+1) &
       ! kvh(i,pver+1) were already initialized to be zero, bprod(i,pver+1)
       ! & sprod(i,pver+1) were direcly calculated from the bflxs and ustar.
       ! Note transport term is assumed to be negligible at STL interfaces.
           
       do k = 2, pver

          if( belongst(i,k) ) then

              turbtype(i,k) = 1    ! STL interfaces
              trma = alph3*alph4exs*ri(i,k) + 2._r8*b1*(alph2-alph4exs*alph5*ri(i,k))
              trmb = (alph3+alph4exs)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
              trmc = ri(i,k)
              det = max(trmb*trmb-4._r8*trma*trmc,0._r8)
              ! Sanity Check
              if( det .lt. 0._r8 ) then
                 errstring = 'The det < 0. for the STL in UW eddy_diff'
                 return
              end if                  
              gh = (-trmb + sqrt(det))/(2._r8*trma)
            ! gh = min(max(gh,-0.28_r8),0.0233_r8)
            ! gh = min(max(gh,-3.5334_r8),0.0233_r8)
              gh = min(max(gh,ghmin),0.0233_r8)
              sh = max(0._r8,alph5/(1._r8+alph3*gh))
              sm = max(0._r8,(alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))

              tke(i,k)   = b1*(leng(i,k)**2)*(-sh*n2(i,k)+sm*s2(i,k))
              tke(i,k)   = min(tke(i,k),tkemax)
              wcap(i,k)  = tke(i,k)/b1
              kvh(i,k)   = leng(i,k) * sqrt(tke(i,k)) * sh
              kvm(i,k)   = leng(i,k) * sqrt(tke(i,k)) * sm
              bprod(i,k) = -kvh(i,k) * n2(i,k)
              sprod(i,k) =  kvm(i,k) * s2(i,k)

          end if

       end do  ! k

       ! --------------------------------------------------- !
       ! End of treatment of Stable Turbulent Regime ( STL ) !
       ! --------------------------------------------------- !

       ! --------------------------------------------------------------- !
       ! Re-computation of eddy diffusivity at the entrainment interface !
       ! assuming that it is purely STL (0<Ri<0.19). Note even Ri>0.19,  !
       ! turbulent can exist at the entrainment interface since 'Sh,Sm'  !
       ! do not necessarily go to zero even when Ri>0.19. Since Ri can   !
       ! be fairly larger than 0.19 at the entrainment interface, I      !
       ! should set minimum value of 'tke' to be 0. in order to prevent  !
       ! sqrt(tke) from being imaginary.                                 !
       ! --------------------------------------------------------------- !

       ! goto 888

         do k = 2, pver

         if( ( turbtype(i,k) .eq. 3 ) .or. ( turbtype(i,k) .eq. 4 ) .or. &
             ( turbtype(i,k) .eq. 5 ) ) then

             trma = alph3*alph4exs*ri(i,k) + 2._r8*b1*(alph2-alph4exs*alph5*ri(i,k))
             trmb = (alph3+alph4exs)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
             trmc = ri(i,k)
             det  = max(trmb*trmb-4._r8*trma*trmc,0._r8)
             gh   = (-trmb + sqrt(det))/(2._r8*trma)
           ! gh   = min(max(gh,-0.28_r8),0.0233_r8)
           ! gh   = min(max(gh,-3.5334_r8),0.0233_r8)
             gh   = min(max(gh,ghmin),0.0233_r8)
             sh   = max(0._r8,alph5/(1._r8+alph3*gh))
             sm   = max(0._r8,(alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))

             lbulk = z(i,k-1) - z(i,k)
             lbulk = min( lbulk, lbulk_max )

             if( choice_tunl .eq. 'rampcl' ) then
                 tunlramp = tunl
             elseif( choice_tunl .eq. 'rampsl' ) then
                 tunlramp = max(1.e-3_r8,ctunl*tunl*exp(-log(ctunl)*ri(i,k)/ricrit))
               ! tunlramp = 0.065_r8 + 0.7_r8*exp(-20._r8*ri(i,k))
             else
                 tunlramp = tunl
             endif
             if( choice_leng .eq. 'origin' ) then
                 leng_imsi = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
               ! leng_imsi = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
             else
                 leng_imsi = min( vk*zi(i,k), tunlramp*lbulk )              
             endif
             leng_imsi = min(leng_max(k), leng_imsi)

             tke_imsi = b1*(leng_imsi**2)*(-sh*n2(i,k)+sm*s2(i,k))
             tke_imsi = min(max(tke_imsi,0._r8),tkemax)
             kvh_imsi = leng_imsi * sqrt(tke_imsi) * sh
             kvm_imsi = leng_imsi * sqrt(tke_imsi) * sm

             if( kvh(i,k) .lt. kvh_imsi ) then 
                 kvh(i,k)   =  kvh_imsi
                 kvm(i,k)   =  kvm_imsi
                 leng(i,k)  = leng_imsi
                 tke(i,k)   =  tke_imsi
                 wcap(i,k)  =  tke_imsi / b1
                 bprod(i,k) = -kvh_imsi * n2(i,k)
                 sprod(i,k) =  kvm_imsi * s2(i,k)
                 turbtype(i,k) = 1          ! This was added on Dec.10.2009 for use in microphysics.
             endif

         end if

         end do

 ! 888   continue 

       ! ------------------------------------------------------------------ !
       ! End of recomputation of eddy diffusivity at entrainment interfaces !
       ! ------------------------------------------------------------------ !

       ! As an option, we can impose a certain minimum back-ground diffusivity.

       ! do k = 1, pver+1
       !    kvh(i,k) = max(0.01_r8,kvh(i,k))
       !    kvm(i,k) = max(0.01_r8,kvm(i,k))
       ! enddo
 
       ! --------------------------------------------------------------------- !
       ! Diagnostic Output                                                     !
       ! Just for diagnostic purpose, calculate stability functions at  each   !
       ! interface including surface. Instead of assuming neutral stability,   !
       ! explicitly calculate stability functions using an reverse procedure   !
       ! starting from tkes(i) similar to the case of SRCL and SBCL in zisocl. !
       ! Note that it is possible to calculate stability functions even when   !
       ! bflxs < 0. Note that this inverse method allows us to define Ri even  !
       ! at the surface. Note also tkes(i) and sprod(i,pver+1) are always      !
       ! positive values by limiters (e.g., ustar_min = 0.01).                 !
       ! Dec.12.2006 : Also just for diagnostic output, re-set                 !
       ! 'bprod(i,pver+1)= bflxs(i)' here. Note that this setting does not     !
       ! influence numerical calculation at all - it is just for diagnostic    !
       ! output.                                                               !
       ! --------------------------------------------------------------------- !

       bprod(i,pver+1) = bflxs(i)
              
       gg = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
       if( abs(alph5-gg*alph3) .le. 1.e-7_r8 ) then
         ! gh = -0.28_r8
           if( bprod(i,pver+1) .gt. 0._r8 ) then
               gh = -3.5334_r8
           else
               gh = ghmin
           endif
       else    
           gh = gg/(alph5-gg*alph3)
       end if 

     ! gh = min(max(gh,-0.28_r8),0.0233_r8)
       if( bprod(i,pver+1) .gt. 0._r8 ) then
           gh = min(max(gh,-3.5334_r8),0.0233_r8)
       else
           gh = min(max(gh,ghmin),0.0233_r8)
       endif

       gh_a(i,pver+1) = gh     
       sh_a(i,pver+1) = max(0._r8,alph5/(1._r8+alph3*gh))
       if( bprod(i,pver+1) .gt. 0._r8 ) then       
           sm_a(i,pver+1) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh))
       else
           sm_a(i,pver+1) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
       endif
       ri_a(i,pver+1)  = -(sm_a(i,pver+1)/sh_a(i,pver+1))*(bprod(i,pver+1)/sprod(i,pver+1))

       do k = 1, pver
          if( ri(i,k) .lt. 0._r8 ) then
              trma = alph3*alph4*ri(i,k) + 2._r8*b1*(alph2-alph4*alph5*ri(i,k))
              trmb = (alph3+alph4)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
              trmc = ri(i,k)
              det  = max(trmb*trmb-4._r8*trma*trmc,0._r8)
              gh   = (-trmb + sqrt(det))/(2._r8*trma)
              gh   = min(max(gh,-3.5334_r8),0.0233_r8)
              gh_a(i,k) = gh
              sh_a(i,k) = max(0._r8,alph5/(1._r8+alph3*gh))
              sm_a(i,k) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh))
              ri_a(i,k) = ri(i,k)
          else
              if( ri(i,k) .gt. ricrit ) then
                  gh_a(i,k) = ghmin
                  sh_a(i,k) = 0._r8
                  sm_a(i,k) = 0._r8
                  ri_a(i,k) = ri(i,k)
              else
                  trma = alph3*alph4exs*ri(i,k) + 2._r8*b1*(alph2-alph4exs*alph5*ri(i,k))
                  trmb = (alph3+alph4exs)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
                  trmc = ri(i,k)
                  det  = max(trmb*trmb-4._r8*trma*trmc,0._r8)
                  gh   = (-trmb + sqrt(det))/(2._r8*trma)
                  gh   = min(max(gh,ghmin),0.0233_r8)
                  gh_a(i,k) = gh
                  sh_a(i,k) = max(0._r8,alph5/(1._r8+alph3*gh))
                  sm_a(i,k) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
                  ri_a(i,k) = ri(i,k)
              endif
          endif

       end do

    end do   ! End of column index loop, i 

    return

    end subroutine caleddy

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

    subroutine exacol( pcols, pver, ncol, ri, bflxs, minpblh, zi, ktop, kbase, ncvfin ) 

    ! ---------------------------------------------------------------------------- !
    ! Object : Find unstable CL regimes and determine the indices                  !
    !          kbase, ktop which delimit these unstable layers :                   !
    !          ri(kbase) > 0 and ri(ktop) > 0, but ri(k) < 0 for ktop < k < kbase. ! 
    ! Author : Chris  Bretherton 08/2000,                                          !
    !          Sungsu Park       08/2006, 11/2008                                  !
    !----------------------------------------------------------------------------- !

    implicit none

    ! --------------- !
    ! Input variables !
    ! --------------- !

    integer,  intent(in) :: pcols                  ! Number of atmospheric columns   
    integer,  intent(in) :: pver                   ! Number of atmospheric vertical layers   
    integer,  intent(in) :: ncol                   ! Number of atmospheric columns   

    real(r8), intent(in) :: ri(pcols,pver)         ! Moist gradient Richardson no.
    real(r8), intent(in) :: bflxs(pcols)           ! Buoyancy flux at surface
    real(r8), intent(in) :: minpblh(pcols)         ! Minimum PBL height based on surface stress
    real(r8), intent(in) :: zi(pcols,pver+1)       ! Interface heights

    ! ---------------- !
    ! Output variables !      
    ! ---------------- !

    integer, intent(out) :: kbase(pcols,ncvmax)    ! External interface index of CL base
    integer, intent(out) :: ktop(pcols,ncvmax)     ! External interface index of CL top
    integer, intent(out) :: ncvfin(pcols)          ! Total number of CLs

    ! --------------- !
    ! Local variables !
    ! --------------- !

    integer              :: i
    integer              :: k
    integer              :: ncv
    real(r8)             :: rimaxentr
    real(r8)             :: riex(pver+1)           ! Column Ri profile extended to surface

    ! ----------------------- !
    ! Main Computation Begins !
    ! ----------------------- !

    do i = 1, ncol
       ncvfin(i) = 0
       do ncv = 1, ncvmax
          ktop(i,ncv)  = 0
          kbase(i,ncv) = 0
       end do
    end do

    ! ------------------------------------------------------ !
    ! Find CL regimes starting from the surface going upward !
    ! ------------------------------------------------------ !
    
    rimaxentr = 0._r8   
    
    do i = 1, ncol

       riex(2:pver) = ri(i,2:pver)

       ! Below allows consistent treatment of surface and other interfaces.
       ! Simply, if surface buoyancy flux is positive, Ri of surface is set to be negative.

       riex(pver+1) = rimaxentr - bflxs(i) 

       ncv = 0
       k   = pver + 1 ! Work upward from surface interface

       do while ( k .gt. ntop_turb + 1 )

        ! Below means that if 'bflxs > 0' (do not contain '=' sign), surface
        ! interface is energetically interior surface. 
       
          if( riex(k) .lt. rimaxentr ) then 

              ! Identify a new CL

              ncv = ncv + 1

              ! First define 'kbase' as the first interface below the lower-most unstable interface
              ! Thus, Richardson number at 'kbase' is positive.

              kbase(i,ncv) = min(k+1,pver+1)

              ! Decrement k until top unstable level

              do while( riex(k) .lt. rimaxentr .and. k .gt. ntop_turb + 1 )
                 k = k - 1
              end do

              ! ktop is the first interface above upper-most unstable interface
              ! Thus, Richardson number at 'ktop' is positive. 

              ktop(i,ncv) = k
             
          else

              ! Search upward for a CL.

              k = k - 1

          end if

       end do ! End of CL regime finding for each atmospheric column

       ncvfin(i) = ncv    

    end do  ! End of atmospheric column do loop

    return 

    end subroutine exacol

    !============================================================================== !
    !                                                                               !
    !============================================================================== !
    
    subroutine zisocl( pcols  , pver  , long ,                                 & 
                       z      , zi    , n2   ,  s2      ,                      & 
                       bprod  , sprod , bflxs,  tkes    ,                      & 
                       ncvfin , kbase , ktop ,  belongcv,                      & 
                       ricl   , ghcl  , shcl ,  smcl    ,                      &
                       lbrk   , wbrk  , ebrk ,  extend  , extend_up, extend_dn,&
                       errstring)

    !------------------------------------------------------------------------ !
    ! Object : This 'zisocl' vertically extends original CLs identified from  !
    !          'exacol' using a merging test based on either 'wint' or 'l2n2' !
    !          and identify new CL regimes. Similar to the case of 'exacol',  !
    !          CL regime index increases with height.  After identifying new  !
    !          CL regimes ( kbase, ktop, ncvfin ),calculate CL internal mean  !
    !          energetics (lbrk : energetic thickness integral, wbrk, ebrk )  !
    !          and stability functions (ricl, ghcl, shcl, smcl) by including  !
    !          surface interfacial layer contribution when bflxs > 0.   Note  !
    !          that there are two options in the treatment of the energetics  !
    !          of surface interfacial layer (use_dw_surf= 'true' or 'false')  !
    ! Author : Sungsu Park 08/2006, 11/2008                                   !
    !------------------------------------------------------------------------ !

    implicit none

    ! --------------- !    
    ! Input variables !
    ! --------------- !

    integer,  intent(in)   :: long                    ! Longitude of the column
    integer,  intent(in)   :: pcols                   ! Number of atmospheric columns   
    integer,  intent(in)   :: pver                    ! Number of atmospheric vertical layers   
    real(r8), intent(in)   :: z(pcols, pver)          ! Layer mid-point height [ m ]
    real(r8), intent(in)   :: zi(pcols, pver+1)       ! Interface height [ m ]
    real(r8), intent(in)   :: n2(pcols, pver)         ! Buoyancy frequency at interfaces except surface [ s-2 ]
    real(r8), intent(in)   :: s2(pcols, pver)         ! Shear frequency at interfaces except surface [ s-2 ]
    real(r8), intent(in)   :: bprod(pcols,pver+1)     ! Buoyancy production [ m2/s3 ]. bprod(i,pver+1) = bflxs 
    real(r8), intent(in)   :: sprod(pcols,pver+1)     ! Shear production [ m2/s3 ]. sprod(i,pver+1) = usta**3/(vk*z(i,pver))
    real(r8), intent(in)   :: bflxs(pcols)            ! Surface buoyancy flux [ m2/s3 ]. bprod(i,pver+1) = bflxs 
    real(r8), intent(in)   :: tkes(pcols)             ! TKE at the surface [ s2/s2 ]

    ! ---------------------- !
    ! Input/output variables !
    ! ---------------------- !

    integer, intent(inout) :: kbase(pcols,ncvmax)     ! Base external interface index of CL
    integer, intent(inout) :: ktop(pcols,ncvmax)      ! Top external interface index of CL
    integer, intent(inout) :: ncvfin(pcols)           ! Total number of CLs

    ! ---------------- !
    ! Output variables !
    ! ---------------- !

    logical,  intent(out) :: belongcv(pcols,pver+1)   ! True if interface is in a CL ( either internal or external )
    real(r8), intent(out) :: ricl(pcols,ncvmax)       ! Mean Richardson number of internal CL
    real(r8), intent(out) :: ghcl(pcols,ncvmax)       ! Half of normalized buoyancy production of internal CL
    real(r8), intent(out) :: shcl(pcols,ncvmax)       ! Galperin instability function of heat-moisture of internal CL
    real(r8), intent(out) :: smcl(pcols,ncvmax)       ! Galperin instability function of momentum of internal CL
    real(r8), intent(out) :: lbrk(pcols,ncvmax)       ! Thickness of (energetically) internal CL ( lint, [m] )
    real(r8), intent(out) :: wbrk(pcols,ncvmax)       ! Mean normalized TKE of internal CL  [ m2/s2 ]
    real(r8), intent(out) :: ebrk(pcols,ncvmax)       ! Mean TKE of internal CL ( b1*wbrk, [m2/s2] )

    character(len=*), intent(out) :: errstring
    ! ------------------ !
    ! Internal variables !
    ! ------------------ !

    logical               :: extend                   ! True when CL is extended in zisocl
    logical               :: extend_up                ! True when CL is extended upward in zisocl
    logical               :: extend_dn                ! True when CL is extended downward in zisocl
    logical               :: bottom                   ! True when CL base is at surface ( kb = pver + 1 )

    integer               :: i                        ! Local index for the longitude
    integer               :: ncv                      ! CL Index increasing with height
    integer               :: incv
    integer               :: k
    integer               :: kb                       ! Local index for kbase
    integer               :: kt                       ! Local index for ktop
    integer               :: ncvinit                  ! Value of ncv at routine entrance 
    integer               :: cntu                     ! Number of merged CLs during upward   extension of individual CL
    integer               :: cntd                     ! Number of merged CLs during downward extension of individual CL
    integer               :: kbinc                    ! Index for incorporating underlying CL
    integer               :: ktinc                    ! Index for incorporating  overlying CL

    real(r8)              :: wint                     ! Normalized TKE of internal CL
    real(r8)              :: dwinc                    ! Normalized TKE of CL external interfaces
    real(r8)              :: dw_surf                  ! Normalized TKE of surface interfacial layer
    real(r8)              :: dzinc
    real(r8)              :: gh
    real(r8)              :: sh
    real(r8)              :: sm
    real(r8)              :: gh_surf                  ! Half of normalized buoyancy production in surface interfacial layer 
    real(r8)              :: sh_surf                  ! Galperin instability function in surface interfacial layer  
    real(r8)              :: sm_surf                  ! Galperin instability function in surface interfacial layer 
    real(r8)              :: l2n2                     ! Vertical integral of 'l^2N^2' over CL. Include thickness product
    real(r8)              :: l2s2                     ! Vertical integral of 'l^2S^2' over CL. Include thickness product
    real(r8)              :: dl2n2                    ! Vertical integration of 'l^2*N^2' of CL external interfaces
    real(r8)              :: dl2s2                    ! Vertical integration of 'l^2*S^2' of CL external interfaces
    real(r8)              :: dl2n2_surf               ! 'dl2n2' defined in the surface interfacial layer
    real(r8)              :: dl2s2_surf               ! 'dl2s2' defined in the surface interfacial layer  
    real(r8)              :: lint                     ! Thickness of (energetically) internal CL
    real(r8)              :: dlint                    ! Interfacial layer thickness of CL external interfaces
    real(r8)              :: dlint_surf               ! Surface interfacial layer thickness 
    real(r8)              :: lbulk                    ! Master Length Scale : Whole CL thickness from top to base external interface
    real(r8)              :: lz                       ! Turbulent length scale
    real(r8)              :: ricll                    ! Mean Richardson number of internal CL 
    real(r8)              :: trma
    real(r8)              :: trmb
    real(r8)              :: trmc
    real(r8)              :: det
    real(r8)              :: zbot                     ! Height of CL base
    real(r8)              :: l2rat                    ! Square of ratio of actual to initial CL (not used)
    real(r8)              :: gg                       ! Intermediate variable used for calculating stability functions of SBCL
    real(r8)              :: tunlramp                 ! Ramping tunl

    ! ----------------------- !
    ! Main Computation Begins !
    ! ----------------------- ! 

    i = long

    ! Initialize main output variables
    
    do k = 1, ncvmax
       ricl(i,k) = 0._r8
       ghcl(i,k) = 0._r8
       shcl(i,k) = 0._r8
       smcl(i,k) = 0._r8
       lbrk(i,k) = 0._r8
       wbrk(i,k) = 0._r8
       ebrk(i,k) = 0._r8
    end do
    extend    = .false.
    extend_up = .false.
    extend_dn = .false.

    ! ----------------------------------------------------------- !
    ! Loop over each CL to see if any of them need to be extended !
    ! ----------------------------------------------------------- !

    ncv = 1

    do while( ncv .le. ncvfin(i) )

       ncvinit = ncv
       cntu    = 0
       cntd    = 0
       kb      = kbase(i,ncv) 
       kt      = ktop(i,ncv)
       
       ! ---------------------------------------------------------------------------- !
       ! Calculation of CL interior energetics including surface before extension     !
       ! ---------------------------------------------------------------------------- !
       ! Note that the contribution of interior interfaces (not surface) to 'wint' is !
       ! accounted by using '-sh*l2n2 + sm*l2s2' while the contribution of surface is !
       ! accounted by using 'dwsurf = tkes/b1' when bflxs > 0. This approach is fully !
       ! reasonable. Another possible alternative,  which seems to be also consistent !
       ! is to calculate 'dl2n2_surf'  and  'dl2s2_surf' of surface interfacial layer !
       ! separately, and this contribution is explicitly added by initializing 'l2n2' !
       ! 'l2s2' not by zero, but by 'dl2n2_surf' and 'ds2n2_surf' below.  At the same !
       ! time, 'dwsurf' should be excluded in 'wint' calculation below. The only diff.!
       ! between two approaches is that in case of the latter approach, contributions !
       ! of surface interfacial layer to the CL mean stability function (ri,gh,sh,sm) !
       ! are explicitly included while the first approach is not. In this sense,  the !
       ! second approach seems to be more conceptually consistent,   but currently, I !
       ! (Sungsu) will keep the first default approach. There is a switch             !
       ! 'use_dw_surf' at the first part of eddy_diff.F90 chosing one of              !
       ! these two options.                                                           !
       ! ---------------------------------------------------------------------------- !
       
       ! ------------------------------------------------------ !   
       ! Step 0: Calculate surface interfacial layer energetics !
       ! ------------------------------------------------------ !

       lbulk      = zi(i,kt) - zi(i,kb)
       lbulk      = min( lbulk, lbulk_max )
       dlint_surf = 0._r8
       dl2n2_surf = 0._r8
       dl2s2_surf = 0._r8
       dw_surf    = 0._r8
       if( kb .eq. pver+1 ) then

           if( bflxs(i) .gt. 0._r8 ) then

               ! Calculate stability functions of surface interfacial layer
               ! from the given 'bprod(i,pver+1)' and 'sprod(i,pver+1)' using
               ! inverse approach. Since alph5>0 and alph3<0, denominator of
               ! gg is always positive if bprod(i,pver+1)>0.               

               gg    = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
               gh    = gg/(alph5-gg*alph3)
             ! gh    = min(max(gh,-0.28_r8),0.0233_r8)
               gh    = min(max(gh,-3.5334_r8),0.0233_r8)
               sh    = alph5/(1._r8+alph3*gh)
               sm    = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
               ricll = min(-(sm/sh)*(bprod(i,pver+1)/sprod(i,pver+1)),ricrit)

               ! Calculate surface interfacial layer contribution to CL internal
               ! energetics. By construction, 'dw_surf = -dl2n2_surf + ds2n2_surf'
               ! is exactly satisfied, which corresponds to assuming turbulent
               ! length scale of surface interfacial layer = vk * z(i,pver). Note
               ! 'dl2n2_surf','dl2s2_surf','dw_surf' include thickness product.   

               dlint_surf = z(i,pver)
               dl2n2_surf = -vk*(z(i,pver)**2)*bprod(i,pver+1)/(sh*sqrt(tkes(i)))
               dl2s2_surf =  vk*(z(i,pver)**2)*sprod(i,pver+1)/(sm*sqrt(tkes(i)))
               dw_surf    = (tkes(i)/b1)*z(i,pver) 

           else

               ! Note that this case can happen when surface is an external 
               ! interface of CL.
               lbulk = zi(i,kt) - z(i,pver)
               lbulk = min( lbulk, lbulk_max )

           end if

       end if
           
       ! ------------------------------------------------------ !   
       ! Step 1: Include surface interfacial layer contribution !
       ! ------------------------------------------------------ !
       
       lint = dlint_surf
       l2n2 = dl2n2_surf
       l2s2 = dl2s2_surf          
       wint = dw_surf
       if( use_dw_surf ) then
           l2n2 = 0._r8
           l2s2 = 0._r8
       else
           wint = 0._r8
       end if    
       
       ! --------------------------------------------------------------------------------- !
       ! Step 2. Include the contribution of 'pure internal interfaces' other than surface !
       ! --------------------------------------------------------------------------------- ! 
       
       if( kt .lt. kb - 1 ) then ! The case of non-SBCL.
                              
           do k = kb - 1, kt + 1, -1       
              if( choice_tunl .eq. 'rampcl' ) then
                ! Modification : I simply used the average tunlramp between the two limits.
                  tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
              elseif( choice_tunl .eq. 'rampsl' ) then
                  tunlramp = ctunl*tunl
                ! tunlramp = 0.765_r8
              else
                  tunlramp = tunl
              endif
              if( choice_leng .eq. 'origin' ) then
                  lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
                ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
              else
                  lz = min( vk*zi(i,k), tunlramp*lbulk )              
              endif
              lz = min(leng_max(k), lz)
              dzinc = z(i,k-1) - z(i,k)
              l2n2  = l2n2 + lz*lz*n2(i,k)*dzinc
              l2s2  = l2s2 + lz*lz*s2(i,k)*dzinc
              lint  = lint + dzinc
           end do

           ! Calculate initial CL stability functions (gh,sh,sm) and net
           ! internal energy of CL including surface contribution if any. 

         ! Modification : It seems that below cannot be applied when ricrit > 0.19.
         !                May need future generalization.

           ricll = min(l2n2/max(l2s2,ntzero),ricrit) ! Mean Ri of internal CL
           trma  = alph3*alph4*ricll+2._r8*b1*(alph2-alph4*alph5*ricll)
           trmb  = ricll*(alph3+alph4)+2._r8*b1*(-alph5*ricll+alph1)
           trmc  = ricll
           det   = max(trmb*trmb-4._r8*trma*trmc,0._r8)
           gh    = (-trmb + sqrt(det))/2._r8/trma
         ! gh    = min(max(gh,-0.28_r8),0.0233_r8)
           gh    = min(max(gh,-3.5334_r8),0.0233_r8)
           sh    = alph5/(1._r8+alph3*gh)
           sm    = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
           wint  = wint - sh*l2n2 + sm*l2s2 

       else ! The case of SBCL
 
           ! If there is no pure internal interface, use only surface interfacial
           ! values. However, re-set surface interfacial values such  that it can
           ! be used in the merging tests (either based on 'wint' or 'l2n2')  and
           ! in such that surface interfacial energy is not double-counted.
           ! Note that regardless of the choise of 'use_dw_surf', below should be
           ! kept as it is below, for consistent merging test of extending SBCL. 
       
           lint = dlint_surf
           l2n2 = dl2n2_surf
           l2s2 = dl2s2_surf 
           wint = dw_surf

           ! Aug.29.2006 : Only for the purpose of merging test of extending SRCL
           ! based on 'l2n2', re-define 'l2n2' of surface interfacial layer using
           ! 'wint'. This part is designed for similar treatment of merging as in
           ! the original 'eddy_diff.F90' code,  where 'l2n2' of SBCL was defined
           ! as 'l2n2 = - wint / sh'. Note that below block is used only when (1)
           ! surface buoyancy production 'bprod(i,pver+1)' is NOT included in the
           ! calculation of surface TKE in the initialization of 'bprod(i,pver+1)'
           ! in the main subroutine ( even though bflxs > 0 ), and (2) to force 
           ! current scheme be similar to the previous scheme in the treatment of  
           ! extending-merging test of SBCL based on 'l2n2'. Otherwise below line
           ! must be commented out. Note at this stage, correct non-zero value of
           ! 'sh' has been already computed.      

           if( choice_tkes .eq. 'ebprod' ) then
               l2n2 = - wint / sh 
           endif
           
       endif
           
       ! Set consistent upper limits on 'l2n2' and 'l2s2'. Below limits are
       ! reasonable since l2n2 of CL interior interface is always negative.

       l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
       l2s2 =  min( l2s2, tkemax*lint/(b1*sm))
       
       ! Note that at this stage, ( gh, sh, sm )  are the values of surface
       ! interfacial layer if there is no pure internal interface, while if
       ! there is pure internal interface, ( gh, sh, sm ) are the values of
       ! pure CL interfaces or the values that include both the CL internal
       ! interfaces and surface interfaces, depending on the 'use_dw_surf'.       
       
       ! ----------------------------------------------------------------------- !
       ! Perform vertical extension-merging process                              !
       ! ----------------------------------------------------------------------- !
       ! During the merging process, we assumed ( lbulk, sh, sm ) of CL external !
       ! interfaces are the same as the ones of the original merging CL. This is !
       ! an inevitable approximation since we don't know  ( sh, sm ) of external !
       ! interfaces at this stage.     Note that current default merging test is !
       ! purely based on buoyancy production without including shear production, !
       ! since we used 'l2n2' instead of 'wint' as a merging parameter. However, !
       ! merging test based on 'wint' maybe conceptually more attractable.       !
       ! Downward CL merging process is identical to the upward merging process, !
       ! but when the base of extended CL reaches to the surface, surface inter  !
       ! facial layer contribution to the energetic of extended CL must be done  !
       ! carefully depending on the sign of surface buoyancy flux. The contribu  !
       ! tion of surface interfacial layer energetic is included to the internal !
       ! energetics of merging CL only when bflxs > 0.                           !
       ! ----------------------------------------------------------------------- !
       
       ! ---------------------------- !
       ! Step 1. Extend the CL upward !
       ! ---------------------------- !
       
       extend = .false.    ! This will become .true. if CL top or base is extended

       ! Calculate contribution of potentially incorporable CL top interface

       if( choice_tunl .eq. 'rampcl' ) then
           tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
       elseif( choice_tunl .eq. 'rampsl' ) then
           tunlramp = ctunl*tunl
         ! tunlramp = 0.765_r8
       else
           tunlramp = tunl
       endif
       if( choice_leng .eq. 'origin' ) then
           lz = ( (vk*zi(i,kt))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
         ! lz = vk*zi(i,kt) / (1._r8+vk*zi(i,kt)/(tunlramp*lbulk))
       else
           lz = min( vk*zi(i,kt), tunlramp*lbulk )              
       endif
       lz = min(leng_max(kt), lz)

       dzinc = z(i,kt-1)-z(i,kt)
       dl2n2 = lz*lz*n2(i,kt)*dzinc
       dl2s2 = lz*lz*s2(i,kt)*dzinc
       dwinc = -sh*dl2n2 + sm*dl2s2

       ! ------------ !
       ! Merging Test !
       ! ------------ !

       ! The part of the below test that involves kt and z has different
       ! effects based on the model top.
       ! If the model top is in the stratosphere, we want the loop to
       ! continue until it either completes normally, or kt is pushed to
       ! the top of the model. The latter case should not happen, so this
       ! causes an error.
       ! If the model top is higher, as in WACCM and WACCM-X, if kt is
       ! pushed close to the model top, this may not represent an error at
       ! all, because of very different and more variable temperature/wind
       ! profiles at the model top. Therefore we simply exit the loop early
       ! and continue with no errors.

     ! do while (  dwinc .gt. ( rinc*dzinc*wint/(lint+(1._r8-rinc)*dzinc)) )  ! Merging test based on wint
     ! do while ( -dl2n2 .gt. (-rinc*dzinc*l2n2/(lint+(1._r8-rinc)*dzinc)) )  ! Merging test based on l2n2 
       do while ( -dl2n2 .gt. (-rinc*l2n2/(1._r8-rinc)) &                     ! Integral merging test
            .and. (kt > ntop_turb+2 .or. z(i,kt) < 50000._r8) )

          ! Add contribution of top external interface to interior energy.
          ! Note even when we chose 'use_dw_surf='true.', the contribution
          ! of surface interfacial layer to 'l2n2' and 'l2s2' are included
          ! here. However it is not double counting of surface interfacial
          ! energy : surface interfacial layer energy is counted in 'wint'
          ! formula and 'l2n2' is just used for performing merging test in
          ! this 'do while' loop.     

          lint = lint + dzinc
          l2n2 = l2n2 + dl2n2
          l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
          l2s2 = l2s2 + dl2s2
          wint = wint + dwinc

          ! Extend top external interface of CL upward after merging

          kt        = kt - 1
          extend    = .true.
          extend_up = .true.
          if( kt .eq. ntop_turb ) then
             errstring = 'zisocl: Error: Tried to extend CL to the model top'
             return
          end if

          ! If the top external interface of extending CL is the same as the 
          ! top interior interface of the overlying CL, overlying CL will be
          ! automatically merged. Then,reduce total number of CL regime by 1. 
          ! and increase 'cntu'(number of merged CLs during upward extension)
          ! by 1.
 
          ktinc = kbase(i,ncv+cntu+1) - 1  ! Lowest interior interface of overlying CL

          if( kt .eq. ktinc ) then

              do k = kbase(i,ncv+cntu+1) - 1, ktop(i,ncv+cntu+1) + 1, -1

                 if( choice_tunl .eq. 'rampcl' ) then
                     tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
                 elseif( choice_tunl .eq. 'rampsl' ) then
                     tunlramp = ctunl*tunl
                   ! tunlramp = 0.765_r8
                 else
                     tunlramp = tunl
                 endif
                 if( choice_leng .eq. 'origin' ) then
                     lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
                   ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
                 else
                     lz = min( vk*zi(i,k), tunlramp*lbulk )              
                 endif
                 lz = min(leng_max(k), lz)

                 dzinc = z(i,k-1)-z(i,k)
                 dl2n2 = lz*lz*n2(i,k)*dzinc
                 dl2s2 = lz*lz*s2(i,k)*dzinc
                 dwinc = -sh*dl2n2 + sm*dl2s2

                 lint = lint + dzinc
                 l2n2 = l2n2 + dl2n2
                 l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
                 l2s2 = l2s2 + dl2s2
                 wint = wint + dwinc

              end do 

              kt        = ktop(i,ncv+cntu+1) 
              ncvfin(i) = ncvfin(i) - 1
              cntu      = cntu + 1
        
          end if

          ! Again, calculate the contribution of potentially incorporatable CL
          ! top external interface of CL regime.

          if( choice_tunl .eq. 'rampcl' ) then
              tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
          elseif( choice_tunl .eq. 'rampsl' ) then
              tunlramp = ctunl*tunl
            ! tunlramp = 0.765_r8
          else
              tunlramp = tunl
          endif
          if( choice_leng .eq. 'origin' ) then
              lz = ( (vk*zi(i,kt))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
            ! lz = vk*zi(i,kt) / (1._r8+vk*zi(i,kt)/(tunlramp*lbulk))
          else
              lz = min( vk*zi(i,kt), tunlramp*lbulk )              
          endif
          lz = min(leng_max(kt), lz)

          dzinc = z(i,kt-1)-z(i,kt)
          dl2n2 = lz*lz*n2(i,kt)*dzinc
          dl2s2 = lz*lz*s2(i,kt)*dzinc
          dwinc = -sh*dl2n2 + sm*dl2s2

       end do   ! End of upward merging test 'do while' loop

       ! Update CL interface indices appropriately if any CL was merged.
       ! Note that below only updated the interface index of merged CL,
       ! not the original merging CL.  Updates of 'kbase' and 'ktop' of 
       ! the original merging CL  will be done after finishing downward
       ! extension also later.

       if( cntu .gt. 0 ) then
           do incv = 1, ncvfin(i) - ncv
              kbase(i,ncv+incv) = kbase(i,ncv+cntu+incv)
              ktop(i,ncv+incv)  = ktop(i,ncv+cntu+incv)
           end do
       end if

       ! ------------------------------ !
       ! Step 2. Extend the CL downward !
       ! ------------------------------ !
       
       if( kb .ne. pver + 1 ) then

           ! Calculate contribution of potentially incorporable CL base interface

           if( choice_tunl .eq. 'rampcl' ) then
               tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
           elseif( choice_tunl .eq. 'rampsl' ) then
               tunlramp = ctunl*tunl
             ! tunlramp = 0.765_r8
           else
               tunlramp = tunl
           endif
           if( choice_leng .eq. 'origin' ) then
               lz = ( (vk*zi(i,kb))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
             ! lz = vk*zi(i,kb) / (1._r8+vk*zi(i,kb)/(tunlramp*lbulk))
           else
               lz = min( vk*zi(i,kb), tunlramp*lbulk )              
           endif
           lz = min(leng_max(kb), lz)

           dzinc = z(i,kb-1)-z(i,kb)
           dl2n2 = lz*lz*n2(i,kb)*dzinc
           dl2s2 = lz*lz*s2(i,kb)*dzinc
           dwinc = -sh*dl2n2 + sm*dl2s2

           ! ------------ ! 
           ! Merging test !
           ! ------------ ! 

           ! In the below merging tests, I must keep '.and.(kb.ne.pver+1)',   
           ! since 'kb' is continuously updated within the 'do while' loop  
           ! whenever CL base is merged.

         ! do while( (  dwinc .gt. ( rinc*dzinc*wint/(lint+(1._r8-rinc)*dzinc)) ) &  ! Merging test based on wint
         ! do while( ( -dl2n2 .gt. (-rinc*dzinc*l2n2/(lint+(1._r8-rinc)*dzinc)) ) &  ! Merging test based on l2n2
         !             .and.(kb.ne.pver+1))
           do while( ( -dl2n2 .gt. (-rinc*l2n2/(1._r8-rinc)) ) &                     ! Integral merging test
                       .and.(kb.ne.pver+1))

              ! Add contributions from interfacial layer kb to CL interior 

              lint = lint + dzinc
              l2n2 = l2n2 + dl2n2
              l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
              l2s2 = l2s2 + dl2s2
              wint = wint + dwinc

              ! Extend the base external interface of CL downward after merging

              kb        =  kb + 1
              extend    = .true.
              extend_dn = .true.

              ! If the base external interface of extending CL is the same as the 
              ! base interior interface of the underlying CL, underlying CL  will
              ! be automatically merged. Then, reduce total number of CL by 1. 
              ! For a consistent treatment with 'upward' extension,  I should use
              ! 'kbinc = kbase(i,ncv-1) - 1' instead of 'ktop(i,ncv-1) + 1' below.
              ! However, it seems that these two methods produce the same results.
              ! Note also that in contrast to upward merging, the decrease of ncv
              ! should be performed here.
              ! Note that below formula correctly works even when upperlying CL 
              ! regime incorporates below SBCL.

              kbinc = 0
              if( ncv .gt. 1 ) kbinc = ktop(i,ncv-1) + 1
              if( kb .eq. kbinc ) then

                  do k =  ktop(i,ncv-1) + 1, kbase(i,ncv-1) - 1

                     if( choice_tunl .eq. 'rampcl' ) then
                         tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
                     elseif( choice_tunl .eq. 'rampsl' ) then
                         tunlramp = ctunl*tunl
                       ! tunlramp = 0.765_r8
                     else
                         tunlramp = tunl
                     endif
                     if( choice_leng .eq. 'origin' ) then
                         lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
                       ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
                     else
                         lz = min( vk*zi(i,k), tunlramp*lbulk )              
                     endif
                     lz = min(leng_max(k), lz)

                     dzinc = z(i,k-1)-z(i,k)
                     dl2n2 = lz*lz*n2(i,k)*dzinc
                     dl2s2 = lz*lz*s2(i,k)*dzinc
                     dwinc = -sh*dl2n2 + sm*dl2s2

                     lint = lint + dzinc
                     l2n2 = l2n2 + dl2n2
                     l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
                     l2s2 = l2s2 + dl2s2
                     wint = wint + dwinc

                  end do 

                  ! We are incorporating interior of CL ncv-1, so merge
                  ! this CL into the current CL.

                  kb        = kbase(i,ncv-1)
                  ncv       = ncv - 1
                  ncvfin(i) = ncvfin(i) -1
                  cntd      = cntd + 1

              end if

              ! Calculate the contribution of potentially incorporatable CL
              ! base external interface. Calculate separately when the base
              ! of extended CL is surface and non-surface.
             
              if( kb .eq. pver + 1 ) then 

                  if( bflxs(i) .gt. 0._r8 ) then 
                      ! Calculate stability functions of surface interfacial layer
                      gg = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
                      gh_surf = gg/(alph5-gg*alph3)
                    ! gh_surf = min(max(gh_surf,-0.28_r8),0.0233_r8)
                      gh_surf = min(max(gh_surf,-3.5334_r8),0.0233_r8)
                      sh_surf = alph5/(1._r8+alph3*gh_surf)
                      sm_surf = (alph1 + alph2*gh_surf)/(1._r8+alph3*gh_surf)/(1._r8+alph4*gh_surf)
                      ! Calculate surface interfacial layer contribution. By construction,
                      ! it exactly becomes 'dw_surf = -dl2n2_surf + ds2n2_surf'  
                      dlint_surf = z(i,pver)
                      dl2n2_surf = -vk*(z(i,pver)**2._r8)*bprod(i,pver+1)/(sh_surf*sqrt(tkes(i)))
                      dl2s2_surf =  vk*(z(i,pver)**2._r8)*sprod(i,pver+1)/(sm_surf*sqrt(tkes(i)))
                      dw_surf = (tkes(i)/b1)*z(i,pver) 
                  else
                      dlint_surf = 0._r8
                      dl2n2_surf = 0._r8
                      dl2s2_surf = 0._r8
                      dw_surf = 0._r8
                  end if
                  ! If (kb.eq.pver+1), updating of CL internal energetics should be 
                  ! performed here inside of 'do while' loop, since 'do while' loop
                  ! contains the constraint of '.and.(kb.ne.pver+1)',so updating of
                  ! CL internal energetics cannot be performed within this do while
                  ! loop when kb.eq.pver+1. Even though I updated all 'l2n2','l2s2',
                  ! 'wint' below, only the updated 'wint' is used in the following
                  ! numerical calculation.                
                  lint = lint + dlint_surf
                  l2n2 = l2n2 + dl2n2_surf
                  l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
                  l2s2 = l2s2 + dl2s2_surf 
                  wint = wint + dw_surf                
                
              else

                  if( choice_tunl .eq. 'rampcl' ) then
                      tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
                  elseif( choice_tunl .eq. 'rampsl' ) then
                      tunlramp = ctunl*tunl
                    ! tunlramp = 0.765_r8
                  else
                      tunlramp = tunl
                  endif
                  if( choice_leng .eq. 'origin' ) then
                      lz = ( (vk*zi(i,kb))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
                    ! lz = vk*zi(i,kb) / (1._r8+vk*zi(i,kb)/(tunlramp*lbulk))
                  else
                      lz = min( vk*zi(i,kb), tunlramp*lbulk )              
                  endif
                  lz = min(leng_max(kb), lz)

                  dzinc = z(i,kb-1)-z(i,kb)
                  dl2n2 = lz*lz*n2(i,kb)*dzinc
                  dl2s2 = lz*lz*s2(i,kb)*dzinc
                  dwinc = -sh*dl2n2 + sm*dl2s2

              end if

          end do ! End of merging test 'do while' loop

          if( (kb.eq.pver+1) .and. (ncv.ne.1) ) then 
               errstring = 'Major mistake zisocl: the CL based at surface is not indexed 1'
               return
          end if

       end if   ! Done with bottom extension of CL 

       ! Update CL interface indices appropriately if any CL was merged.
       ! Note that below only updated the interface index of merged CL,
       ! not the original merging CL.  Updates of 'kbase' and 'ktop' of 
       ! the original merging CL  will be done later below. I should 
       ! check in detail if below index updating is correct or not.   

       if( cntd .gt. 0 ) then
           do incv = 1, ncvfin(i) - ncv
              kbase(i,ncv+incv) = kbase(i,ncvinit+incv)
              ktop(i,ncv+incv)  = ktop(i,ncvinit+incv)
           end do
       end if

       ! Sanity check for positive wint.

       if( wint .lt. 0.01_r8 ) then
           wint = 0.01_r8
       end if

       ! -------------------------------------------------------------------------- !
       ! Finally update CL mean internal energetics including surface contribution  !
       ! after finishing all the CL extension-merging process.  As mentioned above, !
       ! there are two possible ways in the treatment of surface interfacial layer, !
       ! either through 'dw_surf' or 'dl2n2_surf and dl2s2_surf' by setting logical !
       ! variable 'use_dw_surf' =.true. or .false.    In any cases, we should avoid !
       ! double counting of surface interfacial layer and one single consistent way !
       ! should be used throughout the program.                                     !
       ! -------------------------------------------------------------------------- !

       if( extend ) then

           ktop(i,ncv)  = kt
           kbase(i,ncv) = kb

           ! ------------------------------------------------------ !   
           ! Step 1: Include surface interfacial layer contribution !
           ! ------------------------------------------------------ !        
          
           lbulk      = zi(i,kt) - zi(i,kb)
           lbulk      = min( lbulk, lbulk_max )
           dlint_surf = 0._r8
           dl2n2_surf = 0._r8
           dl2s2_surf = 0._r8
           dw_surf    = 0._r8
           if( kb .eq. pver + 1 ) then
               if( bflxs(i) .gt. 0._r8 ) then
                   ! Calculate stability functions of surface interfacial layer
                   gg = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
                   gh = gg/(alph5-gg*alph3)
                 ! gh = min(max(gh,-0.28_r8),0.0233_r8)
                   gh = min(max(gh,-3.5334_r8),0.0233_r8)
                   sh = alph5/(1._r8+alph3*gh)
                   sm = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
                   ! Calculate surface interfacial layer contribution. By construction,
                   ! it exactly becomes 'dw_surf = -dl2n2_surf + ds2n2_surf'  
                   dlint_surf = z(i,pver)
                   dl2n2_surf = -vk*(z(i,pver)**2._r8)*bprod(i,pver+1)/(sh*sqrt(tkes(i)))
                   dl2s2_surf =  vk*(z(i,pver)**2._r8)*sprod(i,pver+1)/(sm*sqrt(tkes(i)))
                   dw_surf    = (tkes(i)/b1)*z(i,pver) 
               else
                   lbulk = zi(i,kt) - z(i,pver)
                   lbulk = min( lbulk, lbulk_max )
               end if
           end if
           lint = dlint_surf
           l2n2 = dl2n2_surf
           l2s2 = dl2s2_surf
           wint = dw_surf
           if( use_dw_surf ) then
               l2n2 = 0._r8
               l2s2 = 0._r8
           else
               wint = 0._r8
           end if   
       
           ! -------------------------------------------------------------- !
           ! Step 2. Include the contribution of 'pure internal interfaces' !
           ! -------------------------------------------------------------- ! 
          
           do k = kt + 1, kb - 1
              if( choice_tunl .eq. 'rampcl' ) then
                  tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
              elseif( choice_tunl .eq. 'rampsl' ) then
                  tunlramp = ctunl*tunl
                ! tunlramp = 0.765_r8
              else
                  tunlramp = tunl
              endif
              if( choice_leng .eq. 'origin' ) then
                  lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
                ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
              else
                  lz = min( vk*zi(i,k), tunlramp*lbulk )              
              endif
              lz = min(leng_max(k), lz)
              dzinc = z(i,k-1) - z(i,k)
              lint = lint + dzinc
              l2n2 = l2n2 + lz*lz*n2(i,k)*dzinc
              l2s2 = l2s2 + lz*lz*s2(i,k)*dzinc
           end do

           ricll = min(l2n2/max(l2s2,ntzero),ricrit)
           trma = alph3*alph4*ricll+2._r8*b1*(alph2-alph4*alph5*ricll)
           trmb = ricll*(alph3+alph4)+2._r8*b1*(-alph5*ricll+alph1)
           trmc = ricll
           det = max(trmb*trmb-4._r8*trma*trmc,0._r8)
           gh = (-trmb + sqrt(det))/2._r8/trma
         ! gh = min(max(gh,-0.28_r8),0.0233_r8)
           gh = min(max(gh,-3.5334_r8),0.0233_r8)
           sh = alph5 / (1._r8+alph3*gh)
           sm = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
           ! Even though the 'wint' after finishing merging was positive, it is 
           ! possible that re-calculated 'wint' here is negative.  In this case,
           ! correct 'wint' to be a small positive number
           wint = max( wint - sh*l2n2 + sm*l2s2, 0.01_r8 )

       end if

       ! ---------------------------------------------------------------------- !
       ! Calculate final output variables of each CL (either has merged or not) !
       ! ---------------------------------------------------------------------- !

       lbrk(i,ncv) = lint
       wbrk(i,ncv) = wint/lint
       ebrk(i,ncv) = b1*wbrk(i,ncv)
       ebrk(i,ncv) = min(ebrk(i,ncv),tkemax)
       ricl(i,ncv) = ricll 
       ghcl(i,ncv) = gh 
       shcl(i,ncv) = sh
       smcl(i,ncv) = sm

       ! Increment counter for next CL. I should check if the increament of 'ncv'
       ! below is reasonable or not, since whenever CL is merged during downward
       ! extension process, 'ncv' is lowered down continuously within 'do' loop.
       ! But it seems that below 'ncv = ncv + 1' is perfectly correct.

       ncv = ncv + 1

    end do                   ! End of loop over each CL regime, ncv.

    ! ---------------------------------------------------------- !
    ! Re-initialize external interface indices which are not CLs !
    ! ---------------------------------------------------------- !

    do ncv = ncvfin(i) + 1, ncvmax
       ktop(i,ncv)  = 0
       kbase(i,ncv) = 0
    end do

    ! ------------------------------------------------ !
    ! Update CL interface identifiers, 'belongcv'      !
    ! CL external interfaces are also identified as CL !
    ! ------------------------------------------------ !

    do k = 1, pver + 1
       belongcv(i,k) = .false.
    end do

    do ncv = 1, ncvfin(i)
       do k = ktop(i,ncv), kbase(i,ncv)
          belongcv(i,k) = .true.
       end do
    end do

    return

    end subroutine zisocl

    real(r8) function compute_cubic(a,b,c)
    ! ------------------------------------------------------------------------- !
    ! Solve canonical cubic : x^3 + a*x^2 + b*x + c = 0,  x = sqrt(e)/sqrt(<e>) !
    ! Set x = max(xmin,x) at the end                                            ! 
    ! ------------------------------------------------------------------------- !
    implicit none
    real(r8), intent(in)     :: a, b, c
    real(r8)  qq, rr, dd, theta, aa, bb, x1, x2, x3
    real(r8), parameter      :: xmin = 1.e-2_r8
    
    qq = (a**2-3._r8*b)/9._r8 
    rr = (2._r8*a**3 - 9._r8*a*b + 27._r8*c)/54._r8
    
    dd = rr**2 - qq**3
    if( dd .le. 0._r8 ) then
        theta = acos(rr/qq**(3._r8/2._r8))
        x1 = -2._r8*sqrt(qq)*cos(theta/3._r8) - a/3._r8
        x2 = -2._r8*sqrt(qq)*cos((theta+2._r8*3.141592_r8)/3._r8) - a/3._r8
        x3 = -2._r8*sqrt(qq)*cos((theta-2._r8*3.141592_r8)/3._r8) - a/3._r8
        compute_cubic = max(max(max(x1,x2),x3),xmin)        
        return
    else
        if( rr .ge. 0._r8 ) then
            aa = -(sqrt(rr**2-qq**3)+rr)**(1._r8/3._r8)
        else
            aa =  (sqrt(rr**2-qq**3)-rr)**(1._r8/3._r8)
        endif
        if( aa .eq. 0._r8 ) then
            bb = 0._r8
        else
            bb = qq/aa
        endif
        compute_cubic = max((aa+bb)-a/3._r8,xmin) 
        return
    endif

    return
    end function compute_cubic

END MODULE eddy_diff