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

              Line data    Source code
           1             :
           2             :   module cldwat2m_macro
           3             :
           4             :   !--------------------------------------------------- !
           5             :   ! Purpose     : CAM Interface for Cloud Macrophysics !
           6             :   ! Author      : Sungsu Park                          !
           7             :   ! Description : Park et al. 2010.                    !
           8             :   ! For questions, contact Sungsu Park                 !
           9             :   !                        e-mail : sungsup@ucar.edu   !
          10             :   !                        phone  : 303-497-1375       !
          11             :   !--------------------------------------------------- !
          12             :
          13             :    use shr_kind_mod,     only: r8=>shr_kind_r8
          14             :    use spmd_utils,       only: masterproc
          15             :    use ppgrid,           only: pcols, pver, pverp
          16             :    use cam_abortutils,   only: endrun
          17             :    use physconst,        only: cpair, latvap, latice, rh2o, gravit, rair
          18             :    use wv_saturation,    only: qsat_water, svp_water, svp_ice, qsat_ice
          19             :    use cam_logfile,      only: iulog
          20             :    use ref_pres,         only: top_lev=>trop_cloud_top_lev
          21             :    use cldfrc2m,         only: astG_PDF_single, astG_PDF, astG_RHU_single, &
          22             :                                astG_RHU, aist_single, aist_vector,         &
          23             :                                rhmini_const, rhmaxi_const
          24             :
          25             :    implicit none
          26             :    private
          27             :    save
          28             :
          29             :    public ::           &
          30             :       ini_macro,       &
          31             :       mmacro_pcond
          32             :
          33             :    ! -------------- !
          34             :    ! Set Parameters !
          35             :    ! -------------- !
          36             :
          37             :    ! ------------------------------------------------------------------------------- !
          38             :    ! Parameter used for selecting generalized critical RH for liquid and ice stratus !
          39             :    ! ------------------------------------------------------------------------------- !
          40             :
          41             :    integer :: i_rhminl ! This is for liquid stratus fraction.
          42             :                        ! If 0 : Original fixed critical RH from the namelist.
          43             :                        ! If 1 : Add convective detrainment effect on the above '0' option.
          44             :                        !        In this case, 'tau_detw' [s] should be specified below.
          45             :                        ! If 2 : Use fully scale-adaptive method.
          46             :                        !        In this case, 'tau_detw' [s] and 'c_aniso' [no unit] should
          47             :                        !        be specified below.
          48             :
          49             :    integer :: i_rhmini ! This is for ice stratus fraction.
          50             :                        ! If 0 : Original fixed critical RH from the namelist.
          51             :                        ! If 1 : Add convective detrainment effect on the above '0' option.
          52             :                        !        In this case, 'tau_deti' [s] should be specified below.
          53             :                        ! If 2 : Use fully scale-adaptive method.
          54             :                        !        In this case, 'tau_deti' [s] and 'c_aniso' [no unit] should
          55             :                        !        be specified below.
          56             :                        ! Note that 'micro_mg_cam' is using below 'rhmini_const', regardless
          57             :                        ! of 'i_rhmini'.  This connection should be built in future.
          58             :
          59             :    real(r8), parameter :: tau_detw =100._r8   ! Dissipation time scale of convective liquid condensate detrained
          60             :                                               !  into the clear portion. [hr]. 0.5-3 hr is possible.
          61             :    real(r8), parameter :: tau_deti =  1._r8   ! Dissipation time scale of convective ice    condensate detrained
          62             :                                               !  into the clear portion. [hr]. 0.5-3 hr is possible.
          63             :    real(r8), parameter :: c_aniso  =  1._r8   ! Inverse of anisotropic factor of PBL turbulence
          64             :
          65             :    ! ----------------------------- !
          66             :    ! Parameters for Liquid Stratus !
          67             :    ! ----------------------------- !
          68             :
          69             :    logical,  parameter  :: CAMstfrac    = .false.    ! If .true. (.false.),
          70             :                                                      ! use Slingo (triangular PDF-based) liquid stratus fraction
          71             :    real(r8), parameter  :: qlst_min     = 2.e-5_r8   ! Minimum in-stratus LWC constraint [ kg/kg ]
          72             :    real(r8), parameter  :: qlst_max     = 3.e-3_r8   ! Maximum in-stratus LWC constraint [ kg/kg ]
          73             :    real(r8), parameter  :: cc           = 0.1_r8     ! For newly formed/dissipated in-stratus CWC ( 0 <= cc <= 1 )
          74             :    integer,  parameter  :: niter        = 2          ! For iterative computation of QQ with 'ramda' below.
          75             :    real(r8), parameter  :: ramda        = 0.5_r8     ! Explicit : ramda = 0, Implicit : ramda = 1 ( 0<= ramda <= 1 )
          76             :    real(r8), private    :: rhminl_const              ! Critical RH for low-level  liquid stratus clouds
          77             :    real(r8), private    :: rhminl_adj_land_const     ! rhminl adjustment for snowfree land
          78             :    real(r8), private    :: rhminh_const              ! Critical RH for high-level liquid stratus clouds
          79             :    real(r8), private    :: premit                    ! Top    height for mid-level liquid stratus fraction
          80             :    real(r8), private    :: premib                    ! Bottom height for mid-level liquid stratus fraction
          81             :
          82             :    real(r8), parameter :: qsmall = 1.e-18_r8         ! Smallest mixing ratio considered in the macrophysics
          83             :
          84             :    contains
          85             :
          86             :    ! -------------- !
          87             :    ! Initialization !
          88             :    ! -------------- !
          89             :
          90           0 :    subroutine ini_macro(rhminl_opt_in, rhmini_opt_in)
          91             :
          92             :    !--------------------------------------------------------------------- !
          93             :    !                                                                      !
          94             :    ! Purpose: Initialize constants for the liquid stratiform macrophysics !
          95             :    !                                                                      !
          96             :    ! Author:  Sungsu Park, Dec.01.2009.                                   !
          97             :    !                                                                      !
          98             :    !--------------------------------------------------------------------- !
          99             :
         100             :    use cloud_fraction, only: cldfrc_getparams
         101             :    use cam_history,    only: addfld
         102             :
         103             :    integer,  intent(in) :: rhminl_opt_in
         104             :    integer,  intent(in) :: rhmini_opt_in
         105             :
         106           0 :    i_rhminl   = rhminl_opt_in
         107           0 :    i_rhmini   = rhmini_opt_in
         108             :
         109             :    call cldfrc_getparams(rhminl_out=rhminl_const, rhminl_adj_land_out=rhminl_adj_land_const,  &
         110           0 :                          rhminh_out=rhminh_const, premit_out=premit, premib_out=premib)
         111             :
         112           0 :    if( masterproc ) then
         113           0 :        write(iulog,*) 'Park Macrophysics Parameters'
         114           0 :        write(iulog,*) '  rhminl          = ', rhminl_const
         115           0 :        write(iulog,*) '  rhminl_adj_land = ', rhminl_adj_land_const
         116           0 :        write(iulog,*) '  rhminh          = ', rhminh_const
         117           0 :        write(iulog,*) '  premit          = ', premit
         118           0 :        write(iulog,*) '  premib          = ', premib
         119           0 :        write(iulog,*) '  i_rhminl        = ', i_rhminl
         120           0 :        write(iulog,*) '  i_rhmini        = ', i_rhmini
         121             :    end if
         122             :
         123             :
         124           0 :    call addfld ('RHMIN_LIQ',     (/ 'lev' /), 'A', 'fraction', 'Default critical RH for liquid-stratus')
         125           0 :    call addfld ('RHMIN_ICE',     (/ 'lev' /), 'A', 'fraction', 'Default critical RH for    ice-stratus')
         126           0 :    call addfld ('DRHMINPBL_LIQ', (/ 'lev' /), 'A', 'fraction', 'Drop of liquid-stratus critical RH by PBL turbulence')
         127           0 :    call addfld ('DRHMINPBL_ICE', (/ 'lev' /), 'A', 'fraction', 'Drop of    ice-stratus critical RH by PBL turbulence')
         128           0 :    call addfld ('DRHMINDET_LIQ', (/ 'lev' /), 'A', 'fraction', 'Drop of liquid-stratus critical RH by convective detrainment')
         129           0 :    call addfld ('DRHMINDET_ICE', (/ 'lev' /), 'A', 'fraction', 'Drop of    ice-stratus critical RH by convective detrainment')
         130             :
         131           0 :    end subroutine ini_macro
         132             :
         133             :    ! ------------------------------ !
         134             :    ! Stratiform Liquid Macrophysics !
         135             :    ! ------------------------------ !
         136             :
         137             :    ! In the version, 'macro --> micro --> advective forcing --> macro...'
         138             :    ! A_...: only 'advective forcing' without 'microphysical tendency'
         139             :    ! C_...: only 'microphysical tendency'
         140             :    ! D_...: only 'detrainment of cumulus condensate'
         141             :    ! So, 'A' and 'C' are exclusive.
         142             :
         143           0 :    subroutine mmacro_pcond( lchnk      , ncol       , dt         , p            , dp         ,              &
         144             :                             T0         , qv0        , ql0        , qi0          , nl0        , ni0        , &
         145             :                             A_T        , A_qv       , A_ql       , A_qi         , A_nl       , A_ni       , &
         146             :                             C_T        , C_qv       , C_ql       , C_qi         , C_nl       , C_ni       , C_qlst, &
         147             :                             D_T        , D_qv       , D_ql       , D_qi         , D_nl       , D_ni       , &
         148             :                             a_cud      , a_cu0      , clrw_old   , clri_old     , landfrac   , snowh      , &
         149             :                             tke        , qtl_flx    , qti_flx    , cmfr_det     , qlr_det    , qir_det    , &
         150             :                             s_tendout  , qv_tendout , ql_tendout , qi_tendout   , nl_tendout , ni_tendout , &
         151             :                             qme        , qvadj      , qladj      , qiadj        , qllim      , qilim      , &
         152             :                             cld        , al_st_star , ai_st_star , ql_st_star   , qi_st_star , do_cldice  )
         153             :
         154           0 :    use constituents,     only : qmin, cnst_get_ind
         155             :    use wv_saturation,    only : findsp_vc
         156             :    use cam_history,      only : outfld, hist_fld_active
         157             :
         158             :    integer   icol
         159             :    integer,  intent(in)    :: lchnk                        ! Chunk number
         160             :    integer,  intent(in)    :: ncol                         ! Number of active columns
         161             :
         162             :    ! Input-Output variables
         163             :
         164             :    real(r8), intent(inout) :: T0(pcols,pver)               ! Temperature [K]
         165             :    real(r8), intent(inout) :: qv0(pcols,pver)              ! Grid-mean water vapor specific humidity [kg/kg]
         166             :    real(r8), intent(inout) :: ql0(pcols,pver)              ! Grid-mean liquid water content [kg/kg]
         167             :    real(r8), intent(inout) :: qi0(pcols,pver)              ! Grid-mean ice water content [kg/kg]
         168             :    real(r8), intent(inout) :: nl0(pcols,pver)              ! Grid-mean number concentration of cloud liquid droplet [#/kg]
         169             :    real(r8), intent(inout) :: ni0(pcols,pver)              ! Grid-mean number concentration of cloud ice    droplet [#/kg]
         170             :
         171             :    ! Input variables
         172             :
         173             :    real(r8), intent(in)    :: dt                           ! Model integration time step [s]
         174             :    real(r8), intent(in)    :: p(pcols,pver)                ! Pressure at the layer mid-point [Pa]
         175             :    real(r8), intent(in)    :: dp(pcols,pver)               ! Pressure thickness [Pa] > 0
         176             :
         177             :    real(r8), intent(in)    :: A_T(pcols,pver)              ! Non-microphysical advective external forcing of T  [K/s]
         178             :    real(r8), intent(in)    :: A_qv(pcols,pver)             ! Non-microphysical advective external forcing of qv [kg/kg/s]
         179             :    real(r8), intent(in)    :: A_ql(pcols,pver)             ! Non-microphysical advective external forcing of ql [kg/kg/s]
         180             :    real(r8), intent(in)    :: A_qi(pcols,pver)             ! Non-microphysical advective external forcing of qi [kg/kg/s]
         181             :    real(r8), intent(in)    :: A_nl(pcols,pver)             ! Non-microphysical advective external forcing of nl [#/kg/s]
         182             :    real(r8), intent(in)    :: A_ni(pcols,pver)             ! Non-microphysical advective external forcing of ni [#/kg/s]
         183             :
         184             :    real(r8), intent(in)    :: C_T(pcols,pver)              ! Microphysical advective external forcing of T  [K/s]
         185             :    real(r8), intent(in)    :: C_qv(pcols,pver)             ! Microphysical advective external forcing of qv [kg/kg/s]
         186             :    real(r8), intent(in)    :: C_ql(pcols,pver)             ! Microphysical advective external forcing of ql [kg/kg/s]
         187             :    real(r8), intent(in)    :: C_qi(pcols,pver)             ! Microphysical advective external forcing of qi [kg/kg/s]
         188             :    real(r8), intent(in)    :: C_nl(pcols,pver)             ! Microphysical advective external forcing of nl [#/kg/s]
         189             :    real(r8), intent(in)    :: C_ni(pcols,pver)             ! Microphysical advective external forcing of ni [#/kg/s]
         190             :    real(r8), intent(in)    :: C_qlst(pcols,pver)           ! Microphysical advective external forcing of ql
         191             :                                                            ! within liquid stratus [kg/kg/s]
         192             :
         193             :    real(r8), intent(in)    :: D_T(pcols,pver)              ! Cumulus detrainment external forcing of T  [K/s]
         194             :    real(r8), intent(in)    :: D_qv(pcols,pver)             ! Cumulus detrainment external forcing of qv [kg/kg/s]
         195             :    real(r8), intent(in)    :: D_ql(pcols,pver)             ! Cumulus detrainment external forcing of ql [kg/kg/s]
         196             :    real(r8), intent(in)    :: D_qi(pcols,pver)             ! Cumulus detrainment external forcing of qi [kg/kg/s]
         197             :    real(r8), intent(in)    :: D_nl(pcols,pver)             ! Cumulus detrainment external forcing of nl [#/kg/s]
         198             :    real(r8), intent(in)    :: D_ni(pcols,pver)             ! Cumulus detrainment external forcing of qi [#/kg/s]
         199             :
         200             :    real(r8), intent(in)    :: a_cud(pcols,pver)            ! Old cumulus fraction before update
         201             :    real(r8), intent(in)    :: a_cu0(pcols,pver)            ! New cumulus fraction after update
         202             :
         203             :    real(r8), intent(in)    :: clrw_old(pcols,pver)         ! Clear sky fraction at the previous time step for liquid stratus process
         204             :    real(r8), intent(in)    :: clri_old(pcols,pver)         ! Clear sky fraction at the previous time step for    ice stratus process
         205             :    real(r8), pointer       :: tke(:,:)                     ! (pcols,pverp) TKE from the PBL scheme
         206             :    real(r8), pointer       :: qtl_flx(:,:)                 ! (pcols,pverp) overbar(w'qtl') from PBL scheme where qtl = qv + ql
         207             :    real(r8), pointer       :: qti_flx(:,:)                 ! (pcols,pverp) overbar(w'qti') from PBL scheme where qti = qv + qi
         208             :    real(r8), pointer       :: cmfr_det(:,:)                ! (pcols,pver)  Detrained mass flux from the convection scheme
         209             :    real(r8), pointer       :: qlr_det(:,:)                 ! (pcols,pver)  Detrained        ql from the convection scheme
         210             :    real(r8), pointer       :: qir_det(:,:)                 ! (pcols,pver)  Detrained        qi from the convection scheme
         211             :
         212             :    real(r8), intent(in)    :: landfrac(pcols)              ! Land fraction
         213             :    real(r8), intent(in)    :: snowh(pcols)                 ! Snow depth (liquid water equivalent)
         214             :    logical,  intent(in)    :: do_cldice                    ! Whether or not cldice should be prognosed
         215             :
         216             :    ! Output variables
         217             :
         218             :    real(r8), intent(out)   :: s_tendout(pcols,pver)        ! Net tendency of grid-mean s  from 'Micro+Macro' processes [J/kg/s]
         219             :    real(r8), intent(out)   :: qv_tendout(pcols,pver)       ! Net tendency of grid-mean qv from 'Micro+Macro' processes [kg/kg/s]
         220             :    real(r8), intent(out)   :: ql_tendout(pcols,pver)       ! Net tendency of grid-mean ql from 'Micro+Macro' processes [kg/kg/s]
         221             :    real(r8), intent(out)   :: qi_tendout(pcols,pver)       ! Net tendency of grid-mean qi from 'Micro+Macro' processes [kg/kg/s]
         222             :    real(r8), intent(out)   :: nl_tendout(pcols,pver)       ! Net tendency of grid-mean nl from 'Micro+Macro' processes [#/kg/s]
         223             :    real(r8), intent(out)   :: ni_tendout(pcols,pver)       ! Net tendency of grid-mean ni from 'Micro+Macro' processes [#/kg/s]
         224             :
         225             :    real(r8), intent(out)   :: qme  (pcols,pver)            ! Net condensation rate [kg/kg/s]
         226             :    real(r8), intent(out)   :: qvadj(pcols,pver)            ! adjustment tendency from "positive_moisture" call (vapor)
         227             :    real(r8), intent(out)   :: qladj(pcols,pver)            ! adjustment tendency from "positive_moisture" call (liquid)
         228             :    real(r8), intent(out)   :: qiadj(pcols,pver)            ! adjustment tendency from "positive_moisture" call (ice)
         229             :    real(r8), intent(out)   :: qllim(pcols,pver)            ! tendency from "instratus_condensate" call (liquid)
         230             :    real(r8), intent(out)   :: qilim(pcols,pver)            ! tendency from "instratus_condensate" call (ice)
         231             :
         232             :    real(r8), intent(out)   :: cld(pcols,pver)              ! Net cloud fraction ( 0 <= cld <= 1 )
         233             :    real(r8), intent(out)   :: al_st_star(pcols,pver)       ! Physical liquid stratus fraction
         234             :    real(r8), intent(out)   :: ai_st_star(pcols,pver)       ! Physical ice stratus fraction
         235             :    real(r8), intent(out)   :: ql_st_star(pcols,pver)       ! In-stratus LWC [kg/kg]
         236             :    real(r8), intent(out)   :: qi_st_star(pcols,pver)       ! In-stratus IWC [kg/kg]
         237             :
         238             :    ! --------------- !
         239             :    ! Local variables !
         240             :    ! --------------- !
         241             :    integer :: ixcldliq, ixcldice
         242             :
         243             :    integer :: i, j, k, iter, ii, jj                        ! Loop indexes
         244             :
         245             :    ! Thermodynamic state variables
         246             :
         247             :    real(r8) T(pcols,pver)                                  ! Temperature of equilibrium reference state
         248             :                                                            ! from which 'Micro & Macro' are computed [K]
         249             :    real(r8) T1(pcols,pver)                                 ! Temperature after 'fice_force' on T01
         250             :    real(r8) T_0(pcols,pver)                                ! Temperature after 'instratus_condensate' on T1
         251             :    real(r8) T_05(pcols,pver)                               ! Temperature after 'advection' on T_0
         252             :    real(r8) T_prime0(pcols,pver)                           ! Temperature after 'Macrophysics (QQ)' on T_05star
         253             :    real(r8) T_dprime(pcols,pver)                           ! Temperature after 'fice_force' on T_prime
         254             :    real(r8) T_star(pcols,pver)                             ! Temperature after 'instratus_condensate' on T_dprime
         255             :
         256             :    real(r8) qv(pcols,pver)                                 ! Grid-mean qv of equilibrium reference state from which
         257             :                                                            ! 'Micro & Macro' are computed [kg/kg]
         258             :    real(r8) qv1(pcols,pver)                                ! Grid-mean qv after 'fice_force' on qv01
         259             :    real(r8) qv_0(pcols,pver)                               ! Grid-mean qv after 'instratus_condensate' on qv1
         260             :    real(r8) qv_05(pcols,pver)                              ! Grid-mean qv after 'advection' on qv_0
         261             :    real(r8) qv_prime0(pcols,pver)                          ! Grid-mean qv after 'Macrophysics (QQ)' on qv_05star
         262             :    real(r8) qv_dprime(pcols,pver)                          ! Grid-mean qv after 'fice_force' on qv_prime
         263             :    real(r8) qv_star(pcols,pver)                            ! Grid-mean qv after 'instratus_condensate' on qv_dprime
         264             :
         265             :    real(r8) ql(pcols,pver)                                 ! Grid-mean ql of equilibrium reference state from which
         266             :                                                            ! 'Micro & Macro' are computed [kg/kg]
         267             :    real(r8) ql1(pcols,pver)                                ! Grid-mean ql after 'fice_force' on ql01
         268             :    real(r8) ql_0(pcols,pver)                               ! Grid-mean ql after 'instratus_condensate' on ql1
         269             :    real(r8) ql_05(pcols,pver)                              ! Grid-mean ql after 'advection' on ql_0
         270             :    real(r8) ql_prime0(pcols,pver)                          ! Grid-mean ql after 'Macrophysics (QQ)' on ql_05star
         271             :    real(r8) ql_dprime(pcols,pver)                          ! Grid-mean ql after 'fice_force' on ql_prime
         272             :    real(r8) ql_star(pcols,pver)                            ! Grid-mean ql after 'instratus_condensate' on ql_dprime
         273             :
         274             :    real(r8) qi(pcols,pver)                                 ! Grid-mean qi of equilibrium reference state from which
         275             :                                                            ! 'Micro & Macro' are computed [kg/kg]
         276             :    real(r8) qi1(pcols,pver)                                ! Grid-mean qi after 'fice_force' on qi01
         277             :    real(r8) qi_0(pcols,pver)                               ! Grid-mean qi after 'instratus_condensate' on qi1
         278             :    real(r8) qi_05(pcols,pver)                              ! Grid-mean qi after 'advection' on qi_0
         279             :    real(r8) qi_prime0(pcols,pver)                          ! Grid-mean qi after 'Macrophysics (QQ)' on qi_05star
         280             :    real(r8) qi_dprime(pcols,pver)                          ! Grid-mean qi after 'fice_force' on qi_prime
         281             :    real(r8) qi_star(pcols,pver)                            ! Grid-mean qi after 'instratus_condensate' on qi_dprime
         282             :
         283             :    real(r8) nl(pcols,pver)                                 ! Grid-mean nl of equilibrium reference state from which
         284             :                                                            ! 'Micro & Macro' are computed [kg/kg]
         285             :    real(r8) nl1(pcols,pver)                                ! Grid-mean nl after 'fice_force' on nl01
         286             :    real(r8) nl_0(pcols,pver)                               ! Grid-mean nl after 'instratus_condensate' on nl1
         287             :    real(r8) nl_05(pcols,pver)                              ! Grid-mean nl after 'advection' on nl_0
         288             :    real(r8) nl_prime0(pcols,pver)                          ! Grid-mean nl after 'Macrophysics (QQ)' on nl_05star
         289             :    real(r8) nl_dprime(pcols,pver)                          ! Grid-mean nl after 'fice_force' on nl_prime
         290             :    real(r8) nl_star(pcols,pver)                            ! Grid-mean nl after 'instratus_condensate' on nl_dprime
         291             :
         292             :    real(r8) ni(pcols,pver)                                 ! Grid-mean ni of equilibrium reference state from which
         293             :                                                            ! 'Micro & Macro' are computed [kg/kg]
         294             :    real(r8) ni1(pcols,pver)                                ! Grid-mean ni after 'fice_force' on ni01
         295             :    real(r8) ni_0(pcols,pver)                               ! Grid-mean ni after 'instratus_condensate' on ni1
         296             :    real(r8) ni_05(pcols,pver)                              ! Grid-mean ni after 'advection' on ni_0
         297             :    real(r8) ni_prime0(pcols,pver)                          ! Grid-mean ni after 'Macrophysics (QQ)' on ni_05star
         298             :    real(r8) ni_dprime(pcols,pver)                          ! Grid-mean ni after 'fice_force' on ni_prime
         299             :    real(r8) ni_star(pcols,pver)                            ! Grid-mean ni after 'instratus_condensate' on ni_dprime
         300             :
         301             :    real(r8) a_st(pcols,pver)                               ! Stratus fraction of equilibrium reference state
         302             :    real(r8) a_st_0(pcols,pver)                             ! Stratus fraction at '_0' state
         303             :    real(r8) a_st_star(pcols,pver)                          ! Stratus fraction at '_star' state
         304             :
         305             :    real(r8) al_st(pcols,pver)                              ! Liquid stratus fraction of equilibrium reference state
         306             :    real(r8) al_st_0(pcols,pver)                            ! Liquid stratus fraction at '_0' state
         307             :    real(r8) al_st_nc(pcols,pver)                           ! Non-physical liquid stratus fraction in the non-cumulus pixels
         308             :
         309             :    real(r8) ai_st(pcols,pver)                              ! Ice stratus fraction of equilibrium reference state
         310             :    real(r8) ai_st_0(pcols,pver)                            ! Ice stratus fraction at '_0' state
         311             :    real(r8) ai_st_nc(pcols,pver)                           ! Non-physical ice stratus fraction in the non-cumulus pixels
         312             :
         313             :    real(r8) ql_st(pcols,pver)                              ! In-stratus LWC of equilibrium reference state [kg/kg]
         314             :    real(r8) ql_st_0(pcols,pver)                            ! In-stratus LWC at '_0' state
         315             :
         316             :    real(r8) qi_st(pcols,pver)                              ! In-stratus IWC of equilibrium reference state [kg/kg]
         317             :    real(r8) qi_st_0(pcols,pver)                            ! In-stratus IWC at '_0' state
         318             :
         319             :  ! Cumulus properties
         320             :
         321             :    real(r8) dacudt(pcols,pver)
         322             :    real(r8) a_cu(pcols,pver)
         323             :
         324             :  ! Adjustment tendency in association with 'positive_moisture'
         325             :
         326             :    real(r8) Tten_pwi1(pcols,pver)                          ! Pre-process T  tendency of input equilibrium state [K/s]
         327             :    real(r8) qvten_pwi1(pcols,pver)                         ! Pre-process qv tendency of input equilibrium state [kg/kg/s]
         328             :    real(r8) qlten_pwi1(pcols,pver)                         ! Pre-process ql tendency of input equilibrium state [kg/kg/s]
         329             :    real(r8) qiten_pwi1(pcols,pver)                         ! Pre-process qi tendency of input equilibrium state [kg/kg/s]
         330             :    real(r8) nlten_pwi1(pcols,pver)                         ! Pre-process nl tendency of input equilibrium state [#/kg/s]
         331             :    real(r8) niten_pwi1(pcols,pver)                         ! Pre-process ni tendency of input equilibrium state [#/kg/s]
         332             :
         333             :    real(r8) Tten_pwi2(pcols,pver)                          ! Post-process T  tendency of provisional equilibrium state [K/s]
         334             :    real(r8) qvten_pwi2(pcols,pver)                         ! Post-process qv tendency of provisional equilibrium state [kg/kg/s]
         335             :    real(r8) qlten_pwi2(pcols,pver)                         ! Post-process ql tendency of provisional equilibrium state [kg/kg/s]
         336             :    real(r8) qiten_pwi2(pcols,pver)                         ! Post-process qi tendency of provisional equilibrium state [kg/kg/s]
         337             :    real(r8) nlten_pwi2(pcols,pver)                         ! Post-process nl tendency of provisoonal equilibrium state [#/kg/s]
         338             :    real(r8) niten_pwi2(pcols,pver)                         ! Post-process ni tendency of provisional equilibrium state [#/kg/s]
         339             :
         340             :    real(r8) A_T_adj(pcols,pver)                            ! After applying external advective forcing [K/s]
         341             :    real(r8) A_qv_adj(pcols,pver)                           ! After applying external advective forcing [kg/kg/s]
         342             :    real(r8) A_ql_adj(pcols,pver)                           ! After applying external advective forcing [kg/kg/s]
         343             :    real(r8) A_qi_adj(pcols,pver)                           ! After applying external advective forcing [kg/kg/s]
         344             :    real(r8) A_nl_adj(pcols,pver)                           ! After applying external advective forcing [#/kg/s]
         345             :    real(r8) A_ni_adj(pcols,pver)                           ! After applying external advective forcing [#/kg/s]
         346             :
         347             :  ! Adjustment tendency in association with 'instratus_condensate'
         348             :
         349             :    real(r8) QQw1(pcols,pver)           ! Effective adjustive condensation into water due to 'instratus_condensate' [kg/kg/s]
         350             :    real(r8) QQi1(pcols,pver)           ! Effective adjustive condensation into ice   due to 'instratus_condensate' [kg/kg/s]
         351             :    real(r8) QQw2(pcols,pver)           ! Effective adjustive condensation into water due to 'instratus_condensate' [kg/kg/s]
         352             :    real(r8) QQi2(pcols,pver)           ! Effective adjustive condensation into ice   due to 'instratus_condensate' [kg/kg/s]
         353             :
         354             :    real(r8) QQnl1(pcols,pver)          ! Tendency of nl associated with QQw1 only when QQw1<0 (net evaporation) [#/kg/s]
         355             :    real(r8) QQni1(pcols,pver)          ! Tendency of ni associated with QQi1 only when QQw1<0 (net evaporation) [#/kg/s]
         356             :    real(r8) QQnl2(pcols,pver)          ! Tendency of nl associated with QQw2 only when QQw2<0 (net evaporation) [#/kg/s]
         357             :    real(r8) QQni2(pcols,pver)          ! Tendency of ni associated with QQi2 only when QQw2<0 (net evaporation) [#/kg/s]
         358             :
         359             :  ! Macrophysical process tendency variables
         360             :
         361             :    real(r8) QQ(pcols,pver)             ! Net condensation rate into water+ice           [kg/kg/s]
         362             :    real(r8) QQw(pcols,pver)            ! Net condensation rate into water               [kg/kg/s]
         363             :    real(r8) QQi(pcols,pver)            ! Net condensation rate into ice                 [kg/kg/s]
         364             :    real(r8) QQnl(pcols,pver)           ! Tendency of nl associated with QQw both for condensation and evaporation [#/kg/s]
         365             :    real(r8) QQni(pcols,pver)           ! Tendency of ni associated with QQi both for condensation and evaporation [#/kg/s]
         366             :    real(r8) ACnl(pcols,pver)           ! Cloud liquid droplet (nl) activation tendency [#/kg/s]
         367             :    real(r8) ACni(pcols,pver)           ! Cloud ice    droplet (ni) activation tendency [#/kg/s]
         368             :
         369             :    real(r8) QQw_prev(pcols,pver)
         370             :    real(r8) QQi_prev(pcols,pver)
         371             :    real(r8) QQnl_prev(pcols,pver)
         372             :    real(r8) QQni_prev(pcols,pver)
         373             :
         374             :    real(r8) QQw_prog(pcols,pver)
         375             :    real(r8) QQi_prog(pcols,pver)
         376             :    real(r8) QQnl_prog(pcols,pver)
         377             :    real(r8) QQni_prog(pcols,pver)
         378             :
         379             :    real(r8) QQ_final(pcols,pver)
         380             :    real(r8) QQw_final(pcols,pver)
         381             :    real(r8) QQi_final(pcols,pver)
         382             :    real(r8) QQn_final(pcols,pver)
         383             :    real(r8) QQnl_final(pcols,pver)
         384             :    real(r8) QQni_final(pcols,pver)
         385             :
         386             :    real(r8) QQ_all(pcols,pver)         ! QQw_all    + QQi_all
         387             :    real(r8) QQw_all(pcols,pver)        ! QQw_final  + QQw1  + QQw2  + qlten_pwi1 + qlten_pwi2 + A_ql_adj [kg/kg/s]
         388             :    real(r8) QQi_all(pcols,pver)        ! QQi_final  + QQi1  + QQi2  + qiten_pwi1 + qiten_pwi2 + A_qi_adj [kg/kg/s]
         389             :    real(r8) QQn_all(pcols,pver)        ! QQnl_all   + QQni_all
         390             :    real(r8) QQnl_all(pcols,pver)       ! QQnl_final + QQnl1 + QQnl2 + nlten_pwi1 + nlten_pwi2 + ACnl [#/kg/s]
         391             :    real(r8) QQni_all(pcols,pver)       ! QQni_final + QQni1 + QQni2 + niten_pwi1 + niten_pwi2 + ACni [#/kg/s]
         392             :
         393             :  ! Coefficient for computing QQ and related processes
         394             :
         395             :    real(r8) U(pcols,pver)                                  ! Grid-mean RH
         396             :    real(r8) U_nc(pcols,pver)                               ! Mean RH of non-cumulus pixels
         397             :    real(r8) G_nc(pcols,pver)                               ! d(U_nc)/d(a_st_nc)
         398             :    real(r8) F_nc(pcols,pver)                               ! A function of second parameter for a_st_nc
         399             :    real(r8) alpha                                          ! = 1/qs
         400             :    real(r8) beta                                           ! = (qv/qs**2)*dqsdT
         401             :    real(r8) betast                                         ! = alpha*dqsdT
         402             :    real(r8) gammal                                         ! = alpha + (latvap/cpair)*beta
         403             :    real(r8) gammai                                         ! = alpha + ((latvap+latice)/cpair)*beta
         404             :    real(r8) gammaQ                                         ! = alpha + (latvap/cpair)*beta
         405             :    real(r8) deltal                                         ! = 1 + a_st*(latvap/cpair)*(betast/alpha)
         406             :    real(r8) deltai                                         ! = 1 + a_st*((latvap+latice)/cpair)*(betast/alpha)
         407             :    real(r8) A_Tc                                           ! Advective external forcing of Tc [K/s]
         408             :    real(r8) A_qt                                           ! Advective external forcing of qt [kg/kg/s]
         409             :    real(r8) C_Tc                                           ! Microphysical forcing of Tc [K/s]
         410             :    real(r8) C_qt                                           ! Microphysical forcing of qt [kg/kg/s]
         411             :    real(r8) dTcdt                                          ! d(Tc)/dt      [K/s]
         412             :    real(r8) dqtdt                                          ! d(qt)/dt      [kg/kg/s]
         413             :    real(r8) dqtstldt                                       ! d(qt_alst)/dt [kg/kg/s]
         414             :    real(r8) dqidt                                          ! d(qi)/dt      [kg/kg/s]
         415             :
         416             :    real(r8) dqlstdt                                        ! d(ql_st)/dt [kg/kg/s]
         417             :    real(r8) dalstdt                                        ! d(al_st)/dt  [1/s]
         418             :    real(r8) dastdt                                         ! d(a_st)/dt  [1/s]
         419             :
         420             :    real(r8) anic                                           ! Fractional area of non-cumulus and non-ice stratus fraction
         421             :    real(r8) GG                                             ! G_nc(i,k)/anic
         422             :
         423             :    real(r8) aa(2,2)
         424             :    real(r8) bb(2,1)
         425             :
         426             :    real(r8) zeros(pcols,pver)
         427             :
         428             :    real(r8) qmin1(pcols,pver)
         429             :    real(r8) qmin2(pcols,pver)
         430             :    real(r8) qmin3(pcols,pver)
         431             :
         432             :    real(r8) esat_a(pcols)                             ! Saturation water vapor pressure [Pa]
         433             :    real(r8) qsat_a(pcols,pver)                        ! Saturation water vapor specific humidity [kg/kg]
         434             :    real(r8) Twb_aw(pcols)                             ! Wet-bulb temperature [K]
         435             :    real(r8) qvwb_aw(pcols,pver)                       ! Wet-bulb water vapor specific humidity [kg/kg]
         436             :
         437             :    real(r8) esat_b(pcols)
         438             :    real(r8) qsat_b(pcols)
         439             :    real(r8) dqsdT_b(pcols)
         440             :
         441             :    logical  land
         442             :    real(r8) tmp
         443             :
         444             :    real(r8) d_rhmin_liq_PBL(pcols,pver)
         445             :    real(r8) d_rhmin_ice_PBL(pcols,pver)
         446             :    real(r8) d_rhmin_liq_det(pcols,pver)
         447             :    real(r8) d_rhmin_ice_det(pcols,pver)
         448             :    real(r8) rhmaxi_arr(pcols,pver)
         449             :    real(r8) rhmini_arr(pcols,pver)
         450             :    real(r8) rhminl_arr(pcols,pver)
         451             :    real(r8) rhminl_adj_land_arr(pcols,pver)
         452             :    real(r8) rhminh_arr(pcols,pver)
         453             :    real(r8) rhmin_liq_diag(pcols,pver)
         454             :    real(r8) rhmin_ice_diag(pcols,pver)
         455             :
         456             :    real(r8) QQmax,QQmin,QQwmin,QQimin                      ! For limiting QQ
         457             :    real(r8) cone                                           ! Number close to but smaller than 1
         458             :
         459           0 :    cone            = 0.999_r8
         460           0 :    zeros(:ncol,:)  = 0._r8
         461             :
         462             :    ! ------------------------------------ !
         463             :    ! Global initialization of main output !
         464             :    ! ------------------------------------ !
         465             :
         466           0 :      s_tendout(:ncol,:)     = 0._r8
         467           0 :      qv_tendout(:ncol,:)    = 0._r8
         468           0 :      ql_tendout(:ncol,:)    = 0._r8
         469           0 :      qi_tendout(:ncol,:)    = 0._r8
         470           0 :      nl_tendout(:ncol,:)    = 0._r8
         471           0 :      ni_tendout(:ncol,:)    = 0._r8
         472             :
         473           0 :      qme(:ncol,:)           = 0._r8
         474             :
         475           0 :      cld(:ncol,:)           = 0._r8
         476           0 :      al_st_star(:ncol,:)    = 0._r8
         477           0 :      ai_st_star(:ncol,:)    = 0._r8
         478           0 :      ql_st_star(:ncol,:)    = 0._r8
         479           0 :      qi_st_star(:ncol,:)    = 0._r8
         480             :
         481             :    ! --------------------------------------- !
         482             :    ! Initialization of internal 2D variables !
         483             :    ! --------------------------------------- !
         484             :
         485           0 :      T(:ncol,:)             = 0._r8
         486           0 :      T1(:ncol,:)            = 0._r8
         487           0 :      T_0(:ncol,:)           = 0._r8
         488           0 :      T_05(:ncol,:)          = 0._r8
         489           0 :      T_prime0(:ncol,:)      = 0._r8
         490           0 :      T_dprime(:ncol,:)      = 0._r8
         491           0 :      T_star(:ncol,:)        = 0._r8
         492             :
         493           0 :      qv(:ncol,:)            = 0._r8
         494           0 :      qv1(:ncol,:)           = 0._r8
         495           0 :      qv_0(:ncol,:)          = 0._r8
         496           0 :      qv_05(:ncol,:)         = 0._r8
         497           0 :      qv_prime0(:ncol,:)     = 0._r8
         498           0 :      qv_dprime(:ncol,:)     = 0._r8
         499           0 :      qv_star(:ncol,:)       = 0._r8
         500             :
         501           0 :      ql(:ncol,:)            = 0._r8
         502           0 :      ql1(:ncol,:)           = 0._r8
         503           0 :      ql_0(:ncol,:)          = 0._r8
         504           0 :      ql_05(:ncol,:)         = 0._r8
         505           0 :      ql_prime0(:ncol,:)     = 0._r8
         506           0 :      ql_dprime(:ncol,:)     = 0._r8
         507           0 :      ql_star(:ncol,:)       = 0._r8
         508             :
         509           0 :      qi(:ncol,:)            = 0._r8
         510           0 :      qi1(:ncol,:)           = 0._r8
         511           0 :      qi_0(:ncol,:)          = 0._r8
         512           0 :      qi_05(:ncol,:)         = 0._r8
         513           0 :      qi_prime0(:ncol,:)     = 0._r8
         514           0 :      qi_dprime(:ncol,:)     = 0._r8
         515           0 :      qi_star(:ncol,:)       = 0._r8
         516             :
         517           0 :      nl(:ncol,:)            = 0._r8
         518           0 :      nl1(:ncol,:)           = 0._r8
         519           0 :      nl_0(:ncol,:)          = 0._r8
         520           0 :      nl_05(:ncol,:)         = 0._r8
         521           0 :      nl_prime0(:ncol,:)     = 0._r8
         522           0 :      nl_dprime(:ncol,:)     = 0._r8
         523           0 :      nl_star(:ncol,:)       = 0._r8
         524             :
         525           0 :      ni(:ncol,:)            = 0._r8
         526           0 :      ni1(:ncol,:)           = 0._r8
         527           0 :      ni_0(:ncol,:)          = 0._r8
         528           0 :      ni_05(:ncol,:)         = 0._r8
         529           0 :      ni_prime0(:ncol,:)     = 0._r8
         530           0 :      ni_dprime(:ncol,:)     = 0._r8
         531           0 :      ni_star(:ncol,:)       = 0._r8
         532             :
         533           0 :      a_st(:ncol,:)          = 0._r8
         534           0 :      a_st_0(:ncol,:)        = 0._r8
         535           0 :      a_st_star(:ncol,:)     = 0._r8
         536             :
         537           0 :      al_st(:ncol,:)         = 0._r8
         538           0 :      al_st_0(:ncol,:)       = 0._r8
         539           0 :      al_st_nc(:ncol,:)      = 0._r8
         540             :
         541           0 :      ai_st(:ncol,:)         = 0._r8
         542           0 :      ai_st_0(:ncol,:)       = 0._r8
         543           0 :      ai_st_nc(:ncol,:)      = 0._r8
         544             :
         545           0 :      ql_st(:ncol,:)         = 0._r8
         546           0 :      ql_st_0(:ncol,:)       = 0._r8
         547             :
         548           0 :      qi_st(:ncol,:)         = 0._r8
         549           0 :      qi_st_0(:ncol,:)       = 0._r8
         550             :
         551             :  ! Cumulus properties
         552             :
         553           0 :      dacudt(:ncol,:)        = 0._r8
         554           0 :      a_cu(:ncol,:)          = 0._r8
         555             :
         556             :  ! Adjustment tendency in association with 'positive_moisture'
         557             :
         558           0 :      Tten_pwi1(:ncol,:)     = 0._r8
         559           0 :      qvten_pwi1(:ncol,:)    = 0._r8
         560           0 :      qlten_pwi1(:ncol,:)    = 0._r8
         561           0 :      qiten_pwi1(:ncol,:)    = 0._r8
         562           0 :      nlten_pwi1(:ncol,:)    = 0._r8
         563           0 :      niten_pwi1(:ncol,:)    = 0._r8
         564             :
         565           0 :      Tten_pwi2(:ncol,:)     = 0._r8
         566           0 :      qvten_pwi2(:ncol,:)    = 0._r8
         567           0 :      qlten_pwi2(:ncol,:)    = 0._r8
         568           0 :      qiten_pwi2(:ncol,:)    = 0._r8
         569           0 :      nlten_pwi2(:ncol,:)    = 0._r8
         570           0 :      niten_pwi2(:ncol,:)    = 0._r8
         571             :
         572           0 :      A_T_adj(:ncol,:)       = 0._r8
         573           0 :      A_qv_adj(:ncol,:)      = 0._r8
         574           0 :      A_ql_adj(:ncol,:)      = 0._r8
         575           0 :      A_qi_adj(:ncol,:)      = 0._r8
         576           0 :      A_nl_adj(:ncol,:)      = 0._r8
         577           0 :      A_ni_adj(:ncol,:)      = 0._r8
         578             :
         579           0 :      qvadj   (:ncol,:)      = 0._r8
         580           0 :      qladj   (:ncol,:)      = 0._r8
         581           0 :      qiadj   (:ncol,:)      = 0._r8
         582             :
         583             :  ! Adjustment tendency in association with 'instratus_condensate'
         584             :
         585           0 :      QQw1(:ncol,:)          = 0._r8
         586           0 :      QQi1(:ncol,:)          = 0._r8
         587           0 :      QQw2(:ncol,:)          = 0._r8
         588           0 :      QQi2(:ncol,:)          = 0._r8
         589             :
         590           0 :      QQnl1(:ncol,:)         = 0._r8
         591           0 :      QQni1(:ncol,:)         = 0._r8
         592           0 :      QQnl2(:ncol,:)         = 0._r8
         593           0 :      QQni2(:ncol,:)         = 0._r8
         594             :
         595           0 :      QQnl(:ncol,:)          = 0._r8
         596           0 :      QQni(:ncol,:)          = 0._r8
         597             :
         598             :  ! Macrophysical process tendency variables
         599             :
         600           0 :      QQ(:ncol,:)            = 0._r8
         601           0 :      QQw(:ncol,:)           = 0._r8
         602           0 :      QQi(:ncol,:)           = 0._r8
         603           0 :      QQnl(:ncol,:)          = 0._r8
         604           0 :      QQni(:ncol,:)          = 0._r8
         605           0 :      ACnl(:ncol,:)          = 0._r8
         606           0 :      ACni(:ncol,:)          = 0._r8
         607             :
         608           0 :      QQw_prev(:ncol,:)      = 0._r8
         609           0 :      QQi_prev(:ncol,:)      = 0._r8
         610           0 :      QQnl_prev(:ncol,:)     = 0._r8
         611           0 :      QQni_prev(:ncol,:)     = 0._r8
         612             :
         613           0 :      QQw_prog(:ncol,:)      = 0._r8
         614           0 :      QQi_prog(:ncol,:)      = 0._r8
         615           0 :      QQnl_prog(:ncol,:)     = 0._r8
         616           0 :      QQni_prog(:ncol,:)     = 0._r8
         617             :
         618           0 :      QQ_final(:ncol,:)      = 0._r8
         619           0 :      QQw_final(:ncol,:)     = 0._r8
         620           0 :      QQi_final(:ncol,:)     = 0._r8
         621           0 :      QQn_final(:ncol,:)     = 0._r8
         622           0 :      QQnl_final(:ncol,:)    = 0._r8
         623           0 :      QQni_final(:ncol,:)    = 0._r8
         624             :
         625           0 :      QQ_all(:ncol,:)        = 0._r8
         626           0 :      QQw_all(:ncol,:)       = 0._r8
         627           0 :      QQi_all(:ncol,:)       = 0._r8
         628           0 :      QQn_all(:ncol,:)       = 0._r8
         629           0 :      QQnl_all(:ncol,:)      = 0._r8
         630           0 :      QQni_all(:ncol,:)      = 0._r8
         631             :
         632             :  ! Coefficient for computing QQ and related processes
         633             :
         634           0 :      U(:ncol,:)             = 0._r8
         635           0 :      U_nc(:ncol,:)          = 0._r8
         636           0 :      G_nc(:ncol,:)          = 0._r8
         637           0 :      F_nc(:ncol,:)          = 0._r8
         638             :
         639             :  ! Other
         640             :
         641           0 :      qmin1(:ncol,:)         = 0._r8
         642           0 :      qmin2(:ncol,:)         = 0._r8
         643           0 :      qmin3(:ncol,:)         = 0._r8
         644             :
         645             :    ! ---------------- !
         646             :    ! Main computation !
         647             :    ! ---------------- !
         648             :
         649             :    ! Compute critical RH for stratus
         650           0 :    rhmaxi_arr(:ncol,:pver) = rhmaxi_const
         651             :    call rhcrit_calc( &
         652             :       ncol, dp, T0, p, &
         653             :       clrw_old, clri_old, tke, qtl_flx, &
         654             :       qti_flx, cmfr_det, qlr_det, qir_det, &
         655             :       rhmaxi_arr, rhmini_arr, rhminl_arr, rhminl_adj_land_arr, rhminh_arr, &
         656           0 :       d_rhmin_liq_PBL, d_rhmin_ice_PBL, d_rhmin_liq_det, d_rhmin_ice_det)
         657             :
         658             :    ! ---------------------------------- !
         659             :    ! Compute cumulus-related properties !
         660             :    ! ---------------------------------- !
         661             :
         662             :    dacudt(:ncol,top_lev:pver) = &
         663           0 :         (a_cu0(:ncol,top_lev:pver) - a_cud(:ncol,top_lev:pver))/dt
         664             :
         665             :    ! ---------------------------------------------------------------------- !
         666             :    ! set to zero for levels above
         667             :    ! ---------------------------------------------------------------------- !
         668           0 :    ql0(:ncol,:top_lev-1) = 0._r8
         669           0 :    qi0(:ncol,:top_lev-1) = 0._r8
         670           0 :    nl0(:ncol,:top_lev-1) = 0._r8
         671           0 :    ni0(:ncol,:top_lev-1) = 0._r8
         672             :
         673             :    ! ---------------------------------------------------------------------- !
         674             :    ! Check if input non-cumulus pixels satisfie a non-negative constraint.  !
         675             :    ! If not, force all water vapor substances to be positive in all layers. !
         676             :    ! We should use 'old' cumulus properties for this routine.               !
         677             :    ! ---------------------------------------------------------------------- !
         678             :
         679           0 :    T1(:ncol,:)    =  T0(:ncol,:)
         680           0 :    qv1(:ncol,:)   = qv0(:ncol,:)
         681           0 :    ql1(:ncol,:)   = ql0(:ncol,:)
         682           0 :    qi1(:ncol,:)   = qi0(:ncol,:)
         683           0 :    nl1(:ncol,:)   = nl0(:ncol,:)
         684           0 :    ni1(:ncol,:)   = ni0(:ncol,:)
         685             :
         686             :
         687           0 :    call cnst_get_ind( 'CLDLIQ', ixcldliq )
         688           0 :    call cnst_get_ind( 'CLDICE', ixcldice )
         689             :
         690             :
         691           0 :    qmin1(:ncol,:) = qmin(1)
         692           0 :    qmin2(:ncol,:) = qmin(ixcldliq)
         693           0 :    qmin3(:ncol,:) = qmin(ixcldice)
         694             :
         695             :    call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp, &
         696             :                            qv1, ql1, qi1, T1, qvten_pwi1, qlten_pwi1, &
         697           0 :                            qiten_pwi1, Tten_pwi1, do_cldice)
         698             :
         699           0 :    do k = top_lev, pver
         700           0 :    do i = 1, ncol
         701           0 :       if( ql1(i,k) .lt. qsmall ) then
         702           0 :           nlten_pwi1(i,k) = -nl1(i,k)/dt
         703           0 :           nl1(i,k)        = 0._r8
         704             :       endif
         705           0 :       if( qi1(i,k) .lt. qsmall ) then
         706           0 :           niten_pwi1(i,k) = -ni1(i,k)/dt
         707           0 :           ni1(i,k)        = 0._r8
         708             :       endif
         709             :    enddo
         710             :    enddo
         711             :
         712             :    ! ------------------------------------------------------------- !
         713             :    ! Impose 'in-stratus condensate amount constraint'              !
         714             :    ! such that it is bounded by two limiting values.               !
         715             :    ! This should also use 'old' cumulus properties since it is     !
         716             :    ! before applying external forcings.                            !
         717             :    ! Below 'QQw1,QQi1' are effective adjustive condensation        !
         718             :    ! Although this process also involves freezing of cloud         !
         719             :    ! liquid into ice, they can be and only can be expressed        !
         720             :    ! in terms of effective condensation.                           !
         721             :    ! ------------------------------------------------------------- !
         722             :
         723           0 :    do k = top_lev, pver
         724             :       call instratus_condensate( lchnk, ncol, k,                                   &
         725             :                                  p(:,k), T1(:,k), qv1(:,k), ql1(:,k), qi1(:,k),    &
         726             :                                  ni1(:,k),                                         &
         727             :                                  a_cud(:,k), zeros(:,k), zeros(:,k),               &
         728             :                                  zeros(:,k), zeros(:,k), zeros(:,k),               &
         729             :                                  landfrac, snowh,                                  &
         730             :                                  rhmaxi_arr(:,k),rhmini_arr(:,k), rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k), &
         731             :                                  T_0(:,k), qv_0(:,k), ql_0(:,k), qi_0(:,k),        &
         732           0 :                                  al_st_0(:,k), ai_st_0(:,k), ql_st_0(:,k), qi_st_0(:,k) )
         733           0 :       a_st_0(:ncol,k) = max(al_st_0(:ncol,k),ai_st_0(:ncol,k))
         734           0 :       QQw1(:ncol,k)   = (ql_0(:ncol,k) - ql1(:ncol,k))/dt
         735           0 :       QQi1(:ncol,k)   = (qi_0(:ncol,k) - qi1(:ncol,k))/dt
         736             :       ! -------------------------------------------------- !
         737             :       ! Reduce droplet concentration if evaporation occurs !
         738             :       ! Set a limit such that negative state not happens.  !
         739             :       ! -------------------------------------------------- !
         740           0 :       do i = 1, ncol
         741           0 :          if( QQw1(i,k) .le. 0._r8 ) then
         742           0 :              if( ql1(i,k) .gt. qsmall ) then
         743           0 :                  QQnl1(i,k) = QQw1(i,k)*nl1(i,k)/ql1(i,k)
         744           0 :                  QQnl1(i,k) = min(0._r8,cone*max(QQnl1(i,k),-nl1(i,k)/dt))
         745             :              else
         746           0 :                  QQnl1(i,k) = 0._r8
         747             :              endif
         748             :          endif
         749           0 :          if( QQi1(i,k) .le. 0._r8 ) then
         750           0 :              if( qi1(i,k) .gt. qsmall ) then
         751           0 :                  QQni1(i,k) = QQi1(i,k)*ni1(i,k)/qi1(i,k)
         752           0 :                  QQni1(i,k) = min(0._r8,cone*max(QQni1(i,k),-ni1(i,k)/dt))
         753             :              else
         754           0 :                  QQni1(i,k) = 0._r8
         755             :              endif
         756             :          endif
         757             :       enddo
         758             :    enddo
         759           0 :    nl_0(:ncol,top_lev:) = max(0._r8,nl1(:ncol,top_lev:)+QQnl1(:ncol,top_lev:)*dt)
         760           0 :    ni_0(:ncol,top_lev:) = max(0._r8,ni1(:ncol,top_lev:)+QQni1(:ncol,top_lev:)*dt)
         761             :
         762             :    ! ----------------------------------------------------------------------------- !
         763             :    ! Check if non-cumulus pixels of '_05' state satisfies non-negative constraint. !
         764             :    ! If not, force all water substances of '_05' state to be positive by imposing  !
         765             :    ! adjustive advection. We should use 'new' cumulus properties for this routine. !
         766             :    ! ----------------------------------------------------------------------------- !
         767             :
         768           0 :    T_05(:ncol,top_lev:)  =  T_0(:ncol,top_lev:) + (  A_T(:ncol,top_lev:) +  C_T(:ncol,top_lev:) ) * dt
         769           0 :    qv_05(:ncol,top_lev:) = qv_0(:ncol,top_lev:) + ( A_qv(:ncol,top_lev:) + C_qv(:ncol,top_lev:) ) * dt
         770           0 :    ql_05(:ncol,top_lev:) = ql_0(:ncol,top_lev:) + ( A_ql(:ncol,top_lev:) + C_ql(:ncol,top_lev:) ) * dt
         771           0 :    qi_05(:ncol,top_lev:) = qi_0(:ncol,top_lev:) + ( A_qi(:ncol,top_lev:) + C_qi(:ncol,top_lev:) ) * dt
         772           0 :    nl_05(:ncol,top_lev:) = max(0._r8, nl_0(:ncol,top_lev:) + ( A_nl(:ncol,top_lev:) + C_nl(:ncol,top_lev:) ) * dt )
         773           0 :    ni_05(:ncol,top_lev:) = max(0._r8, ni_0(:ncol,top_lev:) + ( A_ni(:ncol,top_lev:) + C_ni(:ncol,top_lev:) ) * dt )
         774             :
         775             :    call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp, &
         776             :                            qv_05, ql_05, qi_05, T_05, A_qv_adj, &
         777           0 :                            A_ql_adj, A_qi_adj, A_T_adj, do_cldice)
         778             :
         779             :    ! -------------------------------------------------------------- !
         780             :    ! Define reference state at the first iteration. This will be    !
         781             :    ! continuously updated within the iteration loop below.          !
         782             :    ! While equlibrium state properties are already output from the  !
         783             :    ! 'instratus_condensate', they will be re-computed within the    !
         784             :    ! each iteration process. At the first iteration, they will      !
         785             :    ! produce exactly identical results. Note that except at the     !
         786             :    ! very first iteration iter = 1, we must use updated cumulus     !
         787             :    ! properties at all the other iteration processes. Even at the   !
         788             :    ! first iteration, we should use updated cumulus properties      !
         789             :    ! when computing limiters for (Q,P,E).                           !
         790             :    ! -------------------------------------------------------------- !
         791             :
         792             :    ! -------------------------------------------------------------- !
         793             :    ! Define variables at the reference state of the first iteration !
         794             :    ! -------------------------------------------------------------- !
         795             :
         796           0 :    T(:ncol,top_lev:)     = T_0(:ncol,top_lev:)
         797           0 :    qv(:ncol,top_lev:)    = qv_0(:ncol,top_lev:)
         798           0 :    ql(:ncol,top_lev:)    = ql_0(:ncol,top_lev:)
         799           0 :    qi(:ncol,top_lev:)    = qi_0(:ncol,top_lev:)
         800           0 :    al_st(:ncol,top_lev:) = al_st_0(:ncol,top_lev:)
         801           0 :    ai_st(:ncol,top_lev:) = ai_st_0(:ncol,top_lev:)
         802           0 :    a_st(:ncol,top_lev:)  = a_st_0(:ncol,top_lev:)
         803           0 :    ql_st(:ncol,top_lev:) = ql_st_0(:ncol,top_lev:)
         804           0 :    qi_st(:ncol,top_lev:) = qi_st_0(:ncol,top_lev:)
         805           0 :    nl(:ncol,top_lev:)    = nl_0(:ncol,top_lev:)
         806           0 :    ni(:ncol,top_lev:)    = ni_0(:ncol,top_lev:)
         807             :
         808             :    ! -------------------------- !
         809             :    ! Main iterative computation !
         810             :    ! -------------------------- !
         811             :
         812           0 :    do k = top_lev, pver
         813             :       call findsp_vc(qv_05(:ncol,k), T_05(:ncol,k), p(:ncol,k), .false., &
         814           0 :                      Twb_aw(:ncol), qvwb_aw(:ncol,k))
         815           0 :       call qsat_water(T_05(1:ncol,k), p(1:ncol,k), esat_a(1:ncol), qsat_a(1:ncol,k), ncol)
         816             :    enddo
         817             :
         818           0 :    do iter = 1, niter
         819             :
         820             :       ! ------------------------------------------ !
         821             :       ! Initialize array within the iteration loop !
         822             :       ! ------------------------------------------ !
         823             :
         824           0 :       QQ(:,:)         = 0._r8
         825           0 :       QQw(:,:)        = 0._r8
         826           0 :       QQi(:,:)        = 0._r8
         827           0 :       QQnl(:,:)       = 0._r8
         828           0 :       QQni(:,:)       = 0._r8
         829           0 :       QQw2(:,:)       = 0._r8
         830           0 :       QQi2(:,:)       = 0._r8
         831           0 :       QQnl2(:,:)      = 0._r8
         832           0 :       QQni2(:,:)      = 0._r8
         833           0 :       nlten_pwi2(:,:) = 0._r8
         834           0 :       niten_pwi2(:,:) = 0._r8
         835           0 :       ACnl(:,:)       = 0._r8
         836           0 :       ACni(:,:)       = 0._r8
         837           0 :       aa(:,:)         = 0._r8
         838           0 :       bb(:,:)         = 0._r8
         839             :
         840           0 :       do k = top_lev, pver
         841           0 :          call qsat_water(T(1:ncol,k), p(1:ncol,k), esat_b(1:ncol), qsat_b(1:ncol), ncol, dqsdt=dqsdT_b(1:ncol))
         842           0 :       if( iter .eq. 1 ) then
         843           0 :           a_cu(:ncol,k) = a_cud(:ncol,k)
         844             :       else
         845           0 :           a_cu(:ncol,k) = a_cu0(:ncol,k)
         846             :       endif
         847           0 :       do i = 1, ncol
         848           0 :          U(i,k)    =  qv(i,k)/qsat_b(i)
         849           0 :          U_nc(i,k) =  U(i,k)
         850             :       enddo
         851             :       if( CAMstfrac ) then
         852             :           call astG_RHU(U_nc(:,k),p(:,k),qv(:,k),landfrac(:),snowh(:),al_st_nc(:,k),G_nc(:,k),ncol,&
         853             :                         rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k))
         854             :       else
         855             :           call astG_PDF(U_nc(:,k),p(:,k),qv(:,k),landfrac(:),snowh(:),al_st_nc(:,k),G_nc(:,k),ncol,&
         856           0 :                         rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k))
         857             :       endif
         858             :       call aist_vector(qv(:,k),T(:,k),p(:,k),qi(:,k),ni(:,k),landfrac(:),snowh(:),ai_st_nc(:,k),ncol,&
         859           0 :                        rhmaxi_arr(:,k), rhmini_arr(:,k), rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k))
         860             :
         861           0 :       ai_st(:ncol,k)  =  (1._r8-a_cu(:ncol,k))*ai_st_nc(:ncol,k)
         862           0 :       al_st(:ncol,k)  =  (1._r8-a_cu(:ncol,k))*al_st_nc(:ncol,k)
         863           0 :       a_st(:ncol,k)   =  max(al_st(:ncol,k),ai_st(:ncol,k))
         864             :
         865           0 :       do i = 1, ncol
         866             :
         867             :          ! -------------------------------------------------------- !
         868             :          ! Compute basic thermodynamic coefficients for computing Q !
         869             :          ! -------------------------------------------------------- !
         870             :
         871           0 :          alpha  =  1._r8/qsat_b(i)
         872           0 :          beta   =  dqsdT_b(i)*(qv(i,k)/qsat_b(i)**2)
         873           0 :          betast =  alpha*dqsdT_b(i)
         874           0 :          gammal =  alpha + (latvap/cpair)*beta
         875           0 :          gammai =  alpha + ((latvap+latice)/cpair)*beta
         876           0 :          gammaQ =  alpha + (latvap/cpair)*beta
         877           0 :          deltal =  1._r8 + a_st(i,k)*(latvap/cpair)*(betast/alpha)
         878           0 :          deltai =  1._r8 + a_st(i,k)*((latvap+latice)/cpair)*(betast/alpha)
         879           0 :          A_Tc   =  A_T(i,k)+A_T_adj(i,k)-(latvap/cpair)*(A_ql(i,k)+A_ql_adj(i,k))-((latvap+latice)/cpair)*(A_qi(i,k)+A_qi_adj(i,k))
         880           0 :          A_qt   =  A_qv(i,k) + A_qv_adj(i,k) + A_ql(i,k) + A_ql_adj(i,k) + A_qi(i,k) + A_qi_adj(i,k)
         881           0 :          C_Tc   =  C_T(i,k) - (latvap/cpair)*C_ql(i,k) - ((latvap+latice)/cpair)*C_qi(i,k)
         882           0 :          C_qt   =  C_qv(i,k) + C_ql(i,k) + C_qi(i,k)
         883           0 :          dTcdt  =  A_Tc + C_Tc
         884           0 :          dqtdt  =  A_qt + C_qt
         885             :        ! dqtstldt = A_qt + C_ql(i,k)/max(1.e-2_r8,al_st(i,k))                             ! Original
         886             :        ! dqtstldt = A_qt - A_qi(i,k) - A_qi_adj(i,k) + C_ql(i,k)/max(1.e-2_r8,al_st(i,k)) ! New 1 on Dec.30.2009.
         887           0 :          dqtstldt = A_qt - A_qi(i,k) - A_qi_adj(i,k) + C_qlst(i,k)                        ! New 2 on Dec.30.2009.
         888             :        ! dqtstldt = A_qt + C_qt                                                           ! Original Conservative treatment
         889             :        ! dqtstldt = A_qt - A_qi(i,k) - A_qi_adj(i,k) + C_qt - C_qi(i,k)            ! New Conservative treatment on Dec.30.2009
         890           0 :          dqidt = A_qi(i,k) + A_qi_adj(i,k) + C_qi(i,k)
         891             :
         892           0 :          anic    = max(1.e-8_r8,(1._r8-a_cu(i,k)))
         893           0 :          GG      = G_nc(i,k)/anic
         894           0 :          aa(1,1) = gammal*al_st(i,k)
         895           0 :          aa(1,2) = GG + gammal*cc*ql_st(i,k)
         896           0 :          aa(2,1) = alpha + (latvap/cpair)*betast*al_st(i,k)
         897           0 :          aa(2,2) = (latvap/cpair)*betast*cc*ql_st(i,k)
         898           0 :          bb(1,1) = alpha*dqtdt - beta*dTcdt - gammai*dqidt - GG*al_st_nc(i,k)*dacudt(i,k) + F_nc(i,k)
         899           0 :          bb(2,1) = alpha*dqtstldt - betast*(dTcdt + ((latvap+latice)/cpair)*dqidt)
         900           0 :          call gaussj(aa(1:2,1:2),2,2,bb(1:2,1),1,1)
         901           0 :          dqlstdt = bb(1,1)
         902           0 :          dalstdt = bb(2,1)
         903           0 :          QQ(i,k) = al_st(i,k)*dqlstdt + cc*ql_st(i,k)*dalstdt - ( A_ql(i,k) + A_ql_adj(i,k) + C_ql(i,k) )
         904             :
         905             :        ! ------------------------------------------------------------ !
         906             :        ! Limiter for QQ                                               !
         907             :        ! Here, 'fice' should be from the reference equilibrium state  !
         908             :        ! since QQ itself is computed from the reference state.        !
         909             :        ! From the assumption used for derivation of QQ(i), it must be !
         910             :        ! that QQw(i) = QQ(i)*(1._r8-fice(i)), QQi(i) = QQ(i)*fice(i)  !
         911             :        ! ------------------------------------------------------------ !
         912             :
         913           0 :          if( QQ(i,k) .ge. 0._r8 ) then
         914           0 :              QQmax    = (qv_05(i,k) - qmin(1))/dt ! For ghost cumulus & semi-ghost ice stratus
         915           0 :              QQmax    = max(0._r8,QQmax)
         916           0 :              QQ(i,k)  = min(QQ(i,k),QQmax)
         917           0 :              QQw(i,k) = QQ(i,k)
         918           0 :              QQi(i,k) = 0._r8
         919             :          else
         920           0 :              QQmin  = 0._r8
         921           0 :              if( qv_05(i,k) .lt. qsat_a(i,k) ) QQmin = min(0._r8,cone*(qv_05(i,k)-qvwb_aw(i,k))/dt)
         922           0 :              QQ(i,k)  = max(QQ(i,k),QQmin)
         923           0 :              QQw(i,k) = QQ(i,k)
         924           0 :              QQi(i,k) = 0._r8
         925           0 :              QQwmin   = min(0._r8,-cone*ql_05(i,k)/dt)
         926           0 :              QQimin   = min(0._r8,-cone*qi_05(i,k)/dt)
         927           0 :              QQw(i,k) = min(0._r8,max(QQw(i,k),QQwmin))
         928             :              QQi(i,k) = min(0._r8,max(QQi(i,k),QQimin))
         929             :          endif
         930             :
         931             :        ! -------------------------------------------------- !
         932             :        ! Reduce droplet concentration if evaporation occurs !
         933             :        ! Note 'QQnl1,QQni1' are computed from the reference !
         934             :        ! equilibrium state but limiter is from 'nl_05'.     !
         935             :        ! -------------------------------------------------- !
         936             :
         937           0 :          if( QQw(i,k) .lt. 0._r8 ) then
         938           0 :              if( ql_05(i,k) .gt. qsmall ) then
         939           0 :                  QQnl(i,k) = QQw(i,k)*nl_05(i,k)/ql_05(i,k)
         940           0 :                  QQnl(i,k) = min(0._r8,cone*max(QQnl(i,k),-nl_05(i,k)/dt))
         941             :              else
         942           0 :                  QQnl(i,k) = 0._r8
         943             :              endif
         944             :          endif
         945             :
         946           0 :          if( QQi(i,k) .lt. 0._r8 ) then
         947           0 :              if( qi_05(i,k) .gt. qsmall ) then
         948           0 :                  QQni(i,k) = QQi(i,k)*ni_05(i,k)/qi_05(i,k)
         949           0 :                  QQni(i,k) = min(0._r8,cone*max(QQni(i,k),-ni_05(i,k)/dt))
         950             :              else
         951           0 :                  QQni(i,k) = 0._r8
         952             :              endif
         953             :          endif
         954             :
         955             :       enddo
         956             :       enddo
         957             :
         958             :     ! -------------------------------------------------------------------- !
         959             :     ! Until now, we have finished computing all necessary tendencies       !
         960             :     ! from the equilibrium input state (T_0).                              !
         961             :     ! If ramda = 0 : fully explicit scheme                                 !
         962             :     !    ramda = 1 : fully implicit scheme                                 !
         963             :     ! Note that 'ramda = 0.5 with niter = 2' can mimic                     !
         964             :     ! -------------------------------------------------------------------- !
         965             :
         966           0 :       if( iter .eq. 1 ) then
         967           0 :           QQw_prev(:ncol,top_lev:)  = QQw(:ncol,top_lev:)
         968           0 :           QQi_prev(:ncol,top_lev:)  = QQi(:ncol,top_lev:)
         969           0 :           QQnl_prev(:ncol,top_lev:) = QQnl(:ncol,top_lev:)
         970           0 :           QQni_prev(:ncol,top_lev:) = QQni(:ncol,top_lev:)
         971             :       endif
         972             :
         973           0 :       QQw_prog(:ncol,top_lev:)   = ramda*QQw(:ncol,top_lev:)   + (1._r8-ramda)*QQw_prev(:ncol,top_lev:)
         974           0 :       QQi_prog(:ncol,top_lev:)   = ramda*QQi(:ncol,top_lev:)   + (1._r8-ramda)*QQi_prev(:ncol,top_lev:)
         975           0 :       QQnl_prog(:ncol,top_lev:)  = ramda*QQnl(:ncol,top_lev:)  + (1._r8-ramda)*QQnl_prev(:ncol,top_lev:)
         976           0 :       QQni_prog(:ncol,top_lev:)  = ramda*QQni(:ncol,top_lev:)  + (1._r8-ramda)*QQni_prev(:ncol,top_lev:)
         977             :
         978           0 :       QQw_prev(:ncol,top_lev:)   = QQw_prog(:ncol,top_lev:)
         979           0 :       QQi_prev(:ncol,top_lev:)   = QQi_prog(:ncol,top_lev:)
         980           0 :       QQnl_prev(:ncol,top_lev:)  = QQnl_prog(:ncol,top_lev:)
         981           0 :       QQni_prev(:ncol,top_lev:)  = QQni_prog(:ncol,top_lev:)
         982             :
         983             :     ! -------------------------------------------------------- !
         984             :     ! Compute final prognostic state on which final diagnostic !
         985             :     ! in-stratus condensate adjustment is applied in the below.!
         986             :     ! Important : I must check whether there are any external  !
         987             :     !             advective forcings of 'A_nl(i,k),A_ni(i,k)'. !
         988             :     !             Even they are (i.e., advection of aerosol),  !
         989             :     !             actual droplet activation will be performd   !
         990             :     !             in microphysics, so it will be completely    !
         991             :     !             reasonable to 'A_nl(i,k)=A_ni(i,k)=0'.       !
         992             :     ! -------------------------------------------------------- !
         993             :
         994           0 :     do k = top_lev, pver
         995           0 :     do i = 1, ncol
         996           0 :        T_prime0(i,k)  = T_0(i,k)  + dt*( A_T(i,k)  +  A_T_adj(i,k) +  C_T(i,k) + &
         997           0 :             (latvap*QQw_prog(i,k)+(latvap+latice)*QQi_prog(i,k))/cpair )
         998           0 :        qv_prime0(i,k) = qv_0(i,k) + dt*( A_qv(i,k) + A_qv_adj(i,k) + C_qv(i,k) - QQw_prog(i,k) - QQi_prog(i,k) )
         999           0 :        ql_prime0(i,k) = ql_0(i,k) + dt*( A_ql(i,k) + A_ql_adj(i,k) + C_ql(i,k) + QQw_prog(i,k) )
        1000           0 :        qi_prime0(i,k) = qi_0(i,k) + dt*( A_qi(i,k) + A_qi_adj(i,k) + C_qi(i,k) + QQi_prog(i,k) )
        1001           0 :        nl_prime0(i,k) = max(0._r8,nl_0(i,k) + dt*( A_nl(i,k) + C_nl(i,k) + QQnl_prog(i,k) ))
        1002           0 :        ni_prime0(i,k) = max(0._r8,ni_0(i,k) + dt*( A_ni(i,k) + C_ni(i,k) + QQni_prog(i,k) ))
        1003           0 :        if( ql_prime0(i,k) .lt. qsmall ) nl_prime0(i,k) = 0._r8
        1004           0 :        if( qi_prime0(i,k) .lt. qsmall ) ni_prime0(i,k) = 0._r8
        1005             :     enddo
        1006             :     enddo
        1007             :
        1008             :    ! -------------------------------------------------- !
        1009             :    ! Perform diagnostic 'positive_moisture' constraint. !
        1010             :    ! -------------------------------------------------- !
        1011             :
        1012           0 :    T_dprime(:ncol,top_lev:)  =  T_prime0(:ncol,top_lev:)
        1013           0 :    qv_dprime(:ncol,top_lev:) = qv_prime0(:ncol,top_lev:)
        1014           0 :    ql_dprime(:ncol,top_lev:) = ql_prime0(:ncol,top_lev:)
        1015           0 :    qi_dprime(:ncol,top_lev:) = qi_prime0(:ncol,top_lev:)
        1016           0 :    nl_dprime(:ncol,top_lev:) = nl_prime0(:ncol,top_lev:)
        1017           0 :    ni_dprime(:ncol,top_lev:) = ni_prime0(:ncol,top_lev:)
        1018             :
        1019             :    call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp,          &
        1020             :                            qv_dprime, ql_dprime, qi_dprime, T_dprime,  &
        1021           0 :                            qvten_pwi2, qlten_pwi2, qiten_pwi2, Tten_pwi2, do_cldice)
        1022             :
        1023           0 :    do k = top_lev, pver
        1024           0 :    do i = 1, ncol
        1025           0 :       if( ql_dprime(i,k) .lt. qsmall ) then
        1026           0 :           nlten_pwi2(i,k) = -nl_dprime(i,k)/dt
        1027           0 :           nl_dprime(i,k)   = 0._r8
        1028             :       endif
        1029           0 :       if( qi_dprime(i,k) .lt. qsmall ) then
        1030           0 :           niten_pwi2(i,k) = -ni_dprime(i,k)/dt
        1031           0 :           ni_dprime(i,k)   = 0._r8
        1032             :       endif
        1033             :    enddo
        1034             :    enddo
        1035             :
        1036             :    ! -------------------------------------------------------------- !
        1037             :    ! Add tendency associated with detrainment of cumulus condensate !
        1038             :    ! This tendency is not used in computing Q                       !
        1039             :    ! Since D_ql,D_qi,D_nl,D_ni > 0, don't need to worry about       !
        1040             :    ! negative scalar.                                               !
        1041             :    ! This tendency is not reflected into Fzs2, which is OK.         !
        1042             :    ! -------------------------------------------------------------- !
        1043             :
        1044           0 :    T_dprime(:ncol,top_lev:)   =  T_dprime(:ncol,top_lev:)  + D_T(:ncol,top_lev:) * dt
        1045           0 :    qv_dprime(:ncol,top_lev:)  = qv_dprime(:ncol,top_lev:) + D_qv(:ncol,top_lev:) * dt
        1046           0 :    ql_dprime(:ncol,top_lev:)  = ql_dprime(:ncol,top_lev:) + D_ql(:ncol,top_lev:) * dt
        1047           0 :    qi_dprime(:ncol,top_lev:)  = qi_dprime(:ncol,top_lev:) + D_qi(:ncol,top_lev:) * dt
        1048           0 :    nl_dprime(:ncol,top_lev:)  = nl_dprime(:ncol,top_lev:) + D_nl(:ncol,top_lev:) * dt
        1049           0 :    ni_dprime(:ncol,top_lev:)  = ni_dprime(:ncol,top_lev:) + D_ni(:ncol,top_lev:) * dt
        1050             :
        1051             :    ! ---------------------------------------------------------- !
        1052             :    ! Impose diagnostic upper and lower limits on the in-stratus !
        1053             :    ! condensate amount. This produces a final equilibrium state !
        1054             :    ! at the end of each iterative process.                      !
        1055             :    ! ---------------------------------------------------------- !
        1056             :
        1057           0 :    do k = top_lev, pver
        1058             :       call instratus_condensate( lchnk          , ncol           , k              , p(:,k)        , &
        1059             :                                  T_dprime(:,k)  , qv_dprime(:,k) , ql_dprime(:,k) , qi_dprime(:,k), &
        1060             :                                  ni_dprime(:,k) ,                                                   &
        1061             :                                  a_cu0(:,k)     , zeros(:,k)     , zeros(:,k)     ,                 &
        1062             :                                  zeros(:,k)     , zeros(:,k)     , zeros(:,k)     ,                 &
        1063             :                                  landfrac       , snowh          ,                                  &
        1064             :                                  rhmaxi_arr(:,k),rhmini_arr(:,k), rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k), &
        1065             :                                  T_star(:,k)    , qv_star(:,k)   , ql_star(:,k)   , qi_star(:,k)  , &
        1066           0 :                                  al_st_star(:,k), ai_st_star(:,k), ql_st_star(:,k), qi_st_star(:,k) )
        1067           0 :       a_st_star(:ncol,k)  = max(al_st_star(:ncol,k),ai_st_star(:ncol,k))
        1068           0 :       QQw2(:ncol,k) = (ql_star(:ncol,k) - ql_dprime(:ncol,k))/dt
        1069           0 :       QQi2(:ncol,k) = (qi_star(:ncol,k) - qi_dprime(:ncol,k))/dt
        1070             :       ! -------------------------------------------------- !
        1071             :       ! Reduce droplet concentration if evaporation occurs !
        1072             :       ! -------------------------------------------------- !
        1073           0 :       do i = 1, ncol
        1074           0 :          if( QQw2(i,k) .le. 0._r8 ) then
        1075           0 :              if( ql_dprime(i,k) .ge. qsmall ) then
        1076           0 :                  QQnl2(i,k) = QQw2(i,k)*nl_dprime(i,k)/ql_dprime(i,k)
        1077           0 :                  QQnl2(i,k) = min(0._r8,cone*max(QQnl2(i,k),-nl_dprime(i,k)/dt))
        1078             :              else
        1079           0 :                  QQnl2(i,k) = 0._r8
        1080             :              endif
        1081             :          endif
        1082           0 :          if( QQi2(i,k) .le. 0._r8 ) then
        1083           0 :              if( qi_dprime(i,k) .gt. qsmall ) then
        1084           0 :                  QQni2(i,k) = QQi2(i,k)*ni_dprime(i,k)/qi_dprime(i,k)
        1085           0 :                  QQni2(i,k) = min(0._r8,cone*max(QQni2(i,k),-ni_dprime(i,k)/dt))
        1086             :              else
        1087           0 :                  QQni2(i,k) = 0._r8
        1088             :              endif
        1089             :          endif
        1090             :       enddo
        1091             :    enddo
        1092           0 :    nl_star(:ncol,top_lev:) = max(0._r8,nl_dprime(:ncol,top_lev:)+QQnl2(:ncol,top_lev:)*dt)
        1093           0 :    ni_star(:ncol,top_lev:) = max(0._r8,ni_dprime(:ncol,top_lev:)+QQni2(:ncol,top_lev:)*dt)
        1094             :
        1095             :    ! ------------------------------------------ !
        1096             :    ! Final adjustment of droplet concentration. !
        1097             :    ! Set # to zero if there is no cloud.        !
        1098             :    ! ------------------------------------------ !
        1099             :
        1100           0 :    do k = top_lev, pver
        1101           0 :    do i = 1, ncol
        1102           0 :       if( ql_star(i,k) .lt. qsmall ) then
        1103           0 :           ACnl(i,k) = - nl_star(i,k)/dt
        1104           0 :           nl_star(i,k) = 0._r8
        1105             :       endif
        1106           0 :       if( qi_star(i,k) .lt. qsmall ) then
        1107           0 :           ACni(i,k) = - ni_star(i,k)/dt
        1108           0 :           ni_star(i,k) = 0._r8
        1109             :       endif
        1110             :    enddo
        1111             :    enddo
        1112             :
        1113             :    ! ----------------------------------------------------- !
        1114             :    ! Define equilibrium reference state for next iteration !
        1115             :    ! ----------------------------------------------------- !
        1116             :
        1117           0 :    T(:ncol,top_lev:)     = T_star(:ncol,top_lev:)
        1118           0 :    qv(:ncol,top_lev:)    = qv_star(:ncol,top_lev:)
        1119           0 :    ql(:ncol,top_lev:)    = ql_star(:ncol,top_lev:)
        1120           0 :    qi(:ncol,top_lev:)    = qi_star(:ncol,top_lev:)
        1121           0 :    al_st(:ncol,top_lev:) = al_st_star(:ncol,top_lev:)
        1122           0 :    ai_st(:ncol,top_lev:) = ai_st_star(:ncol,top_lev:)
        1123           0 :    a_st(:ncol,top_lev:)  = a_st_star(:ncol,top_lev:)
        1124           0 :    ql_st(:ncol,top_lev:) = ql_st_star(:ncol,top_lev:)
        1125           0 :    qi_st(:ncol,top_lev:) = qi_st_star(:ncol,top_lev:)
        1126           0 :    nl(:ncol,top_lev:)    = nl_star(:ncol,top_lev:)
        1127           0 :    ni(:ncol,top_lev:)    = ni_star(:ncol,top_lev:)
        1128             :
        1129             :    enddo ! End of 'iter' prognostic iterative computation
        1130             :
        1131             :    ! ------------------------------------------------------------------------ !
        1132             :    ! Compute final tendencies of main output variables and diagnostic outputs !
        1133             :    ! Note that the very input state [T0,qv0,ql0,qi0] are                      !
        1134             :    ! marched to [T_star,qv_star,ql_star,qi_star] with equilibrium             !
        1135             :    ! stratus informations of [a_st_star,ql_st_star,qi_st_star] by             !
        1136             :    ! below final tendencies and [A_T,A_qv,A_ql,A_qi]                          !
        1137             :    ! ------------------------------------------------------------------------ !
        1138             :
        1139             :    ! ------------------ !
        1140             :    ! Process tendencies !
        1141             :    ! ------------------ !
        1142             :
        1143           0 :    QQw_final(:ncol,top_lev:)  = QQw_prog(:ncol,top_lev:)
        1144           0 :    QQi_final(:ncol,top_lev:)  = QQi_prog(:ncol,top_lev:)
        1145           0 :    QQ_final(:ncol,top_lev:)   = QQw_final(:ncol,top_lev:) + QQi_final(:ncol,top_lev:)
        1146             :    QQw_all(:ncol,top_lev:)    = QQw_prog(:ncol,top_lev:)  + QQw1(:ncol,top_lev:) + QQw2(:ncol,top_lev:) + &
        1147           0 :         qlten_pwi1(:ncol,top_lev:) + qlten_pwi2(:ncol,top_lev:) + A_ql_adj(:ncol,top_lev:)
        1148             :    QQi_all(:ncol,top_lev:)    = QQi_prog(:ncol,top_lev:)  + QQi1(:ncol,top_lev:) + QQi2(:ncol,top_lev:) + &
        1149           0 :         qiten_pwi1(:ncol,top_lev:) + qiten_pwi2(:ncol,top_lev:) + A_qi_adj(:ncol,top_lev:)
        1150           0 :    QQ_all(:ncol,top_lev:)     = QQw_all(:ncol,top_lev:)   + QQi_all(:ncol,top_lev:)
        1151           0 :    QQnl_final(:ncol,top_lev:) = QQnl_prog(:ncol,top_lev:)
        1152           0 :    QQni_final(:ncol,top_lev:) = QQni_prog(:ncol,top_lev:)
        1153           0 :    QQn_final(:ncol,top_lev:)  = QQnl_final(:ncol,top_lev:) + QQni_final(:ncol,top_lev:)
        1154             :    QQnl_all(:ncol,top_lev:)   = QQnl_prog(:ncol,top_lev:)  + QQnl1(:ncol,top_lev:) + QQnl2(:ncol,top_lev:) + &
        1155           0 :         nlten_pwi1(:ncol,top_lev:) + nlten_pwi2(:ncol,top_lev:) + ACnl(:ncol,top_lev:) + A_nl_adj(:ncol,top_lev:)
        1156             :    QQni_all(:ncol,top_lev:)   = QQni_prog(:ncol,top_lev:)  + QQni1(:ncol,top_lev:) + QQni2(:ncol,top_lev:) + &
        1157           0 :         niten_pwi1(:ncol,top_lev:) + niten_pwi2(:ncol,top_lev:) + ACni(:ncol,top_lev:) + A_ni_adj(:ncol,top_lev:)
        1158           0 :    QQn_all(:ncol,top_lev:)    = QQnl_all(:ncol,top_lev:)   + QQni_all(:ncol,top_lev:)
        1159           0 :    qme(:ncol,top_lev:)        = QQ_final(:ncol,top_lev:)
        1160           0 :    qvadj(:ncol,top_lev:)      = qvten_pwi1(:ncol,top_lev:) + qvten_pwi2(:ncol,top_lev:) + A_qv_adj(:ncol,top_lev:)
        1161           0 :    qladj(:ncol,top_lev:)      = qlten_pwi1(:ncol,top_lev:) + qlten_pwi2(:ncol,top_lev:) + A_ql_adj(:ncol,top_lev:)
        1162           0 :    qiadj(:ncol,top_lev:)      = qiten_pwi1(:ncol,top_lev:) + qiten_pwi2(:ncol,top_lev:) + A_qi_adj(:ncol,top_lev:)
        1163           0 :    qllim(:ncol,top_lev:)      = QQw1      (:ncol,top_lev:) + QQw2      (:ncol,top_lev:)
        1164           0 :    qilim(:ncol,top_lev:)      = QQi1      (:ncol,top_lev:) + QQi2      (:ncol,top_lev:)
        1165             :
        1166             :    ! ----------------- !
        1167             :    ! Output tendencies !
        1168             :    ! ----------------- !
        1169             :
        1170             :    s_tendout(:ncol,top_lev:)  = cpair*( T_star(:ncol,top_lev:)  -  T0(:ncol,top_lev:) )/dt - &
        1171           0 :         cpair*(A_T(:ncol,top_lev:)+C_T(:ncol,top_lev:))
        1172             :    qv_tendout(:ncol,top_lev:) =    ( qv_star(:ncol,top_lev:) - qv0(:ncol,top_lev:) )/dt - &
        1173           0 :         (A_qv(:ncol,top_lev:)+C_qv(:ncol,top_lev:))
        1174             :    ql_tendout(:ncol,top_lev:) =    ( ql_star(:ncol,top_lev:) - ql0(:ncol,top_lev:) )/dt - &
        1175           0 :         (A_ql(:ncol,top_lev:)+C_ql(:ncol,top_lev:))
        1176             :    qi_tendout(:ncol,top_lev:) =    ( qi_star(:ncol,top_lev:) - qi0(:ncol,top_lev:) )/dt - &
        1177           0 :         (A_qi(:ncol,top_lev:)+C_qi(:ncol,top_lev:))
        1178             :    nl_tendout(:ncol,top_lev:) =    ( nl_star(:ncol,top_lev:) - nl0(:ncol,top_lev:) )/dt - &
        1179           0 :         (A_nl(:ncol,top_lev:)+C_nl(:ncol,top_lev:))
        1180             :    ni_tendout(:ncol,top_lev:) =    ( ni_star(:ncol,top_lev:) - ni0(:ncol,top_lev:) )/dt - &
        1181           0 :         (A_ni(:ncol,top_lev:)+C_ni(:ncol,top_lev:))
        1182             :
        1183           0 :    if (.not. do_cldice) then
        1184           0 :       do k = top_lev, pver
        1185           0 :          do i = 1, ncol
        1186             :
        1187             :             ! Don't want either qi or ni tendencies, but the code above is somewhat convoluted and
        1188             :             ! is trying to adjust both (small numbers). Just force it to zero here.
        1189           0 :             qi_tendout(i,k) = 0._r8
        1190           0 :             ni_tendout(i,k) = 0._r8
        1191             :           end do
        1192             :       end do
        1193             :    end if
        1194             :
        1195             :    ! ------------------ !
        1196             :    ! Net cloud fraction !
        1197             :    ! ------------------ !
        1198             :
        1199           0 :    cld(:ncol,top_lev:) = a_st_star(:ncol,top_lev:) + a_cu0(:ncol,top_lev:)
        1200             :
        1201             :    ! --------------------------------- !
        1202             :    ! Updated grid-mean state variables !
        1203             :    ! --------------------------------- !
        1204             :
        1205           0 :    T0(:ncol,top_lev:)  = T_star(:ncol,top_lev:)
        1206           0 :    qv0(:ncol,top_lev:) = qv_star(:ncol,top_lev:)
        1207           0 :    ql0(:ncol,top_lev:) = ql_star(:ncol,top_lev:)
        1208           0 :    qi0(:ncol,top_lev:) = qi_star(:ncol,top_lev:)
        1209           0 :    nl0(:ncol,top_lev:) = nl_star(:ncol,top_lev:)
        1210           0 :    ni0(:ncol,top_lev:) = ni_star(:ncol,top_lev:)
        1211             :
        1212           0 :    if (hist_fld_active('RHMIN_LIQ')) then
        1213             :       ! Compute default critical RH as a function of height and surface type as in the current code.
        1214           0 :       rhmin_liq_diag(:,:) = 0._r8
        1215           0 :       do k = top_lev, pver
        1216           0 :          do i = 1, ncol
        1217           0 :             land = nint(landfrac(i)) == 1
        1218           0 :             if( p(i,k) .ge. premib ) then
        1219           0 :                if( land .and. (snowh(i).le.0.000001_r8) ) then
        1220           0 :                   rhmin_liq_diag(i,k) = rhminl_const - rhminl_adj_land_const
        1221             :                else
        1222           0 :                   rhmin_liq_diag(i,k) = rhminl_const
        1223             :                endif
        1224           0 :             elseif( p(i,k) .lt. premit ) then
        1225           0 :                rhmin_liq_diag(i,k) = rhminh_const
        1226             :             else
        1227           0 :                tmp = (premib-(max(p(i,k),premit)))/(premib-premit)
        1228           0 :                rhmin_liq_diag(i,k) = rhminh_const*tmp + rhminl_const*(1.0_r8-tmp)
        1229             :             endif
        1230             :          end do
        1231             :       end do
        1232           0 :       call outfld( 'RHMIN_LIQ',      rhmin_liq_diag,  pcols, lchnk )
        1233             :    end if
        1234             :
        1235           0 :    rhmin_ice_diag(:,:) = rhminh_const
        1236           0 :    call outfld( 'RHMIN_ICE',      rhmin_ice_diag,  pcols, lchnk )
        1237             :
        1238           0 :    call outfld( 'DRHMINPBL_LIQ', d_rhmin_liq_PBL,  pcols, lchnk )
        1239           0 :    call outfld( 'DRHMINPBL_ICE', d_rhmin_ice_PBL,  pcols, lchnk )
        1240           0 :    call outfld( 'DRHMINDET_LIQ', d_rhmin_liq_det,  pcols, lchnk )
        1241           0 :    call outfld( 'DRHMINDET_ICE', d_rhmin_ice_det,  pcols, lchnk )
        1242             :
        1243           0 :    end subroutine mmacro_pcond
        1244             :
        1245             :
        1246             : !=======================================================================================================
        1247             :
        1248           0 : subroutine rhcrit_calc( &
        1249             :    ncol, dp, T0, p, &
        1250             :    clrw_old, clri_old, tke, qtl_flx, &
        1251             :    qti_flx, cmfr_det, qlr_det, qir_det, &
        1252             :    rhmaxi_arr, rhmini_arr, rhminl_arr, rhminl_adj_land_arr, rhminh_arr, &
        1253             :    d_rhmin_liq_PBL, d_rhmin_ice_PBL, d_rhmin_liq_det, d_rhmin_ice_det)
        1254             :
        1255             :    ! ------------------------------------------------- !
        1256             :    ! Compute a drop of critical RH for stratus by      !
        1257             :    ! (1) PBL turbulence, and                           !
        1258             :    ! (2) convective detrainment.                       !
        1259             :    ! Note that all of 'd_rhmin...' terms are positive. !
        1260             :    ! ------------------------------------------------- !
        1261             :
        1262             :    integer,  intent(in) :: ncol                         ! Number of active columns
        1263             :    real(r8), intent(in) :: dp(pcols,pver)               ! Pressure thickness [Pa] > 0
        1264             :    real(r8), intent(in) :: T0(pcols,pver)               ! Temperature [K]
        1265             :    real(r8), intent(in) :: p(pcols,pver)                ! Pressure at the layer mid-point [Pa]
        1266             :    real(r8), intent(in) :: clrw_old(pcols,pver)         ! Clear sky fraction at the previous time step for liquid stratus process
        1267             :    real(r8), intent(in) :: clri_old(pcols,pver)         ! Clear sky fraction at the previous time step for    ice stratus process
        1268             :    real(r8), pointer    :: tke(:,:)                     ! (pcols,pverp) TKE from the PBL scheme
        1269             :    real(r8), pointer    :: qtl_flx(:,:)                 ! (pcols,pverp) overbar(w'qtl') from PBL scheme where qtl = qv + ql
        1270             :    real(r8), pointer    :: qti_flx(:,:)                 ! (pcols,pverp) overbar(w'qti') from PBL scheme where qti = qv + qi
        1271             :    real(r8), pointer    :: cmfr_det(:,:)                ! (pcols,pver)  Detrained mass flux from the convection scheme
        1272             :    real(r8), pointer    :: qlr_det(:,:)                 ! (pcols,pver)  Detrained        ql from the convection scheme
        1273             :    real(r8), pointer    :: qir_det(:,:)                 ! (pcols,pver)  Detrained        qi from the convection scheme
        1274             :
        1275             :    real(r8), intent(in)  :: rhmaxi_arr(pcols,pver)
        1276             :
        1277             :    real(r8), intent(out) :: rhmini_arr(pcols,pver)
        1278             :    real(r8), intent(out) :: rhminl_arr(pcols,pver)
        1279             :    real(r8), intent(out) :: rhminl_adj_land_arr(pcols,pver)
        1280             :    real(r8), intent(out) :: rhminh_arr(pcols,pver)
        1281             :    real(r8), intent(out) :: d_rhmin_liq_PBL(pcols,pver)
        1282             :    real(r8), intent(out) :: d_rhmin_ice_PBL(pcols,pver)
        1283             :    real(r8), intent(out) :: d_rhmin_liq_det(pcols,pver)
        1284             :    real(r8), intent(out) :: d_rhmin_ice_det(pcols,pver)
        1285             :
        1286             :    ! local variables
        1287             :
        1288             :    integer :: i, k
        1289             :
        1290             :    real(r8) :: esat_tmp(pcols)          ! Dummy for saturation vapor pressure calc.
        1291             :    real(r8) :: qsat_tmp(pcols)          ! Saturation water vapor specific humidity [kg/kg]
        1292             :    real(r8) :: sig_tmp
        1293             :    !---------------------------------------------------------------------------------------------------
        1294             :
        1295             :
        1296             :
        1297             :    ! ---------------------------------- !
        1298             :    ! Calc critical RH for ice stratus   !
        1299             :    ! ---------------------------------- !
        1300             :
        1301           0 :    rhmini_arr(:,:) = rhmini_const
        1302             :
        1303           0 :    if (i_rhmini > 0) then
        1304             :
        1305             :       ! Compute the drop of critical RH by convective detrainment of cloud condensate
        1306             :
        1307           0 :       do k = top_lev, pver
        1308           0 :          do i = 1, ncol
        1309           0 :             d_rhmin_ice_det(i,k) = tau_deti*(gravit/dp(i,k))*cmfr_det(i,k)*clri_old(i,k)*qir_det(i,k)*3.6e6_r8
        1310           0 :             d_rhmin_ice_det(i,k) = max(0._r8,min(0.5_r8,d_rhmin_ice_det(i,k)))
        1311             :          end do
        1312             :       end do
        1313             :
        1314           0 :       if (i_rhmini == 1) then
        1315           0 :          rhmini_arr(:ncol,:) = rhmini_const - d_rhmin_ice_det(:ncol,:)
        1316             :       end if
        1317             :
        1318             :    end if
        1319             :
        1320           0 :    if (i_rhmini == 2) then
        1321             :
        1322             :       ! Compute the drop of critical RH by the variability induced by PBL turbulence
        1323           0 :       do k = top_lev, pver
        1324           0 :          call qsat_ice(T0(1:ncol,k), p(1:ncol,k), esat_tmp(1:ncol), qsat_tmp(1:ncol), ncol)
        1325           0 :          do i = 1, ncol
        1326           0 :             sig_tmp = 0.5_r8 * ( qti_flx(i,k)   / sqrt(max(qsmall,tke(i,k))) + &
        1327           0 :                                  qti_flx(i,k+1) / sqrt(max(qsmall,tke(i,k+1))) )
        1328           0 :             d_rhmin_ice_PBL(i,k) = c_aniso*sig_tmp/max(qsmall,qsat_tmp(i))
        1329           0 :             d_rhmin_ice_PBL(i,k) = max(0._r8,min(0.5_r8,d_rhmin_ice_PBL(i,k)))
        1330             :
        1331           0 :             rhmini_arr(i,k) = 1._r8 - d_rhmin_ice_PBL(i,k) - d_rhmin_ice_det(i,k)
        1332             :          end do
        1333             :       end do
        1334             :    end if
        1335             :
        1336           0 :    if (i_rhmini > 0) then
        1337           0 :       do k = top_lev, pver
        1338           0 :          do i = 1, ncol
        1339           0 :             rhmini_arr(i,k) = max(0._r8,min(rhmaxi_arr(i,k),rhmini_arr(i,k)))
        1340             :          end do
        1341             :       end do
        1342             :    end if
        1343             :
        1344             :    ! ------------------------------------- !
        1345             :    ! Choose critical RH for liquid stratus !
        1346             :    ! ------------------------------------- !
        1347             :
        1348           0 :    rhminl_arr(:,:)          = rhminl_const
        1349           0 :    rhminl_adj_land_arr(:,:) = rhminl_adj_land_const
        1350           0 :    rhminh_arr(:,:)          = rhminh_const
        1351             :
        1352           0 :    if (i_rhminl > 0) then
        1353             :
        1354             :       ! Compute the drop of critical RH by convective detrainment of cloud condensate
        1355             :
        1356           0 :       do k = top_lev, pver
        1357           0 :          do i = 1, ncol
        1358           0 :             d_rhmin_liq_det(i,k) = tau_detw*(gravit/dp(i,k))*cmfr_det(i,k)*clrw_old(i,k)*qlr_det(i,k)*3.6e6_r8
        1359           0 :             d_rhmin_liq_det(i,k) = max(0._r8,min(0.5_r8,d_rhmin_liq_det(i,k)))
        1360             :          end do
        1361             :       end do
        1362             :
        1363           0 :       if (i_rhminl == 1) then
        1364           0 :          rhminl_arr(:ncol,top_lev:) = rhminl_const - d_rhmin_liq_det(:ncol,top_lev:)
        1365           0 :          rhminh_arr(:ncol,top_lev:) = rhminh_const - d_rhmin_liq_det(:ncol,top_lev:)
        1366             :       end if
        1367             :
        1368             :    end if
        1369             :
        1370           0 :    if (i_rhminl == 2) then
        1371             :
        1372             :       ! Compute the drop of critical RH by the variability induced by PBL turbulence
        1373           0 :       do k = top_lev, pver
        1374           0 :          call qsat_water(T0(1:ncol,k), p(1:ncol,k), esat_tmp(1:ncol), qsat_tmp(1:ncol), ncol)
        1375           0 :          do i = 1, ncol
        1376           0 :             sig_tmp = 0.5_r8 * ( qtl_flx(i,k)   / sqrt(max(qsmall,tke(i,k))) + &
        1377           0 :                                  qtl_flx(i,k+1) / sqrt(max(qsmall,tke(i,k+1))) )
        1378           0 :             d_rhmin_liq_PBL(i,k) = c_aniso*sig_tmp/max(qsmall,qsat_tmp(i))
        1379           0 :             d_rhmin_liq_PBL(i,k) = max(0._r8,min(0.5_r8,d_rhmin_liq_PBL(i,k)))
        1380             :
        1381           0 :             rhminl_arr(i,k) = 1._r8 - d_rhmin_liq_PBL(i,k) - d_rhmin_liq_det(i,k)
        1382           0 :             rhminl_adj_land_arr(i,k) = 0._r8
        1383           0 :             rhminh_arr(i,k) = rhminl_arr(i,k)
        1384             :          end do
        1385             :       end do
        1386             :    end if
        1387             :
        1388           0 :    if (i_rhminl > 0) then
        1389           0 :       do k = top_lev, pver
        1390           0 :          do i = 1, ncol
        1391           0 :             rhminl_arr(i,k) = max(rhminl_adj_land_arr(i,k),min(1._r8,rhminl_arr(i,k)))
        1392           0 :             rhminh_arr(i,k) = max(0._r8,min(1._r8,rhminh_arr(i,k)))
        1393             :          end do
        1394             :       end do
        1395             :    end if
        1396             :
        1397           0 : end subroutine rhcrit_calc
        1398             :
        1399             : !=======================================================================================================
        1400             :
        1401           0 :    subroutine instratus_condensate( lchnk, ncol, k,                      &
        1402             :                                     p_in, T0_in, qv0_in, ql0_in, qi0_in, &
        1403             :                                     ni0_in,                              &
        1404             :                                     a_dc_in, ql_dc_in, qi_dc_in,         &
        1405             :                                     a_sc_in, ql_sc_in, qi_sc_in,         &
        1406             :                                     landfrac, snowh,                     &
        1407             :                                     rhmaxi_in, rhmini_in, rhminl_in, rhminl_adj_land_in, rhminh_in, &
        1408             :                                     T_out, qv_out, ql_out, qi_out,       &
        1409             :                                     al_st_out, ai_st_out, ql_st_out, qi_st_out )
        1410             :
        1411             :    ! ------------------------------------------------------- !
        1412             :    ! Diagnostically force in-stratus condensate to be        !
        1413             :    ! in the range of 'qlst_min < qc_st < qlst_max'           !
        1414             :    ! whenever stratus exists in the equilibrium state        !
        1415             :    ! ------------------------------------------------------- !
        1416             :
        1417             :    integer,  intent(in)  :: lchnk                ! Chunk identifier
        1418             :    integer,  intent(in)  :: ncol                 ! Number of atmospheric columns
        1419             :    integer,  intent(in)  :: k                    ! Layer index
        1420             :
        1421             :    real(r8), intent(in)  :: p_in(pcols)          ! Pressure [Pa]
        1422             :    real(r8), intent(in)  :: T0_in(pcols)         ! Temperature [K]
        1423             :    real(r8), intent(in)  :: qv0_in(pcols)        ! Grid-mean water vapor [kg/kg]
        1424             :    real(r8), intent(in)  :: ql0_in(pcols)        ! Grid-mean LWC [kg/kg]
        1425             :    real(r8), intent(in)  :: qi0_in(pcols)        ! Grid-mean IWC [kg/kg]
        1426             :    real(r8), intent(in)  :: ni0_in(pcols)
        1427             :
        1428             :    real(r8), intent(in)  :: a_dc_in(pcols)       ! Deep cumulus cloud fraction
        1429             :    real(r8), intent(in)  :: ql_dc_in(pcols)      ! In-deep cumulus LWC [kg/kg]
        1430             :    real(r8), intent(in)  :: qi_dc_in(pcols)      ! In-deep cumulus IWC [kg/kg]
        1431             :    real(r8), intent(in)  :: a_sc_in(pcols)       ! Shallow cumulus cloud fraction
        1432             :    real(r8), intent(in)  :: ql_sc_in(pcols)      ! In-shallow cumulus LWC [kg/kg]
        1433             :    real(r8), intent(in)  :: qi_sc_in(pcols)      ! In-shallow cumulus IWC [kg/kg]
        1434             :
        1435             :    real(r8), intent(in)  :: landfrac(pcols)      ! Land fraction
        1436             :    real(r8), intent(in)  :: snowh(pcols)         ! Snow depth (liquid water equivalent)
        1437             :
        1438             :    real(r8), intent(in)  :: rhmaxi_in(pcols)
        1439             :    real(r8), intent(in)  :: rhmini_in(pcols)
        1440             :    real(r8), intent(in)  :: rhminl_in(pcols)
        1441             :    real(r8), intent(in)  :: rhminl_adj_land_in(pcols)
        1442             :    real(r8), intent(in)  :: rhminh_in(pcols)
        1443             :
        1444             :    real(r8), intent(out) :: T_out(pcols)         ! Temperature [K]
        1445             :    real(r8), intent(out) :: qv_out(pcols)        ! Grid-mean water vapor [kg/kg]
        1446             :    real(r8), intent(out) :: ql_out(pcols)        ! Grid-mean LWC [kg/kg]
        1447             :    real(r8), intent(out) :: qi_out(pcols)        ! Grid-mean IWC [kg/kg]
        1448             :
        1449             :    real(r8), intent(out) :: al_st_out(pcols)     ! Liquid stratus fraction
        1450             :    real(r8), intent(out) :: ai_st_out(pcols)     ! Ice stratus fraction
        1451             :    real(r8), intent(out) :: ql_st_out(pcols)     ! In-stratus LWC [kg/kg]
        1452             :    real(r8), intent(out) :: qi_st_out(pcols)     ! In-stratus IWC [kg/kg]
        1453             :
        1454             :    ! Local variables
        1455             :
        1456             :    integer i                                     ! Column    index
        1457             :
        1458             :    real(r8) p
        1459             :    real(r8) T0
        1460             :    real(r8) qv0
        1461             :    real(r8) ql0
        1462             :    real(r8) qi0
        1463             :    real(r8) a_dc
        1464             :    real(r8) ql_dc
        1465             :    real(r8) qi_dc
        1466             :    real(r8) a_sc
        1467             :    real(r8) ql_sc
        1468             :    real(r8) qi_sc
        1469             :    real(r8) esat0
        1470             :    real(r8) qsat0
        1471             :    real(r8) U0
        1472             :    real(r8) U0_nc
        1473             :    real(r8) G0_nc
        1474             :    real(r8) al0_st_nc
        1475             :    real(r8) al0_st
        1476             :    real(r8) ai0_st_nc
        1477             :    real(r8) ai0_st
        1478             :    real(r8) a0_st
        1479             :    real(r8) ql0_nc
        1480             :    real(r8) qi0_nc
        1481             :    real(r8) qc0_nc
        1482             :    real(r8) ql0_st
        1483             :    real(r8) qi0_st
        1484             :    real(r8) qc0_st
        1485             :    real(r8) T
        1486             :    real(r8) qv
        1487             :    real(r8) ql
        1488             :    real(r8) qi
        1489             :    real(r8) ql_st
        1490             :    real(r8) qi_st
        1491             :    real(r8) es
        1492             :    real(r8) qs
        1493             :    real(r8) esat_in(pcols)
        1494             :    real(r8) qsat_in(pcols)
        1495             :    real(r8) U0_in(pcols)
        1496             :    real(r8) al0_st_nc_in(pcols)
        1497             :    real(r8) ai0_st_nc_in(pcols)
        1498             :    real(r8) G0_nc_in(pcols)
        1499             :    integer  idxmod
        1500             :    real(r8) U
        1501             :    real(r8) U_nc
        1502             :    real(r8) al_st_nc
        1503             :    real(r8) ai_st_nc
        1504             :    real(r8) G_nc
        1505             :    real(r8) a_st
        1506             :    real(r8) al_st
        1507             :    real(r8) ai_st
        1508             :    real(r8) Tmin0
        1509             :    real(r8) Tmax0
        1510             :    real(r8) Tmin
        1511             :    real(r8) Tmax
        1512             :    integer caseid
        1513             :
        1514             :    real(r8) rhmaxi
        1515             :    real(r8) rhmini
        1516             :    real(r8) rhminl
        1517             :    real(r8) rhminl_adj_land
        1518             :    real(r8) rhminh
        1519             :
        1520             :    ! ---------------- !
        1521             :    ! Main Computation !
        1522             :    ! ---------------- !
        1523             :
        1524           0 :    call qsat_water(T0_in(1:ncol), p_in(1:ncol), esat_in(1:ncol), qsat_in(1:ncol), ncol)
        1525           0 :    U0_in(:ncol) = qv0_in(:ncol)/qsat_in(:ncol)
        1526           0 :    if( CAMstfrac ) then
        1527             :        call astG_RHU(U0_in(:),p_in(:),qv0_in(:),landfrac(:),snowh(:),al0_st_nc_in(:),G0_nc_in(:),ncol,&
        1528             :                      rhminl_in(:), rhminl_adj_land_in(:), rhminh_in(:))
        1529             :    else
        1530             :        call astG_PDF(U0_in(:),p_in(:),qv0_in(:),landfrac(:),snowh(:),al0_st_nc_in(:),G0_nc_in(:),ncol,&
        1531             :                      rhminl_in(:), rhminl_adj_land_in(:), rhminh_in(:))
        1532             :    endif
        1533             :    call aist_vector(qv0_in(:),T0_in(:),p_in(:),qi0_in(:),ni0_in(:),landfrac(:),snowh(:),ai0_st_nc_in(:),ncol,&
        1534           0 :                     rhmaxi_in(:), rhmini_in(:), rhminl_in(:), rhminl_adj_land_in(:), rhminh_in(:))
        1535             :
        1536           0 :    do i = 1, ncol
        1537             :
        1538             :       ! ---------------------- !
        1539             :       ! Define local variables !
        1540             :       ! ---------------------- !
        1541             :
        1542           0 :       p   = p_in(i)
        1543             :
        1544           0 :       T0  = T0_in(i)
        1545           0 :       qv0 = qv0_in(i)
        1546           0 :       ql0 = ql0_in(i)
        1547           0 :       qi0 = qi0_in(i)
        1548             :
        1549           0 :       a_dc  = a_dc_in(i)
        1550           0 :       ql_dc = ql_dc_in(i)
        1551           0 :       qi_dc = qi_dc_in(i)
        1552             :
        1553           0 :       a_sc  = a_sc_in(i)
        1554           0 :       ql_sc = ql_sc_in(i)
        1555           0 :       qi_sc = qi_sc_in(i)
        1556             :
        1557           0 :       ql_dc = 0._r8
        1558           0 :       qi_dc = 0._r8
        1559           0 :       ql_sc = 0._r8
        1560           0 :       qi_sc = 0._r8
        1561             :
        1562           0 :       es  = esat_in(i)
        1563           0 :       qs  = qsat_in(i)
        1564             :
        1565           0 :       rhmaxi = rhmaxi_in(i)
        1566           0 :       rhmini = rhmini_in(i)
        1567           0 :       rhminl = rhminl_in(i)
        1568           0 :       rhminl_adj_land = rhminl_adj_land_in(i)
        1569           0 :       rhminh = rhminh_in(i)
        1570             :
        1571           0 :       idxmod = 0
        1572           0 :       caseid = -1
        1573             :
        1574             :       ! ------------------------------------------------------------ !
        1575             :       ! Force the grid-mean RH to be smaller than 1 if oversaturated !
        1576             :       ! In order to be compatible with reduced 3x3 QQ, condensation  !
        1577             :       ! should occur only into the liquid in gridmean_RH.            !
        1578             :       ! ------------------------------------------------------------ !
        1579             :
        1580           0 :       if( qv0 .gt. qs ) then
        1581             :           call gridmean_RH( lchnk, i, k, p, T0, qv0, ql0, qi0,      &
        1582             :                             a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, &
        1583             :                             landfrac(i), snowh(i) )
        1584           0 :           call qsat_water(T0, p, esat0, qsat0)
        1585           0 :           U0      = (qv0/qsat0)
        1586           0 :           U0_nc   =  U0
        1587             :           if( CAMstfrac ) then
        1588             :               call astG_RHU_single(U0_nc, p, qv0, landfrac(i), snowh(i), al0_st_nc, G0_nc, &
        1589             :                  rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
        1590             :           else
        1591           0 :               call astG_PDF_single(U0_nc, p, qv0, landfrac(i), snowh(i), al0_st_nc, G0_nc, &
        1592           0 :                  rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
        1593             :           endif
        1594           0 :           call aist_single(qv0,T0,p,qi0,landfrac(i),snowh(i),ai0_st_nc,&
        1595           0 :                            rhmaxi, rhmini, rhminl, rhminl_adj_land, rhminh)
        1596           0 :           ai0_st  = (1._r8-a_dc-a_sc)*ai0_st_nc
        1597           0 :           al0_st  = (1._r8-a_dc-a_sc)*al0_st_nc
        1598           0 :           a0_st   = max(ai0_st,al0_st)
        1599           0 :           idxmod  = 1
        1600             :       else
        1601           0 :           ai0_st  = (1._r8-a_dc-a_sc)*ai0_st_nc_in(i)
        1602           0 :           al0_st  = (1._r8-a_dc-a_sc)*al0_st_nc_in(i)
        1603             :       endif
        1604           0 :       a0_st   = max(ai0_st,al0_st)
        1605             :
        1606             :       ! ----------------------- !
        1607             :       ! Handling of input state !
        1608             :       ! ----------------------- !
        1609             :
        1610           0 :       ql0_nc  = max(0._r8,ql0-a_dc*ql_dc-a_sc*ql_sc)
        1611           0 :       qi0_nc  = max(0._r8,qi0-a_dc*qi_dc-a_sc*qi_sc)
        1612           0 :       qc0_nc  = ql0_nc + qi0_nc
        1613             :
        1614           0 :       Tmin0 = T0 - (latvap/cpair)*ql0
        1615           0 :       Tmax0 = T0 + ((latvap+latice)/cpair)*qv0
        1616             :
        1617             :       ! ------------------------------------------------------------- !
        1618             :       ! Do nothing and just exit if generalized in-stratus condensate !
        1619             :       ! condition is satisfied. This includes the case I.             !
        1620             :       ! For 4x4 liquid stratus, a0_st --> al0_st.                     !
        1621             :       ! ------------------------------------------------------------- !
        1622           0 :       if( ( ql0_nc .ge. qlst_min*al0_st ) .and. ( ql0_nc .le. qlst_max*al0_st ) ) then
        1623             :
        1624             :           ! ------------------ !
        1625             :           ! This is the case I !
        1626             :           ! ------------------ !
        1627           0 :              T = T0
        1628           0 :              qv = qv0
        1629           0 :              ql = ql0
        1630           0 :              qi = qi0
        1631           0 :              caseid = 0
        1632           0 :              goto 10
        1633             :       else
        1634             :          ! ----------------------------- !
        1635             :          ! This is case II : Dense Cloud !
        1636             :          ! ----------------------------- !
        1637           0 :          if( al0_st .eq. 0._r8 .and. ql0_nc .gt. 0._r8 ) then
        1638             :              ! ------------------------------------- !
        1639             :              ! Compute hypothetical full evaporation !
        1640             :              ! ------------------------------------- !
        1641           0 :              T  = Tmin0
        1642           0 :              qv = qv0 + ql0
        1643           0 :              call qsat_water(T, p, es, qs)
        1644           0 :              U  = qv/qs
        1645           0 :              U_nc = U
        1646             :              if( CAMstfrac ) then
        1647             :                  call astG_RHU_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, &
        1648             :                     rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
        1649             :              else
        1650           0 :                  call astG_PDF_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, &
        1651           0 :                     rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
        1652             :              endif
        1653           0 :              al_st = (1._r8-a_dc-a_sc)*al_st_nc
        1654           0 :              caseid = 0
        1655             :
        1656           0 :              if( al_st .eq. 0._r8 ) then
        1657           0 :                  ql = 0._r8
        1658           0 :                  qi = qi0
        1659           0 :                  idxmod = 1
        1660           0 :                  caseid = 1
        1661           0 :                  goto 10
        1662             :              else
        1663             :                  ! ------------------------------------------- !
        1664             :                  ! Evaporate until qc_st decreases to qlst_max !
        1665             :                  ! ------------------------------------------- !
        1666           0 :                  Tmin = Tmin0
        1667           0 :                  Tmax = T0
        1668             :                  call instratus_core( lchnk, i, k, p,                              &
        1669             :                                       T0, qv0, ql0, 0._r8,                         &
        1670             :                                       a_dc, ql_dc, qi_dc,                          &
        1671             :                                       a_sc, ql_sc, qi_sc, ai0_st,                  &
        1672           0 :                                       qlst_max, Tmin, Tmax, landfrac(i), snowh(i), &
        1673             :                                       rhminl, rhminl_adj_land, rhminh,             &
        1674           0 :                                       T, qv, ql, qi )
        1675           0 :                  idxmod = 1
        1676           0 :                  caseid = 2
        1677           0 :                  goto 10
        1678             :              endif
        1679             :          ! ------------------------------ !
        1680             :          ! This is case III : Empty Cloud !
        1681             :          ! ------------------------------ !
        1682           0 :          elseif( al0_st .gt. 0._r8 .and. ql0_nc .eq. 0._r8 ) then
        1683             :               ! ------------------------------------------ !
        1684             :               ! Condense until qc_st increases to qlst_min !
        1685             :               ! ------------------------------------------ !
        1686           0 :               Tmin = Tmin0
        1687           0 :               Tmax = Tmax0
        1688             :               call instratus_core( lchnk, i, k, p,                              &
        1689             :                                    T0, qv0, ql0, 0._r8,                         &
        1690             :                                    a_dc, ql_dc, qi_dc,                          &
        1691             :                                    a_sc, ql_sc, qi_sc, ai0_st,                  &
        1692           0 :                                    qlst_min, Tmin, Tmax, landfrac(i), snowh(i), &
        1693             :                                    rhminl, rhminl_adj_land, rhminh,             &
        1694           0 :                                    T, qv, ql, qi )
        1695           0 :               idxmod = 1
        1696           0 :               caseid = 3
        1697           0 :               goto 10
        1698             :          ! --------------- !
        1699             :          ! This is case IV !
        1700             :          ! --------------- !
        1701           0 :          elseif( al0_st .gt. 0._r8 .and. ql0_nc .gt. 0._r8 ) then
        1702             :
        1703           0 :              if( ql0_nc .gt. qlst_max*al0_st ) then
        1704             :                  ! --------------------------------------- !
        1705             :                  ! Evaporate until qc_st drops to qlst_max !
        1706             :                  ! --------------------------------------- !
        1707           0 :                  Tmin = Tmin0
        1708           0 :                  Tmax = Tmax0
        1709             :                  call instratus_core( lchnk, i, k, p,                              &
        1710             :                                       T0, qv0, ql0, 0._r8,                         &
        1711             :                                       a_dc, ql_dc, qi_dc,                          &
        1712             :                                       a_sc, ql_sc, qi_sc, ai0_st,                  &
        1713           0 :                                       qlst_max, Tmin, Tmax, landfrac(i), snowh(i), &
        1714             :                                       rhminl, rhminl_adj_land, rhminh,             &
        1715           0 :                                       T, qv, ql, qi )
        1716           0 :                  idxmod = 1
        1717           0 :                  caseid = 4
        1718           0 :                  goto 10
        1719           0 :              elseif( ql0_nc .lt. qlst_min*al0_st ) then
        1720             :                  ! -------------------------------------------- !
        1721             :                  ! Condensate until qc_st increases to qlst_min !
        1722             :                  ! -------------------------------------------- !
        1723           0 :                  Tmin = Tmin0
        1724           0 :                  Tmax = Tmax0
        1725             :                  call instratus_core( lchnk, i, k, p,                              &
        1726             :                                       T0, qv0, ql0, 0._r8,                         &
        1727             :                                       a_dc, ql_dc, qi_dc,                          &
        1728             :                                       a_sc, ql_sc, qi_sc, ai0_st,                  &
        1729           0 :                                       qlst_min, Tmin, Tmax, landfrac(i), snowh(i), &
        1730             :                                       rhminl, rhminl_adj_land, rhminh,             &
        1731           0 :                                       T, qv, ql, qi )
        1732           0 :                  idxmod = 1
        1733           0 :                  caseid = 5
        1734           0 :                  goto 10
        1735             :              else
        1736             :                  ! ------------------------------------------------ !
        1737             :                  ! This case should not happen. Issue error message !
        1738             :                  ! ------------------------------------------------ !
        1739           0 :                  write(iulog,*) 'Impossible case1 in instratus_condensate'
        1740           0 :                  call endrun
        1741             :              endif
        1742             :          ! ------------------------------------------------ !
        1743             :          ! This case should not happen. Issue error message !
        1744             :          ! ------------------------------------------------ !
        1745             :          else
        1746           0 :              write(iulog,*) 'Impossible case2 in instratus_condensate'
        1747           0 :              write(iulog,*)  al0_st, a_sc, a_dc
        1748           0 :              write(iulog,*)  1000*ql0_nc, 1000*(ql0+qi0)
        1749           0 :              call endrun
        1750             :          endif
        1751             :       endif
        1752             :
        1753             : 10 continue
        1754             :
        1755             :    ! -------------------------------------------------- !
        1756             :    ! Force final energy-moisture conserving consistency !
        1757             :    ! -------------------------------------------------- !
        1758             :
        1759           0 :      qi = qi0
        1760             :
        1761           0 :      if( idxmod .eq. 1 ) then
        1762           0 :          call aist_single(qv,T,p,qi,landfrac(i),snowh(i),ai_st_nc,&
        1763           0 :                           rhmaxi, rhmini, rhminl, rhminl_adj_land, rhminh)
        1764           0 :          ai_st = (1._r8-a_dc-a_sc)*ai_st_nc
        1765           0 :          call qsat_water(T, p, es, qs)
        1766           0 :          U     = (qv/qs)
        1767           0 :          U_nc  =  U
        1768             :          if( CAMstfrac ) then
        1769             :              call astG_RHU_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, &
        1770             :                 rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
        1771             :          else
        1772           0 :              call astG_PDF_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, &
        1773           0 :                 rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
        1774             :          endif
        1775           0 :          al_st = (1._r8-a_dc-a_sc)*al_st_nc
        1776             :      else
        1777           0 :          ai_st  = (1._r8-a_dc-a_sc)*ai0_st_nc_in(i)
        1778           0 :          al_st  = (1._r8-a_dc-a_sc)*al0_st_nc_in(i)
        1779             :      endif
        1780             :
        1781           0 :      a_st  = max(ai_st,al_st)
        1782             :
        1783           0 :      if( al_st .eq. 0._r8 ) then
        1784             :          ql_st = 0._r8
        1785             :      else
        1786           0 :          ql_st = ql/al_st
        1787           0 :          ql_st = min(qlst_max,max(qlst_min,ql_st)) ! PJR
        1788             :      endif
        1789           0 :      if( ai_st .eq. 0._r8 ) then
        1790             :          qi_st = 0._r8
        1791             :      else
        1792           0 :          qi_st = qi/ai_st
        1793             :      endif
        1794             :
        1795           0 :      qi    = ai_st*qi_st
        1796           0 :      ql    = al_st*ql_st
        1797             :
        1798           0 :      T     = T0 - (latvap/cpair)*(ql0-ql) - ((latvap+latice)/cpair)*(qi0-qi)
        1799           0 :      qv    = qv0 + ql0 - ql + qi0 - qi
        1800             :
        1801             :    ! -------------- !
        1802             :    ! Send to output !
        1803             :    ! -------------- !
        1804             :
        1805           0 :    T_out(i)  = T
        1806           0 :    qv_out(i) = qv
        1807           0 :    ql_out(i) = ql
        1808           0 :    qi_out(i) = qi
        1809           0 :    al_st_out(i) = al_st
        1810           0 :    ai_st_out(i) = ai_st
        1811           0 :    ql_st_out(i) = ql_st
        1812           0 :    qi_st_out(i) = qi_st
        1813             :
        1814             :    enddo
        1815             :
        1816           0 :    return
        1817             :    end subroutine instratus_condensate
        1818             :
        1819             :    ! ----------------- !
        1820             :    ! End of subroutine !
        1821             :    ! ----------------- !
        1822             :
        1823           0 :    subroutine instratus_core( lchnk, icol, k, p,                      &
        1824             :                               T0, qv0, ql0, qi0,                      &
        1825             :                               a_dc, ql_dc, qi_dc,                     &
        1826             :                               a_sc, ql_sc, qi_sc, ai_st,              &
        1827             :                               qcst_crit, Tmin, Tmax, landfrac, snowh, &
        1828             :                               rhminl, rhminl_adj_land, rhminh,        &
        1829             :                               T, qv, ql, qi )
        1830             :
        1831             :    ! ------------------------------------------------------ !
        1832             :    ! Subroutine to find saturation equilibrium state using  !
        1833             :    ! a Newton iteration method, so that 'qc_st = qcst_crit' !
        1834             :    ! is satisfied.                                          !
        1835             :    ! ------------------------------------------------------ !
        1836             :
        1837             :    integer,  intent(in)  :: lchnk      ! Chunk identifier
        1838             :    integer,  intent(in)  :: icol       ! Number of atmospheric columns
        1839             :    integer,  intent(in)  :: k          ! Layer index
        1840             :
        1841             :    real(r8), intent(in)  :: p          ! Pressure [Pa]
        1842             :    real(r8), intent(in)  :: T0         ! Temperature [K]
        1843             :    real(r8), intent(in)  :: qv0        ! Grid-mean water vapor [kg/kg]
        1844             :    real(r8), intent(in)  :: ql0        ! Grid-mean LWC [kg/kg]
        1845             :    real(r8), intent(in)  :: qi0        ! Grid-mean IWC [kg/kg]
        1846             :
        1847             :    real(r8), intent(in)  :: a_dc       ! Deep cumulus cloud fraction
        1848             :    real(r8), intent(in)  :: ql_dc      ! In-deep cumulus LWC [kg/kg]
        1849             :    real(r8), intent(in)  :: qi_dc      ! In-deep cumulus IWC [kg/kg]
        1850             :    real(r8), intent(in)  :: a_sc       ! Shallow cumulus cloud fraction
        1851             :    real(r8), intent(in)  :: ql_sc      ! In-shallow cumulus LWC [kg/kg]
        1852             :    real(r8), intent(in)  :: qi_sc      ! In-shallow cumulus IWC [kg/kg]
        1853             :
        1854             :    real(r8), intent(in)  :: ai_st      ! Ice stratus fraction (fixed)
        1855             :
        1856             :    real(r8), intent(in)  :: Tmin       ! Minimum temperature system can have [K]
        1857             :    real(r8), intent(in)  :: Tmax       ! Maximum temperature system can have [K]
        1858             :    real(r8), intent(in)  :: qcst_crit  ! Critical in-stratus condensate [kg/kg]
        1859             :    real(r8), intent(in)  :: landfrac   ! Land fraction
        1860             :    real(r8), intent(in)  :: snowh      ! Snow depth (liquid water equivalent)
        1861             :
        1862             :    real(r8), intent(in)  :: rhminl
        1863             :    real(r8), intent(in)  :: rhminl_adj_land
        1864             :    real(r8), intent(in)  :: rhminh
        1865             :
        1866             :    real(r8), intent(out) :: T          ! Temperature [K]
        1867             :    real(r8), intent(out) :: qv         ! Grid-mean water vapor [kg/kg]
        1868             :    real(r8), intent(out) :: ql         ! Grid-mean LWC [kg/kg]
        1869             :    real(r8), intent(out) :: qi         ! Grid-mean IWC [kg/kg]
        1870             :
        1871             :    ! Local variables
        1872             :
        1873             :    integer i                           ! Iteration index
        1874             :
        1875             :    real(r8) muQ0, muQ
        1876             :    real(r8) ql_nc0, qi_nc0, qc_nc0, qc_nc
        1877             :    real(r8) fice0, fice
        1878             :    real(r8) ficeg0, ficeg
        1879             :    real(r8) esat0
        1880             :    real(r8) qsat0
        1881             :    real(r8) dqcncdt, dastdt, dUdt
        1882             :    real(r8) alpha, beta
        1883             :    real(r8) U, U_nc
        1884             :    real(r8) al_st_nc, G_nc
        1885             :    real(r8) al_st
        1886             :
        1887             :    ! Variables for root-finding algorithm
        1888             :
        1889             :    integer j
        1890             :    real(r8)  x1, x2
        1891             :    real(r8)  rtsafe
        1892             :    real(r8)  df, dx, dxold, f, fh, fl, temp, xh, xl
        1893             :    real(r8), parameter :: xacc = 1.e-3_r8
        1894             :
        1895             :    ! ---------------- !
        1896             :    ! Main computation !
        1897             :    ! ---------------- !
        1898             :
        1899           0 :    ql_nc0 = max(0._r8,ql0-a_dc*ql_dc-a_sc*ql_sc)
        1900           0 :    qi_nc0 = max(0._r8,qi0-a_dc*qi_dc-a_sc*qi_sc)
        1901           0 :    qc_nc0 = max(0._r8,ql0+qi0-a_dc*(ql_dc+qi_dc)-a_sc*(ql_sc+qi_sc))
        1902           0 :    fice0  = 0._r8
        1903           0 :    ficeg0 = 0._r8
        1904           0 :    muQ0   = 1._r8
        1905             :
        1906             :    ! ------------ !
        1907             :    ! Root finding !
        1908             :    ! ------------ !
        1909             :
        1910           0 :    x1 = Tmin
        1911           0 :    x2 = Tmax
        1912             :    call funcd_instratus( x1, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
        1913             :                          a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
        1914             :                          qcst_crit, landfrac, snowh,                    &
        1915             :                          rhminl, rhminl_adj_land, rhminh,               &
        1916           0 :                          fl, df, qc_nc, fice, al_st )
        1917             :    call funcd_instratus( x2, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
        1918             :                          a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
        1919             :                          qcst_crit, landfrac, snowh,                    &
        1920             :                          rhminl, rhminl_adj_land, rhminh,               &
        1921           0 :                          fh, df, qc_nc, fice, al_st )
        1922           0 :    if((fl > 0._r8 .and. fh > 0._r8) .or. (fl < 0._r8 .and. fh < 0._r8)) then
        1923             :        call funcd_instratus( T0, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
        1924             :                              a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
        1925             :                              qcst_crit, landfrac, snowh,                    &
        1926             :                              rhminl, rhminl_adj_land, rhminh,               &
        1927           0 :                              fl, df, qc_nc, fice, al_st )
        1928           0 :        rtsafe = T0
        1929           0 :        goto 10
        1930             :    endif
        1931           0 :    if( fl == 0._r8) then
        1932           0 :            rtsafe = x1
        1933           0 :            goto 10
        1934           0 :    elseif( fh == 0._r8) then
        1935           0 :            rtsafe = x2
        1936           0 :            goto 10
        1937           0 :    elseif( fl < 0._r8) then
        1938           0 :            xl = x1
        1939           0 :            xh = x2
        1940             :    else
        1941           0 :            xh = x1
        1942           0 :            xl = x2
        1943             :    end if
        1944           0 :    rtsafe = 0.5_r8*(x1+x2)
        1945           0 :    dxold = abs(x2-x1)
        1946           0 :    dx = dxold
        1947             :    call funcd_instratus( rtsafe, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
        1948             :                          a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st,     &
        1949             :                          qcst_crit, landfrac, snowh,                        &
        1950             :                          rhminl, rhminl_adj_land, rhminh,                   &
        1951           0 :                          f, df, qc_nc, fice, al_st )
        1952           0 :    do j = 1, 20
        1953           0 :       if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f) > 0._r8 .or. abs(2.0_r8*f) > abs(dxold*df) ) then
        1954           0 :            dxold = dx
        1955           0 :            dx = 0.5_r8*(xh-xl)
        1956           0 :            rtsafe = xl + dx
        1957           0 :            if(xl == rtsafe) goto 10
        1958             :       else
        1959           0 :            dxold = dx
        1960           0 :            dx = f/df
        1961           0 :            temp = rtsafe
        1962           0 :            rtsafe = rtsafe - dx
        1963           0 :            if (temp == rtsafe) goto 10
        1964             :       end if
        1965             :     ! if(abs(dx) < xacc) goto 10
        1966             :       call funcd_instratus( rtsafe, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
        1967             :                             a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st,     &
        1968             :                             qcst_crit, landfrac, snowh,                        &
        1969             :                             rhminl, rhminl_adj_land, rhminh,                   &
        1970           0 :                             f, df, qc_nc, fice, al_st )
        1971             :     ! Sep.21.2010. Sungsu modified to enhance convergence and guarantee 'qlst_min <  qlst < qlst_max'.
        1972           0 :       if( qcst_crit < 0.5_r8 * ( qlst_min + qlst_max ) ) then
        1973           0 :           if( ( qc_nc*(1._r8-fice) .gt.          qlst_min*al_st .and. &
        1974             :                 qc_nc*(1._r8-fice) .lt. 1.1_r8 * qlst_min*al_st ) ) goto 10
        1975             :       else
        1976           0 :           if( ( qc_nc*(1._r8-fice) .gt. 0.9_r8 * qlst_max*al_st .and. &
        1977             :                 qc_nc*(1._r8-fice) .lt.          qlst_max*al_st ) ) goto 10
        1978             :       endif
        1979           0 :       if(f < 0._r8) then
        1980           0 :           xl = rtsafe
        1981             :       else
        1982           0 :           xh = rtsafe
        1983             :       endif
        1984             :
        1985             :    enddo
        1986             :
        1987             : 10 continue
        1988             :
        1989             :    ! ------------------------------------------- !
        1990             :    ! Final safety check before sending to output !
        1991             :    ! ------------------------------------------- !
        1992             :
        1993           0 :    qc_nc = max(0._r8,qc_nc)
        1994             :
        1995           0 :    T  = rtsafe
        1996           0 :    ql = qc_nc*(1._r8-fice) + a_dc*ql_dc + a_sc*ql_sc
        1997           0 :    qi = qc_nc*fice + a_dc*qi_dc + a_sc*qi_sc
        1998           0 :    qv = qv0 + ql0 + qi0 - (qc_nc + a_dc*(ql_dc+qi_dc) + a_sc*(ql_sc+qi_sc))
        1999           0 :    qv = max(qv,1.e-12_r8)
        2000             :
        2001           0 :    return
        2002             :    end subroutine instratus_core
        2003             :
        2004             :    ! ----------------- !
        2005             :    ! End of subroutine !
        2006             :    ! ----------------- !
        2007             :
        2008           0 :    subroutine funcd_instratus( T, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0,   &
        2009             :                                a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st,  &
        2010             :                                qcst_crit, landfrac, snowh,                     &
        2011             :                                rhminl, rhminl_adj_land, rhminh,                &
        2012             :                                f, fg, qc_nc, fice, al_st )
        2013             :
        2014             :    ! --------------------------------------------------- !
        2015             :    ! Subroutine to find function value and gradient at T !
        2016             :    ! --------------------------------------------------- !
        2017             :
        2018             :    implicit none
        2019             :
        2020             :    real(r8), intent(in)  :: T          ! Iteration temperature [K]
        2021             :
        2022             :    real(r8), intent(in)  :: p          ! Pressure [Pa]
        2023             :    real(r8), intent(in)  :: T0         ! Initial temperature [K]
        2024             :    real(r8), intent(in)  :: qv0        ! Grid-mean water vapor [kg/kg]
        2025             :    real(r8), intent(in)  :: ql0        ! Grid-mean LWC [kg/kg]
        2026             :    real(r8), intent(in)  :: qi0        ! Grid-mean IWC [kg/kg]
        2027             :    real(r8), intent(in)  :: fice0      !
        2028             :    real(r8), intent(in)  :: muQ0       !
        2029             :    real(r8), intent(in)  :: qc_nc0     !
        2030             :
        2031             :    real(r8), intent(in)  :: a_dc       ! Deep cumulus cloud fraction
        2032             :    real(r8), intent(in)  :: ql_dc      ! In-deep cumulus LWC [kg/kg]
        2033             :    real(r8), intent(in)  :: qi_dc      ! In-deep cumulus IWC [kg/kg]
        2034             :    real(r8), intent(in)  :: a_sc       ! Shallow cumulus cloud fraction
        2035             :    real(r8), intent(in)  :: ql_sc      ! In-shallow cumulus LWC [kg/kg]
        2036             :    real(r8), intent(in)  :: qi_sc      ! In-shallow cumulus IWC [kg/kg]
        2037             :
        2038             :    real(r8), intent(in)  :: ai_st      ! Ice stratus fraction (fixed)
        2039             :
        2040             :    real(r8), intent(in)  :: qcst_crit  ! Critical in-stratus condensate [kg/kg]
        2041             :    real(r8), intent(in)  :: landfrac   ! Land fraction
        2042             :    real(r8), intent(in)  :: snowh      ! Snow depth (liquid water equivalent)
        2043             :
        2044             :    real(r8), intent(in)  :: rhminl
        2045             :    real(r8), intent(in)  :: rhminl_adj_land
        2046             :    real(r8), intent(in)  :: rhminh
        2047             :
        2048             :    real(r8), intent(out) :: f          ! Value of minimization function at T
        2049             :    real(r8), intent(out) :: fg         ! Gradient of minimization function
        2050             :    real(r8), intent(out) :: qc_nc      !
        2051             :    real(r8), intent(out) :: al_st      !
        2052             :    real(r8), intent(out) :: fice       !
        2053             :
        2054             :    ! Local variables
        2055             :
        2056             :    real(r8) es
        2057             :    real(r8) qs
        2058             :    real(r8) dqsdT
        2059             :    real(r8) dqcncdt
        2060             :    real(r8) alpha
        2061             :    real(r8) beta
        2062             :    real(r8) U
        2063             :    real(r8) U_nc
        2064             :    real(r8) al_st_nc
        2065             :    real(r8) G_nc
        2066             :    real(r8) dUdt
        2067             :    real(r8) dalstdt
        2068             :    real(r8) qv
        2069             :
        2070             :    ! ---------------- !
        2071             :    ! Main computation !
        2072             :    ! ---------------- !
        2073             :
        2074           0 :    call qsat_water(T, p, es, qs, dqsdt=dqsdT)
        2075             :
        2076           0 :    fice    = fice0
        2077           0 :    qc_nc   = (cpair/latvap)*(T-T0)+muQ0*qc_nc0
        2078           0 :    dqcncdt = (cpair/latvap)
        2079           0 :    qv      = (qv0 + ql0 + qi0 - (qc_nc + a_dc*(ql_dc+qi_dc) + a_sc*(ql_sc+qi_sc)))
        2080           0 :    alpha   = (1._r8/qs)
        2081           0 :    beta    = (qv/qs**2._r8)*dqsdT
        2082             :
        2083           0 :    U      =  (qv/qs)
        2084           0 :    U_nc   =   U
        2085             :    if( CAMstfrac ) then
        2086             :        call astG_RHU_single(U_nc, p, qv, landfrac, snowh, al_st_nc, G_nc, &
        2087             :           rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
        2088             :    else
        2089             :        call astG_PDF_single(U_nc, p, qv, landfrac, snowh, al_st_nc, G_nc, &
        2090           0 :           rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
        2091             :    endif
        2092           0 :    al_st   =  (1._r8-a_dc-a_sc)*al_st_nc
        2093           0 :    dUdt    = -(alpha*dqcncdt+beta)
        2094           0 :    dalstdt =  (1._r8/G_nc)*dUdt
        2095           0 :    if( U_nc .eq. 1._r8 ) dalstdt = 0._r8
        2096             :
        2097           0 :    f  = qc_nc   - qcst_crit*al_st
        2098           0 :    fg = dqcncdt - qcst_crit*dalstdt
        2099             :
        2100           0 :    return
        2101             :    end subroutine funcd_instratus
        2102             :
        2103             :    ! ----------------- !
        2104             :    ! End of subroutine !
        2105             :    ! ----------------- !
        2106             :
        2107           0 :    subroutine gridmean_RH( lchnk, icol, k, p, T, qv, ql, qi,       &
        2108             :                            a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, &
        2109             :                            landfrac, snowh )
        2110             :
        2111             :    ! ------------------------------------------------------------- !
        2112             :    ! Subroutine to force grid-mean RH = 1 when RH > 1              !
        2113             :    ! This is condensation process similar to instratus_condensate. !
        2114             :    ! During condensation, we assume 'fice' is maintained in this   !
        2115             :    ! verison for MG not for RK.                                    !
        2116             :    ! ------------------------------------------------------------- !
        2117             :
        2118             :    integer,  intent(in)    :: lchnk      ! Chunk identifier
        2119             :    integer,  intent(in)    :: icol       ! Number of atmospheric columns
        2120             :    integer,  intent(in)    :: k          ! Layer index
        2121             :
        2122             :    real(r8), intent(in)    :: p          ! Pressure [Pa]
        2123             :    real(r8), intent(inout) :: T          ! Temperature [K]
        2124             :    real(r8), intent(inout) :: qv         ! Grid-mean water vapor [kg/kg]
        2125             :    real(r8), intent(inout) :: ql         ! Grid-mean LWC [kg/kg]
        2126             :    real(r8), intent(inout) :: qi         ! Grid-mean IWC [kg/kg]
        2127             :
        2128             :    real(r8), intent(in)    :: a_dc       ! Deep cumulus cloud fraction
        2129             :    real(r8), intent(in)    :: ql_dc      ! In-deep cumulus LWC [kg/kg]
        2130             :    real(r8), intent(in)    :: qi_dc      ! In-deep cumulus IWC [kg/kg]
        2131             :    real(r8), intent(in)    :: a_sc       ! Shallow cumulus cloud fraction
        2132             :    real(r8), intent(in)    :: ql_sc      ! In-shallow cumulus LWC [kg/kg]
        2133             :    real(r8), intent(in)    :: qi_sc      ! In-shallow cumulus IWC [kg/kg]
        2134             :
        2135             :    real(r8), intent(in)    :: landfrac   ! Land fraction
        2136             :    real(r8), intent(in)    :: snowh      ! Snow depth (liquid water equivalent)
        2137             :
        2138             :    ! Local variables
        2139             :
        2140             :    integer m                             ! Iteration index
        2141             :
        2142             :    real(r8)  ql_nc0, qi_nc0, qc_nc0
        2143             :    real(r8)  Tscale
        2144             :    real(r8)  Tc, qt, qc, dqcdt, qc_nc
        2145             :    real(r8)  es, qs, dqsdT
        2146             :    real(r8)  al_st_nc, G_nc
        2147             :    real(r8)  f, fg
        2148             :    real(r8), parameter :: xacc = 1.e-3_r8
        2149             :
        2150             :    ! ---------------- !
        2151             :    ! Main computation !
        2152             :    ! ---------------- !
        2153             :
        2154           0 :    ql_nc0 = max(0._r8,ql-a_dc*ql_dc-a_sc*ql_sc)
        2155           0 :    qi_nc0 = max(0._r8,qi-a_dc*qi_dc-a_sc*qi_sc)
        2156           0 :    qc_nc0 = max(0._r8,ql+qi-a_dc*(ql_dc+qi_dc)-a_sc*(ql_sc+qi_sc))
        2157           0 :    Tc    = T - (latvap/cpair)*ql
        2158           0 :    qt    = qv + ql
        2159             :
        2160           0 :    do m = 1, 20
        2161           0 :       call qsat_water(T, p, es, qs, dqsdt=dqsdT)
        2162           0 :       Tscale = latvap/cpair
        2163           0 :       qc     = (T-Tc)/Tscale
        2164           0 :       dqcdt  = 1._r8/Tscale
        2165           0 :       f      = qs + qc - qt
        2166           0 :       fg     = dqsdT + dqcdt
        2167           0 :       fg     = sign(1._r8,fg)*max(1.e-10_r8,abs(fg))
        2168             :     ! Sungsu modified convergence criteria to speed up convergence and guarantee RH <= 1.
        2169           0 :       if( qc .ge. 0._r8 .and. ( qt - qc ) .ge. 0.999_r8*qs .and. ( qt - qc ) .le. 1._r8*qs ) then
        2170             :           goto 10
        2171             :       endif
        2172           0 :       T = T - f/fg
        2173             :    enddo
        2174             :  ! write(iulog,*) 'Convergence in gridmean_RH is not reached. RH = ', ( qt - qc ) / qs
        2175             : 10 continue
        2176             :
        2177           0 :    call qsat_water(T, p, es, qs)
        2178             :  ! Sungsu modified 'qv = qs' in consistent with the modified convergence criteria above.
        2179           0 :    qv = min(qt,qs) ! Modified
        2180           0 :    ql = qt - qv
        2181           0 :    T  = Tc + (latvap/cpair)*ql
        2182             :
        2183           0 :    return
        2184             :    end subroutine gridmean_RH
        2185             :
        2186             :    ! ----------------- !
        2187             :    ! End of subroutine !
        2188             :    ! ----------------- !
        2189             :
        2190           0 :    subroutine positive_moisture( ncol, dt, qvmin, qlmin, qimin, dp, &
        2191             :                                  qv,   ql, qi,    t,     qvten, &
        2192             :                                  qlten,    qiten, tten,  do_cldice)
        2193             :
        2194             :    ! ------------------------------------------------------------------------------- !
        2195             :    ! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer,         !
        2196             :    ! force them to be larger than minimum value by (1) condensating water vapor      !
        2197             :    ! into liquid or ice, and (2) by transporting water vapor from the very lower     !
        2198             :    ! layer. '2._r8' is multiplied to the minimum values for safety.                  !
        2199             :    ! Update final state variables and tendencies associated with this correction.    !
        2200             :    ! If any condensation happens, update (s,t) too.                                  !
        2201             :    ! Note that (qv,ql,qi,t,s) are final state variables after applying corresponding !
        2202             :    ! input tendencies.                                                               !
        2203             :    ! Be careful the order of k : '1': top layer, 'pver' : near-surface layer         !
        2204             :    ! ------------------------------------------------------------------------------- !
        2205             :
        2206             :    implicit none
        2207             :    integer,  intent(in)     :: ncol
        2208             :    real(r8), intent(in)     :: dt
        2209             :    real(r8), intent(in)     :: dp(pcols,pver), qvmin(pcols,pver), qlmin(pcols,pver), qimin(pcols,pver)
        2210             :    real(r8), intent(inout)  :: qv(pcols,pver), ql(pcols,pver), qi(pcols,pver), t(pcols,pver)
        2211             :    real(r8), intent(out)    :: qvten(pcols,pver), qlten(pcols,pver), qiten(pcols,pver), tten(pcols,pver)
        2212             :    logical, intent(in)      :: do_cldice
        2213             :    integer   i, k
        2214             :    real(r8)  dql, dqi, dqv, sum, aa, dum
        2215             :
        2216           0 :    tten(:ncol,:pver)  = 0._r8
        2217           0 :    qvten(:ncol,:pver) = 0._r8
        2218           0 :    qlten(:ncol,:pver) = 0._r8
        2219           0 :    qiten(:ncol,:pver) = 0._r8
        2220             :
        2221           0 :    do i = 1, ncol
        2222           0 :       do k = top_lev, pver
        2223           0 :          if( qv(i,k) .lt. qvmin(i,k) .or. ql(i,k) .lt. qlmin(i,k) .or. qi(i,k) .lt. qimin(i,k) ) then
        2224             :              goto 10
        2225             :          endif
        2226             :       enddo
        2227           0 :       goto 11
        2228             :    10 continue
        2229           0 :       do k = top_lev, pver    ! From the top to the 1st (lowest) layer from the surface
        2230           0 :          dql = max(0._r8,1._r8*qlmin(i,k)-ql(i,k))
        2231             :
        2232           0 :          if (do_cldice) then
        2233           0 :          dqi = max(0._r8,1._r8*qimin(i,k)-qi(i,k))
        2234             :          else
        2235             :            dqi = 0._r8
        2236             :          end if
        2237             :
        2238           0 :          qlten(i,k) = qlten(i,k) +  dql/dt
        2239           0 :          qiten(i,k) = qiten(i,k) +  dqi/dt
        2240           0 :          qvten(i,k) = qvten(i,k) - (dql+dqi)/dt
        2241           0 :          tten(i,k)  = tten(i,k)  + (latvap/cpair)*(dql/dt) + ((latvap+latice)/cpair)*(dqi/dt)
        2242           0 :          ql(i,k)    = ql(i,k) + dql
        2243           0 :          qi(i,k)    = qi(i,k) + dqi
        2244           0 :          qv(i,k)    = qv(i,k) - dql - dqi
        2245           0 :          t(i,k)     = t(i,k)  + (latvap * dql + (latvap+latice) * dqi)/cpair
        2246           0 :          dqv        = max(0._r8,1._r8*qvmin(i,k)-qv(i,k))
        2247           0 :          qvten(i,k) = qvten(i,k) + dqv/dt
        2248           0 :          qv(i,k)    = qv(i,k)    + dqv
        2249           0 :          if( k .ne. pver ) then
        2250           0 :              qv(i,k+1)    = qv(i,k+1)    - dqv*dp(i,k)/dp(i,k+1)
        2251           0 :              qvten(i,k+1) = qvten(i,k+1) - dqv*dp(i,k)/dp(i,k+1)/dt
        2252             :          endif
        2253           0 :          qv(i,k) = max(qv(i,k),qvmin(i,k))
        2254           0 :          ql(i,k) = max(ql(i,k),qlmin(i,k))
        2255           0 :          qi(i,k) = max(qi(i,k),qimin(i,k))
        2256             :       end do
        2257             :       ! Extra moisture used to satisfy 'qv(i,pver)=qvmin' is proportionally
        2258             :       ! extracted from all the layers that has 'qv > 2*qvmin'. This fully
        2259             :       ! preserves column moisture.
        2260           0 :       if( dqv .gt. 1.e-20_r8 ) then
        2261             :           sum = 0._r8
        2262           0 :           do k = top_lev, pver
        2263           0 :              if( qv(i,k) .gt. 2._r8*qvmin(i,k) ) sum = sum + qv(i,k)*dp(i,k)
        2264             :           enddo
        2265           0 :           aa = dqv*dp(i,pver)/max(1.e-20_r8,sum)
        2266           0 :           if( aa .lt. 0.5_r8 ) then
        2267           0 :               do k = top_lev, pver
        2268           0 :                  if( qv(i,k) .gt. 2._r8*qvmin(i,k) ) then
        2269           0 :                      dum        = aa*qv(i,k)
        2270           0 :                      qv(i,k)    = qv(i,k) - dum
        2271           0 :                      qvten(i,k) = qvten(i,k) - dum/dt
        2272             :                  endif
        2273             :               enddo
        2274             :           else
        2275           0 :               write(iulog,*) 'Full positive_moisture is impossible in Park Macro'
        2276             :           endif
        2277             :       endif
        2278           0 : 11 continue
        2279             :    enddo
        2280           0 :    return
        2281             :
        2282             :    end subroutine positive_moisture
        2283             :
        2284             :    ! ----------------- !
        2285             :    ! End of subroutine !
        2286             :    ! ----------------- !
        2287             :
        2288           0 :       SUBROUTINE gaussj(a,n,np,b,m,mp)
        2289             :       INTEGER m,mp,n,np,NMAX
        2290             :       real(r8) a(np,np),b(np,mp)
        2291           0 :       real(r8) aa(np,np),bb(np,mp)
        2292             :       PARAMETER (NMAX=50)
        2293             :       INTEGER i,icol,irow,j,k,l,ll,ii,jj,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
        2294             :       real(r8) big,dum,pivinv
        2295             :
        2296           0 :       aa(:,:) = a(:,:)
        2297           0 :       bb(:,:) = b(:,:)
        2298             :
        2299           0 :       do 11 j=1,n
        2300           0 :         ipiv(j)=0
        2301           0 : 11    continue
        2302           0 :       do 22 i=1,n
        2303           0 :         big=0._r8
        2304           0 :         do 13 j=1,n
        2305           0 :           if(ipiv(j).ne.1)then
        2306           0 :             do 12 k=1,n
        2307           0 :               if (ipiv(k).eq.0) then
        2308           0 :                 if (abs(a(j,k)).ge.big)then
        2309           0 :                   big=abs(a(j,k))
        2310           0 :                   irow=j
        2311           0 :                   icol=k
        2312             :                 endif
        2313           0 :               else if (ipiv(k).gt.1) then
        2314           0 :                 write(iulog,*) 'singular matrix in gaussj 1'
        2315           0 :                 do ii = 1, np
        2316           0 :                 do jj = 1, np
        2317           0 :                    write(iulog,*) ii, jj, aa(ii,jj), bb(ii,1)
        2318             :                 end do
        2319             :                 end do
        2320           0 :                 call endrun
        2321             :               endif
        2322           0 : 12          continue
        2323             :           endif
        2324           0 : 13      continue
        2325           0 :         ipiv(icol)=ipiv(icol)+1
        2326           0 :         if (irow.ne.icol) then
        2327           0 :           do 14 l=1,n
        2328           0 :             dum=a(irow,l)
        2329           0 :             a(irow,l)=a(icol,l)
        2330           0 :             a(icol,l)=dum
        2331           0 : 14        continue
        2332           0 :           do 15 l=1,m
        2333           0 :             dum=b(irow,l)
        2334           0 :             b(irow,l)=b(icol,l)
        2335           0 :             b(icol,l)=dum
        2336           0 : 15        continue
        2337             :         endif
        2338           0 :         indxr(i)=irow
        2339           0 :         indxc(i)=icol
        2340           0 :         if (a(icol,icol).eq.0._r8) then
        2341           0 :             write(iulog,*) 'singular matrix in gaussj 2'
        2342           0 :             do ii = 1, np
        2343           0 :             do jj = 1, np
        2344           0 :                write(iulog,*) ii, jj, aa(ii,jj), bb(ii,1)
        2345             :             end do
        2346             :             end do
        2347           0 :             call endrun
        2348             :         endif
        2349           0 :         pivinv=1._r8/a(icol,icol)
        2350           0 :         a(icol,icol)=1._r8
        2351           0 :         do 16 l=1,n
        2352           0 :           a(icol,l)=a(icol,l)*pivinv
        2353           0 : 16      continue
        2354           0 :         do 17 l=1,m
        2355           0 :           b(icol,l)=b(icol,l)*pivinv
        2356           0 : 17      continue
        2357           0 :         do 21 ll=1,n
        2358           0 :           if(ll.ne.icol)then
        2359           0 :             dum=a(ll,icol)
        2360           0 :             a(ll,icol)=0._r8
        2361           0 :             do 18 l=1,n
        2362           0 :               a(ll,l)=a(ll,l)-a(icol,l)*dum
        2363           0 : 18          continue
        2364           0 :             do 19 l=1,m
        2365           0 :               b(ll,l)=b(ll,l)-b(icol,l)*dum
        2366           0 : 19          continue
        2367             :           endif
        2368           0 : 21      continue
        2369           0 : 22    continue
        2370           0 :       do 24 l=n,1,-1
        2371           0 :         if(indxr(l).ne.indxc(l))then
        2372           0 :           do 23 k=1,n
        2373           0 :             dum=a(k,indxr(l))
        2374           0 :             a(k,indxr(l))=a(k,indxc(l))
        2375           0 :             a(k,indxc(l))=dum
        2376           0 : 23        continue
        2377             :         endif
        2378           0 : 24    continue
        2379             :
        2380           0 :       return
        2381             :       end subroutine gaussj
        2382             :
        2383             :    ! ----------------- !
        2384             :    ! End of subroutine !
        2385             :    ! ----------------- !
        2386             :
        2387             : end module cldwat2m_macro

Generated by: LCOV version 1.14