phys/module_cu_camzm_driver.F

   1 MODULE module_cu_camzm_driver
   2
   3   !Roughly based on zm_conv_intr.F90 from CAM
   4
   5   !-------------------------------------------
   6   !Future Modifications and important warning (BSINGH:02/01/2013: Notes from WIG):
   7   !===========================================
   8   !1. Fracis is hardwired to 1 for first 3 constituents (vapor,liq mass,ice mass)
   9   !   for WRF_CHEM simulations and first 5 constituents(liq # and ice #) for WRF
  10   !   physcis ONLY runs. Should be generalized for other constituents (like snow,graupel etc.)
  11   !
  12   !2. In CAM model, ZM scheme only acts up to 40 mb. In WRF, it will function all the way
  13   !   to the model top no matter how high it is. We could do a better job at capping this
  14   !   by checking to see where 40 mb roughly occurs during initialization. It won't change
  15   !   too dramatically over time. Also, WRF often uses a top of 50 mb, so we rarely hit 40 mb.
  16   !   This is a minor issue.
  17   !3. ZM is set to never do convection w/in the PBL. In CAM model, this is dependent on the
  18   !   use of UW shallow. Should this be a conditional setting in WRF or hard-wired as it
  19   !   currently is?
  20   !-----------------------------------------
  21
  22
  23   USE module_cam_support, only: pcnst =>pcnst_runtime, pcols, pver, pverp
  24 #if ( WRF_CHEM == 1 )
  25   USE module_cam_support, only: cam_mam_aerosols
  26 #endif
  27   USE shr_kind_mod,    only: r8 => shr_kind_r8
  28   USE module_cu_camzm, only: convtran, momtran, zm_conv_evap, zm_convr
  29
  30   IMPLICIT NONE
  31
  32   PRIVATE                  !Default to private
  33   PUBLIC :: &              !Public entities
  34        camzm_driver  , &
  35 #if ( WRF_CHEM == 1 )
  36        zm_conv_tend_2, &
  37 #endif
  38        zm_conv_init
  39
  40   !Module level variables which are required for zm_conv_tend_2 subrouine
  41   integer :: ixcldliq, ixcldice, ixnumliq, ixnumice
  42
  43 CONTAINS
  44
  45 !------------------------------------------------------------------------
  46 SUBROUTINE camzm_driver(                                      &
  47               ids,ide, jds,jde, kds,kde                       &
  48              ,ims,ime, jms,jme, kms,kme                       &
  49              ,its,ite, jts,jte, kts,kte                       &
  50              ,itimestep, bl_pbl_physics, sf_sfclay_physics    &
  51              ,th, t_phy, tsk, tke_pbl, ust, qv, qc, qi        &
  52              ,mavail, kpbl, pblh, xland, z                    &
  53              ,z_at_w, dz8w, ht, p, p8w, pi_phy, psfc          &
  54              ,u_phy, v_phy, hfx, qfx, cldfra, cldfra_mp_all   &
  55              ,is_CAMMGMP_used, tpert_camuwpbl                 &
  56              ,dx, dt, stepcu, cudt, curr_secs, adapt_step_flag&
  57              ,cudtacttime                                     &
  58              ,cape_out, mu_out, md_out, zmdt, zmdq            &
  59              ,rliq_out, dlf_out                               &
  60              ,pconvb, pconvt, cubot, cutop, raincv, pratec    &
  61              ,rucuten, rvcuten                                &
  62              ,rthcuten, rqvcuten, rqccuten, rqicuten          &
  63              ,rqcncuten, rqincuten                            &
  64              ,evaptzm, fzsntzm, evsntzm, evapqzm, zmflxprc    &
  65              ,zmflxsnw, zmntprpd, zmntsnpd, zmeiheat          &
  66              ,cmfmc, cmfmcdzm, preccdzm, precz                &
  67              ,zmmtu, zmmtv, zmupgu, zmupgd, zmvpgu, zmvpgd    &
  68              ,zmicuu, zmicud, zmicvu, zmicvd                  &
  69              ,zmdice, zmdliq                                  &
  70              !Following intent-OUT variables are required by wet scavenging code of WRF_CHEM package- Balwinder.Singh@pnnl.gov
  71              ,evapcdp3d, icwmrdp3d, rprddp3d                  &
  72              !Following intent-OUT variables are required by zm_conv_tend_2 call
  73              ,dp3d, du3d, ed3d, eu3d, md3d, mu3d, dsubcld2d   &
  74              ,ideep2d, jt2d, maxg2d, lengath2d                )
  75
  76 ! This routine is based on zm_conv_tend in CAM. It handles the mapping of
  77 ! variables from the WRF to the CAM framework for the Zhang-McFarlane
  78 ! convective parameterization.
  79 !
  80 ! Author: William.Gustafson@pnl.gov, Nov. 2009
  81 ! Last modified: William.Gustafson@pnl.gov, Nov. 2010
  82 !------------------------------------------------------------------------
  83   USE physconst,     only: latice,cpair, gravit
  84   USE wv_saturation, only: fqsatd,epsqs, polysvp
  85   USE wv_saturation,            only: fqsatd,epsqs, polysvp
  86
  87 ! Subroutine arguments...
  88   logical, intent(in)    ::    is_CAMMGMP_used
  89   INTEGER, INTENT(IN   ) ::    ids,ide, jds,jde, kds,kde,  &
  90                                ims,ime, jms,jme, kms,kme,  &
  91                                its,ite, jts,jte, kts,kte,  &
  92                      bl_pbl_physics, & !pbl scheme
  93                           itimestep, & !time step index
  94                   sf_sfclay_physics, & !surface layer scheme
  95                              stepcu    !number of time steps between Cu calls
  96
  97   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: &
  98                              cldfra, & !cloud fraction
  99                       cldfra_mp_all, & !cloud fraction from CAMMGMP microphysics
 100                                dz8w, & !height between interfaces (m)
 101                                   p, & !pressure at mid-level (Pa)
 102                                 p8w, & !pressure at level interface (Pa)
 103                              pi_phy, & !exner function, (p/p0)^(R/cpair) (none)
 104                                  qv, & !water vapor mixing ratio (kg/kg-dry air)
 105                                  th, & !potential temperature (K)
 106                             tke_pbl, & !turbulent kinetic energy from PBL (m2/s2)
 107                               t_phy, & !temperature (K)
 108                               u_phy, & !zonal wind component on T points (m/s)
 109                               v_phy, & !meridional wind component on T points (m/s)
 110                                   z, & !height above sea level at mid-level (m)
 111                              z_at_w    !height above sea level at interface (m)
 112
 113   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN), OPTIONAL :: &
 114                                  qc, & !cloud droplet mixing ratio (kg/kg-dry air)
 115                                  qi    !cloud ice crystal mixing ratio (kg/kg-dry air)
 116
 117   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: &
 118                             dlf_out, & !detraining cloud water tendendcy
 119                             evaptzm, & !temp. tendency - evap/snow prod from ZM (K/s)
 120                             fzsntzm, & !temp. tendency - rain to snow conversion from ZM (K/s)
 121                             evsntzm, & !temp. tendency - snow to rain conversion from ZM (K/s)
 122                             evapqzm, & !spec. humidity tend. - evaporation from ZM (kg/kg/s)
 123                            zmflxprc, & !flux of precipitation from ZM (kg/m2/s)
 124                            zmflxsnw, & !flux of snow from ZM (kg/m2/s)
 125                            zmntprpd, & !net precipitation production from ZM (kg/kg/s)
 126                            zmntsnpd, & !net snow production from ZM (kg/kg/s)
 127                            zmeiheat, & !heating by ice and evaporation from ZM (W/kg)
 128                               cmfmc, & !convective mass flux--m sub c, deep here but ultimately deep+shallow (kg/m2/s)
 129                            cmfmcdzm, & !convection mass flux from ZM deep (kg/m2/s)
 130                              md_out, & !output convective downdraft mass flux (kg/m2/s)
 131                              mu_out, & !output convective updraft mass flux (kg/m2/s)
 132                             rucuten, & !UNcoupled zonal wind tendency due to Cu scheme (m/s2)
 133                             rvcuten, & !UNcoupled meridional wind tendency due to Cu scheme (m/s2)
 134                            rthcuten, & !UNcoupled potential temperature tendendcy due to cu scheme (K/s)
 135                            rqvcuten, & !UNcoupled water vapor mixing ratio tendency due to Cu scheme (kg/kg/s)
 136                                zmdt, & !temp. tendency from moist convection (K/s)
 137                                zmdq, & !spec. humidity tendency from moist convection (kg/kg/s)
 138                               zmmtu, & !U tendency from ZM convective momentum transport (m/s2)
 139                               zmmtv, & !V tendency from ZM convective momentum transport (m/s2)
 140                              zmupgu, & !zonal force from ZM updraft pressure gradient term (m/s2)
 141                              zmupgd, & !zonal force from ZM downdraft pressure gradient term (m/s2)
 142                              zmvpgu, & !meridional force from ZM updraft pressure gradient term (m/s2)
 143                              zmvpgd, & !meridional force from ZM downdraft pressure gradient term (m/s2)
 144                              zmicuu, & !ZM in-cloud U updrafts (m/s)
 145                              zmicud, & !ZM in-cloud U downdrafts (m/s)
 146                              zmicvu, & !ZM in-cloud V updrafts (m/s)
 147                              zmicvd, & !ZM in-cloud V downdrafts (m/s)
 148                              zmdice, & !ZM cloud ice tendency (kg/kg/s)
 149                              zmdliq    !ZM cloud liquid tendency (kg/kg/s)
 150
 151   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT), OPTIONAL :: &
 152                            rqccuten, & !UNcoupled cloud droplet mixing ratio tendency due to Cu scheme (kg/kg/s)
 153                            rqicuten, & !UNcoupled ice crystal mixing ratio tendency due to Cu scheme (kg/kg/s)
 154                           rqcncuten, & !PMA
 155                           rqincuten    !PMA
 156
 157
 158   !variables required by wet scavenging in WRF_CHEM -Balwinder.Singh@pnnl.gov
 159   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT) :: &
 160                           evapcdp3d, & !Evaporation of deep convective precipitation (kg/kg/s)
 161                           icwmrdp3d, & !Deep Convection in-cloud water mixing ratio (kg/m2)
 162                            rprddp3d    !dq/dt due to deep convective rainout (kg/kg/s)
 163
 164   !variables required by zm_conv_tend_2 call
 165   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT) :: &
 166                                dp3d, &
 167                                du3d, &
 168                                ed3d, &
 169                                eu3d, &
 170                                md3d, &
 171                                mu3d
 172
 173
 174   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
 175                                 hfx, & !upward heat flux at sfc (W/m2)
 176                                  ht, & !terrain height (m)
 177                               xland, & !land/water mask, 1=land, 2=water
 178                              mavail, & !soil moisture availability
 179                                pblh, & !planetary boundary layer height (m)
 180                                psfc, & !surface pressure (Pa)
 181                                 qfx, & !upward moisture flux at sfc (kg/m2/s)
 182                      tpert_camuwpbl, & !temperature perturbation from UW CAM PBL
 183                                 tsk, & !skin temperature (K)
 184                                 ust    !u* in similarity theory (m/s)
 185
 186   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
 187                            cape_out, & !convective available potential energy (J/kg)
 188                               cubot, & !level number of base of convection
 189                               cutop, & !level number of top of convection
 190                              pconvb, & !pressure of base of convection (Pa)
 191                              pconvt, & !pressure of top of convection (Pa)
 192                              pratec, & !rain rate returned to WRF for accumulation (mm/s)
 193                            preccdzm, & !convection precipitation rate from ZM deep (m/s)
 194                               precz, & !total precipitation rate from ZM (m/s)
 195                              raincv, & !time-step convective rain amount (mm)
 196                            rliq_out    !vertical integral of reserved cloud water
 197   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: &
 198                            dsubcld2d
 199
 200   REAL, INTENT(IN) :: &
 201                                cudt, & !cumulus time step (min)
 202                           curr_secs, & !current forecast time (s)
 203                                  dt, & !domain time step (s)
 204                                  dx    !grid spacing (m)
 205
 206   REAL, INTENT (INOUT) :: &
 207                         cudtacttime    !for adaptive time stepping, next to to run scheme (s)
 208
 209   INTEGER, DIMENSION( ims:ime, jms:jme), INTENT(IN) :: &
 210                                kpbl    !index of PBL level
 211
 212   INTEGER, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: &
 213                            ideep2d, &
 214                            jt2d,    &
 215                            maxg2d,  &
 216                         lengath2d
 217
 218   LOGICAL, INTENT(IN), OPTIONAL :: &
 219                     adapt_step_flag    !using adaptive time steps?
 220
 221 ! Local variables...
 222   !Variables dimensioned for input to ZM routines
 223   REAL(r8), DIMENSION(pcols, kte+1) ::  &
 224                                mcon, & !convective mass flux--m sub c (sub arg out in CAM)
 225                                pflx, & !scattered precip flux at each level (sub arg out in CAM)
 226                               pint8, & !pressure at layer interface (Pa)
 227                                 zi8, & !height above sea level at mid-level (m)
 228                             flxprec, & !evap outfld: convective-scale flux of precip at interfaces (kg/m2/s)
 229                             flxsnow    !evap outfld: convective-scale flux of snow at interfaces (kg/m2/s)
 230
 231
 232   REAL(r8), DIMENSION(pcols, kte, pcnst) ::  &
 233                                 qh8    !specific humidity (kg/kg-moist air)
 234
 235   REAL(r8), DIMENSION(pcols, kte, 3) ::  &
 236                               cloud, & !holder for cloud water and ice (q in CAM)
 237                            cloudtnd    !holder for cloud tendencies (ptend_loc%q in CAM)
 238
 239 !BSINGH - 02/18/2013: FRACIS is used for qv,qc and qi for WRF_CHEM simulation but for
 240 ! WRF simulations it is used for qnc and qni also
 241 #if ( WRF_CHEM == 1 )
 242  REAL(r8), DIMENSION(pcols, kte, 3) ::  &
 243 #else
 244  REAL(r8), DIMENSION(pcols, kte, 5) ::  &
 245 #endif
 246                              fracis    !fraction of cloud species that are insoluble
 247
 248
 249   REAL(r8), DIMENSION(pcols, kte, 2) ::  &
 250                                icwu, & !in-cloud winds in upraft (m/s)
 251                                icwd, & !in-cloud winds in downdraft (m/s)
 252                              pguall, & !apparent force from updraft pres. gradient (~units?)
 253                              pgdall, & !apparent force from downdraft pres. gradient (~units?)
 254                               winds, & !wind components (m/s)
 255                          wind_tends    !wind component tendencies (m/s2)
 256
 257   REAL(r8), DIMENSION(pcols, kte) ::  &
 258                                cld8, & !cloud fraction
 259                                 cme, & !cmf condensation - evaporation (~units?, sub arg out in CAM)
 260                                 dlf, & !scattered version of the detraining cld h2o tendency (~units?)
 261                          fake_dpdry, & !place holder for dpdry, delta pres. of dry atmos.
 262                             ntprprd, & !evap outfld: net precip production in layer (kg/kg/s)
 263                             ntsnprd, & !evap outfld: net snow production in layer (kg/kg/s)
 264                               pdel8, & !pressure thickness of layer (between interfaces, Pa)
 265                               pmid8, & !pressure at layer middle (Pa)
 266                                 ql8, & !in cloud liquid water (~units?)
 267                              ql8prm, & !cloud liquid water (~units?)
 268                                 qi8, & !cloud ice (~units?)
 269                                  t8, & !temperature (K)
 270                                 zm8, & !height above ground at mid-level (m)
 271                                qtnd, & !specific humidity tendency (kg/kg/s)
 272                                rprd, & !rain production rate (kg/kg/s, comes from pbuf array in CAM, add to restart?~)
 273                                stnd, & !heating rate (dry static energy tendency, W/kg)
 274                       tend_s_snwprd, & !heating rate of snow production (~units?)
 275                     tend_s_snwevmlt, & !heating rate of evap/melting of snow (~units?)
 276                                 zdu, & !detraining mass flux (~units? sub arg out in CAM)
 277                             evapcdp    !Evaporation of deep convective precipitation (kg/kg/s)
 278
 279   REAL(r8), DIMENSION(pcols, kte) ::  &
 280                                 esl, & !saturation vapor pressure
 281                                 qvs, & !saturation specific humidity
 282                           qcten_det, &
 283                           qiten_det, &
 284                          qcnten_det, &
 285                          qinten_det, &
 286                           qsten_det
 287   REAL(r8), DIMENSION(pcols) ::  &
 288                                cape, & !convective available potential energy (J/kg)
 289                               jctop, & !row of top-of-deep-convection indices passed out (sub arg out in CAM)
 290                               jcbot, & !row of base of cloud indices passed out (sub arg out in CAM)
 291                            landfrac, & !land fraction
 292                               pblh8, & !planetary boundary layer height (m)
 293                                phis, & !geopotential at terrain height (m2/s2)
 294                                prec, & !convective-scale precipitation rate from ZM (m/s)
 295                                rliq, & !reserved liquid (not yet in cldliq) for energy integrals (units? sub arg out in CAM)
 296                                snow, & !convective-scale snowfall rate from ZM (m/s)
 297                               tpert    !thermal (convective) temperature excess (K)
 298
 299   REAL(r8), DIMENSION(pcols, kte, (ime-ims+1)*(jme-jms+1)) :: &
 300                                  dp, & !layer pres. thickness between interfaces (mb)
 301                                  du, &
 302                                  ed, &
 303                                  eu, &
 304                                  md, &
 305                                  mu
 306
 307   REAL(r8), DIMENSION(pcols, (ime-ims+1)*(jme-jms+1)) :: &
 308                             dsubcld    ! layer pres. thickness between LCL and maxi (mb)
 309
 310   INTEGER, DIMENSION(pcols, (ime-ims+1)*(jme-jms+1)) :: &
 311                               ideep, & ! holds position of gathered points
 312                                  jt, & ! top-level index of deep cumulus convection
 313                                maxg    ! gathered values of maxi
 314
 315   INTEGER, DIMENSION((ime-ims+1)*(jme-jms+1)) :: &
 316                             lengath    ! number of gathered points
 317
 318   !Other local vars
 319   REAL(r8):: dum1,cudts,hcudts
 320   INTEGER :: i, j, k, kflip, n, ncnst
 321   INTEGER :: lchnk         !chunk identifier, used to map 2-D to 1-D arrays in WRF
 322   INTEGER :: ncol          !number of atmospheric columns in chunk
 323   LOGICAL, DIMENSION(3) :: l_qt    !logical switches for constituent tendency presence
 324   LOGICAL, DIMENSION(2) :: l_windt !logical switches for wind tendency presence
 325   LOGICAL :: run_param , & !flag for handling alternate cumulus call freq.
 326              doing_adapt_dt , decided
 327 !
 328 ! Check to see if this is a convection timestep...
 329 !
 330
 331 !  Initialization for adaptive time step.
 332
 333    if (cudt .eq. 0.) then
 334     cudts = dt
 335     hcudts=cudts*0.5_r8
 336    else
 337     cudts=cudt*60._r8
 338     hcudts=cudts*0.5_r8
 339    end if
 340
 341    doing_adapt_dt = .FALSE.
 342    IF ( PRESENT(adapt_step_flag) ) THEN
 343       IF ( adapt_step_flag ) THEN
 344          doing_adapt_dt = .TRUE.
 345          IF ( cudtacttime .EQ. 0. ) THEN
 346             cudtacttime = curr_secs + cudt*60.
 347          END IF
 348       END IF
 349    END IF
 350
 351 !  Do we run through this scheme or not?
 352
 353 !    Test 1:  If this is the initial model time, then yes.
 354 !                ITIMESTEP=1
 355 !    Test 2:  If the user asked for the cumulus to be run every time step, then yes.
 356 !                CUDT=0 or STEPCU=1
 357 !    Test 3:  If not adaptive dt, and this is on the requested cumulus frequency, then yes.
 358 !                MOD(ITIMESTEP,STEPCU)=0
 359 !    Test 4:  If using adaptive dt and the current time is past the last requested activate cumulus time, then yes.
 360 !                CURR_SECS >= CUDTACTTIME
 361
 362 !  If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
 363 !  to TRUE.  The decided flag says that one of these tests was able to say "yes", run the scheme.
 364 !  We only proceed to other tests if the previous tests all have left decided as FALSE.
 365
 366 !  If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
 367 !  cumulus run.
 368
 369    decided = .FALSE.
 370    run_param = .FALSE.
 371    IF ( ( .NOT. decided ) .AND. &
 372         ( itimestep .EQ. 1 ) ) THEN
 373       run_param   = .TRUE.
 374       decided     = .TRUE.
 375    END IF
 376
 377    IF ( ( .NOT. decided ) .AND. &
 378         ( ( cudt .EQ. 0. ) .OR. ( stepcu .EQ. 1 ) ) ) THEN
 379       run_param   = .TRUE.
 380       decided     = .TRUE.
 381    END IF
 382
 383    IF ( ( .NOT. decided ) .AND. &
 384         ( .NOT. doing_adapt_dt ) .AND. &
 385         ( MOD(itimestep,stepcu) .EQ. 0 ) ) THEN
 386       run_param   = .TRUE.
 387       decided     = .TRUE.
 388    END IF
 389
 390    IF ( ( .NOT. decided ) .AND. &
 391         ( doing_adapt_dt ) .AND. &
 392         ( curr_secs .GE. cudtacttime ) ) THEN
 393       run_param   = .TRUE.
 394       decided     = .TRUE.
 395       cudtacttime = curr_secs + cudt*60
 396    END IF
 397
 398   !Leave the subroutine if it is not yet time to call the cumulus param
 399   if( .not. run_param ) return
 400 !
 401 ! Initialize...
 402 !
 403   ncol  = 1          !chunk size in WRF is 1 since we loop over all columns in a tile
 404
 405   cape_out(its:ite, jts:jte)        = 0.
 406   mu_out(its:ite, kts:kte, jts:jte) = 0.
 407   md_out(its:ite, kts:kte, jts:jte) = 0.
 408   zmdt(its:ite, kts:kte, jts:jte)   = 0.
 409 !
 410 ! Map variables to inputs for zm_convr and call it...
 411 ! Loop over the points in the tile and treat them each as a CAM chunk.
 412 !
 413   do j = jts,jte
 414      do i = its,ite
 415         lchnk = (j-jts)*(ite-its+1) + (i-its+1) !1-D index location from the 2-D tile
 416
 417         !Flip variables on the layer middles
 418         do k = kts,kte
 419            kflip = kte-k+1
 420            if(is_CAMMGMP_used) then
 421               cld8(1,kflip)  = cldfra_mp_all(i,k,j)
 422            else
 423               cld8(1,kflip)  = cldfra(i,k,j)
 424            endif
 425                    if (itimestep .eq. 1) cld8(1,kflip) =0.
 426            pdel8(1,kflip) = p8w(i,k,j) - p8w(i,k+1,j)
 427            pmid8(1,kflip) = p(i,k,j)
 428            qh8(1,kflip,1) = max( qv(i,k,j)/(1.+qv(i,k,j)), 1e-30 ) !values of 0 cause a crash in entropy
 429            if( present(qc) ) then
 430               ql8prm(1,kflip) = qc(i,k,j)/(1.+qv(i,k,j)) !Convert to moist mix ratio
 431            else
 432               ql8prm(1,kflip) = 0.
 433            end if
 434            if( present(qi) ) then
 435               qi8(1,kflip) = qi(i,k,j)/(1.+qv(i,k,j)) !Used in convtran, ditto for conversion
 436            else
 437               qi8(1,kflip) = 0.
 438            end if
 439            t8(1,kflip)    = t_phy(i,k,j)
 440            zm8(1,kflip)   = z(i,k,j) - ht(i,j) !Height above the ground at midlevels
 441
 442            dp(1,kflip,lchnk) = 0.0_r8
 443            du(1,kflip,lchnk) = 0.0_r8
 444            ed(1,kflip,lchnk) = 0.0_r8
 445            eu(1,kflip,lchnk) = 0.0_r8
 446            md(1,kflip,lchnk) = 0.0_r8
 447            mu(1,kflip,lchnk) = 0.0_r8
 448         end do
 449
 450         !Flip variables on the layer interfaces
 451         do k = kts,kte+1
 452            kflip = kte-k+2
 453
 454            pint8(1,kflip) = p8w(i,k,j)
 455            zi8(1,kflip)   = z_at_w(i,k,j) -ht(i,j) !Height above the ground at interfaces
 456         end do
 457
 458         !Other necessary conversions for input to ZM
 459         if( xland(i,j)==2 ) then
 460            landfrac(1) = 1. !land, WRF is all or nothing
 461         else
 462            landfrac(1) = 0. !water
 463         end if
 464         pblh8(1) = pblh(i,j)
 465         phis(1)  = ht(i,j)*gravit
 466
 467         call get_tpert(bl_pbl_physics, sf_sfclay_physics, dx, &
 468              mavail(i,j), kpbl(i,j), pblh(i,j), &
 469              dz8w(i,1,j), psfc(i,j), qv(i,1,j), t_phy(i,1,j), &
 470              th(i,1,j), tsk(i,j), tke_pbl(i,:,j), ust(i,j),   &
 471              u_phy(i,1,j), v_phy(i,1,j), hfx(i,j), qfx(i,j), &
 472              tpert_camuwpbl(i,j), kte, &
 473              tpert(1))
 474
 475         !Call the Zhang-McFarlane (1996) convection parameterization
 476         !PMA: pass in 0.5 x cudt (sec)
 477         !Modified by Balwinder.Singh@pnnl.gov: We are no longer using 0.5*dt
 478
 479         !PMA: pass in 0.5 x cudt (sec)
 480         call zm_convr( lchnk, ncol, &
 481              t8, qh8, prec, jctop, jcbot, &
 482              pblh8, zm8, phis, zi8, qtnd, &
 483              stnd, pmid8, pint8, pdel8, &
 484              hcudts, mcon, cme, cape, &
 485              tpert, dlf, pflx, zdu, rprd, &
 486              mu(:,:,lchnk), md(:,:,lchnk),du(:,:,lchnk),eu(:,:,lchnk),ed(:,:,lchnk), &
 487              dp(:,:,lchnk), dsubcld(:,lchnk), jt(:,lchnk), maxg(:,lchnk), ideep(:,lchnk), &
 488              lengath(lchnk), ql8, rliq, landfrac )
 489
 490         !Start mapping CAM output to WRF output variables. We follow the
 491         !same order as in CAM's zm_conv_tend to ease maintenance...
 492         do k=kts,kte
 493            kflip = kte-k+1
 494            dlf(1,kflip)    = max(min(1.e-6_r8,dlf(1,kflip)),0._r8)
 495            dlf_out(i,k,j)  = dlf(1,kflip)
 496         end do
 497         cape_out(i,j) = cape(1)
 498         rliq_out(i,j) = rliq(1)
 499
 500         !Convert mass flux from reported mb/s to kg/m2/s
 501         mcon(:ncol,:pverp) = mcon(:ncol,:pverp) * 100._r8/gravit
 502
 503         ! Store upward and downward mass fluxes in un-gathered arrays
 504         ! + convert from mb/s to kg/m2/s
 505         do n=1,lengath(lchnk)
 506            do k=kts,kte
 507               kflip = kte-k+1
 508 !              ii = ideep(n,lchnk) <--in WRF this is always 1 because chunk size=1
 509               mu_out(i,k,j) = mu(n,kflip,lchnk) * 100._r8/gravit
 510               md_out(i,k,j) = md(n,kflip,lchnk) * 100._r8/gravit
 511            end do
 512         end do
 513
 514         do k=kts,kte
 515            kflip = kte-k+1
 516            zmdt(i,k,j) = stnd(1,kflip)/cpair
 517            zmdq(i,k,j) = qtnd(1,kflip)
 518         end do
 519
 520         !Top and bottom pressure of convection
 521         pconvt(i,j) = p8w(i,1,j)
 522         pconvb(i,j) = p8w(i,1,j)
 523         do n = 1,lengath(lchnk)
 524            if (maxg(n,lchnk).gt.jt(n,lchnk)) then
 525               pconvt(i,j) = pmid8(ideep(n,lchnk),jt(n,lchnk))  ! gathered array (or jctop ungathered)
 526               pconvb(i,j) = pmid8(ideep(n,lchnk),maxg(n,lchnk))! gathered array
 527            endif
 528         end do
 529         cutop(i,j) = jctop(1)
 530         cubot(i,j) = jcbot(1)
 531
 532         !Add tendency from this process to tendencies arrays. Also,
 533         !increment the local state arrays to reflect additional tendency
 534         !from zm_convr, i.e. the physics update call in CAM. Note that
 535         !we are not readjusting the pressure levels to hydrostatic
 536         !balance for the new virtual temperature like is done in CAM. We
 537         !will let WRF take care of such things later for the total
 538         !tendency.
 539         do k=kts,kte
 540            kflip = kte-k+1
 541
 542            !Convert temperature to potential temperature and
 543            !specific humidity to water vapor mixing ratio
 544            rthcuten(i,k,j) = zmdt(i,k,j)/pi_phy(i,k,j)
 545            rqvcuten(i,k,j) = zmdq(i,k,j)*(1._r8 + qv(i,k,j))**2
 546
 547            t8(1,kflip)    = t8(1,kflip) + zmdt(i,k,j)*cudts   !PMA
 548            qh8(1,kflip,1) = qh8(1,kflip,1) + zmdq(i,k,j)*cudts   !PMA
 549         end do
 550
 551         ! Determine the phase of the precipitation produced and add latent heat of fusion
 552         ! Evaporate some of the precip directly into the environment (Sundqvist)
 553         ! Allow this to use the updated state1 (t8 and qh8 in WRF) and the fresh ptend_loc type
 554         ! heating and specific humidity tendencies produced
 555         qtnd = 0._r8 !re-initialize tendencies (i.e. physics_ptend_init(ptend_loc))
 556         stnd = 0._r8
 557         call zm_conv_evap(ncol, lchnk, &
 558              t8, pmid8, pdel8, qh8(:,:,1), &
 559              stnd, tend_s_snwprd, tend_s_snwevmlt, qtnd, &
 560              rprd, cld8, cudts, &                            !PMA: This subroutine uses delt not 2xdelt, so pass in cudts
 561              prec, snow, ntprprd, ntsnprd , flxprec, flxsnow)
 562         evapcdp(:ncol,:pver) = qtnd(:ncol,:pver) !Added by Balwinder.Singh@pnnl.gov for assisting wetscavenging in WRF_CHEM package
 563
 564         ! Parse output variables from zm_conv_evap
 565         do k=kts,kte
 566            kflip = kte-k+1
 567
 568            evaptzm(i,k,j)  = stnd(1,kflip)/cpair
 569            fzsntzm(i,k,j)  = tend_s_snwprd(1,kflip)/cpair
 570            evsntzm(i,k,j)  = tend_s_snwevmlt(1,kflip)/cpair
 571            evapqzm(i,k,j)  = qtnd(1,kflip)
 572            zmflxprc(i,k,j) = flxprec(1,kflip)
 573            zmflxsnw(i,k,j) = flxsnow(1,kflip)
 574            zmntprpd(i,k,j) = ntprprd(1,kflip)
 575            zmntsnpd(i,k,j) = ntsnprd(1,kflip)
 576            zmeiheat(i,k,j) = stnd(1,kflip) !Do we really need this and evaptzm?
 577            !BSINGH - Following quantities are to be stored at interfaces
 578            !cmfmc(i,k,j)    = mcon(1,kflip) !Set to deep value here, shallow added in UW scheme
 579            !cmfmcdzm(i,k,j) = mcon(1,kflip)
 580            preccdzm(i,j)   = prec(1)       !Rain rate from just deep
 581            precz(i,j)      = prec(1)       !Rain rate for total convection (just deep right now)
 582            pratec(i,j)     = prec(1)*1e3   !Rain rate used in WRF for accumulation (mm/s)
 583            raincv(i,j)     = pratec(i,j)*dt !Rain amount for dynamic time step returned back to WRF   !wig: fixed wrong timestep usage 3-Jun-2016
 584         end do
 585
 586         !BSINGH - storing quantities at interfaces
 587         do k = kts,kte+1
 588            kflip             = kte - k + 2
 589            cmfmc(i,k,j)    = mcon(1,kflip) !Set to deep value here, shallow added in UW scheme
 590            cmfmcdzm(i,k,j) = mcon(1,kflip)
 591         end do
 592         !Add tendency from zm_conv_evap to tendencies arrays. Also,
 593         !increment the local state arrays to reflect additional tendency
 594         !Note that we are not readjusting the pressure levels to hydrostatic
 595         !balance for the new virtual temperature like is done in CAM. We
 596         !will let WRF take care of such things later for the total
 597         !tendency.
 598         do k=kts,kte
 599            kflip = kte-k+1
 600
 601            !Convert temperature to potential temperature and
 602            !specific humidity to water vapor mixing ratio
 603            rthcuten(i,k,j) = rthcuten(i,k,j) + &
 604                              evaptzm(i,k,j)/pi_phy(i,k,j)
 605            rqvcuten(i,k,j) = rqvcuten(i,k,j) + &
 606                              evapqzm(i,k,j)*(1. + qv(i,k,j))**2
 607
 608            t8(1,kflip)    = t8(1,kflip) + evaptzm(i,k,j)*cudts    !PMA
 609            qh8(1,kflip,1) = qh8(1,kflip,1) + evapqzm(i,k,j)*cudts !PMA
 610         end do
 611
 612         ! Momentum transport
 613         stnd       = 0._r8     !Zero relevant tendencies in preparation
 614         wind_tends = 0._r8
 615         do k=kts,kte
 616            kflip = kte-k+1
 617            winds(1,k,1) = u_phy(i,kflip,j)
 618            winds(1,k,2) = v_phy(i,kflip,j)
 619         end do
 620         l_windt(1:2) = .true.
 621
 622         call momtran (lchnk, ncol, &
 623              l_windt, winds, 2, mu(:,:,lchnk), md(:,:,lchnk), &
 624              du(:,:,lchnk), eu(:,:,lchnk), ed(:,:,lchnk), dp(:,:,lchnk), dsubcld(:,lchnk), &
 625              jt(:,lchnk),maxg(:,lchnk), ideep(:,lchnk), 1, lengath(lchnk), &
 626              itimestep, wind_tends, pguall, pgdall, icwu, icwd, hcudts, stnd )  !PMA
 627
 628         !Add tendency from momtran to tendencies arrays. Also,
 629         !increment the local state arrays to reflect additional tendency
 630         !Note that we are not readjusting the pressure levels to hydrostatic
 631         !balance for the new virtual temperature like is done in CAM. We
 632         !will let WRF take care of such things later for the total
 633         !tendency.
 634         do k=kts,kte
 635            kflip = kte-k+1
 636
 637            !Convert temperature to potential temperature and
 638            !specific humidity to water vapor mixing ratio
 639            rucuten(i,k,j)  = wind_tends(1,kflip,1)
 640            rvcuten(i,k,j)  = wind_tends(1,kflip,2)
 641            rthcuten(i,k,j) = rthcuten(i,k,j) + &
 642                                 stnd(1,kflip)/cpair/pi_phy(i,k,j)
 643
 644            t8(1,kflip) = t8(1,kflip) + stnd(1,kflip)/cpair*cudts   !PMA
 645            !winds is not used again so do not bother updating them
 646         end do
 647
 648         !Parse output arrays for momtran
 649         do k=kts,kte
 650            kflip = kte-k+1
 651
 652            zmmtu(i,k,j) = wind_tends(1,kflip,1) !wind tendencies...
 653            zmmtv(i,k,j) = wind_tends(1,kflip,2)
 654
 655            zmupgu(i,k,j) = pguall(1,kflip,1)   !apparent force pres. grads...
 656            zmupgd(i,k,j) = pgdall(1,kflip,1)
 657            zmvpgu(i,k,j) = pguall(1,kflip,2)
 658            zmvpgd(i,k,j) = pgdall(1,kflip,2)
 659
 660            zmicuu(i,k,j) = icwu(1,kflip,1)     !in-cloud winds...
 661            zmicud(i,k,j) = icwd(1,kflip,1)
 662            zmicvu(i,k,j) = icwu(1,kflip,2)
 663            zmicvd(i,k,j) = icwd(1,kflip,2)
 664         end do
 665
 666         !Setup for convective transport of cloud water and ice
 667         !~We should set this up to use an integer pointer instead of
 668         ! hard-coded numbers for each species so that we can easily
 669         ! handle other species. Something to come back to later...
 670         l_qt(1)   = .false.     !do not mix water vapor
 671         l_qt(2:3) = .true.      !do mix cloud water and ice
 672         cloudtnd = 0._r8
 673         fracis(1,:,1:3) = 1._r8 !all cloud liquid & ice
 674 #if ( WRF_CHEM != 1 )
 675         !BSINGH:02/01/2013: For WRF Physics ONLY runs, the liq # and ice # should
 676         !also be transpoted. Please note that liq # and ice # are transported in the
 677         !zm_conv_tend_2 call for WRF_CHEM simulations (ONLY works for MAM aerosols)
 678         !currently.
 679         fracis(1,:,4:5) = 1._r8
 680 #endif
 681         ncnst = 3               !number of constituents in cloud array (including vapor)
 682         fake_dpdry = 0._r8      !delta pres. for dry atmos. (0 since assuming moist mix ratios)
 683         do k=kts,kte
 684            kflip = kte-k+1
 685
 686            cloud(1,kflip,1) = qh8(1,kflip,1)  !We can either use moist mix ratios, as is
 687            cloud(1,kflip,2) = ql8prm(1,kflip)    !done here, or else use dry mix ratios, send
 688            cloud(1,kflip,3) = qi8(1,kflip)    !in appropriate dpdry values, and return the
 689                                               !approp. response from cnst_get_type_byind
 690         end do
 691
 692         call convtran (lchnk, &
 693              l_qt, cloud, ncnst,  mu(:,:,lchnk), md(:,:,lchnk), &
 694              du(:,:,lchnk), eu(:,:,lchnk), ed(:,:,lchnk), dp(:,:,lchnk), dsubcld(:,lchnk), &
 695              jt(:,lchnk), maxg(:,lchnk), ideep(:,lchnk), 1, lengath(lchnk), &
 696              itimestep, fracis, cloudtnd, fake_dpdry)
 697
 698         !Parse output for convtran and increment tendencies
 699 !PMA>
 700         do k=kts,kte
 701            kflip = kte-k+1
 702            esl(1,kflip)     = polysvp(t8(1,kflip),0)
 703            qvs(1,kflip)     = epsqs*esl(1,kflip)/(pmid8(1,kflip)-(1._r8-epsqs)*esl(1,kflip))
 704             if( t8(1,kflip) > 268.15_r8 ) then
 705               dum1 = 0.0_r8
 706             elseif( t8(1,kflip) < 238.15_r8 ) then
 707               dum1 = 1.0_r8
 708             else
 709               dum1 = ( 268.15_r8 - t8(1,kflip) ) / 30._r8
 710             endif
 711            qcten_det(1,kflip) = dlf(1,kflip) * ( 1._r8 - dum1 )
 712            qiten_det(1,kflip) = dlf(1,kflip) * dum1
 713            qcnten_det(1,kflip) = 3._r8 * (dlf(1,kflip)    * ( 1._r8 - dum1 ) ) / (4._r8*3.14159_r8*(8.e-6_r8**3)*997._r8)
 714            qinten_det(1,kflip) = 3._r8 * (dlf(1,kflip)    *  dum1 ) / (4._r8*3.14159_r8*(25.e-6_r8**3)*500._r8)
 715            qsten_det(1,kflip)  =  dlf(1,kflip) * dum1 * latice                    ! liquid to ice heating
 716         end do
 717         do k=kts,kte
 718            kflip = kte-k+1
 719
 720            !For assisting wetscavenging in WRF_CHEM package -Balwinder.Singh@pnnl.gov
 721            evapcdp3d(i,k,j) = evapcdp(1,kflip)
 722            icwmrdp3d(i,k,j) = ql8(1,kflip)
 723            rprddp3d(i,k,j)  = rprd(1,kflip)
 724
 725            zmdice(i,k,j) = cloudtnd(1,kflip,3)+qiten_det(1,kflip)
 726            zmdliq(i,k,j) = cloudtnd(1,kflip,2)+qcten_det(1,kflip)
 727            rthcuten(i,k,j) = rthcuten(i,k,j) + &
 728                             qsten_det(1,kflip)/cpair/pi_phy(i,k,j)
 729
 730            !Convert cloud tendencies from wet to dry mix ratios
 731            if( present(rqccuten) ) then
 732               rqccuten(i,k,j) = (cloudtnd(1,kflip,2)+qcten_det(1,kflip))*(1. + qv(i,k,j))
 733            end if
 734            if( present(rqicuten) ) then
 735               rqicuten(i,k,j) = (cloudtnd(1,kflip,3)+qiten_det(1,kflip))*(1. + qv(i,k,j))
 736            end if
 737            if( present(rqcncuten) ) then
 738               rqcncuten(i,k,j) = qcnten_det(1,kflip)*(1. + qv(i,k,j)) !BSINGH - Added the denominator following qiten_det
 739            end if
 740            if( present(rqincuten) ) then
 741               rqincuten(i,k,j) = qinten_det(1,kflip)*(1. + qv(i,k,j)) !BSINGH - Added the denominator following qiten_det
 742            end if
 743            !Variables required by zm_conv_tend_2 call
 744            dp3d(i,k,j) = dp(1,kflip,lchnk)
 745            du3d(i,k,j) = du(1,kflip,lchnk)
 746            ed3d(i,k,j) = ed(1,kflip,lchnk)
 747            eu3d(i,k,j) = eu(1,kflip,lchnk)
 748            md3d(i,k,j) = md(1,kflip,lchnk)
 749            mu3d(i,k,j) = mu(1,kflip,lchnk)
 750
 751            dsubcld2d(i,j) = dsubcld(1,lchnk)
 752            ideep2d(i,j)   = ideep(1,lchnk)
 753            jt2d(i,j)      = jt(1,lchnk)
 754            maxg2d(i,j)    = maxg(1,lchnk)
 755            lengath2d(i,j) = lengath(lchnk)
 756         end do
 757
 758      end do !i-loop
 759   end do    !j-loop
 760 END SUBROUTINE camzm_driver
 761
 762
 763 !------------------------------------------------------------------------
 764 SUBROUTINE zm_conv_init(DT, DX, rucuten, rvcuten, rthcuten, rqvcuten,   &
 765                      rqccuten, rqicuten, rqcncuten, rqincuten,          &
 766                      p_qc, p_qi, p_qnc, p_qni, param_first_scalar,      &
 767                      restart,                                           &
 768                      ids, ide, jds, jde, kds, kde,                      &
 769                      ims, ime, jms, jme, kms, kme,                      &
 770                      its, ite, jts, jte, kts, kte                    )
 771 ! Initialization routine for Zhang-McFarlane cumulus parameterization
 772 ! from CAM. The routine with this name in CAM is revamped here to give
 773 ! the needed functionality within WRF.
 774 !
 775 ! Author: William.Gustafson@pnl.gov, Nov. 2009
 776 !------------------------------------------------------------------------
 777   USE physconst, only: epsilo, latvap, latice, rh2o, cpair, tmelt
 778   USE module_cu_camzm, only: zm_convi, zmconv_readnl
 779   USE constituents, only: cnst_get_ind
 780
 781   LOGICAL , INTENT(IN)           ::   restart
 782   INTEGER , INTENT(IN)           ::   ids, ide, jds, jde, kds, kde, &
 783                                       ims, ime, jms, jme, kms, kme, &
 784                                       its, ite, jts, jte, kts, kte, &
 785                                       p_qc, p_qi, p_qnc, p_qni,     &
 786                                       param_first_scalar
 787
 788    REAL ,    INTENT(IN)        :: DT, DX
 789
 790
 791   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: &
 792                                                            rucuten, &
 793                                                            rvcuten, &
 794                                                           rthcuten, &
 795                                                           rqvcuten, &
 796                                                           rqccuten, &
 797                                                           rqicuten, &
 798                                                           rqcncuten,&
 799                                                           rqincuten
 800
 801   integer :: i, itf, j, jtf, k, ktf
 802   integer :: limcnv
 803
 804   jtf = min(jte,jde-1)
 805   ktf = min(kte,kde-1)
 806   itf = min(ite,ide-1)
 807
 808   call cnst_get_ind('CLDLIQ', ixcldliq)
 809   call cnst_get_ind('CLDICE', ixcldice)
 810   call cnst_get_ind('NUMLIQ', ixnumliq)
 811   call cnst_get_ind('NUMICE', ixnumice)
 812
 813 !
 814 ! Initialize module_cam_support variables...
 815 ! This could be moved to a master "cam-init" subroutine once multiple CAM
 816 ! parameterizations are implemented in WRF. For now, it doesn't hurt to
 817 ! have these possibly initialized more than once since they are
 818 ! essentially constant.
 819 !
 820   pver  = ktf-kts+1
 821   pverp = pver+1
 822 !~Need code here to set limcnv to max convective level of 40 mb... To
 823 ! properly set an average value for the whole domain would involve doing
 824 ! a collective operation across the whole domain to get pressure averages
 825 ! by level, which is difficult within the WRF framework. So, we will just
 826 ! use the model top for now.
 827 !
 828 ! Technically, limcnv should be calculated for each grid point at each
 829 ! time because the levels in WRF do not have a constant pressure, as
 830 ! opposed to CAM. But, that would involve making this variable local
 831 ! instead of at the module level as is currently done in CAM. For now,
 832 ! we will not worry about this because WRF rarely has domains higher than
 833 ! 40 mb. There is also little variability between grid points near the
 834 ! model top due to the underlying topography so a constant value would
 835 ! be OK across the domain. Also, remember that CAM levels are reversed in
 836 ! the vertical from WRF. So, setting limcnv to 1 means the top of the
 837 ! domain. Note that because of a bug in the parcel_dilute routine, limcnv
 838 ! must be >=2 instead of 1. Otherwise, an array out-of-bounds occurs.
 839   limcnv = 2
 840 !
 841 ! Get the ZM namelist variables (hard-wired for now)...
 842 !
 843   call zmconv_readnl("hard-wired")
 844 !
 845 !~need to determine if convection should happen in PBL and set
 846 ! no_deep_pbl_in accordingly
 847   call zm_convi(DT, DX, limcnv, NO_DEEP_PBL_IN=.true.)
 848
 849 !
 850 ! Set initial values for arrays dependent on restart conditions
 851 !
 852   if(.not.restart)then
 853      do j=jts,jtf
 854         do k=kts,ktf
 855            do i=its,itf
 856               rucuten(i,k,j)  = 0.
 857               rvcuten(i,k,j)  = 0.
 858               rthcuten(i,k,j) = 0.
 859               rqvcuten(i,k,j) = 0.
 860               if( p_qc  > param_first_scalar ) rqccuten(i,k,j)  = 0.
 861               if( p_qi  > param_first_scalar ) rqicuten(i,k,j)  = 0.
 862               if( p_qnc > param_first_scalar ) rqcncuten(i,k,j) = 0.
 863               if( p_qni > param_first_scalar ) rqincuten(i,k,j) = 0.
 864            enddo
 865         enddo
 866      enddo
 867   end if
 868
 869 END SUBROUTINE zm_conv_init
 870
 871
 872 !------------------------------------------------------------------------
 873 SUBROUTINE get_tpert(bl_pbl_physics, sf_sfclay_physics, dx, &
 874      mavail, kpbl, pblh, dzlowest, &
 875      psfc, qvlowest, t_phylowest, thlowest, tsk, tke_pbl, ust,   &
 876      u_phylowest, v_phylowest, hfx, qfx, tpert_camuwpbl, kte, tpert)
 877 ! Calculates the temperature perturbation used to trigger convection.
 878 ! This should only be used if tpert is not already provided by CAM's PBL
 879 ! scheme. Right now, this only works with the Mellor-Yamada boundary
 880 ! layer scheme in combination with the MYJ surface scheme. This is due to
 881 ! the need of TKE for the vertical velocity perturbation. This scheme has
 882 ! not been generalized to handle all schemes that produce TKE at this
 883 ! time, and the non-TKE schemes, e.g. YSU, could probably have an
 884 ! appropriate tpert deduced but we have not taken the time yet to do it.
 885 !
 886 ! Author: William.Gustafson@pnl.gov, Nov. 2009
 887 ! Last updated: William.Gustafson@pnl.gov, Nov. 2010
 888 !------------------------------------------------------------------------
 889   USE module_model_constants, only: cp, ep_1, ep_2, g, r_d, rcp, &
 890                                     svp1, svp2, svp3, svpt0, xlv
 891   USE module_state_description, ONLY : SFCLAYSCHEME              &
 892                                       ,MYJSFCSCHEME              &
 893                                       ,GFSSFCSCHEME              &
 894                                       ,SLABSCHEME                &
 895                                       ,LSMSCHEME                 &
 896                                       ,RUCLSMSCHEME              &
 897                                       ,MYJPBLSCHEME              &
 898                                       ,CAMUWPBLSCHEME
 899 !
 900 ! Subroutine arguments...
 901 !
 902   real, dimension(:), intent(in) :: tke_pbl
 903   real, intent(in) :: dx, dzlowest, hfx, mavail, pblh, psfc, qvlowest, &
 904        tpert_camuwpbl, tsk, t_phylowest, thlowest, ust, u_phylowest, &
 905        v_phylowest
 906   integer, intent(in) :: bl_pbl_physics, kpbl, kte, sf_sfclay_physics
 907   real(r8),intent(inout) :: tpert
 908 !
 909 ! Local vars...
 910 !
 911   real, parameter :: fak      = 8.5   !Constant in surface temperature excess
 912   real, parameter :: tfac     = 1.    !Ratio of 'tpert' to (w't')/wpert
 913   real, parameter :: wfac     = 1.    !Ratio of 'wpert' to sqrt(tke)
 914   real, parameter :: wpertmin = 1.e-6 !Min PBL eddy vertical velocity perturbation
 915   real :: ebrk                        !In CAM, net mean TKE of CL including
 916                                       !entrainment effect (m2/s2). In WRF,
 917                                       !average TKE within the PBL
 918   real :: br2, dthvdz, e1, flux, govrth, psfccmb, qfx, qsfc, rhox, thgb, &
 919        thv, tskv, tvcon, vconv, vsgd, wpert, wspd, za
 920   integer :: k
 921   character(len=250) :: msg
 922   logical :: UnstableOrNeutral
 923 !
 924 ! We can get the perturbation values directly from CAMUWPBLSCHEME if
 925 ! available. Otherwise, we have to calculate them.
 926 !
 927   if( bl_pbl_physics==CAMUWPBLSCHEME ) then
 928      tpert = tpert_camuwpbl
 929 !
 930 !...non-CAMUWPBL. Need MYJ SFC & PBL for now until other schemes
 931 !   get coded to give perturbations too.
 932 ! First, we need to determine if the conditions are stable or unstable.
 933 ! We will do this by mimicing the bulk Richardson calculation from
 934 ! module_sf_sfclay.F because the MYJ scheme does not provide this info.
 935 ! The logic borrowed from the CuP shallow cumulus scheme. Non-MYJ sfc
 936 ! scheme code is commented out for now since we also require MYJ PBL
 937 ! scheme.
 938 !
 939   elseif( bl_pbl_physics==MYJPBLSCHEME ) then
 940
 941      UnstableOrNeutral = .false.
 942      sfclay_case: SELECT CASE (sf_sfclay_physics)
 943      CASE (MYJSFCSCHEME)
 944         ! The MYJ sfc scheme does not already provide a stability criteria.
 945         ! So, we will mimic the bulk Richardson calculation from
 946         ! module_sf_sfclay.F.
 947
 948         if( pblh <= 0. ) call wrf_error_fatal( &
 949              "CAMZMSCHEME needs a PBL height from a PBL scheme.")
 950
 951         za     = 0.5*dzlowest
 952
 953         e1     = svp1*exp(svp2*(tsk-svpt0)/(tsk-svp3))
 954         psfccmb=psfc/1000.  !converts from Pa to cmb
 955         qsfc   = ep_2*e1/(psfccmb-e1)
 956         thgb   = tsk*(100./psfccmb)**rcp
 957         tskv   = thgb*(1.+ep_1*qsfc*mavail)
 958         tvcon  = 1.+ep_1*qvlowest
 959         thv    = thlowest*tvcon
 960         dthvdz = (thv-tskv)
 961
 962         govrth = g/thlowest
 963
 964         rhox   = psfc/(r_d*t_phylowest*tvcon)
 965         flux   = max(hfx/rhox/cp + ep_1*tskv*qfx/rhox,0.)
 966         vconv  = (g/tsk*pblh*flux)**.33
 967         vsgd   = 0.32 * (max(dx/5000.-1.,0.))**.33
 968         wspd   = sqrt(u_phylowest*u_phylowest+v_phylowest*v_phylowest)
 969         wspd   = sqrt(wspd*wspd+vconv*vconv+vsgd*vsgd)
 970         wspd   = max(wspd,0.1)
 971
 972         !And finally, the bulk Richardson number...
 973         br2   = govrth*za*dthvdz/(wspd*wspd)
 974
 975         if( br2 <= 0. ) UnstableOrNeutral = .true.
 976
 977      CASE DEFAULT
 978         call wrf_error_fatal("CAMZMSCHEME requires MYJSFCSCHEME or else CAMUWPBLSCHEME.")
 979
 980      END SELECT sfclay_case
 981 !
 982 ! The perturbation temperature for unstable conditions...
 983 ! The calculation follows the one in caleddy inside eddy_diff.F90 from
 984 ! CAM.
 985 !
 986      !Check that we are using the MJY BL scheme since we are hard-wired to
 987      !use TKE and u* from this scheme. At some point this dependency should
 988      !be broken and a way needs to be found for other schemes to provide
 989      !reasonable tpert values too.
 990      if( bl_pbl_physics /= MYJPBLSCHEME ) &
 991           call wrf_error_fatal("CAMZMSCHEME requires MYJPBLSCHEME or CAMUWPBLSCHEME")
 992
 993
 994      !Recalculate rhox in case a non-MYJ scheme is used to get
 995      !stability and rhox is not calculated above. Right now, this is
 996      !technically reduncant, but cheap.
 997      tvcon = 1.+ep_1*qvlowest
 998      rhox  = psfc/(r_d*t_phylowest*tvcon)
 999
