LCOV - code coverage report
Current view: | top level - physics/cam - eddy_diff.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 896 | 0.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 0 | 7 | 0.0 %

              Line data    Source code
           1             :   module eddy_diff
           2             :
           3             :   !--------------------------------------------------------------------------------- !
           4             :   !                                                                                  !
           5             :   ! The University of Washington Moist Turbulence Scheme to compute eddy diffusion   !
           6             :   ! coefficients associated with dry and moist turbulences in the whole              !
           7             :   ! atmospheric layers.                                                              !
           8             :   !                                                                                  !
           9             :   ! For detailed description of the code and its performances, see                   !
          10             :   !                                                                                  !
          11             :   ! 1.'A new moist turbulence parametrization in the Community Atmosphere Model'     !
          12             :   !    by Christopher S. Bretherton and Sungsu Park. J. Climate. 2009. 22. 3422-3448 !
          13             :   ! 2.'The University of Washington shallow convection and moist turbulence schemes  !
          14             :   !    and their impact on climate simulations with the Community Atmosphere Model'  !
          15             :   !    by Sungsu Park and Christopher S. Bretherton. J. Climate. 2009. 22. 3449-3469 !
          16             :   !                                                                                  !
          17             :   ! For questions on the scheme and code, send an email to                           !
          18             :   !     Sungsu Park      at sungsup@ucar.edu (tel: 303-497-1375)                     !
          19             :   !     Chris Bretherton at breth@washington.edu                                     !
          20             :   !                                                                                  !
          21             :   ! Developed by Chris Bretherton at the University of Washington, Seattle, WA.      !
          22             :   !              Sungsu Park      at the CGD/NCAR, Boulder, CO.                      !
          23             :   ! Last coded on May.2006, Dec.2009 by Sungsu Park.                                 !
          24             :   !                                                                                  !
          25             :   !--------------------------------------------------------------------------------- !
          26             :
          27             :   use wv_saturation,    only: qsat
          28             :
          29             :   implicit none
          30             :   private
          31             :   save
          32             :
          33             :   public :: init_eddy_diff
          34             :   public :: trbintd
          35             :   public :: caleddy
          36             :   public :: ncvmax
          37             :
          38             :   integer,          parameter :: r8 = selected_real_kind(12)    ! 8 byte real
          39             :   integer,          parameter :: i4 = selected_int_kind( 6)     ! 4 byte integer
          40             :   ! --------------------------------- !
          41             :   ! PBL Parameters used in the UW PBL !
          42             :   ! --------------------------------- !
          43             :
          44             :   character,        parameter :: sftype         = 'l'           ! Method for calculating saturation fraction
          45             :
          46             :   character(len=4), parameter :: choice_evhc    = 'maxi'        ! 'orig', 'ramp', 'maxi' : recommended to be used with choice_radf
          47             :   character(len=6), parameter :: choice_radf    = 'maxi'        ! 'orig', 'ramp', 'maxi' : recommended to be used with choice_evhc
          48             :   character(len=6), parameter :: choice_SRCL    = 'nonamb'      ! 'origin', 'remove', 'nonamb'
          49             :
          50             :   character(len=6), parameter :: choice_tunl    = 'rampcl'      ! 'origin', 'rampsl'(Sungsu), 'rampcl'(Chris)
          51             :   real(r8),         parameter :: ctunl          =  2._r8        !  Maximum asympt leng = ctunl*tunl when choice_tunl = 'rampsl(cl)'
          52             :                                                                 ! [ no unit ]
          53             :   character(len=6), parameter :: choice_leng    = 'origin'      ! 'origin', 'takemn'
          54             :   real(r8),         parameter :: cleng          =  3._r8        !  Order of 'leng' when choice_leng = 'origin' [ no unit ]
          55             :   character(len=6), parameter :: choice_tkes    = 'ibprod'      ! 'ibprod' (include tkes in computing bprod), 'ebprod'(exclude)
          56             :
          57             :   real(r8)                    :: lbulk_max      =  40.e3_r8     ! Maximum master length scale designed to address issues in the
          58             :                                                                 ! upper atmosphere where vertical model resolution is coarse [ m ].
          59             :                                                                 ! In order not to disturb turbulence characteristics in the lower
          60             :                                                                 ! troposphere, this should be set at least larger than ~ a few km.
          61             :   real(r8), allocatable       :: leng_max(:)                    ! Maximum length scale designed to address issues in the upper
          62             :                                                                 ! atmosphere.
          63             :
          64             :   ! Parameters for 'sedimentation-entrainment feedback' for liquid stratus
          65             :   ! If .false.,  no sedimentation entrainment feedback ( i.e., use default evhc )
          66             :
          67             :   logical,          parameter :: id_sedfact     = .false.
          68             :   real(r8),         parameter :: ased           =  9._r8        !  Valid only when id_sedfact = .true.
          69             :
          70             :   ! --------------------------------------------------------------------------------------------------- !
          71             :   ! Parameters governing entrainment efficiency A = a1l(i)*evhc, evhc = 1 + a2l * a3l * L * ql / jt2slv !
          72             :   ! Here, 'ql' is cloud-top LWC and 'jt2slv' is the jump in 'slv' across                                !
          73             :   ! the cloud-top entrainment zone ( across two grid layers to consider full mixture )                  !
          74             :   ! --------------------------------------------------------------------------------------------------- !
          75             :
          76             :   real(r8),         parameter :: a1l            =   0.10_r8     ! Dry entrainment efficiency for TKE closure
          77             :                                                                 ! a1l = 0.2*tunl*erat^-1.5,
          78             :                                                                 ! where erat = <e>/wstar^2 for dry CBL =  0.3.
          79             :
          80             :   real(r8),         parameter :: a1i            =   0.2_r8      ! Dry entrainment efficiency for wstar closure
          81             :   real(r8),         parameter :: ccrit          =   0.5_r8      ! Minimum allowable sqrt(tke)/wstar.
          82             :                                                                 ! Used in solving cubic equation for 'ebrk'
          83             :   real(r8),         parameter :: wstar3factcrit =   0.5_r8      ! 1/wstar3factcrit is the maximally allowed enhancement of
          84             :                                                                 ! 'wstar3' due to entrainment.
          85             :
          86             :   real(r8)                    :: a2l                            ! Moist entrainment enhancement param (recommended range : 10~30 )
          87             :   real(r8),         parameter :: a3l            =   0.8_r8      ! Approximation to a complicated thermodynamic parameters
          88             :
          89             :   real(r8),         parameter :: jbumin         =   .001_r8     ! Minimum buoyancy jump at an entrainment jump, [m/s2]
          90             :   real(r8),         parameter :: evhcmax        =   10._r8      ! Upper limit of evaporative enhancement factor
          91             :
          92             :   real(r8),         parameter :: onet           =   1._r8/3._r8 ! 1/3 power in wind gradient expression [ no unit ]
          93             :   integer                     :: ncvmax                         ! Max numbers of CLs (good to set to 'pver')
          94             :   real(r8),         parameter :: qmin           =   1.e-5_r8    ! Minimum grid-mean LWC counted as clouds [kg/kg]
          95             :   real(r8),         parameter :: ntzero         =   1.e-12_r8   ! Not zero (small positive number used in 's2')
          96             :   real(r8),         parameter :: b1             =   5.8_r8      ! TKE dissipation D = e^3/(b1*leng), e = b1*W.
          97             :   real(r8)                    :: b123                           ! b1**(2/3)
          98             :   real(r8),         parameter :: tunl           =   0.085_r8    ! Asympt leng = tunl*(turb lay depth)
          99             :   real(r8),         parameter :: alph1          =   0.5562_r8   ! alph1~alph5 : Galperin instability function parameters
         100             :   real(r8),         parameter :: alph2          =  -4.3640_r8   !               These coefficients are used to calculate
         101             :   real(r8),         parameter :: alph3          = -34.6764_r8   !               'sh' and 'sm' from 'gh'.
         102             :   real(r8),         parameter :: alph4          =  -6.1272_r8   !
         103             :   real(r8),         parameter :: alph5          =   0.6986_r8   !
         104             :   real(r8),         parameter :: ricrit         =   0.19_r8     ! Critical Richardson number for turbulence.
         105             :                                                                 ! Can be any value >= 0.19.
         106             :   real(r8),         parameter :: ae             =   1._r8       ! TKE transport efficiency [no unit]
         107             :   real(r8),         parameter :: rinc           =  -0.04_r8     ! Minimum W/<W> used for CL merging test
         108             :   real(r8),         parameter :: wpertmin       =   1.e-6_r8    ! Minimum PBL eddy vertical velocity perturbation
         109             :   real(r8),         parameter :: wfac           =   1._r8       ! Ratio of 'wpert' to sqrt(tke) for CL.
         110             :   real(r8),         parameter :: tfac           =   1._r8       ! Ratio of 'tpert' to (w't')/wpert for CL.
         111             :                                                                 ! Same ratio also used for q
         112             :   real(r8),         parameter :: fak            =   8.5_r8      ! Constant in surface temperature excess for stable STL.
         113             :                                                                 ! [ no unit ]
         114             :   real(r8),         parameter :: rcapmin        =   0.1_r8      ! Minimum allowable e/<e> in a CL
         115             :   real(r8),         parameter :: rcapmax        =   2.0_r8      ! Maximum allowable e/<e> in a CL
         116             :   real(r8),         parameter :: tkemax         =  20._r8       ! TKE is capped at tkemax [m2/s2]
         117             :
         118             :   logical,          parameter :: use_dw_surf    =  .true.       ! Used in 'zisocl'. Default is 'true'
         119             :                                                                 ! If 'true', surface interfacial energy does not contribute
         120             :                                                                 ! to the CL mean stability functions after finishing merging.
         121             :                                                                 ! For this case, 'dl2n2_surf' is only used for a merging test
         122             :                                                                 ! based on 'l2n2'
         123             :                                                                 ! If 'false',surface interfacial enery explicitly contribute to
         124             :                                                                 ! CL mean stability functions after finishing merging.
         125             :                                                                 ! For this case, 'dl2n2_surf' and 'dl2s2_surf' are directly used
         126             :                                                                 ! for calculating surface interfacial layer energetics
         127             :
         128             :   logical,          parameter :: set_qrlzero    =  .false.      ! .true. ( .false.) : turning-off ( on) radiative-turbulence
         129             :                                                                 ! interaction by setting qrl = 0.
         130             :
         131             :   ! ------------------------------------------------------- !
         132             :   ! PBL constants set using values from other parts of code !
         133             :   ! ------------------------------------------------------- !
         134             :
         135             :   real(r8)                    :: cpair                          ! Specific heat of dry air
         136             :   real(r8)                    :: rair                           ! Gas const for dry air
         137             :   real(r8)                    :: zvir                           ! rh2o/rair - 1
         138             :   real(r8)                    :: latvap                         ! Latent heat of vaporization
         139             :   real(r8)                    :: latice                         ! Latent heat of fusion
         140             :   real(r8)                    :: latsub                         ! Latent heat of sublimation
         141             :   real(r8)                    :: g                              ! Gravitational acceleration
         142             :   real(r8)                    :: vk                             ! Von Karman's constant
         143             :
         144             :   integer                     :: ntop_turb                      ! Top interface level to which turbulent vertical diffusion
         145             :                                                                 ! is applied ( = 1 )
         146             :   integer                     :: nbot_turb                      ! Bottom interface level to which turbulent vertical diff
         147             :                                                                 ! is applied ( = pver )
         148             :
         149             :   CONTAINS
         150             :
         151             :   !============================================================================ !
         152             :   !                                                                             !
         153             :   !============================================================================ !
         154             :
         155           0 :   subroutine init_eddy_diff( pver, gravx, cpairx, rairx, zvirx, &
         156             :                              latvapx, laticex, ntop_eddy, nbot_eddy, vkx, &
         157           0 :                              eddy_lbulk_max, leng_max_in, &
         158             :                              eddy_moist_entrain_a2l, errstring)
         159             :     !---------------------------------------------------------------- !
         160             :     ! Purpose:                                                        !
         161             :     ! Initialize time independent constants/variables of PBL package. !
         162             :     !---------------------------------------------------------------- !
         163             :
         164             :     ! --------- !
         165             :     ! Arguments !
         166             :     ! --------- !
         167             :     integer,  intent(in) :: pver       ! Number of vertical layers
         168             :     integer,  intent(in) :: ntop_eddy  ! Top interface level to which eddy vertical diffusivity is applied ( = 1 )
         169             :     integer,  intent(in) :: nbot_eddy  ! Bottom interface level to which eddy vertical diffusivity is applied ( = pver )
         170             :     real(r8), intent(in) :: gravx      ! Acceleration of gravity
         171             :     real(r8), intent(in) :: cpairx     ! Specific heat of dry air
         172             :     real(r8), intent(in) :: rairx      ! Gas constant for dry air
         173             :     real(r8), intent(in) :: zvirx      ! rh2o/rair - 1
         174             :     real(r8), intent(in) :: latvapx    ! Latent heat of vaporization
         175             :     real(r8), intent(in) :: laticex    ! Latent heat of fusion
         176             :     real(r8), intent(in) :: vkx        ! Von Karman's constant
         177             :     real(r8), intent(in) :: eddy_lbulk_max ! Maximum master length scale
         178             :     real(r8), intent(in) :: leng_max_in(pver) ! Maximum length scale for upper atmosphere
         179             :     real(r8), intent(in) :: eddy_moist_entrain_a2l ! Moist entrainment enhancement param
         180             :
         181             :     character(len=*), intent(out) :: errstring
         182             :
         183             :     integer              :: k          ! Vertical loop index
         184             :
         185           0 :     errstring = ""
         186             :
         187             :     ! --------------- !
         188             :     ! Basic constants !
         189             :     ! --------------- !
         190             :
         191           0 :     ncvmax    = pver
         192             :
         193           0 :     cpair     = cpairx
         194           0 :     rair      = rairx
         195           0 :     g         = gravx
         196           0 :     zvir      = zvirx
         197           0 :     latvap    = latvapx
         198           0 :     latice    = laticex
         199           0 :     latsub    = latvap + latice
         200           0 :     vk        = vkx
         201           0 :     ntop_turb = ntop_eddy
         202           0 :     nbot_turb = nbot_eddy
         203           0 :     b123      = b1**(2._r8/3._r8)
         204           0 :     a2l       = eddy_moist_entrain_a2l
         205             :
         206           0 :     lbulk_max = eddy_lbulk_max
         207             :
         208           0 :     allocate(leng_max(pver))
         209           0 :     leng_max = leng_max_in
         210             :
         211           0 :   end subroutine init_eddy_diff
         212             :
         213             :   !=============================================================================== !
         214             :   !                                                                                !
         215             :   !=============================================================================== !
         216             :
         217           0 :   subroutine sfdiag( pcols   , pver    , ncol    , qt      , ql      , sl      , &
         218           0 :                      pi      , pm      , zi      , cld     , sfi     , sfuh    , &
         219           0 :                      sflh    , slslope , qtslope )
         220             :     !----------------------------------------------------------------------- !
         221             :     !                                                                        !
         222             :     ! Purpose: Interface for calculating saturation fractions  at upper and  !
         223             :     !          lower-half layers, & interfaces for use by turbulence scheme  !
         224             :     !                                                                        !
         225             :     ! Method : Various but 'l' should be chosen for consistency.             !
         226             :     !                                                                        !
         227             :     ! Author : B. Stevens and C. Bretherton (August 2000)                    !
         228             :     !          Sungsu Park. August 2006.                                     !
         229             :     !                       May.   2008.                                     !
         230             :     !                                                                        !
         231             :     ! S.Park : The computed saturation fractions are repeatedly              !
         232             :     !          used to compute buoyancy coefficients in'trbintd' & 'caleddy'.!
         233             :     !----------------------------------------------------------------------- !
         234             :
         235             :     implicit none
         236             :
         237             :     ! --------------- !
         238             :     ! Input arguments !
         239             :     ! --------------- !
         240             :
         241             :     integer,  intent(in)  :: pcols               ! Number of atmospheric columns
         242             :     integer,  intent(in)  :: pver                ! Number of atmospheric layers
         243             :     integer,  intent(in)  :: ncol                ! Number of atmospheric columns
         244             :
         245             :     real(r8), intent(in)  :: sl(pcols,pver)      ! Liquid water static energy [ J/kg ]
         246             :     real(r8), intent(in)  :: qt(pcols,pver)      ! Total water specific humidity [ kg/kg ]
         247             :     real(r8), intent(in)  :: ql(pcols,pver)      ! Liquid water specific humidity [ kg/kg ]
         248             :     real(r8), intent(in)  :: pi(pcols,pver+1)    ! Interface pressures [ Pa ]
         249             :     real(r8), intent(in)  :: pm(pcols,pver)      ! Layer mid-point pressures [ Pa ]
         250             :     real(r8), intent(in)  :: zi(pcols,pver+1)    ! Interface heights [ m ]
         251             :     real(r8), intent(in)  :: cld(pcols,pver)     ! Stratiform cloud fraction [ fraction ]
         252             :     real(r8), intent(in)  :: slslope(pcols,pver) ! Slope of 'sl' in each layer
         253             :     real(r8), intent(in)  :: qtslope(pcols,pver) ! Slope of 'qt' in each layer
         254             :
         255             :     ! ---------------- !
         256             :     ! Output arguments !
         257             :     ! ---------------- !
         258             :
         259             :     real(r8), intent(out) :: sfi(pcols,pver+1)   ! Interfacial layer saturation fraction [ fraction ]
         260             :     real(r8), intent(out) :: sfuh(pcols,pver)    ! Saturation fraction in upper half-layer [ fraction ]
         261             :     real(r8), intent(out) :: sflh(pcols,pver)    ! Saturation fraction in lower half-layer [ fraction ]
         262             :
         263             :     ! --------------- !
         264             :     ! Local Variables !
         265             :     ! --------------- !
         266             :
         267             :     integer               :: i                   ! Longitude index
         268             :     integer               :: k                   ! Vertical index
         269             :     integer               :: km1                 ! k-1
         270             :     integer               :: status              ! Status returned by function calls
         271             :     real(r8)              :: sltop, slbot        ! sl at top/bot of grid layer
         272             :     real(r8)              :: qttop, qtbot        ! qt at top/bot of grid layer
         273             :     real(r8)              :: tltop, tlbot  ! Liquid water temperature at top/bot of grid layer
         274             :     real(r8)              :: qxtop, qxbot        ! Sat excess at top/bot of grid layer
         275             :     real(r8)              :: qxm                 ! Sat excess at midpoint
         276             :     real(r8)              :: es               ! Saturation vapor pressure
         277             :     real(r8)              :: qs               ! Saturation spec. humidity
         278           0 :     real(r8)              :: cldeff(pcols,pver)  ! Effective Cloud Fraction [ fraction ]
         279             :
         280             :     ! ----------------------- !
         281             :     ! Main Computation Begins !
         282             :     ! ----------------------- !
         283             :
         284           0 :     sfi(1:ncol,:)    = 0._r8
         285           0 :     sfuh(1:ncol,:)   = 0._r8
         286           0 :     sflh(1:ncol,:)   = 0._r8
         287           0 :     cldeff(1:ncol,:) = 0._r8
         288             :
         289             :     select case (sftype)
         290             :     case ('d')
         291             :        ! ----------------------------------------------------------------------- !
         292             :        ! Simply use the given stratus fraction ('horizontal' cloud partitioning) !
         293             :        ! ----------------------------------------------------------------------- !
         294             :        do k = ntop_turb + 1, nbot_turb
         295             :           km1 = k - 1
         296             :           do i = 1, ncol
         297             :              sfuh(i,k) = cld(i,k)
         298             :              sflh(i,k) = cld(i,k)
         299             :              sfi(i,k)  = 0.5_r8 * ( sflh(i,km1) + min( sflh(i,km1), sfuh(i,k) ) )
         300             :           end do
         301             :        end do
         302           0 :        do i = 1, ncol
         303             :           sfi(i,pver+1) = sflh(i,pver)
         304             :        end do
         305             :     case ('l')
         306             :        ! ------------------------------------------ !
         307             :        ! Use modified stratus fraction partitioning !
         308             :        ! ------------------------------------------ !
         309           0 :        do k = ntop_turb + 1, nbot_turb
         310           0 :           km1 = k - 1
         311           0 :           do i = 1, ncol
         312           0 :              cldeff(i,k) = cld(i,k)
         313           0 :              sfuh(i,k)   = cld(i,k)
         314           0 :              sflh(i,k)   = cld(i,k)
         315           0 :              if( ql(i,k) .lt. qmin ) then
         316           0 :                  sfuh(i,k) = 0._r8
         317           0 :                  sflh(i,k) = 0._r8
         318             :              end if
         319             :            ! Modification : The contribution of ice should be carefully considered.
         320             :              if( choice_evhc .eq. 'ramp' .or. choice_radf .eq. 'ramp' ) then
         321             :                  cldeff(i,k) = cld(i,k) * min( ql(i,k) / qmin, 1._r8 )
         322             :                  sfuh(i,k)   = cldeff(i,k)
         323             :                  sflh(i,k)   = cldeff(i,k)
         324             :              elseif( choice_evhc .eq. 'maxi' .or. choice_radf .eq. 'maxi' ) then
         325           0 :                  cldeff(i,k) = cld(i,k)
         326           0 :                  sfuh(i,k)   = cldeff(i,k)
         327           0 :                  sflh(i,k)   = cldeff(i,k)
         328             :              endif
         329             :            ! At the stratus top, take the minimum interfacial saturation fraction
         330           0 :              sfi(i,k) = 0.5_r8 * ( sflh(i,km1) + min( sfuh(i,k), sflh(i,km1) ) )
         331             :            ! Modification : Currently sfi at the top and surface interfaces are set to be zero.
         332             :            !                Also, sfuh and sflh in the top model layer is set to be zero.
         333             :            !                However, I may need to set
         334             :            !                         do i = 1, ncol
         335             :            !                            sfi(i,pver+1) = sflh(i,pver)
         336             :            !                         end do
         337             :            !                for treating surface-based fog.
         338             :            ! OK. I added below block similar to the other cases.
         339             :           end do
         340             :        end do
         341           0 :        do i = 1, ncol
         342           0 :           sfi(i,pver+1) = sflh(i,pver)
         343             :        end do
         344             :     case ('u')
         345             :        ! ------------------------------------------------------------------------- !
         346             :        ! Use unsaturated buoyancy - since sfi, sfuh, sflh have already been zeroed !
         347             :        ! nothing more need be done for this case.                                  !
         348             :        ! ------------------------------------------------------------------------- !
         349             :     case ('z')
         350             :        ! ------------------------------------------------------------------------- !
         351             :        ! Calculate saturation fraction based on whether the air just above or just !
         352             :        ! below the interface is saturated, i.e. with vertical cloud partitioning.  !
         353             :        ! The saturation fraction of the interfacial layer between mid-points k and !
         354             :        ! k+1 is computed by averaging the saturation fraction   of the half-layers !
         355             :        ! above and below the interface,  with a special provision   for cloud tops !
         356             :        ! (more cloud in the half-layer below than in the half-layer above).In each !
         357             :        ! half-layer, vertical partitioning of  cloud based on the slopes diagnosed !
         358             :        ! above is used.     Loop down through the layers, computing the saturation !
         359             :        ! fraction in each half-layer (sfuh for upper half, sflh for lower half).   !
         360             :        ! Once sfuh(i,k) is computed, use with sflh(i,k-1) to determine  saturation !
         361             :        ! fraction sfi(i,k) for interfacial layer k-0.5.                            !
         362             :        ! This is 'not' chosen for full consistent treatment of stratus fraction in !
         363             :        ! all physics schemes.                                                      !
         364             :        ! ------------------------------------------------------------------------- !
         365             :        do k = ntop_turb + 1, nbot_turb
         366             :           km1 = k - 1
         367             :           do i = 1, ncol
         368             :            ! Compute saturation excess at the mid-point of layer k
         369             :              sltop    = sl(i,k) + slslope(i,k) * ( pi(i,k) - pm(i,k) )
         370             :              qttop    = qt(i,k) + qtslope(i,k) * ( pi(i,k) - pm(i,k) )
         371             :              tltop = ( sltop - g * zi(i,k) ) / cpair
         372             :              call qsat( tltop, pi(i,k), es, qs)
         373             :              qxtop    = qttop - qs
         374             :              slbot    = sl(i,k) + slslope(i,k) * ( pi(i,k+1) - pm(i,k) )
         375             :              qtbot    = qt(i,k) + qtslope(i,k) * ( pi(i,k+1) - pm(i,k) )
         376             :              tlbot = ( slbot - g * zi(i,k+1) ) / cpair
         377             :              call qsat( tlbot, pi(i,k+1), es, qs)
         378             :              qxbot    = qtbot - qs
         379             :              qxm      = qxtop + ( qxbot - qxtop ) * ( pm(i,k) - pi(i,k) ) / ( pi(i,k+1) - pi(i,k) )
         380             :            ! Find the saturation fraction sfuh(i,k) of the upper half of layer k.
         381             :              if( ( qxtop .lt. 0._r8 ) .and. ( qxm .lt. 0._r8 ) ) then
         382             :                    sfuh(i,k) = 0._r8
         383             :              else if( ( qxtop .gt. 0._r8 ) .and. ( qxm .gt. 0._r8 ) ) then
         384             :                    sfuh(i,k) = 1._r8
         385             :              else ! Either qxm < 0 and qxtop > 0 or vice versa
         386             :                    sfuh(i,k) = max( qxtop, qxm ) / abs( qxtop - qxm )
         387             :              end if
         388             :            ! Combine with sflh(i) (still for layer k-1) to get interfac layer saturation fraction
         389             :              sfi(i,k) = 0.5_r8 * ( sflh(i,k-1) + min( sflh(i,k-1), sfuh(i,k) ) )
         390             :            ! Update sflh to be for the lower half of layer k.
         391             :              if( ( qxbot .lt. 0._r8 ) .and. ( qxm .lt. 0._r8 ) ) then
         392             :                    sflh(i,k) = 0._r8
         393             :              else if( ( qxbot .gt. 0._r8 ) .and. ( qxm .gt. 0._r8 ) ) then
         394             :                    sflh(i,k) = 1._r8
         395             :              else ! Either qxm < 0 and qxbot > 0 or vice versa
         396             :                    sflh(i,k) = max( qxbot, qxm ) / abs( qxbot - qxm )
         397             :              end if
         398             :           end do  ! i
         399             :        end do ! k
         400             :        do i = 1, ncol
         401             :           sfi(i,pver+1) = sflh(i,pver)  ! Saturation fraction in the lowest half-layer.
         402             :        end do
         403             :     end select
         404             :
         405           0 :   return
         406             :   end subroutine sfdiag
         407             :
         408             :   !=============================================================================== !
         409             :   !                                                                                !
         410             :   !=============================================================================== !
         411             :
         412           0 :   subroutine trbintd( pcols   , pver    , ncol    ,                               &
         413           0 :                       z       , u       , v       ,                               &
         414           0 :                       t       , pmid    ,                                         &
         415           0 :                       s2      , n2      , ri      ,                               &
         416           0 :                       zi      , pi      , cld     ,                               &
         417           0 :                       qt      , qv      , ql      , qi      , sfi     , sfuh    , &
         418           0 :                       sflh    , sl      , slv     , slslope , qtslope ,           &
         419           0 :                       chs     , chu     , cms     , cmu     )
         420             :     !----------------------------------------------------------------------- !
         421             :     ! Purpose: Calculate buoyancy coefficients at all interfaces including   !
         422             :     !          surface. Also, computes the profiles of ( sl,qt,n2,s2,ri ).   !
         423             :     !          Note that (n2,s2,ri) are defined at each interfaces except    !
         424             :     !          surface.                                                      !
         425             :     !                                                                        !
         426             :     ! Author: B. Stevens  ( Extracted from pbldiff, August, 2000 )           !
         427             :     !         Sungsu Park ( August 2006, May. 2008 )                         !
         428             :     !----------------------------------------------------------------------- !
         429             :
         430             :     implicit none
         431             :
         432             :     ! --------------- !
         433             :     ! Input arguments !
         434             :     ! --------------- !
         435             :
         436             :     integer,  intent(in)  :: pcols                            ! Number of atmospheric columns
         437             :     integer,  intent(in)  :: pver                             ! Number of atmospheric layers
         438             :     integer,  intent(in)  :: ncol                             ! Number of atmospheric columns
         439             :     real(r8), intent(in)  :: z(pcols,pver)                    ! Layer mid-point height above surface [ m ]
         440             :     real(r8), intent(in)  :: u(pcols,pver)                    ! Layer mid-point u [ m/s ]
         441             :     real(r8), intent(in)  :: v(pcols,pver)                    ! Layer mid-point v [ m/s ]
         442             :     real(r8), intent(in)  :: t(pcols,pver)                    ! Layer mid-point temperature [ K ]
         443             :     real(r8), intent(in)  :: pmid(pcols,pver)                 ! Layer mid-point pressure [ Pa ]
         444             :     real(r8), intent(in)  :: zi(pcols,pver+1)                 ! Interface height [ m ]
         445             :     real(r8), intent(in)  :: pi(pcols,pver+1)                 ! Interface pressure [ Pa ]
         446             :     real(r8), intent(in)  :: cld(pcols,pver)                  ! Stratus fraction
         447             :     real(r8), intent(in)  :: qv(pcols,pver)                   ! Water vapor specific humidity [ kg/kg ]
         448             :     real(r8), intent(in)  :: ql(pcols,pver)                   ! Liquid water specific humidity [ kg/kg ]
         449             :     real(r8), intent(in)  :: qi(pcols,pver)                   ! Ice water specific humidity [ kg/kg ]
         450             :
         451             :     ! ---------------- !
         452             :     ! Output arguments !
         453             :     ! ---------------- !
         454             :
         455             :     real(r8), intent(out) :: s2(pcols,pver)                   ! Interfacial ( except surface ) shear squared [ s-2 ]
         456             :     real(r8), intent(out) :: n2(pcols,pver)                   ! Interfacial ( except surface ) buoyancy frequency [ s-2 ]
         457             :     real(r8), intent(out) :: ri(pcols,pver)                   ! Interfacial ( except surface ) Richardson number, 'n2/s2'
         458             :
         459             :     real(r8), intent(out) :: qt(pcols,pver)                   ! Total specific humidity [ kg/kg ]
         460             :     real(r8), intent(out) :: sfi(pcols,pver+1)                ! Interfacial layer saturation fraction [ fraction ]
         461             :     real(r8), intent(out) :: sfuh(pcols,pver)                 ! Saturation fraction in upper half-layer [ fraction ]
         462             :     real(r8), intent(out) :: sflh(pcols,pver)                 ! Saturation fraction in lower half-layer [ fraction ]
         463             :     real(r8), intent(out) :: sl(pcols,pver)                   ! Liquid water static energy [ J/kg ]
         464             :     real(r8), intent(out) :: slv(pcols,pver)                  ! Liquid water virtual static energy [ J/kg ]
         465             :
         466             :     real(r8), intent(out) :: chu(pcols,pver+1)                ! Heat buoyancy coef for dry states at all interfaces, finally.
         467             :                                                               ! [ unit ? ]
         468             :     real(r8), intent(out) :: chs(pcols,pver+1)                ! heat buoyancy coef for sat states at all interfaces, finally.
         469             :                                                               ! [ unit ? ]
         470             :     real(r8), intent(out) :: cmu(pcols,pver+1)                ! Moisture buoyancy coef for dry states at all interfaces, finally.
         471             :                                                               ! [ unit ? ]
         472             :     real(r8), intent(out) :: cms(pcols,pver+1)                ! Moisture buoyancy coef for sat states at all interfaces, finally.
         473             :                                                               ! [ unit ? ]
         474             :     real(r8), intent(out) :: slslope(pcols,pver)              ! Slope of 'sl' in each layer
         475             :     real(r8), intent(out) :: qtslope(pcols,pver)              ! Slope of 'qt' in each layer
         476             :
         477             :     ! --------------- !
         478             :     ! Local Variables !
         479             :     ! --------------- !
         480             :
         481             :     integer               :: i                                ! Longitude index
         482             :     integer               :: k, km1                           ! Level index
         483             :     integer               :: status                           ! Status returned by function calls
         484             :
         485           0 :     real(r8)              :: qs(pcols,pver)                   ! Saturation specific humidity
         486           0 :     real(r8)              :: es(pcols,pver)                   ! Saturation vapor pressure
         487           0 :     real(r8)              :: gam(pcols,pver)                  ! (l/cp)*(d(qs)/dT)
         488             :     real(r8)              :: rdz                              ! 1 / (delta z) between midpoints
         489             :     real(r8)              :: dsldz                            ! 'delta sl / delta z' at interface
         490             :     real(r8)              :: dqtdz                            ! 'delta qt / delta z' at interface
         491             :     real(r8)              :: ch                               ! 'sfi' weighted ch at the interface
         492             :     real(r8)              :: cm                               ! 'sfi' weighted cm at the interface
         493             :     real(r8)              :: bfact                            ! Buoyancy factor in n2 calculations
         494             :     real(r8)              :: product                          ! Intermediate vars used to find slopes
         495             :     real(r8)              :: dsldp_a, dqtdp_a                 ! Slopes across interface above
         496           0 :     real(r8)              :: dsldp_b(pcols), dqtdp_b(pcols)   ! Slopes across interface below
         497             :
         498             :     ! ----------------------- !
         499             :     ! Main Computation Begins !
         500             :     ! ----------------------- !
         501             :
         502             :     ! Calculate conservative scalars (qt,sl,slv) and buoyancy coefficients at the layer mid-points.
         503             :     ! Note that 'ntop_turb = 1', 'nbot_turb = pver'
         504           0 :     do k = ntop_turb, nbot_turb
         505           0 :        call qsat( t(1:ncol,k), pmid(1:ncol,k), es(1:ncol,k), qs(1:ncol,k), ncol, gam=gam(1:ncol,k))
         506           0 :        do i = 1, ncol
         507           0 :           qt(i,k)  = qv(i,k) + ql(i,k) + qi(i,k)
         508           0 :           sl(i,k)  = cpair * t(i,k) + g * z(i,k) - latvap * ql(i,k) - latsub * qi(i,k)
         509           0 :           slv(i,k) = sl(i,k) * ( 1._r8 + zvir * qt(i,k) )
         510             :         ! Thermodynamic coefficients for buoyancy flux - in this loop these are
         511             :         ! calculated at mid-points; later,  they will be averaged to interfaces,
         512             :         ! where they will ultimately be used.  At the surface, the coefficients
         513             :         ! are taken from the lowest mid point.
         514           0 :           bfact    = g / ( t(i,k) * ( 1._r8 + zvir * qv(i,k) - ql(i,k) - qi(i,k) ) )
         515           0 :           chu(i,k) = ( 1._r8 + zvir * qt(i,k) ) * bfact / cpair
         516           0 :           chs(i,k) = ( ( 1._r8 + ( 1._r8 + zvir ) * gam(i,k) * cpair * t(i,k) / latvap ) / ( 1._r8 + gam(i,k) ) ) * bfact / cpair
         517           0 :           cmu(i,k) = zvir * bfact * t(i,k)
         518           0 :           cms(i,k) = latvap * chs(i,k)  -  bfact * t(i,k)
         519             :        end do
         520             :     end do
         521             :
         522           0 :     do i = 1, ncol
         523           0 :        chu(i,pver+1) = chu(i,pver)
         524           0 :        chs(i,pver+1) = chs(i,pver)
         525           0 :        cmu(i,pver+1) = cmu(i,pver)
         526           0 :        cms(i,pver+1) = cms(i,pver)
         527             :     end do
         528             :
         529             :     ! Compute slopes of conserved variables sl, qt within each layer k.
         530             :     ! 'a' indicates the 'above' gradient from layer k-1 to layer k and
         531             :     ! 'b' indicates the 'below' gradient from layer k   to layer k+1.
         532             :     ! We take a smaller (in absolute value)  of these gradients as the
         533             :     ! slope within layer k. If they have opposite signs,   gradient in
         534             :     ! layer k is taken to be zero. I should re-consider whether   this
         535             :     ! profile reconstruction is the best or not.
         536             :     ! This is similar to the profile reconstruction used in the UWShCu.
         537             :
         538           0 :     do i = 1, ncol
         539             :      ! Slopes at endpoints determined by extrapolation
         540           0 :        slslope(i,pver) = ( sl(i,pver) - sl(i,pver-1) ) / ( pmid(i,pver) - pmid(i,pver-1) )
         541           0 :        qtslope(i,pver) = ( qt(i,pver) - qt(i,pver-1) ) / ( pmid(i,pver) - pmid(i,pver-1) )
         542           0 :        slslope(i,1)    = ( sl(i,2) - sl(i,1) ) / ( pmid(i,2) - pmid(i,1) )
         543           0 :        qtslope(i,1)    = ( qt(i,2) - qt(i,1) ) / ( pmid(i,2) - pmid(i,1) )
         544           0 :        dsldp_b(i)      = slslope(i,1)
         545           0 :        dqtdp_b(i)      = qtslope(i,1)
         546             :     end do
         547             :
         548           0 :     do k = 2, pver - 1
         549           0 :        do i = 1, ncol
         550           0 :           dsldp_a    = dsldp_b(i)
         551           0 :           dqtdp_a    = dqtdp_b(i)
         552           0 :           dsldp_b(i) = ( sl(i,k+1) - sl(i,k) ) / ( pmid(i,k+1) - pmid(i,k) )
         553           0 :           dqtdp_b(i) = ( qt(i,k+1) - qt(i,k) ) / ( pmid(i,k+1) - pmid(i,k) )
         554           0 :           product    = dsldp_a * dsldp_b(i)
         555           0 :           if( product .le. 0._r8 ) then
         556           0 :               slslope(i,k) = 0._r8
         557           0 :           else if( product .gt. 0._r8 .and. dsldp_a .lt. 0._r8 ) then
         558           0 :               slslope(i,k) = max( dsldp_a, dsldp_b(i) )
         559           0 :           else if( product .gt. 0._r8 .and. dsldp_a .gt. 0._r8 ) then
         560           0 :               slslope(i,k) = min( dsldp_a, dsldp_b(i) )
         561             :           end if
         562           0 :           product = dqtdp_a*dqtdp_b(i)
         563           0 :           if( product .le. 0._r8 ) then
         564           0 :               qtslope(i,k) = 0._r8
         565           0 :           else if( product .gt. 0._r8 .and. dqtdp_a .lt. 0._r8 ) then
         566           0 :               qtslope(i,k) = max( dqtdp_a, dqtdp_b(i) )
         567           0 :           else if( product .gt. 0._r8 .and. dqtdp_a .gt. 0._r8 ) then
         568           0 :               qtslope(i,k) = min( dqtdp_a, dqtdp_b(i) )
         569             :           end if
         570             :        end do ! i
         571             :     end do ! k
         572             :
         573             :     !  Compute saturation fraction at the interfacial layers for use in buoyancy
         574             :     !  flux computation.
         575             :
         576             :     call sfdiag( pcols  , pver    , ncol    , qt      , ql      , sl      , &
         577             :                  pi     , pmid    , zi      , cld     , sfi     , sfuh    , &
         578           0 :                  sflh   , slslope , qtslope )
         579             :
         580             :     ! Calculate buoyancy coefficients at all interfaces (1:pver+1) and (n2,s2,ri)
         581             :     ! at all interfaces except surface. Note 'nbot_turb = pver', 'ntop_turb = 1'.
         582             :     ! With the previous definition of buoyancy coefficients at the surface, the
         583             :     ! resulting buoyancy coefficients at the top and surface interfaces becomes
         584             :     ! identical to the buoyancy coefficients at the top and bottom layers. Note
         585             :     ! that even though the dimension of (s2,n2,ri) is 'pver',  they are defined
         586             :     ! at interfaces ( not at the layer mid-points ) except the surface.
         587             :
         588           0 :     do k = nbot_turb, ntop_turb + 1, -1
         589           0 :        km1 = k - 1
         590           0 :        do i = 1, ncol
         591           0 :           rdz      = 1._r8 / ( z(i,km1) - z(i,k) )
         592           0 :           dsldz    = ( sl(i,km1) - sl(i,k) ) * rdz
         593           0 :           dqtdz    = ( qt(i,km1) - qt(i,k) ) * rdz
         594           0 :           chu(i,k) = ( chu(i,km1) + chu(i,k) ) * 0.5_r8
         595           0 :           chs(i,k) = ( chs(i,km1) + chs(i,k) ) * 0.5_r8
         596           0 :           cmu(i,k) = ( cmu(i,km1) + cmu(i,k) ) * 0.5_r8
         597           0 :           cms(i,k) = ( cms(i,km1) + cms(i,k) ) * 0.5_r8
         598           0 :           ch       = chu(i,k) * ( 1._r8 - sfi(i,k) ) + chs(i,k) * sfi(i,k)
         599           0 :           cm       = cmu(i,k) * ( 1._r8 - sfi(i,k) ) + cms(i,k) * sfi(i,k)
         600           0 :           n2(i,k)  = ch * dsldz +  cm * dqtdz
         601           0 :           s2(i,k)  = ( ( u(i,km1) - u(i,k) )**2 + ( v(i,km1) - v(i,k) )**2) * rdz**2
         602           0 :           s2(i,k)  = max( ntzero, s2(i,k) )
         603           0 :           ri(i,k)  = n2(i,k) / s2(i,k)
         604             :        end do
         605             :     end do
         606           0 :     do i = 1, ncol
         607           0 :        n2(i,1) = n2(i,2)
         608           0 :        s2(i,1) = s2(i,2)
         609           0 :        ri(i,1) = ri(i,2)
         610             :     end do
         611             :
         612           0 :   return
         613             :
         614             :   end subroutine trbintd
         615             :
         616             :     ! ---------------------------------------------------------------------------- !
         617             :     !                                                                              !
         618             :     ! The University of Washington Moist Turbulence Scheme                         !
         619             :     !                                                                              !
         620             :     ! Authors : Chris Bretherton at the University of Washington, Seattle, WA      !
         621             :     !           Sungsu Park at the CGD/NCAR, Boulder, CO                           !
         622             :     !                                                                              !
         623             :     ! ---------------------------------------------------------------------------- !
         624             :
         625           0 :     subroutine caleddy( pcols        , pver         , ncol        ,                             &
         626           0 :                         sl           , qt           , ql          , slv        , u            , &
         627           0 :                         v            , pi           , z           , zi         ,                &
         628           0 :                         qflx         , shflx        , slslope     , qtslope    ,                &
         629           0 :                         chu          , chs          , cmu         , cms        , sfuh         , &
         630           0 :                         sflh         , n2           , s2          , ri         , rrho         , &
         631           0 :                         pblh         , ustar        ,                                           &
         632           0 :                         kvh_in       , kvm_in       , kvh         , kvm        ,                &
         633           0 :                         tpert        , qpert        , qrlin       , kvf        , tke          , &
         634           0 :                         wstarent     , bprod        , sprod       , minpblh    , wpert        , &
         635           0 :                         tkes         , went         , turbtype    ,                             &
         636           0 :                         kbase_o      , ktop_o       , ncvfin_o    ,                             &
         637           0 :                         kbase_mg     , ktop_mg      , ncvfin_mg   ,                             &
         638           0 :                         kbase_f      , ktop_f       , ncvfin_f    ,                             &
         639           0 :                         wet_CL       , web_CL       , jtbu_CL     , jbbu_CL    ,                &
         640           0 :                         evhc_CL      , jt2slv_CL    , n2ht_CL     , n2hb_CL    , lwp_CL       , &
         641           0 :                         opt_depth_CL , radinvfrac_CL, radf_CL     , wstar_CL   , wstar3fact_CL, &
         642           0 :                         ebrk         , wbrk         , lbrk        , ricl       , ghcl         , &
         643           0 :                         shcl         , smcl         ,                                           &
         644           0 :                         gh_a         , sh_a         , sm_a        , ri_a       , leng         , &
         645           0 :                         wcap         , pblhp        , cld         , ipbl       , kpblh        , &
         646           0 :                         wsedl        , wsed_CL      , warnstring  , errstring)
         647             :
         648             :     !--------------------------------------------------------------------------------- !
         649             :     !                                                                                  !
         650             :     ! Purpose : This is a driver routine to compute eddy diffusion coefficients        !
         651             :     !           for heat (sl), momentum (u, v), moisture (qt), and other  trace        !
         652             :     !           constituents.   This scheme uses first order closure for stable        !
         653             :     !           turbulent layers (STL). For convective layers (CL), entrainment        !
         654             :     !           closure is used at the CL external interfaces, which is coupled        !
         655             :     !           to the diagnosis of a CL regime mean TKE from the instantaneous        !
         656             :     !           thermodynamic and velocity profiles.   The CLs are diagnosed by        !
         657             :     !           extending original CL layers of moist static instability   into        !
         658             :     !           adjacent weakly stably stratified interfaces,   stopping if the        !
         659             :     !           stability is too strong.   This allows a realistic depiction of        !
         660             :     !           dry convective boundary layers with a downgradient approach.           !
         661             :     !                                                                                  !
         662             :     ! NOTE:     This routine currently assumes ntop_turb = 1, nbot_turb = pver         !
         663             :     !           ( turbulent diffusivities computed at all interior interfaces )        !
         664             :     !           and will require modification to handle a different ntop_turb.         !
         665             :     !                                                                                  !
         666             :     ! Authors:  Sungsu Park and Chris Bretherton. 08/2006, 05/2008.                    !
         667             :     !                                                                                  !
         668             :     ! For details, see                                                                 !
         669             :     !                                                                                  !
         670             :     ! 1. 'A new moist turbulence parametrization in the Community Atmosphere Model'    !
         671             :     !     by Christopher S. Bretherton & Sungsu Park. J. Climate. 22. 3422-3448. 2009. !
         672             :     !                                                                                  !
         673             :     ! 2. 'The University of Washington shallow convection and moist turbulence schemes !
         674             :     !     and their impact on climate simulations with the Community Atmosphere Model' !
         675             :     !     by Sungsu Park & Christopher S. Bretherton. J. Climate. 22. 3449-3469. 2009. !
         676             :     !                                                                                  !
         677             :     ! For questions on the scheme and code, send an email to                           !
         678             :     !     sungsup@ucar.edu or breth@washington.edu                                     !
         679             :     !                                                                                  !
         680             :     !--------------------------------------------------------------------------------- !
         681             :
         682             :     use pbl_utils, only: &
         683             :       compute_radf      ! Subroutine for computing radf
         684             :
         685             :     ! ---------------- !
         686             :     ! Inputs variables !
         687             :     ! ---------------- !
         688             :
         689             :     implicit none
         690             :     integer,  intent(in) :: pcols                     ! Number of atmospheric columns
         691             :     integer,  intent(in) :: pver                      ! Number of atmospheric layers
         692             :     integer,  intent(in) :: ncol                      ! Number of atmospheric columns
         693             :     real(r8), intent(in) :: u(pcols,pver)             ! U wind [ m/s ]
         694             :     real(r8), intent(in) :: v(pcols,pver)             ! V wind [ m/s ]
         695             :     real(r8), intent(in) :: sl(pcols,pver)            ! Liquid water static energy, cp * T + g * z - Lv * ql - Ls * qi [ J/kg ]
         696             :     real(r8), intent(in) :: slv(pcols,pver)           ! Liquid water virtual static energy, sl * ( 1 + 0.608 * qt ) [ J/kg ]
         697             :     real(r8), intent(in) :: qt(pcols,pver)            ! Total speccific humidity  qv + ql + qi [ kg/kg ]
         698             :     real(r8), intent(in) :: ql(pcols,pver)            ! Liquid water specific humidity [ kg/kg ]
         699             :     real(r8), intent(in) :: pi(pcols,pver+1)          ! Interface pressures [ Pa ]
         700             :     real(r8), intent(in) :: z(pcols,pver)             ! Layer midpoint height above surface [ m ]
         701             :     real(r8), intent(in) :: zi(pcols,pver+1)          ! Interface height above surface, i.e., zi(pver+1) = 0 all over the globe
         702             :                                                       ! [ m ]
         703             :     real(r8), intent(in) :: chu(pcols,pver+1)         ! Buoyancy coeffi. unsaturated sl (heat) coef. at all interfaces.
         704             :                                                       ! [ unit ? ]
         705             :     real(r8), intent(in) :: chs(pcols,pver+1)         ! Buoyancy coeffi. saturated sl (heat) coef. at all interfaces.
         706             :                                                       ! [ unit ? ]
         707             :     real(r8), intent(in) :: cmu(pcols,pver+1)         ! Buoyancy coeffi. unsaturated qt (moisture) coef. at all interfaces
         708             :                                                       ! [ unit ? ]
         709             :     real(r8), intent(in) :: cms(pcols,pver+1)         ! Buoyancy coeffi. saturated qt (moisture) coef. at all interfaces
         710             :                                                       ! [ unit ? ]
         711             :     real(r8), intent(in) :: sfuh(pcols,pver)          ! Saturation fraction in upper half-layer [ fraction ]
         712             :     real(r8), intent(in) :: sflh(pcols,pver)          ! Saturation fraction in lower half-layer [ fraction ]
         713             :     real(r8), intent(in) :: n2(pcols,pver)            ! Interfacial (except surface) moist buoyancy frequency [ s-2 ]
         714             :     real(r8), intent(in) :: s2(pcols,pver)            ! Interfacial (except surface) shear frequency [ s-2 ]
         715             :     real(r8), intent(in) :: ri(pcols,pver)            ! Interfacial (except surface) Richardson number
         716             :     real(r8), intent(in) :: qflx(pcols)               ! Kinematic surface constituent ( water vapor ) flux [ kg/m2/s ]
         717             :     real(r8), intent(in) :: shflx(pcols)              ! Kinematic surface heat flux [ unit ? ]
         718             :     real(r8), intent(in) :: slslope(pcols,pver)       ! Slope of 'sl' in each layer [ J/kg/Pa ]
         719             :     real(r8), intent(in) :: qtslope(pcols,pver)       ! Slope of 'qt' in each layer [ kg/kg/Pa ]
         720             :     real(r8), intent(in) :: qrlin(pcols,pver)         ! Input grid-mean LW heating rate : [ K/s ] * cpair * dp = [ W/kg*Pa ]
         721             :     real(r8), intent(in) :: wsedl(pcols,pver)         ! Sedimentation velocity of liquid stratus cloud droplet [ m/s ]
         722             :     real(r8), intent(in) :: ustar(pcols)              ! Surface friction velocity [ m/s ]
         723             :     real(r8), intent(in) :: rrho(pcols)               ! 1./bottom mid-point density. Specific volume [ m3/kg ]
         724             :     real(r8), intent(in) :: kvf(pcols,pver+1)         ! Free atmosphere eddy diffusivity [ m2/s ]
         725             :     logical,  intent(in) :: wstarent                  ! Switch for choosing wstar3 entrainment parameterization
         726             :     real(r8), intent(in) :: minpblh(pcols)            ! Minimum PBL height based on surface stress [ m ]
         727             :     real(r8), intent(in) :: kvh_in(pcols,pver+1)      ! kvh saved from last timestep or last iterative step [ m2/s ]
         728             :     real(r8), intent(in) :: kvm_in(pcols,pver+1)      ! kvm saved from last timestep or last iterative step [ m2/s ]
         729             :     real(r8), intent(in) :: cld(pcols,pver)           ! Stratus Cloud Fraction [ fraction ]
         730             :
         731             :     ! ---------------- !
         732             :     ! Output variables !
         733             :     ! ---------------- !
         734             :
         735             :     real(r8), intent(out) :: kvh(pcols,pver+1)        ! Eddy diffusivity for heat, moisture, and tracers [ m2/s ]
         736             :     real(r8), intent(out) :: kvm(pcols,pver+1)        ! Eddy diffusivity for momentum [ m2/s ]
         737             :     real(r8), intent(out) :: pblh(pcols)              ! PBL top height [ m ]
         738             :     real(r8), intent(out) :: pblhp(pcols)             ! PBL top height pressure [ Pa ]
         739             :     real(r8), intent(out) :: tpert(pcols)             ! Convective temperature excess [ K ]
         740             :     real(r8), intent(out) :: qpert(pcols)             ! Convective humidity excess [ kg/kg ]
         741             :     real(r8), intent(out) :: wpert(pcols)             ! Turbulent velocity excess [ m/s ]
         742             :     real(r8), intent(out) :: tkes(pcols)              ! TKE at surface [ m2/s2 ]
         743             :     real(r8), intent(out) :: went(pcols)              ! Entrainment rate at the PBL top interface [ m/s ]
         744             :     real(r8), intent(out) :: tke(pcols,pver+1)        ! Turbulent kinetic energy [ m2/s2 ], 'tkes' at surface, pver+1.
         745             :     real(r8), intent(out) :: bprod(pcols,pver+1)      ! Buoyancy production [ m2/s3 ],     'bflxs' at surface, pver+1.
         746             :     real(r8), intent(out) :: sprod(pcols,pver+1)      ! Shear production [ m2/s3 ], (ustar(i)**3)/(vk*z(i,pver))
         747             :                                                       ! at surface, pver+1.
         748             :     integer(i4), intent(out) :: turbtype(pcols,pver+1) ! Turbulence type at each interface:
         749             :                                                       ! 0. = Non turbulence interface
         750             :                                                       ! 1. = Stable turbulence interface
         751             :                                                       ! 2. = CL interior interface ( if bflxs > 0, surface is this )
         752             :                                                       ! 3. = Bottom external interface of CL
         753             :                                                       ! 4. = Top external interface of CL.
         754             :                                                       ! 5. = Double entraining CL external interface
         755             :     integer(i4), intent(out) :: ipbl(pcols)           ! If 1, PBL is CL, while if 0, PBL is STL.
         756             :     integer(i4), intent(out) :: kpblh(pcols)          ! Layer index containing PBL within or at the base interface
         757             :     real(r8), intent(out) :: wsed_CL(pcols,ncvmax)    ! Sedimentation velocity at the top of each CL [ m/s ]
         758             :
         759             :     character(len=*), intent(out) :: warnstring
         760             :     character(len=*), intent(out) :: errstring
         761             :
         762             :     ! --------------------------- !
         763             :     ! Diagnostic output variables !
         764             :     ! --------------------------- !
         765             :
         766             :     real(r8) :: kbase_o(pcols,ncvmax)                 ! Original external base interface index of CL just after 'exacol'
         767             :     real(r8) :: ktop_o(pcols,ncvmax)                  ! Original external top  interface index of CL just after 'exacol'
         768             :     real(r8) :: ncvfin_o(pcols)                       ! Original number of CLs just after 'exacol'
         769             :     real(r8) :: kbase_mg(pcols,ncvmax)                ! kbase  just after extending-merging (after 'zisocl') but without SRCL
         770             :     real(r8) :: ktop_mg(pcols,ncvmax)                 ! ktop   just after extending-merging (after 'zisocl') but without SRCL
         771             :     real(r8) :: ncvfin_mg(pcols)                      ! ncvfin just after extending-merging (after 'zisocl') but without SRCL
         772             :     real(r8) :: kbase_f(pcols,ncvmax)                 ! Final kbase  after adding SRCL
         773             :     real(r8) :: ktop_f(pcols,ncvmax)                  ! Final ktop   after adding SRCL
         774             :     real(r8) :: ncvfin_f(pcols)                       ! Final ncvfin after adding SRCL
         775             :     real(r8) :: wet_CL(pcols,ncvmax)                  ! Entrainment rate at the CL top [ m/s ]
         776             :     real(r8) :: web_CL(pcols,ncvmax)                  ! Entrainment rate at the CL base [ m/s ]
         777             :     real(r8) :: jtbu_CL(pcols,ncvmax)                 ! Buoyancy jump across the CL top [ m/s2 ]
         778             :     real(r8) :: jbbu_CL(pcols,ncvmax)                 ! Buoyancy jump across the CL base [ m/s2 ]
         779             :     real(r8) :: evhc_CL(pcols,ncvmax)                 ! Evaporative enhancement factor at the CL top
         780             :     real(r8) :: jt2slv_CL(pcols,ncvmax)               ! Jump of slv ( across two layers ) at CL top for use only in evhc [ J/kg ]
         781             :     real(r8) :: n2ht_CL(pcols,ncvmax)                 ! n2 defined at the CL top  interface
         782             :                                                       ! but using sfuh(kt)   instead of sfi(kt) [ s-2 ]
         783             :     real(r8) :: n2hb_CL(pcols,ncvmax)                 ! n2 defined at the CL base interface
         784             :                                                       ! but using sflh(kb-1) instead of sfi(kb) [ s-2 ]
         785             :     real(r8) :: lwp_CL(pcols,ncvmax)                  ! LWP in the CL top layer [ kg/m2 ]
         786             :     real(r8) :: opt_depth_CL(pcols,ncvmax)            ! Optical depth of the CL top layer
         787             :     real(r8) :: radinvfrac_CL(pcols,ncvmax)           ! Fraction of LW radiative cooling confined in the top portion of CL
         788             :     real(r8) :: radf_CL(pcols,ncvmax)                 ! Buoyancy production at the CL top due to radiative cooling [ m2/s3 ]
         789             :     real(r8) :: wstar_CL(pcols,ncvmax)                ! Convective velocity of CL including entrainment contribution finally [ m/s ]
         790             :     real(r8) :: wstar3fact_CL(pcols,ncvmax)           ! "wstar3fact" of CL. Entrainment enhancement of wstar3 (inverse)
         791             :
         792             :     real(r8) :: gh_a(pcols,pver+1)                    ! Half of normalized buoyancy production, -l2n2/2e. [ no unit ]
         793             :     real(r8) :: sh_a(pcols,pver+1)                    ! Galperin instability function of heat-moisture at all interfaces [ no unit ]
         794             :     real(r8) :: sm_a(pcols,pver+1)                    ! Galperin instability function of momentum      at all interfaces [ no unit ]
         795             :     real(r8) :: ri_a(pcols,pver+1)                    ! Interfacial Richardson number                  at all interfaces [ no unit ]
         796             :
         797             :     real(r8) :: ebrk(pcols,ncvmax)                    ! Net CL mean TKE [ m2/s2 ]
         798             :     real(r8) :: wbrk(pcols,ncvmax)                    ! Net CL mean normalized TKE [ m2/s2 ]
         799             :     real(r8) :: lbrk(pcols,ncvmax)                    ! Net energetic integral thickness of CL [ m ]
         800             :     real(r8) :: ricl(pcols,ncvmax)                    ! Mean Richardson number of CL ( l2n2/l2s2 )
         801             :     real(r8) :: ghcl(pcols,ncvmax)                    ! Half of normalized buoyancy production of CL
         802             :     real(r8) :: shcl(pcols,ncvmax)                    ! Instability function of heat and moisture of CL
         803             :     real(r8) :: smcl(pcols,ncvmax)                    ! Instability function of momentum of CL
         804             :
         805             :     real(r8) :: leng(pcols,pver+1)                    ! Turbulent length scale [ m ], 0 at the surface.
         806             :     real(r8) :: wcap(pcols,pver+1)                    ! Normalized TKE [m2/s2], 'tkes/b1' at the surface and 'tke/b1' at
         807             :                                                       ! the top/bottom entrainment interfaces of CL assuming no transport.
         808             :     ! ------------------------ !
         809             :     ! Local Internal Variables !
         810             :     ! ------------------------ !
         811             :
         812           0 :     logical :: belongcv(pcols,pver+1)                 ! True for interfaces in a CL (both interior and exterior are included)
         813           0 :     logical :: belongst(pcols,pver+1)                 ! True for stable turbulent layer interfaces (STL)
         814             :     logical :: in_CL                                  ! True if interfaces k,k+1 both in same CL.
         815             :     logical :: extend                                 ! True when CL is extended in zisocl
         816             :     logical :: extend_up                              ! True when CL is extended upward in zisocl
         817             :     logical :: extend_dn                              ! True when CL is extended downward in zisocl
         818             :
         819             :     integer :: i                                      ! Longitude index
         820             :     integer :: k                                      ! Vertical index
         821             :     integer :: ks                                     ! Vertical index
         822           0 :     integer :: ncvfin(pcols)                          ! Total number of CL in column
         823             :     integer :: ncvf                                   ! Total number of CL in column prior to adding SRCL
         824             :     integer :: ncv                                    ! Index of current CL
         825             :     integer :: ncvnew                                 ! Index of added SRCL appended after regular CLs from 'zisocl'
         826             :     integer :: ncvsurf                                ! If nonzero, CL index based on surface
         827             :                                                       ! (usually 1, but can be > 1 when SRCL is based at sfc)
         828           0 :     integer :: kbase(pcols,ncvmax)                    ! Vertical index of CL base interface
         829           0 :     integer :: ktop(pcols,ncvmax)                     ! Vertical index of CL top interface
         830             :     integer :: kb, kt                                 ! kbase and ktop for current CL
         831             :     integer :: ktblw                                  ! ktop of the CL located at just below the current CL
         832             :
         833           0 :     integer  :: ktopbl(pcols)                         ! PBL top height or interface index
         834           0 :     real(r8) :: bflxs(pcols)                          ! Surface buoyancy flux [ m2/s3 ]
         835             :     real(r8) :: rcap                                  ! 'tke/ebrk' at all interfaces of CL.
         836             :                                                       ! Set to 1 at the CL entrainment interfaces
         837             :     real(r8) :: jtzm                                  ! Interface layer thickness of CL top interface [ m ]
         838             :     real(r8) :: jtsl                                  ! Jump of s_l across CL top interface [ J/kg ]
         839             :     real(r8) :: jtqt                                  ! Jump of q_t across CL top interface [ kg/kg ]
         840             :     real(r8) :: jtbu                                  ! Jump of buoyancy across CL top interface [ m/s2 ]
         841             :     real(r8) :: jtu                                   ! Jump of u across CL top interface [ m/s ]
         842             :     real(r8) :: jtv                                   ! Jump of v across CL top interface [ m/s ]
         843             :     real(r8) :: jt2slv                                ! Jump of slv ( across two layers ) at CL top for use only in evhc [ J/kg ]
         844             :     real(r8) :: radf                                  ! Buoyancy production at the CL top due to radiative cooling [ m2/s3 ]
         845             :     real(r8) :: jbzm                                  ! Interface layer thickness of CL base interface [ m ]
         846             :     real(r8) :: jbsl                                  ! Jump of s_l across CL base interface [ J/kg ]
         847             :     real(r8) :: jbqt                                  ! Jump of q_t across CL top interface [ kg/kg ]
         848             :     real(r8) :: jbbu                                  ! Jump of buoyancy across CL base interface [ m/s2 ]
         849             :     real(r8) :: jbu                                   ! Jump of u across CL base interface [ m/s ]
         850             :     real(r8) :: jbv                                   ! Jump of v across CL base interface [ m/s ]
         851             :     real(r8) :: ch                                    ! Buoyancy coefficients defined at the CL top and base interfaces
         852             :                                                       ! using CL internal
         853             :     real(r8) :: cm                                    ! sfuh(kt) and sflh(kb-1) instead of sfi(kt) and sfi(kb), respectively.
         854             :                                                       ! These are used for entrainment calculation at CL external interfaces
         855             :                                                       ! and SRCL identification.
         856             :     real(r8) :: n2ht                                  ! n2 defined at the CL top  interface
         857             :                                                       ! but using sfuh(kt)   instead of sfi(kt) [ s-2 ]
         858             :     real(r8) :: n2hb                                  ! n2 defined at the CL base interface
         859             :                                                       ! but using sflh(kb-1) instead of sfi(kb) [ s-2 ]
         860             :     real(r8) :: n2htSRCL                              ! n2 defined at the upper-half layer of SRCL.
         861             :                                                       ! This is used only for identifying SRCL.
         862             :                                                       ! n2htSRCL use SRCL internal slope sl and qt
         863             :                                                       ! as well as sfuh(kt) instead of sfi(kt) [ s-2 ]
         864             :     real(r8) :: gh                                    ! Half of normalized buoyancy production ( -l2n2/2e ) [ no unit ]
         865             :     real(r8) :: sh                                    ! Galperin instability function for heat and moisture
         866             :     real(r8) :: sm                                    ! Galperin instability function for momentum
         867             :     real(r8) :: lbulk                                 ! Depth of turbulent layer, Master length scale (not energetic length)
         868             :     real(r8) :: dzht                                  ! Thickness of top    half-layer [ m ]
         869             :     real(r8) :: dzhb                                  ! Thickness of bottom half-layer [ m ]
         870             :     real(r8) :: rootp                                 ! Sqrt(net CL-mean TKE including entrainment contribution) [ m/s ]
         871             :     real(r8) :: evhc                                  ! Evaporative enhancement factor: (1+E)
         872             :                                                       ! with E = evap. cool. efficiency [ no unit ]
         873             :     real(r8) :: kentr                                 ! Effective entrainment diffusivity 'wet*dz', 'web*dz' [ m2/s ]
         874             :     real(r8) :: lwp                                   ! Liquid water path in the layer kt [ kg/m2 ]
         875             :     real(r8) :: opt_depth                             ! Optical depth of the layer kt [ no unit ]
         876             :     real(r8) :: radinvfrac                            ! Fraction of LW cooling in the layer kt
         877             :                                                       ! concentrated at the CL top [ no unit ]
         878             :     real(r8) :: wet                                   ! CL top entrainment rate [ m/s ]
         879             :     real(r8) :: web                                   ! CL bot entrainment rate [ m/s ]. Set to zero if CL is based at surface.
         880             :     real(r8) :: vyt                                   ! n2ht/n2 at the CL top  interface
         881             :     real(r8) :: vyb                                   ! n2hb/n2 at the CL base interface
         882             :     real(r8) :: vut                                   ! Inverse Ri (=s2/n2) at the CL top  interface
         883             :     real(r8) :: vub                                   ! Inverse Ri (=s2/n2) at the CL base interface
         884             :     real(r8) :: fact                                  ! Factor relating TKE generation to entrainment [ no unit ]
         885             :     real(r8) :: trma                                  ! Intermediate variables used for solving quadratic ( for gh from ri )
         886             :     real(r8) :: trmb                                  ! and cubic equations ( for ebrk: the net CL mean TKE )
         887             :     real(r8) :: trmc                                  !
         888             :     real(r8) :: trmp                                  !
         889             :     real(r8) :: trmq                                  !
         890             :     real(r8) :: qq                                    !
         891             :     real(r8) :: det                                   !
         892             :     real(r8) :: gg                                    ! Intermediate variable used for calculating stability functions of
         893             :                                                       ! SRCL or SBCL based at the surface with bflxs > 0.
         894             :     real(r8) :: dzhb5                                 ! Half thickness of the bottom-most layer of current CL regime
         895             :     real(r8) :: dzht5                                 ! Half thickness of the top-most layer of adjacent CL regime
         896             :                                                       ! just below current CL
         897           0 :     real(r8) :: qrlw(pcols,pver)                      ! Local grid-mean LW heating rate : [K/s] * cpair * dp = [ W/kg*Pa ]
         898             :
         899           0 :     real(r8) :: cldeff(pcols,pver)                    ! Effective stratus fraction
         900             :     real(r8) :: qleff                                 ! Used for computing evhc
         901             :     real(r8) :: tunlramp                              ! Ramping tunl
         902             :     real(r8) :: leng_imsi                             ! For Kv = max(Kv_STL, Kv_entrain)
         903             :     real(r8) :: tke_imsi                              !
         904             :     real(r8) :: kvh_imsi                              !
         905             :     real(r8) :: kvm_imsi                              !
         906             :     real(r8) :: alph4exs                              ! For extended stability function in the stable regime
         907             :     real(r8) :: ghmin                                 !
         908             :
         909             :     real(r8) :: sedfact                               ! For 'sedimentation-entrainment feedback'
         910             :
         911             :     ! Local variables specific for 'wstar' entrainment closure
         912             :
         913             :     real(r8) :: cet                                   ! Proportionality coefficient between wet and wstar3
         914             :     real(r8) :: ceb                                   ! Proportionality coefficient between web and wstar3
         915             :     real(r8) :: wstar                                 ! Convective velocity for CL [ m/s ]
         916             :     real(r8) :: wstar3                                ! Cubed convective velocity for CL [ m3/s3 ]
         917             :     real(r8) :: wstar3fact                            ! 1/(relative change of wstar^3 by entrainment)
         918             :     real(r8) :: rmin                                  ! sqrt(p)
         919             :     real(r8) :: fmin                                  ! f(rmin), where f(r) = r^3 - 3*p*r - 2q
         920             :     real(r8) :: rcrit                                 ! ccrit*wstar
         921             :     real(r8) :: fcrit                                 ! f(rcrit)
         922             :     logical     noroot                                ! True if f(r) has no root r > rcrit
         923             :
         924             :     character(128) :: temp_string
         925             :
         926             :     !-----------------------!
         927             :     ! Start of Main Program !
         928             :     !-----------------------!
         929             :
         930           0 :     warnstring = ""
         931           0 :     errstring = ""
         932             :
         933             :     ! Option: Turn-off LW radiative-turbulence interaction in PBL scheme
         934             :     !         by setting qrlw = 0.  Logical parameter 'set_qrlzero'  was
         935             :     !         defined in the first part of 'eddy_diff.F90' module.
         936             :
         937             :     if( set_qrlzero ) then
         938             :         qrlw(:,:) = 0._r8
         939             :     else
         940           0 :         qrlw(:ncol,:pver) = qrlin(:ncol,:pver)
         941             :     endif
         942             :
         943             :     ! Define effective stratus fraction using the grid-mean ql.
         944             :     ! Modification : The contribution of ice should be carefully considered.
         945             :     !                This should be done in combination with the 'qrlw' and
         946             :     !                overlapping assumption of liquid and ice stratus.
         947             :
         948           0 :     do k = 1, pver
         949           0 :        do i = 1, ncol
         950           0 :           if( choice_evhc .eq. 'ramp' .or. choice_radf .eq. 'ramp' ) then
         951             :               cldeff(i,k) = cld(i,k) * min( ql(i,k) / qmin, 1._r8 )
         952             :           else
         953           0 :               cldeff(i,k) = cld(i,k)
         954             :           endif
         955             :        end do
         956             :     end do
         957             :
         958             :     ! For an extended stability function in the stable regime, re-define
         959             :     ! alph4exe and ghmin. This is for future work.
         960             :
         961             :     if( ricrit .eq. 0.19_r8 ) then
         962           0 :         alph4exs = alph4
         963           0 :         ghmin    = -3.5334_r8
         964             :     elseif( ricrit .gt. 0.19_r8 ) then
         965             :         alph4exs = -2._r8 * b1 * alph2 / ( alph3 - 2._r8 * b1 * alph5 ) / ricrit
         966             :         ghmin    = -1.e10_r8
         967             :     else
         968             :        errstring = 'ricrit should be larger than 0.19 in UW PBL'
         969             :        return
         970             :     endif
         971             :
         972             :     !
         973             :     ! Initialization of Diagnostic Output
         974             :     !
         975             :
         976           0 :     do i = 1, ncol
         977           0 :        went(i)                  = 0._r8
         978           0 :        wet_CL(i,:ncvmax)        = 0._r8
         979           0 :        web_CL(i,:ncvmax)        = 0._r8
         980           0 :        jtbu_CL(i,:ncvmax)       = 0._r8
         981           0 :        jbbu_CL(i,:ncvmax)       = 0._r8
         982           0 :        evhc_CL(i,:ncvmax)       = 0._r8
         983           0 :        jt2slv_CL(i,:ncvmax)     = 0._r8
         984           0 :        n2ht_CL(i,:ncvmax)       = 0._r8
         985           0 :        n2hb_CL(i,:ncvmax)       = 0._r8
         986           0 :        lwp_CL(i,:ncvmax)        = 0._r8
         987           0 :        opt_depth_CL(i,:ncvmax)  = 0._r8
         988           0 :        radinvfrac_CL(i,:ncvmax) = 0._r8
         989           0 :        radf_CL(i,:ncvmax)       = 0._r8
         990           0 :        wstar_CL(i,:ncvmax)      = 0._r8
         991           0 :        wstar3fact_CL(i,:ncvmax) = 0._r8
         992           0 :        ricl(i,:ncvmax)          = 0._r8
         993           0 :        ghcl(i,:ncvmax)          = 0._r8
         994           0 :        shcl(i,:ncvmax)          = 0._r8
         995           0 :        smcl(i,:ncvmax)          = 0._r8
         996           0 :        ebrk(i,:ncvmax)          = 0._r8
         997           0 :        wbrk(i,:ncvmax)          = 0._r8
         998           0 :        lbrk(i,:ncvmax)          = 0._r8
         999           0 :        gh_a(i,:pver+1)          = 0._r8
        1000           0 :        sh_a(i,:pver+1)          = 0._r8
        1001           0 :        sm_a(i,:pver+1)          = 0._r8
        1002           0 :        ri_a(i,:pver+1)          = 0._r8
        1003           0 :        ipbl(i)                  = 0
        1004           0 :        kpblh(i)                 = pver
        1005           0 :        wsed_CL(i,:ncvmax)       = 0._r8
        1006             :     end do
        1007             :
        1008             :     ! kvh and kvm are stored over timesteps in 'vertical_diffusion.F90' and
        1009             :     ! passed in as kvh_in and kvm_in.  However,  at the first timestep they
        1010             :     ! need to be computed and these are done just before calling 'caleddy'.
        1011             :     ! kvm and kvh are also stored over iterative time step in the first part
        1012             :     ! of 'eddy_diff.F90'
        1013             :
        1014             :     ! Initialize kvh and kvm to kvf
        1015           0 :     kvh(:,:) = kvf(:,:)
        1016           0 :     kvm(:,:) = kvf(:,:)
        1017             :     ! Zero diagnostic quantities for the new diffusion step.
        1018           0 :     wcap(:,:) = 0._r8
        1019           0 :     leng(:,:) = 0._r8
        1020           0 :     tke(:,:)  = 0._r8
        1021           0 :     turbtype(:,:) = 0
        1022             :
        1023             :
        1024             :     ! Initialize 'bprod' [ m2/s3 ] and 'sprod' [ m2/s3 ] at all interfaces.
        1025             :     ! Note this initialization is a hybrid initialization since 'n2' [s-2] and 's2' [s-2]
        1026             :     ! are calculated from the given current initial profile, while 'kvh_in' [m2/s] and
        1027             :     ! 'kvm_in' [m2/s] are from the previous iteration or previous time step.
        1028             :     ! This initially guessed 'bprod' and 'sprod' will be updated at the end of this
        1029             :     ! 'caleddy' subroutine for diagnostic output.
        1030             :     ! This computation of 'brpod,sprod' below is necessary for wstar-based entrainment closure.
        1031             :
        1032           0 :     do k = 2, pver
        1033           0 :        do i = 1, ncol
        1034           0 :             bprod(i,k) = -kvh_in(i,k) * n2(i,k)
        1035           0 :             sprod(i,k) =  kvm_in(i,k) * s2(i,k)
        1036             :        end do
        1037             :     end do
        1038             :
        1039             :     ! Set 'bprod' and 'sprod' at top and bottom interface.
        1040             :     ! In calculating 'surface' (actually lowest half-layer) buoyancy flux,
        1041             :     ! 'chu' at surface is defined to be the same as 'chu' at the mid-point
        1042             :     ! of lowest model layer (pver) at the end of 'trbind'. The same is for
        1043             :     ! the other buoyancy coefficients.  'sprod(i,pver+1)'  is defined in a
        1044             :     ! consistent way as the definition of 'tkes' in the original code.
        1045             :     ! ( Important Option ) If I want to isolate surface buoyancy flux from
        1046             :     ! the other parts of CL regimes energetically even though bflxs > 0,
        1047             :     ! all I should do is to re-define 'bprod(i,pver+1)=0' in the below 'do'
        1048             :     ! block. Additionally for merging test of extending SBCL based on 'l2n2'
        1049             :     ! in 'zisocl', I should use 'l2n2 = - wint / sh'  for similar treatment
        1050             :     ! as previous code. All other parts of the code  are fully consistently
        1051             :     ! treated by these change only.
        1052             :     ! My future general convection scheme will use bflxs(i).
        1053             :
        1054           0 :     do i = 1, ncol
        1055           0 :        bprod(i,1) = 0._r8 ! Top interface
        1056           0 :        sprod(i,1) = 0._r8 ! Top interface
        1057           0 :        ch = chu(i,pver+1) * ( 1._r8 - sflh(i,pver) ) + chs(i,pver+1) * sflh(i,pver)
        1058           0 :        cm = cmu(i,pver+1) * ( 1._r8 - sflh(i,pver) ) + cms(i,pver+1) * sflh(i,pver)
        1059           0 :        bflxs(i) = ch * shflx(i) * rrho(i) + cm * qflx(i) * rrho(i)
        1060             :        if( choice_tkes .eq. 'ibprod' ) then
        1061           0 :            bprod(i,pver+1) = bflxs(i)
        1062             :        else
        1063             :            bprod(i,pver+1) = 0._r8
        1064             :        endif
        1065           0 :        sprod(i,pver+1) = (ustar(i)**3)/(vk*z(i,pver))
        1066             :     end do
        1067             :
        1068             :     ! Initially identify CL regimes in 'exacol'
        1069             :     !    ktop  : Interface index of the CL top  external interface
        1070             :     !    kbase : Interface index of the CL base external interface
        1071             :     !    ncvfin: Number of total CLs
        1072             :     ! Note that if surface buoyancy flux is positive ( bflxs = bprod(i,pver+1) > 0 ),
        1073             :     ! surface interface is identified as an internal interface of CL. However, even
        1074             :     ! though bflxs <= 0, if 'pver' interface is a CL internal interface (ri(pver)<0),
        1075             :     ! surface interface is identified as an external interface of CL. If bflxs =< 0
        1076             :     ! and ri(pver) >= 0, then surface interface is identified as a stable turbulent
        1077             :     ! intereface (STL) as shown at the end of 'caleddy'. Even though a 'minpblh' is
        1078             :     ! passed into 'exacol', it is not used in the 'exacol'.
        1079             :
        1080           0 :     call exacol( pcols, pver, ncol, ri, bflxs, minpblh, zi, ktop, kbase, ncvfin )
        1081             :
        1082             :     ! Diagnostic output of CL interface indices before performing 'extending-merging'
        1083             :     ! of CL regimes in 'zisocl'
        1084           0 :     do i = 1, ncol
        1085           0 :     do k = 1, ncvmax
        1086           0 :        kbase_o(i,k) = real(kbase(i,k),r8)
        1087           0 :        ktop_o(i,k)  = real(ktop(i,k),r8)
        1088           0 :        ncvfin_o(i)  = real(ncvfin(i),r8)
        1089             :     end do
        1090             :     end do
        1091             :
        1092             :     ! ----------------------------------- !
        1093             :     ! Perform calculation for each column !
        1094             :     ! ----------------------------------- !
        1095             :
        1096           0 :     do i = 1, ncol
        1097             :
        1098             :        ! Define Surface Interfacial Layer TKE, 'tkes'.
        1099             :        ! In the current code, 'tkes' is used as representing TKE of surface interfacial
        1100             :        ! layer (low half-layer of surface-based grid layer). In the code, when bflxs>0,
        1101             :        ! surface interfacial layer is assumed to be energetically  coupled to the other
        1102             :        ! parts of the CL regime based at the surface. In this sense, it is conceptually
        1103             :        ! more reasonable to include both 'bprod' and 'sprod' in the definition of 'tkes'.
        1104             :        ! Since 'tkes' cannot be negative, it is lower bounded by small positive number.
        1105             :        ! Note that inclusion of 'bprod' in the definition of 'tkes' may increase 'ebrk'
        1106             :        ! and 'wstar3', and eventually, 'wet' at the CL top, especially when 'bflxs>0'.
        1107             :        ! This might help to solve the problem of too shallow PBLH over the overcast Sc
        1108             :        ! regime. If I want to exclude 'bprod(i,pver+1)' in calculating 'tkes' even when
        1109             :        ! bflxs > 0, all I should to do is to set 'bprod(i,pver+1) = 0' in the above
        1110             :        ! initialization 'do' loop (explained above), NOT changing the formulation of
        1111             :        ! tkes(i) in the below block. This is because for consistent treatment in the
        1112             :        ! other parts of the code also.
        1113             :
        1114             :      ! tkes(i) = (b1*vk*z(i,pver)*sprod(i,pver+1))**(2._r8/3._r8)
        1115           0 :        tkes(i) = max(b1*vk*z(i,pver)*(bprod(i,pver+1)+sprod(i,pver+1)), 1.e-7_r8)**(2._r8/3._r8)
        1116           0 :        tkes(i) = min(tkes(i), tkemax)
        1117           0 :        tke(i,pver+1)  = tkes(i)
        1118           0 :        wcap(i,pver+1) = tkes(i)/b1
        1119             :
        1120             :        ! Extend and merge the initially identified CLs, relabel the CLs, and calculate
        1121             :        ! CL internal mean energetics and stability functions in 'zisocl'.
        1122             :        ! The CL nearest to the surface is CL(1) and the CL index, ncv, increases
        1123             :        ! with height. The following outputs are from 'zisocl'. Here, the dimension
        1124             :        ! of below outputs are (pcols,ncvmax) (except the 'ncvfin(pcols)' and
        1125             :        ! 'belongcv(pcols,pver+1)) and 'ncv' goes from 1 to 'ncvfin'.
        1126             :        ! For 'ncv = ncvfin+1, ncvmax', below output are already initialized to be zero.
        1127             :        !      ncvfin       : Total number of CLs
        1128             :        !      kbase(ncv)   : Base external interface index of CL
        1129             :        !      ktop         : Top  external interface index of CL
        1130             :        !      belongcv     : True if the interface (either internal or external) is CL
        1131             :        !      ricl         : Mean Richardson number of internal CL
        1132             :        !      ghcl         : Normalized buoyancy production '-l2n2/2e' [no unit] of internal CL
        1133             :        !      shcl         : Galperin instability function of heat-moisture of internal CL
        1134             :        !      smcl         : Galperin instability function of momentum of internal CL
        1135             :        !      lbrk, <l>int : Thickness of (energetically) internal CL (lint, [m])
        1136             :        !      wbrk, <W>int : Mean normalized TKE of internal CL  ([m2/s2])
        1137             :        !      ebrk, <e>int : Mean TKE of internal CL (b1*wbrk,[m2/s2])
        1138             :        ! The ncvsurf is an identifier saying which CL regime is based at the surface.
        1139             :        ! If 'ncvsurf=1', then the first CL regime is based at the surface. If surface
        1140             :        ! interface is not a part of CL (neither internal nor external), 'ncvsurf = 0'.
        1141             :        ! After identifying and including SRCLs into the normal CL regimes (where newly
        1142             :        ! identified SRCLs are simply appended to the normal CL regimes using regime
        1143             :        ! indices of 'ncvfin+1','ncvfin+2' (as will be shown in the below SRCL part),..
        1144             :        ! where 'ncvfin' is the final CL regime index produced after extending-merging
        1145             :        ! in 'zisocl' but before adding SRCLs), if any newly identified SRCL (e.g.,
        1146             :        ! 'ncvfin+1') is based at surface, then 'ncvsurf = ncvfin+1'. Thus 'ncvsurf' can
        1147             :        ! be 0, 1, or >1. 'ncvsurf' can be a useful diagnostic output.
        1148             :
        1149           0 :        ncvsurf = 0
        1150           0 :        if( ncvfin(i) .gt. 0 ) then
        1151             :            call zisocl( pcols  , pver     , i        ,           &
        1152             :                         z      , zi       , n2       , s2      , &
        1153             :                         bprod  , sprod    , bflxs    , tkes    , &
        1154             :                         ncvfin , kbase    , ktop     , belongcv, &
        1155             :                         ricl   , ghcl     , shcl     , smcl    , &
        1156             :                         lbrk   , wbrk     , ebrk     ,           &
        1157             :                         extend , extend_up, extend_dn, &
        1158           0 :                         errstring)
        1159           0 :            if (trim(errstring) /= "") return
        1160           0 :            if( kbase(i,1) .eq. pver + 1 ) ncvsurf = 1
        1161             :        else
        1162           0 :            belongcv(i,:) = .false.
        1163             :        endif
        1164             :
        1165             :        ! Diagnostic output after finishing extending-merging process in 'zisocl'
        1166             :        ! Since we are adding SRCL additionally, we need to print out these here.
        1167             :
        1168           0 :        do k = 1, ncvmax
        1169           0 :           kbase_mg(i,k) = real(kbase(i,k))
        1170           0 :           ktop_mg(i,k)  = real(ktop(i,k))
        1171           0 :           ncvfin_mg(i)  = real(ncvfin(i))
        1172             :        end do
        1173             :
        1174             :        ! ----------------------- !
        1175             :        ! Identification of SRCLs !
        1176             :        ! ----------------------- !
        1177             :
        1178             :      ! Modification : This cannot identify the 'cirrus' layer due to the condition of
        1179             :      !                ql(i,k) .gt. qmin. This should be modified in future to identify
        1180             :      !                a single thin cirrus layer.
        1181             :      !                Instead of ql, we may use cldn in future, including ice
        1182             :      !                contribution.
        1183             :
        1184             :        ! ------------------------------------------------------------------------------ !
        1185             :        ! Find single-layer radiatively-driven cloud-topped convective layers (SRCLs).   !
        1186             :        ! SRCLs extend through a single model layer k, with entrainment at the top and   !
        1187             :        ! bottom interfaces, unless bottom interface is the surface.                     !
        1188             :        ! The conditions for an SRCL is identified are:                                  !
        1189             :        !                                                                                !
        1190             :        !   1. Cloud in the layer, k : ql(i,k) .gt. qmin = 1.e-5 [ kg/kg ]               !
        1191             :        !   2. No cloud in the above layer (else assuming that some fraction of the LW   !
        1192             :        !      flux divergence in layer k is concentrated at just below top interface    !
        1193             :        !      of layer k is invalid). Then, this condition might be sensitive to the    !
        1194             :        !      vertical resolution of grid.                                              !
        1195             :        !   3. LW radiative cooling (SW heating is assumed uniformly distributed through !
        1196             :        !      layer k, so not relevant to buoyancy production) in the layer k. However, !
        1197             :        !      SW production might also contribute, which may be considered in a future. !
        1198             :        !   4. Internal stratification 'n2ht' of upper-half layer should be unstable.    !
        1199             :        !      The 'n2ht' is pure internal stratification of upper half layer, obtained  !
        1200             :        !      using internal slopes of sl, qt in layer k (in contrast to conventional   !
        1201             :        !      interfacial slope) and saturation fraction in the upper-half layer,       !
        1202             :        !      sfuh(k) (in contrast to sfi(k)).                                          !
        1203             :        !   5. Top and bottom interfaces not both in the same existing convective layer. !
        1204             :        !      If SRCL is within the previouisly identified CL regimes, we don't define  !
        1205             :        !      a new SRCL.                                                               !
        1206             :        !   6. k >= ntop_turb + 1 = 2                                                    !
        1207             :        !   7. Ri at the top interface > ricrit = 0.19 (otherwise turbulent mixing will  !
        1208             :        !      broadly distribute the cloud top in the vertical, preventing localized    !
        1209             :        !      radiative destabilization at the top interface).                          !
        1210             :        !                                                                                !
        1211             :        ! Note if 'k = pver', it identifies a surface-based single fog layer, possibly,  !
        1212             :        ! warm advection fog. Note also the CL regime index of SRCLs itself increases    !
        1213             :        ! with height similar to the regular CLs indices identified from 'zisocl'.       !
        1214             :        ! ------------------------------------------------------------------------------ !
        1215             :
        1216           0 :        ncv  = 1
        1217           0 :        ncvf = ncvfin(i)
        1218             :
        1219             :        if( choice_SRCL .eq. 'remove' ) goto 222
        1220             :
        1221           0 :        do k = nbot_turb, ntop_turb + 1, -1 ! 'k = pver, 2, -1' is a layer index.
        1222             :
        1223           0 :           if( ql(i,k) .gt. qmin .and. ql(i,k-1) .lt. qmin .and. qrlw(i,k) .lt. 0._r8 &
        1224           0 :                                 .and. ri(i,k) .ge. ricrit ) then
        1225             :
        1226             :               ! In order to avoid any confliction with the treatment of ambiguous layer,
        1227             :               ! I need to impose an additional constraint that ambiguous layer cannot be
        1228             :               ! SRCL. So, I added constraint that 'k+1' interface (base interface of k
        1229             :               ! layer) should not be a part of previously identified CL. Since 'belongcv'
        1230             :               ! is even true for external entrainment interfaces, below constraint is
        1231             :               ! fully sufficient.
        1232             :
        1233           0 :               if( choice_SRCL .eq. 'nonamb' .and. belongcv(i,k+1) ) then
        1234             :                   go to 220
        1235             :               endif
        1236             :
        1237           0 :               ch = ( 1._r8 - sfuh(i,k) ) * chu(i,k) + sfuh(i,k) * chs(i,k)
        1238           0 :               cm = ( 1._r8 - sfuh(i,k) ) * cmu(i,k) + sfuh(i,k) * cms(i,k)
        1239             :
        1240           0 :               n2htSRCL = ch * slslope(i,k) + cm * qtslope(i,k)
        1241             :
        1242           0 :               if( n2htSRCL .le. 0._r8 ) then
        1243             :
        1244             :                   ! Test if bottom and top interfaces are part of the pre-existing CL.
        1245             :                   ! If not, find appropriate index for the new SRCL. Note that this
        1246             :                   ! calculation makes use of 'ncv set' obtained from 'zisocl'. The
        1247             :                   ! 'in_CL' is a parameter testing whether the new SRCL is already
        1248             :                   ! within the pre-existing CLs (.true.) or not (.false.).
        1249             :
        1250             :                   in_CL = .false.
        1251             :
        1252           0 :                   do while ( ncv .le. ncvf )
        1253           0 :                      if( ktop(i,ncv) .le. k ) then
        1254           0 :                         if( kbase(i,ncv) .gt. k ) then
        1255             :                             in_CL = .true.
        1256             :                         endif
        1257             :                         exit             ! Exit from 'do while' loop if SRCL is within the CLs.
        1258             :                      else
        1259           0 :                         ncv = ncv + 1    ! Go up one CL
        1260             :                      end if
        1261             :                   end do ! ncv
        1262             :
        1263             :                   if( .not. in_CL ) then ! SRCL is not within the pre-existing CLs.
        1264             :
        1265             :                      ! Identify a new SRCL and add it to the pre-existing CL regime group.
        1266             :
        1267           0 :                      ncvfin(i)       =  ncvfin(i) + 1
        1268           0 :                      ncvnew          =  ncvfin(i)
        1269           0 :                      ktop(i,ncvnew)  =  k
        1270           0 :                      kbase(i,ncvnew) =  k+1
        1271           0 :                      belongcv(i,k)   = .true.
        1272           0 :                      belongcv(i,k+1) = .true.
        1273             :
        1274             :                      ! Calculate internal energy of SRCL. There is no internal energy if
        1275             :                      ! SRCL is elevated from the surface. Also, we simply assume neutral
        1276             :                      ! stability function. Note that this assumption of neutral stability
        1277             :                      ! does not influence numerical calculation- stability functions here
        1278             :                      ! are just for diagnostic output. In general SRCLs other than a SRCL
        1279             :                      ! based at surface with bflxs <= 0, there is no other way but to use
        1280             :                      ! neutral stability function.  However, in case of SRCL based at the
        1281             :                      ! surface,  we can explicitly calculate non-zero stability functions
        1282             :                      ! in a consistent way.   Even though stability functions of SRCL are
        1283             :                      ! just diagnostic outputs not influencing numerical calculations, it
        1284             :                      ! would be informative to write out correct reasonable values rather
        1285             :                      ! than simply assuming neutral stability. I am doing this right now.
        1286             :                      ! Similar calculations were done for the SBCL and when surface inter
        1287             :                      ! facial layer was merged by overlying CL in 'ziscol'.
        1288             :
        1289           0 :                      if( k .lt. pver ) then
        1290             :
        1291           0 :                          wbrk(i,ncvnew) = 0._r8
        1292           0 :                          ebrk(i,ncvnew) = 0._r8
        1293           0 :                          lbrk(i,ncvnew) = 0._r8
        1294           0 :                          ghcl(i,ncvnew) = 0._r8
        1295           0 :                          shcl(i,ncvnew) = 0._r8
        1296           0 :                          smcl(i,ncvnew) = 0._r8
        1297           0 :                          ricl(i,ncvnew) = 0._r8
        1298             :
        1299             :                      else ! Surface-based fog
        1300             :
        1301           0 :                          if( bflxs(i) .gt. 0._r8 ) then    ! Incorporate surface TKE into CL interior energy
        1302             :                                                            ! It is likely that this case cannot exist  since
        1303             :                                                            ! if surface buoyancy flux is positive,  it would
        1304             :                                                            ! have been identified as SBCL in 'zisocl' ahead.
        1305           0 :                              ebrk(i,ncvnew) = tkes(i)
        1306           0 :                              lbrk(i,ncvnew) = z(i,pver)
        1307           0 :                              wbrk(i,ncvnew) = tkes(i) / b1
        1308             :
        1309           0 :                              write(temp_string,*) 'Major mistake in SRCL: bflxs > 0 for surface-based SRCL'
        1310           0 :                              warnstring = trim(warnstring)//temp_string
        1311           0 :                              write(temp_string,*) 'bflxs = ', bflxs(i), &
        1312           0 :                                   'ncvfin_o = ', ncvfin_o(i), &
        1313           0 :                                   'ncvfin_mg = ', ncvfin_mg(i)
        1314           0 :                              warnstring = trim(warnstring)//temp_string
        1315           0 :                              do ks = 1, ncvmax
        1316           0 :                                 write(temp_string,*) 'ncv =', ks, ' ', kbase_o(i,ks), &
        1317           0 :                                      ktop_o(i,ks), kbase_mg(i,ks), ktop_mg(i,ks)
        1318           0 :                                 warnstring = trim(warnstring)//temp_string
        1319             :                              end do
        1320           0 :                              errstring = 'CALEDDY: Major mistake in SRCL: bflxs > 0 for surface-based SRCL'
        1321           0 :                              return
        1322             :                          else                              ! Don't incorporate surface interfacial TKE into CL interior energy
        1323             :
        1324           0 :                              ebrk(i,ncvnew) = 0._r8
        1325           0 :                              lbrk(i,ncvnew) = 0._r8
        1326           0 :                              wbrk(i,ncvnew) = 0._r8
        1327             :
        1328             :                          endif
        1329             :
        1330             :                          ! Calculate stability functions (ghcl, shcl, smcl, ricl) explicitly
        1331             :                          ! using an reverse procedure starting from tkes(i). Note that it is
        1332             :                          ! possible to calculate stability functions even when bflxs < 0.
        1333             :                          ! Previous code just assumed neutral stability functions. Note that
        1334             :                          ! since alph5 = 0.7 > 0, alph3 = -35 < 0, the denominator of gh  is
        1335             :                          ! always positive if bflxs > 0. However, if bflxs < 0,  denominator
        1336             :                          ! can be zero. For this case, we provide a possible maximum negative
        1337             :                          ! value (the most stable state) to gh. Note also tkes(i) is always a
        1338             :                          ! positive value by a limiter. Also, sprod(i,pver+1) > 0 by limiter.
        1339             :
        1340           0 :                          gg = 0.5_r8 * vk * z(i,pver) * bprod(i,pver+1) / ( tkes(i)**(3._r8/2._r8) )
        1341           0 :                          if( abs(alph5-gg*alph3) .le. 1.e-7_r8 ) then
        1342             :                            ! gh = -0.28_r8
        1343             :                            ! gh = -3.5334_r8
        1344             :                              gh = ghmin
        1345             :                          else
        1346           0 :                              gh = gg / ( alph5 - gg * alph3 )
        1347             :                          end if
        1348             :                        ! gh = min(max(gh,-0.28_r8),0.0233_r8)
        1349             :                        ! gh = min(max(gh,-3.5334_r8),0.0233_r8)
        1350           0 :                          gh = min(max(gh,ghmin),0.0233_r8)
        1351           0 :                          ghcl(i,ncvnew) =  gh
        1352           0 :                          shcl(i,ncvnew) =  max(0._r8,alph5/(1._r8+alph3*gh))
        1353           0 :                          smcl(i,ncvnew) =  max(0._r8,(alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
        1354           0 :                          ricl(i,ncvnew) = -(smcl(i,ncvnew)/shcl(i,ncvnew))*(bprod(i,pver+1)/sprod(i,pver+1))
        1355             :
        1356             :                        ! 'ncvsurf' is CL regime index based at the surface. If there is no
        1357             :                        ! such regime, then 'ncvsurf = 0'.
        1358             :
        1359           0 :                          ncvsurf = ncvnew
        1360             :
        1361             :                       end if
        1362             :
        1363             :                   end if
        1364             :
        1365             :               end if
        1366             :
        1367             :           end if
        1368             :
        1369           0 :    220 continue
        1370             :
        1371             :        end do ! End of 'k' loop where 'k' is a grid layer index running from 'pver' to 2
        1372             :
        1373             :    222 continue
        1374             :
        1375             :        ! -------------------------------------------------------------------------- !
        1376             :        ! Up to this point, we identified all kinds of CL regimes :                  !
        1377             :        !   1. A SBCL. By construction, 'bflxs > 0' for SBCL.                        !
        1378             :        !   2. Surface-based CL with multiple layers and 'bflxs =< 0'                !
        1379             :        !   3. Surface-based CL with multiple layers and 'bflxs > 0'                 !
        1380             :        !   4. Regular elevated CL with two entraining interfaces                    !
        1381             :        !   5. SRCLs. If SRCL is based at surface, it will be bflxs < 0.             !
        1382             :        ! '1-4' were identified from 'zisocl' while '5' were identified separately   !
        1383             :        ! after performing 'zisocl'. CL regime index of '1-4' increases with height  !
        1384             :        ! ( e.g., CL = 1 is the CL regime nearest to the surface ) while CL regime   !
        1385             :        ! index of SRCL is simply appended after the final index of CL regimes from  !
        1386             :        ! 'zisocl'. However, CL regime indices of SRCLs itself increases with height !
        1387             :        ! when there are multiple SRCLs, similar to the regular CLs from 'zisocl'.   !
        1388             :        ! -------------------------------------------------------------------------- !
        1389             :
        1390             :        ! Diagnostic output of final CL regimes indices
        1391             :
        1392           0 :        do k = 1, ncvmax
        1393           0 :           kbase_f(i,k) = real(kbase(i,k))
        1394           0 :           ktop_f(i,k)  = real(ktop(i,k))
        1395           0 :           ncvfin_f(i)  = real(ncvfin(i))
        1396             :        end do
        1397             :
        1398             :        ! --------------------------------------------------------------------- !
        1399             :        ! Compute radf for each CL in column by calling subroutine compute_radf !
        1400             :        ! --------------------------------------------------------------------- !
        1401             :        call compute_radf( choice_radf, i, pcols, pver, ncvmax,  ncvfin, ktop, qmin,         &
        1402           0 :                           ql, pi, qrlw, g, cldeff, zi, chs, lwp_CL(i,:), opt_depth_CL(i,:), &
        1403           0 :                           radinvfrac_CL(i,:), radf_CL(i,:) )
        1404             :
        1405             :        ! ---------------------------------------- !
        1406             :        ! Perform do loop for individual CL regime !
        1407             :        ! ---------------------------------------- ! -------------------------------- !
        1408             :        ! For individual CLs, compute                                                 !
        1409             :        !   1. Entrainment rates at the CL top and (if any) base interfaces using     !
        1410             :        !      appropriate entrainment closure (current code use 'wstar' closure).    !
        1411             :        !   2. Net CL mean (i.e., including entrainment contribution) TKE (ebrk)      !
        1412             :        !      and normalized TKE (wbrk).                                             !
        1413             :        !   3. TKE (tke) and normalized TKE (wcap) profiles at all CL interfaces.     !
        1414             :        !   4. ( kvm, kvh ) profiles at all CL interfaces.                            !
        1415             :        !   5. ( bprod, sprod ) profiles at all CL interfaces.                        !
        1416             :        ! Also calculate                                                              !
        1417             :        !   1. PBL height as the top external interface of surface-based CL, if any.  !
        1418             :        !   2. Characteristic excesses of convective 'updraft velocity (wpert)',      !
        1419             :        !      'temperature (tpert)', and 'moisture (qpert)' in the surface-based CL, !
        1420             :        !      if any, for use in the separate convection scheme.                     !
        1421             :        ! If there is no surface-based CL, 'PBL height' and 'convective excesses' are !
        1422             :        ! calculated later from surface-based STL (Stable Turbulent Layer) properties.!
        1423             :        ! --------------------------------------------------------------------------- !
        1424             :
        1425           0 :        ktblw = 0
        1426           0 :        do ncv = 1, ncvfin(i)
        1427             :
        1428           0 :           kt = ktop(i,ncv)
        1429           0 :           kb = kbase(i,ncv)
        1430             :
        1431           0 :           lwp        = lwp_CL(i,ncv)
        1432           0 :           opt_depth  = opt_depth_CL(i,ncv)
        1433           0 :           radinvfrac = radinvfrac_CL(i,ncv)
        1434           0 :           radf       = radf_CL(i, ncv)
        1435             :
        1436             :           ! Check whether surface interface is energetically interior or not.
        1437           0 :           if( kb .eq. (pver+1) .and. bflxs(i) .le. 0._r8 ) then
        1438           0 :               lbulk = zi(i,kt) - z(i,pver)
        1439             :           else
        1440           0 :               lbulk = zi(i,kt) - zi(i,kb)
        1441             :           end if
        1442           0 :           lbulk = min( lbulk, lbulk_max )
        1443             :
        1444             :           ! Calculate 'turbulent length scale (leng)' and 'normalized TKE (wcap)'
        1445             :           ! at all CL interfaces except the surface.  Note that below 'wcap' at
        1446             :           ! external interfaces are not correct. However, it does not influence
        1447             :           ! numerical calculation and correct normalized TKE at the entraining
        1448             :           ! interfaces will be re-calculated at the end of this 'do ncv' loop.
        1449             :
        1450           0 :           do k = min(kb,pver), kt, -1
        1451             :              if( choice_tunl .eq. 'rampcl' ) then
        1452             :                ! In order to treat the case of 'ricl(i,ncv) >> 0' of surface-based SRCL
        1453             :                ! with 'bflxs(i) < 0._r8', I changed ricl(i,ncv) -> min(0._r8,ricl(i,ncv))
        1454             :                ! in the below exponential. This is necessary to prevent the model crash
        1455             :                ! by too large values (e.g., 700) of ricl(i,ncv)
        1456           0 :                  tunlramp = ctunl*tunl*(1._r8-(1._r8-1._r8/ctunl)*exp(min(0._r8,ricl(i,ncv))))
        1457           0 :                  tunlramp = min(max(tunlramp,tunl),ctunl*tunl)
        1458             :              elseif( choice_tunl .eq. 'rampsl' ) then
        1459             :                  tunlramp = ctunl*tunl
        1460             :                ! tunlramp = 0.765_r8
        1461             :              else
        1462             :                  tunlramp = tunl
        1463             :              endif
        1464             :              if( choice_leng .eq. 'origin' ) then
        1465           0 :                  leng(i,k) = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        1466             :                ! leng(i,k) = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
        1467             :              else
        1468             :                  leng(i,k) = min( vk*zi(i,k), tunlramp*lbulk )
        1469             :              endif
        1470           0 :              leng(i,k) = min(leng_max(k), leng(i,k))
        1471           0 :              wcap(i,k) = (leng(i,k)**2) * (-shcl(i,ncv)*n2(i,k)+smcl(i,ncv)*s2(i,k))
        1472             :           end do ! k
        1473             :
        1474             :           ! Calculate basic cross-interface variables ( jump condition ) across the
        1475             :           ! base external interface of CL.
        1476             :
        1477           0 :           if( kb .lt. pver+1 ) then
        1478             :
        1479           0 :               jbzm = z(i,kb-1) - z(i,kb)                                      ! Interfacial layer thickness [m]
        1480           0 :               jbsl = sl(i,kb-1) - sl(i,kb)                                    ! Interfacial jump of 'sl' [J/kg]
        1481           0 :               jbqt = qt(i,kb-1) - qt(i,kb)                                    ! Interfacial jump of 'qt' [kg/kg]
        1482           0 :               jbbu = n2(i,kb) * jbzm                                          ! Interfacial buoyancy jump [m/s2]
        1483             :                                                                               ! considering saturation ( > 0 )
        1484           0 :               jbbu = max(jbbu,jbumin)                                         ! Set minimum buoyancy jump, jbumin = 1.e-3
        1485           0 :               jbu  = u(i,kb-1) - u(i,kb)                                      ! Interfacial jump of 'u' [m/s]
        1486           0 :               jbv  = v(i,kb-1) - v(i,kb)                                      ! Interfacial jump of 'v' [m/s]
        1487           0 :               ch   = (1._r8 -sflh(i,kb-1))*chu(i,kb) + sflh(i,kb-1)*chs(i,kb) ! Buoyancy coefficient just above the base interface
        1488           0 :               cm   = (1._r8 -sflh(i,kb-1))*cmu(i,kb) + sflh(i,kb-1)*cms(i,kb) ! Buoyancy coefficient just above the base interface
        1489           0 :               n2hb = (ch*jbsl + cm*jbqt)/jbzm                                 ! Buoyancy frequency [s-2]
        1490             :                                                                               ! just above the base interface
        1491           0 :               vyb  = n2hb*jbzm/jbbu                                           ! Ratio of 'n2hb/n2' at 'kb' interface
        1492           0 :               vub  = min(1._r8,(jbu**2+jbv**2)/(jbbu*jbzm) )                  ! Ratio of 's2/n2 = 1/Ri' at 'kb' interface
        1493             :
        1494             :           else
        1495             :
        1496             :             ! Below setting is necessary for consistent treatment when 'kb' is at the surface.
        1497             :               jbbu = 0._r8
        1498             :               n2hb = 0._r8
        1499             :               vyb  = 0._r8
        1500             :               vub  = 0._r8
        1501             :               web  = 0._r8
        1502             :
        1503             :           end if
        1504             :
        1505             :           ! Calculate basic cross-interface variables ( jump condition ) across the
        1506             :           ! top external interface of CL. The meanings of variables are similar to
        1507             :           ! the ones at the base interface.
        1508             :
        1509           0 :           jtzm = z(i,kt-1) - z(i,kt)
        1510           0 :           jtsl = sl(i,kt-1) - sl(i,kt)
        1511           0 :           jtqt = qt(i,kt-1) - qt(i,kt)
        1512           0 :           jtbu = n2(i,kt)*jtzm                                      ! Note : 'jtbu' is guaranteed positive by definition of CL top.
        1513           0 :           jtbu = max(jtbu,jbumin)                                   ! But threshold it anyway to be sure.
        1514           0 :           jtu  = u(i,kt-1) - u(i,kt)
        1515           0 :           jtv  = v(i,kt-1) - v(i,kt)
        1516           0 :           ch   = (1._r8 -sfuh(i,kt))*chu(i,kt) + sfuh(i,kt)*chs(i,kt)
        1517           0 :           cm   = (1._r8 -sfuh(i,kt))*cmu(i,kt) + sfuh(i,kt)*cms(i,kt)
        1518           0 :           n2ht = (ch*jtsl + cm*jtqt)/jtzm
        1519           0 :           vyt  = n2ht*jtzm/jtbu
        1520           0 :           vut  = min(1._r8,(jtu**2+jtv**2)/(jtbu*jtzm))
        1521             :
        1522             :           ! Evaporative enhancement factor of entrainment rate at the CL top interface, evhc.
        1523             :           ! We take the full inversion strength to be 'jt2slv = slv(i,kt-2)-slv(i,kt)'
        1524             :           ! where 'kt-1' is in the ambiguous layer. However, for a cloud-topped CL overlain
        1525             :           ! by another CL, it is possible that 'slv(i,kt-2) < slv(i,kt)'. To avoid negative
        1526             :           ! or excessive evhc, we lower-bound jt2slv and upper-bound evhc.  Note 'jtslv' is
        1527             :           ! used only for calculating 'evhc' : when calculating entrainment rate,   we will
        1528             :           ! use normal interfacial buoyancy jump across CL top interface.
        1529             :
        1530           0 :           evhc   = 1._r8
        1531           0 :           jt2slv = 0._r8
        1532             :
        1533             :         ! Modification : I should check whether below 'jbumin' produces reasonable limiting value.
        1534             :         !                In addition, our current formulation does not consider ice contribution.
        1535             :
        1536             :           if( choice_evhc .eq. 'orig' ) then
        1537             :
        1538             :               if( ql(i,kt) .gt. qmin .and. ql(i,kt-1) .lt. qmin ) then
        1539             :                   jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
        1540             :                   jt2slv = max( jt2slv, jbumin*slv(i,kt-1)/g )
        1541             :                   evhc   = 1._r8 + a2l * a3l * latvap * ql(i,kt) / jt2slv
        1542             :                   evhc   = min( evhc, evhcmax )
        1543             :               end if
        1544             :
        1545             :           elseif( choice_evhc .eq. 'ramp' ) then
        1546             :
        1547             :               jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
        1548             :               jt2slv = max( jt2slv, jbumin*slv(i,kt-1)/g )
        1549             :               evhc   = 1._r8 + max(cldeff(i,kt)-cldeff(i,kt-1),0._r8) * a2l * a3l * latvap * ql(i,kt) / jt2slv
        1550             :               evhc   = min( evhc, evhcmax )
        1551             :
        1552             :           elseif( choice_evhc .eq. 'maxi' ) then
        1553             :
        1554           0 :               qleff  = max( ql(i,kt-1), ql(i,kt) )
        1555           0 :               jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
        1556           0 :               jt2slv = max( jt2slv, jbumin*slv(i,kt-1)/g )
        1557           0 :               evhc   = 1._r8 + a2l * a3l * latvap * qleff / jt2slv
        1558           0 :               evhc   = min( evhc, evhcmax )
        1559             :
        1560             :           endif
        1561             :
        1562             :           ! ------------------------------------------------------------------- !
        1563             :           ! Calculate 'wstar3' by summing buoyancy productions within CL from   !
        1564             :           !   1. Interior buoyancy production ( bprod: fcn of TKE )             !
        1565             :           !   2. Cloud-top radiative cooling                                    !
        1566             :           !   3. Surface buoyancy flux contribution only when bflxs > 0.        !
        1567             :           !      Note that master length scale, lbulk, has already been         !
        1568             :           !      corrctly defined at the first part of this 'do ncv' loop       !
        1569             :           !      considering the sign of bflxs.                                 !
        1570             :           ! This 'wstar3' is used for calculation of entrainment rate.          !
        1571             :           ! Note that this 'wstar3' formula does not include shear production   !
        1572             :           ! and the effect of drizzle, which should be included later.          !
        1573             :           ! Q : Strictly speaking, in calculating interior buoyancy production, !
        1574             :           !     the use of 'bprod' is not correct, since 'bprod' is not correct !
        1575             :           !     value but initially guessed value.   More reasonably, we should !
        1576             :           !     use '-leng(i,k)*sqrt(b1*wcap(i,k))*shcl(i,ncv)*n2(i,k)' instead !
        1577             :           !     of 'bprod(i,k)', although this is still an  approximation since !
        1578             :           !     tke(i,k) is not exactly 'b1*wcap(i,k)'  due to a transport term.!
        1579             :           !     However since iterative calculation will be performed after all,!
        1580             :           !     below might also be OK. But I should test this alternative.     !
        1581             :           ! ------------------------------------------------------------------- !
        1582             :
        1583           0 :           dzht   = zi(i,kt)  - z(i,kt)     ! Thickness of CL top half-layer
        1584           0 :           dzhb   = z(i,kb-1) - zi(i,kb)    ! Thickness of CL bot half-layer
        1585           0 :           wstar3 = radf * dzht
        1586           0 :           do k = kt + 1, kb - 1 ! If 'kt = kb - 1', this loop will not be performed.
        1587           0 :                wstar3 =  wstar3 + bprod(i,k) * ( z(i,k-1) - z(i,k) )
        1588             :              ! Below is an alternative which may speed up convergence.
        1589             :              ! However, for interfaces merged into original CL, it can
        1590             :              ! be 'wcap(i,k)<0' since 'n2(i,k)>0'.  Thus, I should use
        1591             :              ! the above original one.
        1592             :              ! wstar3 =  wstar3 - leng(i,k)*sqrt(b1*wcap(i,k))*shcl(i,ncv)*n2(i,k)* &
        1593             :              !                    (z(i,k-1) - z(i,k))
        1594             :           end do
        1595           0 :           if( kb .eq. (pver+1) .and. bflxs(i) .gt. 0._r8 ) then
        1596           0 :              wstar3 = wstar3 + bflxs(i) * dzhb
        1597             :            ! wstar3 = wstar3 + bprod(i,pver+1) * dzhb
        1598             :           end if
        1599           0 :           wstar3 = max( 2.5_r8 * wstar3, 0._r8 )
        1600             :
        1601             :           ! -------------------------------------------------------------- !
        1602             :           ! Below single block is for 'sedimentation-entrainment feedback' !
        1603             :           ! -------------------------------------------------------------- !
        1604             :
        1605             :           if( id_sedfact ) then
        1606             :             ! wsed    = 7.8e5_r8*(ql(i,kt)/ncliq(i,kt))**(2._r8/3._r8)
        1607             :               sedfact = exp(-ased*wsedl(i,kt)/(wstar3**(1._r8/3._r8)+1.e-6_r8))
        1608             :               wsed_CL(i,ncv) = wsedl(i,kt)
        1609             :               if( choice_evhc .eq. 'orig' ) then
        1610             :                   if (ql(i,kt).gt.qmin .and. ql(i,kt-1).lt.qmin) then
        1611             :                       jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
        1612             :                       jt2slv = max(jt2slv, jbumin*slv(i,kt-1)/g)
        1613             :                       evhc = 1._r8+sedfact*a2l*a3l*latvap*ql(i,kt) / jt2slv
        1614             :                       evhc = min(evhc,evhcmax)
        1615             :                   end if
        1616             :               elseif( choice_evhc .eq. 'ramp' ) then
        1617             :                   jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
        1618             :                   jt2slv = max(jt2slv, jbumin*slv(i,kt-1)/g)
        1619             :                   evhc = 1._r8+max(cldeff(i,kt)-cldeff(i,kt-1),0._r8)*sedfact*a2l*a3l*latvap*ql(i,kt) / jt2slv
        1620             :                   evhc = min(evhc,evhcmax)
        1621             :               elseif( choice_evhc .eq. 'maxi' ) then
        1622             :                   qleff  = max(ql(i,kt-1),ql(i,kt))
        1623             :                   jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
        1624             :                   jt2slv = max(jt2slv, jbumin*slv(i,kt-1)/g)
        1625             :                   evhc = 1._r8+sedfact*a2l*a3l*latvap*qleff / jt2slv
        1626             :                   evhc = min(evhc,evhcmax)
        1627             :               endif
        1628             :           endif
        1629             :
        1630             :           ! -------------------------------------------------------------------------- !
        1631             :           ! Now diagnose CL top and bottom entrainment rates (and the contribution of  !
        1632             :           ! top/bottom entrainments to wstar3) using entrainment closures of the form  !
        1633             :           !                                                                            !
        1634             :           !                   wet = cet*wstar3, web = ceb*wstar3                       !
        1635             :           !                                                                            !
        1636             :           ! where cet and ceb depend on the entrainment interface jumps, ql, etc.      !
        1637             :           ! No entrainment is diagnosed unless the wstar3 > 0. Note '1/wstar3fact' is  !
        1638             :           ! a factor indicating the enhancement of wstar3 due to entrainment process.  !
        1639             :           ! Q : Below setting of 'wstar3fact = max(..,0.5)'might prevent the possible  !
        1640             :           !     case when buoyancy consumption by entrainment is  stronger than cloud  !
        1641             :           !     top radiative cooling production. Is that OK ? No.  According to bulk  !
        1642             :           !     modeling study, entrainment buoyancy consumption was always a certain  !
        1643             :           !     fraction of other net productions, rather than a separate sum.  Thus,  !
        1644             :           !     below max limit of wstar3fact is correct.   'wstar3fact = max(.,0.5)'  !
        1645             :           !     prevents unreasonable enhancement of CL entrainment rate by cloud-top  !
        1646             :           !     entrainment instability, CTEI.                                         !
        1647             :           ! Q : Use of the same dry entrainment coefficient, 'a1i' both at the CL  top !
        1648             :           !     and base interfaces may result in too small 'wstar3' and 'ebrk' below, !
        1649             :           !     as was seen in my generalized bulk modeling study. This should be re-  !
        1650             :           !     considered later                                                       !
        1651             :           ! -------------------------------------------------------------------------- !
        1652             :
        1653           0 :           if( wstar3 .gt. 0._r8 ) then
        1654           0 :               cet = a1i * evhc / ( jtbu * lbulk )
        1655           0 :               if( kb .eq. pver + 1 ) then
        1656           0 :                   wstar3fact = max( 1._r8 + 2.5_r8 * cet * n2ht * jtzm * dzht, wstar3factcrit )
        1657             :               else
        1658           0 :                   ceb = a1i / ( jbbu * lbulk )
        1659             :                   wstar3fact = max( 1._r8 + 2.5_r8 * cet * n2ht * jtzm * dzht &
        1660           0 :                                           + 2.5_r8 * ceb * n2hb * jbzm * dzhb, wstar3factcrit )
        1661             :               end if
        1662           0 :               wstar3 = wstar3 / wstar3fact
        1663             :           else ! wstar3 == 0
        1664             :               wstar3fact = 0._r8 ! This is just for dianostic output
        1665             :               cet        = 0._r8
        1666             :               ceb        = 0._r8
        1667             :           end if
        1668             :
        1669             :           ! ---------------------------------------------------------------------------- !
        1670             :           ! Calculate net CL mean TKE including entrainment contribution by solving a    !
        1671             :           ! canonical cubic equation. The solution of cubic equ. is 'rootp**2 = ebrk'    !
        1672             :           ! where 'ebrk' originally (before solving cubic eq.) was interior CL mean TKE, !
        1673             :           ! but after solving cubic equation,  it is replaced by net CL mean TKE in the  !
        1674             :           ! same variable 'ebrk'.                                                        !
        1675             :           ! ---------------------------------------------------------------------------- !
        1676             :           ! Solve cubic equation (canonical form for analytic solution)                  !
        1677             :           !   r^3 - 3*trmp*r - 2*trmq = 0,   r = sqrt<e>                                 !
        1678             :           ! to estimate <e> for CL, derived from layer-mean TKE balance:                 !
        1679             :           !                                                                              !
        1680             :           !   <e>^(3/2)/(b_1*<l>) \approx <B + S>   (*)                                  !
        1681             :           !   <B+S> = (<B+S>_int * l_int + <B+S>_et * dzt + <B+S>_eb * dzb)/lbulk        !
        1682             :           !   <B+S>_int = <e>^(1/2)/(b_1*<l>)*<e>_int                                    !
        1683             :           !   <B+S>_et  = (-vyt+vut)*wet*jtbu + radf                                     !
        1684             :           !   <B+S>_eb  = (-vyb+vub)*web*jbbu                                            !
        1685             :           !                                                                              !
        1686             :           ! where:                                                                       !
        1687             :           !   <> denotes a vertical avg (over the whole CL unless indicated)             !
        1688             :           !   l_int (called lbrk below) is aggregate thickness of interior CL layers     !
        1689             :           !   dzt = zi(i,kt)-z(i,kt)   is thickness of top entrainment layer             !
        1690             :           !   dzb = z(i,kb-1)-zi(i,kb) is thickness of bot entrainment layer             !
        1691             :           !   <e>_int (called ebrk below) is the CL-mean TKE if only interior            !
        1692             :           !                               interfaces contributed.                        !
        1693             :           !   wet, web                  are top. bottom entrainment rates                !
        1694             :           !                                                                              !
        1695             :           ! For a single-level radiatively-driven convective layer, there are no         !
        1696             :           ! interior interfaces so 'ebrk' = 'lbrk' = 0. If the CL goes to the            !
        1697             :           ! surface, 'vyb' and 'vub' are set to zero before and 'ebrk' and 'lbrk'        !
        1698             :           ! have already incorporated the surface interfacial layer contribution,        !
        1699             :           ! so the same formulas still apply.                                            !
        1700             :           !                                                                              !
        1701             :           ! In the original formulation based on TKE,                                    !
        1702             :           !    wet*jtbu = a1l*evhc*<e>^3/2/leng(i,kt)                                    !
        1703             :           !    web*jbbu = a1l*<e>^3/2/leng(i,kt)                                         !
        1704             :           !                                                                              !
        1705             :           ! In the wstar formulation                                                     !
        1706             :           !    wet*jtbu = a1i*evhc*wstar3/lbulk                                          !
        1707             :           !    web*jbbu = a1i*wstar3/lbulk,                                              !
        1708             :           ! ---------------------------------------------------------------------------- !
        1709             :
        1710           0 :           fact = ( evhc * ( -vyt + vut ) * dzht + ( -vyb + vub ) * dzhb * leng(i,kb) / leng(i,kt) ) / lbulk
        1711             :
        1712           0 :           if( wstarent ) then
        1713             :
        1714             :               ! (Option 1) 'wstar' entrainment formulation
        1715             :               ! Here trmq can have either sign, and will usually be nonzero even for non-
        1716             :               ! cloud topped CLs.  If trmq > 0, there will be two positive roots r; we take
        1717             :               ! the larger one. Why ? If necessary, we limit entrainment and wstar to prevent
        1718             :               ! a solution with r < ccrit*wstar ( Why ? ) where we take ccrit = 0.5.
        1719             :
        1720           0 :               trma = 1._r8
        1721           0 :               trmp = ebrk(i,ncv) * ( lbrk(i,ncv) / lbulk ) / 3._r8 + ntzero
        1722           0 :               trmq = 0.5_r8 * b1 * ( leng(i,kt)  / lbulk ) * ( radf * dzht + a1i * fact * wstar3 )
        1723             :
        1724             :               ! Check if there is an acceptable root with r > rcrit = ccrit*wstar.
        1725             :               ! To do this, first find local minimum fmin of the cubic f(r) at sqrt(p),
        1726             :               ! and value fcrit = f(rcrit).
        1727             :
        1728           0 :               rmin  = sqrt(trmp)
        1729           0 :               fmin  = rmin * ( rmin * rmin - 3._r8 * trmp ) - 2._r8 * trmq
        1730           0 :               wstar = wstar3**onet
        1731           0 :               rcrit = ccrit * wstar
        1732           0 :               fcrit = rcrit * ( rcrit * rcrit - 3._r8 * trmp ) - 2._r8 * trmq
        1733             :
        1734             :               ! No acceptable root exists (noroot = .true.) if either:
        1735             :               !    1) rmin < rcrit (in which case cubic is monotone increasing for r > rcrit)
        1736             :               !       and f(rcrit) > 0.
        1737             :               ! or 2) rmin > rcrit (in which case min of f(r) in r > rcrit is at rmin)
        1738             :               !       and f(rmin) > 0.
        1739             :               ! In this case, we reduce entrainment and wstar3 such that r/wstar = ccrit;
        1740             :               ! this changes the coefficients of the cubic.   It might be informative to
        1741             :               ! check when and how many 'noroot' cases occur,  since when 'noroot',   we
        1742             :               ! will impose arbitrary limit on 'wstar3, wet, web, and ebrk' using ccrit.
        1743             :
        1744             :               noroot = ( ( rmin .lt. rcrit ) .and. ( fcrit .gt. 0._r8 ) ) &
        1745           0 :                   .or. ( ( rmin .ge. rcrit ) .and. ( fmin  .gt. 0._r8 ) )
        1746           0 :               if( noroot ) then ! Solve cubic for r
        1747           0 :                   trma = 1._r8 - b1 * ( leng(i,kt) / lbulk ) * a1i * fact / ccrit**3
        1748           0 :                   trma = max( trma, 0.5_r8 )  ! Limit entrainment enhancement of ebrk
        1749           0 :                   trmp = trmp / trma
        1750           0 :                   trmq = 0.5_r8 * b1 * ( leng(i,kt) / lbulk ) * radf * dzht / trma
        1751             :               end if   ! noroot
        1752             :
        1753             :               ! Solve the cubic equation
        1754             :
        1755           0 :               qq = trmq**2 - trmp**3
        1756           0 :               if( qq .ge. 0._r8 ) then
        1757           0 :                   rootp = ( trmq + sqrt(qq) )**(1._r8/3._r8) + ( max( trmq - sqrt(qq), 0._r8 ) )**(1._r8/3._r8)
        1758             :               else
        1759           0 :                   rootp = 2._r8 * sqrt(trmp) * cos( acos( trmq / sqrt(trmp**3) ) / 3._r8 )
        1760             :               end if
        1761             :
        1762             :               ! Adjust 'wstar3' only if there is 'noroot'.
        1763             :               ! And calculate entrainment rates at the top and base interfaces.
        1764             :
        1765           0 :               if( noroot )  wstar3 = ( rootp / ccrit )**3     ! Adjust wstar3
        1766           0 :               wet = cet * wstar3                              ! Find entrainment rates
        1767           0 :               if( kb .lt. pver + 1 ) web = ceb * wstar3       ! When 'kb.eq.pver+1', it was set to web=0.
        1768             :
        1769             :           else !
        1770             :
        1771             :               ! (Option.2) wstarentr = .false. Use original entrainment formulation.
        1772             :               ! trmp > 0 if there are interior interfaces in CL, trmp = 0 otherwise.
        1773             :               ! trmq > 0 if there is cloudtop radiative cooling, trmq = 0 otherwise.
        1774             :
        1775           0 :               trma = 1._r8 - b1 * a1l * fact
        1776           0 :               trma = max( trma, 0.5_r8 )  ! Prevents runaway entrainment instability
        1777           0 :               trmp = ebrk(i,ncv) * ( lbrk(i,ncv) / lbulk ) / ( 3._r8 * trma )
        1778           0 :               trmq = 0.5_r8 * b1 * ( leng(i,kt)  / lbulk ) * radf * dzht / trma
        1779             :
        1780           0 :               qq = trmq**2 - trmp**3
        1781           0 :               if( qq .ge. 0._r8 ) then
        1782           0 :                   rootp = ( trmq + sqrt(qq) )**(1._r8/3._r8) + ( max( trmq - sqrt(qq), 0._r8 ) )**(1._r8/3._r8)
        1783             :               else ! Also part of case 3
        1784           0 :                   rootp = 2._r8 * sqrt(trmp) * cos( acos( trmq / sqrt(trmp**3) ) / 3._r8 )
        1785             :               end if   ! qq
        1786             :
        1787             :              ! Find entrainment rates and limit them by free-entrainment values a1l*sqrt(e)
        1788             :
        1789           0 :               wet = a1l * rootp * min( evhc * rootp**2 / ( leng(i,kt) * jtbu ), 1._r8 )
        1790           0 :               if( kb .lt. pver + 1 ) web = a1l * rootp * min( evhc * rootp**2 / ( leng(i,kb) * jbbu ), 1._r8 )
        1791             :
        1792             :           end if ! wstarentr
        1793             :
        1794             :           ! ---------------------------------------------------- !
        1795             :           ! Finally, get the net CL mean TKE and normalized TKE  !
        1796             :           ! ---------------------------------------------------- !
        1797             :
        1798           0 :           ebrk(i,ncv) = rootp**2
        1799           0 :           ebrk(i,ncv) = min(ebrk(i,ncv),tkemax) ! Limit CL-avg TKE used for entrainment
        1800           0 :           wbrk(i,ncv) = ebrk(i,ncv)/b1
        1801             :
        1802             :           ! The only way ebrk = 0 is for SRCL which are actually radiatively cooled
        1803             :           ! at top interface. In this case, we remove 'convective' label from the
        1804             :           ! interfaces around this layer. This case should now be impossible, so
        1805             :           ! we flag it. Q: I can't understand why this case is impossible now. Maybe,
        1806             :           ! due to various limiting procedures used in solving cubic equation ?
        1807             :           ! In case of SRCL, 'ebrk' should be positive due to cloud top LW radiative
        1808             :           ! cooling contribution, although 'ebrk(internal)' of SRCL before including
        1809             :           ! entrainment contribution (which include LW cooling contribution also) is
        1810             :           ! zero.
        1811             :
        1812           0 :           if( ebrk(i,ncv) .le. 0._r8 ) then
        1813           0 :               write(temp_string,*) 'CALEDDY: Warning, CL with zero TKE, i, kt, kb ', i, kt, kb
        1814           0 :               warnstring = trim(warnstring)//temp_string
        1815           0 :               belongcv(i,kt) = .false.
        1816           0 :               belongcv(i,kb) = .false.
        1817             :           end if
        1818             :
        1819             :           ! ----------------------------------------------------------------------- !
        1820             :           ! Calculate complete TKE profiles at all CL interfaces, capped by tkemax. !
        1821             :           ! We approximate TKE = <e> at entrainment interfaces. However when CL is  !
        1822             :           ! based at surface, correct 'tkes' will be inserted to tke(i,pver+1).     !
        1823             :           ! Note that this approximation at CL external interfaces do not influence !
        1824             :           ! numerical calculation since 'e' at external interfaces are not used  in !
        1825             :           ! actual numerical calculation afterward. In addition in order to extract !
        1826             :           ! correct TKE averaged over the PBL in the cumulus scheme,it is necessary !
        1827             :           ! to set e = <e> at the top entrainment interface.  Since net CL mean TKE !
        1828             :           ! 'ebrk' obtained by solving cubic equation already includes tkes  ( tkes !
        1829             :           ! is included when bflxs > 0 but not when bflxs <= 0 into internal ebrk ),!
        1830             :           ! 'tkes' should be written to tke(i,pver+1)                               !
        1831             :           ! ----------------------------------------------------------------------- !
        1832             :
        1833             :           ! 1. At internal interfaces
        1834           0 :           do k = kb - 1, kt + 1, -1
        1835           0 :              rcap = ( b1 * ae + wcap(i,k) / wbrk(i,ncv) ) / ( b1 * ae + 1._r8 )
        1836           0 :              rcap = min( max(rcap,rcapmin), rcapmax )
        1837           0 :              tke(i,k) = ebrk(i,ncv) * rcap
        1838           0 :              tke(i,k) = min( tke(i,k), tkemax )
        1839           0 :              kvh(i,k) = leng(i,k) * sqrt(tke(i,k)) * shcl(i,ncv)
        1840           0 :              kvm(i,k) = leng(i,k) * sqrt(tke(i,k)) * smcl(i,ncv)
        1841           0 :              bprod(i,k) = -kvh(i,k) * n2(i,k)
        1842           0 :              sprod(i,k) =  kvm(i,k) * s2(i,k)
        1843           0 :              turbtype(i,k) = 2                     ! CL interior interfaces.
        1844             :           end do
        1845             :
        1846             :           ! 2. At CL top entrainment interface
        1847           0 :           kentr = wet * jtzm
        1848           0 :           kvh(i,kt) = kentr
        1849           0 :           kvm(i,kt) = kentr
        1850           0 :           bprod(i,kt) = -kentr * n2ht + radf       ! I must use 'n2ht' not 'n2'
        1851           0 :           sprod(i,kt) =  kentr * s2(i,kt)
        1852           0 :           turbtype(i,kt) = 4                       ! CL top entrainment interface
        1853           0 :           trmp = -b1 * ae / ( 1._r8 + b1 * ae )
        1854           0 :           trmq = -(bprod(i,kt)+sprod(i,kt))*b1*leng(i,kt)/(1._r8+b1*ae)/(ebrk(i,ncv)**(3._r8/2._r8))
        1855           0 :           rcap = compute_cubic(0._r8,trmp,trmq)**2._r8
        1856           0 :           rcap = min( max(rcap,rcapmin), rcapmax )
        1857           0 :           tke(i,kt)  = ebrk(i,ncv) * rcap
        1858           0 :           tke(i,kt)  = min( tke(i,kt), tkemax )
        1859             :
        1860             :           ! 3. At CL base entrainment interface and double entraining interfaces
        1861             :           ! When current CL base is also the top interface of CL regime below,
        1862             :           ! simply add the two contributions for calculating eddy diffusivity
        1863             :           ! and buoyancy/shear production. Below code correctly works because
        1864             :           ! we (CL regime index) always go from surface upward.
        1865             :
        1866           0 :           if( kb .lt. pver + 1 ) then
        1867             :
        1868           0 :               kentr = web * jbzm
        1869             :
        1870           0 :               if( kb .ne. ktblw ) then
        1871             :
        1872           0 :                   kvh(i,kb) = kentr
        1873           0 :                   kvm(i,kb) = kentr
        1874           0 :                   bprod(i,kb) = -kvh(i,kb)*n2hb     ! I must use 'n2hb' not 'n2'
        1875           0 :                   sprod(i,kb) =  kvm(i,kb)*s2(i,kb)
        1876           0 :                   turbtype(i,kb) = 3                ! CL base entrainment interface
        1877           0 :                   trmp = -b1*ae/(1._r8+b1*ae)
        1878           0 :                   trmq = -(bprod(i,kb)+sprod(i,kb))*b1*leng(i,kb)/(1._r8+b1*ae)/(ebrk(i,ncv)**(3._r8/2._r8))
        1879           0 :                   rcap = compute_cubic(0._r8,trmp,trmq)**2._r8
        1880           0 :                   rcap = min( max(rcap,rcapmin), rcapmax )
        1881           0 :                   tke(i,kb)  = ebrk(i,ncv) * rcap
        1882           0 :                   tke(i,kb)  = min( tke(i,kb),tkemax )
        1883             :
        1884             :               else
        1885             :
        1886           0 :                   kvh(i,kb) = kvh(i,kb) + kentr
        1887           0 :                   kvm(i,kb) = kvm(i,kb) + kentr
        1888             :                 ! dzhb5 : Half thickness of the lowest  layer of  current CL regime
        1889             :                 ! dzht5 : Half thickness of the highest layer of adjacent CL regime just below current CL.
        1890           0 :                   dzhb5 = z(i,kb-1) - zi(i,kb)
        1891           0 :                   dzht5 = zi(i,kb) - z(i,kb)
        1892           0 :                   bprod(i,kb) = ( dzht5*bprod(i,kb) - dzhb5*kentr*n2hb )     / ( dzhb5 + dzht5 )
        1893           0 :                   sprod(i,kb) = ( dzht5*sprod(i,kb) + dzhb5*kentr*s2(i,kb) ) / ( dzhb5 + dzht5 )
        1894           0 :                   trmp = -b1*ae/(1._r8+b1*ae)
        1895           0 :                   trmq = -kentr*(s2(i,kb)-n2hb)*b1*leng(i,kb)/(1._r8+b1*ae)/(ebrk(i,ncv)**(3._r8/2._r8))
        1896           0 :                   rcap = compute_cubic(0._r8,trmp,trmq)**2._r8
        1897           0 :                   rcap = min( max(rcap,rcapmin), rcapmax )
        1898           0 :                   tke_imsi = ebrk(i,ncv) * rcap
        1899           0 :                   tke_imsi = min( tke_imsi, tkemax )
        1900           0 :                   tke(i,kb)  = ( dzht5*tke(i,kb) + dzhb5*tke_imsi ) / ( dzhb5 + dzht5 )
        1901           0 :                   tke(i,kb)  = min(tke(i,kb),tkemax)
        1902           0 :                   turbtype(i,kb) = 5                ! CL double entraining interface
        1903             :
        1904             :               end if
        1905             :
        1906             :            else
        1907             :
        1908             :              ! If CL base interface is surface, compute similarly using wcap(i,kb)=tkes/b1
        1909             :              ! Even when bflx < 0, use the same formula in order to impose consistency of
        1910             :              ! tke(i,kb) at bflx = 0._r8
        1911             :
        1912           0 :              rcap = (b1*ae + wcap(i,kb)/wbrk(i,ncv))/(b1*ae + 1._r8)
        1913           0 :              rcap = min( max(rcap,rcapmin), rcapmax )
        1914           0 :              tke(i,kb) = ebrk(i,ncv) * rcap
        1915           0 :              tke(i,kb) = min( tke(i,kb),tkemax )
        1916             :
        1917             :           end if
        1918             :
        1919             :           ! Calculate wcap at all interfaces of CL. Put a  minimum threshold on TKE
        1920             :           ! to prevent possible division by zero.  'wcap' at CL internal interfaces
        1921             :           ! are already calculated in the first part of 'do ncv' loop correctly.
        1922             :           ! When 'kb.eq.pver+1', below formula produces the identical result to the
        1923             :           ! 'tkes(i)/b1' if leng(i,kb) is set to vk*z(i,pver). Note  wcap(i,pver+1)
        1924             :           ! is already defined as 'tkes(i)/b1' at the first part of caleddy.
        1925             :
        1926           0 :           wcap(i,kt) = (bprod(i,kt)+sprod(i,kt))*leng(i,kt)/sqrt(max(tke(i,kt),1.e-6_r8))
        1927           0 :           if( kb .lt. pver + 1 ) then
        1928           0 :               wcap(i,kb) = (bprod(i,kb)+sprod(i,kb))*leng(i,kb)/sqrt(max(tke(i,kb),1.e-6_r8))
        1929             :           end if
        1930             :
        1931             :           ! Save the index of upper external interface of current CL-regime in order to
        1932             :           ! handle the case when this interface is also the lower external interface of
        1933             :           ! CL-regime located just above.
        1934             :
        1935           0 :           ktblw = kt
        1936             :
        1937             :           ! Diagnostic Output
        1938             :
        1939           0 :           wet_CL(i,ncv)        = wet
        1940           0 :           web_CL(i,ncv)        = web
        1941           0 :           jtbu_CL(i,ncv)       = jtbu
        1942           0 :           jbbu_CL(i,ncv)       = jbbu
        1943           0 :           evhc_CL(i,ncv)       = evhc
        1944           0 :           jt2slv_CL(i,ncv)     = jt2slv
        1945           0 :           n2ht_CL(i,ncv)       = n2ht
        1946           0 :           n2hb_CL(i,ncv)       = n2hb
        1947           0 :           wstar_CL(i,ncv)      = wstar
        1948           0 :           wstar3fact_CL(i,ncv) = wstar3fact
        1949             :
        1950             :        end do        ! ncv
        1951             :
        1952             :        ! Calculate PBL height and characteristic cumulus excess for use in the
        1953             :        ! cumulus convection shceme. Also define turbulence type at the surface
        1954             :        ! when the lowest CL is based at the surface. These are just diagnostic
        1955             :        ! outputs, not influencing numerical calculation of current PBL scheme.
        1956             :        ! If the lowest CL is based at the surface, define the PBL depth as the
        1957             :        ! CL top interface. The same rule is applied for all CLs including SRCL.
        1958             :
        1959           0 :        if( ncvsurf .gt. 0 ) then
        1960             :
        1961           0 :            ktopbl(i) = ktop(i,ncvsurf)
        1962           0 :            pblh(i)   = zi(i, ktopbl(i))
        1963           0 :            pblhp(i)  = pi(i, ktopbl(i))
        1964           0 :            wpert(i)  = max(wfac*sqrt(ebrk(i,ncvsurf)),wpertmin)
        1965           0 :            tpert(i)  = max(abs(shflx(i)*rrho(i)/cpair)*tfac/wpert(i),0._r8)
        1966           0 :            qpert(i)  = max(abs(qflx(i)*rrho(i))*tfac/wpert(i),0._r8)
        1967             :
        1968           0 :            if( bflxs(i) .gt. 0._r8 ) then
        1969           0 :                turbtype(i,pver+1) = 2 ! CL interior interface
        1970             :            else
        1971           0 :                turbtype(i,pver+1) = 3 ! CL external base interface
        1972             :            endif
        1973             :
        1974           0 :            ipbl(i)  = 1
        1975           0 :            kpblh(i) = max(ktopbl(i)-1, 1)
        1976           0 :            went(i)  = wet_CL(i,ncvsurf)
        1977             :        end if ! End of the calculationf of te properties of surface-based CL.
        1978             :
        1979             :        ! -------------------------------------------- !
        1980             :        ! Treatment of Stable Turbulent Regime ( STL ) !
        1981             :        ! -------------------------------------------- !
        1982             :
        1983             :        ! Identify top and bottom most (internal) interfaces of STL except surface.
        1984             :        ! Also, calculate 'turbulent length scale (leng)' at each STL interfaces.
        1985             :
        1986           0 :        belongst(i,1) = .false.   ! k = 1 (top interface) is assumed non-turbulent
        1987           0 :        do k = 2, pver            ! k is an interface index
        1988           0 :           belongst(i,k) = ( ri(i,k) .lt. ricrit ) .and. ( .not. belongcv(i,k) )
        1989           0 :           if( belongst(i,k) .and. ( .not. belongst(i,k-1) ) ) then
        1990           0 :               kt = k             ! Top interface index of STL
        1991           0 :           elseif( .not. belongst(i,k) .and. belongst(i,k-1) ) then
        1992           0 :               kb = k - 1         ! Base interface index of STL
        1993           0 :               lbulk = z(i,kt-1) - z(i,kb)
        1994           0 :               lbulk = min( lbulk, lbulk_max )
        1995           0 :               do ks = kt, kb
        1996           0 :                  if( choice_tunl .eq. 'rampcl' ) then
        1997             :                      tunlramp = tunl
        1998             :                  elseif( choice_tunl .eq. 'rampsl' ) then
        1999             :                     tunlramp = max( 1.e-3_r8, ctunl * tunl * exp(-log(ctunl)*ri(i,ks)/ricrit) )
        2000             :                   ! tunlramp = 0.065_r8 + 0.7_r8 * exp(-20._r8*ri(i,ks))
        2001             :                  else
        2002             :                     tunlramp = tunl
        2003             :                  endif
        2004             :                  if( choice_leng .eq. 'origin' ) then
        2005           0 :                      leng(i,ks) = ( (vk*zi(i,ks))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        2006             :                    ! leng(i,ks) = vk*zi(i,ks) / (1._r8+vk*zi(i,ks)/(tunlramp*lbulk))
        2007             :                  else
        2008             :                      leng(i,ks) = min( vk*zi(i,ks), tunlramp*lbulk )
        2009             :                  endif
        2010           0 :                  leng(i,ks) = min(leng_max(ks), leng(i,ks))
        2011             :               end do
        2012             :           end if
        2013             :        end do ! k
        2014             :
        2015             :        ! Now look whether STL extends to ground.  If STL extends to surface,
        2016             :        ! re-define master length scale,'lbulk' including surface interfacial
        2017             :        ! layer thickness, and re-calculate turbulent length scale, 'leng' at
        2018             :        ! all STL interfaces again. Note that surface interface is assumed to
        2019             :        ! always be STL if it is not CL.
        2020             :
        2021           0 :        belongst(i,pver+1) = .not. belongcv(i,pver+1)
        2022             :
        2023           0 :        if( belongst(i,pver+1) ) then     ! kb = pver+1 (surface  STL)
        2024             :
        2025           0 :            turbtype(i,pver+1) = 1        ! Surface is STL interface
        2026             :
        2027           0 :            if( belongst(i,pver) ) then   ! STL includes interior
        2028             :              ! 'kt' already defined above as the top interface of STL
        2029           0 :                lbulk = z(i,kt-1)
        2030             :            else                          ! STL with no interior turbulence
        2031           0 :                kt = pver+1
        2032           0 :                lbulk = z(i,kt-1)
        2033             :            end if
        2034           0 :            lbulk = min( lbulk, lbulk_max )
        2035             :
        2036             :            ! PBL height : Layer mid-point just above the highest STL interface
        2037             :            ! Note in contrast to the surface based CL regime where  PBL height
        2038             :            ! was defined at the top external interface, PBL height of  surface
        2039             :            ! based STL is defined as the layer mid-point.
        2040             :
        2041           0 :            ktopbl(i) = kt - 1
        2042           0 :            pblh(i)   = z(i,ktopbl(i))
        2043           0 :            pblhp(i)  = 0.5_r8 * ( pi(i,ktopbl(i)) + pi(i,ktopbl(i)+1) )
        2044             :
        2045             :            ! Re-calculate turbulent length scale including surface interfacial
        2046             :            ! layer contribution to lbulk.
        2047             :
        2048           0 :            do ks = kt, pver
        2049           0 :               if( choice_tunl .eq. 'rampcl' ) then
        2050             :                   tunlramp = tunl
        2051             :               elseif( choice_tunl .eq. 'rampsl' ) then
        2052             :                   tunlramp = max(1.e-3_r8,ctunl*tunl*exp(-log(ctunl)*ri(i,ks)/ricrit))
        2053             :                 ! tunlramp = 0.065_r8 + 0.7_r8 * exp(-20._r8*ri(i,ks))
        2054             :               else
        2055             :                   tunlramp = tunl
        2056             :               endif
        2057             :               if( choice_leng .eq. 'origin' ) then
        2058           0 :                   leng(i,ks) = ( (vk*zi(i,ks))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        2059             :                 ! leng(i,ks) = vk*zi(i,ks) / (1._r8+vk*zi(i,ks)/(tunlramp*lbulk))
        2060             :               else
        2061             :                   leng(i,ks) = min( vk*zi(i,ks), tunlramp*lbulk )
        2062             :               endif
        2063           0 :              leng(i,ks) = min(leng_max(ks), leng(i,ks))
        2064             :            end do ! ks
        2065             :
        2066             :            ! Characteristic cumulus excess of surface-based STL.
        2067             :            ! We may be able to use ustar for wpert.
        2068             :
        2069           0 :            wpert(i) = 0._r8
        2070           0 :            tpert(i) = max(shflx(i)*rrho(i)/cpair*fak/ustar(i),0._r8) ! CCM stable-layer forms
        2071           0 :            qpert(i) = max(qflx(i)*rrho(i)*fak/ustar(i),0._r8)
        2072             :
        2073           0 :            ipbl(i)  = 0
        2074           0 :            kpblh(i) = ktopbl(i)
        2075             :
        2076             :        end if
        2077             :
        2078             :        ! Calculate stability functions and energetics at the STL interfaces
        2079             :        ! except the surface. Note that tke(i,pver+1) and wcap(i,pver+1) are
        2080             :        ! already calculated in the first part of 'caleddy', kvm(i,pver+1) &
        2081             :        ! kvh(i,pver+1) were already initialized to be zero, bprod(i,pver+1)
        2082             :        ! & sprod(i,pver+1) were direcly calculated from the bflxs and ustar.
        2083             :        ! Note transport term is assumed to be negligible at STL interfaces.
        2084             :
        2085           0 :        do k = 2, pver
        2086             :
        2087           0 :           if( belongst(i,k) ) then
        2088             :
        2089           0 :               turbtype(i,k) = 1    ! STL interfaces
        2090           0 :               trma = alph3*alph4exs*ri(i,k) + 2._r8*b1*(alph2-alph4exs*alph5*ri(i,k))
        2091           0 :               trmb = (alph3+alph4exs)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
        2092           0 :               trmc = ri(i,k)
        2093           0 :               det = max(trmb*trmb-4._r8*trma*trmc,0._r8)
        2094             :               ! Sanity Check
        2095           0 :               if( det .lt. 0._r8 ) then
        2096           0 :                  errstring = 'The det < 0. for the STL in UW eddy_diff'
        2097           0 :                  return
        2098             :               end if
        2099           0 :               gh = (-trmb + sqrt(det))/(2._r8*trma)
        2100             :             ! gh = min(max(gh,-0.28_r8),0.0233_r8)
        2101             :             ! gh = min(max(gh,-3.5334_r8),0.0233_r8)
        2102           0 :               gh = min(max(gh,ghmin),0.0233_r8)
        2103           0 :               sh = max(0._r8,alph5/(1._r8+alph3*gh))
        2104           0 :               sm = max(0._r8,(alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
        2105             :
        2106           0 :               tke(i,k)   = b1*(leng(i,k)**2)*(-sh*n2(i,k)+sm*s2(i,k))
        2107           0 :               tke(i,k)   = min(tke(i,k),tkemax)
        2108           0 :               wcap(i,k)  = tke(i,k)/b1
        2109           0 :               kvh(i,k)   = leng(i,k) * sqrt(tke(i,k)) * sh
        2110           0 :               kvm(i,k)   = leng(i,k) * sqrt(tke(i,k)) * sm
        2111           0 :               bprod(i,k) = -kvh(i,k) * n2(i,k)
        2112           0 :               sprod(i,k) =  kvm(i,k) * s2(i,k)
        2113             :
        2114             :           end if
        2115             :
        2116             :        end do  ! k
        2117             :
        2118             :        ! --------------------------------------------------- !
        2119             :        ! End of treatment of Stable Turbulent Regime ( STL ) !
        2120             :        ! --------------------------------------------------- !
        2121             :
        2122             :        ! --------------------------------------------------------------- !
        2123             :        ! Re-computation of eddy diffusivity at the entrainment interface !
        2124             :        ! assuming that it is purely STL (0<Ri<0.19). Note even Ri>0.19,  !
        2125             :        ! turbulent can exist at the entrainment interface since 'Sh,Sm'  !
        2126             :        ! do not necessarily go to zero even when Ri>0.19. Since Ri can   !
        2127             :        ! be fairly larger than 0.19 at the entrainment interface, I      !
        2128             :        ! should set minimum value of 'tke' to be 0. in order to prevent  !
        2129             :        ! sqrt(tke) from being imaginary.                                 !
        2130             :        ! --------------------------------------------------------------- !
        2131             :
        2132             :        ! goto 888
        2133             :
        2134           0 :          do k = 2, pver
        2135             :
        2136           0 :          if( ( turbtype(i,k) .eq. 3 ) .or. ( turbtype(i,k) .eq. 4 ) .or. &
        2137           0 :              ( turbtype(i,k) .eq. 5 ) ) then
        2138             :
        2139           0 :              trma = alph3*alph4exs*ri(i,k) + 2._r8*b1*(alph2-alph4exs*alph5*ri(i,k))
        2140           0 :              trmb = (alph3+alph4exs)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
        2141           0 :              trmc = ri(i,k)
        2142           0 :              det  = max(trmb*trmb-4._r8*trma*trmc,0._r8)
        2143           0 :              gh   = (-trmb + sqrt(det))/(2._r8*trma)
        2144             :            ! gh   = min(max(gh,-0.28_r8),0.0233_r8)
        2145             :            ! gh   = min(max(gh,-3.5334_r8),0.0233_r8)
        2146           0 :              gh   = min(max(gh,ghmin),0.0233_r8)
        2147           0 :              sh   = max(0._r8,alph5/(1._r8+alph3*gh))
        2148           0 :              sm   = max(0._r8,(alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
        2149             :
        2150           0 :              lbulk = z(i,k-1) - z(i,k)
        2151           0 :              lbulk = min( lbulk, lbulk_max )
        2152             :
        2153           0 :              if( choice_tunl .eq. 'rampcl' ) then
        2154             :                  tunlramp = tunl
        2155             :              elseif( choice_tunl .eq. 'rampsl' ) then
        2156             :                  tunlramp = max(1.e-3_r8,ctunl*tunl*exp(-log(ctunl)*ri(i,k)/ricrit))
        2157             :                ! tunlramp = 0.065_r8 + 0.7_r8*exp(-20._r8*ri(i,k))
        2158             :              else
        2159             :                  tunlramp = tunl
        2160             :              endif
        2161             :              if( choice_leng .eq. 'origin' ) then
        2162           0 :                  leng_imsi = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        2163             :                ! leng_imsi = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
        2164             :              else
        2165             :                  leng_imsi = min( vk*zi(i,k), tunlramp*lbulk )
        2166             :              endif
        2167           0 :              leng_imsi = min(leng_max(k), leng_imsi)
        2168             :
        2169           0 :              tke_imsi = b1*(leng_imsi**2)*(-sh*n2(i,k)+sm*s2(i,k))
        2170           0 :              tke_imsi = min(max(tke_imsi,0._r8),tkemax)
        2171           0 :              kvh_imsi = leng_imsi * sqrt(tke_imsi) * sh
        2172           0 :              kvm_imsi = leng_imsi * sqrt(tke_imsi) * sm
        2173             :
        2174           0 :              if( kvh(i,k) .lt. kvh_imsi ) then
        2175           0 :                  kvh(i,k)   =  kvh_imsi
        2176           0 :                  kvm(i,k)   =  kvm_imsi
        2177           0 :                  leng(i,k)  = leng_imsi
        2178           0 :                  tke(i,k)   =  tke_imsi
        2179           0 :                  wcap(i,k)  =  tke_imsi / b1
        2180           0 :                  bprod(i,k) = -kvh_imsi * n2(i,k)
        2181           0 :                  sprod(i,k) =  kvm_imsi * s2(i,k)
        2182           0 :                  turbtype(i,k) = 1          ! This was added on Dec.10.2009 for use in microphysics.
        2183             :              endif
        2184             :
        2185             :          end if
        2186             :
        2187             :          end do
        2188             :
        2189             :  ! 888   continue
        2190             :
        2191             :        ! ------------------------------------------------------------------ !
        2192             :        ! End of recomputation of eddy diffusivity at entrainment interfaces !
        2193             :        ! ------------------------------------------------------------------ !
        2194             :
        2195             :        ! As an option, we can impose a certain minimum back-ground diffusivity.
        2196             :
        2197             :        ! do k = 1, pver+1
        2198             :        !    kvh(i,k) = max(0.01_r8,kvh(i,k))
        2199             :        !    kvm(i,k) = max(0.01_r8,kvm(i,k))
        2200             :        ! enddo
        2201             :
        2202             :        ! --------------------------------------------------------------------- !
        2203             :        ! Diagnostic Output                                                     !
        2204             :        ! Just for diagnostic purpose, calculate stability functions at  each   !
        2205             :        ! interface including surface. Instead of assuming neutral stability,   !
        2206             :        ! explicitly calculate stability functions using an reverse procedure   !
        2207             :        ! starting from tkes(i) similar to the case of SRCL and SBCL in zisocl. !
        2208             :        ! Note that it is possible to calculate stability functions even when   !
        2209             :        ! bflxs < 0. Note that this inverse method allows us to define Ri even  !
        2210             :        ! at the surface. Note also tkes(i) and sprod(i,pver+1) are always      !
        2211             :        ! positive values by limiters (e.g., ustar_min = 0.01).                 !
        2212             :        ! Dec.12.2006 : Also just for diagnostic output, re-set                 !
        2213             :        ! 'bprod(i,pver+1)= bflxs(i)' here. Note that this setting does not     !
        2214             :        ! influence numerical calculation at all - it is just for diagnostic    !
        2215             :        ! output.                                                               !
        2216             :        ! --------------------------------------------------------------------- !
        2217             :
        2218           0 :        bprod(i,pver+1) = bflxs(i)
        2219             :
        2220           0 :        gg = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
        2221           0 :        if( abs(alph5-gg*alph3) .le. 1.e-7_r8 ) then
        2222             :          ! gh = -0.28_r8
        2223             :            if( bprod(i,pver+1) .gt. 0._r8 ) then
        2224             :                gh = -3.5334_r8
        2225             :            else
        2226             :                gh = ghmin
        2227             :            endif
        2228             :        else
        2229           0 :            gh = gg/(alph5-gg*alph3)
        2230             :        end if
        2231             :
        2232             :      ! gh = min(max(gh,-0.28_r8),0.0233_r8)
        2233           0 :        if( bprod(i,pver+1) .gt. 0._r8 ) then
        2234           0 :            gh = min(max(gh,-3.5334_r8),0.0233_r8)
        2235             :        else
        2236           0 :            gh = min(max(gh,ghmin),0.0233_r8)
        2237             :        endif
        2238             :
        2239           0 :        gh_a(i,pver+1) = gh
        2240           0 :        sh_a(i,pver+1) = max(0._r8,alph5/(1._r8+alph3*gh))
        2241           0 :        if( bprod(i,pver+1) .gt. 0._r8 ) then
        2242           0 :            sm_a(i,pver+1) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh))
        2243             :        else
        2244           0 :            sm_a(i,pver+1) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
        2245             :        endif
        2246           0 :        ri_a(i,pver+1)  = -(sm_a(i,pver+1)/sh_a(i,pver+1))*(bprod(i,pver+1)/sprod(i,pver+1))
        2247             :
        2248           0 :        do k = 1, pver
        2249           0 :           if( ri(i,k) .lt. 0._r8 ) then
        2250           0 :               trma = alph3*alph4*ri(i,k) + 2._r8*b1*(alph2-alph4*alph5*ri(i,k))
        2251           0 :               trmb = (alph3+alph4)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
        2252           0 :               trmc = ri(i,k)
        2253           0 :               det  = max(trmb*trmb-4._r8*trma*trmc,0._r8)
        2254           0 :               gh   = (-trmb + sqrt(det))/(2._r8*trma)
        2255           0 :               gh   = min(max(gh,-3.5334_r8),0.0233_r8)
        2256           0 :               gh_a(i,k) = gh
        2257           0 :               sh_a(i,k) = max(0._r8,alph5/(1._r8+alph3*gh))
        2258           0 :               sm_a(i,k) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh))
        2259           0 :               ri_a(i,k) = ri(i,k)
        2260             :           else
        2261           0 :               if( ri(i,k) .gt. ricrit ) then
        2262           0 :                   gh_a(i,k) = ghmin
        2263           0 :                   sh_a(i,k) = 0._r8
        2264           0 :                   sm_a(i,k) = 0._r8
        2265           0 :                   ri_a(i,k) = ri(i,k)
        2266             :               else
        2267           0 :                   trma = alph3*alph4exs*ri(i,k) + 2._r8*b1*(alph2-alph4exs*alph5*ri(i,k))
        2268           0 :                   trmb = (alph3+alph4exs)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
        2269           0 :                   trmc = ri(i,k)
        2270           0 :                   det  = max(trmb*trmb-4._r8*trma*trmc,0._r8)
        2271           0 :                   gh   = (-trmb + sqrt(det))/(2._r8*trma)
        2272           0 :                   gh   = min(max(gh,ghmin),0.0233_r8)
        2273           0 :                   gh_a(i,k) = gh
        2274           0 :                   sh_a(i,k) = max(0._r8,alph5/(1._r8+alph3*gh))
        2275           0 :                   sm_a(i,k) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
        2276           0 :                   ri_a(i,k) = ri(i,k)
        2277             :               endif
        2278             :           endif
        2279             :
        2280             :        end do
        2281             :
        2282             :     end do   ! End of column index loop, i
        2283             :
        2284             :     return
        2285             :
        2286           0 :     end subroutine caleddy
        2287             :
        2288             :     !============================================================================== !
        2289             :     !                                                                               !
        2290             :     !============================================================================== !
        2291             :
        2292           0 :     subroutine exacol( pcols, pver, ncol, ri, bflxs, minpblh, zi, ktop, kbase, ncvfin )
        2293             :
        2294             :     ! ---------------------------------------------------------------------------- !
        2295             :     ! Object : Find unstable CL regimes and determine the indices                  !
        2296             :     !          kbase, ktop which delimit these unstable layers :                   !
        2297             :     !          ri(kbase) > 0 and ri(ktop) > 0, but ri(k) < 0 for ktop < k < kbase. !
        2298             :     ! Author : Chris  Bretherton 08/2000,                                          !
        2299             :     !          Sungsu Park       08/2006, 11/2008                                  !
        2300             :     !----------------------------------------------------------------------------- !
        2301             :
        2302             :     implicit none
        2303             :
        2304             :     ! --------------- !
        2305             :     ! Input variables !
        2306             :     ! --------------- !
        2307             :
        2308             :     integer,  intent(in) :: pcols                  ! Number of atmospheric columns
        2309             :     integer,  intent(in) :: pver                   ! Number of atmospheric vertical layers
        2310             :     integer,  intent(in) :: ncol                   ! Number of atmospheric columns
        2311             :
        2312             :     real(r8), intent(in) :: ri(pcols,pver)         ! Moist gradient Richardson no.
        2313             :     real(r8), intent(in) :: bflxs(pcols)           ! Buoyancy flux at surface
        2314             :     real(r8), intent(in) :: minpblh(pcols)         ! Minimum PBL height based on surface stress
        2315             :     real(r8), intent(in) :: zi(pcols,pver+1)       ! Interface heights
        2316             :
        2317             :     ! ---------------- !
        2318             :     ! Output variables !
        2319             :     ! ---------------- !
        2320             :
        2321             :     integer, intent(out) :: kbase(pcols,ncvmax)    ! External interface index of CL base
        2322             :     integer, intent(out) :: ktop(pcols,ncvmax)     ! External interface index of CL top
        2323             :     integer, intent(out) :: ncvfin(pcols)          ! Total number of CLs
        2324             :
        2325             :     ! --------------- !
        2326             :     ! Local variables !
        2327             :     ! --------------- !
        2328             :
        2329             :     integer              :: i
        2330             :     integer              :: k
        2331             :     integer              :: ncv
        2332             :     real(r8)             :: rimaxentr
        2333           0 :     real(r8)             :: riex(pver+1)           ! Column Ri profile extended to surface
        2334             :
        2335             :     ! ----------------------- !
        2336             :     ! Main Computation Begins !
        2337             :     ! ----------------------- !
        2338             :
        2339           0 :     do i = 1, ncol
        2340           0 :        ncvfin(i) = 0
        2341           0 :        do ncv = 1, ncvmax
        2342           0 :           ktop(i,ncv)  = 0
        2343           0 :           kbase(i,ncv) = 0
        2344             :        end do
        2345             :     end do
        2346             :
        2347             :     ! ------------------------------------------------------ !
        2348             :     ! Find CL regimes starting from the surface going upward !
        2349             :     ! ------------------------------------------------------ !
        2350             :
        2351             :     rimaxentr = 0._r8
        2352             :
        2353           0 :     do i = 1, ncol
        2354             :
        2355           0 :        riex(2:pver) = ri(i,2:pver)
        2356             :
        2357             :        ! Below allows consistent treatment of surface and other interfaces.
        2358             :        ! Simply, if surface buoyancy flux is positive, Ri of surface is set to be negative.
        2359             :
        2360           0 :        riex(pver+1) = rimaxentr - bflxs(i)
        2361             :
        2362           0 :        ncv = 0
        2363           0 :        k   = pver + 1 ! Work upward from surface interface
        2364             :
        2365           0 :        do while ( k .gt. ntop_turb + 1 )
        2366             :
        2367             :         ! Below means that if 'bflxs > 0' (do not contain '=' sign), surface
        2368             :         ! interface is energetically interior surface.
        2369             :
        2370           0 :           if( riex(k) .lt. rimaxentr ) then
        2371             :
        2372             :               ! Identify a new CL
        2373             :
        2374           0 :               ncv = ncv + 1
        2375             :
        2376             :               ! First define 'kbase' as the first interface below the lower-most unstable interface
        2377             :               ! Thus, Richardson number at 'kbase' is positive.
        2378             :
        2379           0 :               kbase(i,ncv) = min(k+1,pver+1)
        2380             :
        2381             :               ! Decrement k until top unstable level
        2382             :
        2383           0 :               do while( riex(k) .lt. rimaxentr .and. k .gt. ntop_turb + 1 )
        2384           0 :                  k = k - 1
        2385             :               end do
        2386             :
        2387             :               ! ktop is the first interface above upper-most unstable interface
        2388             :               ! Thus, Richardson number at 'ktop' is positive.
        2389             :
        2390           0 :               ktop(i,ncv) = k
        2391             :
        2392             :           else
        2393             :
        2394             :               ! Search upward for a CL.
        2395             :
        2396           0 :               k = k - 1
        2397             :
        2398             :           end if
        2399             :
        2400             :        end do ! End of CL regime finding for each atmospheric column
        2401             :
        2402           0 :        ncvfin(i) = ncv
        2403             :
        2404             :     end do  ! End of atmospheric column do loop
        2405             :
        2406           0 :     return
        2407             :
        2408             :     end subroutine exacol
        2409             :
        2410             :     !============================================================================== !
        2411             :     !                                                                               !
        2412             :     !============================================================================== !
        2413             :
        2414           0 :     subroutine zisocl( pcols  , pver  , long ,                                 &
        2415           0 :                        z      , zi    , n2   ,  s2      ,                      &
        2416           0 :                        bprod  , sprod , bflxs,  tkes    ,                      &
        2417           0 :                        ncvfin , kbase , ktop ,  belongcv,                      &
        2418           0 :                        ricl   , ghcl  , shcl ,  smcl    ,                      &
        2419           0 :                        lbrk   , wbrk  , ebrk ,  extend  , extend_up, extend_dn,&
        2420             :                        errstring)
        2421             :
        2422             :     !------------------------------------------------------------------------ !
        2423             :     ! Object : This 'zisocl' vertically extends original CLs identified from  !
        2424             :     !          'exacol' using a merging test based on either 'wint' or 'l2n2' !
        2425             :     !          and identify new CL regimes. Similar to the case of 'exacol',  !
        2426             :     !          CL regime index increases with height.  After identifying new  !
        2427             :     !          CL regimes ( kbase, ktop, ncvfin ),calculate CL internal mean  !
        2428             :     !          energetics (lbrk : energetic thickness integral, wbrk, ebrk )  !
        2429             :     !          and stability functions (ricl, ghcl, shcl, smcl) by including  !
        2430             :     !          surface interfacial layer contribution when bflxs > 0.   Note  !
        2431             :     !          that there are two options in the treatment of the energetics  !
        2432             :     !          of surface interfacial layer (use_dw_surf= 'true' or 'false')  !
        2433             :     ! Author : Sungsu Park 08/2006, 11/2008                                   !
        2434             :     !------------------------------------------------------------------------ !
        2435             :
        2436             :     implicit none
        2437             :
        2438             :     ! --------------- !
        2439             :     ! Input variables !
        2440             :     ! --------------- !
        2441             :
        2442             :     integer,  intent(in)   :: long                    ! Longitude of the column
        2443             :     integer,  intent(in)   :: pcols                   ! Number of atmospheric columns
        2444             :     integer,  intent(in)   :: pver                    ! Number of atmospheric vertical layers
        2445             :     real(r8), intent(in)   :: z(pcols, pver)          ! Layer mid-point height [ m ]
        2446             :     real(r8), intent(in)   :: zi(pcols, pver+1)       ! Interface height [ m ]
        2447             :     real(r8), intent(in)   :: n2(pcols, pver)         ! Buoyancy frequency at interfaces except surface [ s-2 ]
        2448             :     real(r8), intent(in)   :: s2(pcols, pver)         ! Shear frequency at interfaces except surface [ s-2 ]
        2449             :     real(r8), intent(in)   :: bprod(pcols,pver+1)     ! Buoyancy production [ m2/s3 ]. bprod(i,pver+1) = bflxs
        2450             :     real(r8), intent(in)   :: sprod(pcols,pver+1)     ! Shear production [ m2/s3 ]. sprod(i,pver+1) = usta**3/(vk*z(i,pver))
        2451             :     real(r8), intent(in)   :: bflxs(pcols)            ! Surface buoyancy flux [ m2/s3 ]. bprod(i,pver+1) = bflxs
        2452             :     real(r8), intent(in)   :: tkes(pcols)             ! TKE at the surface [ s2/s2 ]
        2453             :
        2454             :     ! ---------------------- !
        2455             :     ! Input/output variables !
        2456             :     ! ---------------------- !
        2457             :
        2458             :     integer, intent(inout) :: kbase(pcols,ncvmax)     ! Base external interface index of CL
        2459             :     integer, intent(inout) :: ktop(pcols,ncvmax)      ! Top external interface index of CL
        2460             :     integer, intent(inout) :: ncvfin(pcols)           ! Total number of CLs
        2461             :
        2462             :     ! ---------------- !
        2463             :     ! Output variables !
        2464             :     ! ---------------- !
        2465             :
        2466             :     logical,  intent(out) :: belongcv(pcols,pver+1)   ! True if interface is in a CL ( either internal or external )
        2467             :     real(r8), intent(out) :: ricl(pcols,ncvmax)       ! Mean Richardson number of internal CL
        2468             :     real(r8), intent(out) :: ghcl(pcols,ncvmax)       ! Half of normalized buoyancy production of internal CL
        2469             :     real(r8), intent(out) :: shcl(pcols,ncvmax)       ! Galperin instability function of heat-moisture of internal CL
        2470             :     real(r8), intent(out) :: smcl(pcols,ncvmax)       ! Galperin instability function of momentum of internal CL
        2471             :     real(r8), intent(out) :: lbrk(pcols,ncvmax)       ! Thickness of (energetically) internal CL ( lint, [m] )
        2472             :     real(r8), intent(out) :: wbrk(pcols,ncvmax)       ! Mean normalized TKE of internal CL  [ m2/s2 ]
        2473             :     real(r8), intent(out) :: ebrk(pcols,ncvmax)       ! Mean TKE of internal CL ( b1*wbrk, [m2/s2] )
        2474             :
        2475             :     character(len=*), intent(out) :: errstring
        2476             :     ! ------------------ !
        2477             :     ! Internal variables !
        2478             :     ! ------------------ !
        2479             :
        2480             :     logical               :: extend                   ! True when CL is extended in zisocl
        2481             :     logical               :: extend_up                ! True when CL is extended upward in zisocl
        2482             :     logical               :: extend_dn                ! True when CL is extended downward in zisocl
        2483             :     logical               :: bottom                   ! True when CL base is at surface ( kb = pver + 1 )
        2484             :
        2485             :     integer               :: i                        ! Local index for the longitude
        2486             :     integer               :: ncv                      ! CL Index increasing with height
        2487             :     integer               :: incv
        2488             :     integer               :: k
        2489             :     integer               :: kb                       ! Local index for kbase
        2490             :     integer               :: kt                       ! Local index for ktop
        2491             :     integer               :: ncvinit                  ! Value of ncv at routine entrance
        2492             :     integer               :: cntu                     ! Number of merged CLs during upward   extension of individual CL
        2493             :     integer               :: cntd                     ! Number of merged CLs during downward extension of individual CL
        2494             :     integer               :: kbinc                    ! Index for incorporating underlying CL
        2495             :     integer               :: ktinc                    ! Index for incorporating  overlying CL
        2496             :
        2497             :     real(r8)              :: wint                     ! Normalized TKE of internal CL
        2498             :     real(r8)              :: dwinc                    ! Normalized TKE of CL external interfaces
        2499             :     real(r8)              :: dw_surf                  ! Normalized TKE of surface interfacial layer
        2500             :     real(r8)              :: dzinc
        2501             :     real(r8)              :: gh
        2502             :     real(r8)              :: sh
        2503             :     real(r8)              :: sm
        2504             :     real(r8)              :: gh_surf                  ! Half of normalized buoyancy production in surface interfacial layer
        2505             :     real(r8)              :: sh_surf                  ! Galperin instability function in surface interfacial layer
        2506             :     real(r8)              :: sm_surf                  ! Galperin instability function in surface interfacial layer
        2507             :     real(r8)              :: l2n2                     ! Vertical integral of 'l^2N^2' over CL. Include thickness product
        2508             :     real(r8)              :: l2s2                     ! Vertical integral of 'l^2S^2' over CL. Include thickness product
        2509             :     real(r8)              :: dl2n2                    ! Vertical integration of 'l^2*N^2' of CL external interfaces
        2510             :     real(r8)              :: dl2s2                    ! Vertical integration of 'l^2*S^2' of CL external interfaces
        2511             :     real(r8)              :: dl2n2_surf               ! 'dl2n2' defined in the surface interfacial layer
        2512             :     real(r8)              :: dl2s2_surf               ! 'dl2s2' defined in the surface interfacial layer
        2513             :     real(r8)              :: lint                     ! Thickness of (energetically) internal CL
        2514             :     real(r8)              :: dlint                    ! Interfacial layer thickness of CL external interfaces
        2515             :     real(r8)              :: dlint_surf               ! Surface interfacial layer thickness
        2516             :     real(r8)              :: lbulk                    ! Master Length Scale : Whole CL thickness from top to base external interface
        2517             :     real(r8)              :: lz                       ! Turbulent length scale
        2518             :     real(r8)              :: ricll                    ! Mean Richardson number of internal CL
        2519             :     real(r8)              :: trma
        2520             :     real(r8)              :: trmb
        2521             :     real(r8)              :: trmc
        2522             :     real(r8)              :: det
        2523             :     real(r8)              :: zbot                     ! Height of CL base
        2524             :     real(r8)              :: l2rat                    ! Square of ratio of actual to initial CL (not used)
        2525             :     real(r8)              :: gg                       ! Intermediate variable used for calculating stability functions of SBCL
        2526             :     real(r8)              :: tunlramp                 ! Ramping tunl
        2527             :
        2528             :     ! ----------------------- !
        2529             :     ! Main Computation Begins !
        2530             :     ! ----------------------- !
        2531             :
        2532           0 :     i = long
        2533             :
        2534             :     ! Initialize main output variables
        2535             :
        2536           0 :     do k = 1, ncvmax
        2537           0 :        ricl(i,k) = 0._r8
        2538           0 :        ghcl(i,k) = 0._r8
        2539           0 :        shcl(i,k) = 0._r8
        2540           0 :        smcl(i,k) = 0._r8
        2541           0 :        lbrk(i,k) = 0._r8
        2542           0 :        wbrk(i,k) = 0._r8
        2543           0 :        ebrk(i,k) = 0._r8
        2544             :     end do
        2545           0 :     extend    = .false.
        2546           0 :     extend_up = .false.
        2547           0 :     extend_dn = .false.
        2548             :
        2549             :     ! ----------------------------------------------------------- !
        2550             :     ! Loop over each CL to see if any of them need to be extended !
        2551             :     ! ----------------------------------------------------------- !
        2552             :
        2553           0 :     ncv = 1
        2554             :
        2555           0 :     do while( ncv .le. ncvfin(i) )
        2556             :
        2557           0 :        ncvinit = ncv
        2558           0 :        cntu    = 0
        2559           0 :        cntd    = 0
        2560           0 :        kb      = kbase(i,ncv)
        2561           0 :        kt      = ktop(i,ncv)
        2562             :
        2563             :        ! ---------------------------------------------------------------------------- !
        2564             :        ! Calculation of CL interior energetics including surface before extension     !
        2565             :        ! ---------------------------------------------------------------------------- !
        2566             :        ! Note that the contribution of interior interfaces (not surface) to 'wint' is !
        2567             :        ! accounted by using '-sh*l2n2 + sm*l2s2' while the contribution of surface is !
        2568             :        ! accounted by using 'dwsurf = tkes/b1' when bflxs > 0. This approach is fully !
        2569             :        ! reasonable. Another possible alternative,  which seems to be also consistent !
        2570             :        ! is to calculate 'dl2n2_surf'  and  'dl2s2_surf' of surface interfacial layer !
        2571             :        ! separately, and this contribution is explicitly added by initializing 'l2n2' !
        2572             :        ! 'l2s2' not by zero, but by 'dl2n2_surf' and 'ds2n2_surf' below.  At the same !
        2573             :        ! time, 'dwsurf' should be excluded in 'wint' calculation below. The only diff.!
        2574             :        ! between two approaches is that in case of the latter approach, contributions !
        2575             :        ! of surface interfacial layer to the CL mean stability function (ri,gh,sh,sm) !
        2576             :        ! are explicitly included while the first approach is not. In this sense,  the !
        2577             :        ! second approach seems to be more conceptually consistent,   but currently, I !
        2578             :        ! (Sungsu) will keep the first default approach. There is a switch             !
        2579             :        ! 'use_dw_surf' at the first part of eddy_diff.F90 chosing one of              !
        2580             :        ! these two options.                                                           !
        2581             :        ! ---------------------------------------------------------------------------- !
        2582             :
        2583             :        ! ------------------------------------------------------ !
        2584             :        ! Step 0: Calculate surface interfacial layer energetics !
        2585             :        ! ------------------------------------------------------ !
        2586             :
        2587           0 :        lbulk      = zi(i,kt) - zi(i,kb)
        2588           0 :        lbulk      = min( lbulk, lbulk_max )
        2589           0 :        dlint_surf = 0._r8
        2590           0 :        dl2n2_surf = 0._r8
        2591           0 :        dl2s2_surf = 0._r8
        2592           0 :        dw_surf    = 0._r8
        2593           0 :        if( kb .eq. pver+1 ) then
        2594             :
        2595           0 :            if( bflxs(i) .gt. 0._r8 ) then
        2596             :
        2597             :                ! Calculate stability functions of surface interfacial layer
        2598             :                ! from the given 'bprod(i,pver+1)' and 'sprod(i,pver+1)' using
        2599             :                ! inverse approach. Since alph5>0 and alph3<0, denominator of
        2600             :                ! gg is always positive if bprod(i,pver+1)>0.
        2601             :
        2602           0 :                gg    = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
        2603           0 :                gh    = gg/(alph5-gg*alph3)
        2604             :              ! gh    = min(max(gh,-0.28_r8),0.0233_r8)
        2605           0 :                gh    = min(max(gh,-3.5334_r8),0.0233_r8)
        2606           0 :                sh    = alph5/(1._r8+alph3*gh)
        2607           0 :                sm    = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
        2608           0 :                ricll = min(-(sm/sh)*(bprod(i,pver+1)/sprod(i,pver+1)),ricrit)
        2609             :
        2610             :                ! Calculate surface interfacial layer contribution to CL internal
        2611             :                ! energetics. By construction, 'dw_surf = -dl2n2_surf + ds2n2_surf'
        2612             :                ! is exactly satisfied, which corresponds to assuming turbulent
        2613             :                ! length scale of surface interfacial layer = vk * z(i,pver). Note
        2614             :                ! 'dl2n2_surf','dl2s2_surf','dw_surf' include thickness product.
        2615             :
        2616           0 :                dlint_surf = z(i,pver)
        2617           0 :                dl2n2_surf = -vk*(z(i,pver)**2)*bprod(i,pver+1)/(sh*sqrt(tkes(i)))
        2618           0 :                dl2s2_surf =  vk*(z(i,pver)**2)*sprod(i,pver+1)/(sm*sqrt(tkes(i)))
        2619           0 :                dw_surf    = (tkes(i)/b1)*z(i,pver)
        2620             :
        2621             :            else
        2622             :
        2623             :                ! Note that this case can happen when surface is an external
        2624             :                ! interface of CL.
        2625           0 :                lbulk = zi(i,kt) - z(i,pver)
        2626           0 :                lbulk = min( lbulk, lbulk_max )
        2627             :
        2628             :            end if
        2629             :
        2630             :        end if
        2631             :
        2632             :        ! ------------------------------------------------------ !
        2633             :        ! Step 1: Include surface interfacial layer contribution !
        2634             :        ! ------------------------------------------------------ !
        2635             :
        2636           0 :        lint = dlint_surf
        2637             :        l2n2 = dl2n2_surf
        2638             :        l2s2 = dl2s2_surf
        2639           0 :        wint = dw_surf
        2640             :        if( use_dw_surf ) then
        2641           0 :            l2n2 = 0._r8
        2642           0 :            l2s2 = 0._r8
        2643             :        else
        2644             :            wint = 0._r8
        2645             :        end if
        2646             :
        2647             :        ! --------------------------------------------------------------------------------- !
        2648             :        ! Step 2. Include the contribution of 'pure internal interfaces' other than surface !
        2649             :        ! --------------------------------------------------------------------------------- !
        2650             :
        2651           0 :        if( kt .lt. kb - 1 ) then ! The case of non-SBCL.
        2652             :
        2653           0 :            do k = kb - 1, kt + 1, -1
        2654           0 :               if( choice_tunl .eq. 'rampcl' ) then
        2655             :                 ! Modification : I simply used the average tunlramp between the two limits.
        2656             :                   tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
        2657             :               elseif( choice_tunl .eq. 'rampsl' ) then
        2658             :                   tunlramp = ctunl*tunl
        2659             :                 ! tunlramp = 0.765_r8
        2660             :               else
        2661             :                   tunlramp = tunl
        2662             :               endif
        2663             :               if( choice_leng .eq. 'origin' ) then
        2664           0 :                   lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        2665             :                 ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
        2666             :               else
        2667             :                   lz = min( vk*zi(i,k), tunlramp*lbulk )
        2668             :               endif
        2669           0 :               lz = min(leng_max(k), lz)
        2670           0 :               dzinc = z(i,k-1) - z(i,k)
        2671           0 :               l2n2  = l2n2 + lz*lz*n2(i,k)*dzinc
        2672           0 :               l2s2  = l2s2 + lz*lz*s2(i,k)*dzinc
        2673           0 :               lint  = lint + dzinc
        2674             :            end do
        2675             :
        2676             :            ! Calculate initial CL stability functions (gh,sh,sm) and net
        2677             :            ! internal energy of CL including surface contribution if any.
        2678             :
        2679             :          ! Modification : It seems that below cannot be applied when ricrit > 0.19.
        2680             :          !                May need future generalization.
        2681             :
        2682           0 :            ricll = min(l2n2/max(l2s2,ntzero),ricrit) ! Mean Ri of internal CL
        2683           0 :            trma  = alph3*alph4*ricll+2._r8*b1*(alph2-alph4*alph5*ricll)
        2684           0 :            trmb  = ricll*(alph3+alph4)+2._r8*b1*(-alph5*ricll+alph1)
        2685           0 :            trmc  = ricll
        2686           0 :            det   = max(trmb*trmb-4._r8*trma*trmc,0._r8)
        2687           0 :            gh    = (-trmb + sqrt(det))/2._r8/trma
        2688             :          ! gh    = min(max(gh,-0.28_r8),0.0233_r8)
        2689           0 :            gh    = min(max(gh,-3.5334_r8),0.0233_r8)
        2690           0 :            sh    = alph5/(1._r8+alph3*gh)
        2691           0 :            sm    = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
        2692           0 :            wint  = wint - sh*l2n2 + sm*l2s2
        2693             :
        2694             :        else ! The case of SBCL
        2695             :
        2696             :            ! If there is no pure internal interface, use only surface interfacial
        2697             :            ! values. However, re-set surface interfacial values such  that it can
        2698             :            ! be used in the merging tests (either based on 'wint' or 'l2n2')  and
        2699             :            ! in such that surface interfacial energy is not double-counted.
        2700             :            ! Note that regardless of the choise of 'use_dw_surf', below should be
        2701             :            ! kept as it is below, for consistent merging test of extending SBCL.
        2702             :
        2703             :            lint = dlint_surf
        2704             :            l2n2 = dl2n2_surf
        2705             :            l2s2 = dl2s2_surf
        2706             :            wint = dw_surf
        2707             :
        2708             :            ! Aug.29.2006 : Only for the purpose of merging test of extending SRCL
        2709             :            ! based on 'l2n2', re-define 'l2n2' of surface interfacial layer using
        2710             :            ! 'wint'. This part is designed for similar treatment of merging as in
        2711             :            ! the original 'eddy_diff.F90' code,  where 'l2n2' of SBCL was defined
        2712             :            ! as 'l2n2 = - wint / sh'. Note that below block is used only when (1)
        2713             :            ! surface buoyancy production 'bprod(i,pver+1)' is NOT included in the
        2714             :            ! calculation of surface TKE in the initialization of 'bprod(i,pver+1)'
        2715             :            ! in the main subroutine ( even though bflxs > 0 ), and (2) to force
        2716             :            ! current scheme be similar to the previous scheme in the treatment of
        2717             :            ! extending-merging test of SBCL based on 'l2n2'. Otherwise below line
        2718             :            ! must be commented out. Note at this stage, correct non-zero value of
        2719             :            ! 'sh' has been already computed.
        2720             :
        2721             :            if( choice_tkes .eq. 'ebprod' ) then
        2722             :                l2n2 = - wint / sh
        2723             :            endif
        2724             :
        2725             :        endif
        2726             :
        2727             :        ! Set consistent upper limits on 'l2n2' and 'l2s2'. Below limits are
        2728             :        ! reasonable since l2n2 of CL interior interface is always negative.
        2729             :
        2730           0 :        l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
        2731           0 :        l2s2 =  min( l2s2, tkemax*lint/(b1*sm))
        2732             :
        2733             :        ! Note that at this stage, ( gh, sh, sm )  are the values of surface
        2734             :        ! interfacial layer if there is no pure internal interface, while if
        2735             :        ! there is pure internal interface, ( gh, sh, sm ) are the values of
        2736             :        ! pure CL interfaces or the values that include both the CL internal
        2737             :        ! interfaces and surface interfaces, depending on the 'use_dw_surf'.
        2738             :
        2739             :        ! ----------------------------------------------------------------------- !
        2740             :        ! Perform vertical extension-merging process                              !
        2741             :        ! ----------------------------------------------------------------------- !
        2742             :        ! During the merging process, we assumed ( lbulk, sh, sm ) of CL external !
        2743             :        ! interfaces are the same as the ones of the original merging CL. This is !
        2744             :        ! an inevitable approximation since we don't know  ( sh, sm ) of external !
        2745             :        ! interfaces at this stage.     Note that current default merging test is !
        2746             :        ! purely based on buoyancy production without including shear production, !
        2747             :        ! since we used 'l2n2' instead of 'wint' as a merging parameter. However, !
        2748             :        ! merging test based on 'wint' maybe conceptually more attractable.       !
        2749             :        ! Downward CL merging process is identical to the upward merging process, !
        2750             :        ! but when the base of extended CL reaches to the surface, surface inter  !
        2751             :        ! facial layer contribution to the energetic of extended CL must be done  !
        2752             :        ! carefully depending on the sign of surface buoyancy flux. The contribu  !
        2753             :        ! tion of surface interfacial layer energetic is included to the internal !
        2754             :        ! energetics of merging CL only when bflxs > 0.                           !
        2755             :        ! ----------------------------------------------------------------------- !
        2756             :
        2757             :        ! ---------------------------- !
        2758             :        ! Step 1. Extend the CL upward !
        2759             :        ! ---------------------------- !
        2760             :
        2761           0 :        extend = .false.    ! This will become .true. if CL top or base is extended
        2762             :
        2763             :        ! Calculate contribution of potentially incorporable CL top interface
        2764             :
        2765           0 :        if( choice_tunl .eq. 'rampcl' ) then
        2766             :            tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
        2767             :        elseif( choice_tunl .eq. 'rampsl' ) then
        2768             :            tunlramp = ctunl*tunl
        2769             :          ! tunlramp = 0.765_r8
        2770             :        else
        2771             :            tunlramp = tunl
        2772             :        endif
        2773             :        if( choice_leng .eq. 'origin' ) then
        2774           0 :            lz = ( (vk*zi(i,kt))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        2775             :          ! lz = vk*zi(i,kt) / (1._r8+vk*zi(i,kt)/(tunlramp*lbulk))
        2776             :        else
        2777             :            lz = min( vk*zi(i,kt), tunlramp*lbulk )
        2778             :        endif
        2779           0 :        lz = min(leng_max(kt), lz)
        2780             :
        2781           0 :        dzinc = z(i,kt-1)-z(i,kt)
        2782           0 :        dl2n2 = lz*lz*n2(i,kt)*dzinc
        2783           0 :        dl2s2 = lz*lz*s2(i,kt)*dzinc
        2784           0 :        dwinc = -sh*dl2n2 + sm*dl2s2
        2785             :
        2786             :        ! ------------ !
        2787             :        ! Merging Test !
        2788             :        ! ------------ !
        2789             :
        2790             :        ! The part of the below test that involves kt and z has different
        2791             :        ! effects based on the model top.
        2792             :        ! If the model top is in the stratosphere, we want the loop to
        2793             :        ! continue until it either completes normally, or kt is pushed to
        2794             :        ! the top of the model. The latter case should not happen, so this
        2795             :        ! causes an error.
        2796             :        ! If the model top is higher, as in WACCM and WACCM-X, if kt is
        2797             :        ! pushed close to the model top, this may not represent an error at
        2798             :        ! all, because of very different and more variable temperature/wind
        2799             :        ! profiles at the model top. Therefore we simply exit the loop early
        2800             :        ! and continue with no errors.
        2801             :
        2802             :      ! do while (  dwinc .gt. ( rinc*dzinc*wint/(lint+(1._r8-rinc)*dzinc)) )  ! Merging test based on wint
        2803             :      ! do while ( -dl2n2 .gt. (-rinc*dzinc*l2n2/(lint+(1._r8-rinc)*dzinc)) )  ! Merging test based on l2n2
        2804             :        do while ( -dl2n2 .gt. (-rinc*l2n2/(1._r8-rinc)) &                     ! Integral merging test
        2805           0 :             .and. (kt > ntop_turb+2 .or. z(i,kt) < 50000._r8) )
        2806             :
        2807             :           ! Add contribution of top external interface to interior energy.
        2808             :           ! Note even when we chose 'use_dw_surf='true.', the contribution
        2809             :           ! of surface interfacial layer to 'l2n2' and 'l2s2' are included
        2810             :           ! here. However it is not double counting of surface interfacial
        2811             :           ! energy : surface interfacial layer energy is counted in 'wint'
        2812             :           ! formula and 'l2n2' is just used for performing merging test in
        2813             :           ! this 'do while' loop.
        2814             :
        2815           0 :           lint = lint + dzinc
        2816           0 :           l2n2 = l2n2 + dl2n2
        2817           0 :           l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
        2818           0 :           l2s2 = l2s2 + dl2s2
        2819           0 :           wint = wint + dwinc
        2820             :
        2821             :           ! Extend top external interface of CL upward after merging
        2822             :
        2823           0 :           kt        = kt - 1
        2824           0 :           extend    = .true.
        2825           0 :           extend_up = .true.
        2826           0 :           if( kt .eq. ntop_turb ) then
        2827           0 :              errstring = 'zisocl: Error: Tried to extend CL to the model top'
        2828           0 :              return
        2829             :           end if
        2830             :
        2831             :           ! If the top external interface of extending CL is the same as the
        2832             :           ! top interior interface of the overlying CL, overlying CL will be
        2833             :           ! automatically merged. Then,reduce total number of CL regime by 1.
        2834             :           ! and increase 'cntu'(number of merged CLs during upward extension)
        2835             :           ! by 1.
        2836             :
        2837           0 :           ktinc = kbase(i,ncv+cntu+1) - 1  ! Lowest interior interface of overlying CL
        2838             :
        2839           0 :           if( kt .eq. ktinc ) then
        2840             :
        2841           0 :               do k = kbase(i,ncv+cntu+1) - 1, ktop(i,ncv+cntu+1) + 1, -1
        2842             :
        2843           0 :                  if( choice_tunl .eq. 'rampcl' ) then
        2844             :                      tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
        2845             :                  elseif( choice_tunl .eq. 'rampsl' ) then
        2846             :                      tunlramp = ctunl*tunl
        2847             :                    ! tunlramp = 0.765_r8
        2848             :                  else
        2849             :                      tunlramp = tunl
        2850             :                  endif
        2851             :                  if( choice_leng .eq. 'origin' ) then
        2852           0 :                      lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        2853             :                    ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
        2854             :                  else
        2855             :                      lz = min( vk*zi(i,k), tunlramp*lbulk )
        2856             :                  endif
        2857           0 :                  lz = min(leng_max(k), lz)
        2858             :
        2859           0 :                  dzinc = z(i,k-1)-z(i,k)
        2860           0 :                  dl2n2 = lz*lz*n2(i,k)*dzinc
        2861           0 :                  dl2s2 = lz*lz*s2(i,k)*dzinc
        2862           0 :                  dwinc = -sh*dl2n2 + sm*dl2s2
        2863             :
        2864           0 :                  lint = lint + dzinc
        2865           0 :                  l2n2 = l2n2 + dl2n2
        2866           0 :                  l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
        2867           0 :                  l2s2 = l2s2 + dl2s2
        2868           0 :                  wint = wint + dwinc
        2869             :
        2870             :               end do
        2871             :
        2872           0 :               kt        = ktop(i,ncv+cntu+1)
        2873           0 :               ncvfin(i) = ncvfin(i) - 1
        2874           0 :               cntu      = cntu + 1
        2875             :
        2876             :           end if
        2877             :
        2878             :           ! Again, calculate the contribution of potentially incorporatable CL
        2879             :           ! top external interface of CL regime.
        2880             :
        2881           0 :           if( choice_tunl .eq. 'rampcl' ) then
        2882             :               tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
        2883             :           elseif( choice_tunl .eq. 'rampsl' ) then
        2884             :               tunlramp = ctunl*tunl
        2885             :             ! tunlramp = 0.765_r8
        2886             :           else
        2887             :               tunlramp = tunl
        2888             :           endif
        2889             :           if( choice_leng .eq. 'origin' ) then
        2890           0 :               lz = ( (vk*zi(i,kt))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        2891             :             ! lz = vk*zi(i,kt) / (1._r8+vk*zi(i,kt)/(tunlramp*lbulk))
        2892             :           else
        2893             :               lz = min( vk*zi(i,kt), tunlramp*lbulk )
        2894             :           endif
        2895           0 :           lz = min(leng_max(kt), lz)
        2896             :
        2897           0 :           dzinc = z(i,kt-1)-z(i,kt)
        2898           0 :           dl2n2 = lz*lz*n2(i,kt)*dzinc
        2899           0 :           dl2s2 = lz*lz*s2(i,kt)*dzinc
        2900           0 :           dwinc = -sh*dl2n2 + sm*dl2s2
        2901             :
        2902             :        end do   ! End of upward merging test 'do while' loop
        2903             :
        2904             :        ! Update CL interface indices appropriately if any CL was merged.
        2905             :        ! Note that below only updated the interface index of merged CL,
        2906             :        ! not the original merging CL.  Updates of 'kbase' and 'ktop' of
        2907             :        ! the original merging CL  will be done after finishing downward
        2908             :        ! extension also later.
        2909             :
        2910           0 :        if( cntu .gt. 0 ) then
        2911           0 :            do incv = 1, ncvfin(i) - ncv
        2912           0 :               kbase(i,ncv+incv) = kbase(i,ncv+cntu+incv)
        2913           0 :               ktop(i,ncv+incv)  = ktop(i,ncv+cntu+incv)
        2914             :            end do
        2915             :        end if
        2916             :
        2917             :        ! ------------------------------ !
        2918             :        ! Step 2. Extend the CL downward !
        2919             :        ! ------------------------------ !
        2920             :
        2921           0 :        if( kb .ne. pver + 1 ) then
        2922             :
        2923             :            ! Calculate contribution of potentially incorporable CL base interface
        2924             :
        2925           0 :            if( choice_tunl .eq. 'rampcl' ) then
        2926             :                tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
        2927             :            elseif( choice_tunl .eq. 'rampsl' ) then
        2928             :                tunlramp = ctunl*tunl
        2929             :              ! tunlramp = 0.765_r8
        2930             :            else
        2931             :                tunlramp = tunl
        2932             :            endif
        2933             :            if( choice_leng .eq. 'origin' ) then
        2934           0 :                lz = ( (vk*zi(i,kb))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        2935             :              ! lz = vk*zi(i,kb) / (1._r8+vk*zi(i,kb)/(tunlramp*lbulk))
        2936             :            else
        2937             :                lz = min( vk*zi(i,kb), tunlramp*lbulk )
        2938             :            endif
        2939           0 :            lz = min(leng_max(kb), lz)
        2940             :
        2941           0 :            dzinc = z(i,kb-1)-z(i,kb)
        2942           0 :            dl2n2 = lz*lz*n2(i,kb)*dzinc
        2943           0 :            dl2s2 = lz*lz*s2(i,kb)*dzinc
        2944           0 :            dwinc = -sh*dl2n2 + sm*dl2s2
        2945             :
        2946             :            ! ------------ !
        2947             :            ! Merging test !
        2948             :            ! ------------ !
        2949             :
        2950             :            ! In the below merging tests, I must keep '.and.(kb.ne.pver+1)',
        2951             :            ! since 'kb' is continuously updated within the 'do while' loop
        2952             :            ! whenever CL base is merged.
        2953             :
        2954             :          ! do while( (  dwinc .gt. ( rinc*dzinc*wint/(lint+(1._r8-rinc)*dzinc)) ) &  ! Merging test based on wint
        2955             :          ! do while( ( -dl2n2 .gt. (-rinc*dzinc*l2n2/(lint+(1._r8-rinc)*dzinc)) ) &  ! Merging test based on l2n2
        2956             :          !             .and.(kb.ne.pver+1))
        2957             :            do while( ( -dl2n2 .gt. (-rinc*l2n2/(1._r8-rinc)) ) &                     ! Integral merging test
        2958           0 :                        .and.(kb.ne.pver+1))
        2959             :
        2960             :               ! Add contributions from interfacial layer kb to CL interior
        2961             :
        2962           0 :               lint = lint + dzinc
        2963           0 :               l2n2 = l2n2 + dl2n2
        2964           0 :               l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
        2965           0 :               l2s2 = l2s2 + dl2s2
        2966           0 :               wint = wint + dwinc
        2967             :
        2968             :               ! Extend the base external interface of CL downward after merging
        2969             :
        2970           0 :               kb        =  kb + 1
        2971           0 :               extend    = .true.
        2972           0 :               extend_dn = .true.
        2973             :
        2974             :               ! If the base external interface of extending CL is the same as the
        2975             :               ! base interior interface of the underlying CL, underlying CL  will
        2976             :               ! be automatically merged. Then, reduce total number of CL by 1.
        2977             :               ! For a consistent treatment with 'upward' extension,  I should use
        2978             :               ! 'kbinc = kbase(i,ncv-1) - 1' instead of 'ktop(i,ncv-1) + 1' below.
        2979             :               ! However, it seems that these two methods produce the same results.
        2980             :               ! Note also that in contrast to upward merging, the decrease of ncv
        2981             :               ! should be performed here.
        2982             :               ! Note that below formula correctly works even when upperlying CL
        2983             :               ! regime incorporates below SBCL.
        2984             :
        2985           0 :               kbinc = 0
        2986           0 :               if( ncv .gt. 1 ) kbinc = ktop(i,ncv-1) + 1
        2987           0 :               if( kb .eq. kbinc ) then
        2988             :
        2989           0 :                   do k =  ktop(i,ncv-1) + 1, kbase(i,ncv-1) - 1
        2990             :
        2991           0 :                      if( choice_tunl .eq. 'rampcl' ) then
        2992             :                          tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
        2993             :                      elseif( choice_tunl .eq. 'rampsl' ) then
        2994             :                          tunlramp = ctunl*tunl
        2995             :                        ! tunlramp = 0.765_r8
        2996             :                      else
        2997             :                          tunlramp = tunl
        2998             :                      endif
        2999             :                      if( choice_leng .eq. 'origin' ) then
        3000           0 :                          lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        3001             :                        ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
        3002             :                      else
        3003             :                          lz = min( vk*zi(i,k), tunlramp*lbulk )
        3004             :                      endif
        3005           0 :                      lz = min(leng_max(k), lz)
        3006             :
        3007           0 :                      dzinc = z(i,k-1)-z(i,k)
        3008           0 :                      dl2n2 = lz*lz*n2(i,k)*dzinc
        3009           0 :                      dl2s2 = lz*lz*s2(i,k)*dzinc
        3010           0 :                      dwinc = -sh*dl2n2 + sm*dl2s2
        3011             :
        3012           0 :                      lint = lint + dzinc
        3013           0 :                      l2n2 = l2n2 + dl2n2
        3014           0 :                      l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
        3015           0 :                      l2s2 = l2s2 + dl2s2
        3016           0 :                      wint = wint + dwinc
        3017             :
        3018             :                   end do
        3019             :
        3020             :                   ! We are incorporating interior of CL ncv-1, so merge
        3021             :                   ! this CL into the current CL.
        3022             :
        3023           0 :                   kb        = kbase(i,ncv-1)
        3024           0 :                   ncv       = ncv - 1
        3025           0 :                   ncvfin(i) = ncvfin(i) -1
        3026           0 :                   cntd      = cntd + 1
        3027             :
        3028             :               end if
        3029             :
        3030             :               ! Calculate the contribution of potentially incorporatable CL
        3031             :               ! base external interface. Calculate separately when the base
        3032             :               ! of extended CL is surface and non-surface.
        3033             :
        3034           0 :               if( kb .eq. pver + 1 ) then
        3035             :
        3036           0 :                   if( bflxs(i) .gt. 0._r8 ) then
        3037             :                       ! Calculate stability functions of surface interfacial layer
        3038           0 :                       gg = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
        3039           0 :                       gh_surf = gg/(alph5-gg*alph3)
        3040             :                     ! gh_surf = min(max(gh_surf,-0.28_r8),0.0233_r8)
        3041           0 :                       gh_surf = min(max(gh_surf,-3.5334_r8),0.0233_r8)
        3042           0 :                       sh_surf = alph5/(1._r8+alph3*gh_surf)
        3043           0 :                       sm_surf = (alph1 + alph2*gh_surf)/(1._r8+alph3*gh_surf)/(1._r8+alph4*gh_surf)
        3044             :                       ! Calculate surface interfacial layer contribution. By construction,
        3045             :                       ! it exactly becomes 'dw_surf = -dl2n2_surf + ds2n2_surf'
        3046           0 :                       dlint_surf = z(i,pver)
        3047           0 :                       dl2n2_surf = -vk*(z(i,pver)**2._r8)*bprod(i,pver+1)/(sh_surf*sqrt(tkes(i)))
        3048           0 :                       dl2s2_surf =  vk*(z(i,pver)**2._r8)*sprod(i,pver+1)/(sm_surf*sqrt(tkes(i)))
        3049           0 :                       dw_surf = (tkes(i)/b1)*z(i,pver)
        3050             :                   else
        3051             :                       dlint_surf = 0._r8
        3052             :                       dl2n2_surf = 0._r8
        3053             :                       dl2s2_surf = 0._r8
        3054             :                       dw_surf = 0._r8
        3055             :                   end if
        3056             :                   ! If (kb.eq.pver+1), updating of CL internal energetics should be
        3057             :                   ! performed here inside of 'do while' loop, since 'do while' loop
        3058             :                   ! contains the constraint of '.and.(kb.ne.pver+1)',so updating of
        3059             :                   ! CL internal energetics cannot be performed within this do while
        3060             :                   ! loop when kb.eq.pver+1. Even though I updated all 'l2n2','l2s2',
        3061             :                   ! 'wint' below, only the updated 'wint' is used in the following
        3062             :                   ! numerical calculation.
        3063           0 :                   lint = lint + dlint_surf
        3064           0 :                   l2n2 = l2n2 + dl2n2_surf
        3065           0 :                   l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
        3066           0 :                   l2s2 = l2s2 + dl2s2_surf
        3067           0 :                   wint = wint + dw_surf
        3068             :
        3069             :               else
        3070             :
        3071           0 :                   if( choice_tunl .eq. 'rampcl' ) then
        3072             :                       tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
        3073             :                   elseif( choice_tunl .eq. 'rampsl' ) then
        3074             :                       tunlramp = ctunl*tunl
        3075             :                     ! tunlramp = 0.765_r8
        3076             :                   else
        3077             :                       tunlramp = tunl
        3078             :                   endif
        3079             :                   if( choice_leng .eq. 'origin' ) then
        3080           0 :                       lz = ( (vk*zi(i,kb))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        3081             :                     ! lz = vk*zi(i,kb) / (1._r8+vk*zi(i,kb)/(tunlramp*lbulk))
        3082             :                   else
        3083             :                       lz = min( vk*zi(i,kb), tunlramp*lbulk )
        3084             :                   endif
        3085           0 :                   lz = min(leng_max(kb), lz)
        3086             :
        3087           0 :                   dzinc = z(i,kb-1)-z(i,kb)
        3088           0 :                   dl2n2 = lz*lz*n2(i,kb)*dzinc
        3089           0 :                   dl2s2 = lz*lz*s2(i,kb)*dzinc
        3090           0 :                   dwinc = -sh*dl2n2 + sm*dl2s2
        3091             :
        3092             :               end if
        3093             :
        3094             :           end do ! End of merging test 'do while' loop
        3095             :
        3096           0 :           if( (kb.eq.pver+1) .and. (ncv.ne.1) ) then
        3097           0 :                errstring = 'Major mistake zisocl: the CL based at surface is not indexed 1'
        3098           0 :                return
        3099             :           end if
        3100             :
        3101             :        end if   ! Done with bottom extension of CL
        3102             :
        3103             :        ! Update CL interface indices appropriately if any CL was merged.
        3104             :        ! Note that below only updated the interface index of merged CL,
        3105             :        ! not the original merging CL.  Updates of 'kbase' and 'ktop' of
        3106             :        ! the original merging CL  will be done later below. I should
        3107             :        ! check in detail if below index updating is correct or not.
        3108             :
        3109           0 :        if( cntd .gt. 0 ) then
        3110           0 :            do incv = 1, ncvfin(i) - ncv
        3111           0 :               kbase(i,ncv+incv) = kbase(i,ncvinit+incv)
        3112           0 :               ktop(i,ncv+incv)  = ktop(i,ncvinit+incv)
        3113             :            end do
        3114             :        end if
        3115             :
        3116             :        ! Sanity check for positive wint.
        3117             :
        3118           0 :        if( wint .lt. 0.01_r8 ) then
        3119           0 :            wint = 0.01_r8
        3120             :        end if
        3121             :
        3122             :        ! -------------------------------------------------------------------------- !
        3123             :        ! Finally update CL mean internal energetics including surface contribution  !
        3124             :        ! after finishing all the CL extension-merging process.  As mentioned above, !
        3125             :        ! there are two possible ways in the treatment of surface interfacial layer, !
        3126             :        ! either through 'dw_surf' or 'dl2n2_surf and dl2s2_surf' by setting logical !
        3127             :        ! variable 'use_dw_surf' =.true. or .false.    In any cases, we should avoid !
        3128             :        ! double counting of surface interfacial layer and one single consistent way !
        3129             :        ! should be used throughout the program.                                     !
        3130             :        ! -------------------------------------------------------------------------- !
        3131             :
        3132           0 :        if( extend ) then
        3133             :
        3134           0 :            ktop(i,ncv)  = kt
        3135           0 :            kbase(i,ncv) = kb
        3136             :
        3137             :            ! ------------------------------------------------------ !
        3138             :            ! Step 1: Include surface interfacial layer contribution !
        3139             :            ! ------------------------------------------------------ !
        3140             :
        3141           0 :            lbulk      = zi(i,kt) - zi(i,kb)
        3142           0 :            lbulk      = min( lbulk, lbulk_max )
        3143           0 :            dlint_surf = 0._r8
        3144           0 :            dl2n2_surf = 0._r8
        3145           0 :            dl2s2_surf = 0._r8
        3146           0 :            dw_surf    = 0._r8
        3147           0 :            if( kb .eq. pver + 1 ) then
        3148           0 :                if( bflxs(i) .gt. 0._r8 ) then
        3149             :                    ! Calculate stability functions of surface interfacial layer
        3150           0 :                    gg = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
        3151           0 :                    gh = gg/(alph5-gg*alph3)
        3152             :                  ! gh = min(max(gh,-0.28_r8),0.0233_r8)
        3153           0 :                    gh = min(max(gh,-3.5334_r8),0.0233_r8)
        3154           0 :                    sh = alph5/(1._r8+alph3*gh)
        3155           0 :                    sm = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
        3156             :                    ! Calculate surface interfacial layer contribution. By construction,
        3157             :                    ! it exactly becomes 'dw_surf = -dl2n2_surf + ds2n2_surf'
        3158           0 :                    dlint_surf = z(i,pver)
        3159           0 :                    dl2n2_surf = -vk*(z(i,pver)**2._r8)*bprod(i,pver+1)/(sh*sqrt(tkes(i)))
        3160           0 :                    dl2s2_surf =  vk*(z(i,pver)**2._r8)*sprod(i,pver+1)/(sm*sqrt(tkes(i)))
        3161           0 :                    dw_surf    = (tkes(i)/b1)*z(i,pver)
        3162             :                else
        3163           0 :                    lbulk = zi(i,kt) - z(i,pver)
        3164           0 :                    lbulk = min( lbulk, lbulk_max )
        3165             :                end if
        3166             :            end if
        3167           0 :            lint = dlint_surf
        3168             :            l2n2 = dl2n2_surf
        3169             :            l2s2 = dl2s2_surf
        3170           0 :            wint = dw_surf
        3171             :            if( use_dw_surf ) then
        3172           0 :                l2n2 = 0._r8
        3173           0 :                l2s2 = 0._r8
        3174             :            else
        3175             :                wint = 0._r8
        3176             :            end if
        3177             :
        3178             :            ! -------------------------------------------------------------- !
        3179             :            ! Step 2. Include the contribution of 'pure internal interfaces' !
        3180             :            ! -------------------------------------------------------------- !
        3181             :
        3182           0 :            do k = kt + 1, kb - 1
        3183           0 :               if( choice_tunl .eq. 'rampcl' ) then
        3184             :                   tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
        3185             :               elseif( choice_tunl .eq. 'rampsl' ) then
        3186             :                   tunlramp = ctunl*tunl
        3187             :                 ! tunlramp = 0.765_r8
        3188             :               else
        3189             :                   tunlramp = tunl
        3190             :               endif
        3191             :               if( choice_leng .eq. 'origin' ) then
        3192           0 :                   lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
        3193             :                 ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
        3194             :               else
        3195             :                   lz = min( vk*zi(i,k), tunlramp*lbulk )
        3196             :               endif
        3197           0 :               lz = min(leng_max(k), lz)
        3198           0 :               dzinc = z(i,k-1) - z(i,k)
        3199           0 :               lint = lint + dzinc
        3200           0 :               l2n2 = l2n2 + lz*lz*n2(i,k)*dzinc
        3201           0 :               l2s2 = l2s2 + lz*lz*s2(i,k)*dzinc
        3202             :            end do
        3203             :
        3204           0 :            ricll = min(l2n2/max(l2s2,ntzero),ricrit)
        3205           0 :            trma = alph3*alph4*ricll+2._r8*b1*(alph2-alph4*alph5*ricll)
        3206           0 :            trmb = ricll*(alph3+alph4)+2._r8*b1*(-alph5*ricll+alph1)
        3207           0 :            trmc = ricll
        3208           0 :            det = max(trmb*trmb-4._r8*trma*trmc,0._r8)
        3209           0 :            gh = (-trmb + sqrt(det))/2._r8/trma
        3210             :          ! gh = min(max(gh,-0.28_r8),0.0233_r8)
        3211           0 :            gh = min(max(gh,-3.5334_r8),0.0233_r8)
        3212           0 :            sh = alph5 / (1._r8+alph3*gh)
        3213           0 :            sm = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
        3214             :            ! Even though the 'wint' after finishing merging was positive, it is
        3215             :            ! possible that re-calculated 'wint' here is negative.  In this case,
        3216             :            ! correct 'wint' to be a small positive number
        3217           0 :            wint = max( wint - sh*l2n2 + sm*l2s2, 0.01_r8 )
        3218             :
        3219             :        end if
        3220             :
        3221             :        ! ---------------------------------------------------------------------- !
        3222             :        ! Calculate final output variables of each CL (either has merged or not) !
        3223             :        ! ---------------------------------------------------------------------- !
        3224             :
        3225           0 :        lbrk(i,ncv) = lint
        3226           0 :        wbrk(i,ncv) = wint/lint
        3227           0 :        ebrk(i,ncv) = b1*wbrk(i,ncv)
        3228           0 :        ebrk(i,ncv) = min(ebrk(i,ncv),tkemax)
        3229           0 :        ricl(i,ncv) = ricll
        3230           0 :        ghcl(i,ncv) = gh
        3231           0 :        shcl(i,ncv) = sh
        3232           0 :        smcl(i,ncv) = sm
        3233             :
        3234             :        ! Increment counter for next CL. I should check if the increament of 'ncv'
        3235             :        ! below is reasonable or not, since whenever CL is merged during downward
        3236             :        ! extension process, 'ncv' is lowered down continuously within 'do' loop.
        3237             :        ! But it seems that below 'ncv = ncv + 1' is perfectly correct.
        3238             :
        3239           0 :        ncv = ncv + 1
        3240             :
        3241             :     end do                   ! End of loop over each CL regime, ncv.
        3242             :
        3243             :     ! ---------------------------------------------------------- !
        3244             :     ! Re-initialize external interface indices which are not CLs !
        3245             :     ! ---------------------------------------------------------- !
        3246             :
        3247           0 :     do ncv = ncvfin(i) + 1, ncvmax
        3248           0 :        ktop(i,ncv)  = 0
        3249           0 :        kbase(i,ncv) = 0
        3250             :     end do
        3251             :
        3252             :     ! ------------------------------------------------ !
        3253             :     ! Update CL interface identifiers, 'belongcv'      !
        3254             :     ! CL external interfaces are also identified as CL !
        3255             :     ! ------------------------------------------------ !
        3256             :
        3257           0 :     do k = 1, pver + 1
        3258           0 :        belongcv(i,k) = .false.
        3259             :     end do
        3260             :
        3261           0 :     do ncv = 1, ncvfin(i)
        3262           0 :        do k = ktop(i,ncv), kbase(i,ncv)
        3263           0 :           belongcv(i,k) = .true.
        3264             :        end do
        3265             :     end do
        3266             :
        3267             :     return
        3268             :
        3269           0 :     end subroutine zisocl
        3270             :
        3271           0 :     real(r8) function compute_cubic(a,b,c)
        3272             :     ! ------------------------------------------------------------------------- !
        3273             :     ! Solve canonical cubic : x^3 + a*x^2 + b*x + c = 0,  x = sqrt(e)/sqrt(<e>) !
        3274             :     ! Set x = max(xmin,x) at the end                                            !
        3275             :     ! ------------------------------------------------------------------------- !
        3276             :     implicit none
        3277             :     real(r8), intent(in)     :: a, b, c
        3278             :     real(r8)  qq, rr, dd, theta, aa, bb, x1, x2, x3
        3279             :     real(r8), parameter      :: xmin = 1.e-2_r8
        3280             :
        3281           0 :     qq = (a**2-3._r8*b)/9._r8
        3282           0 :     rr = (2._r8*a**3 - 9._r8*a*b + 27._r8*c)/54._r8
        3283             :
        3284           0 :     dd = rr**2 - qq**3
        3285           0 :     if( dd .le. 0._r8 ) then
        3286           0 :         theta = acos(rr/qq**(3._r8/2._r8))
        3287           0 :         x1 = -2._r8*sqrt(qq)*cos(theta/3._r8) - a/3._r8
        3288           0 :         x2 = -2._r8*sqrt(qq)*cos((theta+2._r8*3.141592_r8)/3._r8) - a/3._r8
        3289           0 :         x3 = -2._r8*sqrt(qq)*cos((theta-2._r8*3.141592_r8)/3._r8) - a/3._r8
        3290           0 :         compute_cubic = max(max(max(x1,x2),x3),xmin)
        3291           0 :         return
        3292             :     else
        3293           0 :         if( rr .ge. 0._r8 ) then
        3294           0 :             aa = -(sqrt(rr**2-qq**3)+rr)**(1._r8/3._r8)
        3295             :         else
        3296           0 :             aa =  (sqrt(rr**2-qq**3)-rr)**(1._r8/3._r8)
        3297             :         endif
        3298           0 :         if( aa .eq. 0._r8 ) then
        3299             :             bb = 0._r8
        3300             :         else
        3301           0 :             bb = qq/aa
        3302             :         endif
        3303           0 :         compute_cubic = max((aa+bb)-a/3._r8,xmin)
        3304           0 :         return
        3305             :     endif
        3306             :
        3307             :     return
        3308             :     end function compute_cubic
        3309             :
        3310             : END MODULE eddy_diff

Generated by: LCOV version 1.14