phys/module_cu_camzm.F

   1 #define WRF_PORT
   2 #define MODAL_AERO
   3 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
   4 module module_cu_camzm
   5
   6 !---------------------------------------------------------------------------------
   7 ! Purpose:
   8 !
   9 ! Interface from Zhang-McFarlane convection scheme, includes evaporation of convective
  10 ! precip from the ZM scheme
  11 !
  12 ! Apr 2006: RBN: Code added to perform a dilute ascent for closure of the CM mass flux
  13 !                based on an entraining plume a la Raymond and Blythe (1992)
  14 !
  15 ! Author: Byron Boville, from code in tphysbc
  16 !
  17 !---------------------------------------------------------------------------------
  18   use shr_kind_mod,    only: r8 => shr_kind_r8
  19 #ifndef WRF_PORT
  20   use spmd_utils,      only: masterproc
  21   use ppgrid,          only: pcols, pver, pverp
  22 #else
  23   use module_cam_support, only: masterproc, pcols, pver, pverp
  24 #endif
  25   use cldwat,          only: cldwat_fice
  26   use physconst,       only: cpair, epsilo, gravit, latice, latvap, tmelt, rair, &
  27                              cpwv, cpliq, rh2o
  28 #ifndef WRF_PORT
  29   use abortutils,      only: endrun
  30   use cam_logfile,     only: iulog
  31 #else
  32   use module_cam_support,only: endrun, iulog
  33 #endif
  34
  35   implicit none
  36
  37   save
  38   private                         ! Make default type private to the module
  39 !
  40 ! PUBLIC: interfaces
  41 !
  42   public zmconv_readnl            ! read zmconv_nl namelist
  43   public zm_convi                 ! ZM schemea
  44   public zm_convr                 ! ZM schemea
  45   public zm_conv_evap             ! evaporation of precip from ZM schemea
  46   public convtran                 ! convective transport
  47   public momtran                  ! convective momentum transport
  48
  49 !
  50 ! Private data
  51 !
  52    real(r8), parameter :: unset_r8 = huge(1.0_r8)
  53    real(r8) :: zmconv_c0_lnd = unset_r8
  54    real(r8) :: zmconv_c0_ocn = unset_r8
  55    real(r8) :: zmconv_ke     = unset_r8
  56
  57    real(r8) rl         ! wg latent heat of vaporization.
  58    real(r8) cpres      ! specific heat at constant pressure in j/kg-degk.
  59    real(r8), parameter :: capelmt = 70._r8  ! threshold value for cape for deep convection.
  60    real(r8) :: ke           ! Tunable evaporation efficiency set from namelist input zmconv_ke
  61    real(r8) :: c0_lnd       ! set from namelist input zmconv_c0_lnd
  62    real(r8) :: c0_ocn       ! set from namelist input zmconv_c0_ocn
  63    real(r8) tau   ! convective time scale
  64    real(r8),parameter :: a = 21.656_r8
  65    real(r8),parameter :: b = 5418._r8
  66    real(r8),parameter :: c1 = 6.112_r8
  67    real(r8),parameter :: c2 = 17.67_r8
  68    real(r8),parameter :: c3 = 243.5_r8
  69    real(r8) :: tfreez
  70    real(r8) :: eps1
  71
  72
  73    logical :: no_deep_pbl ! default = .false.
  74                           ! no_deep_pbl = .true. eliminates deep convection entirely within PBL
  75
  76
  77 !moved from moistconvection.F90
  78    real(r8) :: rgrav       ! reciprocal of grav
  79    real(r8) :: rgas        ! gas constant for dry air
  80    real(r8) :: grav        ! = gravit
  81    real(r8) :: cp          ! = cpres = cpair
  82
  83    integer  limcnv       ! top interface level limit for convection
  84
  85    real(r8),parameter ::  tiedke_add = 0.5_r8
  86
  87 contains
  88
  89 subroutine zmconv_readnl(nlfile)
  90 #ifndef WRF_PORT
  91    use namelist_utils,  only: find_group_name
  92    use units,           only: getunit, freeunit
  93    use mpishorthand
  94 #endif
  95
  96    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
  97 #ifndef WRF_PORT
  98    ! Local variables
  99    integer :: unitn, ierr
 100    character(len=*), parameter :: subname = 'zmconv_readnl'
 101
 102    namelist /zmconv_nl/ zmconv_c0_lnd, zmconv_c0_ocn, zmconv_ke
 103    !-----------------------------------------------------------------------------
 104
 105    if (masterproc) then
 106       unitn = getunit()
 107       open( unitn, file=trim(nlfile), status='old' )
 108       call find_group_name(unitn, 'zmconv_nl', status=ierr)
 109       if (ierr == 0) then
 110          read(unitn, zmconv_nl, iostat=ierr)
 111          if (ierr /= 0) then
 112             call endrun(subname // ':: ERROR reading namelist')
 113          end if
 114       end if
 115       close(unitn)
 116       call freeunit(unitn)
 117
 118       ! set local variables
 119       c0_lnd = zmconv_c0_lnd
 120       c0_ocn = zmconv_c0_ocn
 121       ke = zmconv_ke
 122
 123    end if
 124
 125 #ifdef SPMD
 126    ! Broadcast namelist variables
 127    call mpibcast(c0_lnd,            1, mpir8,  0, mpicom)
 128    call mpibcast(c0_ocn,            1, mpir8,  0, mpicom)
 129    call mpibcast(ke,                1, mpir8,  0, mpicom)
 130 #endif
 131 #else
 132 ! WRF_PORT currently uses hard-wired values for the namelist input. The
 133 ! values could easily be setup to come from the Registry in the future.
 134 ! The hard-wired values are the defaults for the fv core. They should be
 135 ! verified by somebody knowledgable on the matter.
 136    c0_lnd = 0.0059D0
 137    c0_ocn = 0.0450D0
 138    ke     = 1.0E-6
 139 #endif
 140 end subroutine zmconv_readnl
 141
 142
 143 #ifndef WRF_PORT
 144    subroutine zm_convi(limcnv_in, no_deep_pbl_in)
 145    use dycore,       only: dycore_is, get_resolution
 146 #else
 147    subroutine zm_convi(DT,DX,limcnv_in, no_deep_pbl_in)
 148 #endif
 149
 150
 151    integer, intent(in)           :: limcnv_in       ! top interface level limit for convection
 152    logical, intent(in), optional :: no_deep_pbl_in  ! no_deep_pbl = .true. eliminates ZM convection entirely within PBL
 153 #ifdef WRF_PORT
 154    REAL ,    INTENT(IN)        :: DT, DX
 155 #endif
 156    real(r8) :: DTT
 157
 158
 159
 160    ! local variables
 161    character(len=32)   :: hgrid           ! horizontal grid specifier
 162
 163 #ifdef WRF_PORT
 164    real(r8) :: deltax     ! model delta x
 165    real(r8) :: ref_dx     ! reference delta x (CAM5 1.9 x 2.5 deg)
 166    real(r8) :: taumin     ! minimum allowable tau
 167    real(r8) :: taumax     ! maximum allowable tau
 168 #endif
 169
 170    ! Initialization of ZM constants
 171    limcnv = limcnv_in
 172    tfreez = tmelt
 173    eps1   = epsilo
 174    rl     = latvap
 175    cpres  = cpair
 176    rgrav  = 1.0_r8/gravit
 177    rgas   = rair
 178    grav   = gravit
 179    cp     = cpres
 180 #ifdef WRF_PORT
 181    !PMA sets parameters
 182    ref_dx = 275000._r8
 183    taumin = 600._r8
 184    taumax = 3600._r8
 185    deltax=dx
 186 #endif
 187
 188    if ( present(no_deep_pbl_in) )  then
 189       no_deep_pbl = no_deep_pbl_in
 190    else
 191       no_deep_pbl = .true.
 192    endif
 193
 194    DTT=DT
 195
 196    ! tau=4800. were used in canadian climate center. however, in echam3 t42,
 197    ! convection is too weak, thus adjusted to 2400.
 198 #ifndef WRF_PORT
 199    hgrid = get_resolution()
 200    tau = 3600._r8
 201 !  PMA: Standard 2-deg CAM5.1 uses 1-hr relaxation time scale which is too long for the
 202 !  mesoscale model to remove CAPE quickly. As a result, the model generates grid-scale
 203 !  storms and becomes unstable sometimes (Williamson, 2012). Hence, we sets tau as a function
 204 !  of the length of time step. The caveat of this is that the model behavior
 205 !  of precipitation will be different to CAM5. We acknowledge that reducing dynamical
 206 !  time step can also make a stable model simulation. March 12, 2013
 207 !
 208    tau = max(min(2._r8*DTT,1200._r8),120._r8)
 209 #else
 210
 211 !PMA: The 1.9x2.5 deg CAM5 uses tau = 3600s. With higher resolution and shorter timestep in CAM,
 212 !     tau is set to a smaller value. Williamson (2012) found that without reducing tau in a high-res
 213 !     simulation, 'grid-scale storms' can form because the ZM scheme fails to remove CAPE quickly enough
 214 !     so that the resolved scheme (micro+macro) has to work very hard. WRF usually is run at much higher
 215 !     resolution than CAM and this grid-scale storm problem is apparent, which can cause the model to
 216 !     crash. Hence, we make tau to be a function of spatial resolution.
 217 !     For WRF, one needs to make dt very small if tau=3600s for the model to run stably.
 218
 219 !   tau = 3600._r8
 220
 221 !PMA formula: tau = 3600s when dx = 275km ~ 2.5deg
 222 !             tau =  600s when dx = 10km
 223
 224 !  tau = 11.32_r8 * deltax + 489.68_r8
 225
 226 !PMA adopts PJR formula
 227
 228    tau = max(taumin, taumax*min(1._r8,deltax/ref_dx))
 229
 230 #endif
 231    if ( masterproc ) then
 232       write(iulog,*) 'delta X =',deltax
 233 #ifdef WRF_PORT
 234       call wrf_debug(1,iulog)
 235 #endif
 236       write(iulog,*) 'Convective relaxation time scale (tau) is a tunable parameter in CAM and is a function of spatial resolution.'
 237       write(iulog,*) 'Users are encouraged to consult with the PNNL WRF-CAM5 development team if they want to change tau.'
 238 #ifdef WRF_PORT
 239       call wrf_debug(1,iulog)
 240 #endif
 241       write(iulog,*) 'tuning parameters zm_convi: tau',tau
 242 #ifdef WRF_PORT
 243       call wrf_message(iulog)
 244 #endif
 245       write(iulog,*) 'Standard 2-deg CAM5.1 sets tau=3600s and is reduced for the mesoscale model WRF'
 246 #ifdef WRF_PORT
 247       call wrf_debug(1,iulog)
 248 #endif
 249       write(iulog,*) 'tuning parameters zm_convi: c0_lnd',c0_lnd, ', c0_ocn', c0_ocn
 250 #ifdef WRF_PORT
 251       call wrf_debug(1,iulog)
 252 #endif
 253       write(iulog,*) 'tuning parameters zm_convi: ke',ke
 254 #ifdef WRF_PORT
 255       call wrf_debug(1,iulog)
 256 #endif
 257       write(iulog,*) 'tuning parameters zm_convi: no_deep_pbl',no_deep_pbl
 258 #ifdef WRF_PORT
 259       call wrf_debug(1,iulog)
 260 #endif
 261    endif
 262
 263    if (masterproc) write(iulog,*)'**** ZM: DILUTE Buoyancy Calculation ****'
 264 #ifdef WRF_PORT
 265       call wrf_debug(1,iulog)
 266 #endif
 267
 268 end subroutine zm_convi
 269
 270
 271
 272 subroutine zm_convr(lchnk   ,ncol    , &
 273                     t       ,qh      ,prec    ,jctop   ,jcbot   , &
 274                     pblh    ,zm      ,geos    ,zi      ,qtnd    , &
 275                     heat    ,pap     ,paph    ,dpp     , &
 276                     delt    ,mcon    ,cme     ,cape    , &
 277                     tpert   ,dlf     ,pflx    ,zdu     ,rprd    , &
 278                     mu      ,md      ,du      ,eu      ,ed      , &
 279                     dp      ,dsubcld ,jt      ,maxg    ,ideep   , &
 280                     lengath ,ql      ,rliq    ,landfrac)
 281 !-----------------------------------------------------------------------
 282 !
 283 ! Purpose:
 284 ! Main driver for zhang-mcfarlane convection scheme
 285 !
 286 ! Method:
 287 ! performs deep convective adjustment based on mass-flux closure
 288 ! algorithm.
 289 !
 290 ! Author:guang jun zhang, m.lazare, n.mcfarlane. CAM Contact: P. Rasch
 291 !
 292 ! This is contributed code not fully standardized by the CAM core group.
 293 ! All variables have been typed, where most are identified in comments
 294 ! The current procedure will be reimplemented in a subsequent version
 295 ! of the CAM where it will include a more straightforward formulation
 296 ! and will make use of the standard CAM nomenclature
 297 !
 298 !-----------------------------------------------------------------------
 299 #ifdef WRF_PORT
 300   use module_cam_support, only: pcnst =>pcnst_runtime
 301 #else
 302    use constituents, only: pcnst
 303    use phys_control, only: cam_physpkg_is
 304 #endif
 305
 306 !
 307 ! ************************ index of variables **********************
 308 !
 309 !  wg * alpha    array of vertical differencing used (=1. for upstream).
 310 !  w  * cape     convective available potential energy.
 311 !  wg * capeg    gathered convective available potential energy.
 312 !  c  * capelmt  threshold value for cape for deep convection.
 313 !  ic  * cpres    specific heat at constant pressure in j/kg-degk.
 314 !  i  * dpp
 315 !  ic  * delt     length of model time-step in seconds.
 316 !  wg * dp       layer thickness in mbs (between upper/lower interface).
 317 !  wg * dqdt     mixing ratio tendency at gathered points.
 318 !  wg * dsdt     dry static energy ("temp") tendency at gathered points.
 319 !  wg * dudt     u-wind tendency at gathered points.
 320 !  wg * dvdt     v-wind tendency at gathered points.
 321 !  wg * dsubcld  layer thickness in mbs between lcl and maxi.
 322 !  ic  * grav     acceleration due to gravity in m/sec2.
 323 !  wg * du       detrainment in updraft. specified in mid-layer
 324 !  wg * ed       entrainment in downdraft.
 325 !  wg * eu       entrainment in updraft.
 326 !  wg * hmn      moist static energy.
 327 !  wg * hsat     saturated moist static energy.
 328 !  w  * ideep    holds position of gathered points vs longitude index.
 329 !  ic  * pver     number of model levels.
 330 !  wg * j0       detrainment initiation level index.
 331 !  wg * jd       downdraft   initiation level index.
 332 !  ic  * jlatpr   gaussian latitude index for printing grids (if needed).
 333 !  wg * jt       top  level index of deep cumulus convection.
 334 !  w  * lcl      base level index of deep cumulus convection.
 335 !  wg * lclg     gathered values of lcl.
 336 !  w  * lel      index of highest theoretical convective plume.
 337 !  wg * lelg     gathered values of lel.
 338 !  w  * lon      index of onset level for deep convection.
 339 !  w  * maxi     index of level with largest moist static energy.
 340 !  wg * maxg     gathered values of maxi.
 341 !  wg * mb       cloud base mass flux.
 342 !  wg * mc       net upward (scaled by mb) cloud mass flux.
 343 !  wg * md       downward cloud mass flux (positive up).
 344 !  wg * mu       upward   cloud mass flux (positive up). specified
 345 !                at interface
 346 !  ic  * msg      number of missing moisture levels at the top of model.
 347 !  w  * p        grid slice of ambient mid-layer pressure in mbs.
 348 !  i  * pblt     row of pbl top indices.
 349 !  w  * pcpdh    scaled surface pressure.
 350 !  w  * pf       grid slice of ambient interface pressure in mbs.
 351 !  wg * pg       grid slice of gathered values of p.
 352 !  w  * q        grid slice of mixing ratio.
 353 !  wg * qd       grid slice of mixing ratio in downdraft.
 354 !  wg * qg       grid slice of gathered values of q.
 355 !  i/o * qh       grid slice of specific humidity.
 356 !  w  * qh0      grid slice of initial specific humidity.
 357 !  wg * qhat     grid slice of upper interface mixing ratio.
 358 !  wg * ql       grid slice of cloud liquid water.
 359 !  wg * qs       grid slice of saturation mixing ratio.
 360 !  w  * qstp     grid slice of parcel temp. saturation mixing ratio.
 361 !  wg * qstpg    grid slice of gathered values of qstp.
 362 !  wg * qu       grid slice of mixing ratio in updraft.
 363 !  ic  * rgas     dry air gas constant.
 364 !  wg * rl       latent heat of vaporization.
 365 !  w  * s        grid slice of scaled dry static energy (t+gz/cp).
 366 !  wg * sd       grid slice of dry static energy in downdraft.
 367 !  wg * sg       grid slice of gathered values of s.
 368 !  wg * shat     grid slice of upper interface dry static energy.
 369 !  wg * su       grid slice of dry static energy in updraft.
 370 !  i/o * t
 371 !  o  * jctop    row of top-of-deep-convection indices passed out.
 372 !  O  * jcbot    row of base of cloud indices passed out.
 373 !  wg * tg       grid slice of gathered values of t.
 374 !  w  * tl       row of parcel temperature at lcl.
 375 !  wg * tlg      grid slice of gathered values of tl.
 376 !  w  * tp       grid slice of parcel temperatures.
 377 !  wg * tpg      grid slice of gathered values of tp.
 378 !  i/o * u        grid slice of u-wind (real).
 379 !  wg * ug       grid slice of gathered values of u.
 380 !  i/o * utg      grid slice of u-wind tendency (real).
 381 !  i/o * v        grid slice of v-wind (real).
 382 !  w  * va       work array re-used by called subroutines.
 383 !  wg * vg       grid slice of gathered values of v.
 384 !  i/o * vtg      grid slice of v-wind tendency (real).
 385 !  i  * w        grid slice of diagnosed large-scale vertical velocity.
 386 !  w  * z        grid slice of ambient mid-layer height in metres.
 387 !  w  * zf       grid slice of ambient interface height in metres.
 388 !  wg * zfg      grid slice of gathered values of zf.
 389 !  wg * zg       grid slice of gathered values of z.
 390 !
 391 !-----------------------------------------------------------------------
 392 !
 393 ! multi-level i/o fields:
 394 !  i      => input arrays.
 395 !  i/o    => input/output arrays.
 396 !  w      => work arrays.
 397 !  wg     => work arrays operating only on gathered points.
 398 !  ic     => input data constants.
 399 !  c      => data constants pertaining to subroutine itself.
 400 !
 401 ! input arguments
 402 !
 403    integer, intent(in) :: lchnk                   ! chunk identifier
 404    integer, intent(in) :: ncol                    ! number of atmospheric columns
 405
 406    real(r8), intent(in) :: t(pcols,pver)          ! grid slice of temperature at mid-layer.
 407    real(r8), intent(in) :: qh(pcols,pver,pcnst)   ! grid slice of specific humidity.
 408    real(r8), intent(in) :: pap(pcols,pver)
 409    real(r8), intent(in) :: paph(pcols,pver+1)
 410    real(r8), intent(in) :: dpp(pcols,pver)        ! local sigma half-level thickness (i.e. dshj).
 411    real(r8), intent(in) :: zm(pcols,pver)
 412    real(r8), intent(in) :: geos(pcols)
 413    real(r8), intent(in) :: zi(pcols,pver+1)
 414    real(r8), intent(in) :: pblh(pcols)
 415    real(r8), intent(in) :: tpert(pcols)
 416    real(r8), intent(in) :: landfrac(pcols) ! RBN Landfrac
 417 !
 418 ! output arguments
 419 !
 420    real(r8), intent(out) :: qtnd(pcols,pver)           ! specific humidity tendency (kg/kg/s)
 421    real(r8), intent(out) :: heat(pcols,pver)           ! heating rate (dry static energy tendency, W/kg)
 422    real(r8), intent(out) :: mcon(pcols,pverp)
 423    real(r8), intent(out) :: dlf(pcols,pver)    ! scattrd version of the detraining cld h2o tend
 424    real(r8), intent(out) :: pflx(pcols,pverp)  ! scattered precip flux at each level
 425    real(r8), intent(out) :: cme(pcols,pver)
 426    real(r8), intent(out) :: cape(pcols)        ! w  convective available potential energy.
 427    real(r8), intent(out) :: zdu(pcols,pver)
 428    real(r8), intent(out) :: rprd(pcols,pver)     ! rain production rate
 429 ! move these vars from local storage to output so that convective
 430 ! transports can be done in outside of conv_cam.
 431    real(r8), intent(out) :: mu(pcols,pver)
 432    real(r8), intent(out) :: eu(pcols,pver)
 433    real(r8), intent(out) :: du(pcols,pver)
 434    real(r8), intent(out) :: md(pcols,pver)
 435    real(r8), intent(out) :: ed(pcols,pver)
 436    real(r8), intent(out) :: dp(pcols,pver)       ! wg layer thickness in mbs (between upper/lower interface).
 437    real(r8), intent(out) :: dsubcld(pcols)       ! wg layer thickness in mbs between lcl and maxi.
 438    real(r8), intent(out) :: jctop(pcols)  ! o row of top-of-deep-convection indices passed out.
 439    real(r8), intent(out) :: jcbot(pcols)  ! o row of base of cloud indices passed out.
 440    real(r8), intent(out) :: prec(pcols)
 441    real(r8), intent(out) :: rliq(pcols) ! reserved liquid (not yet in cldliq) for energy integrals
 442
 443    real(r8) zs(pcols)
 444    real(r8) dlg(pcols,pver)    ! gathrd version of the detraining cld h2o tend
 445    real(r8) pflxg(pcols,pverp) ! gather precip flux at each level
 446    real(r8) cug(pcols,pver)    ! gathered condensation rate
 447    real(r8) evpg(pcols,pver)   ! gathered evap rate of rain in downdraft
 448    real(r8) mumax(pcols)
 449    integer jt(pcols)                          ! wg top  level index of deep cumulus convection.
 450    integer maxg(pcols)                        ! wg gathered values of maxi.
 451    integer ideep(pcols)                       ! w holds position of gathered points vs longitude index.
 452    integer lengath
 453 !     diagnostic field used by chem/wetdep codes
 454    real(r8) ql(pcols,pver)                    ! wg grid slice of cloud liquid water.
 455 !
 456    real(r8) pblt(pcols)           ! i row of pbl top indices.
 457
 458
 459
 460
 461 !
 462 !-----------------------------------------------------------------------
 463 !
 464 ! general work fields (local variables):
 465 !
 466    real(r8) q(pcols,pver)              ! w  grid slice of mixing ratio.
 467    real(r8) p(pcols,pver)              ! w  grid slice of ambient mid-layer pressure in mbs.
 468    real(r8) z(pcols,pver)              ! w  grid slice of ambient mid-layer height in metres.
 469    real(r8) s(pcols,pver)              ! w  grid slice of scaled dry static energy (t+gz/cp).
 470    real(r8) tp(pcols,pver)             ! w  grid slice of parcel temperatures.
 471    real(r8) zf(pcols,pver+1)           ! w  grid slice of ambient interface height in metres.
 472    real(r8) pf(pcols,pver+1)           ! w  grid slice of ambient interface pressure in mbs.
 473    real(r8) qstp(pcols,pver)           ! w  grid slice of parcel temp. saturation mixing ratio.
 474
 475    real(r8) tl(pcols)                  ! w  row of parcel temperature at lcl.
 476
 477    integer lcl(pcols)                  ! w  base level index of deep cumulus convection.
 478    integer lel(pcols)                  ! w  index of highest theoretical convective plume.
 479    integer lon(pcols)                  ! w  index of onset level for deep convection.
 480    integer maxi(pcols)                 ! w  index of level with largest moist static energy.
 481    integer index(pcols)
 482    real(r8) precip
 483 !
 484 ! gathered work fields:
 485 !
 486    real(r8) qg(pcols,pver)             ! wg grid slice of gathered values of q.
 487    real(r8) tg(pcols,pver)             ! w  grid slice of temperature at interface.
 488    real(r8) pg(pcols,pver)             ! wg grid slice of gathered values of p.
 489    real(r8) zg(pcols,pver)             ! wg grid slice of gathered values of z.
 490    real(r8) sg(pcols,pver)             ! wg grid slice of gathered values of s.
 491    real(r8) tpg(pcols,pver)            ! wg grid slice of gathered values of tp.
 492    real(r8) zfg(pcols,pver+1)          ! wg grid slice of gathered values of zf.
 493    real(r8) qstpg(pcols,pver)          ! wg grid slice of gathered values of qstp.
 494    real(r8) ug(pcols,pver)             ! wg grid slice of gathered values of u.
 495    real(r8) vg(pcols,pver)             ! wg grid slice of gathered values of v.
 496    real(r8) cmeg(pcols,pver)
 497
 498    real(r8) rprdg(pcols,pver)           ! wg gathered rain production rate
 499    real(r8) capeg(pcols)               ! wg gathered convective available potential energy.
 500    real(r8) tlg(pcols)                 ! wg grid slice of gathered values of tl.
 501    real(r8) landfracg(pcols)            ! wg grid slice of landfrac
 502
 503    integer lclg(pcols)       ! wg gathered values of lcl.
 504    integer lelg(pcols)
 505 !
 506 ! work fields arising from gathered calculations.
 507 !
 508    real(r8) dqdt(pcols,pver)           ! wg mixing ratio tendency at gathered points.
 509    real(r8) dsdt(pcols,pver)           ! wg dry static energy ("temp") tendency at gathered points.
 510 !      real(r8) alpha(pcols,pver)      ! array of vertical differencing used (=1. for upstream).
 511    real(r8) sd(pcols,pver)             ! wg grid slice of dry static energy in downdraft.
 512    real(r8) qd(pcols,pver)             ! wg grid slice of mixing ratio in downdraft.
 513    real(r8) mc(pcols,pver)             ! wg net upward (scaled by mb) cloud mass flux.
 514    real(r8) qhat(pcols,pver)           ! wg grid slice of upper interface mixing ratio.
 515    real(r8) qu(pcols,pver)             ! wg grid slice of mixing ratio in updraft.
 516    real(r8) su(pcols,pver)             ! wg grid slice of dry static energy in updraft.
 517    real(r8) qs(pcols,pver)             ! wg grid slice of saturation mixing ratio.
 518    real(r8) shat(pcols,pver)           ! wg grid slice of upper interface dry static energy.
 519    real(r8) hmn(pcols,pver)            ! wg moist static energy.
 520    real(r8) hsat(pcols,pver)           ! wg saturated moist static energy.
 521    real(r8) qlg(pcols,pver)
 522    real(r8) dudt(pcols,pver)           ! wg u-wind tendency at gathered points.
 523    real(r8) dvdt(pcols,pver)           ! wg v-wind tendency at gathered points.
 524 !      real(r8) ud(pcols,pver)
 525 !      real(r8) vd(pcols,pver)
 526
 527    real(r8) mb(pcols)                  ! wg cloud base mass flux.
 528
 529    integer jlcl(pcols)
 530    integer j0(pcols)                 ! wg detrainment initiation level index.
 531    integer jd(pcols)                 ! wg downdraft initiation level index.
 532
 533    real(r8) delt                     ! length of model time-step in seconds.
 534
 535    integer i
 536    integer ii
 537    integer k
 538    integer msg                      !  ic number of missing moisture levels at the top of model.
 539    real(r8) qdifr
 540    real(r8) sdifr
 541
 542 !
 543 !--------------------------Data statements------------------------------
 544 !
 545 ! Set internal variable "msg" (convection limit) to "limcnv-1"
 546 !
 547    msg = limcnv - 1
 548 !
 549 ! initialize necessary arrays.
 550 ! zero out variables not used in cam
 551 !
 552    qtnd(:,:) = 0._r8
 553    heat(:,:) = 0._r8
 554    mcon(:,:) = 0._r8
 555    rliq(:ncol)   = 0._r8
 556 !
 557 ! initialize convective tendencies
 558 !
 559    prec(:ncol) = 0._r8
 560    do k = 1,pver
 561       do i = 1,ncol
 562          dqdt(i,k)  = 0._r8
 563          dsdt(i,k)  = 0._r8
 564          dudt(i,k)  = 0._r8
 565          dvdt(i,k)  = 0._r8
 566          pflx(i,k)  = 0._r8
 567          pflxg(i,k) = 0._r8
 568          cme(i,k)   = 0._r8
 569          rprd(i,k)  = 0._r8
 570          zdu(i,k)   = 0._r8
 571          ql(i,k)    = 0._r8
 572          qlg(i,k)   = 0._r8
 573          dlf(i,k)   = 0._r8
 574          dlg(i,k)   = 0._r8
 575       end do
 576    end do
 577    do i = 1,ncol
 578       pflx(i,pverp) = 0
 579       pflxg(i,pverp) = 0
 580    end do
 581 !
 582    do i = 1,ncol
 583       pblt(i) = pver
 584       dsubcld(i) = 0._r8
 585
 586       jctop(i) = pver
 587       jcbot(i) = 1
 588
 589    end do
 590 !
 591 ! calculate local pressure (mbs) and height (m) for both interface
 592 ! and mid-layer locations.
 593 !
 594    do i = 1,ncol
 595       zs(i) = geos(i)*rgrav
 596       pf(i,pver+1) = paph(i,pver+1)*0.01_r8
 597       zf(i,pver+1) = zi(i,pver+1) + zs(i)
 598    end do
 599    do k = 1,pver
 600       do i = 1,ncol
 601          p(i,k) = pap(i,k)*0.01_r8
 602          pf(i,k) = paph(i,k)*0.01_r8
 603          z(i,k) = zm(i,k) + zs(i)
 604          zf(i,k) = zi(i,k) + zs(i)
 605       end do
 606    end do
 607 !
 608    do k = pver - 1,msg + 1,-1
 609       do i = 1,ncol
 610          if (abs(z(i,k)-zs(i)-pblh(i)) < (zf(i,k)-zf(i,k+1))*0.5_r8) pblt(i) = k
 611       end do
 612    end do
 613 !
 614 ! store incoming specific humidity field for subsequent calculation
 615 ! of precipitation (through change in storage).
 616 ! define dry static energy (normalized by cp).
 617 !
 618    do k = 1,pver
 619       do i = 1,ncol
 620          q(i,k) = qh(i,k,1)
 621          s(i,k) = t(i,k) + (grav/cpres)*z(i,k)
 622          tp(i,k)=0.0_r8
 623          shat(i,k) = s(i,k)
 624          qhat(i,k) = q(i,k)
 625       end do
 626    end do
 627
 628    do i = 1,ncol
 629       capeg(i) = 0._r8
 630       lclg(i) = 1
 631       lelg(i) = pver
 632       maxg(i) = 1
 633       tlg(i) = 400._r8
 634       dsubcld(i) = 0._r8
 635    end do
 636 #ifndef WRF_PORT
 637    if( cam_physpkg_is('cam3')) then
 638
 639       !  For cam3 physics package, call non-dilute
 640
 641       call buoyan(lchnk   ,ncol    , &
 642                   q       ,t       ,p       ,z       ,pf       , &
 643                   tp      ,qstp    ,tl      ,rl      ,cape     , &
 644                   pblt    ,lcl     ,lel     ,lon     ,maxi     , &
 645                   rgas    ,grav    ,cpres   ,msg     , &
 646                   tpert   )
 647    else
 648 #else
 649
 650       !  Evaluate Tparcel, qsat(Tparcel), buoyancy and CAPE,
 651       !     lcl, lel, parcel launch level at index maxi()=hmax
 652
 653       call buoyan_dilute(lchnk   ,ncol    , &
 654                   q       ,t       ,p       ,z       ,pf       , &
 655                   tp      ,qstp    ,tl      ,rl      ,cape     , &
 656                   pblt    ,lcl     ,lel     ,lon     ,maxi     , &
 657                   rgas    ,grav    ,cpres   ,msg     , &
 658                   tpert   )
 659 #endif
 660 #ifndef WRF_PORT
 661    end if
 662 #endif
 663
 664 !
 665 ! determine whether grid points will undergo some deep convection
 666 ! (ideep=1) or not (ideep=0), based on values of cape,lcl,lel
 667 ! (require cape.gt. 0 and lel<lcl as minimum conditions).
 668 !
 669    lengath = 0
 670    do i=1,ncol
 671       if (cape(i) > capelmt) then
 672          lengath = lengath + 1
 673          index(lengath) = i
 674       end if
 675    end do
 676
 677    if (lengath.eq.0) return
 678    do ii=1,lengath
 679       i=index(ii)
 680       ideep(ii)=i
 681    end do
 682 !
 683 ! obtain gathered arrays necessary for ensuing calculations.
 684 !
 685    do k = 1,pver
 686       do i = 1,lengath
 687          dp(i,k) = 0.01_r8*dpp(ideep(i),k)
 688          qg(i,k) = q(ideep(i),k)
 689          tg(i,k) = t(ideep(i),k)
 690          pg(i,k) = p(ideep(i),k)
 691          zg(i,k) = z(ideep(i),k)
 692          sg(i,k) = s(ideep(i),k)
 693          tpg(i,k) = tp(ideep(i),k)
 694          zfg(i,k) = zf(ideep(i),k)
 695          qstpg(i,k) = qstp(ideep(i),k)
 696          ug(i,k) = 0._r8
 697          vg(i,k) = 0._r8
 698       end do
 699    end do
 700 !
 701    do i = 1,lengath
 702       zfg(i,pver+1) = zf(ideep(i),pver+1)
 703    end do
 704    do i = 1,lengath
 705       capeg(i) = cape(ideep(i))
 706       lclg(i) = lcl(ideep(i))
 707       lelg(i) = lel(ideep(i))
 708       maxg(i) = maxi(ideep(i))
 709       tlg(i) = tl(ideep(i))
 710       landfracg(i) = landfrac(ideep(i))
 711    end do
 712 !
 713 ! calculate sub-cloud layer pressure "thickness" for use in
 714 ! closure and tendency routines.
 715 !
 716    do k = msg + 1,pver
 717       do i = 1,lengath
 718          if (k >= maxg(i)) then
 719             dsubcld(i) = dsubcld(i) + dp(i,k)
 720          end if
 721       end do
 722    end do
 723 !
 724 ! define array of factors (alpha) which defines interfacial
 725 ! values, as well as interfacial values for (q,s) used in
 726 ! subsequent routines.
 727 !
 728    do k = msg + 2,pver
 729       do i = 1,lengath
 730 !            alpha(i,k) = 0.5
 731          sdifr = 0._r8
 732          qdifr = 0._r8
 733          if (sg(i,k) > 0._r8 .or. sg(i,k-1) > 0._r8) &
 734             sdifr = abs((sg(i,k)-sg(i,k-1))/max(sg(i,k-1),sg(i,k)))
 735          if (qg(i,k) > 0._r8 .or. qg(i,k-1) > 0._r8) &
 736             qdifr = abs((qg(i,k)-qg(i,k-1))/max(qg(i,k-1),qg(i,k)))
 737          if (sdifr > 1.E-6_r8) then
 738             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))
 739          else
 740             shat(i,k) = 0.5_r8* (sg(i,k)+sg(i,k-1))
 741          end if
 742          if (qdifr > 1.E-6_r8) then
 743             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))
 744          else
 745             qhat(i,k) = 0.5_r8* (qg(i,k)+qg(i,k-1))
 746          end if
 747       end do
 748    end do
 749 !
 750 ! obtain cloud properties.
 751 !
 752
 753    call cldprp(lchnk   , &
 754                qg      ,tg      ,ug      ,vg      ,pg      , &
 755                zg      ,sg      ,mu      ,eu      ,du      , &
 756                md      ,ed      ,sd      ,qd      ,mc      , &
 757                qu      ,su      ,zfg     ,qs      ,hmn     , &
 758                hsat    ,shat    ,qlg     , &
 759                cmeg    ,maxg    ,lelg    ,jt      ,jlcl    , &
 760                maxg    ,j0      ,jd      ,rl      ,lengath , &
 761                rgas    ,grav    ,cpres   ,msg     , &
 762                pflxg   ,evpg    ,cug     ,rprdg   ,limcnv  ,landfracg)
 763 !
 764 ! convert detrainment from units of "1/m" to "1/mb".
 765 !
 766    do k = msg + 1,pver
 767       do i = 1,lengath
 768          du   (i,k) = du   (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
 769          eu   (i,k) = eu   (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
 770          ed   (i,k) = ed   (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
 771          cug  (i,k) = cug  (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
 772          cmeg (i,k) = cmeg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
 773          rprdg(i,k) = rprdg(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
 774          evpg (i,k) = evpg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
 775       end do
 776    end do
 777
 778    call closure(lchnk   , &
 779                 qg      ,tg      ,pg      ,zg      ,sg      , &
 780                 tpg     ,qs      ,qu      ,su      ,mc      , &
 781                 du      ,mu      ,md      ,qd      ,sd      , &
 782                 qhat    ,shat    ,dp      ,qstpg   ,zfg     , &
 783                 qlg     ,dsubcld ,mb      ,capeg   ,tlg     , &
 784                 lclg    ,lelg    ,jt      ,maxg    ,1       , &
 785                 lengath ,rgas    ,grav    ,cpres   ,rl      , &
 786                 msg     ,capelmt    )
 787 !
 788 ! limit cloud base mass flux to theoretical upper bound.
 789 !
 790    do i=1,lengath
 791       mumax(i) = 0
 792    end do
 793    do k=msg + 2,pver
 794       do i=1,lengath
 795         mumax(i) = max(mumax(i), mu(i,k)/dp(i,k))
 796       end do
 797    end do
 798
 799    do i=1,lengath
 800       if (mumax(i) > 0._r8) then
 801          mb(i) = min(mb(i),0.5_r8/(delt*mumax(i)))
 802       else
 803          mb(i) = 0._r8
 804       endif
 805    end do
 806    ! If no_deep_pbl = .true., don't allow convection entirely
 807    ! within PBL (suggestion of Bjorn Stevens, 8-2000)
 808
 809    if (no_deep_pbl) then
 810       do i=1,lengath
 811          if (zm(ideep(i),jt(i)) < pblh(ideep(i))) mb(i) = 0
 812       end do
 813    end if
 814
 815
 816    do k=msg+1,pver
 817       do i=1,lengath
 818          mu   (i,k)  = mu   (i,k)*mb(i)
 819          md   (i,k)  = md   (i,k)*mb(i)
 820          mc   (i,k)  = mc   (i,k)*mb(i)
 821          du   (i,k)  = du   (i,k)*mb(i)
 822          eu   (i,k)  = eu   (i,k)*mb(i)
 823          ed   (i,k)  = ed   (i,k)*mb(i)
 824          cmeg (i,k)  = cmeg (i,k)*mb(i)
 825          rprdg(i,k)  = rprdg(i,k)*mb(i)
 826          cug  (i,k)  = cug  (i,k)*mb(i)
 827          evpg (i,k)  = evpg (i,k)*mb(i)
 828          pflxg(i,k+1)= pflxg(i,k+1)*mb(i)*100._r8/grav
 829       end do
 830    end do
 831 !
 832 ! compute temperature and moisture changes due to convection.
 833 !
 834    call q1q2_pjr(lchnk   , &
 835                  dqdt    ,dsdt    ,qg      ,qs      ,qu      , &
 836                  su      ,du      ,qhat    ,shat    ,dp      , &
 837                  mu      ,md      ,sd      ,qd      ,qlg     , &
 838                  dsubcld ,jt      ,maxg    ,1       ,lengath , &
 839                  cpres   ,rl      ,msg     ,          &
 840                  dlg     ,evpg    ,cug     )
 841 !
 842 ! gather back temperature and mixing ratio.
 843 !
 844    do k = msg + 1,pver
 845 !DIR$ CONCURRENT
 846       do i = 1,lengath
 847 !
 848 ! q is updated to compute net precip.
 849 !
 850          q(ideep(i),k) = qh(ideep(i),k,1) + 2._r8*delt*dqdt(i,k)
 851          qtnd(ideep(i),k) = dqdt (i,k)
 852          cme (ideep(i),k) = cmeg (i,k)
 853          rprd(ideep(i),k) = rprdg(i,k)
 854          zdu (ideep(i),k) = du   (i,k)
 855          mcon(ideep(i),k) = mc   (i,k)
 856          heat(ideep(i),k) = dsdt (i,k)*cpres
 857          dlf (ideep(i),k) = dlg  (i,k)
 858          pflx(ideep(i),k) = pflxg(i,k)
 859          ql  (ideep(i),k) = qlg  (i,k)
 860       end do
 861    end do
 862 !
 863 !DIR$ CONCURRENT
 864    do i = 1,lengath
 865       jctop(ideep(i)) = jt(i)
 866 !++bee
 867       jcbot(ideep(i)) = maxg(i)
 868 !--bee
 869       pflx(ideep(i),pverp) = pflxg(i,pverp)
 870    end do
 871
 872 ! Compute precip by integrating change in water vapor minus detrained cloud water
 873    do k = pver,msg + 1,-1
 874       do i = 1,ncol
 875          prec(i) = prec(i) - dpp(i,k)* (q(i,k)-qh(i,k,1)) - dpp(i,k)*dlf(i,k)*2._r8*delt
 876       end do
 877    end do
 878
 879 ! obtain final precipitation rate in m/s.
 880    do i = 1,ncol
 881       prec(i) = rgrav*max(prec(i),0._r8)/ (2._r8*delt)/1000._r8
 882    end do
 883
 884 ! Compute reserved liquid (not yet in cldliq) for energy integrals.
 885 ! Treat rliq as flux out bottom, to be added back later.
 886    do k = 1, pver
 887       do i = 1, ncol
 888          rliq(i) = rliq(i) + dlf(i,k)*dpp(i,k)/gravit
 889       end do
 890    end do
 891    rliq(:ncol) = rliq(:ncol) /1000._r8
 892
 893    return
 894 end subroutine zm_convr
 895
 896 !===============================================================================
 897 subroutine zm_conv_evap(ncol,lchnk, &
 898      t,pmid,pdel,q, &
 899      tend_s, tend_s_snwprd, tend_s_snwevmlt, tend_q, &
 900      prdprec, cldfrc, deltat,  &
 901      prec, snow, ntprprd, ntsnprd, flxprec, flxsnow )
 902
 903 !-----------------------------------------------------------------------
 904 ! Compute tendencies due to evaporation of rain from ZM scheme
 905 !--
 906 ! Compute the total precipitation and snow fluxes at the surface.
 907 ! Add in the latent heat of fusion for snow formation and melt, since it not dealt with
 908 ! in the Zhang-MacFarlane parameterization.
 909 ! Evaporate some of the precip directly into the environment using a Sundqvist type algorithm
 910 !-----------------------------------------------------------------------
 911
 912     use wv_saturation,  only: aqsat
 913 #ifndef WRF_PORT
 914     use phys_grid, only: get_rlat_all_p
 915 #endif
 916
 917 !------------------------------Arguments--------------------------------
 918     integer,intent(in) :: ncol, lchnk             ! number of columns and chunk index
 919     real(r8),intent(in), dimension(pcols,pver) :: t          ! temperature (K)
 920     real(r8),intent(in), dimension(pcols,pver) :: pmid       ! midpoint pressure (Pa)
 921     real(r8),intent(in), dimension(pcols,pver) :: pdel       ! layer thickness (Pa)
 922     real(r8),intent(in), dimension(pcols,pver) :: q          ! water vapor (kg/kg)
 923     real(r8),intent(inout), dimension(pcols,pver) :: tend_s     ! heating rate (J/kg/s)
 924     real(r8),intent(inout), dimension(pcols,pver) :: tend_q     ! water vapor tendency (kg/kg/s)
 925     real(r8),intent(out  ), dimension(pcols,pver) :: tend_s_snwprd ! Heating rate of snow production
 926     real(r8),intent(out  ), dimension(pcols,pver) :: tend_s_snwevmlt ! Heating rate of evap/melting of snow
 927
 928
 929
 930     real(r8), intent(in   ) :: prdprec(pcols,pver)! precipitation production (kg/ks/s)
 931     real(r8), intent(in   ) :: cldfrc(pcols,pver) ! cloud fraction
 932     real(r8), intent(in   ) :: deltat             ! time step
 933
 934     real(r8), intent(inout) :: prec(pcols)        ! Convective-scale preciptn rate
 935     real(r8), intent(out)   :: snow(pcols)        ! Convective-scale snowfall rate
 936 !
 937 !---------------------------Local storage-------------------------------
 938
 939     real(r8) :: est    (pcols,pver)    ! Saturation vapor pressure
 940     real(r8) :: fice   (pcols,pver)    ! ice fraction in precip production
 941     real(r8) :: fsnow_conv(pcols,pver) ! snow fraction in precip production
 942     real(r8) :: qsat   (pcols,pver)    ! saturation specific humidity
 943     real(r8),intent(out) :: flxprec(pcols,pverp)   ! Convective-scale flux of precip at interfaces (kg/m2/s)
 944     real(r8),intent(out) :: flxsnow(pcols,pverp)   ! Convective-scale flux of snow   at interfaces (kg/m2/s)
 945     real(r8),intent(out) :: ntprprd(pcols,pver)    ! net precip production in layer
 946     real(r8),intent(out) :: ntsnprd(pcols,pver)    ! net snow production in layer
 947     real(r8) :: work1                  ! temp variable (pjr)
 948     real(r8) :: work2                  ! temp variable (pjr)
 949
 950     real(r8) :: evpvint(pcols)         ! vertical integral of evaporation
 951     real(r8) :: evpprec(pcols)         ! evaporation of precipitation (kg/kg/s)
 952     real(r8) :: evpsnow(pcols)         ! evaporation of snowfall (kg/kg/s)
 953     real(r8) :: snowmlt(pcols)         ! snow melt tendency in layer
 954     real(r8) :: flxsntm(pcols)         ! flux of snow into layer, after melting
 955
 956     real(r8) :: evplimit               ! temp variable for evaporation limits
 957     real(r8) :: rlat(pcols)
 958
 959     integer :: i,k                     ! longitude,level indices
 960
 961
 962 !-----------------------------------------------------------------------
 963
 964 ! convert input precip to kg/m2/s
 965     prec(:ncol) = prec(:ncol)*1000._r8
 966
 967 ! determine saturation vapor pressure
 968     call aqsat (t    ,pmid  ,est    ,qsat    ,pcols   , &
 969          ncol ,pver  ,1       ,pver    )
 970
 971 ! determine ice fraction in rain production (use cloud water parameterization fraction at present)
 972     call cldwat_fice(ncol, t, fice, fsnow_conv)
 973
 974 ! zero the flux integrals on the top boundary
 975     flxprec(:ncol,1) = 0._r8
 976     flxsnow(:ncol,1) = 0._r8
 977     evpvint(:ncol)   = 0._r8
 978
 979     do k = 1, pver
 980        do i = 1, ncol
 981
 982 ! Melt snow falling into layer, if necessary.
 983           if (t(i,k) > tmelt) then
 984              flxsntm(i) = 0._r8
 985              snowmlt(i) = flxsnow(i,k) * gravit/ pdel(i,k)
 986           else
 987              flxsntm(i) = flxsnow(i,k)
 988              snowmlt(i) = 0._r8
 989           end if
 990
 991 ! relative humidity depression must be > 0 for evaporation
 992           evplimit = max(1._r8 - q(i,k)/qsat(i,k), 0._r8)
 993
 994 ! total evaporation depends on flux in the top of the layer
 995 ! flux prec is the net production above layer minus evaporation into environmet
 996           evpprec(i) = ke * (1._r8 - cldfrc(i,k)) * evplimit * sqrt(flxprec(i,k))
 997 !**********************************************************
 998 !!          evpprec(i) = 0.    ! turn off evaporation for now
 999 !**********************************************************
1000
1001 ! Don't let evaporation supersaturate layer (approx). Layer may already be saturated.
1002 ! Currently does not include heating/cooling change to qsat
1003           evplimit   = max(0._r8, (qsat(i,k)-q(i,k)) / deltat)
1004
1005 ! Don't evaporate more than is falling into the layer - do not evaporate rain formed
1006 ! in this layer but if precip production is negative, remove from the available precip
1007 ! Negative precip production occurs because of evaporation in downdrafts.
1008 !!$          evplimit   = flxprec(i,k) * gravit / pdel(i,k) + min(prdprec(i,k), 0.)
1009           evplimit   = min(evplimit, flxprec(i,k) * gravit / pdel(i,k))
1010
1011 ! Total evaporation cannot exceed input precipitation
1012           evplimit   = min(evplimit, (prec(i) - evpvint(i)) * gravit / pdel(i,k))
1013
1014           evpprec(i) = min(evplimit, evpprec(i))
1015
1016 ! evaporation of snow depends on snow fraction of total precipitation in the top after melting
1017           if (flxprec(i,k) > 0._r8) then
1018 !            evpsnow(i) = evpprec(i) * flxsntm(i) / flxprec(i,k)
1019 !            prevent roundoff problems
1020              work1 = min(max(0._r8,flxsntm(i)/flxprec(i,k)),1._r8)
1021              evpsnow(i) = evpprec(i) * work1
1022           else
1023              evpsnow(i) = 0._r8
1024           end if
1025
1026 ! vertically integrated evaporation
1027           evpvint(i) = evpvint(i) + evpprec(i) * pdel(i,k)/gravit
1028
1029 ! net precip production is production - evaporation
1030           ntprprd(i,k) = prdprec(i,k) - evpprec(i)
1031 ! net snow production is precip production * ice fraction - evaporation - melting
1032 !pjrworks ntsnprd(i,k) = prdprec(i,k)*fice(i,k) - evpsnow(i) - snowmlt(i)
1033 !pjrwrks2 ntsnprd(i,k) = prdprec(i,k)*fsnow_conv(i,k) - evpsnow(i) - snowmlt(i)
1034 ! the small amount added to flxprec in the work1 expression has been increased from
1035 ! 1e-36 to 8.64e-11 (1e-5 mm/day).  This causes the temperature based partitioning
1036 ! scheme to be used for small flxprec amounts.  This is to address error growth problems.
1037 #ifdef PERGRO
1038           work1 = min(max(0._r8,flxsnow(i,k)/(flxprec(i,k)+8.64e-11_r8)),1._r8)
1039 #else
1040           if (flxprec(i,k).gt.0._r8) then
1041              work1 = min(max(0._r8,flxsnow(i,k)/flxprec(i,k)),1._r8)
1042           else
1043              work1 = 0._r8
1044           endif
1045 #endif
1046           work2 = max(fsnow_conv(i,k), work1)
1047           if (snowmlt(i).gt.0._r8) work2 = 0._r8
1048 !         work2 = fsnow_conv(i,k)
1049           ntsnprd(i,k) = prdprec(i,k)*work2 - evpsnow(i) - snowmlt(i)
1050           tend_s_snwprd  (i,k) = prdprec(i,k)*work2*latice
1051           tend_s_snwevmlt(i,k) = - ( evpsnow(i) + snowmlt(i) )*latice
1052
1053 ! precipitation fluxes
1054           flxprec(i,k+1) = flxprec(i,k) + ntprprd(i,k) * pdel(i,k)/gravit
1055           flxsnow(i,k+1) = flxsnow(i,k) + ntsnprd(i,k) * pdel(i,k)/gravit
1056
1057 ! protect against rounding error
1058           flxprec(i,k+1) = max(flxprec(i,k+1), 0._r8)
1059           flxsnow(i,k+1) = max(flxsnow(i,k+1), 0._r8)
1060 ! more protection (pjr)
1061 !         flxsnow(i,k+1) = min(flxsnow(i,k+1), flxprec(i,k+1))
1062
1063 ! heating (cooling) and moistening due to evaporation
1064 ! - latent heat of vaporization for precip production has already been accounted for
1065 ! - snow is contained in prec
1066           tend_s(i,k)   =-evpprec(i)*latvap + ntsnprd(i,k)*latice
1067           tend_q(i,k) = evpprec(i)
1068        end do
1069     end do
1070
1071 ! set output precipitation rates (m/s)
1072     prec(:ncol) = flxprec(:ncol,pver+1) / 1000._r8
1073     snow(:ncol) = flxsnow(:ncol,pver+1) / 1000._r8
1074
1075 !**********************************************************
1076 !!$    tend_s(:ncol,:)   = 0.      ! turn heating off
1077 !**********************************************************
1078
1079   end subroutine zm_conv_evap
1080
1081
1082
1083 subroutine convtran(lchnk   , &
1084                     doconvtran,q       ,ncnst   ,mu      ,md      , &
1085                     du      ,eu      ,ed      ,dp      ,dsubcld , &
1086                     jt      ,mx      ,ideep   ,il1g    ,il2g    , &
1087                     nstep   ,fracis  ,dqdt    ,dpdry   )
1088 !-----------------------------------------------------------------------
1089 !
1090 ! Purpose:
1091 ! Convective transport of trace species
1092 !
1093 ! Mixing ratios may be with respect to either dry or moist air
1094 !
1095 ! Method:
1096 ! <Describe the algorithm(s) used in the routine.>
1097 ! <Also include any applicable external references.>
1098 !
1099 ! Author: P. Rasch
1100 !
1101 !-----------------------------------------------------------------------
1102    use shr_kind_mod, only: r8 => shr_kind_r8
1103    use constituents,    only: cnst_get_type_byind
1104 #ifndef WRF_PORT
1105    use ppgrid
1106    use abortutils, only: endrun
1107 #endif
1108
1109    implicit none
1110 !-----------------------------------------------------------------------
1111 !
1112 ! Input arguments
1113 !
1114    integer, intent(in) :: lchnk                 ! chunk identifier
1115    integer, intent(in) :: ncnst                 ! number of tracers to transport
1116    logical, intent(in) :: doconvtran(ncnst)     ! flag for doing convective transport
1117    real(r8), intent(in) :: q(pcols,pver,ncnst)  ! Tracer array including moisture
1118    real(r8), intent(in) :: mu(pcols,pver)       ! Mass flux up
1119    real(r8), intent(in) :: md(pcols,pver)       ! Mass flux down
1120    real(r8), intent(in) :: du(pcols,pver)       ! Mass detraining from updraft
1121    real(r8), intent(in) :: eu(pcols,pver)       ! Mass entraining from updraft
1122    real(r8), intent(in) :: ed(pcols,pver)       ! Mass entraining from downdraft
1123    real(r8), intent(in) :: dp(pcols,pver)       ! Delta pressure between interfaces
1124    real(r8), intent(in) :: dsubcld(pcols)       ! Delta pressure from cloud base to sfc
1125    real(r8), intent(in) :: fracis(pcols,pver,ncnst) ! fraction of tracer that is insoluble
1126
1127    integer, intent(in) :: jt(pcols)         ! Index of cloud top for each column
1128    integer, intent(in) :: mx(pcols)         ! Index of cloud top for each column
1129    integer, intent(in) :: ideep(pcols)      ! Gathering array
1130    integer, intent(in) :: il1g              ! Gathered min lon indices over which to operate
1131    integer, intent(in) :: il2g              ! Gathered max lon indices over which to operate
1132    integer, intent(in) :: nstep             ! Time step index
1133
1134    real(r8), intent(in) :: dpdry(pcols,pver)       ! Delta pressure between interfaces
1135
1136
1137 ! input/output
1138
1139    real(r8), intent(out) :: dqdt(pcols,pver,ncnst)  ! Tracer tendency array
1140
1141 !--------------------------Local Variables------------------------------
1142
1143    integer i                 ! Work index
1144    integer k                 ! Work index
1145    integer kbm               ! Highest altitude index of cloud base
1146    integer kk                ! Work index
1147    integer kkp1              ! Work index
1148    integer km1               ! Work index
1149    integer kp1               ! Work index
1150    integer ktm               ! Highest altitude index of cloud top
1151    integer m                 ! Work index
1152
1153    real(r8) cabv                 ! Mix ratio of constituent above
1154    real(r8) cbel                 ! Mix ratio of constituent below
1155    real(r8) cdifr                ! Normalized diff between cabv and cbel
1156    real(r8) chat(pcols,pver)     ! Mix ratio in env at interfaces
1157    real(r8) cond(pcols,pver)     ! Mix ratio in downdraft at interfaces
1158    real(r8) const(pcols,pver)    ! Gathered tracer array
1159    real(r8) fisg(pcols,pver)     ! gathered insoluble fraction of tracer
1160    real(r8) conu(pcols,pver)     ! Mix ratio in updraft at interfaces
1161    real(r8) dcondt(pcols,pver)   ! Gathered tend array
1162    real(r8) small                ! A small number
1163    real(r8) mbsth                ! Threshold for mass fluxes
1164    real(r8) mupdudp              ! A work variable
1165    real(r8) minc                 ! A work variable
1166    real(r8) maxc                 ! A work variable
1167    real(r8) fluxin               ! A work variable
1168    real(r8) fluxout              ! A work variable
1169    real(r8) netflux              ! A work variable
1170
1171    real(r8) dutmp(pcols,pver)       ! Mass detraining from updraft
1172    real(r8) eutmp(pcols,pver)       ! Mass entraining from updraft
1173    real(r8) edtmp(pcols,pver)       ! Mass entraining from downdraft
1174    real(r8) dptmp(pcols,pver)    ! Delta pressure between interfaces
1175 !-----------------------------------------------------------------------
1176 !
1177    small = 1.e-36_r8
1178 ! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
1179    mbsth = 1.e-15_r8
1180
1181 ! Find the highest level top and bottom levels of convection
1182    ktm = pver
1183    kbm = pver
1184    do i = il1g, il2g
1185       ktm = min(ktm,jt(i))
1186       kbm = min(kbm,mx(i))
1187    end do
1188
1189 ! Loop ever each constituent
1190    do m = 2, ncnst
1191       if (doconvtran(m)) then
1192
1193          if (cnst_get_type_byind(m).eq.'dry') then
1194             do k = 1,pver
1195                do i =il1g,il2g
1196                   dptmp(i,k) = dpdry(i,k)
1197                   dutmp(i,k) = du(i,k)*dp(i,k)/dpdry(i,k)
1198                   eutmp(i,k) = eu(i,k)*dp(i,k)/dpdry(i,k)
1199                   edtmp(i,k) = ed(i,k)*dp(i,k)/dpdry(i,k)
1200                end do
1201             end do
1202          else
1203             do k = 1,pver
1204                do i =il1g,il2g
1205                   dptmp(i,k) = dp(i,k)
1206                   dutmp(i,k) = du(i,k)
1207                   eutmp(i,k) = eu(i,k)
1208                   edtmp(i,k) = ed(i,k)
1209                end do
1210             end do
1211          endif
1212 !        dptmp = dp
1213
1214 ! Gather up the constituent and set tend to zero
1215          do k = 1,pver
1216             do i =il1g,il2g
1217                const(i,k) = q(ideep(i),k,m)
1218                fisg(i,k) = fracis(ideep(i),k,m)
1219             end do
1220          end do
1221
1222 ! From now on work only with gathered data
1223
1224 ! Interpolate environment tracer values to interfaces
1225          do k = 1,pver
1226             km1 = max(1,k-1)
1227             do i = il1g, il2g
1228                minc = min(const(i,km1),const(i,k))
1229                maxc = max(const(i,km1),const(i,k))
1230                if (minc < 0) then
1231                   cdifr = 0._r8
1232                else
1233                   cdifr = abs(const(i,k)-const(i,km1))/max(maxc,small)
1234                endif
1235
1236 ! If the two layers differ significantly use a geometric averaging
1237 ! procedure
1238                if (cdifr > 1.E-6_r8) then
1239                   cabv = max(const(i,km1),maxc*1.e-12_r8)
1240                   cbel = max(const(i,k),maxc*1.e-12_r8)
1241                   chat(i,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel
1242
1243                else             ! Small diff, so just arithmetic mean
1244                   chat(i,k) = 0.5_r8* (const(i,k)+const(i,km1))
1245                end if
1246
1247 ! Provisional up and down draft values
1248                conu(i,k) = chat(i,k)
1249                cond(i,k) = chat(i,k)
1250
1251 !              provisional tends
1252                dcondt(i,k) = 0._r8
1253
1254             end do
1255          end do
1256
1257 ! Do levels adjacent to top and bottom
1258          k = 2
1259          km1 = 1
1260          kk = pver
1261          do i = il1g,il2g
1262             mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
1263             if (mupdudp > mbsth) then
1264                conu(i,kk) = (+eutmp(i,kk)*fisg(i,kk)*const(i,kk)*dptmp(i,kk))/mupdudp
1265             endif
1266             if (md(i,k) < -mbsth) then
1267                cond(i,k) =  (-edtmp(i,km1)*fisg(i,km1)*const(i,km1)*dptmp(i,km1))/md(i,k)
1268             endif
1269          end do
1270
1271 ! Updraft from bottom to top
1272          do kk = pver-1,1,-1
1273             kkp1 = min(pver,kk+1)
1274             do i = il1g,il2g
1275                mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
1276                if (mupdudp > mbsth) then
1277                   conu(i,kk) = (  mu(i,kkp1)*conu(i,kkp1)+eutmp(i,kk)*fisg(i,kk)* &
1278                                   const(i,kk)*dptmp(i,kk) )/mupdudp
1279                endif
1280             end do
1281          end do
1282
1283 ! Downdraft from top to bottom
1284          do k = 3,pver
1285             km1 = max(1,k-1)
1286             do i = il1g,il2g
1287                if (md(i,k) < -mbsth) then
1288                   cond(i,k) =  (  md(i,km1)*cond(i,km1)-edtmp(i,km1)*fisg(i,km1)*const(i,km1) &
1289                                   *dptmp(i,km1) )/md(i,k)
1290                endif
1291             end do
1292          end do
1293
1294
1295          do k = ktm,pver
1296             km1 = max(1,k-1)
1297             kp1 = min(pver,k+1)
1298             do i = il1g,il2g
1299
1300 ! version 1 hard to check for roundoff errors
1301 !               dcondt(i,k) =
1302 !     $                  +(+mu(i,kp1)* (conu(i,kp1)-chat(i,kp1))
1303 !     $                    -mu(i,k)*   (conu(i,k)-chat(i,k))
1304 !     $                    +md(i,kp1)* (cond(i,kp1)-chat(i,kp1))
1305 !     $                    -md(i,k)*   (cond(i,k)-chat(i,k))
1306 !     $                   )/dp(i,k)
1307
1308 ! version 2 hard to limit fluxes
1309 !               fluxin =  mu(i,kp1)*conu(i,kp1) + mu(i,k)*chat(i,k)
1310 !     $                 -(md(i,k)  *cond(i,k)   + md(i,kp1)*chat(i,kp1))
1311 !               fluxout = mu(i,k)*conu(i,k)     + mu(i,kp1)*chat(i,kp1)
1312 !     $                 -(md(i,kp1)*cond(i,kp1) + md(i,k)*chat(i,k))
1313
1314 ! version 3 limit fluxes outside convection to mass in appropriate layer
1315 ! these limiters are probably only safe for positive definite quantitities
1316 ! it assumes that mu and md already satify a courant number limit of 1
1317                fluxin =  mu(i,kp1)*conu(i,kp1)+ mu(i,k)*min(chat(i,k),const(i,km1)) &
1318                          -(md(i,k)  *cond(i,k) + md(i,kp1)*min(chat(i,kp1),const(i,kp1)))
1319                fluxout = mu(i,k)*conu(i,k) + mu(i,kp1)*min(chat(i,kp1),const(i,k)) &
1320                          -(md(i,kp1)*cond(i,kp1) + md(i,k)*min(chat(i,k),const(i,k)))
1321
1322                netflux = fluxin - fluxout
1323                if (abs(netflux) < max(fluxin,fluxout)*1.e-12_r8) then
1324                   netflux = 0._r8
1325                endif
1326                dcondt(i,k) = netflux/dptmp(i,k)
1327             end do
1328          end do
1329 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1330 !
1331 !DIR$ NOINTERCHANGE
1332          do k = kbm,pver
1333             km1 = max(1,k-1)
1334             do i = il1g,il2g
1335                if (k == mx(i)) then
1336
1337 ! version 1
1338 !                  dcondt(i,k) = (1./dsubcld(i))*
1339 !     $              (-mu(i,k)*(conu(i,k)-chat(i,k))
1340 !     $               -md(i,k)*(cond(i,k)-chat(i,k))
1341 !     $              )
1342
1343 ! version 2
1344 !                  fluxin =  mu(i,k)*chat(i,k) - md(i,k)*cond(i,k)
1345 !                  fluxout = mu(i,k)*conu(i,k) - md(i,k)*chat(i,k)
1346 ! version 3
1347                   fluxin =  mu(i,k)*min(chat(i,k),const(i,km1)) - md(i,k)*cond(i,k)
1348                   fluxout = mu(i,k)*conu(i,k) - md(i,k)*min(chat(i,k),const(i,k))
1349
1350                   netflux = fluxin - fluxout
1351                   if (abs(netflux) < max(fluxin,fluxout)*1.e-12_r8) then
1352                      netflux = 0._r8
1353                   endif
1354 !                  dcondt(i,k) = netflux/dsubcld(i)
1355                   dcondt(i,k) = netflux/dptmp(i,k)
1356                else if (k > mx(i)) then
1357 !                  dcondt(i,k) = dcondt(i,k-1)
1358                   dcondt(i,k) = 0._r8
1359                end if
1360             end do
1361          end do
1362
1363 ! Initialize to zero everywhere, then scatter tendency back to full array
1364          dqdt(:,:,m) = 0._r8
1365          do k = 1,pver
1366             kp1 = min(pver,k+1)
1367 !DIR$ CONCURRENT
1368             do i = il1g,il2g
1369                dqdt(ideep(i),k,m) = dcondt(i,k)
1370             end do
1371          end do
1372
1373       end if      ! for doconvtran
1374
1375    end do
1376
1377    return
1378 end subroutine convtran
1379
1380 !=========================================================================================
1381
1382 subroutine momtran(lchnk, ncol, &
1383                     domomtran,q       ,ncnst   ,mu      ,md    , &
1384                     du      ,eu      ,ed      ,dp      ,dsubcld , &
1385                     jt      ,mx      ,ideep   ,il1g    ,il2g    , &
1386                     nstep   ,dqdt    ,pguall     ,pgdall, icwu, icwd, dt, seten    )
1387 !-----------------------------------------------------------------------
1388 !
1389 ! Purpose:
1390 ! Convective transport of momentum
1391 !
1392 ! Mixing ratios may be with respect to either dry or moist air
1393 !
1394 ! Method:
1395 ! Based on the convtran subroutine by P. Rasch
1396 ! <Also include any applicable external references.>
1397 !
1398 ! Author: J. Richter and P. Rasch
1399 !
1400 !-----------------------------------------------------------------------
1401    use shr_kind_mod, only: r8 => shr_kind_r8
1402 #ifndef WRF_PORT
1403    use constituents,    only: cnst_get_type_byind
1404    use ppgrid
1405    use abortutils, only: endrun
1406 #endif
1407
1408    implicit none
1409 !-----------------------------------------------------------------------
1410 !
1411 ! Input arguments
1412 !
1413    integer, intent(in) :: lchnk                 ! chunk identifier
1414    integer, intent(in) :: ncol                  ! number of atmospheric columns
1415    integer, intent(in) :: ncnst                 ! number of tracers to transport
1416    logical, intent(in) :: domomtran(ncnst)      ! flag for doing convective transport
1417    real(r8), intent(in) :: q(pcols,pver,ncnst)  ! Wind array
1418    real(r8), intent(in) :: mu(pcols,pver)       ! Mass flux up
1419    real(r8), intent(in) :: md(pcols,pver)       ! Mass flux down
1420    real(r8), intent(in) :: du(pcols,pver)       ! Mass detraining from updraft
1421    real(r8), intent(in) :: eu(pcols,pver)       ! Mass entraining from updraft
1422    real(r8), intent(in) :: ed(pcols,pver)       ! Mass entraining from downdraft
1423    real(r8), intent(in) :: dp(pcols,pver)       ! Delta pressure between interfaces
1424    real(r8), intent(in) :: dsubcld(pcols)       ! Delta pressure from cloud base to sfc
1425    real(r8), intent(in)    :: dt                       !  time step in seconds : 2*delta_t
1426
1427    integer, intent(in) :: jt(pcols)         ! Index of cloud top for each column
1428    integer, intent(in) :: mx(pcols)         ! Index of cloud top for each column
1429    integer, intent(in) :: ideep(pcols)      ! Gathering array
1430    integer, intent(in) :: il1g              ! Gathered min lon indices over which to operate
1431    integer, intent(in) :: il2g              ! Gathered max lon indices over which to operate
1432    integer, intent(in) :: nstep             ! Time step index
1433
1434
1435
1436 ! input/output
1437
1438    real(r8), intent(out) :: dqdt(pcols,pver,ncnst)  ! Tracer tendency array
1439
1440 !--------------------------Local Variables------------------------------
1441
1442    integer i                 ! Work index
1443    integer k                 ! Work index
1444    integer kbm               ! Highest altitude index of cloud base
1445    integer kk                ! Work index
1446    integer kkp1              ! Work index
1447    integer kkm1              ! Work index
1448    integer km1               ! Work index
1449    integer kp1               ! Work index
1450    integer ktm               ! Highest altitude index of cloud top
1451    integer m                 ! Work index
1452    integer ii                 ! Work index
1453
1454    real(r8) cabv                 ! Mix ratio of constituent above
1455    real(r8) cbel                 ! Mix ratio of constituent below
1456    real(r8) cdifr                ! Normalized diff between cabv and cbel
1457    real(r8) chat(pcols,pver)     ! Mix ratio in env at interfaces
1458    real(r8) cond(pcols,pver)     ! Mix ratio in downdraft at interfaces
1459    real(r8) const(pcols,pver)    ! Gathered wind array
1460    real(r8) conu(pcols,pver)     ! Mix ratio in updraft at interfaces
1461    real(r8) dcondt(pcols,pver)   ! Gathered tend array
1462    real(r8) small                ! A small number
1463    real(r8) mbsth                ! Threshold for mass fluxes
1464    real(r8) mupdudp              ! A work variable
1465    real(r8) minc                 ! A work variable
1466    real(r8) maxc                 ! A work variable
1467    real(r8) fluxin               ! A work variable
1468    real(r8) fluxout              ! A work variable
1469    real(r8) netflux              ! A work variable
1470
1471    real(r8) momcu                ! constant for updraft pressure gradient term
1472    real(r8) momcd                ! constant for downdraft pressure gradient term
1473    real(r8) sum                  ! sum
1474    real(r8) sum2                  ! sum2
1475
1476    real(r8) mududp(pcols,pver) ! working variable
1477    real(r8) mddudp(pcols,pver)     ! working variable
1478
1479    real(r8) pgu(pcols,pver)      ! Pressure gradient term for updraft
1480    real(r8) pgd(pcols,pver)      ! Pressure gradient term for downdraft
1481
1482    real(r8),intent(out) ::  pguall(pcols,pver,ncnst)      ! Apparent force from  updraft PG
1483    real(r8),intent(out) ::  pgdall(pcols,pver,ncnst)      ! Apparent force from  downdraft PG
1484
1485    real(r8),intent(out) ::  icwu(pcols,pver,ncnst)      ! In-cloud winds in updraft
1486    real(r8),intent(out) ::  icwd(pcols,pver,ncnst)      ! In-cloud winds in downdraft
1487
1488    real(r8),intent(out) ::  seten(pcols,pver) ! Dry static energy tendency
1489    real(r8)                 gseten(pcols,pver) ! Gathered dry static energy tendency
1490
1491    real(r8)  mflux(pcols,pverp,ncnst)   ! Gathered momentum flux
1492
1493    real(r8)  wind0(pcols,pver,ncnst)       !  gathered  wind before time step
1494    real(r8)  windf(pcols,pver,ncnst)       !  gathered  wind after time step
1495    real(r8) fkeb, fket, ketend_cons, ketend, utop, ubot, vtop, vbot, gset2
1496
1497
1498 !-----------------------------------------------------------------------
1499 !
1500
1501 ! Initialize outgoing fields
1502    pguall(:,:,:)     = 0.0_r8
1503    pgdall(:,:,:)     = 0.0_r8
1504 ! Initialize in-cloud winds to environmental wind
1505    icwu(:ncol,:,:)       = q(:ncol,:,:)
1506    icwd(:ncol,:,:)       = q(:ncol,:,:)
1507
1508 ! Initialize momentum flux and  final winds
1509    mflux(:,:,:)       = 0.0_r8
1510    wind0(:,:,:)         = 0.0_r8
1511    windf(:,:,:)         = 0.0_r8
1512
1513 ! Initialize dry static energy
1514
1515    seten(:,:)         = 0.0_r8
1516    gseten(:,:)         = 0.0_r8
1517
1518 ! Define constants for parameterization
1519
1520    momcu = 0.4_r8
1521    momcd = 0.4_r8
1522
1523    small = 1.e-36_r8
1524 ! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
1525    mbsth = 1.e-15_r8
1526
1527 ! Find the highest level top and bottom levels of convection
1528    ktm = pver
1529    kbm = pver
1530    do i = il1g, il2g
1531       ktm = min(ktm,jt(i))
1532       kbm = min(kbm,mx(i))
1533    end do
1534
1535 ! Loop ever each wind component
1536    do m = 1, ncnst                    !start at m = 1 to transport momentum
1537       if (domomtran(m)) then
1538
1539 ! Gather up the winds and set tend to zero
1540          do k = 1,pver
1541             do i =il1g,il2g
1542                const(i,k) = q(ideep(i),k,m)
1543                 wind0(i,k,m) = const(i,k)
1544             end do
1545          end do
1546
1547
1548 ! From now on work only with gathered data
1549
1550 ! Interpolate winds to interfaces
1551
1552          do k = 1,pver
1553             km1 = max(1,k-1)
1554             do i = il1g, il2g
1555
1556                ! use arithmetic mean
1557                chat(i,k) = 0.5_r8* (const(i,k)+const(i,km1))
1558
1559 ! Provisional up and down draft values
1560                conu(i,k) = chat(i,k)
1561                cond(i,k) = chat(i,k)
1562
1563 !              provisional tends
1564                dcondt(i,k) = 0._r8
1565
1566             end do
1567          end do
1568
1569
1570 !
1571 ! Pressure Perturbation Term
1572 !
1573
1574       !Top boundary:  assume mu is zero
1575
1576          k=1
1577          pgu(:il2g,k) = 0.0_r8
1578          pgd(:il2g,k) = 0.0_r8
1579
1580          do k=2,pver-1
1581             km1 = max(1,k-1)
1582             kp1 = min(pver,k+1)
1583             do i = il1g,il2g
1584
1585                !interior points
1586
1587                mududp(i,k) =  ( mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
1588                            +  mu(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
1589
1590                pgu(i,k) = - momcu * 0.5_r8 * mududp(i,k)
1591
1592
1593                mddudp(i,k) =  ( md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
1594                            +  md(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
1595
1596                pgd(i,k) = - momcd * 0.5_r8 * mddudp(i,k)
1597
1598
1599             end do
1600          end do
1601
1602        ! bottom boundary
1603        k = pver
1604        km1 = max(1,k-1)
1605        do i=il1g,il2g
1606
1607           mududp(i,k) =   mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
1608           pgu(i,k) = - momcu *  mududp(i,k)
1609
1610           mddudp(i,k) =   md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
1611
1612           pgd(i,k) = - momcd * mddudp(i,k)
1613
1614        end do
1615
1616
1617 !
1618 ! In-cloud velocity calculations
1619 !
1620
1621 ! Do levels adjacent to top and bottom
1622          k = 2
1623          km1 = 1
1624          kk = pver
1625          kkm1 = max(1,kk-1)
1626          do i = il1g,il2g
1627             mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
1628             if (mupdudp > mbsth) then
1629
1630                conu(i,kk) = (+eu(i,kk)*const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp
1631             endif
1632             if (md(i,k) < -mbsth) then
1633                cond(i,k) =  (-ed(i,km1)*const(i,km1)*dp(i,km1))-pgd(i,km1)*dp(i,km1)/md(i,k)
1634             endif
1635
1636
1637          end do
1638
1639
1640
1641 ! Updraft from bottom to top
1642          do kk = pver-1,1,-1
1643             kkm1 = max(1,kk-1)
1644             kkp1 = min(pver,kk+1)
1645             do i = il1g,il2g
1646                mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
1647                if (mupdudp > mbsth) then
1648
1649                   conu(i,kk) = (  mu(i,kkp1)*conu(i,kkp1)+eu(i,kk)* &
1650                                   const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp
1651                endif
1652             end do
1653
1654          end do
1655
1656
1657 ! Downdraft from top to bottom
1658          do k = 3,pver
1659             km1 = max(1,k-1)
1660             do i = il1g,il2g
1661                if (md(i,k) < -mbsth) then
1662
1663                   cond(i,k) =  (  md(i,km1)*cond(i,km1)-ed(i,km1)*const(i,km1) &
1664                                   *dp(i,km1)-pgd(i,km1)*dp(i,km1) )/md(i,k)
1665
1666                endif
1667             end do
1668          end do
1669
1670
1671          sum = 0._r8
1672          sum2 = 0._r8
1673
1674
1675          do k = ktm,pver
1676             km1 = max(1,k-1)
1677             kp1 = min(pver,k+1)
1678             do i = il1g,il2g
1679                ii = ideep(i)
1680
1681 ! version 1 hard to check for roundoff errors
1682                dcondt(i,k) =  &
1683                            +(mu(i,kp1)* (conu(i,kp1)-chat(i,kp1)) &
1684                            -mu(i,k)*   (conu(i,k)-chat(i,k))      &
1685                            +md(i,kp1)* (cond(i,kp1)-chat(i,kp1)) &
1686                            -md(i,k)*   (cond(i,k)-chat(i,k)) &
1687                           )/dp(i,k)
1688
1689             end do
1690          end do
1691
1692   ! dcont for bottom layer
1693           !
1694           !DIR$ NOINTERCHANGE
1695           do k = kbm,pver
1696              km1 = max(1,k-1)
1697              do i = il1g,il2g
1698                 if (k == mx(i)) then
1699
1700                    ! version 1
1701                    dcondt(i,k) = (1./dp(i,k))*   &
1702                         (-mu(i,k)*(conu(i,k)-chat(i,k)) &
1703                         -md(i,k)*(cond(i,k)-chat(i,k)) &
1704                         )
1705                 end if
1706              end do
1707           end do
1708
1709 ! Initialize to zero everywhere, then scatter tendency back to full array
1710          dqdt(:,:,m) = 0._r8
1711
1712          do k = 1,pver
1713             do i = il1g,il2g
1714                ii = ideep(i)
1715                dqdt(ii,k,m) = dcondt(i,k)
1716     ! Output apparent force on the mean flow from pressure gradient
1717                pguall(ii,k,m) = -pgu(i,k)
1718                pgdall(ii,k,m) = -pgd(i,k)
1719                icwu(ii,k,m)   =  conu(i,k)
1720                icwd(ii,k,m)   =  cond(i,k)
1721             end do
1722          end do
1723
1724           ! Calculate momentum flux in units of mb*m/s2
1725
1726           do k = ktm,pver
1727              do i = il1g,il2g
1728                 ii = ideep(i)
1729                 mflux(i,k,m) = &
1730                      -mu(i,k)*   (conu(i,k)-chat(i,k))      &
1731                      -md(i,k)*   (cond(i,k)-chat(i,k))
1732              end do
1733           end do
1734
1735
1736           ! Calculate winds at the end of the time step
1737
1738           do k = ktm,pver
1739              do i = il1g,il2g
1740                 ii = ideep(i)
1741                 km1 = max(1,k-1)
1742                 kp1 = k+1
1743                 windf(i,k,m) = const(i,k)    -   (mflux(i,kp1,m) - mflux(i,k,m)) * dt /dp(i,k)
1744
1745              end do
1746           end do
1747
1748        end if      ! for domomtran
1749    end do
1750
1751  ! Need to add an energy fix to account for the dissipation of kinetic energy
1752     ! Formulation follows from Boville and Bretherton (2003)
1753     ! formulation by PJR
1754
1755     do k = ktm,pver
1756        km1 = max(1,k-1)
1757        kp1 = min(pver,k+1)
1758        do i = il1g,il2g
1759
1760           ii = ideep(i)
1761
1762           ! calculate the KE fluxes at top and bot of layer
1763           ! based on a discrete approximation to b&b eq(35) F_KE = u*F_u + v*F_v at interface
1764           utop = (wind0(i,k,1)+wind0(i,km1,1))/2.
1765           vtop = (wind0(i,k,2)+wind0(i,km1,2))/2.
1766           ubot = (wind0(i,kp1,1)+wind0(i,k,1))/2.
1767           vbot = (wind0(i,kp1,2)+wind0(i,k,2))/2.
1768           fket = utop*mflux(i,k,1)   + vtop*mflux(i,k,2)    ! top of layer
1769           fkeb = ubot*mflux(i,k+1,1) + vbot*mflux(i,k+1,2)  ! bot of layer
1770
1771           ! divergence of these fluxes should give a conservative redistribution of KE
1772           ketend_cons = (fket-fkeb)/dp(i,k)
1773
1774           ! tendency in kinetic energy resulting from the momentum transport
1775           ketend = ((windf(i,k,1)**2 + windf(i,k,2)**2) - (wind0(i,k,1)**2 + wind0(i,k,2)**2))*0.5/dt
1776
1777           ! the difference should be the dissipation
1778           gset2 = ketend_cons - ketend
1779           gseten(i,k) = gset2
1780
1781        end do
1782
1783     end do
1784
1785     ! Scatter dry static energy to full array
1786     do k = 1,pver
1787        do i = il1g,il2g
1788           ii = ideep(i)
1789           seten(ii,k) = gseten(i,k)
1790
1791        end do
1792     end do
1793
1794    return
1795 end subroutine momtran
1796
1797 !=========================================================================================
1798
1799 subroutine buoyan(lchnk   ,ncol    , &
1800                   q       ,t       ,p       ,z       ,pf      , &
1801                   tp      ,qstp    ,tl      ,rl      ,cape    , &
1802                   pblt    ,lcl     ,lel     ,lon     ,mx      , &
1803                   rd      ,grav    ,cp      ,msg     , &
1804                   tpert   )
1805 !-----------------------------------------------------------------------
1806 !
1807 ! Purpose:
1808 ! <Say what the routine does>
1809 !
1810 ! Method:
1811 ! <Describe the algorithm(s) used in the routine.>
1812 ! <Also include any applicable external references.>
1813 !
1814 ! Author:
1815 ! This is contributed code not fully standardized by the CCM core group.
1816 ! The documentation has been enhanced to the degree that we are able.
1817 ! Reviewed:          P. Rasch, April 1996
1818 !
1819 !-----------------------------------------------------------------------
1820    implicit none
1821 !-----------------------------------------------------------------------
1822 !
1823 ! input arguments
1824 !
1825    integer, intent(in) :: lchnk                 ! chunk identifier
1826    integer, intent(in) :: ncol                  ! number of atmospheric columns
1827
1828    real(r8), intent(in) :: q(pcols,pver)        ! spec. humidity
1829    real(r8), intent(in) :: t(pcols,pver)        ! temperature
1830    real(r8), intent(in) :: p(pcols,pver)        ! pressure
1831    real(r8), intent(in) :: z(pcols,pver)        ! height
1832    real(r8), intent(in) :: pf(pcols,pver+1)     ! pressure at interfaces
1833    real(r8), intent(in) :: pblt(pcols)          ! index of pbl depth
1834    real(r8), intent(in) :: tpert(pcols)         ! perturbation temperature by pbl processes
1835
1836 !
1837 ! output arguments
1838 !
1839    real(r8), intent(out) :: tp(pcols,pver)       ! parcel temperature
1840    real(r8), intent(out) :: qstp(pcols,pver)     ! saturation mixing ratio of parcel
1841    real(r8), intent(out) :: tl(pcols)            ! parcel temperature at lcl
1842    real(r8), intent(out) :: cape(pcols)          ! convective aval. pot. energy.
1843    integer lcl(pcols)        !
1844    integer lel(pcols)        !
1845    integer lon(pcols)        ! level of onset of deep convection
1846    integer mx(pcols)         ! level of max moist static energy
1847 !
1848 !--------------------------Local Variables------------------------------
1849 !
1850    real(r8) capeten(pcols,5)     ! provisional value of cape
1851    real(r8) tv(pcols,pver)       !
1852    real(r8) tpv(pcols,pver)      !
1853    real(r8) buoy(pcols,pver)
1854
1855    real(r8) a1(pcols)
1856    real(r8) a2(pcols)
1857    real(r8) estp(pcols)
1858    real(r8) pl(pcols)
1859    real(r8) plexp(pcols)
1860    real(r8) hmax(pcols)
1861    real(r8) hmn(pcols)
1862    real(r8) y(pcols)
1863
1864    logical plge600(pcols)
1865    integer knt(pcols)
1866    integer lelten(pcols,5)
1867
1868    real(r8) cp
1869    real(r8) e
1870    real(r8) grav
1871
1872    integer i
1873    integer k
1874    integer msg
1875    integer n
1876
1877    real(r8) rd
1878    real(r8) rl
1879 #ifdef PERGRO
1880    real(r8) rhd
1881 #endif
1882 !
1883 !-----------------------------------------------------------------------
1884 !
1885    do n = 1,5
1886       do i = 1,ncol
1887          lelten(i,n) = pver
1888          capeten(i,n) = 0._r8
1889       end do
1890    end do
1891 !
1892    do i = 1,ncol
1893       lon(i) = pver
1894       knt(i) = 0
1895       lel(i) = pver
1896       mx(i) = lon(i)
1897       cape(i) = 0._r8
1898       hmax(i) = 0._r8
1899    end do
1900
1901    tp(:ncol,:) = t(:ncol,:)
1902    qstp(:ncol,:) = q(:ncol,:)
1903
1904 !!! RBN - Initialize tv and buoy for output.
1905 !!! tv=tv : tpv=tpv : qstp=q : buoy=0.
1906    tv(:ncol,:) = t(:ncol,:) *(1._r8+1.608*q(:ncol,:))/ (1._r8+q(:ncol,:))
1907    tpv(:ncol,:) = tv(:ncol,:)
1908    buoy(:ncol,:) = 0._r8
1909
1910 !
1911 ! set "launching" level(mx) to be at maximum moist static energy.
1912 ! search for this level stops at planetary boundary layer top.
1913 !
1914 #ifdef PERGRO
1915    do k = pver,msg + 1,-1
1916       do i = 1,ncol
1917          hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
1918 !
1919 ! Reset max moist static energy level when relative difference exceeds 1.e-4
1920 !
1921          rhd = (hmn(i) - hmax(i))/(hmn(i) + hmax(i))
1922          if (k >= nint(pblt(i)) .and. k <= lon(i) .and. rhd > -1.e-4_r8) then
1923             hmax(i) = hmn(i)
1924             mx(i) = k
1925          end if
1926       end do
1927    end do
1928 #else
1929    do k = pver,msg + 1,-1
1930       do i = 1,ncol
1931          hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
1932          if (k >= nint(pblt(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then
1933             hmax(i) = hmn(i)
1934             mx(i) = k
1935          end if
1936       end do
1937    end do
1938 #endif
1939 !
1940    do i = 1,ncol
1941       lcl(i) = mx(i)
1942       e = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i)))
1943       tl(i) = 2840._r8/ (3.5_r8*log(t(i,mx(i)))-log(e)-4.805_r8) + 55._r8
1944       if (tl(i) < t(i,mx(i))) then
1945          plexp(i) = (1._r8/ (0.2854_r8* (1._r8-0.28_r8*q(i,mx(i)))))
1946          pl(i) = p(i,mx(i))* (tl(i)/t(i,mx(i)))**plexp(i)
1947       else
1948          tl(i) = t(i,mx(i))
1949          pl(i) = p(i,mx(i))
1950       end if
1951    end do
1952
1953 !
1954 ! calculate lifting condensation level (lcl).
1955 !
1956    do k = pver,msg + 2,-1
1957       do i = 1,ncol
1958          if (k <= mx(i) .and. (p(i,k) > pl(i) .and. p(i,k-1) <= pl(i))) then
1959             lcl(i) = k - 1
1960          end if
1961       end do
1962    end do
1963 !
1964 ! if lcl is above the nominal level of non-divergence (600 mbs),
1965 ! no deep convection is permitted (ensuing calculations
1966 ! skipped and cape retains initialized value of zero).
1967 !
1968    do i = 1,ncol
1969       plge600(i) = pl(i).ge.600._r8
1970    end do
1971 !
1972 ! initialize parcel properties in sub-cloud layer below lcl.
1973 !
1974    do k = pver,msg + 1,-1
1975       do i=1,ncol
1976          if (k > lcl(i) .and. k <= mx(i) .and. plge600(i)) then
1977             tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
1978             qstp(i,k) = q(i,mx(i))
1979             tp(i,k) = t(i,mx(i))* (p(i,k)/p(i,mx(i)))**(0.2854_r8* (1._r8-0.28_r8*q(i,mx(i))))
1980 !
1981 ! buoyancy is increased by 0.5 k as in tiedtke
1982 !
1983 !-jjh          tpv (i,k)=tp(i,k)*(1.+1.608*q(i,mx(i)))/
1984 !-jjh     1                     (1.+q(i,mx(i)))
1985             tpv(i,k) = (tp(i,k)+tpert(i))*(1._r8+1.608_r8*q(i,mx(i)))/ (1._r8+q(i,mx(i)))
1986             buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
1987          end if
1988       end do
1989    end do
1990
1991 !
1992 ! define parcel properties at lcl (i.e. level immediately above pl).
1993 !
1994    do k = pver,msg + 1,-1
1995       do i=1,ncol
1996          if (k == lcl(i) .and. plge600(i)) then
1997             tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
1998             qstp(i,k) = q(i,mx(i))
1999             tp(i,k) = tl(i)* (p(i,k)/pl(i))**(0.2854_r8* (1._r8-0.28_r8*qstp(i,k)))
2000 !              estp(i)  =exp(a-b/tp(i,k))
2001 ! use of different formulas for est has about 1 g/kg difference
2002 ! in qs at t= 300k, and 0.02 g/kg at t=263k, with the formula
2003 ! above giving larger qs.
2004 !
2005             estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/((tp(i,k)-tfreez)+c3))
2006
2007             qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i))
2008             a1(i) = cp / rl + qstp(i,k) * (1._r8+ qstp(i,k) / eps1) * rl * eps1 / &
2009                     (rd * tp(i,k) ** 2)
2010             a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/eps1*qstp(i,k))* &
2011                     (1._r8+qstp(i,k)/eps1)*eps1**2*rl*rl/ &
2012                     (rd**2*tp(i,k)**4)-qstp(i,k)* &
2013                     (1._r8+qstp(i,k)/eps1)*2._r8*eps1*rl/ &
2014                     (rd*tp(i,k)**3))
2015             a1(i) = 1._r8/a1(i)
2016             a2(i) = -a2(i)*a1(i)**3
2017             y(i) = q(i,mx(i)) - qstp(i,k)
2018             tp(i,k) = tp(i,k) + a1(i)*y(i) + a2(i)*y(i)**2
2019 !          estp(i)  =exp(a-b/tp(i,k))
2020             estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/ ((tp(i,k)-tfreez)+c3))
2021
2022             qstp(i,k) = eps1*estp(i) / (p(i,k)-estp(i))
2023 !
2024 ! buoyancy is increased by 0.5 k in cape calculation.
2025 ! dec. 9, 1994
2026 !-jjh          tpv(i,k) =tp(i,k)*(1.+1.608*qstp(i,k))/(1.+q(i,mx(i)))
2027 !
2028             tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k)) / (1._r8+q(i,mx(i)))
2029             buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
2030          end if
2031       end do
2032    end do
2033 !
2034 ! main buoyancy calculation.
2035 !
2036    do k = pver - 1,msg + 1,-1
2037       do i=1,ncol
2038          if (k < lcl(i) .and. plge600(i)) then
2039             tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
2040             qstp(i,k) = qstp(i,k+1)
2041             tp(i,k) = tp(i,k+1)* (p(i,k)/p(i,k+1))**(0.2854_r8* (1._r8-0.28_r8*qstp(i,k)))
2042 !          estp(i) = exp(a-b/tp(i,k))
2043             estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/((tp(i,k)-tfreez)+c3))
2044
2045             qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i))
2046             a1(i) = cp/rl + qstp(i,k)* (1._r8+qstp(i,k)/eps1)*rl*eps1/ (rd*tp(i,k)**2)
2047             a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/eps1*qstp(i,k))* &
2048                     (1._r8+qstp(i,k)/eps1)*eps1**2*rl*rl/ &
2049                     (rd**2*tp(i,k)**4)-qstp(i,k)* &
2050                     (1._r8+qstp(i,k)/eps1)*2._r8*eps1*rl/ &
2051                     (rd*tp(i,k)**3))
2052             a1(i) = 1._r8/a1(i)
2053             a2(i) = -a2(i)*a1(i)**3
2054             y(i) = qstp(i,k+1) - qstp(i,k)
2055             tp(i,k) = tp(i,k) + a1(i)*y(i) + a2(i)*y(i)**2
2056 !          estp(i)  =exp(a-b/tp(i,k))
2057             estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/ ((tp(i,k)-tfreez)+c3))
2058
2059             qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i))
2060 !-jjh          tpv(i,k) =tp(i,k)*(1.+1.608*qstp(i,k))/
2061 !jt            (1.+q(i,mx(i)))
2062             tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k))/(1._r8+q(i,mx(i)))
2063             buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
2064          end if
2065       end do
2066    end do
2067
2068 !
2069    do k = msg + 2,pver
2070       do i = 1,ncol
2071          if (k < lcl(i) .and. plge600(i)) then
2072             if (buoy(i,k+1) > 0._r8 .and. buoy(i,k) <= 0._r8) then
2073                knt(i) = min(5,knt(i) + 1)
2074                lelten(i,knt(i)) = k
2075             end if
2076          end if
2077       end do
2078    end do
2079 !
2080 ! calculate convective available potential energy (cape).
2081 !
2082    do n = 1,5
2083       do k = msg + 1,pver
2084          do i = 1,ncol
2085             if (plge600(i) .and. k <= mx(i) .and. k > lelten(i,n)) then
2086                capeten(i,n) = capeten(i,n) + rd*buoy(i,k)*log(pf(i,k+1)/pf(i,k))
2087             end if
2088          end do
2089       end do
2090    end do
2091 !
2092 ! find maximum cape from all possible tentative capes from
2093 ! one sounding,
2094 ! and use it as the final cape, april 26, 1995
2095 !
2096    do n = 1,5
2097       do i = 1,ncol
2098          if (capeten(i,n) > cape(i)) then
2099             cape(i) = capeten(i,n)
2100             lel(i) = lelten(i,n)
2101          end if
2102       end do
2103    end do
2104 !
2105 ! put lower bound on cape for diagnostic purposes.
2106 !
2107    do i = 1,ncol
2108       cape(i) = max(cape(i), 0._r8)
2109    end do
2110 !
2111    return
2112 end subroutine buoyan
2113
2114 subroutine cldprp(lchnk   , &
2115                   q       ,t       ,u       ,v       ,p       , &
2116                   z       ,s       ,mu      ,eu      ,du      , &
2117                   md      ,ed      ,sd      ,qd      ,mc      , &
2118                   qu      ,su      ,zf      ,qst     ,hmn     , &
2119                   hsat    ,shat    ,ql      , &
2120                   cmeg    ,jb      ,lel     ,jt      ,jlcl    , &
2121                   mx      ,j0      ,jd      ,rl      ,il2g    , &
2122                   rd      ,grav    ,cp      ,msg     , &
2123                   pflx    ,evp     ,cu      ,rprd    ,limcnv  ,landfrac)
2124 !-----------------------------------------------------------------------
2125 !
2126 ! Purpose:
2127 ! <Say what the routine does>
2128 !
2129 ! Method:
2130 ! may 09/91 - guang jun zhang, m.lazare, n.mcfarlane.
2131 !             original version cldprop.
2132 !
2133 ! Author: See above, modified by P. Rasch
2134 ! This is contributed code not fully standardized by the CCM core group.
2135 !
2136 ! this code is very much rougher than virtually anything else in the CCM
2137 ! there are debug statements left strewn about and code segments disabled
2138 ! these are to facilitate future development. We expect to release a
2139 ! cleaner code in a future release
2140 !
2141 ! the documentation has been enhanced to the degree that we are able
2142 !
2143 !-----------------------------------------------------------------------
2144
2145    implicit none
2146
2147 !------------------------------------------------------------------------------
2148 !
2149 ! Input arguments
2150 !
2151    integer, intent(in) :: lchnk                  ! chunk identifier
2152
2153    real(r8), intent(in) :: q(pcols,pver)         ! spec. humidity of env
2154    real(r8), intent(in) :: t(pcols,pver)         ! temp of env
2155    real(r8), intent(in) :: p(pcols,pver)         ! pressure of env
2156    real(r8), intent(in) :: z(pcols,pver)         ! height of env
2157    real(r8), intent(in) :: s(pcols,pver)         ! normalized dry static energy of env
2158    real(r8), intent(in) :: zf(pcols,pverp)       ! height of interfaces
2159    real(r8), intent(in) :: u(pcols,pver)         ! zonal velocity of env
2160    real(r8), intent(in) :: v(pcols,pver)         ! merid. velocity of env
2161
2162    real(r8), intent(in) :: landfrac(pcols) ! RBN Landfrac
2163
2164    integer, intent(in) :: jb(pcols)              ! updraft base level
2165    integer, intent(in) :: lel(pcols)             ! updraft launch level
2166    integer, intent(out) :: jt(pcols)              ! updraft plume top
2167    integer, intent(out) :: jlcl(pcols)            ! updraft lifting cond level
2168    integer, intent(in) :: mx(pcols)              ! updraft base level (same is jb)
2169    integer, intent(out) :: j0(pcols)              ! level where updraft begins detraining
2170    integer, intent(out) :: jd(pcols)              ! level of downdraft
2171    integer, intent(in) :: limcnv                 ! convection limiting level
2172    integer, intent(in) :: il2g                   !CORE GROUP REMOVE
2173    integer, intent(in) :: msg                    ! missing moisture vals (always 0)
2174    real(r8), intent(in) :: rl                    ! latent heat of vap
2175    real(r8), intent(in) :: shat(pcols,pver)      ! interface values of dry stat energy
2176 !
2177 ! output
2178 !
2179    real(r8), intent(out) :: rprd(pcols,pver)     ! rate of production of precip at that layer
2180    real(r8), intent(out) :: du(pcols,pver)       ! detrainement rate of updraft
2181    real(r8), intent(out) :: ed(pcols,pver)       ! entrainment rate of downdraft
2182    real(r8), intent(out) :: eu(pcols,pver)       ! entrainment rate of updraft
2183    real(r8), intent(out) :: hmn(pcols,pver)      ! moist stat energy of env
2184    real(r8), intent(out) :: hsat(pcols,pver)     ! sat moist stat energy of env
2185    real(r8), intent(out) :: mc(pcols,pver)       ! net mass flux
2186    real(r8), intent(out) :: md(pcols,pver)       ! downdraft mass flux
2187    real(r8), intent(out) :: mu(pcols,pver)       ! updraft mass flux
2188    real(r8), intent(out) :: pflx(pcols,pverp)    ! precipitation flux thru layer
2189    real(r8), intent(out) :: qd(pcols,pver)       ! spec humidity of downdraft
2190    real(r8), intent(out) :: ql(pcols,pver)       ! liq water of updraft
2191    real(r8), intent(out) :: qst(pcols,pver)      ! saturation spec humidity of env.
2192    real(r8), intent(out) :: qu(pcols,pver)       ! spec hum of updraft
2193    real(r8), intent(out) :: sd(pcols,pver)       ! normalized dry stat energy of downdraft
2194    real(r8), intent(out) :: su(pcols,pver)       ! normalized dry stat energy of updraft
2195
2196
2197    real(r8) rd                   ! gas constant for dry air
2198    real(r8) grav                 ! gravity
2199    real(r8) cp                   ! heat capacity of dry air
2200
2201 !
2202 ! Local workspace
2203 !
2204    real(r8) gamma(pcols,pver)
2205    real(r8) dz(pcols,pver)
2206    real(r8) iprm(pcols,pver)
2207    real(r8) hu(pcols,pver)
2208    real(r8) hd(pcols,pver)
2209    real(r8) eps(pcols,pver)
2210    real(r8) f(pcols,pver)
2211    real(r8) k1(pcols,pver)
2212    real(r8) i2(pcols,pver)
2213    real(r8) ihat(pcols,pver)
2214    real(r8) i3(pcols,pver)
2215    real(r8) idag(pcols,pver)
2216    real(r8) i4(pcols,pver)
2217    real(r8) qsthat(pcols,pver)
2218    real(r8) hsthat(pcols,pver)
2219    real(r8) gamhat(pcols,pver)
2220    real(r8) cu(pcols,pver)
2221    real(r8) evp(pcols,pver)
2222    real(r8) cmeg(pcols,pver)
2223    real(r8) qds(pcols,pver)
2224 ! RBN For c0mask
2225    real(r8) c0mask(pcols)
2226
2227    real(r8) hmin(pcols)
2228    real(r8) expdif(pcols)
2229    real(r8) expnum(pcols)
2230    real(r8) ftemp(pcols)
2231    real(r8) eps0(pcols)
2232    real(r8) rmue(pcols)
2233    real(r8) zuef(pcols)
2234    real(r8) zdef(pcols)
2235    real(r8) epsm(pcols)
2236    real(r8) ratmjb(pcols)
2237    real(r8) est(pcols)
2238    real(r8) totpcp(pcols)
2239    real(r8) totevp(pcols)
2240    real(r8) alfa(pcols)
2241    real(r8) ql1
2242    real(r8) tu
2243    real(r8) estu
2244    real(r8) qstu
2245
2246    real(r8) small
2247    real(r8) mdt
2248
2249    integer khighest
2250    integer klowest
2251    integer kount
2252    integer i,k
2253
2254    logical doit(pcols)
2255    logical done(pcols)
2256 !
2257 !------------------------------------------------------------------------------
2258 !
2259    do i = 1,il2g
2260       ftemp(i) = 0._r8
2261       expnum(i) = 0._r8
2262       expdif(i) = 0._r8
2263       c0mask(i)  = c0_ocn * (1._r8-landfrac(i)) +   c0_lnd * landfrac(i)
2264    end do
2265 !
2266 !jr Change from msg+1 to 1 to prevent blowup
2267 !
2268    do k = 1,pver
2269       do i = 1,il2g
2270          dz(i,k) = zf(i,k) - zf(i,k+1)
2271       end do
2272    end do
2273
2274 !
2275 ! initialize many output and work variables to zero
2276 !
2277    pflx(:il2g,1) = 0
2278
2279    do k = 1,pver
2280       do i = 1,il2g
2281          k1(i,k) = 0._r8
2282          i2(i,k) = 0._r8
2283          i3(i,k) = 0._r8
2284          i4(i,k) = 0._r8
2285          mu(i,k) = 0._r8
2286          f(i,k) = 0._r8
2287          eps(i,k) = 0._r8
2288          eu(i,k) = 0._r8
2289          du(i,k) = 0._r8
2290          ql(i,k) = 0._r8
2291          cu(i,k) = 0._r8
2292          evp(i,k) = 0._r8
2293          cmeg(i,k) = 0._r8
2294          qds(i,k) = q(i,k)
2295          md(i,k) = 0._r8
2296          ed(i,k) = 0._r8
2297          sd(i,k) = s(i,k)
2298          qd(i,k) = q(i,k)
2299          mc(i,k) = 0._r8
2300          qu(i,k) = q(i,k)
2301          su(i,k) = s(i,k)
2302 !        est(i)=exp(a-b/t(i,k))
2303          est(i) = c1*exp((c2* (t(i,k)-tfreez))/((t(i,k)-tfreez)+c3))
2304 !++bee
2305          if ( p(i,k)-est(i) > 0._r8 ) then
2306             qst(i,k) = eps1*est(i)/ (p(i,k)-est(i))
2307          else
2308             qst(i,k) = 1.0_r8
2309          end if
2310 !--bee
2311          gamma(i,k) = qst(i,k)*(1._r8 + qst(i,k)/eps1)*eps1*rl/(rd*t(i,k)**2)*rl/cp
2312          hmn(i,k) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
2313          hsat(i,k) = cp*t(i,k) + grav*z(i,k) + rl*qst(i,k)
2314          hu(i,k) = hmn(i,k)
2315          hd(i,k) = hmn(i,k)
2316          rprd(i,k) = 0._r8
2317       end do
2318    end do
2319 !
2320 !jr Set to zero things which make this routine blow up
2321 !
2322    do k=1,msg
2323       do i=1,il2g
2324          rprd(i,k) = 0._r8
2325       end do
2326    end do
2327 !
2328 ! interpolate the layer values of qst, hsat and gamma to
2329 ! layer interfaces
2330 !
2331    do i = 1,il2g
2332       hsthat(i,msg+1) = hsat(i,msg+1)
2333       qsthat(i,msg+1) = qst(i,msg+1)
2334       gamhat(i,msg+1) = gamma(i,msg+1)
2335       totpcp(i) = 0._r8
2336       totevp(i) = 0._r8
2337    end do
2338    do k = msg + 2,pver
2339       do i = 1,il2g
2340          if (abs(qst(i,k-1)-qst(i,k)) > 1.E-6_r8) then
2341             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))
2342          else
2343             qsthat(i,k) = qst(i,k)
2344          end if
2345          hsthat(i,k) = cp*shat(i,k) + rl*qsthat(i,k)
2346          if (abs(gamma(i,k-1)-gamma(i,k)) > 1.E-6_r8) then
2347             gamhat(i,k) = log(gamma(i,k-1)/gamma(i,k))*gamma(i,k-1)*gamma(i,k)/ &
2348                                 (gamma(i,k-1)-gamma(i,k))
2349          else
2350             gamhat(i,k) = gamma(i,k)
2351          end if
2352       end do
2353    end do
2354 !
2355 ! initialize cloud top to highest plume top.
2356 !jr changed hard-wired 4 to limcnv+1 (not to exceed pver)
2357 !
2358    jt(:) = pver
2359    do i = 1,il2g
2360       jt(i) = max(lel(i),limcnv+1)
2361       jt(i) = min(jt(i),pver)
2362       jd(i) = pver
2363       jlcl(i) = lel(i)
2364       hmin(i) = 1.E6_r8
2365    end do
2366 !
2367 ! find the level of minimum hsat, where detrainment starts
2368 !
2369
2370    do k = msg + 1,pver
2371       do i = 1,il2g
2372          if (hsat(i,k) <= hmin(i) .and. k >= jt(i) .and. k <= jb(i)) then
2373             hmin(i) = hsat(i,k)
2374             j0(i) = k
2375          end if
2376       end do
2377    end do
2378    do i = 1,il2g
2379       j0(i) = min(j0(i),jb(i)-2)
2380       j0(i) = max(j0(i),jt(i)+2)
2381 !
2382 ! Fix from Guang Zhang to address out of bounds array reference
2383 !
2384       j0(i) = min(j0(i),pver)
2385    end do
2386 !
2387 ! Initialize certain arrays inside cloud
2388 !
2389    do k = msg + 1,pver
2390       do i = 1,il2g
2391          if (k >= jt(i) .and. k <= jb(i)) then
2392             hu(i,k) = hmn(i,mx(i)) + cp*tiedke_add
2393             su(i,k) = s(i,mx(i)) + tiedke_add
2394          end if
2395       end do
2396    end do
2397 !
2398 ! *********************************************************
2399 ! compute taylor series for approximate eps(z) below
2400 ! *********************************************************
2401 !
2402    do k = pver - 1,msg + 1,-1
2403       do i = 1,il2g
2404          if (k < jb(i) .and. k >= jt(i)) then
2405             k1(i,k) = k1(i,k+1) + (hmn(i,mx(i))-hmn(i,k))*dz(i,k)
2406             ihat(i,k) = 0.5_r8* (k1(i,k+1)+k1(i,k))
2407             i2(i,k) = i2(i,k+1) + ihat(i,k)*dz(i,k)
2408             idag(i,k) = 0.5_r8* (i2(i,k+1)+i2(i,k))
2409             i3(i,k) = i3(i,k+1) + idag(i,k)*dz(i,k)
2410             iprm(i,k) = 0.5_r8* (i3(i,k+1)+i3(i,k))
2411             i4(i,k) = i4(i,k+1) + iprm(i,k)*dz(i,k)
2412          end if
2413       end do
2414    end do
2415 !
2416 ! re-initialize hmin array for ensuing calculation.
2417 !
2418    do i = 1,il2g
2419       hmin(i) = 1.E6_r8
2420    end do
2421    do k = msg + 1,pver
2422       do i = 1,il2g
2423          if (k >= j0(i) .and. k <= jb(i) .and. hmn(i,k) <= hmin(i)) then
2424             hmin(i) = hmn(i,k)
2425             expdif(i) = hmn(i,mx(i)) - hmin(i)
2426          end if
2427       end do
2428    end do
2429 !
2430 ! *********************************************************
2431 ! compute approximate eps(z) using above taylor series
2432 ! *********************************************************
2433 !
2434    do k = msg + 2,pver
2435       do i = 1,il2g
2436          expnum(i) = 0._r8
2437          ftemp(i) = 0._r8
2438          if (k < jt(i) .or. k >= jb(i)) then
2439             k1(i,k) = 0._r8
2440             expnum(i) = 0._r8
2441          else
2442             expnum(i) = hmn(i,mx(i)) - (hsat(i,k-1)*(zf(i,k)-z(i,k)) + &
2443                         hsat(i,k)* (z(i,k-1)-zf(i,k)))/(z(i,k-1)-z(i,k))
2444          end if
2445          if ((expdif(i) > 100._r8 .and. expnum(i) > 0._r8) .and. &
2446              k1(i,k) > expnum(i)*dz(i,k)) then
2447             ftemp(i) = expnum(i)/k1(i,k)
2448             f(i,k) = ftemp(i) + i2(i,k)/k1(i,k)*ftemp(i)**2 + &
2449                      (2._r8*i2(i,k)**2-k1(i,k)*i3(i,k))/k1(i,k)**2* &
2450                      ftemp(i)**3 + (-5._r8*k1(i,k)*i2(i,k)*i3(i,k)+ &
2451                      5._r8*i2(i,k)**3+k1(i,k)**2*i4(i,k))/ &
2452                      k1(i,k)**3*ftemp(i)**4
2453             f(i,k) = max(f(i,k),0._r8)
2454             f(i,k) = min(f(i,k),0.0002_r8)
2455          end if
2456       end do
2457    end do
2458    do i = 1,il2g
2459       if (j0(i) < jb(i)) then
2460          if (f(i,j0(i)) < 1.E-6_r8 .and. f(i,j0(i)+1) > f(i,j0(i))) j0(i) = j0(i) + 1
2461       end if
2462    end do
2463    do k = msg + 2,pver
2464       do i = 1,il2g
2465          if (k >= jt(i) .and. k <= j0(i)) then
2466             f(i,k) = max(f(i,k),f(i,k-1))
2467          end if
2468       end do
2469    end do
2470    do i = 1,il2g
2471       eps0(i) = f(i,j0(i))
2472       eps(i,jb(i)) = eps0(i)
2473    end do
2474 !
2475 ! This is set to match the Rasch and Kristjansson paper
2476 !
2477    do k = pver,msg + 1,-1
2478       do i = 1,il2g
2479          if (k >= j0(i) .and. k <= jb(i)) then
2480             eps(i,k) = f(i,j0(i))
2481          end if
2482       end do
2483    end do
2484    do k = pver,msg + 1,-1
2485       do i = 1,il2g
2486          if (k < j0(i) .and. k >= jt(i)) eps(i,k) = f(i,k)
2487       end do
2488    end do
2489 !
2490 ! specify the updraft mass flux mu, entrainment eu, detrainment du
2491 ! and moist static energy hu.
2492 ! here and below mu, eu,du, md and ed are all normalized by mb
2493 !
2494    do i = 1,il2g
2495       if (eps0(i) > 0._r8) then
2496          mu(i,jb(i)) = 1._r8
2497          eu(i,jb(i)) = mu(i,jb(i))/dz(i,jb(i))
2498       end if
2499    end do
2500    do k = pver,msg + 1,-1
2501       do i = 1,il2g
2502          if (eps0(i) > 0._r8 .and. (k >= jt(i) .and. k < jb(i))) then
2503             zuef(i) = zf(i,k) - zf(i,jb(i))
2504             rmue(i) = (1._r8/eps0(i))* (exp(eps(i,k+1)*zuef(i))-1._r8)/zuef(i)
2505             mu(i,k) = (1._r8/eps0(i))* (exp(eps(i,k  )*zuef(i))-1._r8)/zuef(i)
2506             eu(i,k) = (rmue(i)-mu(i,k+1))/dz(i,k)
2507             du(i,k) = (rmue(i)-mu(i,k))/dz(i,k)
2508          end if
2509       end do
2510    end do
2511 !
2512    khighest = pverp
2513    klowest = 1
2514    do i=1,il2g
2515       khighest = min(khighest,lel(i))
2516       klowest = max(klowest,jb(i))
2517    end do
2518    do k = klowest-1,khighest,-1
2519 !cdir$ ivdep
2520       do i = 1,il2g
2521          if (k <= jb(i)-1 .and. k >= lel(i) .and. eps0(i) > 0._r8) then
2522             if (mu(i,k) < 0.01_r8) then
2523                hu(i,k) = hu(i,jb(i))
2524                mu(i,k) = 0._r8
2525                eu(i,k) = 0._r8
2526                du(i,k) = mu(i,k+1)/dz(i,k)
2527             else
2528                hu(i,k) = mu(i,k+1)/mu(i,k)*hu(i,k+1) + &
2529                          dz(i,k)/mu(i,k)* (eu(i,k)*hmn(i,k)- du(i,k)*hsat(i,k))
2530             end if
2531          end if
2532       end do
2533    end do
2534 !
2535 ! reset cloud top index beginning from two layers above the
2536 ! cloud base (i.e. if cloud is only one layer thick, top is not reset
2537 !
2538    do i=1,il2g
2539       doit(i) = .true.
2540    end do
2541    do k=klowest-2,khighest-1,-1
2542       do i=1,il2g
2543          if (doit(i) .and. k <= jb(i)-2 .and. k >= lel(i)-1) then
2544            if (hu(i,k) <= hsthat(i,k) .and. hu(i,k+1) > hsthat(i,k+1) &
2545                .and. mu(i,k) >= 0.02_r8) then
2546                if (hu(i,k)-hsthat(i,k) < -2000._r8) then
2547                   jt(i) = k + 1
2548                   doit(i) = .false.
2549                else
2550                   jt(i) = k
2551                   if (eps0(i) <= 0._r8) doit(i) = .false.
2552                end if
2553             else if (hu(i,k) > hu(i,jb(i)) .or. mu(i,k) < 0.01_r8) then
2554                jt(i) = k + 1
2555                doit(i) = .false.
2556             end if
2557          end if
2558       end do
2559    end do
2560    do k = pver,msg + 1,-1
2561 !cdir$ ivdep
2562       do i = 1,il2g
2563          if (k >= lel(i) .and. k <= jt(i) .and. eps0(i) > 0._r8) then
2564             mu(i,k) = 0._r8
2565             eu(i,k) = 0._r8
2566             du(i,k) = 0._r8
2567             hu(i,k) = hu(i,jb(i))
2568          end if
2569          if (k == jt(i) .and. eps0(i) > 0._r8) then
2570             du(i,k) = mu(i,k+1)/dz(i,k)
2571             eu(i,k) = 0._r8
2572             mu(i,k) = 0._r8
2573          end if
2574       end do
2575    end do
2576 !
2577 ! specify downdraft properties (no downdrafts if jd.ge.jb).
2578 ! scale down downward mass flux profile so that net flux
2579 ! (up-down) at cloud base in not negative.
2580 !
2581    do i = 1,il2g
2582 !
2583 ! in normal downdraft strength run alfa=0.2.  In test4 alfa=0.1
2584 !
2585       alfa(i) = 0.1_r8
2586       jt(i) = min(jt(i),jb(i)-1)
2587       jd(i) = max(j0(i),jt(i)+1)
2588       jd(i) = min(jd(i),jb(i))
2589       hd(i,jd(i)) = hmn(i,jd(i)-1)
2590       if (jd(i) < jb(i) .and. eps0(i) > 0._r8) then
2591          epsm(i) = eps0(i)
2592          md(i,jd(i)) = -alfa(i)*epsm(i)/eps0(i)
2593       end if
2594    end do
2595    do k = msg + 1,pver
2596       do i = 1,il2g
2597          if ((k > jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8) then
2598             zdef(i) = zf(i,jd(i)) - zf(i,k)
2599             md(i,k) = -alfa(i)/ (2._r8*eps0(i))*(exp(2._r8*epsm(i)*zdef(i))-1._r8)/zdef(i)
2600          end if
2601       end do
2602    end do
2603    do k = msg + 1,pver
2604       do i = 1,il2g
2605          if ((k >= jt(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8 .and. jd(i) < jb(i)) then
2606             ratmjb(i) = min(abs(mu(i,jb(i))/md(i,jb(i))),1._r8)
2607             md(i,k) = md(i,k)*ratmjb(i)
2608          end if
2609       end do
2610    end do
2611
2612    small = 1.e-20_r8
2613    do k = msg + 1,pver
2614       do i = 1,il2g
2615          if ((k >= jt(i) .and. k <= pver) .and. eps0(i) > 0._r8) then
2616             ed(i,k-1) = (md(i,k-1)-md(i,k))/dz(i,k-1)
2617             mdt = min(md(i,k),-small)
2618             hd(i,k) = (md(i,k-1)*hd(i,k-1) - dz(i,k-1)*ed(i,k-1)*hmn(i,k-1))/mdt
2619          end if
2620       end do
2621    end do
2622 !
2623 ! calculate updraft and downdraft properties.
2624 !
2625    do k = msg + 2,pver
2626       do i = 1,il2g
2627          if ((k >= jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8 .and. jd(i) < jb(i)) then
2628             qds(i,k) = qsthat(i,k) + gamhat(i,k)*(hd(i,k)-hsthat(i,k))/ &
2629                (rl*(1._r8 + gamhat(i,k)))
2630          end if
2631       end do
2632    end do
2633 !
2634    do i = 1,il2g
2635       done(i) = .false.
2636    end do
2637    kount = 0
2638    do k = pver,msg + 2,-1
2639       do i = 1,il2g
2640          if (( .not. done(i) .and. k > jt(i) .and. k < jb(i)) .and. eps0(i) > 0._r8) then
2641             su(i,k) = mu(i,k+1)/mu(i,k)*su(i,k+1) + &
2642                       dz(i,k)/mu(i,k)* (eu(i,k)-du(i,k))*s(i,k)
2643             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)- &
2644                             du(i,k)*qst(i,k))
2645             tu = su(i,k) - grav/cp*zf(i,k)
2646             estu = c1*exp((c2* (tu-tfreez))/ ((tu-tfreez)+c3))
2647             qstu = eps1*estu/ ((p(i,k)+p(i,k-1))/2._r8-estu)
2648             if (qu(i,k) >= qstu) then
2649                jlcl(i) = k
2650                kount = kount + 1
2651                done(i) = .true.
2652             end if
2653          end if
2654       end do
2655       if (kount >= il2g) goto 690
2656    end do
2657 690 continue
2658    do k = msg + 2,pver
2659       do i = 1,il2g
2660          if (k == jb(i) .and. eps0(i) > 0._r8) then
2661             qu(i,k) = q(i,mx(i))
2662             su(i,k) = (hu(i,k)-rl*qu(i,k))/cp
2663          end if
2664          if ((k > jt(i) .and. k <= jlcl(i)) .and. eps0(i) > 0._r8) then
2665             su(i,k) = shat(i,k) + (hu(i,k)-hsthat(i,k))/(cp* (1._r8+gamhat(i,k)))
2666             qu(i,k) = qsthat(i,k) + gamhat(i,k)*(hu(i,k)-hsthat(i,k))/ &
2667                      (rl* (1._r8+gamhat(i,k)))
2668          end if
2669       end do
2670    end do
2671
2672 ! compute condensation in updraft
2673    do k = pver,msg + 2,-1
2674       do i = 1,il2g
2675          if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0._r8) then
2676             cu(i,k) = ((mu(i,k)*su(i,k)-mu(i,k+1)*su(i,k+1))/ &
2677                       dz(i,k)- (eu(i,k)-du(i,k))*s(i,k))/(rl/cp)
2678             if (k == jt(i)) cu(i,k) = 0._r8
2679             cu(i,k) = max(0._r8,cu(i,k))
2680          end if
2681       end do
2682    end do
2683
2684 ! compute condensed liquid, rain production rate
2685 ! accumulate total precipitation (condensation - detrainment of liquid)
2686 ! Note ql1 = ql(k) + rprd(k)*dz(k)/mu(k)
2687 ! The differencing is somewhat strange (e.g. du(i,k)*ql(i,k+1)) but is
2688 ! consistently applied.
2689 !    mu, ql are interface quantities
2690 !    cu, du, eu, rprd are midpoint quantites
2691    do k = pver,msg + 2,-1
2692       do i = 1,il2g
2693          rprd(i,k) = 0._r8
2694          if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0._r8 .and. mu(i,k) >= 0.0_r8) then
2695             if (mu(i,k) > 0._r8) then
2696                ql1 = 1._r8/mu(i,k)* (mu(i,k+1)*ql(i,k+1)- &
2697                      dz(i,k)*du(i,k)*ql(i,k+1)+dz(i,k)*cu(i,k))
2698                ql(i,k) = ql1/ (1._r8+dz(i,k)*c0mask(i))
2699             else
2700                ql(i,k) = 0._r8
2701             end if
2702             totpcp(i) = totpcp(i) + dz(i,k)*(cu(i,k)-du(i,k)*ql(i,k+1))
2703             rprd(i,k) = c0mask(i)*mu(i,k)*ql(i,k)
2704          end if
2705       end do
2706    end do
2707 !
2708    do i = 1,il2g
2709       qd(i,jd(i)) = qds(i,jd(i))
2710       sd(i,jd(i)) = (hd(i,jd(i)) - rl*qd(i,jd(i)))/cp
2711    end do
2712 !
2713    do k = msg + 2,pver
2714       do i = 1,il2g
2715          if (k >= jd(i) .and. k < jb(i) .and. eps0(i) > 0._r8) then
2716             qd(i,k+1) = qds(i,k+1)
2717             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)
2718             evp(i,k) = max(evp(i,k),0._r8)
2719             mdt = min(md(i,k+1),-small)
2720             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
2721             totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k)
2722          end if
2723       end do
2724    end do
2725    do i = 1,il2g
2726 !*guang         totevp(i) = totevp(i) + md(i,jd(i))*q(i,jd(i)-1) -
2727       totevp(i) = totevp(i) + md(i,jd(i))*qd(i,jd(i)) - md(i,jb(i))*qd(i,jb(i))
2728    end do
2729 !!$   if (.true.) then
2730    if (.false.) then
2731       do i = 1,il2g
2732          k = jb(i)
2733          if (eps0(i) > 0._r8) then
2734             evp(i,k) = -ed(i,k)*q(i,k) + (md(i,k)*qd(i,k))/dz(i,k)
2735             evp(i,k) = max(evp(i,k),0._r8)
2736             totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k)
2737          end if
2738       end do
2739    endif
2740
2741    do i = 1,il2g
2742       totpcp(i) = max(totpcp(i),0._r8)
2743       totevp(i) = max(totevp(i),0._r8)
2744    end do
2745 !
2746    do k = msg + 2,pver
2747       do i = 1,il2g
2748          if (totevp(i) > 0._r8 .and. totpcp(i) > 0._r8) then
2749             md(i,k)  = md (i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i)))
2750             ed(i,k)  = ed (i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i)))
2751             evp(i,k) = evp(i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i)))
2752          else
2753             md(i,k) = 0._r8
2754             ed(i,k) = 0._r8
2755             evp(i,k) = 0._r8
2756          end if
2757 ! cmeg is the cloud water condensed - rain water evaporated
2758 ! rprd is the cloud water converted to rain - (rain evaporated)
2759          cmeg(i,k) = cu(i,k) - evp(i,k)
2760          rprd(i,k) = rprd(i,k)-evp(i,k)
2761       end do
2762    end do
2763
2764 ! compute the net precipitation flux across interfaces
2765    pflx(:il2g,1) = 0._r8
2766    do k = 2,pverp
2767       do i = 1,il2g
2768          pflx(i,k) = pflx(i,k-1) + rprd(i,k-1)*dz(i,k-1)
2769       end do
2770    end do
2771 !
2772    do k = msg + 1,pver
2773       do i = 1,il2g
2774          mc(i,k) = mu(i,k) + md(i,k)
2775       end do
2776    end do
2777 !
2778    return
2779 end subroutine cldprp
2780
2781 subroutine closure(lchnk   , &
2782                    q       ,t       ,p       ,z       ,s       , &
2783                    tp      ,qs      ,qu      ,su      ,mc      , &
2784                    du      ,mu      ,md      ,qd      ,sd      , &
2785                    qhat    ,shat    ,dp      ,qstp    ,zf      , &
2786                    ql      ,dsubcld ,mb      ,cape    ,tl      , &
2787                    lcl     ,lel     ,jt      ,mx      ,il1g    , &
2788                    il2g    ,rd      ,grav    ,cp      ,rl      , &
2789                    msg     ,capelmt )
2790 !-----------------------------------------------------------------------
2791 !
2792 ! Purpose:
2793 ! <Say what the routine does>
2794 !
2795 ! Method:
2796 ! <Describe the algorithm(s) used in the routine.>
2797 ! <Also include any applicable external references.>
2798 !
2799 ! Author: G. Zhang and collaborators. CCM contact:P. Rasch
2800 ! This is contributed code not fully standardized by the CCM core group.
2801 !
2802 ! this code is very much rougher than virtually anything else in the CCM
2803 ! We expect to release cleaner code in a future release
2804 !
2805 ! the documentation has been enhanced to the degree that we are able
2806 !
2807 !-----------------------------------------------------------------------
2808 #ifndef WRF_PORT
2809    use dycore,    only: dycore_is, get_resolution
2810 #endif
2811
2812    implicit none
2813
2814 !
2815 !-----------------------------Arguments---------------------------------
2816 !
2817    integer, intent(in) :: lchnk                 ! chunk identifier
2818
2819    real(r8), intent(inout) :: q(pcols,pver)        ! spec humidity
2820    real(r8), intent(inout) :: t(pcols,pver)        ! temperature
2821    real(r8), intent(inout) :: p(pcols,pver)        ! pressure (mb)
2822    real(r8), intent(inout) :: mb(pcols)            ! cloud base mass flux
2823    real(r8), intent(in) :: z(pcols,pver)        ! height (m)
2824    real(r8), intent(in) :: s(pcols,pver)        ! normalized dry static energy
2825    real(r8), intent(in) :: tp(pcols,pver)       ! parcel temp
2826    real(r8), intent(in) :: qs(pcols,pver)       ! sat spec humidity
2827    real(r8), intent(in) :: qu(pcols,pver)       ! updraft spec. humidity
2828    real(r8), intent(in) :: su(pcols,pver)       ! normalized dry stat energy of updraft
2829    real(r8), intent(in) :: mc(pcols,pver)       ! net convective mass flux
2830    real(r8), intent(in) :: du(pcols,pver)       ! detrainment from updraft
2831    real(r8), intent(in) :: mu(pcols,pver)       ! mass flux of updraft
2832    real(r8), intent(in) :: md(pcols,pver)       ! mass flux of downdraft
2833    real(r8), intent(in) :: qd(pcols,pver)       ! spec. humidity of downdraft
2834    real(r8), intent(in) :: sd(pcols,pver)       ! dry static energy of downdraft
2835    real(r8), intent(in) :: qhat(pcols,pver)     ! environment spec humidity at interfaces
2836    real(r8), intent(in) :: shat(pcols,pver)     ! env. normalized dry static energy at intrfcs
2837    real(r8), intent(in) :: dp(pcols,pver)       ! pressure thickness of layers
2838    real(r8), intent(in) :: qstp(pcols,pver)     ! spec humidity of parcel
2839    real(r8), intent(in) :: zf(pcols,pver+1)     ! height of interface levels
2840    real(r8), intent(in) :: ql(pcols,pver)       ! liquid water mixing ratio
2841
2842    real(r8), intent(in) :: cape(pcols)          ! available pot. energy of column
2843    real(r8), intent(in) :: tl(pcols)
2844    real(r8), intent(in) :: dsubcld(pcols)       ! thickness of subcloud layer
2845
2846    integer, intent(in) :: lcl(pcols)        ! index of lcl
2847    integer, intent(in) :: lel(pcols)        ! index of launch leve
2848    integer, intent(in) :: jt(pcols)         ! top of updraft
2849    integer, intent(in) :: mx(pcols)         ! base of updraft
2850 !
2851 !--------------------------Local variables------------------------------
2852 !
2853    real(r8) dtpdt(pcols,pver)
2854    real(r8) dqsdtp(pcols,pver)
2855    real(r8) dtmdt(pcols,pver)
2856    real(r8) dqmdt(pcols,pver)
2857    real(r8) dboydt(pcols,pver)
2858    real(r8) thetavp(pcols,pver)
2859    real(r8) thetavm(pcols,pver)
2860
2861    real(r8) dtbdt(pcols),dqbdt(pcols),dtldt(pcols)
2862    real(r8) beta
2863    real(r8) capelmt
2864    real(r8) cp
2865    real(r8) dadt(pcols)
2866    real(r8) debdt
2867    real(r8) dltaa
2868    real(r8) eb
2869    real(r8) grav
2870
2871    integer i
2872    integer il1g
2873    integer il2g
2874    integer k, kmin, kmax
2875    integer msg
2876
2877    real(r8) rd
2878    real(r8) rl
2879 ! change of subcloud layer properties due to convection is
2880 ! related to cumulus updrafts and downdrafts.
2881 ! mc(z)=f(z)*mb, mub=betau*mb, mdb=betad*mb are used
2882 ! to define betau, betad and f(z).
2883 ! note that this implies all time derivatives are in effect
2884 ! time derivatives per unit cloud-base mass flux, i.e. they
2885 ! have units of 1/mb instead of 1/sec.
2886 !
2887    do i = il1g,il2g
2888       mb(i) = 0._r8
2889       eb = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i)))
2890       dtbdt(i) = (1._r8/dsubcld(i))* (mu(i,mx(i))*(shat(i,mx(i))-su(i,mx(i)))+ &
2891                   md(i,mx(i))* (shat(i,mx(i))-sd(i,mx(i))))
2892       dqbdt(i) = (1._r8/dsubcld(i))* (mu(i,mx(i))*(qhat(i,mx(i))-qu(i,mx(i)))+ &
2893                  md(i,mx(i))* (qhat(i,mx(i))-qd(i,mx(i))))
2894       debdt = eps1*p(i,mx(i))/ (eps1+q(i,mx(i)))**2*dqbdt(i)
2895       dtldt(i) = -2840._r8* (3.5_r8/t(i,mx(i))*dtbdt(i)-debdt/eb)/ &
2896                  (3.5_r8*log(t(i,mx(i)))-log(eb)-4.805_r8)**2
2897    end do
2898 !
2899 !   dtmdt and dqmdt are cumulus heating and drying.
2900 !
2901    do k = msg + 1,pver
2902       do i = il1g,il2g
2903          dtmdt(i,k) = 0._r8
2904          dqmdt(i,k) = 0._r8
2905       end do
2906    end do
2907 !
2908    do k = msg + 1,pver - 1
2909       do i = il1g,il2g
2910          if (k == jt(i)) then
2911             dtmdt(i,k) = (1._r8/dp(i,k))*(mu(i,k+1)* (su(i,k+1)-shat(i,k+1)- &
2912                           rl/cp*ql(i,k+1))+md(i,k+1)* (sd(i,k+1)-shat(i,k+1)))
2913             dqmdt(i,k) = (1._r8/dp(i,k))*(mu(i,k+1)* (qu(i,k+1)- &
2914                          qhat(i,k+1)+ql(i,k+1))+md(i,k+1)*(qd(i,k+1)-qhat(i,k+1)))
2915          end if
2916       end do
2917    end do
2918 !
2919    beta = 0._r8
2920    do k = msg + 1,pver - 1
2921       do i = il1g,il2g
2922          if (k > jt(i) .and. k < mx(i)) then
2923             dtmdt(i,k) = (mc(i,k)* (shat(i,k)-s(i,k))+mc(i,k+1)* (s(i,k)-shat(i,k+1)))/ &
2924                          dp(i,k) - rl/cp*du(i,k)*(beta*ql(i,k)+ (1-beta)*ql(i,k+1))
2925 !          dqmdt(i,k)=(mc(i,k)*(qhat(i,k)-q(i,k))
2926 !     1                +mc(i,k+1)*(q(i,k)-qhat(i,k+1)))/dp(i,k)
2927 !     2                +du(i,k)*(qs(i,k)-q(i,k))
2928 !     3                +du(i,k)*(beta*ql(i,k)+(1-beta)*ql(i,k+1))
2929
2930             dqmdt(i,k) = (mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)+cp/rl* (su(i,k+1)-s(i,k)))- &
2931                           mu(i,k)* (qu(i,k)-qhat(i,k)+cp/rl*(su(i,k)-s(i,k)))+md(i,k+1)* &
2932                          (qd(i,k+1)-qhat(i,k+1)+cp/rl*(sd(i,k+1)-s(i,k)))-md(i,k)* &
2933                          (qd(i,k)-qhat(i,k)+cp/rl*(sd(i,k)-s(i,k))))/dp(i,k) + &
2934                           du(i,k)* (beta*ql(i,k)+(1-beta)*ql(i,k+1))
2935          end if
2936       end do
2937    end do
2938 !
2939    do k = msg + 1,pver
2940       do i = il1g,il2g
2941          if (k >= lel(i) .and. k <= lcl(i)) then
2942             thetavp(i,k) = tp(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+1.608_r8*qstp(i,k)-q(i,mx(i)))
2943             thetavm(i,k) = t(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,k))
2944             dqsdtp(i,k) = qstp(i,k)* (1._r8+qstp(i,k)/eps1)*eps1*rl/(rd*tp(i,k)**2)
2945 !
2946 ! dtpdt is the parcel temperature change due to change of
2947 ! subcloud layer properties during convection.
2948 !
2949             dtpdt(i,k) = tp(i,k)/ (1._r8+rl/cp* (dqsdtp(i,k)-qstp(i,k)/tp(i,k)))* &
2950                         (dtbdt(i)/t(i,mx(i))+rl/cp* (dqbdt(i)/tl(i)-q(i,mx(i))/ &
2951                          tl(i)**2*dtldt(i)))
2952 !
2953 ! dboydt is the integrand of cape change.
2954 !
2955             dboydt(i,k) = ((dtpdt(i,k)/tp(i,k)+1._r8/(1._r8+1.608_r8*qstp(i,k)-q(i,mx(i)))* &
2956                           (1.608_r8 * dqsdtp(i,k) * dtpdt(i,k) -dqbdt(i))) - (dtmdt(i,k)/t(i,k)+0.608_r8/ &
2957                           (1._r8+0.608_r8*q(i,k))*dqmdt(i,k)))*grav*thetavp(i,k)/thetavm(i,k)
2958          end if
2959       end do
2960    end do
2961 !
2962    do k = msg + 1,pver
2963       do i = il1g,il2g
2964          if (k > lcl(i) .and. k < mx(i)) then
2965             thetavp(i,k) = tp(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,mx(i)))
2966             thetavm(i,k) = t(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,k))
2967 !
2968 ! dboydt is the integrand of cape change.
2969 !
2970             dboydt(i,k) = (dtbdt(i)/t(i,mx(i))+0.608_r8/ (1._r8+0.608_r8*q(i,mx(i)))*dqbdt(i)- &
2971                           dtmdt(i,k)/t(i,k)-0.608_r8/ (1._r8+0.608_r8*q(i,k))*dqmdt(i,k))* &
2972                           grav*thetavp(i,k)/thetavm(i,k)
2973          end if
2974       end do
2975    end do
2976
2977 !
2978 ! buoyant energy change is set to 2/3*excess cape per 3 hours
2979 !
2980    dadt(il1g:il2g)  = 0._r8
2981    kmin = minval(lel(il1g:il2g))
2982    kmax = maxval(mx(il1g:il2g)) - 1
2983    do k = kmin, kmax
2984       do i = il1g,il2g
2985          if ( k >= lel(i) .and. k <= mx(i) - 1) then
2986             dadt(i) = dadt(i) + dboydt(i,k)* (zf(i,k)-zf(i,k+1))
2987          endif
2988       end do
2989    end do
2990    do i = il1g,il2g
2991       dltaa = -1._r8* (cape(i)-capelmt)
2992       if (dadt(i) /= 0._r8) mb(i) = max(dltaa/tau/dadt(i),0._r8)
2993    end do
2994 !
2995    return
2996 end subroutine closure
2997
2998 subroutine q1q2_pjr(lchnk   , &
2999                     dqdt    ,dsdt    ,q       ,qs      ,qu      , &
3000                     su      ,du      ,qhat    ,shat    ,dp      , &
3001                     mu      ,md      ,sd      ,qd      ,ql      , &
3002                     dsubcld ,jt      ,mx      ,il1g    ,il2g    , &
3003                     cp      ,rl      ,msg     ,          &
3004                     dl      ,evp     ,cu      )
3005
3006
3007    implicit none
3008
3009 !-----------------------------------------------------------------------
3010 !
3011 ! Purpose:
3012 ! <Say what the routine does>
3013 !
3014 ! Method:
3015 ! <Describe the algorithm(s) used in the routine.>
3016 ! <Also include any applicable external references.>
3017 !
3018 ! Author: phil rasch dec 19 1995
3019 !
3020 !-----------------------------------------------------------------------
3021
3022
3023    real(r8), intent(in) :: cp
3024
3025    integer, intent(in) :: lchnk             ! chunk identifier
3026    integer, intent(in) :: il1g
3027    integer, intent(in) :: il2g
3028    integer, intent(in) :: msg
3029
3030    real(r8), intent(in) :: q(pcols,pver)
3031    real(r8), intent(in) :: qs(pcols,pver)
3032    real(r8), intent(in) :: qu(pcols,pver)
3033    real(r8), intent(in) :: su(pcols,pver)
3034    real(r8), intent(in) :: du(pcols,pver)
3035    real(r8), intent(in) :: qhat(pcols,pver)
3036    real(r8), intent(in) :: shat(pcols,pver)
3037    real(r8), intent(in) :: dp(pcols,pver)
3038    real(r8), intent(in) :: mu(pcols,pver)
3039    real(r8), intent(in) :: md(pcols,pver)
3040    real(r8), intent(in) :: sd(pcols,pver)
3041    real(r8), intent(in) :: qd(pcols,pver)
3042    real(r8), intent(in) :: ql(pcols,pver)
3043    real(r8), intent(in) :: evp(pcols,pver)
3044    real(r8), intent(in) :: cu(pcols,pver)
3045    real(r8), intent(in) :: dsubcld(pcols)
3046
3047    real(r8),intent(out) :: dqdt(pcols,pver),dsdt(pcols,pver)
3048    real(r8),intent(out) :: dl(pcols,pver)
3049    integer kbm
3050    integer ktm
3051    integer jt(pcols)
3052    integer mx(pcols)
3053 !
3054 ! work fields:
3055 !
3056    integer i
3057    integer k
3058
3059    real(r8) emc
3060    real(r8) rl
3061 !-------------------------------------------------------------------
3062    do k = msg + 1,pver
3063       do i = il1g,il2g
3064          dsdt(i,k) = 0._r8
3065          dqdt(i,k) = 0._r8
3066          dl(i,k) = 0._r8
3067       end do
3068    end do
3069 !
3070 ! find the highest level top and bottom levels of convection
3071 !
3072    ktm = pver
3073    kbm = pver
3074    do i = il1g, il2g
3075       ktm = min(ktm,jt(i))
3076       kbm = min(kbm,mx(i))
3077    end do
3078
3079    do k = ktm,pver-1
3080       do i = il1g,il2g
3081          emc = -cu (i,k)               &         ! condensation in updraft
3082                +evp(i,k)                         ! evaporating rain in downdraft
3083
3084          dsdt(i,k) = -rl/cp*emc &
3085                      + (+mu(i,k+1)* (su(i,k+1)-shat(i,k+1)) &
3086                         -mu(i,k)*   (su(i,k)-shat(i,k)) &
3087                         +md(i,k+1)* (sd(i,k+1)-shat(i,k+1)) &
3088                         -md(i,k)*   (sd(i,k)-shat(i,k)) &
3089                        )/dp(i,k)
3090
3091          dqdt(i,k) = emc + &
3092                     (+mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)) &
3093                      -mu(i,k)*   (qu(i,k)-qhat(i,k)) &
3094                      +md(i,k+1)* (qd(i,k+1)-qhat(i,k+1)) &
3095                      -md(i,k)*   (qd(i,k)-qhat(i,k)) &
3096                     )/dp(i,k)
3097
3098          dl(i,k) = du(i,k)*ql(i,k+1)
3099
3100       end do
3101    end do
3102
3103 !
3104 !DIR$ NOINTERCHANGE!
3105    do k = kbm,pver
3106       do i = il1g,il2g
3107          if (k == mx(i)) then
3108             dsdt(i,k) = (1._r8/dsubcld(i))* &
3109                         (-mu(i,k)* (su(i,k)-shat(i,k)) &
3110                          -md(i,k)* (sd(i,k)-shat(i,k)) &
3111                         )
3112             dqdt(i,k) = (1._r8/dsubcld(i))* &
3113                         (-mu(i,k)*(qu(i,k)-qhat(i,k)) &
3114                          -md(i,k)*(qd(i,k)-qhat(i,k)) &
3115                         )
3116          else if (k > mx(i)) then
3117             dsdt(i,k) = dsdt(i,k-1)
3118             dqdt(i,k) = dqdt(i,k-1)
3119          end if
3120       end do
3121    end do
3122 !
3123    return
3124 end subroutine q1q2_pjr
3125
3126 subroutine buoyan_dilute(lchnk   ,ncol    , &
3127                   q       ,t       ,p       ,z       ,pf      , &
3128                   tp      ,qstp    ,tl      ,rl      ,cape    , &
3129                   pblt    ,lcl     ,lel     ,lon     ,mx      , &
3130                   rd      ,grav    ,cp      ,msg     , &
3131                   tpert   )
3132 !-----------------------------------------------------------------------
3133 !
3134 ! Purpose:
3135 ! Calculates CAPE the lifting condensation level and the convective top
3136 ! where buoyancy is first -ve.
3137 !
3138 ! Method: Calculates the parcel temperature based on a simple constant
3139 ! entraining plume model. CAPE is integrated from buoyancy.
3140 ! 09/09/04 - Simplest approach using an assumed entrainment rate for
3141 !            testing (dmpdp).
3142 ! 08/04/05 - Swap to convert dmpdz to dmpdp
3143 !
3144 ! SCAM Logical Switches - DILUTE:RBN - Now Disabled
3145 ! ---------------------
3146 ! switch(1) = .T. - Uses the dilute parcel calculation to obtain tendencies.
3147 ! switch(2) = .T. - Includes entropy/q changes due to condensate loss and freezing.
3148 ! switch(3) = .T. - Adds the PBL Tpert for the parcel temperature at all levels.
3149 !
3150 ! References:
3151 ! Raymond and Blythe (1992) JAS
3152 !
3153 ! Author:
3154 ! Richard Neale - September 2004
3155 !
3156 !-----------------------------------------------------------------------
3157    implicit none
3158 !-----------------------------------------------------------------------
3159 !
3160 ! input arguments
3161 !
3162    integer, intent(in) :: lchnk                 ! chunk identifier
3163    integer, intent(in) :: ncol                  ! number of atmospheric columns
3164
3165    real(r8), intent(in) :: q(pcols,pver)        ! spec. humidity
3166    real(r8), intent(in) :: t(pcols,pver)        ! temperature
3167    real(r8), intent(in) :: p(pcols,pver)        ! pressure
3168    real(r8), intent(in) :: z(pcols,pver)        ! height
3169    real(r8), intent(in) :: pf(pcols,pver+1)     ! pressure at interfaces
3170    real(r8), intent(in) :: pblt(pcols)          ! index of pbl depth
3171    real(r8), intent(in) :: tpert(pcols)         ! perturbation temperature by pbl processes
3172
3173 !
3174 ! output arguments
3175 !
3176    real(r8), intent(out) :: tp(pcols,pver)       ! parcel temperature
3177    real(r8), intent(out) :: qstp(pcols,pver)     ! saturation mixing ratio of parcel (only above lcl, just q below).
3178    real(r8), intent(out) :: tl(pcols)            ! parcel temperature at lcl
3179    real(r8), intent(out) :: cape(pcols)          ! convective aval. pot. energy.
3180    integer lcl(pcols)        !
3181    integer lel(pcols)        !
3182    integer lon(pcols)        ! level of onset of deep convection
3183    integer mx(pcols)         ! level of max moist static energy
3184 !
3185 !--------------------------Local Variables------------------------------
3186 !
3187    real(r8) capeten(pcols,5)     ! provisional value of cape
3188    real(r8) tv(pcols,pver)       !
3189    real(r8) tpv(pcols,pver)      !
3190    real(r8) buoy(pcols,pver)
3191
3192    real(r8) a1(pcols)
3193    real(r8) a2(pcols)
3194    real(r8) estp(pcols)
3195    real(r8) pl(pcols)
3196    real(r8) plexp(pcols)
3197    real(r8) hmax(pcols)
3198    real(r8) hmn(pcols)
3199    real(r8) y(pcols)
3200
3201    logical plge600(pcols)
3202    integer knt(pcols)
3203    integer lelten(pcols,5)
3204
3205    real(r8) cp
3206    real(r8) e
3207    real(r8) grav
3208
3209    integer i
3210    integer k
3211    integer msg
3212    integer n
3213
3214    real(r8) rd
3215    real(r8) rl
3216 #ifdef PERGRO
3217    real(r8) rhd
3218 #endif
3219 !
3220 !-----------------------------------------------------------------------
3221 !
3222    do n = 1,5
3223       do i = 1,ncol
3224          lelten(i,n) = pver
3225          capeten(i,n) = 0._r8
3226       end do
3227    end do
3228 !
3229    do i = 1,ncol
3230       lon(i) = pver
3231       knt(i) = 0
3232       lel(i) = pver
3233       mx(i) = lon(i)
3234       cape(i) = 0._r8
3235       hmax(i) = 0._r8
3236    end do
3237
3238    tp(:ncol,:) = t(:ncol,:)
3239    qstp(:ncol,:) = q(:ncol,:)
3240
3241 !!! RBN - Initialize tv and buoy for output.
3242 !!! tv=tv : tpv=tpv : qstp=q : buoy=0.
3243    tv(:ncol,:) = t(:ncol,:) *(1._r8+1.608_r8*q(:ncol,:))/ (1._r8+q(:ncol,:))
3244    tpv(:ncol,:) = tv(:ncol,:)
3245    buoy(:ncol,:) = 0._r8
3246
3247 !
3248 ! set "launching" level(mx) to be at maximum moist static energy.
3249 ! search for this level stops at planetary boundary layer top.
3250 !
3251 #ifdef PERGRO
3252    do k = pver,msg + 1,-1
3253       do i = 1,ncol
3254          hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
3255 !
3256 ! Reset max moist static energy level when relative difference exceeds 1.e-4
3257 !
3258          rhd = (hmn(i) - hmax(i))/(hmn(i) + hmax(i))
3259          if (k >= nint(pblt(i)) .and. k <= lon(i) .and. rhd > -1.e-4_r8) then
3260             hmax(i) = hmn(i)
3261             mx(i) = k
3262          end if
3263       end do
3264    end do
3265 #else
3266    do k = pver,msg + 1,-1
3267       do i = 1,ncol
3268          hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
3269          if (k >= nint(pblt(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then
3270             hmax(i) = hmn(i)
3271             mx(i) = k
3272          end if
3273       end do
3274    end do
3275 #endif
3276
3277 ! LCL dilute calculation - initialize to mx(i)
3278 ! Determine lcl in parcel_dilute and get pl,tl after parcel_dilute
3279 ! Original code actually sets LCL as level above wher condensate forms.
3280 ! Therefore in parcel_dilute lcl(i) will be at first level where qsmix < qtmix.
3281
3282    do i = 1,ncol ! Initialise LCL variables.
3283       lcl(i) = mx(i)
3284       tl(i) = t(i,mx(i))
3285       pl(i) = p(i,mx(i))
3286    end do
3287
3288 !
3289 ! main buoyancy calculation.
3290 !
3291 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3292 !!! DILUTE PLUME CALCULATION USING ENTRAINING PLUME !!!
3293 !!!   RBN 9/9/04   !!!
3294
3295    call parcel_dilute(lchnk, ncol, msg, mx, p, t, q, tpert, tp, tpv, qstp, pl, tl, lcl)
3296
3297
3298 ! If lcl is above the nominal level of non-divergence (600 mbs),
3299 ! no deep convection is permitted (ensuing calculations
3300 ! skipped and cape retains initialized value of zero).
3301 !
3302    do i = 1,ncol
3303       plge600(i) = pl(i).ge.600._r8 ! Just change to always allow buoy calculation.
3304    end do
3305
3306 !
3307 ! Main buoyancy calculation.
3308 !
3309    do k = pver,msg + 1,-1
3310       do i=1,ncol
3311          if (k <= mx(i) .and. plge600(i)) then   ! Define buoy from launch level to cloud top.
3312             tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
3313             buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add  ! +0.5K or not?
3314          else
3315             qstp(i,k) = q(i,k)
3316             tp(i,k)   = t(i,k)
3317             tpv(i,k)  = tv(i,k)
3318          endif
3319       end do
3320    end do
3321
3322
3323
3324 !-------------------------------------------------------------------------------
3325
3326 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3327
3328
3329
3330 !
3331    do k = msg + 2,pver
3332       do i = 1,ncol
3333          if (k < lcl(i) .and. plge600(i)) then
3334             if (buoy(i,k+1) > 0. .and. buoy(i,k) <= 0._r8) then
3335                knt(i) = min(5,knt(i) + 1)
3336                lelten(i,knt(i)) = k
3337             end if
3338          end if
3339       end do
3340    end do
3341 !
3342 ! calculate convective available potential energy (cape).
3343 !
3344    do n = 1,5
3345       do k = msg + 1,pver
3346          do i = 1,ncol
3347             if (plge600(i) .and. k <= mx(i) .and. k > lelten(i,n)) then
3348                capeten(i,n) = capeten(i,n) + rd*buoy(i,k)*log(pf(i,k+1)/pf(i,k))
3349             end if
3350          end do
3351       end do
3352    end do
3353 !
3354 ! find maximum cape from all possible tentative capes from
3355 ! one sounding,
3356 ! and use it as the final cape, april 26, 1995
3357 !
3358    do n = 1,5
3359       do i = 1,ncol
3360          if (capeten(i,n) > cape(i)) then
3361             cape(i) = capeten(i,n)
3362             lel(i) = lelten(i,n)
3363          end if
3364       end do
3365    end do
3366 !
3367 ! put lower bound on cape for diagnostic purposes.
3368 !
3369    do i = 1,ncol
3370       cape(i) = max(cape(i), 0._r8)
3371    end do
3372 !
3373    return
3374 end subroutine buoyan_dilute
3375
3376 subroutine parcel_dilute (lchnk, ncol, msg, klaunch, p, t, q, tpert, tp, tpv, qstp, pl, tl, lcl)
3377
3378 ! Routine  to determine
3379 !   1. Tp   - Parcel temperature
3380 !   2. qstp - Saturated mixing ratio at the parcel temperature.
3381
3382 !--------------------
3383 implicit none
3384 !--------------------
3385
3386 integer, intent(in) :: lchnk
3387 integer, intent(in) :: ncol
3388 integer, intent(in) :: msg
3389
3390 integer, intent(in), dimension(pcols) :: klaunch(pcols)
3391
3392 real(r8), intent(in), dimension(pcols,pver) :: p
3393 real(r8), intent(in), dimension(pcols,pver) :: t
3394 real(r8), intent(in), dimension(pcols,pver) :: q
3395 real(r8), intent(in), dimension(pcols) :: tpert ! PBL temperature perturbation.
3396
3397 real(r8), intent(inout), dimension(pcols,pver) :: tp    ! Parcel temp.
3398 real(r8), intent(inout), dimension(pcols,pver) :: qstp  ! Parcel water vapour (sat value above lcl).
3399 real(r8), intent(inout), dimension(pcols) :: tl         ! Actual temp of LCL.
3400 real(r8), intent(inout), dimension(pcols) :: pl          ! Actual pressure of LCL.
3401
3402 integer, intent(inout), dimension(pcols) :: lcl ! Lifting condesation level (first model level with saturation).
3403
3404 real(r8), intent(out), dimension(pcols,pver) :: tpv   ! Define tpv within this routine.
3405
3406 !--------------------
3407
3408 ! Have to be careful as s is also dry static energy.
3409
3410
3411 ! If we are to retain the fact that CAM loops over grid-points in the internal
3412 ! loop then we need to dimension sp,atp,mp,xsh2o with ncol.
3413
3414
3415 real(r8) tmix(pcols,pver)        ! Tempertaure of the entraining parcel.
3416 real(r8) qtmix(pcols,pver)       ! Total water of the entraining parcel.
3417 real(r8) qsmix(pcols,pver)       ! Saturated mixing ratio at the tmix.
3418 real(r8) smix(pcols,pver)        ! Entropy of the entraining parcel.
3419 real(r8) xsh2o(pcols,pver)       ! Precipitate lost from parcel.
3420 real(r8) ds_xsh2o(pcols,pver)    ! Entropy change due to loss of condensate.
3421 real(r8) ds_freeze(pcols,pver)   ! Entropy change sue to freezing of precip.
3422
3423 real(r8) mp(pcols)    ! Parcel mass flux.
3424 real(r8) qtp(pcols)   ! Parcel total water.
3425 real(r8) sp(pcols)    ! Parcel entropy.
3426
3427 real(r8) sp0(pcols)    ! Parcel launch entropy.
3428 real(r8) qtp0(pcols)   ! Parcel launch total water.
3429 real(r8) mp0(pcols)    ! Parcel launch relative mass flux.
3430
3431 real(r8) lwmax      ! Maximum condesate that can be held in cloud before rainout.
3432 real(r8) dmpdp      ! Parcel fractional mass entrainment rate (/mb).
3433 !real(r8) dmpdpc     ! In cloud parcel mass entrainment rate (/mb).
3434 real(r8) dmpdz      ! Parcel fractional mass entrainment rate (/m)
3435 real(r8) dpdz,dzdp  ! Hydrstatic relation and inverse of.
3436 real(r8) senv       ! Environmental entropy at each grid point.
3437 real(r8) qtenv      ! Environmental total water "   "   ".
3438 real(r8) penv       ! Environmental total pressure "   "   ".
3439 real(r8) tenv       ! Environmental total temperature "   "   ".
3440 real(r8) new_s      ! Hold value for entropy after condensation/freezing adjustments.
3441 real(r8) new_q      ! Hold value for total water after condensation/freezing adjustments.
3442 real(r8) dp         ! Layer thickness (center to center)
3443 real(r8) tfguess    ! First guess for entropy inversion - crucial for efficiency!
3444 real(r8) tscool     ! Super cooled temperature offset (in degC) (eg -35).
3445
3446 real(r8) qxsk, qxskp1        ! LCL excess water (k, k+1)
3447 real(r8) dsdp, dqtdp, dqxsdp ! LCL s, qt, p gradients (k, k+1)
3448 real(r8) slcl,qtlcl,qslcl    ! LCL s, qt, qs values.
3449
3450 integer rcall       ! Number of ientropy call for errors recording
3451 integer nit_lheat     ! Number of iterations for condensation/freezing loop.
3452 integer i,k,ii   ! Loop counters.
3453
3454 !======================================================================
3455 !    SUMMARY
3456 !
3457 !  9/9/04 - Assumes parcel is initiated from level of maxh (klaunch)
3458 !           and entrains at each level with a specified entrainment rate.
3459 !
3460 ! 15/9/04 - Calculates lcl(i) based on k where qsmix is first < qtmix.
3461 !
3462 !======================================================================
3463 !
3464 ! Set some values that may be changed frequently.
3465 !
3466
3467 nit_lheat = 2 ! iterations for ds,dq changes from condensation freezing.
3468 dmpdz=-1.e-3_r8        ! Entrainment rate. (-ve for /m)
3469 !dmpdpc = 3.e-2_r8   ! In cloud entrainment rate (/mb).
3470 lwmax = 1.e-3_r8    ! Need to put formula in for this.
3471 tscool = 0.0_r8   ! Temp at which water loading freezes in the cloud.
3472
3473 qtmix=0._r8
3474 smix=0._r8
3475
3476 qtenv = 0._r8
3477 senv = 0._r8
3478 tenv = 0._r8
3479 penv = 0._r8
3480
3481 qtp0 = 0._r8
3482 sp0  = 0._r8
3483 mp0 = 0._r8
3484
3485 qtp = 0._r8
3486 sp = 0._r8
3487 mp = 0._r8
3488
3489 new_q = 0._r8
3490 new_s = 0._r8
3491
3492 ! **** Begin loops ****
3493
3494 do k = pver, msg+1, -1
3495    do i=1,ncol
3496
3497 ! Initialize parcel values at launch level.
3498
3499       if (k == klaunch(i)) then
3500          qtp0(i) = q(i,k)   ! Parcel launch total water (assuming subsaturated) - OK????.
3501          sp0(i)  = entropy(t(i,k),p(i,k),qtp0(i))  ! Parcel launch entropy.
3502          mp0(i)  = 1._r8       ! Parcel launch relative mass (i.e. 1 parcel stays 1 parcel for dmpdp=0, undilute).
3503          smix(i,k)  = sp0(i)
3504          qtmix(i,k) = qtp0(i)
3505          tfguess = t(i,k)
3506          rcall = 1
3507          call ientropy (rcall,i,lchnk,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess)
3508       end if
3509
3510 ! Entraining levels
3511
3512       if (k < klaunch(i)) then
3513
3514 ! Set environmental values for this level.
3515
3516          dp = (p(i,k)-p(i,k+1)) ! In -ve mb as p decreasing with height - difference between center of layers.
3517          qtenv = 0.5_r8*(q(i,k)+q(i,k+1))         ! Total water of environment.
3518          tenv  = 0.5_r8*(t(i,k)+t(i,k+1))
3519          penv  = 0.5_r8*(p(i,k)+p(i,k+1))
3520
3521          senv  = entropy(tenv,penv,qtenv)  ! Entropy of environment.
3522
3523 ! Determine fractional entrainment rate /pa given value /m.
3524
3525          dpdz = -(penv*grav)/(rgas*tenv) ! in mb/m since  p in mb.
3526          dzdp = 1._r8/dpdz                  ! in m/mb
3527          dmpdp = dmpdz*dzdp              ! /mb Fractional entrainment
3528
3529 ! Sum entrainment to current level
3530 ! entrains q,s out of intervening dp layers, in which linear variation is assumed
3531 ! so really it entrains the mean of the 2 stored values.
3532
3533          sp(i)  = sp(i)  - dmpdp*dp*senv
3534          qtp(i) = qtp(i) - dmpdp*dp*qtenv
3535          mp(i)  = mp(i)  - dmpdp*dp
3536
3537 ! Entrain s and qt to next level.
3538
3539          smix(i,k)  = (sp0(i)  +  sp(i)) / (mp0(i) + mp(i))
3540          qtmix(i,k) = (qtp0(i) + qtp(i)) / (mp0(i) + mp(i))
3541
3542 ! Invert entropy from s and q to determine T and saturation-capped q of mixture.
3543 ! t(i,k) used as a first guess so that it converges faster.
3544
3545          tfguess = tmix(i,k+1)
3546          rcall = 2
3547          call ientropy(rcall,i,lchnk,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess)
3548
3549 !
3550 ! Determine if this is lcl of this column if qsmix <= qtmix.
3551 ! FIRST LEVEL where this happens on ascending.
3552
3553          if (qsmix(i,k) <= qtmix(i,k) .and. qsmix(i,k+1) > qtmix(i,k+1)) then
3554             lcl(i) = k
3555             qxsk   = qtmix(i,k) - qsmix(i,k)
3556             qxskp1 = qtmix(i,k+1) - qsmix(i,k+1)
3557             dqxsdp = (qxsk - qxskp1)/dp
3558             pl(i)  = p(i,k+1) - qxskp1/dqxsdp    ! pressure level of actual lcl.
3559             dsdp   = (smix(i,k)  - smix(i,k+1))/dp
3560             dqtdp  = (qtmix(i,k) - qtmix(i,k+1))/dp
3561             slcl   = smix(i,k+1)  +  dsdp* (pl(i)-p(i,k+1))
3562             qtlcl  = qtmix(i,k+1) +  dqtdp*(pl(i)-p(i,k+1))
3563
3564             tfguess = tmix(i,k)
3565             rcall = 3
3566             call ientropy (rcall,i,lchnk,slcl,pl(i),qtlcl,tl(i),qslcl,tfguess)
3567
3568 !            write(iulog,*)' '
3569 !            write(iulog,*)' p',p(i,k+1),pl(i),p(i,lcl(i))
3570 !            write(iulog,*)' t',tmix(i,k+1),tl(i),tmix(i,lcl(i))
3571 !            write(iulog,*)' s',smix(i,k+1),slcl,smix(i,lcl(i))
3572 !            write(iulog,*)'qt',qtmix(i,k+1),qtlcl,qtmix(i,lcl(i))
3573 !            write(iulog,*)'qs',qsmix(i,k+1),qslcl,qsmix(i,lcl(i))
3574
3575          endif
3576 !
3577       end if !  k < klaunch
3578
3579
3580    end do ! Levels loop
3581 end do ! Columns loop
3582
3583 !!!!!!!!!!!!!!!!!!!!!!!!!!END ENTRAINMENT LOOP!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3584
3585 !! Could stop now and test with this as it will provide some estimate of buoyancy
3586 !! without the effects of freezing/condensation taken into account for tmix.
3587
3588 !! So we now have a profile of entropy and total water of the entraining parcel
3589 !! Varying with height from the launch level klaunch parcel=environment. To the
3590 !! top allowed level for the existence of convection.
3591
3592 !! Now we have to adjust these values such that the water held in vaopor is < or
3593 !! = to qsmix. Therefore, we assume that the cloud holds a certain amount of
3594 !! condensate (lwmax) and the rest is rained out (xsh2o). This, obviously
3595 !! provides latent heating to the mixed parcel and so this has to be added back
3596 !! to it. But does this also increase qsmix as well? Also freezing processes
3597
3598
3599 xsh2o = 0._r8
3600 ds_xsh2o = 0._r8
3601 ds_freeze = 0._r8
3602
3603 !!!!!!!!!!!!!!!!!!!!!!!!!PRECIPITATION/FREEZING LOOP!!!!!!!!!!!!!!!!!!!!!!!!!!
3604 !! Iterate solution twice for accuracy
3605
3606
3607
3608 do k = pver, msg+1, -1
3609    do i=1,ncol
3610
3611 ! Initialize variables at k=klaunch
3612
3613       if (k == klaunch(i)) then
3614
3615 ! Set parcel values at launch level assume no liquid water.
3616
3617          tp(i,k)    = tmix(i,k)
3618          qstp(i,k)  = q(i,k)
3619          tpv(i,k)   =  (tp(i,k) + tpert(i)) * (1._r8+1.608_r8*qstp(i,k)) / (1._r8+qstp(i,k))
3620
3621       end if
3622
3623       if (k < klaunch(i)) then
3624
3625 ! Initiaite loop if switch(2) = .T. - RBN:DILUTE - TAKEN OUT BUT COULD BE RETURNED LATER.
3626
3627 ! Iterate nit_lheat times for s,qt changes.
3628
3629          do ii=0,nit_lheat-1
3630
3631 ! Rain (xsh2o) is excess condensate, bar LWMAX (Accumulated loss from qtmix).
3632
3633             xsh2o(i,k) = max (0._r8, qtmix(i,k) - qsmix(i,k) - lwmax)
3634
3635 ! Contribution to ds from precip loss of condensate (Accumulated change from smix).(-ve)
3636
3637             ds_xsh2o(i,k) = ds_xsh2o(i,k+1) - cpliq * log (tmix(i,k)/tfreez) * max(0._r8,(xsh2o(i,k)-xsh2o(i,k+1)))
3638 !
3639 ! Entropy of freezing: latice times amount of water involved divided by T.
3640 !
3641
3642             if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) == 0._r8) then ! One off freezing of condensate.
3643                ds_freeze(i,k) = (latice/tmix(i,k)) * max(0._r8,qtmix(i,k)-qsmix(i,k)-xsh2o(i,k)) ! Gain of LH
3644             end if
3645
3646             if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) /= 0._r8) then ! Continual freezing of additional condensate.
3647                ds_freeze(i,k) = ds_freeze(i,k+1)+(latice/tmix(i,k)) * max(0._r8,(qsmix(i,k+1)-qsmix(i,k)))
3648             end if
3649
3650 ! Adjust entropy and accordingly to sum of ds (be careful of signs).
3651
3652             new_s = smix(i,k) + ds_xsh2o(i,k) + ds_freeze(i,k)
3653
3654 ! Adjust liquid water and accordingly to xsh2o.
3655
3656             new_q = qtmix(i,k) - xsh2o(i,k)
3657
3658 ! Invert entropy to get updated Tmix and qsmix of parcel.
3659
3660             tfguess = tmix(i,k)
3661             rcall =4
3662             call ientropy (rcall,i,lchnk,new_s, p(i,k), new_q, tmix(i,k), qsmix(i,k), tfguess)
3663
3664          end do  ! Iteration loop for freezing processes.
3665
3666 ! tp  - Parcel temp is temp of mixture.
3667 ! tpv - Parcel v. temp should be density temp with new_q total water.
3668
3669          tp(i,k)    = tmix(i,k)
3670
3671 ! tpv = tprho in the presence of condensate (i.e. when new_q > qsmix)
3672
3673          if (new_q > qsmix(i,k)) then  ! Super-saturated so condensate present - reduces buoyancy.
3674             qstp(i,k) = qsmix(i,k)
3675          else                          ! Just saturated/sub-saturated - no condensate virtual effects.
3676             qstp(i,k) = new_q
3677          end if
3678
3679          tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k)) / (1._r8+ new_q)
3680
3681       end if ! k < klaunch
3682
3683    end do ! Loop for columns
3684
3685 end do  ! Loop for vertical levels.
3686
3687
3688 return
3689 end subroutine parcel_dilute
3690
3691 !-----------------------------------------------------------------------------------------
3692 real(r8) function entropy(TK,p,qtot)
3693 !-----------------------------------------------------------------------------------------
3694 !
3695 ! TK(K),p(mb),qtot(kg/kg)
3696 ! from Raymond and Blyth 1992
3697 !
3698      real(r8), intent(in) :: p,qtot,TK
3699      real(r8) :: qv,qsat,e,esat,L,eref,pref
3700
3701 pref = 1000.0_r8           ! mb
3702 eref = 6.106_r8            ! sat p at tfreez (mb)
3703
3704 L = rl - (cpliq - cpwv)*(TK-tfreez)         ! T IN CENTIGRADE
3705
3706 ! Replace call to satmixutils.
3707
3708 esat = c1*exp(c2*(TK-tfreez)/(c3+TK-tfreez))       ! esat(T) in mb
3709 qsat=eps1*esat/(p-esat)                      ! Sat. mixing ratio (in kg/kg).
3710
3711 qv = min(qtot,qsat)                         ! Partition qtot into vapor part only.
3712 e = qv*p / (eps1 +qv)
3713
3714 entropy = (cpres + qtot*cpliq)*log( TK/tfreez) - rgas*log( (p-e)/pref ) + &
3715         L*qv/TK - qv*rh2o*log(qv/qsat)
3716 !
3717 return
3718 end FUNCTION entropy
3719
3720 !
3721 !-----------------------------------------------------------------------------------------
3722    SUBROUTINE ientropy (rcall,icol,lchnk,s,p,qt,T,qsat,Tfg)
3723 !-----------------------------------------------------------------------------------------
3724 !
3725 ! p(mb), Tfg/T(K), qt/qv(kg/kg), s(J/kg).
3726 ! Inverts entropy, pressure and total water qt
3727 ! for T and saturated vapor mixing ratio
3728 !
3729 #ifndef WRF_PORT
3730      use phys_grid, only: get_rlon_p, get_rlat_p
3731 #endif
3732
3733      integer, intent(in) :: icol, lchnk, rcall
3734      real(r8), intent(in)  :: s, p, Tfg, qt
3735      real(r8), intent(out) :: qsat, T
3736      real(r8) :: qv,Ts,dTs,fs1,fs2,esat
3737      real(r8) :: pref,eref,L,e
3738      real(r8) :: this_lat,this_lon
3739      integer :: LOOPMAX,i
3740
3741 LOOPMAX = 100                   !* max number of iteration loops
3742
3743 ! Values for entropy
3744 pref = 1000.0_r8           ! mb ref pressure.
3745 eref = 6.106_r8           ! sat p at tfreez (mb)
3746
3747 ! Invert the entropy equation -- use Newton's method
3748
3749 Ts = Tfg                  ! Better first guess based on Tprofile from conv.
3750
3751 converge: do i=0, LOOPMAX
3752
3753    L = rl - (cpliq - cpwv)*(Ts-tfreez)
3754
3755    esat = c1*exp(c2*(Ts-tfreez)/(c3+Ts-tfreez)) ! Bolton (eq. 10)
3756    qsat = eps1*esat/(p-esat)
3757    qv = min(qt,qsat)
3758    e = qv*p / (eps1 +qv)  ! Bolton (eq. 16)
3759    fs1 = (cpres + qt*cpliq)*log( Ts/tfreez ) - rgas*log( (p-e)/pref ) + &
3760         L*qv/Ts - qv*rh2o*log(qv/qsat) - s
3761
3762    L = rl - (cpliq - cpwv)*(Ts-1._r8-tfreez)
3763
3764    esat = c1*exp(c2*(Ts-1._r8-tfreez)/(c3+Ts-1._r8-tfreez))
3765    qsat = eps1*esat/(p-esat)
3766    qv = min(qt,qsat)
3767    e = qv*p / (eps1 +qv)
3768    fs2 = (cpres + qt*cpliq)*log( (Ts-1._r8)/tfreez ) - rgas*log( (p-e)/pref ) + &
3769         L*qv/(Ts-1._r8) - qv*rh2o*log(qv/qsat) - s
3770
3771    dTs = fs1/(fs2 - fs1)
3772    Ts  = Ts+dTs
3773    if (abs(dTs).lt.0.001_r8) exit converge
3774    if (i .eq. LOOPMAX - 1) then
3775 #ifndef WRF_PORT
3776       this_lat = get_rlat_p(lchnk, icol)*57.296_r8
3777       this_lon = get_rlon_p(lchnk, icol)*57.296_r8
3778 #else
3779 !Do not worry about the specific lat/lon in WRF
3780       this_lat = 0.
3781       this_lon = 0.
3782 #endif
3783       write(iulog,*) '*** ZM_CONV: IENTROPY: Failed and about to exit, info follows ****'
3784 #ifdef WRF_PORT
3785       call wrf_debug(1,iulog)
3786 #endif
3787       write(iulog,100) 'ZM_CONV: IENTROPY. Details: call#,lchnk,icol= ',rcall,lchnk,icol, &
3788        ' lat: ',this_lat,' lon: ',this_lon, &
3789        ' P(mb)= ', p, ' Tfg(K)= ', Tfg, ' qt(g/kg) = ', 1000._r8*qt, &
3790        ' qsat(g/kg) = ', 1000._r8*qsat,', s(J/kg) = ',s
3791 #ifdef WRF_PORT
3792       call wrf_debug(1,iulog)
3793 #endif
3794       write(iulog,*) '*** Please report this crash to Po-Lun.Ma@pnnl.gov ***'
3795 #ifdef WRF_PORT
3796       call wrf_debug(1,iulog)
3797 #endif
3798       call endrun('**** ZM_CONV IENTROPY: Tmix did not converge ****')
3799    end if
3800 enddo converge
3801
3802 ! Replace call to satmixutils.
3803
3804 esat = c1*exp(c2*(Ts-tfreez)/(c3+Ts-tfreez))
3805 qsat=eps1*esat/(p-esat)
3806
3807 qv = min(qt,qsat)                             !       /* check for saturation */
3808 T = Ts
3809
3810  100    format (A,I1,I4,I4,7(A,F6.2))
3811
3812 return
3813 end SUBROUTINE ientropy
3814
3815 end module module_cu_camzm