1000      if( UnstableOrNeutral ) then
1001         !1st, get avg TKE w/n the PBL as a proxy for ebrk variable in CAM
1002         ebrk = 0.
1003         do k=1,min(kpbl,kte+1) !BSINGH - Changed KTE to KTE+1 as tke_pbl is @ interfaces
1004            ebrk = ebrk + tke_pbl(k)
1005         end do
1006         ebrk = ebrk/real(kpbl)
1007
1008         wpert = max( wfac*sqrt(ebrk), wpertmin )
1009         tpert = max( abs(hfx/rhox/cp)*tfac/wpert, 0. )
1010 !
1011 ! Or, the perturbation temperature for stable conditions...
1012 !
1013      else !Stable conditions
1014         tpert = max( hfx/rhox/cp*fak/ust, 0. )
1015      end if
1016
1017   else
1018      call wrf_error_fatal("CAMZMSCHEME requires MYJPBLSCHEME or CAMUWPBLSCHEME")
1019
1020   end if !PBL choice
1021
1022 END SUBROUTINE get_tpert
1023
1024
1025 !=========================================================================================
1026
1027 #if ( WRF_CHEM == 1 )
1028 subroutine zm_conv_tend_2(itimestep, dt, p8w, fracis3d, dp3d, du3d, ed3d, eu3d, md3d, mu3d, &
1029      dsubcld2d, ideep2d, jt2d, maxg2d, lengath2d, moist, scalar, chem,                      &
1030      ids,ide, jds,jde, kds,kde, &
1031      ims,ime, jms,jme, kms,kme, &
1032      its,ite, jts,jte, kts,kte)
1033
1034   use module_state_description,  only: num_moist, num_scalar, num_chem, param_first_scalar, &
1035        P_QV, P_QC, P_QI, P_QNC, P_QNI, CBMZ_CAM_MAM3_NOAQ, CBMZ_CAM_MAM3_AQ,                &
1036        CBMZ_CAM_MAM7_NOAQ,CBMZ_CAM_MAM7_AQ
1037   use module_data_cam_mam_asect, only: lptr_chem_to_q, factconv_chem_to_q
1038   use module_mp_cammgmp_driver,  only: physics_update, physics_ptend_init
1039
1040   implicit none
1041
1042   ! Arguments
1043
1044   !intent-ins
1045   integer, intent(in) :: itimestep
1046   integer, intent(in) :: ids,ide, jds,jde, kds,kde
1047   integer, intent(in) :: ims,ime, jms,jme, kms,kme
1048   integer, intent(in) :: its,ite, jts,jte, kts,kte
1049
1050   integer, dimension(ims:ime,jms:jme), intent(in) :: ideep2d, jt2d, maxg2d, lengath2d
1051
1052   real, intent(in) :: dt
1053
1054   real, dimension(ims:ime,jms:jme), intent(in) :: dsubcld2d
1055
1056   real, dimension(ims:ime,kms:kme,jms:jme), intent(in) :: p8w, dp3d, du3d, ed3d, eu3d, md3d
1057   real, dimension(ims:ime,kms:kme,jms:jme), intent(in) :: mu3d
1058
1059   real, dimension(ims:ime,kms:kme,jms:jme,pcnst), intent(in) :: fracis3d
1060
1061   !intent-inouts
1062   real, dimension(ims:ime,kms:kme,jms:jme,num_moist),  intent(inout) :: moist
1063   real, dimension(ims:ime,kms:kme,jms:jme,num_scalar), intent(inout) :: scalar
1064   real, dimension(ims:ime,kms:kme,jms:jme,num_chem),   intent(inout) :: chem
1065
1066
1067   ! Local variables
1068   character(len=1000) :: msg
1069   character*24 :: ptend_name            !ptend%name in CAM5 - used in parameterization
1070   logical      :: ptend_ls              !ptend%ls   in CAM5 - used for calling physics_update subroutine
1071   logical      :: ptend_lq(pcnst)       !ptend%lq   in CAM5
1072
1073   integer :: iw, jw, kw, kflip, l, l2
1074   integer :: i, lchnk, istat
1075   integer :: nstep
1076
1077   real(r8) :: delta_p, multFrc, dtime
1078   real(r8) :: dpdry(pcols,kte)
1079   real(r8) :: state_pdel(pcols,kte)                ! Pressure difference (Pa)
1080   real(r8) :: state_pdeldry(pcols,kte)             ! dry pressure difference (1/Pa)
1081   real(r8) :: state_q(pcols,kte,pcnst)
1082   real(r8) :: state_s(pcols,kte)
1083   real(r8) :: ptend_s(pcols,kte)                   !Dummy arguments for physics_update call
1084   real(r8) :: ptend_q(pcols,kte,pcnst)
1085
1086   real(r8) :: dp(pcols, kte, (ime-ims+1)*(jme-jms+1)) !layer pres. thickness between interfaces (mb)
1087   real(r8) :: du(pcols, kte, (ime-ims+1)*(jme-jms+1))
1088   real(r8) :: ed(pcols, kte, (ime-ims+1)*(jme-jms+1))
1089   real(r8) :: eu(pcols, kte, (ime-ims+1)*(jme-jms+1))
1090   real(r8) :: md(pcols, kte, (ime-ims+1)*(jme-jms+1))
1091   real(r8) :: mu(pcols, kte, (ime-ims+1)*(jme-jms+1))
1092
1093   real(r8) :: dsubcld(pcols, (ime-ims+1)*(jme-jms+1)) ! layer pres. thickness between LCL and maxi (mb)
1094
1095   integer :: ideep(pcols, (ime-ims+1)*(jme-jms+1)) ! holds position of gathered points
1096   integer :: jt(pcols, (ime-ims+1)*(jme-jms+1))    ! top-level index of deep cumulus convection
1097   integer :: maxg(pcols, (ime-ims+1)*(jme-jms+1))  ! gathered values of maxi
1098
1099   integer :: lengath((ime-ims+1)*(jme-jms+1))    ! number of gathered points
1100
1101
1102   ! physics buffer fields
1103   integer itim, ifld
1104   real(r8), dimension(pcols,kte,pcnst) :: fracis  ! fraction of transported species that are insoluble
1105
1106   !Time step is stored in the (r8) format in dtime
1107   dtime = dt
1108
1109   ! Map variables for CAM subrouine call
1110   ! Loop over the points in the tile and treat them each as a CAM chunk.
1111   !
1112
1113   !BSINGH:02/01/2013: Sanity check for Non-MAM simulations
1114   !This subroutine will not be called for chem_opt=0, so we dont need to worry about chem_opt = 0
1115   if(.NOT.cam_mam_aerosols) then
1116      write(msg,*)'CAMZM DRIVER - zm_conv_tend_2 is valid for only MAM aerosols ', &
1117           '(chem_opts:',CBMZ_CAM_MAM3_NOAQ,CBMZ_CAM_MAM3_AQ,CBMZ_CAM_MAM7_NOAQ,CBMZ_CAM_MAM7_AQ ,')'
1118      call wrf_error_fatal( msg )
1119   endif
1120
1121   do jw = jts,jte
1122      do iw = its,ite
1123         lchnk = (jw-jts)*(ite-its+1) + (iw-its+1) !1-D index location from the 2-D tile
1124
1125         !Flip variables on the layer middles
1126         do kw = kts,kte
1127            kflip = kte-kw+1
1128
1129            !Following three formulas are obtained from ported CAM's ZM cumulus scheme
1130            !Values of 0 cause a crash in entropy
1131            multFrc              = 1._r8/(1._r8 + moist(iw,kw,jw,P_QV))
1132            state_q(1,kflip,1)   = moist(iw,kw,jw,P_QV)*multFrc      !Specific humidity
1133            state_q(1,kflip,ixcldliq) = moist(iw,kw,jw,P_QC)*multFrc !Convert to moist mix ratio-cloud liquid
1134            state_q(1,kflip,ixcldice) = moist(iw,kw,jw,P_QI)*multFrc !cloud ice
1135            state_q(1,kflip,ixnumliq) = scalar(iw,kw,jw,P_QNC)*multFrc
1136            state_q(1,kflip,ixnumice) = scalar(iw,kw,jw,P_QNI)*multFrc
1137            do l = 1, 5
1138               fracis(1,kflip,l)      = fracis3d(iw,kw,jw,l)
1139            enddo
1140            !populate state_q and qqcw arrays
1141            !Following Do-Loop is obtained from chem/module_cam_mam_aerchem_driver.F
1142            do l = param_first_scalar, num_chem
1143               l2 = lptr_chem_to_q(l)
1144               if ((l2 >= 1) .and. (l2 <= pcnst)) then
1145                  fracis(1,kflip,l2)  = fracis3d(iw,kw,jw,l2)
1146                  state_q(1,kflip,l2) = chem(iw,kw,jw,l)*factconv_chem_to_q(l)
1147               end if
1148            end do ! l
1149
1150            delta_p                 = p8w(iw,kw,jw) - p8w(iw,kw+1,jw) !Change in pressure (Pa)
1151            state_pdel(  1,kflip)   = delta_p
1152            state_pdeldry(1,kflip)  = state_pdel(1,kflip)*(1._r8-state_q(1,kflip,1))
1153            !Variables required by zm_conv_tend_2 call
1154            dp(1,kflip,lchnk) = dp3d(iw,kw,jw)
1155            du(1,kflip,lchnk) = du3d(iw,kw,jw)
1156            ed(1,kflip,lchnk) = ed3d(iw,kw,jw)
1157            eu(1,kflip,lchnk) = eu3d(iw,kw,jw)
1158            md(1,kflip,lchnk) = md3d(iw,kw,jw)
1159            mu(1,kflip,lchnk) = mu3d(iw,kw,jw)
1160
1161            dsubcld(1,lchnk)  = dsubcld2d(iw,jw)
1162            ideep(1,lchnk)    = ideep2d(iw,jw)
1163            jt(1,lchnk)       = jt2d(iw,jw)
1164            maxg(1,lchnk)     = maxg2d(iw,jw)
1165            lengath(lchnk)    = lengath2d(iw,jw)
1166         enddo
1167
1168         !
1169         ! Initialize
1170         !
1171         call physics_ptend_init(ptend_name,ptend_q,ptend_s,ptend_lq,ptend_ls,pcnst)  !! Initialize local ptend type
1172
1173         !
1174         ! Transport all constituents except cloud water and ice
1175         !
1176
1177         nstep = itimestep
1178
1179         !
1180         !     Convective transport of all trace species except cloud liquid
1181         !     and cloud ice done here because we need to do the scavenging first
1182         !     to determine the interstitial fraction.
1183         !
1184         ptend_name  = 'convtran2'
1185         ptend_lq(:) = .true.
1186         ptend_lq(1)        = .false.
1187         ptend_lq(ixcldice) = .false.
1188         ptend_lq(ixcldliq) = .false.
1189
1190         ! initialize dpdry for call to convtran
1191         ! it is used for tracers of dry mixing ratio type
1192         dpdry = 0._r8
1193         do i = 1,lengath(lchnk)
1194            dpdry(i,:) = state_pdeldry(ideep(i,lchnk),:)/100._r8
1195         end do
1196
1197         call convtran (lchnk,                                        &
1198              ptend_lq,state_q, pcnst,  mu(:,:,lchnk), md(:,:,lchnk),   &
1199              du(:,:,lchnk), eu(:,:,lchnk), ed(:,:,lchnk), dp(:,:,lchnk), dsubcld(:,lchnk),  &
1200              jt(:,lchnk),maxg(:,lchnk),ideep(:,lchnk), 1, lengath(lchnk),  &
1201              nstep,   fracis,  ptend_q, dpdry)
1202
1203         !Update chem array and state constiuents
1204         !populate state_s, ptend_s, ptend_ls with dummy values (zeros) for physics update call
1205         state_s(:,:) = 0.0_r8
1206         ptend_s(:,:) = 0.0_r8
1207         ptend_ls     = .FALSE.
1208
1209         call physics_update(lchnk,dtime,state_q,ptend_q,state_s,ptend_s,ptend_name,ptend_lq,ptend_ls,pcnst)
1210
1211         !Post processing of the output to update WRF arrays
1212           do kw=kts,kte
1213
1214              kflip = kte-kw+1
1215
1216              !Following equation are derived following UWPBL and CAMZM schemes
1217              moist(iw,kw,jw,P_QV)   = state_q(1,kflip,1) / (1.0_r8 - state_q(1,kflip,1))
1218              multFrc                = 1._r8 + moist(iw,kw,jw,P_QV)
1219
1220              moist(iw,kw,jw,P_QC)   = state_q(1,kflip,2) * multFrc
1221              moist(iw,kw,jw,P_QI)   = state_q(1,kflip,3) * multFrc
1222
1223              scalar(iw,kw,jw,P_QNC) = state_q(1,kflip,4) * multFrc
1224              scalar(iw,kw,jw,P_QNI) = state_q(1,kflip,5) * multFrc
1225
1226              do l = param_first_scalar, num_chem
1227                 l2 = lptr_chem_to_q(l)
1228                 if ((l2 >= 1) .and. (l2 <= pcnst)) then
1229                    chem(iw,kw,jw,l) = state_q(1,kflip,l2)/factconv_chem_to_q(l)
1230                 end if
1231              end do ! l
1232           end do
1233
1234
1235
1236      enddo
1237   enddo
1238 end subroutine zm_conv_tend_2
1239 #endif
1240 END MODULE module_cu_camzm_driver