LCOV - code coverage report
Current view: | top level - physics/cam - hk_conv.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 13 | 323 | 4.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 1 | 3 | 33.3 %

              Line data    Source code
           1             :
           2             : module hk_conv
           3             : !
           4             : ! Moist convection. Primarily data used by both Zhang-McFarlane convection
           5             : ! and Hack shallow convective schemes.
           6             : !
           7             : ! $Id$
           8             : !
           9             :    use shr_kind_mod,     only: r8 => shr_kind_r8
          10             :    use cam_logfile,      only: iulog
          11             :    use spmd_utils,       only: masterproc
          12             :    use cam_abortutils,   only: endrun
          13             :    implicit none
          14             :
          15             :    private
          16             :    save
          17             : !
          18             : ! Public interfaces
          19             : !
          20             :    public mfinti   !  Initialization of data for moist convection
          21             :    public cmfmca   !  Hack shallow convection
          22             :    public hkconv_readnl ! read hkconv_nl namelist
          23             :
          24             : !
          25             : ! Private data used for Hack shallow convection
          26             : !
          27             :    real(r8), parameter :: unset_r8 = huge(1.0_r8)
          28             :
          29             :   ! Namelist variables
          30             :    real(r8) :: hkconv_c0 = unset_r8
          31             :    real(r8) :: hkconv_cmftau = unset_r8
          32             :
          33             :    real(r8) :: hlat        ! latent heat of vaporization
          34             :    real(r8) :: c0          ! rain water autoconversion coefficient set from namelist input hkconv_c0
          35             :    real(r8) :: betamn      ! minimum overshoot parameter
          36             :    real(r8) :: rhlat       ! reciprocal of hlat
          37             :    real(r8) :: rcp         ! reciprocal of cp
          38             :    real(r8) :: cmftau      ! characteristic adjustment time scale set from namelist input hkconv_cmftau
          39             :    real(r8) :: rhoh2o      ! density of liquid water (STP)
          40             :    real(r8) :: dzmin       ! minimum convective depth for precipitation
          41             :    real(r8) :: tiny        ! arbitrary small num used in transport estimates
          42             :    real(r8) :: eps         ! convergence criteria (machine dependent)
          43             :    real(r8) :: tpmax       ! maximum acceptable t perturbation (degrees C)
          44             :    real(r8) :: shpmax      ! maximum acceptable q perturbation (g/g)
          45             :
          46             :    integer :: iloc         ! longitude location for diagnostics
          47             :    integer :: jloc         ! latitude  location for diagnostics
          48             :    integer :: nsloc        ! nstep for which to produce diagnostics
          49             : !
          50             :    logical :: rlxclm       ! logical to relax column versus cloud triplet
          51             :
          52             :    real(r8) cp          ! specific heat of dry air
          53             :    real(r8) grav        ! gravitational constant
          54             :    real(r8) rgrav       ! reciprocal of grav
          55             :    real(r8) rgas        ! gas constant for dry air
          56             :    integer limcnv          ! top interface level limit for convection
          57             :
          58             :
          59             :
          60             :
          61             : contains
          62        1536 : subroutine hkconv_readnl(nlfile)
          63             :
          64             :    use namelist_utils,  only: find_group_name
          65             :    use units,           only: getunit, freeunit
          66             :    use mpishorthand
          67             :
          68             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          69             :
          70             :    ! Local variables
          71             :    integer :: unitn, ierr
          72             :    character(len=*), parameter :: subname = 'hkconv_readnl'
          73             :
          74             :    namelist /hkconv_nl/ hkconv_cmftau, hkconv_c0
          75             :    !-----------------------------------------------------------------------------
          76             :
          77        1536 :    if (masterproc) then
          78           2 :       unitn = getunit()
          79           2 :       open( unitn, file=trim(nlfile), status='old' )
          80           2 :       call find_group_name(unitn, 'hkconv_nl', status=ierr)
          81           2 :       if (ierr == 0) then
          82           0 :          read(unitn, hkconv_nl, iostat=ierr)
          83           0 :          if (ierr /= 0) then
          84           0 :             call endrun(subname // ':: ERROR reading namelist')
          85             :          end if
          86             :       end if
          87           2 :       close(unitn)
          88           2 :       call freeunit(unitn)
          89             :
          90             :       ! set local variables
          91           2 :       cmftau = hkconv_cmftau
          92           2 :       c0     = hkconv_c0
          93             :
          94             :    end if
          95             :
          96             : #ifdef SPMD
          97             :    ! Broadcast namelist variables
          98        1536 :    call mpibcast(cmftau,            1, mpir8,  0, mpicom)
          99        1536 :    call mpibcast(c0,                1, mpir8,  0, mpicom)
         100             : #endif
         101             :
         102        1536 : end subroutine hkconv_readnl
         103             :
         104             : !================================================================================================
         105             :
         106           0 : subroutine mfinti (rair    ,cpair   ,gravit  ,latvap  ,rhowtr,limcnv_in )
         107             : !-----------------------------------------------------------------------
         108             : !
         109             : ! Purpose:
         110             : ! Initialize moist convective mass flux procedure common block, cmfmca
         111             : !
         112             : ! Method:
         113             : ! <Describe the algorithm(s) used in the routine.>
         114             : ! <Also include any applicable external references.>
         115             : !
         116             : ! Author: J. Hack
         117             : !
         118             : !-----------------------------------------------------------------------
         119             :    use spmd_utils, only: masterproc
         120             : !------------------------------Arguments--------------------------------
         121             : !
         122             : ! Input arguments
         123             : !
         124             :    real(r8), intent(in) :: rair              ! gas constant for dry air
         125             :    real(r8), intent(in) :: cpair             ! specific heat of dry air
         126             :    real(r8), intent(in) :: gravit            ! acceleration due to gravity
         127             :    real(r8), intent(in) :: latvap            ! latent heat of vaporization
         128             :    real(r8), intent(in) :: rhowtr            ! density of liquid water (STP)
         129             :    integer,  intent(in) :: limcnv_in         ! top interface level limit for convection
         130             :
         131             :    ! local variables
         132             :    character(len=32)    :: hgrid             ! horizontal grid specifier
         133             : !
         134             : !-----------------------------------------------------------------------
         135             : !
         136             : ! Initialize physical constants for moist convective mass flux procedure
         137             : !
         138           0 :    cp     = cpair         ! specific heat of dry air
         139           0 :    hlat   = latvap        ! latent heat of vaporization
         140           0 :    grav   = gravit        ! gravitational constant
         141           0 :    rgas   = rair          ! gas constant for dry air
         142           0 :    rhoh2o = rhowtr        ! density of liquid water (STP)
         143             :
         144           0 :    limcnv = limcnv_in
         145             :
         146             :    ! Initialize free parameters for moist convective mass flux procedure
         147             :    ! cmftau - characteristic adjustment time scale
         148             :    ! c0     - rain water autoconversion coeff (1/m)
         149             :
         150           0 :    if (masterproc) then
         151           0 :       write(iulog,*) 'tuning parameters hk_conv: cmftau',cmftau
         152           0 :       write(iulog,*) 'tuning parameters hk_conv: c0',c0
         153             :    endif
         154           0 :    dzmin  = 0.0_r8           ! minimum cloud depth to precipitate (m)
         155           0 :    betamn = 0.10_r8          ! minimum overshoot parameter
         156             :
         157             :
         158           0 :    tpmax  = 1.50_r8          ! maximum acceptable t perturbation (deg C)
         159           0 :    shpmax = 1.50e-3_r8       ! maximum acceptable q perturbation (g/g)
         160           0 :    rlxclm = .true.        ! logical variable to specify that relaxation
         161             : !                                time scale should applied to column as
         162             : !                                opposed to triplets individually
         163             : !
         164             : ! Initialize miscellaneous (frequently used) constants
         165             : !
         166           0 :    rhlat  = 1.0_r8/hlat      ! reciprocal latent heat of vaporization
         167           0 :    rcp    = 1.0_r8/cp        ! reciprocal specific heat of dry air
         168           0 :    rgrav  = 1.0_r8/grav      ! reciprocal gravitational constant
         169             : !
         170             : ! Initialize diagnostic location information for moist convection scheme
         171             : !
         172           0 :    iloc   = 1             ! longitude point for diagnostic info
         173           0 :    jloc   = 1             ! latitude  point for diagnostic info
         174           0 :    nsloc  = 1             ! nstep value at which to begin diagnostics
         175             : !
         176             : ! Initialize other miscellaneous parameters
         177             : !
         178           0 :    tiny   = 1.0e-36_r8       ! arbitrary small number (scalar transport)
         179           0 :    eps    = 1.0e-13_r8       ! convergence criteria (machine dependent)
         180             : !
         181           0 :    return
         182             : end subroutine mfinti
         183             :
         184           0 : subroutine cmfmca(lchnk   ,ncol    , &
         185             :                   nstep   ,ztodt     ,pmid    ,pdel    , &
         186             :                   rpdel   ,zm      ,tpert   ,qpert   ,phis    , &
         187             :                   pblh    ,t       ,q       ,cmfdt   ,dq      , &
         188             :                   cmfmc   ,cmfdqr  ,cmfsl   ,cmflq   ,precc   , &
         189             :                   qc      ,cnt     ,cnb     ,icwmr   ,rliq    , &
         190             :                   pmiddry ,pdeldry ,rpdeldry)
         191             : !-----------------------------------------------------------------------
         192             : !
         193             : ! Purpose:
         194             : ! Moist convective mass flux procedure:
         195             : !
         196             : ! Method:
         197             : ! If stratification is unstable to nonentraining parcel ascent,
         198             : ! complete an adjustment making successive use of a simple cloud model
         199             : ! consisting of three layers (sometimes referred to as a triplet)
         200             : !
         201             : ! Code generalized to allow specification of parcel ("updraft")
         202             : ! properties, as well as convective transport of an arbitrary
         203             : ! number of passive constituents (see q array).  The code
         204             : ! is written so the water vapor field is passed independently
         205             : ! in the calling list from the block of other transported
         206             : ! constituents, even though as currently designed, it is the
         207             : ! first component in the constituents field.
         208             : !
         209             : ! Author: J. Hack
         210             : !
         211             : ! BAB: changed code to report tendencies in cmfdt and dq, instead of
         212             : ! updating profiles. Cmfdq contains water only, made it a local variable
         213             : ! made dq (all constituents) the argument.
         214             : !
         215             : !-----------------------------------------------------------------------
         216             :
         217             : !#######################################################################
         218             : !#                                                                     #
         219             : !# Debugging blocks are marked this way for easy identification        #
         220             : !#                                                                     #
         221             : !#######################################################################
         222             :    use constituents,  only: pcnst
         223             :    use constituents,    only: cnst_get_type_byind
         224             :    use ppgrid,    only: pcols, pver, pverp
         225             : #if ( defined DIAGNS )
         226             :    use phys_grid, only: get_lat_all_p, get_lon_all_p
         227             : #endif
         228             :    use wv_saturation, only: qsat
         229             :
         230             :    real(r8) ssfac               ! supersaturation bound (detrained air)
         231             :    parameter (ssfac = 1.001_r8)
         232             :
         233             : !------------------------------Arguments--------------------------------
         234             : !
         235             : ! Input arguments
         236             : !
         237             :    integer, intent(in) :: lchnk                ! chunk identifier
         238             :    integer, intent(in) :: ncol                 ! number of atmospheric columns
         239             :    integer, intent(in) :: nstep                ! current time step index
         240             :
         241             :    real(r8), intent(in) :: ztodt               ! 2 delta-t (seconds)
         242             :    real(r8), intent(in) :: pmid(pcols,pver)    ! pressure
         243             :    real(r8), intent(in) :: pdel(pcols,pver)    ! delta-p
         244             :    real(r8), intent(in) :: pmiddry(pcols,pver)    ! pressure
         245             :    real(r8), intent(in) :: pdeldry(pcols,pver)    ! delta-p
         246             :    real(r8), intent(in) :: rpdel(pcols,pver)   ! 1./pdel
         247             :    real(r8), intent(in) :: rpdeldry(pcols,pver)   ! 1./pdel
         248             :    real(r8), intent(in) :: zm(pcols,pver)      ! height abv sfc at midpoints
         249             :    real(r8), intent(in) :: tpert(pcols)        ! PBL perturbation theta
         250             :    real(r8), intent(in) :: qpert(pcols,pcnst)  ! PBL perturbation specific humidity
         251             :    real(r8), intent(in) :: phis(pcols)         ! surface geopotential
         252             :    real(r8), intent(in) :: pblh(pcols)         ! PBL height (provided by PBL routine)
         253             :    real(r8), intent(in) :: t(pcols,pver)       ! temperature (t bar)
         254             :    real(r8), intent(in) :: q(pcols,pver,pcnst) ! specific humidity (sh bar)
         255             : !
         256             : ! Output arguments
         257             : !
         258             :    real(r8), intent(out) :: cmfdt(pcols,pver)   ! dt/dt due to moist convection
         259             :    real(r8), intent(out) :: cmfmc(pcols,pverp)  ! moist convection cloud mass flux
         260             :    real(r8), intent(out) :: cmfdqr(pcols,pver)  ! dq/dt due to convective rainout
         261             :    real(r8), intent(out) :: cmfsl(pcols,pver )  ! convective lw static energy flux
         262             :    real(r8), intent(out) :: cmflq(pcols,pver )  ! convective total water flux
         263             :    real(r8), intent(out) :: precc(pcols)        ! convective precipitation rate
         264             : ! JJH mod to explicitly export cloud water
         265             :    real(r8), intent(out) :: qc(pcols,pver)      ! dq/dt due to export of cloud water
         266             :    real(r8), intent(out) :: cnt(pcols)          ! top level of convective activity
         267             :    real(r8), intent(out) :: cnb(pcols)          ! bottom level of convective activity
         268             :    real(r8), intent(out) :: dq(pcols,pver,pcnst) ! constituent tendencies
         269             :    real(r8), intent(out) :: icwmr(pcols,pver)
         270             :    real(r8), intent(out) :: rliq(pcols)
         271             : !
         272             : !---------------------------Local workspace-----------------------------
         273             : !
         274             :    real(r8) pm(pcols,pver)    ! pressure
         275             :    real(r8) pd(pcols,pver)    ! delta-p
         276             :    real(r8) rpd(pcols,pver)   ! 1./pdel
         277             :
         278             :    real(r8) cmfdq(pcols,pver)   ! dq/dt due to moist convection
         279             :    real(r8) gam(pcols,pver)     ! 1/cp (d(qsat)/dT)
         280             :    real(r8) sb(pcols,pver)      ! dry static energy (s bar)
         281             :    real(r8) hb(pcols,pver)      ! moist static energy (h bar)
         282             :    real(r8) shbs(pcols,pver)    ! sat. specific humidity (sh bar star)
         283             :    real(r8) hbs(pcols,pver)     ! sat. moist static energy (h bar star)
         284             :    real(r8) shbh(pcols,pverp)   ! specific humidity on interfaces
         285             :    real(r8) sbh(pcols,pverp)    ! s bar on interfaces
         286             :    real(r8) hbh(pcols,pverp)    ! h bar on interfaces
         287             :    real(r8) cmrh(pcols,pverp)   ! interface constituent mixing ratio
         288             :    real(r8) prec(pcols)         ! instantaneous total precipitation
         289             :    real(r8) dzcld(pcols)        ! depth of convective layer (m)
         290             :    real(r8) beta(pcols)         ! overshoot parameter (fraction)
         291             :    real(r8) betamx(pcols)       ! local maximum on overshoot
         292             :    real(r8) eta(pcols)          ! convective mass flux (kg/m^2 s)
         293             :    real(r8) etagdt(pcols)       ! eta*grav*dt
         294             :    real(r8) cldwtr(pcols)       ! cloud water (mass)
         295             :    real(r8) rnwtr(pcols)        ! rain water  (mass)
         296             : !  JJH extension to facilitate export of cloud liquid water
         297             :    real(r8) totcond(pcols)      ! total condensate; mix of precip and cloud water (mass)
         298             :    real(r8) sc  (pcols)         ! dry static energy   ("in-cloud")
         299             :    real(r8) shc (pcols)         ! specific humidity   ("in-cloud")
         300             :    real(r8) hc  (pcols)         ! moist static energy ("in-cloud")
         301             :    real(r8) cmrc(pcols)         ! constituent mix rat ("in-cloud")
         302             :    real(r8) dq1(pcols)          ! shb  convective change (lower lvl)
         303             :    real(r8) dq2(pcols)          ! shb  convective change (mid level)
         304             :    real(r8) dq3(pcols)          ! shb  convective change (upper lvl)
         305             :    real(r8) ds1(pcols)          ! sb   convective change (lower lvl)
         306             :    real(r8) ds2(pcols)          ! sb   convective change (mid level)
         307             :    real(r8) ds3(pcols)          ! sb   convective change (upper lvl)
         308             :    real(r8) dcmr1(pcols)        ! q convective change (lower lvl)
         309             :    real(r8) dcmr2(pcols)        ! q convective change (mid level)
         310             :    real(r8) dcmr3(pcols)        ! q convective change (upper lvl)
         311             :    real(r8) estemp(pcols,pver)  ! saturation vapor pressure (scratch)
         312             :    real(r8) vtemp1(2*pcols)     ! intermediate scratch vector
         313             :    real(r8) vtemp2(2*pcols)     ! intermediate scratch vector
         314             :    real(r8) vtemp3(2*pcols)     ! intermediate scratch vector
         315             :    real(r8) vtemp4(2*pcols)     ! intermediate scratch vector
         316             :    real(r8) vtemp5(2*pcols)     ! intermediate scratch vector
         317             :    integer indx1(pcols)     ! longitude indices for condition true
         318             :    logical etagt0           ! true if eta > 0.0
         319             :    real(r8) sh1                 ! dummy arg in qhalf statement func.
         320             :    real(r8) sh2                 ! dummy arg in qhalf statement func.
         321             :    real(r8) shbs1               ! dummy arg in qhalf statement func.
         322             :    real(r8) shbs2               ! dummy arg in qhalf statement func.
         323             :    real(r8) cats                ! modified characteristic adj. time
         324             :    real(r8) rtdt                ! 1./ztodt
         325             :    real(r8) qprime              ! modified specific humidity pert.
         326             :    real(r8) tprime              ! modified thermal perturbation
         327             :    real(r8) pblhgt              ! bounded pbl height (max[pblh,1m])
         328             :    real(r8) fac1                ! intermediate scratch variable
         329             :    real(r8) shprme              ! intermediate specific humidity pert.
         330             :    real(r8) qsattp              ! sat mix rat for thermally pert PBL parcels
         331             :    real(r8) dz                  ! local layer depth
         332             :    real(r8) temp1               ! intermediate scratch variable
         333             :    real(r8) b1                  ! bouyancy measure in detrainment lvl
         334             :    real(r8) b2                  ! bouyancy measure in condensation lvl
         335             :    real(r8) temp2               ! intermediate scratch variable
         336             :    real(r8) temp3               ! intermediate scratch variable
         337             :    real(r8) g                   ! bounded vertical gradient of hb
         338             :    real(r8) tmass               ! total mass available for convective exch
         339             :    real(r8) denom               ! intermediate scratch variable
         340             :    real(r8) qtest1              ! used in negative q test (middle lvl)
         341             :    real(r8) qtest2              ! used in negative q test (lower lvl)
         342             :    real(r8) fslkp               ! flux lw static energy (bot interface)
         343             :    real(r8) fslkm               ! flux lw static energy (top interface)
         344             :    real(r8) fqlkp               ! flux total water (bottom interface)
         345             :    real(r8) fqlkm               ! flux total water (top interface)
         346             :    real(r8) botflx              ! bottom constituent mixing ratio flux
         347             :    real(r8) topflx              ! top constituent mixing ratio flux
         348             :    real(r8) efac1               ! ratio q to convectively induced chg (btm lvl)
         349             :    real(r8) efac2               ! ratio q to convectively induced chg (mid lvl)
         350             :    real(r8) efac3               ! ratio q to convectively induced chg (top lvl)
         351             :    real(r8) tb(pcols,pver)      ! working storage for temp (t bar)
         352             :    real(r8) shb(pcols,pver)     ! working storage for spec hum (sh bar)
         353             :    real(r8) adjfac              ! adjustment factor (relaxation related)
         354             :    real(r8) rktp
         355             :    real(r8) rk
         356             : #if ( defined DIAGNS )
         357             : !
         358             : !  Following 7 real variables are used in diagnostics calculations
         359             : !
         360             :    real(r8) rh                  ! relative humidity
         361             :    real(r8) es                  ! sat vapor pressure
         362             :    real(r8) hsum1               ! moist static energy integral
         363             :    real(r8) qsum1               ! total water integral
         364             :    real(r8) hsum2               ! final moist static energy integral
         365             :    real(r8) qsum2               ! final total water integral
         366             :    real(r8) fac                 ! intermediate scratch variable
         367             : #endif
         368             :    integer i,k              ! longitude, level indices
         369             :    integer ii               ! index on "gathered" vectors
         370             :    integer len1             ! vector length of "gathered" vectors
         371             :    integer m                ! constituent index
         372             :    integer ktp              ! tmp indx used to track top of convective layer
         373             : #if ( defined DIAGNS )
         374             :    integer n                ! vertical index     (diagnostics)
         375             :    integer kp               ! vertical index     (diagnostics)
         376             :    integer kpp              ! index offset, kp+1 (diagnostics)
         377             :    integer kpm1             ! index offset, kp-1 (diagnostics)
         378             :    integer lat(pcols)       ! latitude indices
         379             :    integer lon(pcols)       ! longitude indices
         380             : #endif
         381             : !
         382             : !---------------------------Statement functions-------------------------
         383             : !
         384             :    real(r8) qhalf
         385             :    qhalf(sh1,sh2,shbs1,shbs2) = min(max(sh1,sh2),(shbs2*sh1 + shbs1*sh2)/(shbs1+shbs2))
         386             : !
         387             : !-----------------------------------------------------------------------
         388             :
         389             : !** BAB initialize output tendencies here
         390             : !       copy q to dq; use dq below for passive tracer transport
         391           0 :    cmfdt(:ncol,:)  = 0._r8
         392           0 :    cmfdq(:ncol,:)  = 0._r8
         393           0 :    dq(:ncol,:,2:)  = q(:ncol,:,2:)
         394           0 :    cmfmc(:ncol,:)  = 0._r8
         395           0 :    cmfdqr(:ncol,:) = 0._r8
         396           0 :    cmfsl(:ncol,:)  = 0._r8
         397           0 :    cmflq(:ncol,:)  = 0._r8
         398           0 :    qc(:ncol,:)     = 0._r8
         399           0 :    rliq(:ncol)     = 0._r8
         400             : !
         401             : #if ( defined DIAGNS )
         402             : ! Determine chunk latitudes and longitudes
         403             :    call get_lat_all_p(lchnk, ncol, lat)
         404             :    call get_lon_all_p(lchnk, ncol, lon)
         405             : #endif
         406             : !
         407             : ! Ensure that characteristic adjustment time scale (cmftau) assumed
         408             : ! in estimate of eta isn't smaller than model time scale (ztodt)
         409             : ! The time over which the convection is assumed to act (the adjustment
         410             : ! time scale) can be applied with each application of the three-level
         411             : ! cloud model, or applied to the column tendencies after a "hard"
         412             : ! adjustment (i.e., on a 2-delta t time scale) is evaluated
         413             : !
         414           0 :    if (rlxclm) then
         415           0 :       cats   = ztodt             ! relaxation applied to column
         416           0 :       adjfac = ztodt/(max(ztodt,cmftau))
         417             :    else
         418           0 :       cats   = max(ztodt,cmftau) ! relaxation applied to triplet
         419           0 :       adjfac = 1.0_r8
         420             :    endif
         421           0 :    rtdt = 1.0_r8/ztodt
         422             : !
         423             : ! Move temperature and moisture into working storage
         424             : !
         425           0 :    do k=limcnv,pver
         426           0 :       do i=1,ncol
         427           0 :          tb (i,k) = t(i,k)
         428           0 :          shb(i,k) = q(i,k,1)
         429             :       end do
         430             :    end do
         431           0 :    do k=1,pver
         432           0 :       do i=1,ncol
         433           0 :          icwmr(i,k) = 0._r8
         434             :       end do
         435             :    end do
         436             : !
         437             : ! Compute sb,hb,shbs,hbs
         438             : !
         439           0 :    do k = limcnv,pver
         440             :       call qsat(tb(1:ncol,k), pmid(1:ncol,k), &
         441             :            estemp(1:ncol,k), shbs(1:ncol,k), ncol, &
         442           0 :            gam=gam(1:ncol,k))
         443             :    end do
         444             : !
         445           0 :    do k=limcnv,pver
         446           0 :       do i=1,ncol
         447           0 :          sb (i,k) = cp*tb(i,k) + zm(i,k)*grav + phis(i)
         448           0 :          hb (i,k) = sb(i,k) + hlat*shb(i,k)
         449           0 :          hbs(i,k) = sb(i,k) + hlat*shbs(i,k)
         450             :       end do
         451             :    end do
         452             : !
         453             : ! Compute sbh, shbh
         454             : !
         455           0 :    do k=limcnv+1,pver
         456           0 :       do i=1,ncol
         457           0 :          sbh (i,k) = 0.5_r8*(sb(i,k-1) + sb(i,k))
         458           0 :          shbh(i,k) = qhalf(shb(i,k-1),shb(i,k),shbs(i,k-1),shbs(i,k))
         459           0 :          hbh (i,k) = sbh(i,k) + hlat*shbh(i,k)
         460             :       end do
         461             :    end do
         462             : !
         463             : ! Specify properties at top of model (not used, but filling anyway)
         464             : !
         465           0 :    do i=1,ncol
         466           0 :       sbh (i,limcnv) = sb(i,limcnv)
         467           0 :       shbh(i,limcnv) = shb(i,limcnv)
         468           0 :       hbh (i,limcnv) = hb(i,limcnv)
         469             :    end do
         470             : !
         471             : ! Zero vertically independent control, tendency & diagnostic arrays
         472             : !
         473           0 :    do i=1,ncol
         474           0 :       prec(i)  = 0.0_r8
         475           0 :       dzcld(i) = 0.0_r8
         476           0 :       cnb(i)   = 0.0_r8
         477           0 :       cnt(i)   = real(pver+1,r8)
         478             :    end do
         479             : #if ( defined DIAGNS )
         480             : !#######################################################################
         481             : !#                                                                     #
         482             : !#    output initial thermodynamic profile if debug diagnostics        #
         483             : !#                                                                     #
         484             :    do i=1,ncol
         485             :      if ((lat(i).eq.jloc) .and. (lon(i).eq.iloc) &
         486             :          .and. (nstep.ge.nsloc)) then
         487             : !#                                                                     #
         488             : !#       approximate vertical integral of moist static energy          #
         489             : !#       and total preciptable water                                   #
         490             : !#                                                                     #
         491             :       hsum1 = 0.0_r8
         492             :       qsum1 = 0.0_r8
         493             :       do k=limcnv,pver
         494             :          hsum1 = hsum1 + pdel(i,k)*rgrav*hb(i,k)
         495             :          qsum1 = qsum1 + pdel(i,k)*rgrav*shb(i,k)
         496             :       end do
         497             : !#                                                                     #
         498             :       write(iulog,8010)
         499             :       fac = grav*864._r8
         500             :       do k=limcnv,pver
         501             :          rh = shb(i,k)/shbs(i,k)
         502             :          write(iulog,8020) shbh(i,k),sbh(i,k),hbh(i,k),fac*cmfmc(i,k),cmfsl(i,k), cmflq(i,k)
         503             :          write(iulog,8040) tb(i,k),shb(i,k),rh,sb(i,k),hb(i,k),hbs(i,k),ztodt*cmfdt(i,k), &
         504             :                        ztodt*cmfdq(i,k),ztodt*cmfdqr(i,k)
         505             :       end do
         506             :       write(iulog, 8000) prec(i)
         507             :      end if
         508             :    enddo
         509             : #endif
         510             : !#                                                                     #
         511             : !#                                                                     #
         512             : !#######################################################################
         513             : !
         514             : ! Begin moist convective mass flux adjustment procedure.
         515             : ! Formalism ensures that negative cloud liquid water can never occur
         516             : !
         517           0 :    do 70 k=pver-1,limcnv+1,-1
         518           0 :       do 10 i=1,ncol
         519           0 :          etagdt(i) = 0.0_r8
         520           0 :          eta   (i) = 0.0_r8
         521           0 :          beta  (i) = 0.0_r8
         522           0 :          ds1   (i) = 0.0_r8
         523           0 :          ds2   (i) = 0.0_r8
         524           0 :          ds3   (i) = 0.0_r8
         525           0 :          dq1   (i) = 0.0_r8
         526           0 :          dq2   (i) = 0.0_r8
         527           0 :          dq3   (i) = 0.0_r8
         528             : !
         529             : ! Specification of "cloud base" conditions
         530             : !
         531           0 :          qprime    = 0.0_r8
         532           0 :          tprime    = 0.0_r8
         533             : !
         534             : ! Assign tprime within the PBL to be proportional to the quantity
         535             : ! tpert (which will be bounded by tpmax), passed to this routine by
         536             : ! the PBL routine.  Don't allow perturbation to produce a dry
         537             : ! adiabatically unstable parcel.  Assign qprime within the PBL to be
         538             : ! an appropriately modified value of the quantity qpert (which will be
         539             : ! bounded by shpmax) passed to this routine by the PBL routine.  The
         540             : ! quantity qprime should be less than the local saturation value
         541             : ! (qsattp=qsat[t+tprime,p]).  In both cases, tpert and qpert are
         542             : ! linearly reduced toward zero as the PBL top is approached.
         543             : !
         544           0 :          pblhgt = max(pblh(i),1.0_r8)
         545           0 :          if ( (zm(i,k+1) <= pblhgt) .and. dzcld(i) == 0.0_r8 ) then
         546           0 :             fac1   = max(0.0_r8,1.0_r8-zm(i,k+1)/pblhgt)
         547           0 :             tprime = min(tpert(i),tpmax)*fac1
         548           0 :             qsattp = shbs(i,k+1) + cp*rhlat*gam(i,k+1)*tprime
         549           0 :             shprme = min(min(qpert(i,1),shpmax)*fac1,max(qsattp-shb(i,k+1),0.0_r8))
         550           0 :             qprime = max(qprime,shprme)
         551             :          else
         552             :             tprime = 0.0_r8
         553             :             qprime = 0.0_r8
         554             :          end if
         555             : !
         556             : ! Specify "updraft" (in-cloud) thermodynamic properties
         557             : !
         558           0 :          sc (i)    = sb (i,k+1) + cp*tprime
         559           0 :          shc(i)    = shb(i,k+1) + qprime
         560           0 :          hc (i)    = sc (i    ) + hlat*shc(i)
         561           0 :          vtemp4(i) = hc(i) - hbs(i,k)
         562           0 :          dz        = pdel(i,k)*rgas*tb(i,k)*rgrav/pmid(i,k)
         563           0 :          if (vtemp4(i) > 0.0_r8) then
         564           0 :             dzcld(i) = dzcld(i) + dz
         565             :          else
         566           0 :             dzcld(i) = 0.0_r8
         567             :          end if
         568           0 : 10       continue
         569             : #if ( defined DIAGNS )
         570             : !#######################################################################
         571             : !#                                                                     #
         572             : !#    output thermodynamic perturbation information                    #
         573             : !#                                                                     #
         574             :          do i=1,ncol
         575             :            if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
         576             :             write(iulog,8090) k+1,sc(iloc),shc(iloc),hc(iloc)
         577             :            end if
         578             :          enddo
         579             : !#                                                                     #
         580             : !#######################################################################
         581             : #endif
         582             : !
         583             : ! Check on moist convective instability
         584             : ! Build index vector of points where instability exists
         585             : !
         586             :          len1 = 0
         587           0 :          do i=1,ncol
         588           0 :             if (vtemp4(i) > 0.0_r8) then
         589           0 :                len1 = len1 + 1
         590           0 :                indx1(len1) = i
         591             :             end if
         592             :          end do
         593           0 :          if (len1 <= 0) go to 70
         594             : !
         595             : ! Current level just below top level => no overshoot
         596             : !
         597           0 :          if (k <= limcnv+1) then
         598           0 :             do ii=1,len1
         599           0 :                i = indx1(ii)
         600           0 :                temp1     = vtemp4(i)/(1.0_r8 + gam(i,k))
         601           0 :                cldwtr(i) = max(0.0_r8,(sb(i,k) - sc(i) + temp1))
         602           0 :                beta(i)   = 0.0_r8
         603           0 :                vtemp3(i) = (1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k))
         604             :             end do
         605             :          else
         606             : !
         607             : ! First guess at overshoot parameter using crude buoyancy closure
         608             : ! 10% overshoot assumed as a minimum and 1-c0*dz maximum to start
         609             : ! If pre-existing supersaturation in detrainment layer, beta=0
         610             : ! cldwtr is temporarily equal to hlat*l (l=> liquid water)
         611             : !
         612           0 :             do ii=1,len1
         613           0 :                i = indx1(ii)
         614           0 :                temp1     = vtemp4(i)/(1.0_r8 + gam(i,k))
         615           0 :                cldwtr(i) = max(0.0_r8,(sb(i,k)-sc(i)+temp1))
         616           0 :                betamx(i) = 1.0_r8 - c0*max(0.0_r8,(dzcld(i)-dzmin))
         617           0 :                b1        = (hc(i) - hbs(i,k-1))*pdel(i,k-1)
         618           0 :                b2        = (hc(i) - hbs(i,k  ))*pdel(i,k  )
         619           0 :                beta(i)   = max(betamn,min(betamx(i), 1.0_r8 + b1/b2))
         620           0 :                if (hbs(i,k-1) <= hb(i,k-1)) beta(i) = 0.0_r8
         621             : !
         622             : ! Bound maximum beta to ensure physically realistic solutions
         623             : !
         624             : ! First check constrains beta so that eta remains positive
         625             : ! (assuming that eta is already positive for beta equal zero)
         626             : !
         627           0 :                vtemp1(i) = -(hbh(i,k+1) - hc(i))*pdel(i,k)*rpdel(i,k+1)+ &
         628           0 :                            (1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k+1) + cldwtr(i))
         629           0 :                vtemp2(i) = (1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k))
         630           0 :                vtemp3(i) = vtemp2(i)
         631           0 :                if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._r8) then
         632           0 :                   betamx(i) = 0.99_r8*(vtemp1(i)/vtemp2(i))
         633           0 :                   beta(i)   = max(0.0_r8,min(betamx(i),beta(i)))
         634             :                end if
         635             :             end do
         636             : !
         637             : ! Second check involves supersaturation of "detrainment layer"
         638             : ! small amount of supersaturation acceptable (by ssfac factor)
         639             : !
         640           0 :             do ii=1,len1
         641           0 :                i = indx1(ii)
         642           0 :                if (hb(i,k-1) < hbs(i,k-1)) then
         643           0 :                   vtemp1(i) = vtemp1(i)*rpdel(i,k)
         644           0 :                   temp2 = gam(i,k-1)*(sbh(i,k) - sc(i) + cldwtr(i)) -  &
         645           0 :                           hbh(i,k) + hc(i) - sc(i) + sbh(i,k)
         646           0 :                   temp3 = vtemp3(i)*rpdel(i,k)
         647             :                   vtemp2(i) = (ztodt/cats)*(hc(i) - hbs(i,k))*temp2/ &
         648           0 :                               (pdel(i,k-1)*(hbs(i,k-1) - hb(i,k-1))) + temp3
         649           0 :                   if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._r8) then
         650           0 :                      betamx(i) = ssfac*(vtemp1(i)/vtemp2(i))
         651           0 :                      beta(i)   = max(0.0_r8,min(betamx(i),beta(i)))
         652             :                   end if
         653             :                else
         654           0 :                   beta(i) = 0.0_r8
         655             :                end if
         656             :             end do
         657             : !
         658             : ! Third check to avoid introducing 2 delta x thermodynamic
         659             : ! noise in the vertical ... constrain adjusted h (or theta e)
         660             : ! so that the adjustment doesn't contribute to "kinks" in h
         661             : !
         662           0 :             do ii=1,len1
         663           0 :                i = indx1(ii)
         664           0 :                g = min(0.0_r8,hb(i,k) - hb(i,k-1))
         665           0 :                temp1 = (hb(i,k) - hb(i,k-1) - g)*(cats/ztodt)/(hc(i) - hbs(i,k))
         666           0 :                vtemp1(i) = temp1*vtemp1(i) + (hc(i) - hbh(i,k+1))*rpdel(i,k)
         667           0 :                vtemp2(i) = temp1*vtemp3(i)*rpdel(i,k) + (hc(i) - hbh(i,k) - cldwtr(i))* &
         668           0 :                            (rpdel(i,k) + rpdel(i,k+1))
         669           0 :                if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._r8) then
         670           0 :                   if (vtemp2(i) /= 0.0_r8) then
         671           0 :                      betamx(i) = vtemp1(i)/vtemp2(i)
         672             :                   else
         673           0 :                      betamx(i) = 0.0_r8
         674             :                   end if
         675           0 :                   beta(i) = max(0.0_r8,min(betamx(i),beta(i)))
         676             :                end if
         677             :             end do
         678             :          end if
         679             : !
         680             : ! Calculate mass flux required for stabilization.
         681             : !
         682             : ! Ensure that the convective mass flux, eta, is positive by
         683             : ! setting negative values of eta to zero..
         684             : ! Ensure that estimated mass flux cannot move more than the
         685             : ! minimum of total mass contained in either layer k or layer k+1.
         686             : ! Also test for other pathological cases that result in non-
         687             : ! physical states and adjust eta accordingly.
         688             : !
         689           0 :          do ii=1,len1
         690           0 :             i = indx1(ii)
         691           0 :             beta(i) = max(0.0_r8,beta(i))
         692           0 :             temp1 = hc(i) - hbs(i,k)
         693           0 :             temp2 = ((1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k+1) + cldwtr(i)) - &
         694           0 :                       beta(i)*vtemp3(i))*rpdel(i,k) - (hbh(i,k+1) - hc(i))*rpdel(i,k+1)
         695           0 :             eta(i) = temp1/(temp2*grav*cats)
         696           0 :             tmass = min(pdel(i,k),pdel(i,k+1))*rgrav
         697           0 :             if (eta(i) > tmass*rtdt .or. eta(i) <= 0.0_r8) eta(i) = 0.0_r8
         698             : !
         699             : ! Check on negative q in top layer (bound beta)
         700             : !
         701           0 :             if (shc(i)-shbh(i,k) < 0.0_r8 .and. beta(i)*eta(i) /= 0.0_r8) then
         702           0 :                denom = eta(i)*grav*ztodt*(shc(i) - shbh(i,k))*rpdel(i,k-1)
         703           0 :                beta(i) = max(0.0_r8,min(-0.999_r8*shb(i,k-1)/denom,beta(i)))
         704             :             end if
         705             : !
         706             : ! Check on negative q in middle layer (zero eta)
         707             : !
         708             :             qtest1 = shb(i,k) + eta(i)*grav*ztodt*((shc(i) - shbh(i,k+1)) - &
         709             :                      (1.0_r8 - beta(i))*cldwtr(i)*rhlat - beta(i)*(shc(i) - shbh(i,k)))* &
         710           0 :                      rpdel(i,k)
         711           0 :             if (qtest1 <= 0.0_r8) eta(i) = 0.0_r8
         712             : !
         713             : ! Check on negative q in lower layer (bound eta)
         714             : !
         715           0 :             fac1 = -(shbh(i,k+1) - shc(i))*rpdel(i,k+1)
         716           0 :             qtest2 = shb(i,k+1) - eta(i)*grav*ztodt*fac1
         717           0 :             if (qtest2 < 0.0_r8) then
         718           0 :                eta(i) = 0.99_r8*shb(i,k+1)/(grav*ztodt*fac1)
         719             :             end if
         720           0 :             etagdt(i) = eta(i)*grav*ztodt
         721             :          end do
         722             : !
         723             : #if ( defined DIAGNS )
         724             : !#######################################################################
         725             : !#                                                                     #
         726             :          do i=1,ncol
         727             :            if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep >= nsloc)) then
         728             :             write(iulog,8080) beta(iloc), eta(iloc)
         729             :            end if
         730             :          enddo
         731             : !#                                                                     #
         732             : !#######################################################################
         733             : #endif
         734             : !
         735             : ! Calculate cloud water, rain water, and thermodynamic changes
         736             : !
         737           0 :          do 30 ii=1,len1
         738           0 :             i = indx1(ii)
         739           0 :             icwmr(i,k) = cldwtr(i)*rhlat
         740           0 :             cldwtr(i) = etagdt(i)*cldwtr(i)*rhlat*rgrav
         741             : ! JJH changes to facilitate export of cloud liquid water --------------------------------
         742           0 :             totcond(i) = (1.0_r8 - beta(i))*cldwtr(i)
         743           0 :             rnwtr(i) = min(totcond(i),c0*(dzcld(i)-dzmin)*cldwtr(i))
         744           0 :             ds1(i) = etagdt(i)*(sbh(i,k+1) - sc(i))*rpdel(i,k+1)
         745           0 :             dq1(i) = etagdt(i)*(shbh(i,k+1) - shc(i))*rpdel(i,k+1)
         746             :             ds2(i) = (etagdt(i)*(sc(i) - sbh(i,k+1)) +  &
         747           0 :                      hlat*grav*cldwtr(i) - beta(i)*etagdt(i)*(sc(i) - sbh(i,k)))*rpdel(i,k)
         748             : ! JJH change for export of cloud liquid water; must use total condensate
         749             : ! since rainwater no longer represents total condensate
         750             :             dq2(i) = (etagdt(i)*(shc(i) - shbh(i,k+1)) - grav*totcond(i) - beta(i)* &
         751           0 :                      etagdt(i)*(shc(i) - shbh(i,k)))*rpdel(i,k)
         752             :             ds3(i) = beta(i)*(etagdt(i)*(sc(i) - sbh(i,k)) - hlat*grav*cldwtr(i))* &
         753           0 :                      rpdel(i,k-1)
         754           0 :             dq3(i) = beta(i)*etagdt(i)*(shc(i) - shbh(i,k))*rpdel(i,k-1)
         755             : !
         756             : ! Isolate convective fluxes for later diagnostics
         757             : !
         758           0 :             fslkp = eta(i)*(sc(i) - sbh(i,k+1))
         759           0 :             fslkm = beta(i)*(eta(i)*(sc(i) - sbh(i,k)) - hlat*cldwtr(i)*rtdt)
         760           0 :             fqlkp = eta(i)*(shc(i) - shbh(i,k+1))
         761           0 :             fqlkm = beta(i)*eta(i)*(shc(i) - shbh(i,k))
         762             : !
         763             : ! Update thermodynamic profile (update sb, hb, & hbs later)
         764             : !
         765           0 :             tb (i,k+1) = tb(i,k+1)  + ds1(i)*rcp
         766           0 :             tb (i,k  ) = tb(i,k  )  + ds2(i)*rcp
         767           0 :             tb (i,k-1) = tb(i,k-1)  + ds3(i)*rcp
         768           0 :             shb(i,k+1) = shb(i,k+1) + dq1(i)
         769           0 :             shb(i,k  ) = shb(i,k  ) + dq2(i)
         770           0 :             shb(i,k-1) = shb(i,k-1) + dq3(i)
         771             : !
         772             : ! ** Update diagnostic information for final budget **
         773             : ! Tracking precipitation, temperature & specific humidity tendencies,
         774             : ! rainout term, convective mass flux, convective liquid
         775             : ! water static energy flux, and convective total water flux
         776             : ! The variable afac makes the necessary adjustment to the
         777             : ! diagnostic fluxes to account for adjustment time scale based on
         778             : ! how relaxation time scale is to be applied (column vs. triplet)
         779             : !
         780           0 :             prec(i)    = prec(i) + (rnwtr(i)/rhoh2o)*adjfac
         781             : !
         782             : ! The following variables have units of "units"/second
         783             : !
         784           0 :             cmfdt (i,k+1) = cmfdt (i,k+1) + ds1(i)*rtdt*adjfac
         785           0 :             cmfdt (i,k  ) = cmfdt (i,k  ) + ds2(i)*rtdt*adjfac
         786           0 :             cmfdt (i,k-1) = cmfdt (i,k-1) + ds3(i)*rtdt*adjfac
         787           0 :             cmfdq (i,k+1) = cmfdq (i,k+1) + dq1(i)*rtdt*adjfac
         788           0 :             cmfdq (i,k  ) = cmfdq (i,k  ) + dq2(i)*rtdt*adjfac
         789           0 :             cmfdq (i,k-1) = cmfdq (i,k-1) + dq3(i)*rtdt*adjfac
         790             : ! JJH changes to export cloud liquid water --------------------------------
         791           0 :             qc    (i,k  ) = (grav*(totcond(i)-rnwtr(i))*rpdel(i,k))*rtdt*adjfac
         792           0 :             cmfdqr(i,k  ) = cmfdqr(i,k  ) + (grav*rnwtr(i)*rpdel(i,k))*rtdt*adjfac
         793           0 :             cmfmc (i,k+1) = cmfmc (i,k+1) + eta(i)*adjfac
         794           0 :             cmfmc (i,k  ) = cmfmc (i,k  ) + beta(i)*eta(i)*adjfac
         795             : !
         796             : ! The following variables have units of w/m**2
         797             : !
         798           0 :             cmfsl (i,k+1) = cmfsl (i,k+1) + fslkp*adjfac
         799           0 :             cmfsl (i,k  ) = cmfsl (i,k  ) + fslkm*adjfac
         800           0 :             cmflq (i,k+1) = cmflq (i,k+1) + hlat*fqlkp*adjfac
         801           0 :             cmflq (i,k  ) = cmflq (i,k  ) + hlat*fqlkm*adjfac
         802           0 : 30          continue
         803             : !
         804             : ! Next, convectively modify passive constituents
         805             : ! For now, when applying relaxation time scale to thermal fields after
         806             : ! entire column has undergone convective overturning, constituents will
         807             : ! be mixed using a "relaxed" value of the mass flux determined above
         808             : ! Although this will be inconsistant with the treatment of the thermal
         809             : ! fields, it's computationally much cheaper, no more-or-less justifiable,
         810             : ! and consistent with how the history tape mass fluxes would be used in
         811             : ! an off-line mode (i.e., using an off-line transport model)
         812             : !
         813           0 :             do 50 m=2,pcnst    ! note: indexing assumes water is first field
         814           0 :                if (cnst_get_type_byind(m).eq.'dry') then
         815           0 :                   pd(:ncol,:) = pdeldry(:ncol,:)
         816           0 :                   rpd(:ncol,:) = rpdeldry(:ncol,:)
         817           0 :                   pm(:ncol,:) = pmiddry(:ncol,:)
         818             :                else
         819           0 :                   pd(:ncol,:) = pdel(:ncol,:)
         820           0 :                   rpd(:ncol,:) = rpdel(:ncol,:)
         821           0 :                   pm(:ncol,:) = pmid(:ncol,:)
         822             :                endif
         823           0 :                do 40 ii=1,len1
         824           0 :                   i = indx1(ii)
         825             : !
         826             : ! If any of the reported values of the constituent is negative in
         827             : ! the three adjacent levels, nothing will be done to the profile
         828             : !
         829           0 :                   if ((dq(i,k+1,m) < 0.0_r8) .or. (dq(i,k,m) < 0.0_r8) .or. (dq(i,k-1,m) < 0.0_r8)) go to 40
         830             : !
         831             : ! Specify constituent interface values (linear interpolation)
         832             : !
         833           0 :                   cmrh(i,k  ) = 0.5_r8*(dq(i,k-1,m) + dq(i,k  ,m))
         834           0 :                   cmrh(i,k+1) = 0.5_r8*(dq(i,k  ,m) + dq(i,k+1,m))
         835             : !
         836             : ! Specify perturbation properties of constituents in PBL
         837             : !
         838           0 :                   pblhgt = max(pblh(i),1.0_r8)
         839           0 :                   if ( (zm(i,k+1) <= pblhgt) .and. dzcld(i) == 0.0_r8 ) then
         840           0 :                      fac1 = max(0.0_r8,1.0_r8-zm(i,k+1)/pblhgt)
         841           0 :                      cmrc(i) = dq(i,k+1,m) + qpert(i,m)*fac1
         842             :                   else
         843           0 :                      cmrc(i) = dq(i,k+1,m)
         844             :                   end if
         845             : !
         846             : ! Determine fluxes, flux divergence => changes due to convection
         847             : ! Logic must be included to avoid producing negative values. A bit
         848             : ! messy since there are no a priori assumptions about profiles.
         849             : ! Tendency is modified (reduced) when pending disaster detected.
         850             : !
         851           0 :                   botflx   = etagdt(i)*(cmrc(i) - cmrh(i,k+1))*adjfac
         852           0 :                   topflx   = beta(i)*etagdt(i)*(cmrc(i)-cmrh(i,k))*adjfac
         853           0 :                   dcmr1(i) = -botflx*rpd(i,k+1)
         854           0 :                   efac1    = 1.0_r8
         855           0 :                   efac2    = 1.0_r8
         856           0 :                   efac3    = 1.0_r8
         857             : !
         858           0 :                   if (dq(i,k+1,m)+dcmr1(i) < 0.0_r8) then
         859           0 :                      if ( abs(dcmr1(i)) > 1.e-300_r8 ) then
         860           0 :                         efac1 = max(tiny,abs(dq(i,k+1,m)/dcmr1(i)) - eps)
         861             :                      else
         862           0 :                         efac1 = tiny
         863             :                      endif
         864             :                   end if
         865             : !
         866           0 :                   if (efac1 == tiny .or. efac1 > 1.0_r8) efac1 = 0.0_r8
         867           0 :                   dcmr1(i) = -efac1*botflx*rpd(i,k+1)
         868           0 :                   dcmr2(i) = (efac1*botflx - topflx)*rpd(i,k)
         869             : !
         870           0 :                   if (dq(i,k,m)+dcmr2(i) < 0.0_r8) then
         871           0 :                      if ( abs(dcmr2(i)) > 1.e-300_r8 ) then
         872           0 :                         efac2 = max(tiny,abs(dq(i,k  ,m)/dcmr2(i)) - eps)
         873             :                      else
         874             :                         efac2 = tiny
         875             :                      endif
         876             :                   end if
         877             : !
         878           0 :                   if (efac2 == tiny .or. efac2 > 1.0_r8) efac2 = 0.0_r8
         879           0 :                   dcmr2(i) = (efac1*botflx - efac2*topflx)*rpd(i,k)
         880           0 :                   dcmr3(i) = efac2*topflx*rpd(i,k-1)
         881             : !
         882           0 :                   if ( (dq(i,k-1,m)+dcmr3(i) < 0.0_r8 ) ) then
         883           0 :                      if  ( abs(dcmr3(i)) > 1.e-300_r8 ) then
         884           0 :                         efac3 = max(tiny,abs(dq(i,k-1,m)/dcmr3(i)) - eps)
         885             :                      else
         886             :                         efac3 = tiny
         887             :                      endif
         888             :                   end if
         889             : !
         890           0 :                   if (efac3 == tiny .or. efac3 > 1.0_r8) efac3 = 0.0_r8
         891           0 :                   efac3    = min(efac2,efac3)
         892           0 :                   dcmr2(i) = (efac1*botflx - efac3*topflx)*rpd(i,k)
         893           0 :                   dcmr3(i) = efac3*topflx*rpd(i,k-1)
         894             : !
         895           0 :                   dq(i,k+1,m) = dq(i,k+1,m) + dcmr1(i)
         896           0 :                   dq(i,k  ,m) = dq(i,k  ,m) + dcmr2(i)
         897           0 :                   dq(i,k-1,m) = dq(i,k-1,m) + dcmr3(i)
         898           0 : 40                continue
         899           0 : 50                continue                ! end of m=2,pcnst loop
         900             : !
         901             : ! Constituent modifications complete
         902             : !
         903           0 :                   if (k == limcnv+1) go to 60
         904             : !
         905             : ! Complete update of thermodynamic structure at integer levels
         906             : ! gather/scatter points that need new values of shbs and gamma
         907             : !
         908           0 :                   do ii=1,len1
         909           0 :                      i = indx1(ii)
         910           0 :                      vtemp1(ii     ) = tb(i,k)
         911           0 :                      vtemp1(ii+len1) = tb(i,k-1)
         912           0 :                      vtemp2(ii     ) = pmid(i,k)
         913           0 :                      vtemp2(ii+len1) = pmid(i,k-1)
         914             :                   end do
         915           0 :                   call qsat(vtemp1(1:2*len1), vtemp2(1:2*len1), &
         916           0 :                             vtemp5(1:2*len1), vtemp3(1:2*len1), 2*len1, gam=vtemp4(1:2*len1))
         917           0 :                   do ii=1,len1
         918           0 :                      i = indx1(ii)
         919           0 :                      shbs(i,k  ) = vtemp3(ii     )
         920           0 :                      shbs(i,k-1) = vtemp3(ii+len1)
         921           0 :                      gam(i,k  ) = vtemp4(ii     )
         922           0 :                      gam(i,k-1) = vtemp4(ii+len1)
         923           0 :                      sb (i,k  ) = sb(i,k  ) + ds2(i)
         924           0 :                      sb (i,k-1) = sb(i,k-1) + ds3(i)
         925           0 :                      hb (i,k  ) = sb(i,k  ) + hlat*shb(i,k  )
         926           0 :                      hb (i,k-1) = sb(i,k-1) + hlat*shb(i,k-1)
         927           0 :                      hbs(i,k  ) = sb(i,k  ) + hlat*shbs(i,k  )
         928           0 :                      hbs(i,k-1) = sb(i,k-1) + hlat*shbs(i,k-1)
         929             :                   end do
         930             : !
         931             : ! Update thermodynamic information at half (i.e., interface) levels
         932             : !
         933           0 :                   do ii=1,len1
         934           0 :                      i = indx1(ii)
         935           0 :                      sbh (i,k) = 0.5_r8*(sb(i,k) + sb(i,k-1))
         936           0 :                      shbh(i,k) = qhalf(shb(i,k-1),shb(i,k),shbs(i,k-1),shbs(i,k))
         937           0 :                      hbh (i,k) = sbh(i,k) + hlat*shbh(i,k)
         938           0 :                      sbh (i,k-1) = 0.5_r8*(sb(i,k-1) + sb(i,k-2))
         939           0 :                      shbh(i,k-1) = qhalf(shb(i,k-2),shb(i,k-1),shbs(i,k-2),shbs(i,k-1))
         940           0 :                      hbh (i,k-1) = sbh(i,k-1) + hlat*shbh(i,k-1)
         941             :                   end do
         942             : !
         943             : #if ( defined DIAGNS )
         944             : !#######################################################################
         945             : !#                                                                     #
         946             : !#    this update necessary, only if debugging diagnostics requested   #
         947             : !#                                                                     #
         948             :                   do i=1,ncol
         949             :                      if (lat(i) == jloc .and. nstep >= nsloc) then
         950             :                         call qsat(tb(i,k+1), pmid(i,k+1), &
         951             :                              es, shbs(i,k+1), gam=gam(i,k+1))
         952             :                         sb (i,k+1) = sb(i,k+1) + ds1(i)
         953             :                         hb (i,k+1) = sb(i,k+1) + hlat*shb(i,k+1)
         954             :                         hbs(i,k+1) = sb(i,k+1) + hlat*shbs(i,k+1)
         955             :                         kpp = k + 2
         956             :                         if (k+1 == pver) kpp = k + 1
         957             :                         do kp=k+1,kpp
         958             :                            kpm1 = kp-1
         959             :                            sbh(i,kp)  = 0.5_r8*(sb(i,kpm1) + sb(i,kp))
         960             :                            shbh(i,kp) = qhalf(shb(i,kpm1),shb(i,kp),shbs(i,kpm1),shbs(i,kp))
         961             :                            hbh(i,kp)  = sbh(i,kp) + hlat*shbh(i,kp)
         962             :                         end do
         963             :                      end if
         964             :                   end do
         965             : !#                                                                     #
         966             : !#          diagnostic output                                          #
         967             : !#                                                                     #
         968             :                   do i=1,ncol
         969             :                     if ((lat(i)== jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
         970             :                      write(iulog, 8060) k
         971             :                      fac = grav*864._r8
         972             :                      do n=limcnv,pver
         973             :                         rh  = shb(i,n)/shbs(i,n)
         974             :                         write(iulog,8020)shbh(i,n),sbh(i,n),hbh(i,n),fac*cmfmc(i,n), &
         975             :                                      cmfsl(i,n), cmflq(i,n)
         976             : !--------------write(iulog, 8050)
         977             : !--------------write(iulog, 8030) fac*cmfmc(i,n),cmfsl(i,n), cmflq(i,n)
         978             :                         write(iulog, 8040) tb(i,n),shb(i,n),rh,sb(i,n),hb(i,n), &
         979             :                                        hbs(i,n), ztodt*cmfdt(i,n),ztodt*cmfdq(i,n), &
         980             :                                        ztodt*cmfdqr(i,n)
         981             :                      end do
         982             :                      write(iulog, 8000) prec(i)
         983             :                     end if
         984             :                   end do
         985             : !#                                                                     #
         986             : !#                                                                     #
         987             : !#######################################################################
         988             : #endif
         989             : !
         990             : ! Ensure that dzcld is reset if convective mass flux zero
         991             : ! specify the current vertical extent of the convective activity
         992             : ! top of convective layer determined by size of overshoot param.
         993             : !
         994           0 : 60                do i=1,ncol
         995           0 :                      etagt0 = eta(i).gt.0.0_r8
         996           0 :                      if ( .not. etagt0) dzcld(i) = 0.0_r8
         997           0 :                      if (etagt0 .and. beta(i) > betamn) then
         998           0 :                         ktp = k-1
         999             :                      else
        1000             :                         ktp = k
        1001             :                      end if
        1002           0 :                      if (etagt0) then
        1003           0 :                         rk=k
        1004           0 :                         rktp=ktp
        1005           0 :                         cnt(i) = min(cnt(i),rktp)
        1006           0 :                         cnb(i) = max(cnb(i),rk)
        1007             :                      end if
        1008             :                   end do
        1009           0 : 70                continue                  ! end of k loop
        1010             : !
        1011             : ! ** apply final thermodynamic tendencies **
        1012             : !
        1013             : !**BAB don't update input profiles
        1014             : !!$                  do k=limcnv,pver
        1015             : !!$                     do i=1,ncol
        1016             : !!$                        t (i,k) = t (i,k) + cmfdt(i,k)*ztodt
        1017             : !!$                        q(i,k,1) = q(i,k,1) + cmfdq(i,k)*ztodt
        1018             : !!$                     end do
        1019             : !!$                  end do
        1020             : ! Set output q tendencies
        1021           0 :       dq(:ncol,:,1 ) = cmfdq(:ncol,:)
        1022           0 :       dq(:ncol,:,2:) = (dq(:ncol,:,2:) - q(:ncol,:,2:))/ztodt
        1023             : !
        1024             : ! Kludge to prevent cnb-cnt from being zero (in the event
        1025             : ! someone decides that they want to divide by this quantity)
        1026             : !
        1027           0 :                   do i=1,ncol
        1028           0 :                      if (cnb(i) /= 0.0_r8 .and. cnb(i) == cnt(i)) then
        1029           0 :                         cnt(i) = cnt(i) - 1.0_r8
        1030             :                      end if
        1031             :                   end do
        1032             : !
        1033           0 :                   do i=1,ncol
        1034           0 :                      precc(i) = prec(i)*rtdt
        1035             :                   end do
        1036             : !
        1037             : ! Compute reserved liquid (not yet in cldliq) for energy integrals.
        1038             : ! Treat rliq as flux out bottom, to be added back later.
        1039           0 :    do k = 1, pver
        1040           0 :       do i = 1, ncol
        1041           0 :          rliq(i) = rliq(i) + qc(i,k)*pdel(i,k)/grav
        1042             :       end do
        1043             :    end do
        1044           0 :    rliq(:ncol) = rliq(:ncol) /1000._r8
        1045             :
        1046             : #if ( defined DIAGNS )
        1047             : !#######################################################################
        1048             : !#                                                                     #
        1049             : !#    we're done ... show final result if debug diagnostics requested  #
        1050             : !#                                                                     #
        1051             :                   do i=1,ncol
        1052             :                     if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
        1053             :                      fac = grav*864._r8
        1054             :                      write(iulog, 8010)
        1055             :                      do k=limcnv,pver
        1056             :                         rh = shb(i,k)/shbs(i,k)
        1057             :                         write(iulog, 8020) shbh(i,k),sbh(i,k),hbh(i,k),fac*cmfmc(i,k), &
        1058             :                                        cmfsl(i,k), cmflq(i,k)
        1059             :                         write(iulog, 8040) tb(i,k),shb(i,k),rh   ,sb(i,k),hb(i,k), &
        1060             :                                        hbs(i,k), ztodt*cmfdt(i,k),ztodt*cmfdq(i,k), &
        1061             :                                        ztodt*cmfdqr(i,k)
        1062             :                      end do
        1063             :                      write(iulog, 8000) prec(i)
        1064             : !#                                                                     #
        1065             : !#       approximate vertical integral of moist static energy and      #
        1066             : !#       total preciptable water after adjustment and output changes   #
        1067             : !#                                                                     #
        1068             :                      hsum2 = 0.0_r8
        1069             :                      qsum2 = 0.0_r8
        1070             :                      do k=limcnv,pver
        1071             :                         hsum2 = hsum2 + pdel(i,k)*rgrav*hb(i,k)
        1072             :                         qsum2 = qsum2 + pdel(i,k)*rgrav*shb(i,k)
        1073             :                      end do
        1074             : !#                                                                     #
        1075             :                      write(iulog,8070) hsum1, hsum2, abs(hsum2-hsum1)/hsum2, &
        1076             :                                     qsum1, qsum2, abs(qsum2-qsum1)/qsum2
        1077             :                     end if
        1078             :                   enddo
        1079             : !#                                                                     #
        1080             : !#                                                                     #
        1081             : !#######################################################################
        1082             : #endif
        1083           0 :                   return                 ! we're all done ... return to calling procedure
        1084             : #if ( defined DIAGNS )
        1085             : !
        1086             : ! Formats
        1087             : !
        1088             : 8000              format(///,10x,'PREC = ',3pf12.6,/)
        1089             : 8010              format('1**        TB      SHB      RH       SB', &
        1090             :                         '       HB      HBS      CAH      CAM       PRECC ', &
        1091             :                         '     ETA      FSL       FLQ     **', /)
        1092             : 8020              format(' ----- ',     9x,3p,f7.3,2x,2p,     9x,-3p,f7.3,2x, &
        1093             :                         f7.3, 37x, 0p,2x,f8.2,  0p,2x,f8.2,2x,f8.2, ' ----- ')
        1094             : 8030              format(' ----- ',  0p,82x,f8.2,  0p,2x,f8.2,2x,f8.2, &
        1095             :                          ' ----- ')
        1096             : 8040              format(' - - - ',f7.3,2x,3p,f7.3,2x,2p,f7.3,2x,-3p,f7.3,2x, &
        1097             :                         f7.3, 2x,f8.3,2x,0p,f7.3,3p,2x,f7.3,2x,f7.3,30x, &
        1098             :                          ' - - - ')
        1099             : 8050              format(' ----- ',110x,' ----- ')
        1100             : 8060              format('1 K =>',  i4,/, &
        1101             :                            '           TB      SHB      RH       SB', &
        1102             :                            '       HB      HBS      CAH      CAM       PREC ', &
        1103             :                            '     ETA      FSL       FLQ', /)
        1104             : 8070              format(' VERTICALLY INTEGRATED MOIST STATIC ENERGY BEFORE, AFTER', &
        1105             :                         ' AND PERCENTAGE DIFFERENCE => ',1p,2e15.7,2x,2p,f7.3,/, &
        1106             :                         ' VERTICALLY INTEGRATED MOISTURE            BEFORE, AFTER' &,
        1107             :                         ' AND PERCENTAGE DIFFERENCE => ',1p,2e15_r8.7_r8,2x,2p,f7.3,/)
        1108             : 8080              format(' BETA, ETA => ', 1p,2e12.3)
        1109             : 8090              format (' k+1, sc, shc, hc => ', 1x, i2, 1p, 3e12.4)
        1110             : #endif
        1111             : !
        1112             : end subroutine cmfmca
        1113             : end module hk_conv

Generated by: LCOV version 1.14