LCOV - code coverage report
Current view: | top level - atmos_phys/schemes/zhang_mcfarlane - zm_convr.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 992 | 1027 | 96.6 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 10 | 10 | 100.0 %

              Line data    Source code
           1             : module zm_convr
           2             :
           3             :   use ccpp_kinds, only:  kind_phys
           4             :
           5             :   implicit none
           6             :
           7             :   save
           8             :   private                         ! Make default type private to the module
           9             : !
          10             : ! PUBLIC: interfaces
          11             : !
          12             :   public zm_convr_init                 ! ZM schemea
          13             :   public zm_convr_run                  ! ZM schemea
          14             :
          15             :    real(kind_phys) rl         ! wg latent heat of vaporization.
          16             :    real(kind_phys) cpres      ! specific heat at constant pressure in j/kg-degk.
          17             :    real(kind_phys) :: capelmt ! namelist configurable:
          18             :                        ! threshold value for cape for deep convection.
          19             :    real(kind_phys) :: ke           ! Tunable evaporation efficiency set from namelist input zmconv_ke
          20             :    real(kind_phys) :: ke_lnd
          21             :    real(kind_phys) :: c0_lnd       ! set from namelist input zmconv_c0_lnd
          22             :    real(kind_phys) :: c0_ocn       ! set from namelist input zmconv_c0_ocn
          23             :    integer  :: num_cin      ! set from namelist input zmconv_num_cin
          24             :                             ! The number of negative buoyancy regions that are allowed
          25             :                             ! before the convection top and CAPE calculations are completed.
          26             :    logical  :: zm_org
          27             :    real(kind_phys) tau   ! convective time scale
          28             :    real(kind_phys) :: tfreez
          29             :    real(kind_phys) :: eps1
          30             :    real(kind_phys) :: momcu
          31             :    real(kind_phys) :: momcd
          32             :
          33             :    logical :: no_deep_pbl ! default = .false.
          34             :                           ! no_deep_pbl = .true. eliminates deep convection entirely within PBL
          35             :
          36             :
          37             :    real(kind_phys) :: rgrav       ! reciprocal of grav
          38             :    real(kind_phys) :: rgas        ! gas constant for dry air
          39             :    real(kind_phys) :: grav        ! = gravit
          40             :    real(kind_phys) :: cp          ! = cpres = cpair
          41             :
          42             :    integer  limcnv       ! top interface level limit for convection
          43             :
          44             :    logical :: lparcel_pbl     ! Switch to turn on mixing of parcel MSE air, and picking launch level to be the top of the PBL.
          45             :
          46             :
          47             :    real(kind_phys) :: tiedke_add      ! namelist configurable
          48             :    real(kind_phys) :: dmpdz_param     ! namelist configurable
          49             :
          50             : contains
          51             :
          52             :
          53             :
          54             : !===============================================================================
          55             : !> \section arg_table_zm_convr_init Argument Table
          56             : !! \htmlinclude zm_convr_init.html
          57             : !!
          58        1536 : subroutine zm_convr_init(cpair, epsilo, gravit, latvap, tmelt, rair, &
          59             :                     limcnv_in, zmconv_c0_lnd, zmconv_c0_ocn, zmconv_ke, zmconv_ke_lnd, &
          60             :                     zmconv_momcu, zmconv_momcd, zmconv_num_cin, zmconv_org, &
          61             :                     no_deep_pbl_in, zmconv_tiedke_add, &
          62             :                     zmconv_capelmt, zmconv_dmpdz, zmconv_parcel_pbl, zmconv_tau, &
          63        1536 :                     masterproc, iulog, errmsg, errflg)
          64             :
          65             :    real(kind_phys), intent(in)   :: cpair           ! specific heat of dry air (J K-1 kg-1)
          66             :    real(kind_phys), intent(in)   :: epsilo          ! ratio of h2o to dry air molecular weights
          67             :    real(kind_phys), intent(in)   :: gravit          ! gravitational acceleration (m s-2)
          68             :    real(kind_phys), intent(in)   :: latvap          ! Latent heat of vaporization (J kg-1)
          69             :    real(kind_phys), intent(in)   :: tmelt           ! Freezing point of water (K)
          70             :    real(kind_phys), intent(in)   :: rair            ! Dry air gas constant     (J K-1 kg-1)
          71             :    integer, intent(in)           :: limcnv_in       ! top interface level limit for convection
          72             :    integer, intent(in)           :: zmconv_num_cin  ! Number negative buoyancy regions that are allowed
          73             :                                                     ! before the convection top and CAPE calculations are completed.
          74             :    real(kind_phys),intent(in)           :: zmconv_c0_lnd
          75             :    real(kind_phys),intent(in)           :: zmconv_c0_ocn
          76             :    real(kind_phys),intent(in)           :: zmconv_ke
          77             :    real(kind_phys),intent(in)           :: zmconv_ke_lnd
          78             :    real(kind_phys),intent(in)           :: zmconv_momcu
          79             :    real(kind_phys),intent(in)           :: zmconv_momcd
          80             :    logical, intent(in)           :: zmconv_org
          81             :    logical, intent(in)           :: no_deep_pbl_in  ! no_deep_pbl = .true. eliminates ZM convection entirely within PBL
          82             :    real(kind_phys),intent(in)           :: zmconv_tiedke_add
          83             :    real(kind_phys),intent(in)           :: zmconv_capelmt
          84             :    real(kind_phys),intent(in)           :: zmconv_dmpdz
          85             :    logical, intent(in)           :: zmconv_parcel_pbl ! Should the parcel properties include PBL mixing?
          86             :    real(kind_phys),intent(in)           :: zmconv_tau
          87             :    logical, intent(in)                  :: masterproc
          88             :    integer, intent(in)                  :: iulog
          89             :    character(len=512), intent(out)      :: errmsg
          90             :    integer, intent(out)                 :: errflg
          91             :
          92        1536 :    errmsg =''
          93        1536 :    errflg = 0
          94             :
          95             :    ! Initialization of ZM constants
          96        1536 :    limcnv = limcnv_in
          97        1536 :    tfreez = tmelt
          98        1536 :    eps1   = epsilo
          99        1536 :    rl     = latvap
         100        1536 :    cpres  = cpair
         101        1536 :    rgrav  = 1.0_kind_phys/gravit
         102        1536 :    rgas   = rair
         103        1536 :    grav   = gravit
         104        1536 :    cp     = cpres
         105             :
         106        1536 :    c0_lnd  = zmconv_c0_lnd
         107        1536 :    c0_ocn  = zmconv_c0_ocn
         108        1536 :    num_cin = zmconv_num_cin
         109        1536 :    ke      = zmconv_ke
         110        1536 :    ke_lnd  = zmconv_ke_lnd
         111        1536 :    zm_org  = zmconv_org
         112        1536 :    momcu   = zmconv_momcu
         113        1536 :    momcd   = zmconv_momcd
         114             :
         115        1536 :    tiedke_add = zmconv_tiedke_add
         116        1536 :    capelmt = zmconv_capelmt
         117        1536 :    dmpdz_param = zmconv_dmpdz
         118        1536 :    no_deep_pbl = no_deep_pbl_in
         119        1536 :    lparcel_pbl = zmconv_parcel_pbl
         120             :
         121        1536 :    tau = zmconv_tau
         122             :
         123        1536 :    if ( masterproc ) then
         124           2 :       write(iulog,*) 'tuning parameters zm_convr_init: tau',tau
         125           2 :       write(iulog,*) 'tuning parameters zm_convr_init: c0_lnd',c0_lnd, ', c0_ocn', c0_ocn
         126           2 :       write(iulog,*) 'tuning parameters zm_convr_init: num_cin', num_cin
         127           2 :       write(iulog,*) 'tuning parameters zm_convr_init: ke',ke
         128           2 :       write(iulog,*) 'tuning parameters zm_convr_init: no_deep_pbl',no_deep_pbl
         129           2 :       write(iulog,*) 'tuning parameters zm_convr_init: zm_capelmt', capelmt
         130           2 :       write(iulog,*) 'tuning parameters zm_convr_init: zm_dmpdz', dmpdz_param
         131           2 :       write(iulog,*) 'tuning parameters zm_convr_init: zm_tiedke_add', tiedke_add
         132           2 :       write(iulog,*) 'tuning parameters zm_convr_init: zm_parcel_pbl', lparcel_pbl
         133             :    endif
         134             :
         135        1536 :    if (masterproc) write(iulog,*)'**** ZM: DILUTE Buoyancy Calculation ****'
         136             :
         137        1536 : end subroutine zm_convr_init
         138             :
         139             :
         140             : !===============================================================================
         141             : !> \section arg_table_zm_convr_run Argument Table
         142             : !! \htmlinclude zm_convr_run.html
         143             : !!
         144     1495368 : subroutine zm_convr_run(     ncol    ,pver    , &
         145             :                     pverp,   gravit  ,latice  ,cpwv    ,cpliq   , rh2o, &
         146     1495368 :                     t       ,qh      ,prec    , &
         147     1495368 :                     pblh    ,zm      ,geos    ,zi      ,qtnd    , &
         148     2990736 :                     heat    ,pap     ,paph    ,dpp     , &
         149     2990736 :                     delt    ,mcon    ,cme     ,cape    , &
         150     2990736 :                     tpert   ,dlf     ,zdu     ,rprd    , &
         151     1495368 :                     mu      ,md      ,du      ,eu      ,ed      , &
         152     2990736 :                     dp      ,dsubcld ,jt      ,maxg    ,ideep   , &
         153     4486104 :                     ql      ,rliq    ,landfrac,                   &
         154     2990736 :                     org     ,orgt    ,org2d   ,  &
         155     1495368 :                     dif     ,dnlf    ,dnif    , &
         156     1495368 :                     rice   ,errmsg  ,errflg)
         157             : !-----------------------------------------------------------------------
         158             : !
         159             : ! Purpose:
         160             : ! Main driver for zhang-mcfarlane convection scheme
         161             : !
         162             : ! Method:
         163             : ! performs deep convective adjustment based on mass-flux closure
         164             : ! algorithm.
         165             : !
         166             : ! Author:guang jun zhang, m.lazare, n.mcfarlane. CAM Contact: P. Rasch
         167             : !
         168             : ! This is contributed code not fully standardized by the CAM core group.
         169             : ! All variables have been typed, where most are identified in comments
         170             : ! The current procedure will be reimplemented in a subsequent version
         171             : ! of the CAM where it will include a more straightforward formulation
         172             : ! and will make use of the standard CAM nomenclature
         173             : !
         174             : !-----------------------------------------------------------------------
         175             : !
         176             : ! ************************ index of variables **********************
         177             : !
         178             : !  wg * alpha    array of vertical differencing used (=1. for upstream).
         179             : !  w  * cape     convective available potential energy.
         180             : !  wg * capeg    gathered convective available potential energy.
         181             : !  c  * capelmt  threshold value for cape for deep convection.
         182             : !  ic  * cpres    specific heat at constant pressure in j/kg-degk.
         183             : !  i  * dpp
         184             : !  ic  * delt     length of model time-step in seconds.
         185             : !  wg * dp       layer thickness in mbs (between upper/lower interface).
         186             : !  wg * dqdt     mixing ratio tendency at gathered points.
         187             : !  wg * dsdt     dry static energy ("temp") tendency at gathered points.
         188             : !  wg * dudt     u-wind tendency at gathered points.
         189             : !  wg * dvdt     v-wind tendency at gathered points.
         190             : !  wg * dsubcld  layer thickness in mbs between lcl and maxi.
         191             : !  ic  * grav     acceleration due to gravity in m/sec2.
         192             : !  wg * du       detrainment in updraft. specified in mid-layer
         193             : !  wg * ed       entrainment in downdraft.
         194             : !  wg * eu       entrainment in updraft.
         195             : !  wg * hmn      moist static energy.
         196             : !  wg * hsat     saturated moist static energy.
         197             : !  w  * ideep    holds position of gathered points vs longitude index.
         198             : !  ic  * pver     number of model levels.
         199             : !  wg * j0       detrainment initiation level index.
         200             : !  wg * jd       downdraft   initiation level index.
         201             : !  ic  * jlatpr   gaussian latitude index for printing grids (if needed).
         202             : !  wg * jt       top  level index of deep cumulus convection.
         203             : !  w  * lcl      base level index of deep cumulus convection.
         204             : !  wg * lclg     gathered values of lcl.
         205             : !  w  * lel      index of highest theoretical convective plume.
         206             : !  wg * lelg     gathered values of lel.
         207             : !  w  * lon      index of onset level for deep convection.
         208             : !  w  * maxi     index of level with largest moist static energy.
         209             : !  wg * maxg     gathered values of maxi.
         210             : !  wg * mb       cloud base mass flux.
         211             : !  wg * mc       net upward (scaled by mb) cloud mass flux.
         212             : !  wg * md       downward cloud mass flux (positive up).
         213             : !  wg * mu       upward   cloud mass flux (positive up). specified
         214             : !                at interface
         215             : !  ic  * msg      number of missing moisture levels at the top of model.
         216             : !  w  * p        grid slice of ambient mid-layer pressure in mbs.
         217             : !  i  * pblt     row of pbl top indices.
         218             : !  w  * pcpdh    scaled surface pressure.
         219             : !  w  * pf       grid slice of ambient interface pressure in mbs.
         220             : !  wg * pg       grid slice of gathered values of p.
         221             : !  w  * q        grid slice of mixing ratio.
         222             : !  wg * qd       grid slice of mixing ratio in downdraft.
         223             : !  wg * qg       grid slice of gathered values of q.
         224             : !  i/o * qh       grid slice of specific humidity.
         225             : !  w  * qh0      grid slice of initial specific humidity.
         226             : !  wg * qhat     grid slice of upper interface mixing ratio.
         227             : !  wg * ql       grid slice of cloud liquid water.
         228             : !  wg * qs       grid slice of saturation mixing ratio.
         229             : !  w  * qstp     grid slice of parcel temp. saturation mixing ratio.
         230             : !  wg * qstpg    grid slice of gathered values of qstp.
         231             : !  wg * qu       grid slice of mixing ratio in updraft.
         232             : !  ic  * rgas     dry air gas constant.
         233             : !  wg * rl       latent heat of vaporization.
         234             : !  w  * s        grid slice of scaled dry static energy (t+gz/cp).
         235             : !  wg * sd       grid slice of dry static energy in downdraft.
         236             : !  wg * sg       grid slice of gathered values of s.
         237             : !  wg * shat     grid slice of upper interface dry static energy.
         238             : !  wg * su       grid slice of dry static energy in updraft.
         239             : !  i/o * t
         240             : !  wg * tg       grid slice of gathered values of t.
         241             : !  w  * tl       row of parcel temperature at lcl.
         242             : !  wg * tlg      grid slice of gathered values of tl.
         243             : !  w  * tp       grid slice of parcel temperatures.
         244             : !  wg * tpg      grid slice of gathered values of tp.
         245             : !  i/o * u        grid slice of u-wind (real).
         246             : !  wg * ug       grid slice of gathered values of u.
         247             : !  i/o * utg      grid slice of u-wind tendency (real).
         248             : !  i/o * v        grid slice of v-wind (real).
         249             : !  w  * va       work array re-used by called subroutines.
         250             : !  wg * vg       grid slice of gathered values of v.
         251             : !  i/o * vtg      grid slice of v-wind tendency (real).
         252             : !  i  * w        grid slice of diagnosed large-scale vertical velocity.
         253             : !  w  * z        grid slice of ambient mid-layer height in metres.
         254             : !  w  * zf       grid slice of ambient interface height in metres.
         255             : !  wg * zfg      grid slice of gathered values of zf.
         256             : !  wg * zg       grid slice of gathered values of z.
         257             : !
         258             : !-----------------------------------------------------------------------
         259             : !
         260             : ! multi-level i/o fields:
         261             : !  i      => input arrays.
         262             : !  i/o    => input/output arrays.
         263             : !  w      => work arrays.
         264             : !  wg     => work arrays operating only on gathered points.
         265             : !  ic     => input data constants.
         266             : !  c      => data constants pertaining to subroutine itself.
         267             : !
         268             : ! input arguments
         269             : !
         270             :    integer, intent(in) :: ncol                    ! number of atmospheric columns
         271             :    integer, intent(in) :: pver, pverp
         272             :
         273             :    real(kind_phys), intent(in) :: gravit          ! gravitational acceleration (m s-2)
         274             :    real(kind_phys), intent(in) :: latice          ! Latent heat of fusion (J kg-1)
         275             :    real(kind_phys), intent(in) :: cpwv            ! specific heat of water vapor (J K-1 kg-1)
         276             :    real(kind_phys), intent(in) :: cpliq           ! specific heat of fresh h2o (J K-1 kg-1)
         277             :    real(kind_phys), intent(in) :: rh2o            ! Water vapor gas constant (J K-1 kg-1)
         278             :
         279             :    real(kind_phys), intent(in) :: t(:,:)          ! grid slice of temperature at mid-layer.           (ncol,pver)
         280             :    real(kind_phys), intent(in) :: qh(:,:)         ! grid slice of specific humidity.                  (ncol,pver)
         281             :    real(kind_phys), intent(in) :: pap(:,:)        !                                                   (ncol,pver)
         282             :    real(kind_phys), intent(in) :: paph(:,:)     !                                                     (ncol,pver+1)
         283             :    real(kind_phys), intent(in) :: dpp(:,:)        ! local sigma half-level thickness (i.e. dshj).     (ncol,pver)
         284             :    real(kind_phys), intent(in) :: zm(:,:)       !                                                     (ncol,pver)
         285             :    real(kind_phys), intent(in) :: geos(:)       !                                                     (ncol)
         286             :    real(kind_phys), intent(in) :: zi(:,:)      !                                                      (ncol,pver+1)
         287             :    real(kind_phys), intent(in) :: pblh(:)     !                                                       (ncol)
         288             :    real(kind_phys), intent(in) :: tpert(:)      !                                                     (ncol)
         289             :    real(kind_phys), intent(in) :: landfrac(:) ! RBN Landfrac                                          (ncol)
         290             :
         291             : ! output arguments
         292             : !
         293             :    real(kind_phys), intent(out) :: qtnd(:,:)           ! specific humidity tendency (kg/kg/s)             (ncol,pver)
         294             :    real(kind_phys), intent(out) :: heat(:,:)           ! heating rate (dry static energy tendency, W/kg)  (ncol,pver)
         295             :    real(kind_phys), intent(out) :: mcon(:,:)  !   (ncol,pverp)
         296             :    real(kind_phys), intent(out) :: dlf(:,:)    ! scattrd version of the detraining cld h2o tend (ncol,pver)
         297             :    real(kind_phys), intent(out) :: cme(:,:)    !                                                          (ncol,pver)
         298             :    real(kind_phys), intent(out) :: cape(:)        ! w  convective available potential energy.             (ncol)
         299             :    real(kind_phys), intent(out) :: zdu(:,:)    ! (ncol,pver)
         300             :    real(kind_phys), intent(out) :: rprd(:,:)     ! rain production rate (ncol,pver)
         301             :    real(kind_phys), intent(out) :: dif(:,:)        ! detrained convective cloud ice mixing ratio.         (ncol,pver)
         302             :    real(kind_phys), intent(out) :: dnlf(:,:)       ! detrained convective cloud water num concen.         (ncol,pver)
         303             :    real(kind_phys), intent(out) :: dnif(:,:)       ! detrained convective cloud ice num concen.           (ncol,pver)
         304             :
         305             : ! move these vars from local storage to output so that convective
         306             : ! transports can be done in outside of conv_cam.
         307             :    real(kind_phys), intent(out) :: mu(:,:)  !                                                                 (ncol,pver)
         308             :    real(kind_phys), intent(out) :: eu(:,:)  !                                                                 (ncol,pver)
         309             :    real(kind_phys), intent(out) :: du(:,:)  !                                                                 (ncol,pver)
         310             :    real(kind_phys), intent(out) :: md(:,:)  !                                                                 (ncol,pver)
         311             :    real(kind_phys), intent(out) :: ed(:,:)  !                                                                 (ncol,pver)
         312             :    real(kind_phys), intent(out) :: dp(:,:)       ! wg layer thickness in mbs (between upper/lower interface). (ncol,pver)
         313             :    real(kind_phys), intent(out) :: dsubcld(:)       ! wg layer thickness in mbs between lcl and maxi.         (ncol)
         314             :    real(kind_phys), intent(out) :: prec(:)  !                                                                 (ncol)
         315             :    real(kind_phys), intent(out) :: rliq(:) ! reserved liquid (not yet in cldliq) for energy integrals         (ncol)
         316             :    real(kind_phys), intent(out) :: rice(:) ! reserved ice (not yet in cldce) for energy integrals             (ncol)
         317             :
         318             :    integer,  intent(out) :: ideep(:)  ! column indices of gathered points                              (ncol)
         319             :
         320             :    character(len=512), intent(out)      :: errmsg
         321             :    integer, intent(out)                 :: errflg
         322             :
         323             :
         324             :    real(kind_phys), intent(in) :: org(:,:)     ! Only used if zm_org is true  ! in
         325             :    real(kind_phys), intent(out) :: orgt(:,:)   ! Only used if zm_org is true   ! out
         326             :    real(kind_phys), intent(out) :: org2d(:,:)  ! Only used if zm_org is true   ! out
         327             :
         328             :    ! Local variables
         329             :
         330     2990736 :    real(kind_phys) zs(ncol)
         331     2990736 :    real(kind_phys) dlg(ncol,pver)    ! gathrd version of the detraining cld h2o tend
         332     2990736 :    real(kind_phys) cug(ncol,pver)    ! gathered condensation rate
         333             :
         334     2990736 :    real(kind_phys) evpg(ncol,pver)   ! gathered evap rate of rain in downdraft
         335     2990736 :    real(kind_phys) orgavg(ncol)
         336     2990736 :    real(kind_phys) dptot(ncol)
         337             :
         338     2990736 :    real(kind_phys) mumax(ncol)
         339             :
         340             : !CACNOTE - Figure out real intent for jt and maxg
         341             :    integer, intent(inout) :: jt(ncol)                          ! wg top  level index of deep cumulus convection.
         342             :    integer, intent(inout) :: maxg(ncol)                        ! wg gathered values of maxi.
         343             :
         344             :    integer lengath
         345             : !     diagnostic field used by chem/wetdep codes
         346             :
         347             : !CACNOTE - Figure out real intent for ql
         348             :    real(kind_phys),intent(inout):: ql(ncol,pver)                    ! wg grid slice of cloud liquid water.
         349             : !
         350     2990736 :    real(kind_phys) pblt(ncol)           ! i row of pbl top indices.
         351             :
         352             :
         353             :
         354             :
         355             : !
         356             : !-----------------------------------------------------------------------
         357             : !
         358             : ! general work fields (local variables):
         359             : !
         360     2990736 :    real(kind_phys) q(ncol,pver)              ! w  grid slice of mixing ratio.
         361     2990736 :    real(kind_phys) p(ncol,pver)              ! w  grid slice of ambient mid-layer pressure in mbs.
         362     2990736 :    real(kind_phys) z(ncol,pver)              ! w  grid slice of ambient mid-layer height in metres.
         363     2990736 :    real(kind_phys) s(ncol,pver)              ! w  grid slice of scaled dry static energy (t+gz/cp).
         364     2990736 :    real(kind_phys) tp(ncol,pver)             ! w  grid slice of parcel temperatures.
         365     2990736 :    real(kind_phys) zf(ncol,pver+1)           ! w  grid slice of ambient interface height in metres.
         366     2990736 :    real(kind_phys) pf(ncol,pver+1)           ! w  grid slice of ambient interface pressure in mbs.
         367     2990736 :    real(kind_phys) qstp(ncol,pver)           ! w  grid slice of parcel temp. saturation mixing ratio.
         368             :
         369     2990736 :    real(kind_phys) tl(ncol)                  ! w  row of parcel temperature at lcl.
         370             :
         371     2990736 :    integer lcl(ncol)                  ! w  base level index of deep cumulus convection.
         372     2990736 :    integer lel(ncol)                  ! w  index of highest theoretical convective plume.
         373     2990736 :    integer lon(ncol)                  ! w  index of onset level for deep convection.
         374     2990736 :    integer maxi(ncol)                 ! w  index of level with largest moist static energy.
         375             :
         376             :    real(kind_phys) precip
         377             : !
         378             : ! gathered work fields:
         379             : !
         380     2990736 :    real(kind_phys) qg(ncol,pver)             ! wg grid slice of gathered values of q.
         381     2990736 :    real(kind_phys) tg(ncol,pver)             ! w  grid slice of temperature at interface.
         382     2990736 :    real(kind_phys) pg(ncol,pver)             ! wg grid slice of gathered values of p.
         383     2990736 :    real(kind_phys) zg(ncol,pver)             ! wg grid slice of gathered values of z.
         384     2990736 :    real(kind_phys) sg(ncol,pver)             ! wg grid slice of gathered values of s.
         385     2990736 :    real(kind_phys) tpg(ncol,pver)            ! wg grid slice of gathered values of tp.
         386     2990736 :    real(kind_phys) zfg(ncol,pver+1)          ! wg grid slice of gathered values of zf.
         387     2990736 :    real(kind_phys) qstpg(ncol,pver)          ! wg grid slice of gathered values of qstp.
         388     2990736 :    real(kind_phys) ug(ncol,pver)             ! wg grid slice of gathered values of u.
         389     2990736 :    real(kind_phys) vg(ncol,pver)             ! wg grid slice of gathered values of v.
         390     2990736 :    real(kind_phys) cmeg(ncol,pver)
         391             :
         392     2990736 :    real(kind_phys) rprdg(ncol,pver)           ! wg gathered rain production rate
         393     2990736 :    real(kind_phys) capeg(ncol)               ! wg gathered convective available potential energy.
         394     2990736 :    real(kind_phys) tlg(ncol)                 ! wg grid slice of gathered values of tl.
         395     2990736 :    real(kind_phys) landfracg(ncol)            ! wg grid slice of landfrac
         396             :
         397     2990736 :    integer lclg(ncol)       ! wg gathered values of lcl.
         398     2990736 :    integer lelg(ncol)
         399             : !
         400             : ! work fields arising from gathered calculations.
         401             : !
         402     2990736 :    real(kind_phys) dqdt(ncol,pver)           ! wg mixing ratio tendency at gathered points.
         403     2990736 :    real(kind_phys) dsdt(ncol,pver)           ! wg dry static energy ("temp") tendency at gathered points.
         404     2990736 :    real(kind_phys) sd(ncol,pver)             ! wg grid slice of dry static energy in downdraft.
         405     2990736 :    real(kind_phys) qd(ncol,pver)             ! wg grid slice of mixing ratio in downdraft.
         406     2990736 :    real(kind_phys) mc(ncol,pver)             ! wg net upward (scaled by mb) cloud mass flux.
         407     2990736 :    real(kind_phys) qhat(ncol,pver)           ! wg grid slice of upper interface mixing ratio.
         408     2990736 :    real(kind_phys) qu(ncol,pver)             ! wg grid slice of mixing ratio in updraft.
         409     2990736 :    real(kind_phys) su(ncol,pver)             ! wg grid slice of dry static energy in updraft.
         410     2990736 :    real(kind_phys) qs(ncol,pver)             ! wg grid slice of saturation mixing ratio.
         411     2990736 :    real(kind_phys) shat(ncol,pver)           ! wg grid slice of upper interface dry static energy.
         412     2990736 :    real(kind_phys) hmn(ncol,pver)            ! wg moist static energy.
         413     2990736 :    real(kind_phys) hsat(ncol,pver)           ! wg saturated moist static energy.
         414     2990736 :    real(kind_phys) qlg(ncol,pver)
         415     2990736 :    real(kind_phys) dudt(ncol,pver)           ! wg u-wind tendency at gathered points.
         416     2990736 :    real(kind_phys) dvdt(ncol,pver)           ! wg v-wind tendency at gathered points.
         417             :
         418     2990736 :    real(kind_phys) qldeg(ncol,pver)        ! cloud liquid water mixing ratio for detrainment (kg/kg)
         419     2990736 :    real(kind_phys) mb(ncol)                ! wg cloud base mass flux.
         420             :
         421     2990736 :    integer jlcl(ncol)
         422     2990736 :    integer j0(ncol)                 ! wg detrainment initiation level index.
         423     1495368 :    integer jd(ncol)                 ! wg downdraft initiation level index.
         424             :
         425             :    real(kind_phys),intent(in):: delt                     ! length of model time-step in seconds.
         426             :
         427             :    integer i
         428             :    integer ii
         429             :    integer k, kk, l, m
         430             :
         431             :    integer msg                      !  ic number of missing moisture levels at the top of model.
         432             :    real(kind_phys) qdifr
         433             :    real(kind_phys) sdifr
         434             :
         435             :    real(kind_phys), parameter :: dcon  = 25.e-6_kind_phys
         436             :    real(kind_phys), parameter :: mucon = 5.3_kind_phys
         437             :    real(kind_phys) negadq
         438             :    logical doliq
         439             :
         440             :
         441             : !
         442             : !--------------------------Data statements------------------------------
         443             :
         444     1495368 :    errmsg = ''
         445     1495368 :    errflg = 0
         446             : !
         447             : ! Set internal variable "msg" (convection limit) to "limcnv-1"
         448             : !
         449     1495368 :    msg = limcnv - 1
         450             : !
         451             : ! initialize necessary arrays.
         452             : ! zero out variables not used in cam
         453             : !
         454             :
         455     1495368 :    if (zm_org) then
         456           0 :       orgt(:,:) = 0._kind_phys
         457             :    end if
         458             :
         459  2323627992 :    qtnd(:,:) = 0._kind_phys
         460  2323627992 :    heat(:,:) = 0._kind_phys
         461  2348597160 :    mcon(:,:) = 0._kind_phys
         462    24969168 :    rliq(:ncol)   = 0._kind_phys
         463    24969168 :    rice(:ncol)   = 0._kind_phys
         464             :
         465             : !
         466             : ! initialize convective tendencies
         467             : !
         468    24969168 :    prec(:ncol) = 0._kind_phys
         469   140564592 :    do k = 1,pver
         470  2323627992 :       do i = 1,ncol
         471  2183063400 :          dqdt(i,k)  = 0._kind_phys
         472  2183063400 :          dsdt(i,k)  = 0._kind_phys
         473  2183063400 :          dudt(i,k)  = 0._kind_phys
         474  2183063400 :          dvdt(i,k)  = 0._kind_phys
         475  2183063400 :          cme(i,k)   = 0._kind_phys
         476  2183063400 :          rprd(i,k)  = 0._kind_phys
         477  2183063400 :          zdu(i,k)   = 0._kind_phys
         478  2183063400 :          ql(i,k)    = 0._kind_phys
         479  2183063400 :          qlg(i,k)   = 0._kind_phys
         480  2183063400 :          dlf(i,k)   = 0._kind_phys
         481  2183063400 :          dlg(i,k)   = 0._kind_phys
         482  2183063400 :          qldeg(i,k) = 0._kind_phys
         483             :
         484  2183063400 :          dif(i,k)   = 0._kind_phys
         485  2183063400 :          dnlf(i,k)  = 0._kind_phys
         486  2322132624 :          dnif(i,k)  = 0._kind_phys
         487             :
         488             :       end do
         489             :    end do
         490             :
         491    24969168 :    do i = 1,ncol
         492    23473800 :       pblt(i) = pver
         493    24969168 :       dsubcld(i) = 0._kind_phys
         494             :
         495             :    end do
         496             :
         497     1495368 :   if (zm_org) then
         498             : ! compute vertical average here
         499           0 :       orgavg(:) = 0._kind_phys
         500           0 :       dptot(:) = 0._kind_phys
         501             :
         502           0 :       do k = 1, pver
         503           0 :         do i = 1,ncol
         504           0 :           if (org(i,k) .gt. 0) then
         505           0 :             orgavg(i) = orgavg(i)+dpp(i,k)*org(i,k)
         506           0 :             dptot(i) = dptot(i)+dpp(i,k)
         507             :           endif
         508             :         enddo
         509             :       enddo
         510             :
         511           0 :       do i = 1,ncol
         512           0 :         if (dptot(i) .gt. 0) then
         513           0 :           orgavg(i) = orgavg(i)/dptot(i)
         514             :         endif
         515             :       enddo
         516             :
         517           0 :       do k = 1, pver
         518           0 :         do i = 1, ncol
         519           0 :            org2d(i,k) = orgavg(i)
         520             :         enddo
         521             :       enddo
         522             :
         523             :    endif
         524             :
         525             : !
         526             : ! calculate local pressure (mbs) and height (m) for both interface
         527             : ! and mid-layer locations.
         528             : !
         529    24969168 :    do i = 1,ncol
         530    23473800 :       zs(i) = geos(i)*rgrav
         531    23473800 :       pf(i,pver+1) = paph(i,pver+1)*0.01_kind_phys
         532    24969168 :       zf(i,pver+1) = zi(i,pver+1) + zs(i)
         533             :    end do
         534   140564592 :    do k = 1,pver
         535  2323627992 :       do i = 1,ncol
         536  2183063400 :          p(i,k) = pap(i,k)*0.01_kind_phys
         537  2183063400 :          pf(i,k) = paph(i,k)*0.01_kind_phys
         538  2183063400 :          z(i,k) = zm(i,k) + zs(i)
         539  2322132624 :          zf(i,k) = zi(i,k) + zs(i)
         540             :       end do
         541             :    end do
         542             : !
         543    89722080 :    do k = pver - 1,msg + 1,-1
         544  1474676280 :       do i = 1,ncol
         545  1473180912 :          if (abs(z(i,k)-zs(i)-pblh(i)) < (zf(i,k)-zf(i,k+1))*0.5_kind_phys) pblt(i) = k
         546             :       end do
         547             :    end do
         548             : !
         549             : ! store incoming specific humidity field for subsequent calculation
         550             : ! of precipitation (through change in storage).
         551             : ! define dry static energy (normalized by cp).
         552             : !
         553   140564592 :    do k = 1,pver
         554  2323627992 :       do i = 1,ncol
         555  2183063400 :          q(i,k) = qh(i,k)
         556  2183063400 :          s(i,k) = t(i,k) + (grav/cpres)*z(i,k)
         557  2183063400 :          tp(i,k)=0.0_kind_phys
         558  2183063400 :          shat(i,k) = s(i,k)
         559  2322132624 :          qhat(i,k) = q(i,k)
         560             :       end do
         561             :    end do
         562             :
         563    24969168 :    do i = 1,ncol
         564    23473800 :       capeg(i) = 0._kind_phys
         565    23473800 :       lclg(i) = 1
         566    23473800 :       lelg(i) = pver
         567    23473800 :       maxg(i) = 1
         568    23473800 :       tlg(i) = 400._kind_phys
         569    24969168 :       dsubcld(i) = 0._kind_phys
         570             :    end do
         571             :
         572             :
         573             :    !  Evaluate Tparcel, qs(Tparcel), buoyancy and CAPE,
         574             :    !     lcl, lel, parcel launch level at index maxi()=hmax
         575             :
         576             :    call buoyan_dilute(  ncol   ,pver     , &
         577             :                cpliq   ,latice  ,cpwv    ,rh2o    ,&
         578             :                q       ,t       ,p       ,z       ,pf       , &
         579             :                tp      ,qstp    ,tl      ,rl      ,cape     , &
         580             :                pblt    ,lcl     ,lel     ,lon     ,maxi     , &
         581             :                rgas    ,grav    ,cpres   ,msg     , &
         582             :                zi      ,zs      ,tpert   , org2d  , landfrac,&
         583  1096185888 :                errmsg  ,errflg)
         584             :
         585             : !
         586             : ! determine whether grid points will undergo some deep convection
         587             : ! (ideep=1) or not (ideep=0), based on values of cape,lcl,lel
         588             : ! (require cape.gt. 0 and lel<lcl as minimum conditions).
         589             : !
         590     1495368 :    lengath = 0
         591    24969168 :    ideep   = 0
         592    24969168 :    do i=1,ncol
         593    24969168 :       if (cape(i) > capelmt) then
         594     4224820 :          lengath = lengath + 1
         595     4224820 :          ideep(lengath) = i
         596             :       end if
         597             :    end do
         598             :
         599     1495368 :    if (lengath.eq.0) return
         600             : !
         601             : ! obtain gathered arrays necessary for ensuing calculations.
         602             : !
         603    54383700 :    do k = 1,pver
         604   447291960 :       do i = 1,lengath
         605   392908260 :          dp(i,k) = 0.01_kind_phys*dpp(ideep(i),k)
         606   392908260 :          qg(i,k) = q(ideep(i),k)
         607   392908260 :          tg(i,k) = t(ideep(i),k)
         608   392908260 :          pg(i,k) = p(ideep(i),k)
         609   392908260 :          zg(i,k) = z(ideep(i),k)
         610   392908260 :          sg(i,k) = s(ideep(i),k)
         611   392908260 :          tpg(i,k) = tp(ideep(i),k)
         612   392908260 :          zfg(i,k) = zf(ideep(i),k)
         613   392908260 :          qstpg(i,k) = qstp(ideep(i),k)
         614   392908260 :          ug(i,k) = 0._kind_phys
         615   446713410 :          vg(i,k) = 0._kind_phys
         616             :       end do
         617             :    end do
         618             :
         619             : !
         620     4803370 :    do i = 1,lengath
         621     4803370 :       zfg(i,pver+1) = zf(ideep(i),pver+1)
         622             :    end do
         623     4803370 :    do i = 1,lengath
         624     4224820 :       capeg(i) = cape(ideep(i))
         625     4224820 :       lclg(i) = lcl(ideep(i))
         626     4224820 :       lelg(i) = lel(ideep(i))
         627     4224820 :       maxg(i) = maxi(ideep(i))
         628     4224820 :       tlg(i) = tl(ideep(i))
         629     4803370 :       landfracg(i) = landfrac(ideep(i))
         630             :    end do
         631             : !
         632             : ! calculate sub-cloud layer pressure "thickness" for use in
         633             : ! closure and tendency routines.
         634             : !
         635    35291550 :    do k = msg + 1,pver
         636   288780750 :       do i = 1,lengath
         637   288202200 :          if (k >= maxg(i)) then
         638    23218179 :             dsubcld(i) = dsubcld(i) + dp(i,k)
         639             :          end if
         640             :       end do
         641             :    end do
         642             : !
         643             : ! define array of factors (alpha) which defines interfacial
         644             : ! values, as well as interfacial values for (q,s) used in
         645             : ! subsequent routines.
         646             : !
         647    34713000 :    do k = msg + 2,pver
         648   283977380 :       do i = 1,lengath
         649   249264380 :          sdifr = 0._kind_phys
         650   249264380 :          qdifr = 0._kind_phys
         651   249264380 :          if (sg(i,k) > 0._kind_phys .or. sg(i,k-1) > 0._kind_phys) &
         652   249264380 :             sdifr = abs((sg(i,k)-sg(i,k-1))/max(sg(i,k-1),sg(i,k)))
         653   249264380 :          if (qg(i,k) > 0._kind_phys .or. qg(i,k-1) > 0._kind_phys) &
         654   249264380 :             qdifr = abs((qg(i,k)-qg(i,k-1))/max(qg(i,k-1),qg(i,k)))
         655   249264380 :          if (sdifr > 1.E-6_kind_phys) then
         656   249166037 :             shat(i,k) = log(sg(i,k-1)/sg(i,k))*sg(i,k-1)*sg(i,k)/(sg(i,k-1)-sg(i,k))
         657             :          else
         658       98343 :             shat(i,k) = 0.5_kind_phys* (sg(i,k)+sg(i,k-1))
         659             :          end if
         660   283398830 :          if (qdifr > 1.E-6_kind_phys) then
         661   249261709 :             qhat(i,k) = log(qg(i,k-1)/qg(i,k))*qg(i,k-1)*qg(i,k)/(qg(i,k-1)-qg(i,k))
         662             :          else
         663        2671 :             qhat(i,k) = 0.5_kind_phys* (qg(i,k)+qg(i,k-1))
         664             :          end if
         665             :       end do
         666             :    end do
         667             : !
         668             : ! obtain cloud properties.
         669             : !
         670             :
         671             :    call cldprp(ncol   ,pver    ,pverp   ,cpliq  , &
         672             :                latice  ,cpwv    ,rh2o    ,&
         673             :                qg      ,tg      ,ug      ,vg      ,pg      , &
         674             :                zg      ,sg      ,mu      ,eu      ,du      , &
         675             :                md      ,ed      ,sd      ,qd      ,mc      , &
         676             :                qu      ,su      ,zfg     ,qs      ,hmn     , &
         677             :                hsat    ,shat    ,qlg     , &
         678             :                cmeg    ,maxg    ,lelg    ,jt      ,jlcl    , &
         679             :                maxg    ,j0      ,jd      ,rl      ,lengath , &
         680             :                rgas    ,grav    ,cpres   ,msg     , &
         681             :                evpg    ,cug     ,rprdg   ,limcnv  ,landfracg , &
         682  1045949235 :                qldeg    ,qhat    )
         683             :
         684             :
         685             : !
         686             : ! convert detrainment from units of "1/m" to "1/mb".
         687             : !
         688             :
         689    35291550 :    do k = msg + 1,pver
         690   288780750 :       do i = 1,lengath
         691   253489200 :          du   (i,k) = du   (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         692   253489200 :          eu   (i,k) = eu   (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         693   253489200 :          ed   (i,k) = ed   (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         694   253489200 :          cug  (i,k) = cug  (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         695   253489200 :          cmeg (i,k) = cmeg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         696   253489200 :          rprdg(i,k) = rprdg(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         697   288202200 :          evpg (i,k) = evpg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         698             :       end do
         699             :    end do
         700             :
         701             :    call closure(ncol   ,pver    , &
         702             :                 qg      ,tg      ,pg      ,zg      ,sg      , &
         703             :                 tpg     ,qs      ,qu      ,su      ,mc      , &
         704             :                 du      ,mu      ,md      ,qd      ,sd      , &
         705             :                 qhat    ,shat    ,dp      ,qstpg   ,zfg     , &
         706             :                 qlg     ,dsubcld ,mb      ,capeg   ,tlg     , &
         707             :                 lclg    ,lelg    ,jt      ,maxg    ,1       , &
         708             :                 lengath ,rgas    ,grav    ,cpres   ,rl      , &
         709   836875098 :                 msg     ,capelmt    )
         710             : !
         711             : ! limit cloud base mass flux to theoretical upper bound.
         712             : !
         713     4803370 :    do i=1,lengath
         714     4803370 :       mumax(i) = 0
         715             :    end do
         716    34713000 :    do k=msg + 2,pver
         717   283977380 :       do i=1,lengath
         718   283398830 :         mumax(i) = max(mumax(i), mu(i,k)/dp(i,k))
         719             :       end do
         720             :    end do
         721             :
         722     4803370 :    do i=1,lengath
         723     4803370 :       if (mumax(i) > 0._kind_phys) then
         724     4208860 :          mb(i) = min(mb(i),1._kind_phys/(delt*mumax(i)))
         725             :       else
         726       15960 :          mb(i) = 0._kind_phys
         727             :       endif
         728             :    end do
         729             :    ! If no_deep_pbl = .true., don't allow convection entirely
         730             :    ! within PBL (suggestion of Bjorn Stevens, 8-2000)
         731             :
         732      578550 :    if (no_deep_pbl) then
         733           0 :       do i=1,lengath
         734           0 :          if (zm(ideep(i),jt(i)) < pblh(ideep(i))) mb(i) = 0
         735             :       end do
         736             :    end if
         737             :
         738    35291550 :    do k=msg+1,pver
         739   288780750 :       do i=1,lengath
         740   253489200 :          mu   (i,k)  = mu   (i,k)*mb(i)
         741   253489200 :          md   (i,k)  = md   (i,k)*mb(i)
         742   253489200 :          mc   (i,k)  = mc   (i,k)*mb(i)
         743   253489200 :          du   (i,k)  = du   (i,k)*mb(i)
         744   253489200 :          eu   (i,k)  = eu   (i,k)*mb(i)
         745   253489200 :          ed   (i,k)  = ed   (i,k)*mb(i)
         746   253489200 :          cmeg (i,k)  = cmeg (i,k)*mb(i)
         747   253489200 :          rprdg(i,k)  = rprdg(i,k)*mb(i)
         748   253489200 :          cug  (i,k)  = cug  (i,k)*mb(i)
         749   288202200 :          evpg (i,k)  = evpg (i,k)*mb(i)
         750             :
         751             :       end do
         752             :    end do
         753             : !
         754             : ! compute temperature and moisture changes due to convection.
         755             : !
         756             :    call q1q2_pjr(ncol   ,pver    ,latice  , &
         757             :                  dqdt    ,dsdt    ,qg      ,qs      ,qu      , &
         758             :                  su      ,du      ,qhat    ,shat    ,dp      , &
         759             :                  mu      ,md      ,sd      ,qd      ,qldeg   , &
         760             :                  dsubcld ,jt      ,maxg    ,1       ,lengath , &
         761             :                  cpres   ,rl      ,msg     ,          &
         762   836875098 :                  dlg     ,evpg    ,cug)
         763             :
         764             : !
         765             : ! gather back temperature and mixing ratio.
         766             : !
         767             :
         768    35291550 :    do k = msg + 1,pver
         769   288780750 :       do i = 1,lengath
         770             : !
         771             : ! q is updated to compute net precip.
         772             : !
         773   253489200 :          q(ideep(i),k) = qh(ideep(i),k) + delt*dqdt(i,k)
         774   253489200 :          qtnd(ideep(i),k) = dqdt (i,k)
         775   253489200 :          cme (ideep(i),k) = cmeg (i,k)
         776   253489200 :          rprd(ideep(i),k) = rprdg(i,k)
         777   253489200 :          zdu (ideep(i),k) = du   (i,k)
         778   253489200 :          mcon(ideep(i),k) = mc   (i,k)
         779   253489200 :          heat(ideep(i),k) = dsdt (i,k)*cpres
         780   253489200 :          dlf (ideep(i),k) = dlg  (i,k)
         781   288202200 :          ql  (ideep(i),k) = qlg  (i,k)
         782             :       end do
         783             :    end do
         784             :
         785             : ! Compute precip by integrating change in water vapor minus detrained cloud water
         786    35291550 :    do k = pver,msg + 1,-1
         787   581488650 :       do i = 1,ncol
         788   580910100 :           prec(i) = prec(i) - dpp(i,k)* (q(i,k)-qh(i,k)) - dpp(i,k)*(dlf(i,k)+dif(i,k))*delt
         789             :       end do
         790             :    end do
         791             :
         792             : ! obtain final precipitation rate in m/s.
         793     9681835 :    do i = 1,ncol
         794     9681835 :       prec(i) = rgrav*max(prec(i),0._kind_phys)/ delt/1000._kind_phys
         795             :    end do
         796             :
         797             : ! Compute reserved liquid (not yet in cldliq) for energy integrals.
         798             : ! Treat rliq as flux out bottom, to be added back later.
         799    54383700 :    do k = 1, pver
         800   900989205 :       do i = 1, ncol
         801   846605505 :           rliq(i) = rliq(i) + (dlf(i,k)+dif(i,k))*dpp(i,k)/gravit
         802   900410655 :           rice(i) = rice(i) + dif(i,k)*dpp(i,k)/gravit
         803             :       end do
         804             :    end do
         805     9681835 :    rliq(:ncol) = rliq(:ncol) /1000._kind_phys
         806     9681835 :    rice(:ncol) = rice(:ncol) /1000._kind_phys
         807             :
         808             :    return
         809             : end subroutine zm_convr_run
         810             :
         811             : subroutine zm_convr_finalize
         812             : end subroutine zm_convr_finalize
         813             :
         814             : !=========================================================================================
         815             :
         816     1495368 : subroutine buoyan_dilute(  ncol   ,pver    , &
         817             :                   cpliq   ,latice  ,cpwv    ,rh2o    ,&
         818     1495368 :                   q       ,t       ,p       ,z       ,pf      , &
         819     1495368 :                   tp      ,qstp    ,tl      ,rl      ,cape    , &
         820     1495368 :                   pblt    ,lcl     ,lel     ,lon     ,mx      , &
         821             :                   rd      ,grav    ,cp      ,msg     , &
         822     1495368 :                   zi      ,zs      ,tpert    ,org    , landfrac,&
         823           0 :                   errmsg  ,errflg)
         824             : !-----------------------------------------------------------------------
         825             : !
         826             : ! Purpose:
         827             : ! Calculates CAPE the lifting condensation level and the convective top
         828             : ! where buoyancy is first -ve.
         829             : !
         830             : ! Method: Calculates the parcel temperature based on a simple constant
         831             : ! entraining plume model. CAPE is integrated from buoyancy.
         832             : ! 09/09/04 - Simplest approach using an assumed entrainment rate for
         833             : !            testing (dmpdp).
         834             : ! 08/04/05 - Swap to convert dmpdz to dmpdp
         835             : !
         836             : ! SCAM Logical Switches - DILUTE:RBN - Now Disabled
         837             : ! ---------------------
         838             : ! switch(1) = .T. - Uses the dilute parcel calculation to obtain tendencies.
         839             : ! switch(2) = .T. - Includes entropy/q changes due to condensate loss and freezing.
         840             : ! switch(3) = .T. - Adds the PBL Tpert for the parcel temperature at all levels.
         841             : !
         842             : ! References:
         843             : ! Raymond and Blythe (1992) JAS
         844             : !
         845             : ! Author:
         846             : ! Richard Neale - September 2004
         847             : !
         848             : !-----------------------------------------------------------------------
         849             :    implicit none
         850             : !-----------------------------------------------------------------------
         851             : !
         852             : ! input arguments
         853             : !
         854             :    integer, intent(in) :: ncol                  ! number of atmospheric columns
         855             :    integer, intent(in) :: pver
         856             :    real(kind_phys), intent(in) :: cpliq
         857             :    real(kind_phys), intent(in) :: latice
         858             :    real(kind_phys), intent(in) :: cpwv
         859             :    real(kind_phys), intent(in) :: rh2o
         860             :
         861             :    real(kind_phys), intent(in) :: q(ncol,pver)        ! spec. humidity
         862             :    real(kind_phys), intent(in) :: t(ncol,pver)        ! temperature
         863             :    real(kind_phys), intent(in) :: p(ncol,pver)        ! pressure
         864             :    real(kind_phys), intent(in) :: z(ncol,pver)        ! height
         865             :    real(kind_phys), intent(in) :: pf(ncol,pver+1)     ! pressure at interfaces
         866             :    real(kind_phys), intent(in) :: pblt(ncol)          ! index of pbl depth
         867             :    real(kind_phys), intent(in) :: tpert(ncol)         ! perturbation temperature by pbl processes
         868             :
         869             : ! Use z interface/surface relative values for PBL parcel calculations.
         870             :    real(kind_phys), intent(in) :: zi(ncol,pver+1)
         871             :    real(kind_phys), intent(in) :: zs(ncol)
         872             :
         873             : !
         874             : ! output arguments
         875             : !
         876             :
         877             :    real(kind_phys), intent(out) :: tp(ncol,pver)       ! parcel temperature
         878             :    real(kind_phys), intent(out) :: qstp(ncol,pver)     ! saturation mixing ratio of parcel (only above lcl, just q below).
         879             :    real(kind_phys), intent(out) :: tl(ncol)            ! parcel temperature at lcl
         880             :    real(kind_phys), intent(out) :: cape(ncol)          ! convective aval. pot. energy.
         881             :    integer lcl(ncol)        !
         882             :    integer lel(ncol)        !
         883             :    integer lon(ncol)        ! level of onset of deep convection
         884             :    integer mx(ncol)         ! level of max moist static energy
         885             :
         886             :    real(kind_phys)  :: org(:,:)      ! organization parameter
         887             :    real(kind_phys), intent(in) :: landfrac(ncol)
         888             :    character(len=512), intent(out)      :: errmsg
         889             :    integer, intent(out)                 :: errflg
         890             :
         891             : !
         892             : !--------------------------Local Variables------------------------------
         893             : !
         894     2990736 :    real(kind_phys) capeten(ncol,5)     ! provisional value of cape
         895     2990736 :    real(kind_phys) tv(ncol,pver)       !
         896     2990736 :    real(kind_phys) tpv(ncol,pver)      !
         897     2990736 :    real(kind_phys) buoy(ncol,pver)
         898             :
         899             :    real(kind_phys) a1(ncol)
         900             :    real(kind_phys) a2(ncol)
         901             :    real(kind_phys) estp(ncol)
         902     2990736 :    real(kind_phys) pl(ncol)
         903             :    real(kind_phys) plexp(ncol)
         904     2990736 :    real(kind_phys) hmax(ncol)
         905     2990736 :    real(kind_phys) hmn(ncol)
         906             :    real(kind_phys) y(ncol)
         907             :
         908     2990736 :    logical plge600(ncol)
         909     2990736 :    integer knt(ncol)
         910     2990736 :    integer lelten(ncol,5)
         911             :
         912             :
         913             :
         914             :
         915             : ! Parcel property variables
         916             :
         917     2990736 :   real(kind_phys)           :: hmn_lev(ncol,pver)  ! Vertical profile of moist static energy for each column
         918     2990736 :   real(kind_phys)           :: dp_lev(ncol,pver)   ! Level dpressure between interfaces
         919     2990736 :   real(kind_phys)           :: hmn_zdp(ncol,pver)  ! Integrals of hmn_lev*dp_lev at each level
         920     2990736 :   real(kind_phys)           :: q_zdp(ncol,pver)    ! Integrals of q*dp_lev at each level
         921             :   real(kind_phys)           :: dp_zfrac             ! Fraction of vertical grid box below mixing top (usually pblt)
         922     2990736 :   real(kind_phys)           :: parcel_dz(ncol)     ! Depth of parcel mixing (usually parcel_hscale*parcel_dz)
         923     2990736 :   real(kind_phys)           :: parcel_ztop(ncol)   ! Height of parcel mixing (usually parcel_ztop+zm(nlev))
         924     2990736 :   real(kind_phys)           :: parcel_dp(ncol)     ! Pressure integral over parcel mixing depth (usually pblt)
         925     2990736 :   real(kind_phys)           :: parcel_hdp(ncol)    ! Pressure*MSE integral over parcel mixing depth (usually pblt)
         926     2990736 :   real(kind_phys)           :: parcel_qdp(ncol)    ! Pressure*q integral over parcel mixing depth (usually pblt)
         927     2990736 :   real(kind_phys)           :: pbl_dz(ncol)        ! Previously diagnosed PBL height
         928     2990736 :   real(kind_phys)           :: hpar(ncol)          ! Initial MSE of the parcel
         929     2990736 :   real(kind_phys)           :: qpar(ncol)          ! Initial humidity of the parcel
         930     1495368 :   real(kind_phys)           :: ql(ncol)          ! Initial parcel humidity (for ientropy routine)
         931             :   integer            :: ipar ! Index for top of parcel mixing/launch level.
         932             :
         933             :
         934             :
         935             :
         936             :    real(kind_phys) cp
         937             :    real(kind_phys) e
         938             :    real(kind_phys) grav
         939             :
         940             :    integer i
         941             :    integer k
         942             :    integer msg
         943             :    integer n
         944             :
         945             :    real(kind_phys) rd
         946             :    real(kind_phys) rl
         947             :
         948             :
         949             : ! Scaling of PBL height to give parcel mixing length for lparcel_pbl=True
         950             :
         951             :    real(kind_phys), parameter :: parcel_hscale  = 0.5_kind_phys
         952             :
         953             :
         954             : !
         955             : !-----------------------------------------------------------------------
         956             : !
         957     8972208 :    do n = 1,5
         958   126341208 :       do i = 1,ncol
         959   117369000 :          lelten(i,n) = pver
         960   124845840 :          capeten(i,n) = 0._kind_phys
         961             :       end do
         962             :    end do
         963             : !
         964    24969168 :    do i = 1,ncol
         965    23473800 :       lon(i) = pver
         966    23473800 :       knt(i) = 0
         967    23473800 :       lel(i) = pver
         968    23473800 :       mx(i) = lon(i)
         969    23473800 :       cape(i) = 0._kind_phys
         970    23473800 :       hmax(i) = 0._kind_phys
         971    23473800 :       pbl_dz(i) = z(i,nint(pblt(i)))-zs(i) ! mid-point z (zm) reference to PBL depth
         972    23473800 :       parcel_dz(i) = max(zi(i,pver),parcel_hscale*pbl_dz(i)) ! PBL mixing depth [parcel_hscale*Boundary, but no thinner than zi(i,pver)]
         973    23473800 :       parcel_ztop(i) = parcel_dz(i)+zs(i) ! PBL mixing height ztop this is wrt zs=0
         974    23473800 :       parcel_hdp(i) = 0._kind_phys
         975    23473800 :       parcel_dp(i) = 0._kind_phys
         976    23473800 :       parcel_qdp(i) = 0._kind_phys
         977    23473800 :       hpar(i) = 0._kind_phys
         978    24969168 :       qpar(i) = 0._kind_phys
         979             :    end do
         980             :
         981  2323627992 :    tp(:ncol,:) = t(:ncol,:)
         982  2323627992 :    qstp(:ncol,:) = q(:ncol,:)
         983  2323627992 :    hmn_lev(:ncol,:) = 0._kind_phys
         984             :
         985             :
         986             :
         987             : !!! Initialize tv and buoy for output.
         988             : !!! tv=tv : tpv=tpv : qstp=q : buoy=0.
         989  2323627992 :    tv(:ncol,:) = t(:ncol,:) *(1._kind_phys+1.608_kind_phys*q(:ncol,:))/ (1._kind_phys+q(:ncol,:))
         990  2323627992 :    tpv(:ncol,:) = tv(:ncol,:)
         991  2323627992 :    buoy(:ncol,:) = 0._kind_phys
         992             :
         993             :
         994             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         995             : ! Mix the parcel over a certain dp or dz and take the launch level as the top level
         996             : ! of this mixing region and the parcel properties as this mixed value
         997             : ! Should be well mixed by other processes in the very near PBL.
         998             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         999             :
        1000             :
        1001             :
        1002     1495368 : if (lparcel_pbl) then
        1003             :
        1004             : ! Vertical profile of MSE and pressure weighted of the same.
        1005  2323627992 :    hmn_lev(:ncol,1:pver) = cp*t(:ncol,1:pver) + grav*z(:ncol,1:pver) + rl*q(:ncol,1:pver)
        1006  2323627992 :    dp_lev(:ncol,1:pver) = pf(:ncol,2:pver+1)-pf(:ncol,1:pver)
        1007  2323627992 :    hmn_zdp(:ncol,1:pver) = hmn_lev(:ncol,1:pver)*dp_lev(:ncol,1:pver)
        1008  2323627992 :    q_zdp(:ncol,1:pver) = q(:ncol,1:pver)*dp_lev(:ncol,1:pver)
        1009             :
        1010             :
        1011             : ! Mix profile over vertical length scale of 0.5*PBLH.
        1012             :
        1013    24969168 :    do i = 1,ncol ! Loop columns
        1014  1431901800 :       do k = pver,msg + 1,-1
        1015             :
        1016  1431901800 :          if (zi(i,k+1)<= parcel_dz(i)) then ! Has to be relative to near-surface layer center elevation
        1017   124110348 :             ipar = k
        1018             :
        1019   124110348 :             if (k == pver) then ! Always at least the full depth of lowest model layer.
        1020             :                dp_zfrac = 1._kind_phys
        1021             :             else
        1022             :                ! Fraction of grid cell depth (mostly 1, except when parcel_ztop is in between levels.
        1023   100636548 :                dp_zfrac =  min(1._kind_phys,(parcel_dz(i)-zi(i,k+1))/(zi(i,k)-zi(i,k+1)))
        1024             :             end if
        1025             :
        1026   124110348 :             parcel_hdp(i) = parcel_hdp(i)+hmn_zdp(i,k)*dp_zfrac ! Sum parcel profile up to a certain level.
        1027   124110348 :             parcel_qdp(i) = parcel_qdp(i)+q_zdp(i,k)*dp_zfrac ! Sum parcel profile up to a certain level.
        1028   124110348 :             parcel_dp(i)  = parcel_dp(i)+dp_lev(i,k)*dp_zfrac ! SUM dp's for weighting of parcel_hdp
        1029             :
        1030             :          end if
        1031             :       end do
        1032    23473800 :       hpar(i) = parcel_hdp(i)/parcel_dp(i)
        1033    23473800 :       qpar(i) = parcel_qdp(i)/parcel_dp(i)
        1034    24969168 :       mx(i) = ipar
        1035             :    end do
        1036             :
        1037             : else ! Default method finding level of MSE maximum (nlev sensitive though)
        1038             :     !
        1039             :     ! set "launching" level(mx) to be at maximum moist static energy.
        1040             :     ! search for this level stops at planetary boundary layer top.
        1041             :     !
        1042           0 :     do k = pver,msg + 1,-1
        1043           0 :        do i = 1,ncol
        1044           0 :           hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
        1045           0 :           if (k >= nint(pblt(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then
        1046           0 :              hmax(i) = hmn(i)
        1047           0 :              mx(i) = k
        1048             :           end if
        1049             :        end do
        1050             :     end do
        1051             :
        1052             : end if ! Default method of determining parcel launch properties.
        1053             :
        1054             :
        1055             :
        1056             :
        1057             :
        1058             : ! LCL dilute calculation - initialize to mx(i)
        1059             : ! Determine lcl in parcel_dilute and get pl,tl after parcel_dilute
        1060             : ! Original code actually sets LCL as level above wher condensate forms.
        1061             : ! Therefore in parcel_dilute lcl(i) will be at first level where qsmix < qtmix.
        1062             :
        1063     1495368 : if (lparcel_pbl) then
        1064             :
        1065             : ! For parcel dilute need to invert hpar and qpar.
        1066             : ! Now need to supply ql(i) as it is mixed parcel version, just q(i,max(i)) in default
        1067             :
        1068    24969168 :    do i = 1,ncol             ! Initialise LCL variables.
        1069    23473800 :       lcl(i) = mx(i)
        1070    23473800 :       tl(i) = (hpar(i)-rl*qpar(i)-grav*parcel_ztop(i))/cp
        1071    23473800 :       ql(i) = qpar(i)
        1072    24969168 :       pl(i) = p(i,mx(i))
        1073             :    end do
        1074             :
        1075             : else
        1076             :
        1077           0 :    do i = 1,ncol
        1078           0 :       lcl(i) = mx(i)
        1079           0 :       tl(i) = t(i,mx(i))
        1080           0 :       ql(i) = q(i,mx(i))
        1081           0 :       pl(i) = p(i,mx(i))
        1082             :    end do
        1083             :
        1084             : end if ! Mixed parcel properties
        1085             :
        1086             :
        1087             :
        1088             : !
        1089             : ! main buoyancy calculation.
        1090             : !
        1091             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        1092             : !!! DILUTE PLUME CALCULATION USING ENTRAINING PLUME !!!
        1093             : !!!   RBN 9/9/04   !!!
        1094             :
        1095             :    call parcel_dilute(ncol, pver, cpliq, cpwv, rh2o, latice, msg, mx, p, t, q, &
        1096             :    tpert, tp, tpv, qstp, pl, tl, ql, lcl, &
        1097     1495368 :    org, landfrac, errmsg, errflg)
        1098             :
        1099             :
        1100             : ! If lcl is above the nominal level of non-divergence (600 mbs),
        1101             : ! no deep convection is permitted (ensuing calculations
        1102             : ! skipped and cape retains initialized value of zero).
        1103             : !
        1104    24969168 :    do i = 1,ncol
        1105    24969168 :       plge600(i) = pl(i).ge.600._kind_phys ! Just change to always allow buoy calculation.
        1106             :    end do
        1107             :
        1108             : !
        1109             : ! Main buoyancy calculation.
        1110             : !
        1111    91217448 :    do k = pver,msg + 1,-1
        1112  1499645448 :       do i=1,ncol
        1113  1498150080 :          if (k <= mx(i) .and. plge600(i)) then   ! Define buoy from launch level to cloud top.
        1114  1245471862 :             tv(i,k) = t(i,k)* (1._kind_phys+1.608_kind_phys*q(i,k))/ (1._kind_phys+q(i,k))
        1115  1245471862 :             buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add  ! +0.5K or not?
        1116             :          else
        1117   162956138 :             qstp(i,k) = q(i,k)
        1118   162956138 :             tp(i,k)   = t(i,k)
        1119   162956138 :             tpv(i,k)  = tv(i,k)
        1120             :          endif
        1121             :       end do
        1122             :    end do
        1123             :
        1124             :
        1125             :
        1126             : !-------------------------------------------------------------------------------
        1127             :
        1128             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        1129             :
        1130             :
        1131             :
        1132             : !
        1133    89722080 :    do k = msg + 2,pver
        1134  1474676280 :       do i = 1,ncol
        1135  1473180912 :          if (k < lcl(i) .and. plge600(i)) then
        1136  1119120109 :             if (buoy(i,k+1) > 0._kind_phys .and. buoy(i,k) <= 0._kind_phys) then
        1137    22839231 :                knt(i) = min(num_cin,knt(i) + 1)
        1138    22839231 :                lelten(i,knt(i)) = k
        1139             :             end if
        1140             :          end if
        1141             :       end do
        1142             :    end do
        1143             : !
        1144             : ! calculate convective available potential energy (cape).
        1145             : !
        1146     2990736 :    do n = 1,num_cin
        1147    92712816 :       do k = msg + 1,pver
        1148  1499645448 :          do i = 1,ncol
        1149  1498150080 :             if (plge600(i) .and. k <= mx(i) .and. k > lelten(i,n)) then
        1150   179903275 :                capeten(i,n) = capeten(i,n) + rd*buoy(i,k)*log(pf(i,k+1)/pf(i,k))
        1151             :             end if
        1152             :          end do
        1153             :       end do
        1154             :    end do
        1155             : !
        1156             : ! find maximum cape from all possible tentative capes from
        1157             : ! one sounding,
        1158             : ! and use it as the final cape, april 26, 1995
        1159             : !
        1160     2990736 :    do n = 1,num_cin
        1161    26464536 :       do i = 1,ncol
        1162    24969168 :          if (capeten(i,n) > cape(i)) then
        1163    16731523 :             cape(i) = capeten(i,n)
        1164    16731523 :             lel(i) = lelten(i,n)
        1165             :          end if
        1166             :       end do
        1167             :    end do
        1168             : !
        1169             : ! put lower bound on cape for diagnostic purposes.
        1170             : !
        1171    24969168 :    do i = 1,ncol
        1172    24969168 :       cape(i) = max(cape(i), 0._kind_phys)
        1173             :    end do
        1174             : !
        1175     1495368 :    return
        1176             : end subroutine buoyan_dilute
        1177             :
        1178     1495368 : subroutine parcel_dilute (ncol, pver, cpliq, cpwv, rh2o, latice, msg, klaunch, p, t, q, &
        1179     2990736 :   tpert, tp, tpv, qstp, pl, tl, ql, lcl, &
        1180     1495368 :   org, landfrac,errmsg,errflg)
        1181             :
        1182             : ! Routine  to determine
        1183             : !   1. Tp   - Parcel temperature
        1184             : !   2. qstp - Saturated mixing ratio at the parcel temperature.
        1185             :
        1186             : !--------------------
        1187             : implicit none
        1188             : !--------------------
        1189             :
        1190             : integer, intent(in) :: ncol
        1191             : integer, intent(in) :: pver
        1192             : real(kind_phys), intent(in) :: cpliq
        1193             : real(kind_phys), intent(in) :: cpwv
        1194             : real(kind_phys), intent(in) :: rh2o
        1195             : real(kind_phys), intent(in) :: latice
        1196             : integer, intent(in) :: msg
        1197             :
        1198             : integer, intent(in), dimension(ncol) :: klaunch(ncol)
        1199             :
        1200             : real(kind_phys), intent(in), dimension(ncol,pver) :: p
        1201             : real(kind_phys), intent(in), dimension(ncol,pver) :: t
        1202             : real(kind_phys), intent(in), dimension(ncol,pver) :: q
        1203             : real(kind_phys), intent(in), dimension(ncol) :: tpert ! PBL temperature perturbation.
        1204             :
        1205             : real(kind_phys), intent(inout), dimension(ncol,pver) :: tp    ! Parcel temp.
        1206             : real(kind_phys), intent(inout), dimension(ncol,pver) :: qstp  ! Parcel water vapour (sat value above lcl).
        1207             : real(kind_phys), intent(inout), dimension(ncol) :: tl         ! Actual temp of LCL.
        1208             : real(kind_phys), intent(inout), dimension(ncol) :: ql ! Actual humidity of LCL
        1209             : real(kind_phys), intent(inout), dimension(ncol) :: pl          ! Actual pressure of LCL.
        1210             :
        1211             : integer, intent(inout), dimension(ncol) :: lcl ! Lifting condesation level (first model level with saturation).
        1212             :
        1213             : real(kind_phys), intent(out), dimension(ncol,pver) :: tpv   ! Define tpv within this routine.
        1214             :
        1215             : character(len=512), intent(out)      :: errmsg
        1216             : integer, intent(out)                 :: errflg
        1217             :
        1218             :
        1219             :
        1220             : real(kind_phys), dimension(:,:) :: org
        1221             : real(kind_phys), intent(in), dimension(ncol) :: landfrac
        1222             : !--------------------
        1223             :
        1224             : ! Have to be careful as s is also dry static energy.
        1225             :
        1226             :
        1227             : ! If we are to retain the fact that CAM loops over grid-points in the internal
        1228             : ! loop then we need to dimension sp,atp,mp,xsh2o with ncol.
        1229             :
        1230             :
        1231     2990736 : real(kind_phys) tmix(ncol,pver)        ! Tempertaure of the entraining parcel.
        1232     2990736 : real(kind_phys) qtmix(ncol,pver)       ! Total water of the entraining parcel.
        1233     2990736 : real(kind_phys) qsmix(ncol,pver)       ! Saturated mixing ratio at the tmix.
        1234     2990736 : real(kind_phys) smix(ncol,pver)        ! Entropy of the entraining parcel.
        1235     2990736 : real(kind_phys) xsh2o(ncol,pver)       ! Precipitate lost from parcel.
        1236     2990736 : real(kind_phys) ds_xsh2o(ncol,pver)    ! Entropy change due to loss of condensate.
        1237     2990736 : real(kind_phys) ds_freeze(ncol,pver)   ! Entropy change sue to freezing of precip.
        1238             : real(kind_phys) dmpdz2d(ncol,pver)     ! variable detrainment rate
        1239             :
        1240     2990736 : real(kind_phys) mp(ncol)    ! Parcel mass flux.
        1241     2990736 : real(kind_phys) qtp(ncol)   ! Parcel total water.
        1242     2990736 : real(kind_phys) sp(ncol)    ! Parcel entropy.
        1243             :
        1244     2990736 : real(kind_phys) sp0(ncol)    ! Parcel launch entropy.
        1245     2990736 : real(kind_phys) qtp0(ncol)   ! Parcel launch total water.
        1246     1495368 : real(kind_phys) mp0(ncol)    ! Parcel launch relative mass flux.
        1247             :
        1248             : real(kind_phys) lwmax      ! Maximum condesate that can be held in cloud before rainout.
        1249             : real(kind_phys) dmpdp      ! Parcel fractional mass entrainment rate (/mb).
        1250             : real(kind_phys) dmpdz      ! Parcel fractional mass entrainment rate (/m)
        1251             : real(kind_phys) dpdz,dzdp  ! Hydrstatic relation and inverse of.
        1252             : real(kind_phys) senv       ! Environmental entropy at each grid point.
        1253             : real(kind_phys) qtenv      ! Environmental total water "   "   ".
        1254             : real(kind_phys) penv       ! Environmental total pressure "   "   ".
        1255             : real(kind_phys) tenv       ! Environmental total temperature "   "   ".
        1256             : real(kind_phys) new_s      ! Hold value for entropy after condensation/freezing adjustments.
        1257             : real(kind_phys) new_q      ! Hold value for total water after condensation/freezing adjustments.
        1258             : real(kind_phys) dp         ! Layer thickness (center to center)
        1259             : real(kind_phys) tfguess    ! First guess for entropy inversion - crucial for efficiency!
        1260             : real(kind_phys) tscool     ! Super cooled temperature offset (in degC) (eg -35).
        1261             :
        1262             : real(kind_phys) qxsk, qxskp1        ! LCL excess water (k, k+1)
        1263             : real(kind_phys) dsdp, dqtdp, dqxsdp ! LCL s, qt, p gradients (k, k+1)
        1264             : real(kind_phys) slcl,qtlcl,qslcl    ! LCL s, qt, qs values.
        1265             : real(kind_phys) org2rkm, org2Tpert
        1266             : real(kind_phys) dmpdz_lnd, dmpdz_mask
        1267             :
        1268             : integer rcall       ! Number of ientropy call for errors recording
        1269             : integer nit_lheat     ! Number of iterations for condensation/freezing loop.
        1270             : integer i,k,ii   ! Loop counters.
        1271             :
        1272             : !======================================================================
        1273             : !    SUMMARY
        1274             : !
        1275             : !  9/9/04 - Assumes parcel is initiated from level of maxh (klaunch)
        1276             : !           and entrains at each level with a specified entrainment rate.
        1277             : !
        1278             : ! 15/9/04 - Calculates lcl(i) based on k where qsmix is first < qtmix.
        1279             : !
        1280             : !======================================================================
        1281             : !
        1282             : ! Set some values that may be changed frequently.
        1283             : !
        1284             :
        1285             : if (zm_org) then
        1286             :    org2rkm = 10._kind_phys
        1287             :    org2Tpert = 0._kind_phys
        1288             : endif
        1289     1495368 : nit_lheat = 2 ! iterations for ds,dq changes from condensation freezing.
        1290     1495368 : dmpdz=dmpdz_param       ! Entrainment rate. (-ve for /m)
        1291     1495368 : dmpdz_lnd=-1.e-3_kind_phys
        1292     1495368 : lwmax = 1.e-3_kind_phys    ! Need to put formula in for this.
        1293     1495368 : tscool = 0.0_kind_phys   ! Temp at which water loading freezes in the cloud.
        1294             :
        1295  2323627992 : qtmix=0._kind_phys
        1296  2323627992 : smix=0._kind_phys
        1297             :
        1298     1495368 : qtenv = 0._kind_phys
        1299     1495368 : senv = 0._kind_phys
        1300     1495368 : tenv = 0._kind_phys
        1301     1495368 : penv = 0._kind_phys
        1302             :
        1303    24969168 : qtp0 = 0._kind_phys
        1304    24969168 : sp0  = 0._kind_phys
        1305    24969168 : mp0 = 0._kind_phys
        1306             :
        1307    24969168 : qtp = 0._kind_phys
        1308    24969168 : sp = 0._kind_phys
        1309    24969168 : mp = 0._kind_phys
        1310             :
        1311     1495368 : new_q = 0._kind_phys
        1312     1495368 : new_s = 0._kind_phys
        1313             :
        1314             : ! **** Begin loops ****
        1315             :
        1316    91217448 : do k = pver, msg+1, -1
        1317  1499645448 :    do i=1,ncol
        1318             :
        1319             : ! Initialize parcel values at launch level.
        1320             :
        1321  1408428000 :       if (k == klaunch(i)) then
        1322             :
        1323    23473800 :          if (lparcel_pbl) then ! Modifcations to parcel properties if lparcel_pbl set.
        1324             :
        1325    23473800 :             qtp0(i) = ql(i)     ! Parcel launch q (PBL mixed value).
        1326    23473800 :             sp0(i)  = entropy(tl(i),pl(i),qtp0(i),cpliq,cpwv,rh2o) ! Parcel launch entropy could be a mixed parcel.
        1327             :
        1328             :          else
        1329             :
        1330           0 :             qtp0(i) = q(i,k)    ! Parcel launch total water (assuming subsaturated)
        1331           0 :             sp0(i)  = entropy(t(i,k),p(i,k),qtp0(i),cpliq,cpwv,rh2o) ! Parcel launch entropy.
        1332             :
        1333             :          end if
        1334             :
        1335    23473800 :          mp0(i)  = 1._kind_phys       ! Parcel launch relative mass (i.e. 1 parcel stays 1 parcel for dmpdp=0, undilute).
        1336    23473800 :          smix(i,k)  = sp0(i)
        1337    23473800 :          qtmix(i,k) = qtp0(i)
        1338    23473800 :          tfguess = t(i,k)
        1339    23473800 :          rcall = 1
        1340             :          call ientropy (rcall,i,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess,cpliq,cpwv,rh2o,errmsg,errflg)
        1341             :       end if
        1342             :
        1343             : ! Entraining levels
        1344             :
        1345  1498150080 :       if (k < klaunch(i)) then
        1346             :
        1347             : ! Set environmental values for this level.
        1348             :
        1349  1284317652 :          dp = (p(i,k)-p(i,k+1)) ! In -ve mb as p decreasing with height - difference between center of layers.
        1350  1284317652 :          qtenv = 0.5_kind_phys*(q(i,k)+q(i,k+1))         ! Total water of environment.
        1351  1284317652 :          tenv  = 0.5_kind_phys*(t(i,k)+t(i,k+1))
        1352  1284317652 :          penv  = 0.5_kind_phys*(p(i,k)+p(i,k+1))
        1353             :
        1354  1284317652 :          senv  = entropy(tenv,penv,qtenv,cpliq,cpwv,rh2o)  ! Entropy of environment.
        1355             :
        1356             : ! Determine fractional entrainment rate /pa given value /m.
        1357             :
        1358  1284317652 :          dpdz = -(penv*grav)/(rgas*tenv) ! in mb/m since  p in mb.
        1359  1284317652 :          dzdp = 1._kind_phys/dpdz                  ! in m/mb
        1360  1284317652 :          if (zm_org) then
        1361           0 :             dmpdz_mask = landfrac(i) * dmpdz_lnd + (1._kind_phys - landfrac(i)) * dmpdz
        1362           0 :             dmpdp = (dmpdz_mask/(1._kind_phys+org(i,k)*org2rkm))*dzdp              ! /mb Fractional entrainment
        1363             :          else
        1364  1284317652 :             dmpdp = dmpdz*dzdp
        1365             :          endif
        1366             :
        1367             : ! Sum entrainment to current level
        1368             : ! entrains q,s out of intervening dp layers, in which linear variation is assumed
        1369             : ! so really it entrains the mean of the 2 stored values.
        1370             :
        1371  1284317652 :          sp(i)  = sp(i)  - dmpdp*dp*senv
        1372  1284317652 :          qtp(i) = qtp(i) - dmpdp*dp*qtenv
        1373  1284317652 :          mp(i)  = mp(i)  - dmpdp*dp
        1374             :
        1375             : ! Entrain s and qt to next level.
        1376             :
        1377  1284317652 :          smix(i,k)  = (sp0(i)  +  sp(i)) / (mp0(i) + mp(i))
        1378  1284317652 :          qtmix(i,k) = (qtp0(i) + qtp(i)) / (mp0(i) + mp(i))
        1379             :
        1380             : ! Invert entropy from s and q to determine T and saturation-capped q of mixture.
        1381             : ! t(i,k) used as a first guess so that it converges faster.
        1382             :
        1383  1284317652 :          tfguess = tmix(i,k+1)
        1384  1284317652 :          rcall = 2
        1385             :          call ientropy(rcall,i,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess,cpliq,cpwv,rh2o,errmsg,errflg)
        1386             :
        1387             : !
        1388             : ! Determine if this is lcl of this column if qsmix <= qtmix.
        1389             : ! FIRST LEVEL where this happens on ascending.
        1390             :
        1391  1284317652 :          if (qsmix(i,k) <= qtmix(i,k) .and. qsmix(i,k+1) > qtmix(i,k+1)) then
        1392    21501361 :             lcl(i) = k
        1393    21501361 :             qxsk   = qtmix(i,k) - qsmix(i,k)
        1394    21501361 :             qxskp1 = qtmix(i,k+1) - qsmix(i,k+1)
        1395    21501361 :             dqxsdp = (qxsk - qxskp1)/dp
        1396    21501361 :             pl(i)  = p(i,k+1) - qxskp1/dqxsdp    ! pressure level of actual lcl.
        1397    21501361 :             dsdp   = (smix(i,k)  - smix(i,k+1))/dp
        1398    21501361 :             dqtdp  = (qtmix(i,k) - qtmix(i,k+1))/dp
        1399    21501361 :             slcl   = smix(i,k+1)  +  dsdp* (pl(i)-p(i,k+1))
        1400    21501361 :             qtlcl  = qtmix(i,k+1) +  dqtdp*(pl(i)-p(i,k+1))
        1401             :
        1402    21501361 :             tfguess = tmix(i,k)
        1403    21501361 :             rcall = 3
        1404    21501361 :             call ientropy (rcall,i,slcl,pl(i),qtlcl,tl(i),qslcl,tfguess,cpliq,cpwv,rh2o,errmsg,errflg)
        1405             :
        1406             : !            write(iulog,*)' '
        1407             : !            write(iulog,*)' p',p(i,k+1),pl(i),p(i,lcl(i))
        1408             : !            write(iulog,*)' t',tmix(i,k+1),tl(i),tmix(i,lcl(i))
        1409             : !            write(iulog,*)' s',smix(i,k+1),slcl,smix(i,lcl(i))
        1410             : !            write(iulog,*)'qt',qtmix(i,k+1),qtlcl,qtmix(i,lcl(i))
        1411             : !            write(iulog,*)'qs',qsmix(i,k+1),qslcl,qsmix(i,lcl(i))
        1412             :
        1413             :          endif
        1414             : !
        1415             :       end if !  k < klaunch
        1416             :
        1417             :
        1418             :    end do ! Levels loop
        1419             : end do ! Columns loop
        1420             :
        1421             : !!!!!!!!!!!!!!!!!!!!!!!!!!END ENTRAINMENT LOOP!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        1422             :
        1423             : !! Could stop now and test with this as it will provide some estimate of buoyancy
        1424             : !! without the effects of freezing/condensation taken into account for tmix.
        1425             :
        1426             : !! So we now have a profile of entropy and total water of the entraining parcel
        1427             : !! Varying with height from the launch level klaunch parcel=environment. To the
        1428             : !! top allowed level for the existence of convection.
        1429             :
        1430             : !! Now we have to adjust these values such that the water held in vaopor is < or
        1431             : !! = to qsmix. Therefore, we assume that the cloud holds a certain amount of
        1432             : !! condensate (lwmax) and the rest is rained out (xsh2o). This, obviously
        1433             : !! provides latent heating to the mixed parcel and so this has to be added back
        1434             : !! to it. But does this also increase qsmix as well? Also freezing processes
        1435             :
        1436             :
        1437  2323627992 : xsh2o = 0._kind_phys
        1438  2323627992 : ds_xsh2o = 0._kind_phys
        1439  2323627992 : ds_freeze = 0._kind_phys
        1440             :
        1441             : !!!!!!!!!!!!!!!!!!!!!!!!!PRECIPITATION/FREEZING LOOP!!!!!!!!!!!!!!!!!!!!!!!!!!
        1442             : !! Iterate solution twice for accuracy
        1443             :
        1444             :
        1445             :
        1446    91217448 : do k = pver, msg+1, -1
        1447  1499645448 :    do i=1,ncol
        1448             :
        1449             : ! Initialize variables at k=klaunch
        1450             :
        1451  1408428000 :       if (k == klaunch(i)) then
        1452             :
        1453             : ! Set parcel values at launch level assume no liquid water.
        1454             :
        1455    23473800 :          tp(i,k)    = tmix(i,k)
        1456    23473800 :          qstp(i,k)  = q(i,k)
        1457    23473800 :          if (zm_org) then
        1458           0 :             tpv(i,k)   =  (tp(i,k) + (org2Tpert*org(i,k)+tpert(i))) * (1._kind_phys+1.608_kind_phys*qstp(i,k)) / (1._kind_phys+qstp(i,k))
        1459             :          else
        1460    23473800 :             tpv(i,k)   =  (tp(i,k) + tpert(i)) * (1._kind_phys+1.608_kind_phys*qstp(i,k)) / (1._kind_phys+qstp(i,k))
        1461             :          endif
        1462             :
        1463             :       end if
        1464             :
        1465  1498150080 :       if (k < klaunch(i)) then
        1466             :
        1467             : ! Initiaite loop if switch(2) = .T. - RBN:DILUTE - TAKEN OUT BUT COULD BE RETURNED LATER.
        1468             :
        1469             : ! Iterate nit_lheat times for s,qt changes.
        1470             :
        1471  3852952956 :          do ii=0,nit_lheat-1
        1472             :
        1473             : ! Rain (xsh2o) is excess condensate, bar LWMAX (Accumulated loss from qtmix).
        1474             :
        1475  2568635304 :             xsh2o(i,k) = max (0._kind_phys, qtmix(i,k) - qsmix(i,k) - lwmax)
        1476             :
        1477             : ! Contribution to ds from precip loss of condensate (Accumulated change from smix).(-ve)
        1478             :
        1479  2568635304 :             ds_xsh2o(i,k) = ds_xsh2o(i,k+1) - cpliq * log (tmix(i,k)/tfreez) * max(0._kind_phys,(xsh2o(i,k)-xsh2o(i,k+1)))
        1480             : !
        1481             : ! Entropy of freezing: latice times amount of water involved divided by T.
        1482             : !
        1483             :
        1484  2568635304 :             if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) == 0._kind_phys) then ! One off freezing of condensate.
        1485   106642733 :                ds_freeze(i,k) = (latice/tmix(i,k)) * max(0._kind_phys,qtmix(i,k)-qsmix(i,k)-xsh2o(i,k)) ! Gain of LH
        1486             :             end if
        1487             :
        1488  2568635304 :             if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) /= 0._kind_phys) then ! Continual freezing of additional condensate.
        1489  2001007784 :                ds_freeze(i,k) = ds_freeze(i,k+1)+(latice/tmix(i,k)) * max(0._kind_phys,(qsmix(i,k+1)-qsmix(i,k)))
        1490             :             end if
        1491             :
        1492             : ! Adjust entropy and accordingly to sum of ds (be careful of signs).
        1493             :
        1494  2568635304 :             new_s = smix(i,k) + ds_xsh2o(i,k) + ds_freeze(i,k)
        1495             :
        1496             : ! Adjust liquid water and accordingly to xsh2o.
        1497             :
        1498  2568635304 :             new_q = qtmix(i,k) - xsh2o(i,k)
        1499             :
        1500             : ! Invert entropy to get updated Tmix and qsmix of parcel.
        1501             :
        1502  2568635304 :             tfguess = tmix(i,k)
        1503  2568635304 :             rcall =4
        1504  3852952956 :             call ientropy (rcall,i,new_s, p(i,k), new_q, tmix(i,k), qsmix(i,k), tfguess,cpliq,cpwv,rh2o,errmsg,errflg)
        1505             :
        1506             :          end do  ! Iteration loop for freezing processes.
        1507             :
        1508             : ! tp  - Parcel temp is temp of mixture.
        1509             : ! tpv - Parcel v. temp should be density temp with new_q total water.
        1510             :
        1511  1284317652 :          tp(i,k)    = tmix(i,k)
        1512             :
        1513             : ! tpv = tprho in the presence of condensate (i.e. when new_q > qsmix)
        1514             :
        1515  1284317652 :          if (new_q > qsmix(i,k)) then  ! Super-saturated so condensate present - reduces buoyancy.
        1516  1206901008 :             qstp(i,k) = qsmix(i,k)
        1517             :          else                          ! Just saturated/sub-saturated - no condensate virtual effects.
        1518    77416644 :             qstp(i,k) = new_q
        1519             :          end if
        1520             :
        1521  1284317652 :          if (zm_org) then
        1522           0 :             tpv(i,k) = (tp(i,k)+(org2Tpert*org(i,k)+tpert(i)))* (1._kind_phys+1.608_kind_phys*qstp(i,k)) / (1._kind_phys+ new_q)
        1523             :          else
        1524  1284317652 :             tpv(i,k) = (tp(i,k)+tpert(i))* (1._kind_phys+1.608_kind_phys*qstp(i,k)) / (1._kind_phys+ new_q)
        1525             :          endif
        1526             :
        1527             :       end if ! k < klaunch
        1528             :
        1529             :    end do ! Loop for columns
        1530             :
        1531             : end do  ! Loop for vertical levels.
        1532             :
        1533             :
        1534     1495368 : return
        1535             : end subroutine parcel_dilute
        1536             :
        1537             : !-----------------------------------------------------------------------------------------
        1538 24718872406 : real(kind_phys) function entropy(TK,p,qtot,cpliq,cpwv,rh2o)
        1539             : !-----------------------------------------------------------------------------------------
        1540             : !
        1541             : ! TK(K),p(mb),qtot(kg/kg)
        1542             : ! from Raymond and Blyth 1992
        1543             : !
        1544             :      real(kind_phys), intent(in) :: p,qtot,TK
        1545             :      real(kind_phys), intent(in) :: cpliq
        1546             :      real(kind_phys), intent(in) :: cpwv
        1547             :      real(kind_phys), intent(in) :: rh2o
        1548             :
        1549             :      real(kind_phys) :: qv,qst,e,est,L
        1550             :      real(kind_phys), parameter :: pref = 1000._kind_phys
        1551             :
        1552 24718872406 : L = rl - (cpliq - cpwv)*(TK-tfreez)         ! T IN CENTIGRADE
        1553             :
        1554 24718872406 : call qsat_hPa(TK, p, est, qst)
        1555             :
        1556 24718872406 : qv = min(qtot,qst)                         ! Partition qtot into vapor part only.
        1557 24718872406 : e = qv*p / (eps1 +qv)
        1558             :
        1559             : entropy = (cpres + qtot*cpliq)*log( TK/tfreez) - rgas*log( (p-e)/pref ) + &
        1560 24718872406 :         L*qv/TK - qv*rh2o*log(qv/qst)
        1561             :
        1562 24718872406 : end FUNCTION entropy
        1563             :
        1564             : !
        1565             : !-----------------------------------------------------------------------------------------
        1566  3897928117 : SUBROUTINE ientropy (rcall,icol,s,p,qt,T,qst,Tfg,cpliq,cpwv,rh2o,errmsg,errflg)
        1567             : !-----------------------------------------------------------------------------------------
        1568             : !
        1569             : ! p(mb), Tfg/T(K), qt/qv(kg/kg), s(J/kg).
        1570             : ! Inverts entropy, pressure and total water qt
        1571             : ! for T and saturated vapor mixing ratio
        1572             : !
        1573             :
        1574             : ! CACNOTE - Remove this when pass in lat/lon or pass out lchnk,icol (Note: lchnk will not exist in CAM-SIMA)
        1575             : !  use phys_grid, only: get_rlon_p, get_rlat_p
        1576             :
        1577             :   integer, intent(in) :: icol, rcall
        1578             :   real(kind_phys), intent(in)  :: s, p, Tfg, qt
        1579             :   real(kind_phys), intent(in) :: cpliq
        1580             :   real(kind_phys), intent(in) :: cpwv
        1581             :   real(kind_phys), intent(in) :: rh2o
        1582             :   real(kind_phys), intent(out) :: qst, T
        1583             :   character(len=512), intent(out)      :: errmsg
        1584             :   integer, intent(out)                 :: errflg
        1585             :
        1586             :   real(kind_phys) :: est, this_lat,this_lon
        1587             :   real(kind_phys) :: a,b,c,d,ebr,fa,fb,fc,pbr,qbr,rbr,sbr,tol1,xm,tol
        1588             :   integer :: i
        1589             :
        1590             :   logical :: converged
        1591             :
        1592             :   ! Max number of iteration loops.
        1593             :   integer, parameter :: LOOPMAX = 100
        1594             :   real(kind_phys), parameter :: EPS = 3.e-8_kind_phys
        1595             :
        1596  3897928117 :   converged = .false.
        1597             :
        1598             :   ! Invert the entropy equation -- use Brent's method
        1599             :   ! Brent, R. P. Ch. 3-4 in Algorithms for Minimization Without Derivatives. Englewood Cliffs, NJ: Prentice-Hall, 1973.
        1600             :
        1601  3897928117 :   T = Tfg                  ! Better first guess based on Tprofile from conv.
        1602             :
        1603  3897928117 :   a = Tfg-10    !low bracket
        1604  3897928117 :   b = Tfg+10    !high bracket
        1605             :
        1606  3897928117 :   fa = entropy(a, p, qt,cpliq,cpwv,rh2o) - s
        1607  3897928117 :   fb = entropy(b, p, qt,cpliq,cpwv,rh2o) - s
        1608             :
        1609  3897928117 :   c=b
        1610  3897928117 :   fc=fb
        1611  3897928117 :   tol=0.001_kind_phys
        1612             :
        1613 19513152837 :   converge: do i=0, LOOPMAX
        1614 19513152837 :      if ((fb > 0.0_kind_phys .and. fc > 0.0_kind_phys) .or. &
        1615             :           (fb < 0.0_kind_phys .and. fc < 0.0_kind_phys)) then
        1616 12782668392 :         c=a
        1617 12782668392 :         fc=fa
        1618 12782668392 :         d=b-a
        1619 12782668392 :         ebr=d
        1620             :      end if
        1621 19513152837 :      if (abs(fc) < abs(fb)) then
        1622  4641732623 :         a=b
        1623  4641732623 :         b=c
        1624  4641732623 :         c=a
        1625  4641732623 :         fa=fb
        1626  4641732623 :         fb=fc
        1627  4641732623 :         fc=fa
        1628             :      end if
        1629             :
        1630 19513152837 :      tol1=2.0_kind_phys*EPS*abs(b)+0.5_kind_phys*tol
        1631 19513152837 :      xm=0.5_kind_phys*(c-b)
        1632 19513152837 :      converged = (abs(xm) <= tol1 .or. fb == 0.0_kind_phys)
        1633 19513152837 :      if (converged) exit converge
        1634             :
        1635 15615224720 :      if (abs(ebr) >= tol1 .and. abs(fa) > abs(fb)) then
        1636 15615224720 :         sbr=fb/fa
        1637 15615224720 :         if (a == c) then
        1638  8884834670 :            pbr=2.0_kind_phys*xm*sbr
        1639  8884834670 :            qbr=1.0_kind_phys-sbr
        1640             :         else
        1641  6730390050 :            qbr=fa/fc
        1642  6730390050 :            rbr=fb/fc
        1643  6730390050 :            pbr=sbr*(2.0_kind_phys*xm*qbr*(qbr-rbr)-(b-a)*(rbr-1.0_kind_phys))
        1644  6730390050 :            qbr=(qbr-1.0_kind_phys)*(rbr-1.0_kind_phys)*(sbr-1.0_kind_phys)
        1645             :         end if
        1646 15615224720 :         if (pbr > 0.0_kind_phys) qbr=-qbr
        1647 15615224720 :         pbr=abs(pbr)
        1648 15615224720 :         if (2.0_kind_phys*pbr  <  min(3.0_kind_phys*xm*qbr-abs(tol1*qbr),abs(ebr*qbr))) then
        1649 15613826358 :            ebr=d
        1650 15613826358 :            d=pbr/qbr
        1651             :         else
        1652             :            d=xm
        1653             :            ebr=d
        1654             :         end if
        1655             :      else
        1656             :         d=xm
        1657             :         ebr=d
        1658             :      end if
        1659 15615224720 :      a=b
        1660 15615224720 :      fa=fb
        1661 15615224720 :      b=b+merge(d,sign(tol1,xm), abs(d) > tol1 )
        1662             :
        1663 19513152837 :      fb = entropy(b, p, qt,cpliq,cpwv,rh2o) - s
        1664             :
        1665             :   end do converge
        1666             :
        1667  3897928117 :   T = b
        1668  3897928117 :   call qsat_hPa(T, p, est, qst)
        1669             :
        1670  3897928117 :   if (.not. converged) then
        1671             : !CACNOTE - Revisit this with Jesse
        1672             : !     this_lat = get_rlat_p(lchnk, icol)*57.296_kind_phys
        1673             : !     this_lon = get_rlon_p(lchnk, icol)*57.296_kind_phys
        1674             : !     write(errmsg,100) 'ZM_CONV: IENTROPY. Details: call#,lchnk,icol= ',rcall,lchnk,icol, &
        1675             : !          ' lat: ',this_lat,' lon: ',this_lon, &
        1676             : !          ' P(mb)= ', p, ' Tfg(K)= ', Tfg, ' qt(g/kg) = ', 1000._kind_phys*qt, &
        1677             : !          ' qst(g/kg) = ', 1000._kind_phys*qst,', s(J/kg) = ',s
        1678           0 :      errflg=1
        1679             :   end if
        1680             :
        1681             : 100 format (A,I1,I4,I4,7(A,F6.2))
        1682             :
        1683  3897928117 : end SUBROUTINE ientropy
        1684             :
        1685      578550 : subroutine cldprp(ncol   ,pver    ,pverp   ,cpliq   , &
        1686             :                   latice  ,cpwv    ,rh2o    ,&
        1687      578550 :                   q       ,t       ,u       ,v       ,p       , &
        1688      578550 :                   z       ,s       ,mu      ,eu      ,du      , &
        1689      578550 :                   md      ,ed      ,sd      ,qd      ,mc      , &
        1690      578550 :                   qu      ,su      ,zf      ,qst     ,hmn     , &
        1691      578550 :                   hsat    ,shat    ,ql      , &
        1692      578550 :                   cmeg    ,jb      ,lel     ,jt      ,jlcl    , &
        1693      578550 :                   mx      ,j0      ,jd      ,rl      ,il2g    , &
        1694             :                   rd      ,grav    ,cp      ,msg     , &
        1695      578550 :                   evp     ,cu      ,rprd    ,limcnv  ,landfrac, &
        1696      578550 :                   qcde     ,qhat  )
        1697             :
        1698             : !-----------------------------------------------------------------------
        1699             : !CACNOTE - fill in documentation
        1700             : !
        1701             : ! Purpose:
        1702             : ! <Say what the routine does>
        1703             : !
        1704             : ! Method:
        1705             : ! may 09/91 - guang jun zhang, m.lazare, n.mcfarlane.
        1706             : !             original version cldprop.
        1707             : !
        1708             : ! Author: See above, modified by P. Rasch
        1709             : ! This is contributed code not fully standardized by the CCM core group.
        1710             : !
        1711             : ! this code is very much rougher than virtually anything else in the CCM
        1712             : ! there are debug statements left strewn about and code segments disabled
        1713             : ! these are to facilitate future development. We expect to release a
        1714             : ! cleaner code in a future release
        1715             : !
        1716             : ! the documentation has been enhanced to the degree that we are able
        1717             : !
        1718             : !-----------------------------------------------------------------------
        1719             :
        1720             :    implicit none
        1721             :
        1722             : !------------------------------------------------------------------------------
        1723             : !
        1724             : ! Input arguments
        1725             : !
        1726             :    integer, intent(in) :: ncol
        1727             :    integer, intent(in) :: pver
        1728             :    integer, intent(in) :: pverp
        1729             :
        1730             :    real(kind_phys), intent(in) :: cpliq
        1731             :    real(kind_phys), intent(in) :: latice
        1732             :    real(kind_phys), intent(in) :: cpwv
        1733             :    real(kind_phys), intent(in) :: rh2o
        1734             :
        1735             :    real(kind_phys), intent(in) :: q(ncol,pver)         ! spec. humidity of env
        1736             :    real(kind_phys), intent(in) :: t(ncol,pver)         ! temp of env
        1737             :    real(kind_phys), intent(in) :: p(ncol,pver)         ! pressure of env
        1738             :    real(kind_phys), intent(in) :: z(ncol,pver)         ! height of env
        1739             :    real(kind_phys), intent(in) :: s(ncol,pver)         ! normalized dry static energy of env
        1740             :    real(kind_phys), intent(in) :: zf(ncol,pverp)       ! height of interfaces
        1741             :    real(kind_phys), intent(in) :: u(ncol,pver)         ! zonal velocity of env
        1742             :    real(kind_phys), intent(in) :: v(ncol,pver)         ! merid. velocity of env
        1743             :
        1744             :    real(kind_phys), intent(in) :: landfrac(ncol) ! RBN Landfrac
        1745             :
        1746             :    integer, intent(in) :: jb(ncol)              ! updraft base level
        1747             :    integer, intent(in) :: lel(ncol)             ! updraft launch level
        1748             :    integer, intent(in) :: mx(ncol)              ! updraft base level (same is jb)
        1749             :    integer, intent(out) :: jt(ncol)              ! updraft plume top
        1750             :    integer, intent(out) :: jlcl(ncol)            ! updraft lifting cond level
        1751             :    integer, intent(out) :: j0(ncol)              ! level where updraft begins detraining
        1752             :    integer, intent(out) :: jd(ncol)              ! level of downdraft
        1753             :    integer, intent(in) :: limcnv                 ! convection limiting level
        1754             :    integer, intent(in) :: il2g                   !CORE GROUP REMOVE
        1755             :    integer, intent(in) :: msg                    ! missing moisture vals (always 0)
        1756             :    real(kind_phys), intent(in) :: rl                    ! latent heat of vap
        1757             :    real(kind_phys), intent(in) :: shat(ncol,pver)      ! interface values of dry stat energy
        1758             :    real(kind_phys), intent(in) :: qhat(ncol,pver)      ! wg grid slice of upper interface mixing ratio.
        1759             :
        1760             : !
        1761             : ! output
        1762             : !
        1763             :
        1764             :    real(kind_phys), intent(out) :: rprd(ncol,pver)     ! rate of production of precip at that layer
        1765             :    real(kind_phys), intent(out) :: du(ncol,pver)       ! detrainement rate of updraft
        1766             :    real(kind_phys), intent(out) :: ed(ncol,pver)       ! entrainment rate of downdraft
        1767             :    real(kind_phys), intent(out) :: eu(ncol,pver)       ! entrainment rate of updraft
        1768             :    real(kind_phys), intent(out) :: hmn(ncol,pver)      ! moist stat energy of env
        1769             :    real(kind_phys), intent(out) :: hsat(ncol,pver)     ! sat moist stat energy of env
        1770             :    real(kind_phys), intent(out) :: mc(ncol,pver)       ! net mass flux
        1771             :    real(kind_phys), intent(out) :: md(ncol,pver)       ! downdraft mass flux
        1772             :    real(kind_phys), intent(out) :: mu(ncol,pver)       ! updraft mass flux
        1773             :    real(kind_phys), intent(out) :: qd(ncol,pver)       ! spec humidity of downdraft
        1774             :    real(kind_phys), intent(out) :: ql(ncol,pver)       ! liq water of updraft
        1775             :    real(kind_phys), intent(out) :: qst(ncol,pver)      ! saturation mixing ratio of env.
        1776             :    real(kind_phys), intent(out) :: qu(ncol,pver)       ! spec hum of updraft
        1777             :    real(kind_phys), intent(out) :: sd(ncol,pver)       ! normalized dry stat energy of downdraft
        1778             :    real(kind_phys), intent(out) :: su(ncol,pver)       ! normalized dry stat energy of updraft
        1779             :    real(kind_phys), intent(out) :: qcde(ncol,pver)     ! cloud water mixing ratio for detrainment (kg/kg)
        1780             :
        1781             :    real(kind_phys) rd                   ! gas constant for dry air
        1782             :    real(kind_phys) grav                 ! gravity
        1783             :    real(kind_phys) cp                   ! heat capacity of dry air
        1784             :
        1785             : !
        1786             : ! Local workspace
        1787             : !
        1788     1157100 :    real(kind_phys) gamma(ncol,pver)
        1789     1157100 :    real(kind_phys) dz(ncol,pver)
        1790     1157100 :    real(kind_phys) iprm(ncol,pver)
        1791     1157100 :    real(kind_phys) hu(ncol,pver)
        1792     1157100 :    real(kind_phys) hd(ncol,pver)
        1793     1157100 :    real(kind_phys) eps(ncol,pver)
        1794     1157100 :    real(kind_phys) f(ncol,pver)
        1795     1157100 :    real(kind_phys) k1(ncol,pver)
        1796     1157100 :    real(kind_phys) i2(ncol,pver)
        1797     1157100 :    real(kind_phys) ihat(ncol,pver)
        1798     1157100 :    real(kind_phys) i3(ncol,pver)
        1799     1157100 :    real(kind_phys) idag(ncol,pver)
        1800     1157100 :    real(kind_phys) i4(ncol,pver)
        1801     1157100 :    real(kind_phys) qsthat(ncol,pver)
        1802     1157100 :    real(kind_phys) hsthat(ncol,pver)
        1803     1157100 :    real(kind_phys) gamhat(ncol,pver)
        1804             :    real(kind_phys) cu(ncol,pver)
        1805             :    real(kind_phys) evp(ncol,pver)
        1806             :    real(kind_phys) cmeg(ncol,pver)
        1807     1157100 :    real(kind_phys) qds(ncol,pver)
        1808     1157100 :    real(kind_phys) c0mask(ncol)
        1809             :
        1810     1157100 :    real(kind_phys) hmin(ncol)
        1811     1157100 :    real(kind_phys) expdif(ncol)
        1812     1157100 :    real(kind_phys) expnum(ncol)
        1813     1157100 :    real(kind_phys) ftemp(ncol)
        1814     1157100 :    real(kind_phys) eps0(ncol)
        1815     1157100 :    real(kind_phys) rmue(ncol)
        1816     1157100 :    real(kind_phys) zuef(ncol)
        1817     1157100 :    real(kind_phys) zdef(ncol)
        1818     1157100 :    real(kind_phys) epsm(ncol)
        1819     1157100 :    real(kind_phys) ratmjb(ncol)
        1820     1157100 :    real(kind_phys) est(ncol)
        1821     1157100 :    real(kind_phys) totpcp(ncol)
        1822     1157100 :    real(kind_phys) totevp(ncol)
        1823     1157100 :    real(kind_phys) alfa(ncol)
        1824             :    real(kind_phys) ql1
        1825             :    real(kind_phys) tu
        1826             :    real(kind_phys) estu
        1827             :    real(kind_phys) qstu
        1828             :
        1829             :    real(kind_phys) small
        1830             :    real(kind_phys) mdt
        1831             :
        1832     1157100 :    real(kind_phys) fice(ncol,pver)        ! ice fraction in precip production
        1833     1157100 :    real(kind_phys) tug(ncol,pver)
        1834             :
        1835     1157100 :    real(kind_phys) tvuo(ncol,pver)        ! updraft virtual T w/o freezing heating
        1836     1157100 :    real(kind_phys) tvu(ncol,pver)         ! updraft virtual T with freezing heating
        1837     1157100 :    real(kind_phys) totfrz(ncol)
        1838     1157100 :    real(kind_phys) frz (ncol,pver)        ! rate of freezing
        1839     1157100 :    integer  jto(ncol)              ! updraft plume old top
        1840     1157100 :    integer  tmplel(ncol)
        1841             :
        1842             :    integer  iter, itnum
        1843             :    integer  m
        1844             :
        1845             :    integer khighest
        1846             :    integer klowest
        1847             :    integer kount
        1848             :    integer i,k
        1849             :
        1850     1157100 :    logical doit(ncol)
        1851     1157100 :    logical done(ncol)
        1852             : !
        1853             : !------------------------------------------------------------------------------
        1854             : !
        1855             :
        1856     4803370 :    do i = 1,il2g
        1857     4224820 :       ftemp(i) = 0._kind_phys
        1858     4224820 :       expnum(i) = 0._kind_phys
        1859     4224820 :       expdif(i) = 0._kind_phys
        1860     4803370 :       c0mask(i)  = c0_ocn * (1._kind_phys-landfrac(i)) +   c0_lnd * landfrac(i)
        1861             :    end do
        1862             : !
        1863             : !jr Change from msg+1 to 1 to prevent blowup
        1864             : !
        1865    54383700 :    do k = 1,pver
        1866   447291960 :       do i = 1,il2g
        1867   446713410 :          dz(i,k) = zf(i,k) - zf(i,k+1)
        1868             :       end do
        1869             :    end do
        1870             :
        1871             : !
        1872             : ! initialize many output and work variables to zero
        1873             : !
        1874    54383700 :    do k = 1,pver
        1875   447291960 :       do i = 1,il2g
        1876   392908260 :          k1(i,k) = 0._kind_phys
        1877   392908260 :          i2(i,k) = 0._kind_phys
        1878   392908260 :          i3(i,k) = 0._kind_phys
        1879   392908260 :          i4(i,k) = 0._kind_phys
        1880   392908260 :          mu(i,k) = 0._kind_phys
        1881   392908260 :          f(i,k) = 0._kind_phys
        1882   392908260 :          eps(i,k) = 0._kind_phys
        1883   392908260 :          eu(i,k) = 0._kind_phys
        1884   392908260 :          du(i,k) = 0._kind_phys
        1885   392908260 :          ql(i,k) = 0._kind_phys
        1886   392908260 :          cu(i,k) = 0._kind_phys
        1887   392908260 :          evp(i,k) = 0._kind_phys
        1888   392908260 :          cmeg(i,k) = 0._kind_phys
        1889   392908260 :          qds(i,k) = q(i,k)
        1890   392908260 :          md(i,k) = 0._kind_phys
        1891   392908260 :          ed(i,k) = 0._kind_phys
        1892   392908260 :          sd(i,k) = s(i,k)
        1893   392908260 :          qd(i,k) = q(i,k)
        1894   392908260 :          mc(i,k) = 0._kind_phys
        1895   392908260 :          qu(i,k) = q(i,k)
        1896   392908260 :          su(i,k) = s(i,k)
        1897   392908260 :          call qsat_hPa(t(i,k), p(i,k), est(i), qst(i,k))
        1898             :
        1899   392908260 :          if ( p(i,k)-est(i) <= 0._kind_phys ) then
        1900    33307924 :             qst(i,k) = 1.0_kind_phys
        1901             :          end if
        1902             :
        1903   392908260 :          gamma(i,k) = qst(i,k)*(1._kind_phys + qst(i,k)/eps1)*eps1*rl/(rd*t(i,k)**2)*rl/cp
        1904   392908260 :          hmn(i,k) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
        1905   392908260 :          hsat(i,k) = cp*t(i,k) + grav*z(i,k) + rl*qst(i,k)
        1906   392908260 :          hu(i,k) = hmn(i,k)
        1907   392908260 :          hd(i,k) = hmn(i,k)
        1908   392908260 :          rprd(i,k) = 0._kind_phys
        1909             :
        1910   392908260 :          fice(i,k) = 0._kind_phys
        1911   392908260 :          tug(i,k)  = 0._kind_phys
        1912   392908260 :          qcde(i,k)   = 0._kind_phys
        1913   392908260 :          tvuo(i,k) = (shat(i,k) - grav/cp*zf(i,k))*(1._kind_phys + 0.608_kind_phys*qhat(i,k))
        1914   392908260 :          tvu(i,k) = tvuo(i,k)
        1915   839621670 :          frz(i,k)  = 0._kind_phys
        1916             :
        1917             :       end do
        1918             :    end do
        1919             :
        1920             : !
        1921             : !jr Set to zero things which make this routine blow up
        1922             : !
        1923    19670700 :    do k=1,msg
        1924   159089760 :       do i=1,il2g
        1925   158511210 :          rprd(i,k) = 0._kind_phys
        1926             :       end do
        1927             :    end do
        1928             : !
        1929             : ! interpolate the layer values of qst, hsat and gamma to
        1930             : ! layer interfaces
        1931             : !
        1932    20249250 :    do k = 1, msg+1
        1933   163893130 :       do i = 1,il2g
        1934   143643880 :          hsthat(i,k) = hsat(i,k)
        1935   143643880 :          qsthat(i,k) = qst(i,k)
        1936   163314580 :          gamhat(i,k) = gamma(i,k)
        1937             :       end do
        1938             :    end do
        1939     4803370 :    do i = 1,il2g
        1940     4224820 :       totpcp(i) = 0._kind_phys
        1941     4803370 :       totevp(i) = 0._kind_phys
        1942             :    end do
        1943    34713000 :    do k = msg + 2,pver
        1944   283977380 :       do i = 1,il2g
        1945   249264380 :          if (abs(qst(i,k-1)-qst(i,k)) > 1.E-6_kind_phys) then
        1946   238641064 :             qsthat(i,k) = log(qst(i,k-1)/qst(i,k))*qst(i,k-1)*qst(i,k)/ (qst(i,k-1)-qst(i,k))
        1947             :          else
        1948    10623316 :             qsthat(i,k) = qst(i,k)
        1949             :          end if
        1950   249264380 :          hsthat(i,k) = cp*shat(i,k) + rl*qsthat(i,k)
        1951   283398830 :          if (abs(gamma(i,k-1)-gamma(i,k)) > 1.E-6_kind_phys) then
        1952             :             gamhat(i,k) = log(gamma(i,k-1)/gamma(i,k))*gamma(i,k-1)*gamma(i,k)/ &
        1953   249230721 :                                 (gamma(i,k-1)-gamma(i,k))
        1954             :          else
        1955       33659 :             gamhat(i,k) = gamma(i,k)
        1956             :          end if
        1957             :       end do
        1958             :    end do
        1959             : !
        1960             : ! initialize cloud top to highest plume top.
        1961             : !jr changed hard-wired 4 to limcnv+1 (not to exceed pver)
        1962             : !
        1963     9681835 :    jt(:) = pver
        1964     4803370 :    do i = 1,il2g
        1965     4224820 :       jt(i) = max(lel(i),limcnv+1)
        1966     4224820 :       jt(i) = min(jt(i),pver)
        1967     4224820 :       jd(i) = pver
        1968     4224820 :       jlcl(i) = lel(i)
        1969     4803370 :       hmin(i) = 1.E6_kind_phys
        1970             :    end do
        1971             : !
        1972             : ! find the level of minimum hsat, where detrainment starts
        1973             : !
        1974             :
        1975    35291550 :    do k = msg + 1,pver
        1976   288780750 :       do i = 1,il2g
        1977   288202200 :          if (hsat(i,k) <= hmin(i) .and. k >= jt(i) .and. k <= jb(i)) then
        1978    18825375 :             hmin(i) = hsat(i,k)
        1979    18825375 :             j0(i) = k
        1980             :          end if
        1981             :       end do
        1982             :    end do
        1983     4803370 :    do i = 1,il2g
        1984     4224820 :       j0(i) = min(j0(i),jb(i)-2)
        1985     4224820 :       j0(i) = max(j0(i),jt(i)+2)
        1986             : !
        1987             : ! Fix from Guang Zhang to address out of bounds array reference
        1988             : !
        1989     4803370 :       j0(i) = min(j0(i),pver)
        1990             :    end do
        1991             : !
        1992             : ! Initialize certain arrays inside cloud
        1993             : !
        1994    35291550 :    do k = msg + 1,pver
        1995   288780750 :       do i = 1,il2g
        1996   288202200 :          if (k >= jt(i) .and. k <= jb(i)) then
        1997    73449496 :             hu(i,k) = hmn(i,mx(i)) + cp*tiedke_add
        1998    73449496 :             su(i,k) = s(i,mx(i)) + tiedke_add
        1999             :          end if
        2000             :       end do
        2001             :    end do
        2002             : !
        2003             : ! *********************************************************
        2004             : ! compute taylor series for approximate eps(z) below
        2005             : ! *********************************************************
        2006             : !
        2007    34713000 :    do k = pver - 1,msg + 1,-1
        2008   283977380 :       do i = 1,il2g
        2009   283398830 :          if (k < jb(i) .and. k >= jt(i)) then
        2010    69224676 :             k1(i,k) = k1(i,k+1) + (hmn(i,mx(i))-hmn(i,k))*dz(i,k)
        2011    69224676 :             ihat(i,k) = 0.5_kind_phys* (k1(i,k+1)+k1(i,k))
        2012    69224676 :             i2(i,k) = i2(i,k+1) + ihat(i,k)*dz(i,k)
        2013    69224676 :             idag(i,k) = 0.5_kind_phys* (i2(i,k+1)+i2(i,k))
        2014    69224676 :             i3(i,k) = i3(i,k+1) + idag(i,k)*dz(i,k)
        2015    69224676 :             iprm(i,k) = 0.5_kind_phys* (i3(i,k+1)+i3(i,k))
        2016    69224676 :             i4(i,k) = i4(i,k+1) + iprm(i,k)*dz(i,k)
        2017             :          end if
        2018             :       end do
        2019             :    end do
        2020             : !
        2021             : ! re-initialize hmin array for ensuing calculation.
        2022             : !
        2023     4803370 :    do i = 1,il2g
        2024     4803370 :       hmin(i) = 1.E6_kind_phys
        2025             :    end do
        2026    35291550 :    do k = msg + 1,pver
        2027   288780750 :       do i = 1,il2g
        2028   288202200 :          if (k >= j0(i) .and. k <= jb(i) .and. hmn(i,k) <= hmin(i)) then
        2029     6423900 :             hmin(i) = hmn(i,k)
        2030     6423900 :             expdif(i) = hmn(i,mx(i)) - hmin(i)
        2031             :          end if
        2032             :       end do
        2033             :    end do
        2034             : !
        2035             : ! *********************************************************
        2036             : ! compute approximate eps(z) using above taylor series
        2037             : ! *********************************************************
        2038             : !
        2039    34713000 :    do k = msg + 2,pver
        2040   283977380 :       do i = 1,il2g
        2041   249264380 :          expnum(i) = 0._kind_phys
        2042   249264380 :          ftemp(i) = 0._kind_phys
        2043   249264380 :          if (k < jt(i) .or. k >= jb(i)) then
        2044   180039704 :             k1(i,k) = 0._kind_phys
        2045   180039704 :             expnum(i) = 0._kind_phys
        2046             :          else
        2047   207674028 :             expnum(i) = hmn(i,mx(i)) - (hsat(i,k-1)*(zf(i,k)-z(i,k)) + &
        2048   276898704 :                         hsat(i,k)* (z(i,k-1)-zf(i,k)))/(z(i,k-1)-z(i,k))
        2049             :          end if
        2050   249264380 :          if ((expdif(i) > 100._kind_phys .and. expnum(i) > 0._kind_phys) .and. &
        2051    34134450 :             k1(i,k) > expnum(i)*dz(i,k)) then
        2052    55843352 :             ftemp(i) = expnum(i)/k1(i,k)
        2053             :             f(i,k) = ftemp(i) + i2(i,k)/k1(i,k)*ftemp(i)**2 + &
        2054             :                      (2._kind_phys*i2(i,k)**2-k1(i,k)*i3(i,k))/k1(i,k)**2* &
        2055             :                      ftemp(i)**3 + (-5._kind_phys*k1(i,k)*i2(i,k)*i3(i,k)+ &
        2056             :                      5._kind_phys*i2(i,k)**3+k1(i,k)**2*i4(i,k))/ &
        2057    55843352 :                      k1(i,k)**3*ftemp(i)**4
        2058    55843352 :             f(i,k) = max(f(i,k),0._kind_phys)
        2059    55843352 :             f(i,k) = min(f(i,k),0.0002_kind_phys)
        2060             :          end if
        2061             :       end do
        2062             :    end do
        2063     4803370 :    do i = 1,il2g
        2064     4803370 :       if (j0(i) < jb(i)) then
        2065     4224820 :          if (f(i,j0(i)) < 1.E-6_kind_phys .and. f(i,j0(i)+1) > f(i,j0(i))) j0(i) = j0(i) + 1
        2066             :       end if
        2067             :    end do
        2068    34713000 :    do k = msg + 2,pver
        2069   283977380 :       do i = 1,il2g
        2070   283398830 :          if (k >= jt(i) .and. k <= j0(i)) then
        2071    23284594 :             f(i,k) = max(f(i,k),f(i,k-1))
        2072             :          end if
        2073             :       end do
        2074             :    end do
        2075     4803370 :    do i = 1,il2g
        2076     4224820 :       eps0(i) = f(i,j0(i))
        2077     4803370 :       eps(i,jb(i)) = eps0(i)
        2078             :    end do
        2079             : !
        2080             : ! This is set to match the Rasch and Kristjansson paper
        2081             : !
        2082    35291550 :    do k = pver,msg + 1,-1
        2083   288780750 :       do i = 1,il2g
        2084   288202200 :          if (k >= j0(i) .and. k <= jb(i)) then
        2085    54389722 :             eps(i,k) = f(i,j0(i))
        2086             :          end if
        2087             :       end do
        2088             :    end do
        2089    35291550 :    do k = pver,msg + 1,-1
        2090   288780750 :       do i = 1,il2g
        2091   288202200 :          if (k < j0(i) .and. k >= jt(i)) eps(i,k) = f(i,k)
        2092             :       end do
        2093             :    end do
        2094             :
        2095             :    itnum = 1
        2096     1157100 :    do iter=1, itnum
        2097             :
        2098             : !
        2099             : ! specify the updraft mass flux mu, entrainment eu, detrainment du
        2100             : ! and moist static energy hu.
        2101             : ! here and below mu, eu,du, md and ed are all normalized by mb
        2102             : !
        2103     4803370 :       do i = 1,il2g
        2104     4224820 :          if (eps0(i) > 0._kind_phys) then
        2105     4208860 :             mu(i,jb(i)) = 1._kind_phys
        2106     4208860 :             eu(i,jb(i)) = mu(i,jb(i))/dz(i,jb(i))
        2107             :          end if
        2108             :
        2109     4803370 :          tmplel(i) = jt(i)
        2110             :       end do
        2111    35291550 :       do k = pver,msg + 1,-1
        2112   288780750 :          do i = 1,il2g
        2113   288202200 :             if (eps0(i) > 0._kind_phys .and. (k >= tmplel(i) .and. k < jb(i))) then
        2114    68848238 :                zuef(i) = zf(i,k) - zf(i,jb(i))
        2115    68848238 :                rmue(i) = (1._kind_phys/eps0(i))* (exp(eps(i,k+1)*zuef(i))-1._kind_phys)/zuef(i)
        2116    68848238 :                mu(i,k) = (1._kind_phys/eps0(i))* (exp(eps(i,k  )*zuef(i))-1._kind_phys)/zuef(i)
        2117    68848238 :                eu(i,k) = (rmue(i)-mu(i,k+1))/dz(i,k)
        2118    68848238 :                du(i,k) = (rmue(i)-mu(i,k))/dz(i,k)
        2119             :             end if
        2120             :          end do
        2121             :       end do
        2122             :
        2123             :       khighest = pverp
        2124             :       klowest = 1
        2125     4803370 :       do i=1,il2g
        2126     4224820 :          khighest = min(khighest,lel(i))
        2127     4803370 :          klowest = max(klowest,jb(i))
        2128             :       end do
        2129    11472496 :       do k = klowest-1,khighest,-1
        2130    97097036 :          do i = 1,il2g
        2131    96518486 :             if (k <= jb(i)-1 .and. k >= lel(i) .and. eps0(i) > 0._kind_phys) then
        2132    68848238 :                if (mu(i,k) < 0.02_kind_phys) then
        2133     2248398 :                   hu(i,k) = hmn(i,k)
        2134     2248398 :                   mu(i,k) = 0._kind_phys
        2135     2248398 :                   eu(i,k) = 0._kind_phys
        2136     2248398 :                   du(i,k) = mu(i,k+1)/dz(i,k)
        2137             :                else
        2138    66599840 :                   hu(i,k) = mu(i,k+1)/mu(i,k)*hu(i,k+1) + &
        2139   133199680 :                             dz(i,k)/mu(i,k)* (eu(i,k)*hmn(i,k)- du(i,k)*hsat(i,k))
        2140             :                end if
        2141             :             end if
        2142             :          end do
        2143             :       end do
        2144             : !
        2145             : ! reset cloud top index beginning from two layers above the
        2146             : ! cloud base (i.e. if cloud is only one layer thick, top is not reset
        2147             : !
        2148     4803370 :       do i=1,il2g
        2149     4224820 :          doit(i) = .true.
        2150     4224820 :          totfrz(i)= 0._kind_phys
        2151   258292570 :          do k = pver,msg + 1,-1
        2152   257714020 :             totfrz(i)= totfrz(i)+ frz(i,k)*dz(i,k)
        2153             :          end do
        2154             :       end do
        2155    11472496 :       do k=klowest-2,khighest-1,-1
        2156    97097036 :          do i=1,il2g
        2157    96518486 :             if (doit(i) .and. k <= jb(i)-2 .and. k >= lel(i)-1) then
        2158   126640866 :                if (hu(i,k) <= hsthat(i,k) .and. hu(i,k+1) > hsthat(i,k+1) &
        2159   189961299 :                   .and. mu(i,k) >= 0.02_kind_phys) then
        2160       96390 :                   if (hu(i,k)-hsthat(i,k) < -2000._kind_phys) then
        2161       34006 :                      jt(i) = k + 1
        2162       34006 :                      doit(i) = .false.
        2163             :                   else
        2164       62384 :                      jt(i) = k
        2165       62384 :                      doit(i) = .false.
        2166             :                   end if
        2167    63224043 :                else if ( (hu(i,k) > hu(i,jb(i)) .and. totfrz(i)<=0._kind_phys) .or. mu(i,k) < 0.02_kind_phys) then
        2168     4128430 :                   jt(i) = k + 1
        2169     4128430 :                   doit(i) = .false.
        2170             :                end if
        2171             :             end if
        2172             :          end do
        2173             :       end do
        2174             :
        2175     9681835 :       if (iter == 1)  jto(:) = jt(:)
        2176             :
        2177    35291550 :       do k = pver,msg + 1,-1
        2178   288780750 :          do i = 1,il2g
        2179   253489200 :             if (k >= lel(i) .and. k <= jt(i) .and. eps0(i) > 0._kind_phys) then
        2180     9690241 :                mu(i,k) = 0._kind_phys
        2181     9690241 :                eu(i,k) = 0._kind_phys
        2182     9690241 :                du(i,k) = 0._kind_phys
        2183     9690241 :                hu(i,k) = hmn(i,k)
        2184             :             end if
        2185   288202200 :             if (k == jt(i) .and. eps0(i) > 0._kind_phys) then
        2186     4208860 :                du(i,k) = mu(i,k+1)/dz(i,k)
        2187     4208860 :                eu(i,k) = 0._kind_phys
        2188     4208860 :                mu(i,k) = 0._kind_phys
        2189             :             end if
        2190             :          end do
        2191             :       end do
        2192             :
        2193     4803370 :       do i = 1,il2g
        2194     4803370 :          done(i) = .false.
        2195             :       end do
        2196             :       kount = 0
        2197     7149672 :       do k = pver,msg + 2,-1
        2198    59663894 :          do i = 1,il2g
        2199    52552984 :             if (k == jb(i) .and. eps0(i) > 0._kind_phys) then
        2200     4208860 :                qu(i,k) = q(i,mx(i))
        2201     4208860 :                su(i,k) = (hu(i,k)-rl*qu(i,k))/cp
        2202             :             end if
        2203    59663894 :             if (( .not. done(i) .and. k > jt(i) .and. k < jb(i)) .and. eps0(i) > 0._kind_phys) then
        2204    21047964 :                su(i,k) = mu(i,k+1)/mu(i,k)*su(i,k+1) + &
        2205    21047964 :                          dz(i,k)/mu(i,k)* (eu(i,k)-du(i,k))*s(i,k)
        2206             :                qu(i,k) = mu(i,k+1)/mu(i,k)*qu(i,k+1) + dz(i,k)/mu(i,k)* (eu(i,k)*q(i,k)- &
        2207    10523982 :                                du(i,k)*qst(i,k))
        2208    10523982 :                tu = su(i,k) - grav/cp*zf(i,k)
        2209    10523982 :                call qsat_hPa(tu, (p(i,k)+p(i,k-1))/2._kind_phys, estu, qstu)
        2210    10523982 :                if (qu(i,k) >= qstu) then
        2211     4141002 :                   jlcl(i) = k
        2212     4141002 :                   kount = kount + 1
        2213     4141002 :                   done(i) = .true.
        2214             :                end if
        2215             :             end if
        2216             :          end do
        2217     7149672 :          if (kount >= il2g) goto 690
        2218             :       end do
        2219             : 690   continue
        2220    34713000 :       do k = msg + 2,pver
        2221   283977380 :          do i = 1,il2g
        2222   283398830 :             if ((k > jt(i) .and. k <= jlcl(i)) .and. eps0(i) > 0._kind_phys) then
        2223    52775017 :                su(i,k) = shat(i,k) + (hu(i,k)-hsthat(i,k))/(cp* (1._kind_phys+gamhat(i,k)))
        2224             :                qu(i,k) = qsthat(i,k) + gamhat(i,k)*(hu(i,k)-hsthat(i,k))/ &
        2225    52775017 :                         (rl* (1._kind_phys+gamhat(i,k)))
        2226             :             end if
        2227             :          end do
        2228             :       end do
        2229             :
        2230             : ! compute condensation in updraft
        2231     4803370 :       tmplel(:il2g) = jb(:il2g)
        2232             :
        2233    34713000 :       do k = pver,msg + 2,-1
        2234   283977380 :          do i = 1,il2g
        2235   283398830 :              if (k >= jt(i) .and. k < tmplel(i) .and. eps0(i) > 0._kind_phys) then
        2236             :
        2237   126733714 :                cu(i,k) = ((mu(i,k)*su(i,k)-mu(i,k+1)*su(i,k+1))/ &
        2238   126733714 :                          dz(i,k)- (eu(i,k)-du(i,k))*s(i,k))/(rl/cp)
        2239    63366857 :                if (k == jt(i)) cu(i,k) = 0._kind_phys
        2240    63366857 :                cu(i,k) = max(0._kind_phys,cu(i,k))
        2241             :             end if
        2242             :          end do
        2243             :       end do
        2244             :
        2245             :
        2246             : ! compute condensed liquid, rain production rate
        2247             : ! accumulate total precipitation (condensation - detrainment of liquid)
        2248             : ! Note ql1 = ql(k) + rprd(k)*dz(k)/mu(k)
        2249             : ! The differencing is somewhat strange (e.g. du(i,k)*ql(i,k+1)) but is
        2250             : ! consistently applied.
        2251             : !    mu, ql are interface quantities
        2252             : !    cu, du, eu, rprd are midpoint quantites
        2253             :
        2254    35291550 :          do k = pver,msg + 2,-1
        2255   283977380 :             do i = 1,il2g
        2256   249264380 :                rprd(i,k) = 0._kind_phys
        2257   283398830 :                if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0._kind_phys .and. mu(i,k) >= 0.0_kind_phys) then
        2258    63366857 :                   if (mu(i,k) > 0._kind_phys) then
        2259    59157997 :                      ql1 = 1._kind_phys/mu(i,k)* (mu(i,k+1)*ql(i,k+1)- &
        2260    59157997 :                            dz(i,k)*du(i,k)*ql(i,k+1)+dz(i,k)*cu(i,k))
        2261    59157997 :                      ql(i,k) = ql1/ (1._kind_phys+dz(i,k)*c0mask(i))
        2262             :                   else
        2263     4208860 :                      ql(i,k) = 0._kind_phys
        2264             :                   end if
        2265    63366857 :                   totpcp(i) = totpcp(i) + dz(i,k)*(cu(i,k)-du(i,k)*ql(i,k+1))
        2266    63366857 :                   rprd(i,k) = c0mask(i)*mu(i,k)*ql(i,k)
        2267    63366857 :                   qcde(i,k) = ql(i,k)
        2268             :
        2269             :
        2270             :                end if
        2271             :             end do
        2272             :          end do
        2273             : !
        2274             :    end do   !iter
        2275             : !
        2276             : ! specify downdraft properties (no downdrafts if jd.ge.jb).
        2277             : ! scale down downward mass flux profile so that net flux
        2278             : ! (up-down) at cloud base in not negative.
        2279             : !
        2280     4803370 :    do i = 1,il2g
        2281             : !
        2282             : ! in normal downdraft strength run alfa=0.2.  In test4 alfa=0.1
        2283             : !
        2284     4224820 :       alfa(i) = 0.1_kind_phys
        2285     4224820 :       jt(i) = min(jt(i),jb(i)-1)
        2286     4224820 :       jd(i) = max(j0(i),jt(i)+1)
        2287     4224820 :       jd(i) = min(jd(i),jb(i))
        2288     4224820 :       hd(i,jd(i)) = hmn(i,jd(i)-1)
        2289     4803370 :       if (jd(i) < jb(i) .and. eps0(i) > 0._kind_phys) then
        2290     4152989 :          epsm(i) = eps0(i)
        2291     4152989 :          md(i,jd(i)) = -alfa(i)*epsm(i)/eps0(i)
        2292             :       end if
        2293             :    end do
        2294    35291550 :    do k = msg + 1,pver
        2295   288780750 :       do i = 1,il2g
        2296   288202200 :          if ((k > jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._kind_phys) then
        2297    48167037 :             zdef(i) = zf(i,jd(i)) - zf(i,k)
        2298    48167037 :             md(i,k) = -alfa(i)/ (2._kind_phys*eps0(i))*(exp(2._kind_phys*epsm(i)*zdef(i))-1._kind_phys)/zdef(i)
        2299             :          end if
        2300             :       end do
        2301             :    end do
        2302             :
        2303    35291550 :    do k = msg + 1,pver
        2304   288780750 :       do i = 1,il2g
        2305   288202200 :          if ((k >= jt(i) .and. k <= jb(i)) .and. eps0(i) > 0._kind_phys .and. jd(i) < jb(i)) then
        2306    67463975 :             ratmjb(i) = min(abs(mu(i,jb(i))/md(i,jb(i))),1._kind_phys)
        2307    67463975 :             md(i,k) = md(i,k)*ratmjb(i)
        2308             :          end if
        2309             :       end do
        2310             :    end do
        2311             :
        2312             :    small = 1.e-20_kind_phys
        2313    35291550 :    do k = msg + 1,pver
        2314   288780750 :       do i = 1,il2g
        2315   288202200 :          if ((k >= jt(i) .and. k <= pver) .and. eps0(i) > 0._kind_phys) then
        2316    86502733 :             ed(i,k-1) = (md(i,k-1)-md(i,k))/dz(i,k-1)
        2317    86502733 :             mdt = min(md(i,k),-small)
        2318    86502733 :             hd(i,k) = (md(i,k-1)*hd(i,k-1) - dz(i,k-1)*ed(i,k-1)*hmn(i,k-1))/mdt
        2319             :          end if
        2320             :       end do
        2321             :    end do
        2322             : !
        2323             : ! calculate updraft and downdraft properties.
        2324             : !
        2325    34713000 :    do k = msg + 2,pver
        2326   283977380 :       do i = 1,il2g
        2327   283398830 :          if ((k >= jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._kind_phys .and. jd(i) < jb(i)) then
        2328    52320026 :             qds(i,k) = qsthat(i,k) + gamhat(i,k)*(hd(i,k)-hsthat(i,k))/ &
        2329   104640052 :                (rl*(1._kind_phys + gamhat(i,k)))
        2330             :          end if
        2331             :       end do
        2332             :    end do
        2333             :
        2334     4803370 :    do i = 1,il2g
        2335     4224820 :       qd(i,jd(i)) = qds(i,jd(i))
        2336     4803370 :       sd(i,jd(i)) = (hd(i,jd(i)) - rl*qd(i,jd(i)))/cp
        2337             :    end do
        2338             : !
        2339    34713000 :    do k = msg + 2,pver
        2340   283977380 :       do i = 1,il2g
        2341   283398830 :          if (k >= jd(i) .and. k < jb(i) .and. eps0(i) > 0._kind_phys) then
        2342    48167037 :             qd(i,k+1) = qds(i,k+1)
        2343    48167037 :             evp(i,k) = -ed(i,k)*q(i,k) + (md(i,k)*qd(i,k)-md(i,k+1)*qd(i,k+1))/dz(i,k)
        2344    48167037 :             evp(i,k) = max(evp(i,k),0._kind_phys)
        2345    48167037 :             mdt = min(md(i,k+1),-small)
        2346    48167037 :             sd(i,k+1) = ((rl/cp*evp(i,k)-ed(i,k)*s(i,k))*dz(i,k) + md(i,k)*sd(i,k))/mdt
        2347    48167037 :             totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k)
        2348             :          end if
        2349             :       end do
        2350             :    end do
        2351     4803370 :    do i = 1,il2g
        2352     4803370 :       totevp(i) = totevp(i) + md(i,jd(i))*qd(i,jd(i)) - md(i,jb(i))*qd(i,jb(i))
        2353             :    end do
        2354             : !!$   if (.true.) then
        2355             :    if (.false.) then
        2356             :       do i = 1,il2g
        2357             :          k = jb(i)
        2358             :          if (eps0(i) > 0._kind_phys) then
        2359             :             evp(i,k) = -ed(i,k)*q(i,k) + (md(i,k)*qd(i,k))/dz(i,k)
        2360             :             evp(i,k) = max(evp(i,k),0._kind_phys)
        2361             :             totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k)
        2362             :          end if
        2363             :       end do
        2364             :    endif
        2365             :
        2366     4803370 :    do i = 1,il2g
        2367     4224820 :       totpcp(i) = max(totpcp(i),0._kind_phys)
        2368     4803370 :       totevp(i) = max(totevp(i),0._kind_phys)
        2369             :    end do
        2370             : !
        2371    34713000 :    do k = msg + 2,pver
        2372   283977380 :       do i = 1,il2g
        2373   249264380 :          if (totevp(i) > 0._kind_phys .and. totpcp(i) > 0._kind_phys) then
        2374   244844690 :             md(i,k)  = md (i,k)*min(1._kind_phys, totpcp(i)/(totevp(i)+totpcp(i)))
        2375   244844690 :             ed(i,k)  = ed (i,k)*min(1._kind_phys, totpcp(i)/(totevp(i)+totpcp(i)))
        2376   244844690 :             evp(i,k) = evp(i,k)*min(1._kind_phys, totpcp(i)/(totevp(i)+totpcp(i)))
        2377             :          else
        2378     4419690 :             md(i,k) = 0._kind_phys
        2379     4419690 :             ed(i,k) = 0._kind_phys
        2380     4419690 :             evp(i,k) = 0._kind_phys
        2381             :          end if
        2382             : ! cmeg is the cloud water condensed - rain water evaporated
        2383             : ! rprd is the cloud water converted to rain - (rain evaporated)
        2384   249264380 :          cmeg(i,k) = cu(i,k) - evp(i,k)
        2385   283398830 :          rprd(i,k) = rprd(i,k)-evp(i,k)
        2386             :       end do
        2387             :    end do
        2388             :
        2389             : !
        2390    35291550 :    do k = msg + 1,pver
        2391   288780750 :       do i = 1,il2g
        2392   288202200 :          mc(i,k) = mu(i,k) + md(i,k)
        2393             :       end do
        2394             :    end do
        2395             : !
        2396      578550 :    return
        2397             : end subroutine cldprp
        2398             :
        2399      578550 : subroutine closure(ncol   ,pver, &
        2400      578550 :                    q       ,t       ,p       ,z       ,s       , &
        2401      578550 :                    tp      ,qs      ,qu      ,su      ,mc      , &
        2402      578550 :                    du      ,mu      ,md      ,qd      ,sd      , &
        2403      578550 :                    qhat    ,shat    ,dp      ,qstp    ,zf      , &
        2404      578550 :                    ql      ,dsubcld ,mb      ,cape    ,tl      , &
        2405      578550 :                    lcl     ,lel     ,jt      ,mx      ,il1g    , &
        2406             :                    il2g    ,rd      ,grav    ,cp      ,rl      , &
        2407             :                    msg     ,capelmt )
        2408             : !-----------------------------------------------------------------------
        2409             : !CACNOTE - fill in documentation
        2410             : !
        2411             : ! Purpose:
        2412             : ! <Say what the routine does>
        2413             : !
        2414             : ! Method:
        2415             : ! <Describe the algorithm(s) used in the routine.>
        2416             : ! <Also include any applicable external references.>
        2417             : !
        2418             : ! Author: G. Zhang and collaborators. CCM contact:P. Rasch
        2419             : ! This is contributed code not fully standardized by the CCM core group.
        2420             : !
        2421             : ! this code is very much rougher than virtually anything else in the CCM
        2422             : ! We expect to release cleaner code in a future release
        2423             : !
        2424             : ! the documentation has been enhanced to the degree that we are able
        2425             : !
        2426             : !-----------------------------------------------------------------------
        2427             :
        2428             : !
        2429             : !-----------------------------Arguments---------------------------------
        2430             : !
        2431             :    integer, intent(in) :: ncol
        2432             :    integer, intent(in) :: pver
        2433             :
        2434             :    real(kind_phys), intent(inout) :: q(ncol,pver)        ! spec humidity
        2435             :    real(kind_phys), intent(inout) :: t(ncol,pver)        ! temperature
        2436             :    real(kind_phys), intent(inout) :: p(ncol,pver)        ! pressure (mb)
        2437             :    real(kind_phys), intent(inout) :: mb(ncol)            ! cloud base mass flux
        2438             :    real(kind_phys), intent(in) :: z(ncol,pver)        ! height (m)
        2439             :    real(kind_phys), intent(in) :: s(ncol,pver)        ! normalized dry static energy
        2440             :    real(kind_phys), intent(in) :: tp(ncol,pver)       ! parcel temp
        2441             :    real(kind_phys), intent(in) :: qs(ncol,pver)       ! sat spec humidity
        2442             :    real(kind_phys), intent(in) :: qu(ncol,pver)       ! updraft spec. humidity
        2443             :    real(kind_phys), intent(in) :: su(ncol,pver)       ! normalized dry stat energy of updraft
        2444             :    real(kind_phys), intent(in) :: mc(ncol,pver)       ! net convective mass flux
        2445             :    real(kind_phys), intent(in) :: du(ncol,pver)       ! detrainment from updraft
        2446             :    real(kind_phys), intent(in) :: mu(ncol,pver)       ! mass flux of updraft
        2447             :    real(kind_phys), intent(in) :: md(ncol,pver)       ! mass flux of downdraft
        2448             :    real(kind_phys), intent(in) :: qd(ncol,pver)       ! spec. humidity of downdraft
        2449             :    real(kind_phys), intent(in) :: sd(ncol,pver)       ! dry static energy of downdraft
        2450             :    real(kind_phys), intent(in) :: qhat(ncol,pver)     ! environment spec humidity at interfaces
        2451             :    real(kind_phys), intent(in) :: shat(ncol,pver)     ! env. normalized dry static energy at intrfcs
        2452             :    real(kind_phys), intent(in) :: dp(ncol,pver)       ! pressure thickness of layers
        2453             :    real(kind_phys), intent(in) :: qstp(ncol,pver)     ! spec humidity of parcel
        2454             :    real(kind_phys), intent(in) :: zf(ncol,pver+1)     ! height of interface levels
        2455             :    real(kind_phys), intent(in) :: ql(ncol,pver)       ! liquid water mixing ratio
        2456             :
        2457             :    real(kind_phys), intent(in) :: cape(ncol)          ! available pot. energy of column
        2458             :    real(kind_phys), intent(in) :: tl(ncol)
        2459             :    real(kind_phys), intent(in) :: dsubcld(ncol)       ! thickness of subcloud layer
        2460             :
        2461             :    integer, intent(in) :: lcl(ncol)        ! index of lcl
        2462             :    integer, intent(in) :: lel(ncol)        ! index of launch leve
        2463             :    integer, intent(in) :: jt(ncol)         ! top of updraft
        2464             :    integer, intent(in) :: mx(ncol)         ! base of updraft
        2465             : !
        2466             : !--------------------------Local variables------------------------------
        2467             : !
        2468     1157100 :    real(kind_phys) dtpdt(ncol,pver)
        2469     1157100 :    real(kind_phys) dqsdtp(ncol,pver)
        2470     1157100 :    real(kind_phys) dtmdt(ncol,pver)
        2471     1157100 :    real(kind_phys) dqmdt(ncol,pver)
        2472     1157100 :    real(kind_phys) dboydt(ncol,pver)
        2473     1157100 :    real(kind_phys) thetavp(ncol,pver)
        2474     1157100 :    real(kind_phys) thetavm(ncol,pver)
        2475             :
        2476     1157100 :    real(kind_phys) dtbdt(ncol),dqbdt(ncol),dtldt(ncol)
        2477             :    real(kind_phys) beta
        2478             :    real(kind_phys) capelmt
        2479             :    real(kind_phys) cp
        2480     1157100 :    real(kind_phys) dadt(ncol)
        2481             :    real(kind_phys) debdt
        2482             :    real(kind_phys) dltaa
        2483             :    real(kind_phys) eb
        2484             :    real(kind_phys) grav
        2485             :
        2486             :    integer i
        2487             :    integer il1g
        2488             :    integer il2g
        2489             :    integer k, kmin, kmax
        2490             :    integer msg
        2491             :
        2492             :    real(kind_phys) rd
        2493             :    real(kind_phys) rl
        2494             : ! change of subcloud layer properties due to convection is
        2495             : ! related to cumulus updrafts and downdrafts.
        2496             : ! mc(z)=f(z)*mb, mub=betau*mb, mdb=betad*mb are used
        2497             : ! to define betau, betad and f(z).
        2498             : ! note that this implies all time derivatives are in effect
        2499             : ! time derivatives per unit cloud-base mass flux, i.e. they
        2500             : ! have units of 1/mb instead of 1/sec.
        2501             : !
        2502     4803370 :    do i = il1g,il2g
        2503     4224820 :       mb(i) = 0._kind_phys
        2504     4224820 :       eb = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i)))
        2505             :       dtbdt(i) = (1._kind_phys/dsubcld(i))* (mu(i,mx(i))*(shat(i,mx(i))-su(i,mx(i)))+ &
        2506     4224820 :                   md(i,mx(i))* (shat(i,mx(i))-sd(i,mx(i))))
        2507             :       dqbdt(i) = (1._kind_phys/dsubcld(i))* (mu(i,mx(i))*(qhat(i,mx(i))-qu(i,mx(i)))+ &
        2508     4224820 :                  md(i,mx(i))* (qhat(i,mx(i))-qd(i,mx(i))))
        2509     4224820 :       debdt = eps1*p(i,mx(i))/ (eps1+q(i,mx(i)))**2*dqbdt(i)
        2510             :       dtldt(i) = -2840._kind_phys* (3.5_kind_phys/t(i,mx(i))*dtbdt(i)-debdt/eb)/ &
        2511     4803370 :                  (3.5_kind_phys*log(t(i,mx(i)))-log(eb)-4.805_kind_phys)**2
        2512             :    end do
        2513             : !
        2514             : !   dtmdt and dqmdt are cumulus heating and drying.
        2515             : !
        2516    35291550 :    do k = msg + 1,pver
        2517   288780750 :       do i = il1g,il2g
        2518   253489200 :          dtmdt(i,k) = 0._kind_phys
        2519   288202200 :          dqmdt(i,k) = 0._kind_phys
        2520             :       end do
        2521             :    end do
        2522             : !
        2523    34713000 :    do k = msg + 1,pver - 1
        2524   283977380 :       do i = il1g,il2g
        2525   283398830 :          if (k == jt(i)) then
        2526     8449640 :             dtmdt(i,k) = (1._kind_phys/dp(i,k))*(mu(i,k+1)* (su(i,k+1)-shat(i,k+1)- &
        2527     8449640 :                           rl/cp*ql(i,k+1))+md(i,k+1)* (sd(i,k+1)-shat(i,k+1)))
        2528             :             dqmdt(i,k) = (1._kind_phys/dp(i,k))*(mu(i,k+1)* (qu(i,k+1)- &
        2529     4224820 :                          qhat(i,k+1)+ql(i,k+1))+md(i,k+1)*(qd(i,k+1)-qhat(i,k+1)))
        2530             :          end if
        2531             :       end do
        2532             :    end do
        2533             : !
        2534      578550 :    beta = 0._kind_phys
        2535    34713000 :    do k = msg + 1,pver - 1
        2536   283977380 :       do i = il1g,il2g
        2537   283398830 :          if (k > jt(i) .and. k < mx(i)) then
        2538   118315994 :             dtmdt(i,k) = (mc(i,k)* (shat(i,k)-s(i,k))+mc(i,k+1)* (s(i,k)-shat(i,k+1)))/ &
        2539   118315994 :                          dp(i,k) - rl/cp*du(i,k)*(beta*ql(i,k)+ (1-beta)*ql(i,k+1))
        2540             :
        2541             :             dqmdt(i,k) = (mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)+cp/rl* (su(i,k+1)-s(i,k)))- &
        2542             :                           mu(i,k)* (qu(i,k)-qhat(i,k)+cp/rl*(su(i,k)-s(i,k)))+md(i,k+1)* &
        2543             :                          (qd(i,k+1)-qhat(i,k+1)+cp/rl*(sd(i,k+1)-s(i,k)))-md(i,k)* &
        2544             :                          (qd(i,k)-qhat(i,k)+cp/rl*(sd(i,k)-s(i,k))))/dp(i,k) + &
        2545    59157997 :                           du(i,k)* (beta*ql(i,k)+(1-beta)*ql(i,k+1))
        2546             :          end if
        2547             :       end do
        2548             :    end do
        2549             : !
        2550    35291550 :    do k = msg + 1,pver
        2551   288780750 :       do i = il1g,il2g
        2552   288202200 :          if (k >= lel(i) .and. k <= lcl(i)) then
        2553    63991323 :             thetavp(i,k) = tp(i,k)* (1000._kind_phys/p(i,k))** (rd/cp)*(1._kind_phys+1.608_kind_phys*qstp(i,k)-q(i,mx(i)))
        2554    63991323 :             thetavm(i,k) = t(i,k)* (1000._kind_phys/p(i,k))** (rd/cp)*(1._kind_phys+0.608_kind_phys*q(i,k))
        2555    63991323 :             dqsdtp(i,k) = qstp(i,k)* (1._kind_phys+qstp(i,k)/eps1)*eps1*rl/(rd*tp(i,k)**2)
        2556             : !
        2557             : ! dtpdt is the parcel temperature change due to change of
        2558             : ! subcloud layer properties during convection.
        2559             : !
        2560             :             dtpdt(i,k) = tp(i,k)/ (1._kind_phys+rl/cp* (dqsdtp(i,k)-qstp(i,k)/tp(i,k)))* &
        2561             :                         (dtbdt(i)/t(i,mx(i))+rl/cp* (dqbdt(i)/tl(i)-q(i,mx(i))/ &
        2562    63991323 :                          tl(i)**2*dtldt(i)))
        2563             : !
        2564             : ! dboydt is the integrand of cape change.
        2565             : !
        2566             :             dboydt(i,k) = ((dtpdt(i,k)/tp(i,k)+1._kind_phys/(1._kind_phys+1.608_kind_phys*qstp(i,k)-q(i,mx(i)))* &
        2567             :                           (1.608_kind_phys * dqsdtp(i,k) * dtpdt(i,k) -dqbdt(i))) - (dtmdt(i,k)/t(i,k)+0.608_kind_phys/ &
        2568    63991323 :                           (1._kind_phys+0.608_kind_phys*q(i,k))*dqmdt(i,k)))*grav*thetavp(i,k)/thetavm(i,k)
        2569             :          end if
        2570             :       end do
        2571             :    end do
        2572             : !
        2573    35291550 :    do k = msg + 1,pver
        2574   288780750 :       do i = il1g,il2g
        2575   288202200 :          if (k > lcl(i) .and. k < mx(i)) then
        2576     5455621 :             thetavp(i,k) = tp(i,k)* (1000._kind_phys/p(i,k))** (rd/cp)*(1._kind_phys+0.608_kind_phys*q(i,mx(i)))
        2577     5455621 :             thetavm(i,k) = t(i,k)* (1000._kind_phys/p(i,k))** (rd/cp)*(1._kind_phys+0.608_kind_phys*q(i,k))
        2578             : !
        2579             : ! dboydt is the integrand of cape change.
        2580             : !
        2581             :             dboydt(i,k) = (dtbdt(i)/t(i,mx(i))+0.608_kind_phys/ (1._kind_phys+0.608_kind_phys*q(i,mx(i)))*dqbdt(i)- &
        2582             :                           dtmdt(i,k)/t(i,k)-0.608_kind_phys/ (1._kind_phys+0.608_kind_phys*q(i,k))*dqmdt(i,k))* &
        2583     5455621 :                           grav*thetavp(i,k)/thetavm(i,k)
        2584             :          end if
        2585             :       end do
        2586             :    end do
        2587             :
        2588             : !
        2589             : ! buoyant energy change is set to 2/3*excess cape per 3 hours
        2590             : !
        2591     4803370 :    dadt(il1g:il2g)  = 0._kind_phys
        2592     4803370 :    kmin = minval(lel(il1g:il2g))
        2593     4803370 :    kmax = maxval(mx(il1g:il2g)) - 1
        2594    11472496 :    do k = kmin, kmax
        2595    97097036 :       do i = il1g,il2g
        2596    96518486 :          if ( k >= lel(i) .and. k <= mx(i) - 1) then
        2597    69224676 :             dadt(i) = dadt(i) + dboydt(i,k)* (zf(i,k)-zf(i,k+1))
        2598             :          endif
        2599             :       end do
        2600             :    end do
        2601     4803370 :    do i = il1g,il2g
        2602     4224820 :       dltaa = -1._kind_phys* (cape(i)-capelmt)
        2603     4803370 :       if (dadt(i) /= 0._kind_phys) mb(i) = max(dltaa/tau/dadt(i),0._kind_phys)
        2604             :    end do
        2605             : !
        2606      578550 :    return
        2607             : end subroutine closure
        2608             :
        2609      578550 : subroutine q1q2_pjr(ncol   ,pver    ,latice  ,&
        2610      578550 :                     dqdt    ,dsdt    ,q       ,qs      ,qu      , &
        2611      578550 :                     su      ,du      ,qhat    ,shat    ,dp      , &
        2612      578550 :                     mu      ,md      ,sd      ,qd      ,ql      , &
        2613      578550 :                     dsubcld ,jt      ,mx      ,il1g    ,il2g    , &
        2614             :                     cp      ,rl      ,msg     ,          &
        2615      578550 :                     dl      ,evp     ,cu)
        2616             :
        2617             :    implicit none
        2618             :
        2619             : !-----------------------------------------------------------------------
        2620             : ! CACNOTE -fill in documentation
        2621             : !
        2622             : ! Purpose:
        2623             : ! <Say what the routine does>
        2624             : !
        2625             : ! Method:
        2626             : ! <Describe the algorithm(s) used in the routine.>
        2627             : ! <Also include any applicable external references.>
        2628             : !
        2629             : ! Author: phil rasch dec 19 1995
        2630             : !
        2631             : !-----------------------------------------------------------------------
        2632             :
        2633             :
        2634             :    real(kind_phys), intent(in) :: cp
        2635             :
        2636             :    integer, intent(in) :: ncol
        2637             :    integer, intent(in) :: pver
        2638             :    real(kind_phys), intent(in) :: latice
        2639             :    integer, intent(in) :: il1g
        2640             :    integer, intent(in) :: il2g
        2641             :    integer, intent(in) :: msg
        2642             :
        2643             :    real(kind_phys), intent(in) :: q(ncol,pver)
        2644             :    real(kind_phys), intent(in) :: qs(ncol,pver)
        2645             :    real(kind_phys), intent(in) :: qu(ncol,pver)
        2646             :    real(kind_phys), intent(in) :: su(ncol,pver)
        2647             :    real(kind_phys), intent(in) :: du(ncol,pver)
        2648             :    real(kind_phys), intent(in) :: qhat(ncol,pver)
        2649             :    real(kind_phys), intent(in) :: shat(ncol,pver)
        2650             :    real(kind_phys), intent(in) :: dp(ncol,pver)
        2651             :    real(kind_phys), intent(in) :: mu(ncol,pver)
        2652             :    real(kind_phys), intent(in) :: md(ncol,pver)
        2653             :    real(kind_phys), intent(in) :: sd(ncol,pver)
        2654             :    real(kind_phys), intent(in) :: qd(ncol,pver)
        2655             :    real(kind_phys), intent(in) :: ql(ncol,pver)
        2656             :    real(kind_phys), intent(in) :: evp(ncol,pver)
        2657             :    real(kind_phys), intent(in) :: cu(ncol,pver)
        2658             :    real(kind_phys), intent(in) :: dsubcld(ncol)
        2659             :
        2660             :    real(kind_phys),intent(out) :: dqdt(ncol,pver),dsdt(ncol,pver)
        2661             :    real(kind_phys),intent(out) :: dl(ncol,pver)
        2662             :
        2663             :    integer kbm
        2664             :    integer ktm
        2665             :    integer jt(ncol)
        2666             :    integer mx(ncol)
        2667             : !
        2668             : ! work fields:
        2669             : !
        2670             :    integer i
        2671             :    integer k
        2672             :
        2673             :    real(kind_phys) emc
        2674             :    real(kind_phys) rl
        2675             : !-------------------------------------------------------------------
        2676    35291550 :    do k = msg + 1,pver
        2677   288780750 :       do i = il1g,il2g
        2678   253489200 :          dsdt(i,k) = 0._kind_phys
        2679   253489200 :          dqdt(i,k) = 0._kind_phys
        2680   288202200 :          dl(i,k) = 0._kind_phys
        2681             :       end do
        2682             :    end do
        2683             :
        2684             : !
        2685             : ! find the highest level top and bottom levels of convection
        2686             : !
        2687      578550 :    ktm = pver
        2688      578550 :    kbm = pver
        2689     4803370 :    do i = il1g, il2g
        2690     4224820 :       ktm = min(ktm,jt(i))
        2691     4803370 :       kbm = min(kbm,mx(i))
        2692             :    end do
        2693             :
        2694    12905130 :    do k = ktm,pver-1
        2695   109984991 :       do i = il1g,il2g
        2696   194159722 :          emc = -cu (i,k)               &         ! condensation in updraft
        2697   194159722 :                +evp(i,k)                         ! evaporating rain in downdraft
        2698             :
        2699             :          dsdt(i,k) = -rl/cp*emc &
        2700    97079861 :                      + (+mu(i,k+1)* (su(i,k+1)-shat(i,k+1)) &
        2701             :                         -mu(i,k)*   (su(i,k)-shat(i,k)) &
        2702             :                         +md(i,k+1)* (sd(i,k+1)-shat(i,k+1)) &
        2703             :                         -md(i,k)*   (sd(i,k)-shat(i,k)) &
        2704    97079861 :                        )/dp(i,k)
        2705             :
        2706             :          dqdt(i,k) = emc + &
        2707             :                     (+mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)) &
        2708             :                      -mu(i,k)*   (qu(i,k)-qhat(i,k)) &
        2709             :                      +md(i,k+1)* (qd(i,k+1)-qhat(i,k+1)) &
        2710             :                      -md(i,k)*   (qd(i,k)-qhat(i,k)) &
        2711    97079861 :                     )/dp(i,k)
        2712             :
        2713   109406441 :          dl(i,k) = du(i,k)*ql(i,k+1)
        2714             :
        2715             :       end do
        2716             :    end do
        2717             :
        2718             : !
        2719     4123638 :    do k = kbm,pver
        2720    30999945 :       do i = il1g,il2g
        2721    30421395 :          if (k == mx(i)) then
        2722     4224820 :             dsdt(i,k) = (1._kind_phys/dsubcld(i))* &
        2723             :                         (-mu(i,k)* (su(i,k)-shat(i,k)) &
        2724             :                          -md(i,k)* (sd(i,k)-shat(i,k)) &
        2725     4224820 :                         )
        2726             :             dqdt(i,k) = (1._kind_phys/dsubcld(i))* &
        2727             :                         (-mu(i,k)*(qu(i,k)-qhat(i,k)) &
        2728             :                          -md(i,k)*(qd(i,k)-qhat(i,k)) &
        2729     4224820 :                         )
        2730    22651487 :          else if (k > mx(i)) then
        2731    18993359 :             dsdt(i,k) = dsdt(i,k-1)
        2732    18993359 :             dqdt(i,k) = dqdt(i,k-1)
        2733             :          end if
        2734             :       end do
        2735             :    end do
        2736             : !
        2737      578550 :    return
        2738             : end subroutine q1q2_pjr
        2739             :
        2740             :
        2741             : ! Wrapper for qsat_water that does translation between Pa and hPa
        2742             : ! qsat_water uses Pa internally, so get it right, need to pass in Pa.
        2743             : ! Afterward, set es back to hPa.
        2744 29020232765 : subroutine qsat_hPa(t, p, es, qm)
        2745             : !CACNOTE - Need to figure out how to handle this
        2746             :   use wv_saturation, only: qsat_water
        2747             :
        2748             :   ! Inputs
        2749             :   real(kind_phys), intent(in) :: t    ! Temperature (K)
        2750             :   real(kind_phys), intent(in) :: p    ! Pressure (hPa)
        2751             :   ! Outputs
        2752             :   real(kind_phys), intent(out) :: es  ! Saturation vapor pressure (hPa)
        2753             :   real(kind_phys), intent(out) :: qm  ! Saturation mass mixing ratio
        2754             :                                ! (vapor mass over dry mass, kg/kg)
        2755             :
        2756 29020232765 :   call qsat_water(t, p*100._kind_phys, es, qm)
        2757             :
        2758 29020232765 :   es = es*0.01_kind_phys
        2759             :
        2760 29020232765 : end subroutine qsat_hPa
        2761             :
        2762             : end module zm_convr

Generated by: LCOV version 1.14