phys/module_mp_wdm6.F

   1 #if ( RWORDSIZE == 4 )
   2 #  define VREC vsrec
   3 #  define VSQRT vssqrt
   4 #else
   5 #  define VREC vrec
   6 #  define VSQRT vsqrt
   7 #endif
   8 !-------------------------------------------------------------------------------
   9    module module_mp_wdm6
  10 !-------------------------------------------------------------------------------
  11 !
  12    use module_mp_radar
  13    use module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG
  14 !
  15 !
  16 !-------------------------------------------------------------------------------
  17 !
  18 !  constant   :
  19 !
  20 !    dtcldcr       - maximum time step for minor loops
  21 !    n0r           - intercept parameter rain
  22 !    n0s           - intercept parameter snow
  23 !    n0g           - intercept parameter graupel
  24 !    n0smax        - maximum n0s (t=-90C unlimited)
  25 !    alpha         - 0.122 exponen factor for n0s
  26 !    avtr, bvtr    - a constant for terminal velocity of rain
  27 !    avts, bvts    - a constant for terminal velocity of snow
  28 !    avtg, bvtg    - a constant for terminal velocity of graupel
  29 !    lamdarmax     - limited maximum value for slope parameter of rain
  30 !    lamdarmin     - limited minimum value for slope parameter of rain
  31 !    lamdasmax     - limited maximum value for slope parameter of snow
  32 !    lamdagmax     - limited maximum value for slope parameter of graupel
  33 !    r0            - 8 microm  in contrast to 10 micro m
  34 !    peaut         - collection efficiency
  35 !    xncr          - maritime cloud in contrast to 3.e8 in tc80
  36 !    xmyu          - the dynamic viscosity kg m-1 s-1
  37 !    dicon         - constant for the cloud-ice diamter
  38 !    dimax         - limited maximum value for the cloud-ice diamter
  39 !    pfrz1, pfrz2  - constant in Biggs freezing
  40 !    qcrmin        - minimum values for qr and qs
  41 !    ncmin         - minimum value for nc
  42 !    nrmin         - minimum value for nr
  43 !    eacrc         - now/cloud-water collection efficiency
  44 !    qs0           - threshold amount for aggretion to occur
  45 !    satmax        - maximum saturation value for CCN activation
  46 !                    1.008 for maritime /1.0048 for conti
  47 !    actk          - parameter for the CCN activation
  48 !    actr          - radius of activated CCN drops
  49 !    ncrk1, ncrk2  - Long's collection kernel coefficient
  50 !    di100, di600, di2000  - parameter related with accretion and collection
  51 !                            of cloud drops
  52 !    di82          - dimater related with raindrops evaporation
  53 !    di15          - auto conversion takes place beyond this diameter
  54 !
  55 !-------------------------------------------------------------------------------
  56    real, parameter, private :: dtcldcr = 120.
  57    real, parameter, private :: n0r     = 8.e6
  58    real, parameter, private :: n0s     = 2.e6
  59 !   real, parameter, private :: n0g     = 4.e6
  60    real, parameter, private :: n0smax  =  1.e11
  61    real, parameter, private :: dens    = 100.0
  62    real, parameter, private :: alpha   = .12
  63    real, parameter, private :: avtr    = 841.9
  64    real, parameter, private :: bvtr    = 0.8
  65    real, parameter, private :: avts    = 11.72
  66    real, parameter, private :: bvts    = .41
  67 !   real, parameter, private :: avtg    = 330.
  68 !   real, parameter, private :: bvtg    = 0.8
  69    real, parameter, private :: lamdacmax = 5.e5
  70    real, parameter, private :: lamdacmin = 2.e4
  71    real, parameter, private :: lamdarmax = 5.e4
  72    real, parameter, private :: lamdarmin = 2.e3
  73    real, parameter, private :: lamdasmax = 1.e5
  74 !   real, parameter, private :: lamdagmax = 6.e4
  75    real, parameter, private :: r0      = .8e-5
  76    real, parameter, private :: peaut   = .55
  77    real, parameter, private :: xncr    = 3.e8
  78    real, parameter, private :: xncr0   = 5.e7
  79    real, parameter, private :: xncr1   = 5.e8
  80    real, parameter, private :: xmyu    = 1.718e-5
  81    real, parameter, private :: dicon   = 11.9
  82    real, parameter, private :: dimax   = 500.e-6
  83    real, parameter, private :: pfrz1   = 100.
  84    real, parameter, private :: pfrz2   = 0.66
  85    real, parameter, private :: qcrmin  = 1.e-9
  86    real, parameter, private :: ncmin   = 1.e1
  87    real, parameter, private :: nrmin   = 1.e-2
  88    real, parameter, private :: eacrc   = 1.0
  89    real, parameter, private :: qs0     = 6.e-4
  90    real, parameter, private :: satmax  = 1.0048
  91    real, parameter, private :: actk    = 0.6
  92    real, parameter, private :: actr    = 1.5
  93    real, parameter, private :: ncrk1   = 3.03e3
  94    real, parameter, private :: ncrk2   = 2.59e15
  95    real, parameter, private :: di100   = 1.e-4
  96    real, parameter, private :: di600   = 6.e-4
  97    real, parameter, private :: di2000  = 2000.e-6
  98    real, parameter, private :: di82    = 82.e-6
  99    real, parameter, private :: di15    = 15.e-6
 100    real, save ::                                                               &
 101              qc0,qc1,qck1,pidnc,bvtr1,bvtr2,bvtr3,bvtr4,bvtr5,                 &
 102              bvtr6,bvtr7, bvtr2o5,bvtr3o5,                                     &
 103              g1pbr,g2pbr,g3pbr,g4pbr,g5pbr,g6pbr,g7pbr,                        &
 104              g5pbro2,g7pbro2,pi,                                               &
 105              pvtr,pvtrn,eacrr,pacrr,pidn0r,pidnr,                              &
 106              precr1,precr2,xmmax,roqimax,bvts1,bvts2,                          &
 107              bvts3,bvts4,g1pbs,g3pbs,g4pbs,g5pbso2,                            &
 108              pvts,pacrs,precs1,precs2,pidn0s,xlv1,pacrc,                       &
 109              bvtg1,bvtg2,bvtg3,bvtg4,g1pbg,g3pbg,g4pbg,                        &
 110              g5pbgo2,pvtg,pacrg,precg1,precg2,pidn0g,                          &
 111              n0g,avtg,bvtg,deng,lamdagmax,                                     &
 112              rslopecmax,rslopec2max,rslopec3max,                               &
 113              rslopermax,rslopesmax,rslopegmax,                                 &
 114              rsloperbmax,rslopesbmax,rslopegbmax,                              &
 115              rsloper2max,rslopes2max,rslopeg2max,                              &
 116              rsloper3max,rslopes3max,rslopeg3max
 117 !
 118     contains
 119 !===================================================================
 120 !
 121     subroutine wdm6(th, q, qc, qr, qi, qs, qg,             &
 122                     nn, nc, nr,                            &
 123                     den, pii, p, delz,                     &
 124                     delt,g, cpd, cpv, ccn0, rd, rv, t0c,   &
 125                     ep1, ep2, qmin,                        &
 126                     xls, xlv0, xlf0, den0, denr,           &
 127                     cliq,cice,psat,                        &
 128                     xland,                                 &
 129                     rain, rainncv,                         &
 130                     snow, snowncv,                         &
 131                     sr,                                    &
 132                     refl_10cm, diagflag, do_radar_ref,     &
 133                     graupel, graupelncv,                   &
 134                     itimestep,                             &
 135                     has_reqc, has_reqi, has_reqs,          &  ! for radiation
 136                     re_cloud, re_ice,   re_snow,           &  ! for radiation
 137                     ids,ide, jds,jde, kds,kde,             &
 138                     ims,ime, jms,jme, kms,kme,             &
 139                     its,ite, jts,jte, kts,kte              &
 140                                                            )
 141 !-------------------------------------------------------------------
 142    implicit none
 143 !-------------------------------------------------------------------
 144    integer                                   , intent(in   ) :: itimestep
 145    integer                                   , intent(in   ) ::                &
 146                                                     ids,ide, jds,jde, kds,kde, &
 147                                                     ims,ime, jms,jme, kms,kme, &
 148                                                     its,ite, jts,jte, kts,kte
 149    real                                      , intent(in   ) :: delt, g, ccn0, &
 150                                                                 rd, rv, t0c,   &
 151                                                                 cpd, cpv,      &
 152                                                                 den0, qmin,    &
 153                                                                 ep1, ep2, xls, &
 154                                                                 xlv0, xlf0,    &
 155                                                                 cliq, cice,    &
 156                                                                 psat, denr
 157    real, dimension(ims:ime,jms:jme)          , intent(in   ) :: xland
 158    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(in   ) :: den
 159    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(in   ) :: pii
 160    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(in   ) :: p
 161    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(in   ) :: delz
 162    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(inout) :: th
 163    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(inout) :: q
 164    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(inout) :: qc
 165    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(inout) :: qi
 166    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(inout) :: qr
 167    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(inout) :: qs
 168    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(inout) :: qg
 169    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(inout) :: nn
 170    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(inout) :: nc
 171    real, dimension(ims:ime,kms:kme,jms:jme)  , intent(inout) :: nr
 172    real, dimension(ims:ime,jms:jme)          , intent(inout) :: rain
 173    real, dimension(ims:ime,jms:jme)          , intent(inout) :: rainncv
 174    real, dimension(ims:ime,jms:jme)          , intent(inout) :: sr
 175    real, dimension(ims:ime,jms:jme), optional, intent(inout) :: snow
 176    real, dimension(ims:ime,jms:jme), optional, intent(inout) :: snowncv
 177    real, dimension(ims:ime,jms:jme), optional, intent(inout) :: graupel
 178    real, dimension(ims:ime,jms:jme), optional, intent(inout) :: graupelncv
 179 !
 180 ! for radiation connecting
 181 !
 182    integer                                 , intent(in   ) :: has_reqc
 183    integer                                 , intent(in   ) :: has_reqi
 184    integer                                 , intent(in   ) :: has_reqs
 185    real, dimension(ims:ime,kms:kme,jms:jme), intent(inout) :: re_cloud
 186    real, dimension(ims:ime,kms:kme,jms:jme), intent(inout) :: re_ice
 187    real, dimension(ims:ime,kms:kme,jms:jme), intent(inout) :: re_snow
 188 !
 189 ! for reflectivity
 190 !
 191    real, dimension(ims:ime,kms:kme,jms:jme), intent(inout) :: refl_10cm
 192    logical                       , optional, intent(in   ) :: diagflag
 193    integer                       , optional, intent(in   ) :: do_radar_ref
 194 !
 195 ! local variables
 196 !
 197    real, dimension(its:ite,kts:kte)   :: t
 198    real, dimension(its:ite,kts:kte,2) :: qci
 199    real, dimension(its:ite,kts:kte,3) :: qrs, ncr
 200    integer                            :: i,j,k
 201 !
 202 ! for reflectivity
 203 !
 204    real, dimension(kts:kte) :: qv1d
 205    real, dimension(kts:kte) :: t1d
 206    real, dimension(kts:kte) :: p1d
 207    real, dimension(kts:kte) :: qr1d
 208    real, dimension(kts:kte) :: nr1d
 209    real, dimension(kts:kte) :: qs1d
 210    real, dimension(kts:kte) :: qg1d
 211    real, dimension(kts:kte) :: dBZ
 212 !
 213 ! to calculate effective radius for radiation
 214 !
 215    real, dimension(kts:kte) :: qc1d, nc1d, den1d
 216    real, dimension(kts:kte) :: qi1d
 217    real, dimension(kts:kte) :: re_qc, re_qi, re_qs
 218 !
 219    if (itimestep .eq. 1) then
 220      do j = jms,jme
 221        do k = kms,kme
 222          do i = ims,ime
 223            nn(i,k,j) = ccn0
 224          enddo
 225        enddo
 226      enddo
 227    endif
 228 !
 229    do j = jts,jte
 230      do k = kts,kte
 231        do i = its,ite
 232          t(i,k) = th(i,k,j)*pii(i,k,j)
 233          qci(i,k,1) = qc(i,k,j)
 234          qci(i,k,2) = qi(i,k,j)
 235          qrs(i,k,1) = qr(i,k,j)
 236          qrs(i,k,2) = qs(i,k,j)
 237          qrs(i,k,3) = qg(i,k,j)
 238          ncr(i,k,1) = nn(i,k,j)
 239          ncr(i,k,2) = nc(i,k,j)
 240          ncr(i,k,3) = nr(i,k,j)
 241        enddo
 242      enddo
 243      !
 244      !  Sending array starting locations of optional variables may cause
 245      !  troubles, so we explicitly change the call.
 246      !
 247      call wdm62D(t, q(ims,kms,j), qci, qrs, ncr               &
 248                 ,den(ims,kms,j)                               &
 249                 ,p(ims,kms,j), delz(ims,kms,j)                &
 250                 ,delt,g, cpd, cpv, ccn0, rd, rv, t0c          &
 251                 ,ep1, ep2, qmin                               &
 252                 ,xls, xlv0, xlf0, den0, denr                  &
 253                 ,cliq,cice,psat                               &
 254                 ,j                                            &
 255                 ,xland(ims,j)                                 &
 256                 ,rain(ims,j),rainncv(ims,j)                   &
 257                 ,sr(ims,j)                                    &
 258                 ,ids,ide, jds,jde, kds,kde                    &
 259                 ,ims,ime, jms,jme, kms,kme                    &
 260                 ,its,ite, jts,jte, kts,kte                    &
 261                 ,snow(ims,j),snowncv(ims,j)                   &
 262                 ,graupel(ims,j),graupelncv(ims,j)             &
 263                                                                )
 264 !
 265      do k = kts,kte
 266        do i = its,ite
 267          th(i,k,j) = t(i,k)/pii(i,k,j)
 268          qc(i,k,j) = qci(i,k,1)
 269          qi(i,k,j) = qci(i,k,2)
 270          qr(i,k,j) = qrs(i,k,1)
 271          qs(i,k,j) = qrs(i,k,2)
 272          qg(i,k,j) = qrs(i,k,3)
 273          nn(i,k,j) = ncr(i,k,1)
 274          nc(i,k,j) = ncr(i,k,2)
 275          nr(i,k,j) = ncr(i,k,3)
 276        enddo
 277      enddo
 278 !
 279      if ( present (diagflag) ) then
 280        if (diagflag .and. do_radar_ref == 1) then
 281          do i = its,ite
 282            do k = kts,kte
 283              t1d(k)  = th(i,k,j)*pii(i,k,j)
 284              p1d(k)  = p(i,k,j)
 285              qv1d(k) = q(i,k,j)
 286              qr1d(k) = qr(i,k,j)
 287              nr1d(k) = nr(i,k,j)
 288              qs1d(k) = qs(i,k,j)
 289              qg1d(k) = qg(i,k,j)
 290            enddo
 291            call refl10cm_wdm6 (qv1d, qr1d, nr1d, qs1d, qg1d,        &
 292                                t1d, p1d, dBZ, kts, kte, i, j)
 293            do k = kts, kte
 294              refl_10cm(i,k,j) = max(-35., dBZ(k))
 295            enddo
 296          enddo
 297        endif
 298      endif
 299 !
 300 ! calculate effective radius of cloud, ice, and snow
 301 !
 302      if (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) then
 303        do i = its,ite
 304          do k = kts,kte
 305            re_qc(k) = RE_QC_BG
 306            re_qi(k) = RE_QI_BG
 307            re_qs(k) = RE_QS_BG
 308 !
 309            t1d(k)  = th(i,k,j)*pii(i,k,j)
 310            den1d(k)= den(i,k,j)
 311            qc1d(k) = qc(i,k,j)
 312            qi1d(k) = qi(i,k,j)
 313            qs1d(k) = qs(i,k,j)
 314            nc1d(k) = nc(i,k,j)
 315          enddo
 316 !
 317          call effectRad_wdm6(t1d, qc1d, nc1d, qi1d, qs1d, den1d,   &
 318                              qmin, t0c, re_qc, re_qi, re_qs,       &
 319                              kts, kte, i, j)
 320 !
 321          do k = kts,kte
 322            re_cloud(i,k,j) = max(RE_QC_BG, min(re_qc(k),  50.e-6))
 323            re_ice(i,k,j)   = max(RE_QI_BG, min(re_qi(k), 125.e-6))
 324            re_snow(i,k,j)  = max(RE_QS_BG, min(re_qs(k), 999.e-6))
 325          enddo
 326        enddo
 327      endif
 328    enddo
 329 !
 330    end subroutine wdm6
 331 !
 332 !------------------------------------------------------------------------------
 333 !
 334    subroutine wdm62D(t, q, qci, qrs, ncr, den, p, delz            &
 335                    ,delt,g, cpd, cpv, ccn0, rd, rv, t0c           &
 336                    ,ep1, ep2, qmin                                &
 337                    ,xls, xlv0, xlf0, den0, denr                   &
 338                    ,cliq,cice,psat                                &
 339                    ,lat                                           &
 340                    ,slmsk                                         &
 341                    ,rain,rainncv                                  &
 342                    ,sr                                            &
 343                    ,ids,ide, jds,jde, kds,kde                     &
 344                    ,ims,ime, jms,jme, kms,kme                     &
 345                    ,its,ite, jts,jte, kts,kte                     &
 346                    ,snow,snowncv                                  &
 347                    ,graupel,graupelncv                            &
 348                                                                   )
 349 !------------------------------------------------------------------------------
 350     implicit none
 351 !-------------------------------------------------------------------------------
 352 !
 353 !  This code is a WRF double-moment 6-class GRAUPEL phase
 354 !  microphyiscs scheme (wdm6). The WDM microphysics scheme predicts
 355 !  number concentrations for warm rain species including clouds and
 356 !  rain. cloud condensation nuclei (ccn) is also predicted.
 357 !  The cold rain species including ice, snow, graupel follow the
 358 !  WRF single-moment 6-class microphysics (wsm6, Hong and Lim 2006)
 359 !  in which theoretical background for WSM ice phase microphysics is
 360 !  based on Hong et al. (2004). A new mixed-phase terminal velocity
 361 !  for precipitating ice is introduced in WSM6 (Dudhia et al. 2008).
 362 !  The WDM scheme is described in Lim and Hong (2010).
 363 !
 364 !  wdm6 cloud scheme
 365 !
 366 !  Coded by Kyo-Sun Lim and Song-You Hong (Yonsei Univ.) Fall 2008
 367 !
 368 !  Implemented by Kyo-Sun Lim and Jimy Dudhia (NCAR) Winter 2008
 369 !
 370 !  further modifications :
 371 !        semi-lagrangian sedimentation (JH,2010),hong, aug 2009
 372 !        ==> higher accuracy and efficient at lower resolutions
 373 !        reflectivity computation from greg thompson, lim, jun 2011
 374 !        ==> only dignostic, but with removal of too large drops
 375 !        add hail option from afwa, aug 2014
 376 !        ==> switch graupel or hail by changing no, den, fall vel.
 377 !        effetive radius of hydrometeors, bae from kiaps, jan 2015
 378 !        ==> consistency in solar insolation of rrtmg radiation
 379 !        bug fix in melting terms, bae from kiaps, nov 2015
 380 !        ==> density of air is divided, which has not been
 381 !        bug fix in sedimentation of rain, bae from kiaps, mar 2017
 382 !        change autoconversion rate to equation, bae from kiaps, oct 2017
 383 !        nccn > 100 cm-3, bae from kiaps, oct 2017
 384 !        bug fix hydrometeros update before condensation process of rain,
 385 !           bae from kiaps, oct2017
 386 !        bug fix in snow melting and autoconversion.
 387 !        ==> Revisions enchance melting and autoconverison.
 388 !           hong from noaa, july 2023
 389 !        remedy for tiny nc in the presence of qc (nc = xncr, 300/cm**3)
 390 !        ==> prevents the crash in computing autoconversion
 391 !           hong from noaa, july 2023
 392 !
 393 !  history log :
 394 !
 395 !  Reference) Lim and Hong (LH, 2010) Mon. Wea. Rev.
 396 !             Juang and Hong (JH, 2010) Mon. Wea. Rev.
 397 !             Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
 398 !             Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
 399 !             Cohard and Pinty (CP, 2000) Quart. J. Roy. Meteor. Soc.
 400 !             Khairoutdinov and Kogan (KK, 2000) Mon. Wea. Rev.
 401 !             Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
 402 !
 403 !             Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
 404 !             Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
 405 !             Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
 406 !
 407 !-------------------------------------------------------------------------------
 408 !
 409 !  input :
 410 !    deltim                       - timestep
 411 !    g, cpd, cpv, t0c, den0,      - constant
 412 !    rd, rv, ep1, ep2, qmin,
 413 !    xls, xlv0, xlf0, denr,
 414 !    cliq, cice, psat
 415 !    ids, ide, jds, jde, kds, kde - dimension
 416 !    ims, ime, jms, jme, kms, kme
 417 !    its, ite, jts, jte, kts, kte
 418 !    ncloud                       - number of hydrometeor
 419 !    p                            - pressure, Pa
 420 !    delz                         - depth of model layer, m
 421 !   * with chemistry
 422 !    slmsk                        - land/sea mask, 0: sea, 1: land
 423 !    so4                          - mixing ratio of sulfate
 424 !    ss                           - mixing ratio of sea salt
 425 !    oc                           - mixing ratio of organic carbon
 426 !
 427 !  inout :
 428 !    t1                           - temperautre
 429 !    q1                           - specific humidity
 430 !    q2                           - mixing ratio of cloud water, rain, ice,
 431 !                                   snow, and graupel
 432 !                                   qc, qr, qi, qs, qg
 433 !    n2                           - number concentration of ccn, cloud droplet,
 434 !                                   and raindrop
 435 !                                   nn, nc, nr
 436 !    rain, rainncv                - precipitation
 437 !    sr                           - ratio of snow to rain
 438 !
 439 !  all units are in m.k.s. and source/sink terms in kg kg-1 s-1.
 440 !-------------------------------------------------------------------------------
 441    integer                                , intent(in   ) ::                  &
 442                                                    ids,ide, jds,jde, kds,kde, &
 443                                                    ims,ime, jms,jme, kms,kme, &
 444                                                    its,ite, jts,jte, kts,kte
 445    integer                                , intent(in   ) :: lat
 446    real                                   , intent(in   ) :: delt
 447    real                                   , intent(in   ) :: g, cpd, cpv, t0c, &
 448                                                              den0, rd, rv,     &
 449                                                              ep1, ep2, qmin,   &
 450                                                              xls, xlv0, xlf0,  &
 451                                                              cliq, cice, psat, &
 452                                                              denr
 453    real                                   , intent(in   ) :: ccn0
 454    real, dimension(ims:ime)               , intent(in   ) :: slmsk
 455    real, dimension(ims:ime,kms:kme)       , intent(in   ) :: p
 456    real, dimension(ims:ime,kms:kme)       , intent(in   ) :: delz
 457    real, dimension(ims:ime,kms:kme)       , intent(in   ) :: den
 458    real, dimension(its:ite,kts:kte)       , intent(inout) :: t
 459    real, dimension(its:ite,kts:kte,2)     , intent(inout) :: qci
 460    real, dimension(its:ite,kts:kte,3)     , intent(inout) :: qrs
 461    real, dimension(its:ite,kts:kte,3)     , intent(inout) :: ncr
 462    real, dimension(ims:ime,kms:kme)       , intent(inout) :: q
 463    real, dimension(ims:ime)               , intent(inout) :: rain
 464    real, dimension(ims:ime)               , intent(inout) :: rainncv
 465    real, dimension(ims:ime)               , intent(inout) :: sr
 466    real, dimension(ims:ime), optional     , intent(inout) :: snow
 467    real, dimension(ims:ime), optional     , intent(inout) :: snowncv
 468    real, dimension(ims:ime), optional     , intent(inout) :: graupel
 469    real, dimension(ims:ime), optional     , intent(inout) :: graupelncv
 470 !
 471 !  local variables
 472 !
 473    real, dimension(its:ite,kts:kte)   :: dend
 474    real, dimension(its:ite,kts:kte)   :: qcr
 475    real, dimension(its:ite,kts:kte,3) :: rh
 476    real, dimension(its:ite,kts:kte,3) :: qs
 477    real, dimension(its:ite,kts:kte,3) :: rslope, rslope2, rslope3, rslopeb
 478    real, dimension(its:ite,kts:kte,3) :: falk, fall
 479    real, dimension(its:ite,kts:kte,3) :: work1
 480    real, dimension(its:ite,kts:kte,3) :: qrs_tmp
 481    real, dimension(its:ite,kts:kte,2) :: avedia
 482    real, dimension(its:ite,kts:kte)   :: rslopec, rslopec2,rslopec3
 483    real, dimension(its:ite,kts:kte)   :: workn, falln, falkn
 484    real, dimension(its:ite,kts:kte)   :: worka, workr
 485    real, dimension(its:ite,kts:kte)   :: den_tmp, delz_tmp, ncr_tmp
 486    real, dimension(its:ite,kts:kte)   :: lamdr_tmp
 487    real, dimension(its:ite,kts:kte)   :: lamdc_tmp
 488    real, dimension(its:ite,kts:kte)   :: falkc, work1c, work2c, fallc
 489    real, dimension(its:ite,kts:kte)   :: dqr,dnr
 490    real, dimension(its:ite,kts:kte)   :: pcact, prevp, psdep, pgdep, praut,    &
 491                                          psaut, pgaut, pracw, psacw, pgacw,    &
 492                                          pgacr, pgacs, psaci, pgmlt, praci,    &
 493                                          piacr, pracs, psacr, pgaci, pseml,    &
 494                                          pgeml, paacw, pigen, pidep, pcond,    &
 495                                          psmlt, psevp, pgevp
 496    real, dimension(its:ite,kts:kte)   :: nraut, nracw, ncevp, nccol, nrcol,    &
 497                                          nsacw, ngacw, niacr, nsacr, ngacr,    &
 498                                          naacw, nseml, ngeml, ncact
 499    real, dimension(its:ite,kts:kte)   :: xl, cpm, work2, denfac, xni, n0sfac,  &
 500                                          qsum
 501    real, dimension(its:ite,kts:kte)   :: denqrs1, denqrs2, denqrs3
 502    real, dimension(its:ite,kts:kte)   :: denqr1, denncr3, denqci
 503    real, dimension(its:ite)           :: delqrs1, delqrs2, delqrs3, delncr3
 504    real, dimension(its:ite)           :: delqi
 505    real, dimension(its:ite)           :: tstepsnow, tstepgraup
 506    real                               :: gfac, sfac
 507 !
 508 ! top level for computing of cloud microphysical process
 509 !
 510    integer                            :: kte_in
 511 !
 512 ! variables for optimization
 513 !
 514    real,    dimension(its:ite)        :: tvec1
 515    integer, dimension(its:ite)        :: mnstep, numndt
 516    integer, dimension(its:ite)        :: mstep, numdt
 517    logical, dimension(its:ite)        :: flgcld
 518    real                               :: temp
 519    real                               :: vt2ave
 520    real                               :: holdc, holdci
 521    real                               :: cpmcal, xlcal, lamdac, diffus,        &
 522                                          viscos, xka, venfac, conden, diffac,  &
 523                                          x, y, z, a, b, c, d, e,               &
 524                                          ndt, qdt, holdrr, holdrs, holdrg,     &
 525                                          supcol, supcolt, pvt, coeres, supsat, &
 526                                          dtcld, xmi, eacrs, satdt, qimax,      &
 527                                          diameter, xni0, roqi0,                &
 528                                          fallsum, fallsum_qsi, fallsum_qg,     &
 529                                          vt2i, vt2r, vt2s, vt2g, acrfac,       &
 530                                          egs, egi, coecol,                     &
 531                                          xlwork2, factor, source, value,       &
 532                                          nfrzdtr, nfrzdtc,                     &
 533                                          taucon, lencon, lenconcr,             &
 534                                          xlf, pfrzdtc, pfrzdtr, supice,        &
 535                                          alpha2, delta2, delta3
 536 !
 537    integer                            :: i, j, k, mstepmax,                    &
 538                                          iprt, latd, lond, loop, loops, ifsat, &
 539                                          n, idim, kdim
 540 !
 541 ! Temporaries used for inlining fpvs function
 542 !
 543    real                               :: dldti, xb, xai, tr, xbi, xa, hvap,    &
 544                                          cvap, hsub, dldt, ttp
 545 !-------------------------------------------------------------------------------
 546 !
 547 ! compute internal functions
 548 !
 549    cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
 550    xlcal(x) = xlv0-xlv1*(x-t0c)
 551 !
 552 ! size distributions: (x=mixing ratio, y=air density):
 553 ! valid for mixing ratio > 1.e-9 kg/kg.
 554 !
 555 ! optimizatin : A**B => exp(log(A)*(B))
 556 !
 557    lamdac(x,y,z)= exp(log(((pidnc*z)/(x*y)))*((.33333333)))
 558 !
 559 ! diffus: diffusion coefficient of the water vapor
 560 ! viscos: kinematic viscosity(m2s-1)
 561 !
 562    diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y   ! 8.794e-5*x**1.81/y
 563    viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y  ! 1.496e-6*x**1.5/(x+120.)/y
 564    xka(x,y)    = 1.414e3*viscos(x,y)*y
 565    diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
 566    venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333)))            &
 567                    /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
 568    conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
 569 !
 570    idim = ite-its+1
 571    kdim = kte-kts+1
 572 !
 573    kte_in = kte
 574 !
 575 !  paddint 0 for negative values generated by dynamics
 576 !
 577    do k = kts, kte
 578      do i = its, ite
 579        qci(i,k,1) = max(qci(i,k,1),0.0)
 580        qrs(i,k,1) = max(qrs(i,k,1),0.0)
 581        qci(i,k,2) = max(qci(i,k,2),0.0)
 582        qrs(i,k,2) = max(qrs(i,k,2),0.0)
 583        qrs(i,k,3) = max(qrs(i,k,3),0.0)
 584        ncr(i,k,1) = min(max(ncr(i,k,1),1.e8),2.e10)
 585        ncr(i,k,2) = max(ncr(i,k,2),0.0)
 586        ncr(i,k,3) = max(ncr(i,k,3),0.0)
 587      enddo
 588    enddo
 589 ! imsi
 590    do k = kts,kte
 591      do i = its,ite
 592        dend(i,k) = den(i,k)
 593      enddo
 594    enddo
 595 !
 596 !  latent heat for phase changes and heat capacity. neglect the
 597 !  changes during microphysical process calculation
 598 !  emanuel(1994)
 599 !
 600    do k = kts,kte_in
 601      do i = its,ite
 602        cpm(i,k) = cpmcal(q(i,k))
 603        xl(i,k) = xlcal(t(i,k))
 604      enddo
 605    enddo
 606 !
 607    qcr(:,:) = 0.0
 608    do i = its,ite
 609      if(slmsk(i).eq.2) then      ! water
 610        qcr(i,:) = qc0
 611      else
 612        qcr(i,:) = qc1
 613      endif
 614    enddo
 615 !
 616    do k = kts,kte
 617      do i = its,ite
 618        delz_tmp(i,k) = delz(i,k)
 619        den_tmp(i,k) = den(i,k)
 620      enddo
 621    enddo
 622 !
 623 ! initialize the surface rain, snow, graupel
 624 !
 625    do i = its,ite
 626      rainncv(i) = 0.
 627      if(present(snowncv) .and. present(snow)) snowncv(i) = 0.
 628      if(present (graupelncv) .and. present(graupel)) graupelncv(i) = 0.
 629      sr(i) = 0.
 630 !
 631 ! new local array to catch step snow and graupel
 632 !
 633      tstepsnow(i) = 0.
 634      tstepgraup(i) = 0.
 635    enddo
 636 !
 637 ! compute the minor time steps.
 638 !
 639    loops = max(nint(delt/dtcldcr),1)
 640    dtcld = delt/loops
 641    if(delt.le.dtcldcr) dtcld = delt
 642 !
 643    do loop = 1,loops
 644 !
 645 ! initialize the large scale variables
 646 !
 647      do i = its,ite
 648        mstep(i) = 1
 649        mnstep(i) = 1
 650        flgcld(i) = .true.
 651      enddo
 652 !
 653      do k = kts, kte
 654        call VREC( tvec1(its), den(its,k), ite-its+1)
 655        do i = its, ite
 656          tvec1(i) = tvec1(i)*den0
 657        enddo
 658        call VSQRT( denfac(its,k), tvec1(its), ite-its+1)
 659      enddo
 660 !
 661 ! Inline expansion for fpvs
 662 !   qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
 663 !   qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
 664 !
 665      hsub = xls
 666      hvap = xlv0
 667      cvap = cpv
 668      ttp = t0c+0.01
 669      dldt = cvap-cliq
 670      xa = -dldt/rv
 671      xb = xa+hvap/(rv*ttp)
 672      dldti = cvap-cice
 673      xai = -dldti/rv
 674      xbi = xai+hsub/(rv*ttp)
 675 !
 676      do k = kts,kte_in
 677        do i = its,ite
 678          tr = ttp/t(i,k)
 679          qs(i,k,1) = psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
 680          qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
 681          qs(i,k,1) = ep2*qs(i,k,1)/(p(i,k)-qs(i,k,1))
 682          qs(i,k,1) = max(qs(i,k,1),qmin)
 683          rh(i,k,1) = max(q(i,k)/qs(i,k,1),qmin)
 684          tr=ttp/t(i,k)
 685          if(t(i,k).lt.ttp) then
 686            qs(i,k,2) = psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
 687          else
 688            qs(i,k,2) = psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
 689          endif
 690          qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
 691          qs(i,k,2) = ep2*qs(i,k,2)/(p(i,k)-qs(i,k,2))
 692          qs(i,k,2) = max(qs(i,k,2),qmin)
 693          rh(i,k,2) = max(q(i,k)/qs(i,k,2),qmin)
 694        enddo
 695      enddo
 696 !
 697 !----------------------------------------------------------------
 698 !     initialize the variables for microphysical physics
 699 !
 700 !
 701       do k = kts, kte
 702         do i = its, ite
 703           prevp(i,k) = 0.
 704           psdep(i,k) = 0.
 705           pgdep(i,k) = 0.
 706           praut(i,k) = 0.
 707           psaut(i,k) = 0.
 708           pgaut(i,k) = 0.
 709           pracw(i,k) = 0.
 710           praci(i,k) = 0.
 711           piacr(i,k) = 0.
 712           psaci(i,k) = 0.
 713           psacw(i,k) = 0.
 714           pracs(i,k) = 0.
 715           psacr(i,k) = 0.
 716           pgacw(i,k) = 0.
 717           paacw(i,k) = 0.
 718           pgaci(i,k) = 0.
 719           pgacr(i,k) = 0.
 720           pgacs(i,k) = 0.
 721           pigen(i,k) = 0.
 722           pidep(i,k) = 0.
 723           pcond(i,k) = 0.
 724           psmlt(i,k) = 0.
 725           pgmlt(i,k) = 0.
 726           pseml(i,k) = 0.
 727           pgeml(i,k) = 0.
 728           psevp(i,k) = 0.
 729           pgevp(i,k) = 0.
 730           pcact(i,k) = 0.
 731           falk(i,k,1) = 0.
 732           falk(i,k,2) = 0.
 733           falk(i,k,3) = 0.
 734           fall(i,k,1) = 0.
 735           fall(i,k,2) = 0.
 736           fall(i,k,3) = 0.
 737           fallc(i,k) = 0.
 738           falkc(i,k) = 0.
 739           falln(i,k) =0.
 740           falkn(i,k) =0.
 741           xni(i,k) = 1.e3
 742           nsacw(i,k) = 0.
 743           ngacw(i,k) = 0.
 744           naacw(i,k) = 0.
 745           niacr(i,k) = 0.
 746           nsacr(i,k) = 0.
 747           ngacr(i,k) = 0.
 748           nseml(i,k) = 0.
 749           ngeml(i,k) = 0.
 750           nracw(i,k) = 0.
 751           nccol(i,k) = 0.
 752           nrcol(i,k) = 0.
 753           ncact(i,k) = 0.
 754           nraut(i,k) = 0.
 755           ncevp(i,k) = 0.
 756         enddo
 757       enddo
 758 !
 759       do k = kts, kte
 760         do i = its, ite
 761           if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin ) then
 762             rslopec(i,k) = rslopecmax
 763             rslopec2(i,k) = rslopec2max
 764             rslopec3(i,k) = rslopec3max
 765           else
 766             rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
 767             rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
 768             rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
 769           endif
 770 !-------------------------------------------------------------
 771 ! Ni: ice crystal number concentraiton   [HDC 5c]
 772 !-------------------------------------------------------------
 773           temp = (den(i,k)*max(qci(i,k,2),qmin))
 774           temp = sqrt(sqrt(temp*temp*temp))
 775           xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
 776         enddo
 777       enddo
 778 !----------------------------------------------------------------
 779 !     compute the fallout term:
 780 !     first, vertical terminal velosity for minor loops
 781 !----------------------------------------------------------------
 782       do k = kts, kte
 783         do i = its, ite
 784           qrs_tmp(i,k,1) = qrs(i,k,1)
 785           qrs_tmp(i,k,2) = qrs(i,k,2)
 786           qrs_tmp(i,k,3) = qrs(i,k,3)
 787           ncr_tmp(i,k) = ncr(i,k,3)
 788         enddo
 789       enddo
 790       call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
 791                      rslope3,work1,workn,its,ite,kts,kte)
 792 !
 793 ! vt update for qr and nr
 794 !
 795       mstepmax = 1
 796       numdt = 1
 797       do k = kte, kts, -1
 798         do i = its, ite
 799           work1(i,k,1) = work1(i,k,1)/delz(i,k)
 800           workn(i,k) = workn(i,k)/delz(i,k)
 801           numdt(i) = max(nint(max(work1(i,k,1),workn(i,k))*dtcld+.5),1)
 802           if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
 803         enddo
 804       enddo
 805       do i = its, ite
 806         if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
 807       enddo
 808 !
 809       do n = 1, mstepmax
 810         k = kte
 811         do i = its, ite
 812          if(n.le.mstep(i)) then
 813            falk(i,k,1) = dend(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
 814            falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
 815            fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
 816            falln(i,k) = falln(i,k)+falkn(i,k)
 817            qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/dend(i,k),0.)
 818            ncr(i,k,3) = max(ncr(i,k,3)-falkn(i,k)*dtcld,0.)
 819          endif
 820        enddo
 821 !
 822        do k = kte_in-1,kts,-1
 823          do i = its,ite
 824            if(n.le.mstep(i)) then
 825              falk(i,k,1) = dend(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
 826              falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
 827              fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
 828              falln(i,k) = falln(i,k)+falkn(i,k)
 829              dqr(i,k) = min(falk(i,k,1)*dtcld/dend(i,k),qrs(i,k,1))
 830              dqr(i,k+1) = min(falk(i,k+1,1)*delz(i,k+1)/delz(i,k)              &
 831                           *dtcld/dend(i,k),qrs(i,k+1,1))
 832              dnr(i,k) = min(falkn(i,k)*dtcld,ncr(i,k,3))
 833              dnr(i,k+1) = min(falkn(i,k+1)*delz(i,k+1)/delz(i,k)*dtcld,        &
 834                           ncr(i,k+1,3))
 835              qrs(i,k,1) = max(qrs(i,k,1)-dqr(i,k)+dqr(i,k+1),0.)
 836              ncr(i,k,3) = max(ncr(i,k,3)-dnr(i,k)+dnr(i,k+1),0.)
 837             endif
 838           enddo
 839         enddo
 840         do k = kts, kte
 841           do i = its, ite
 842             qrs_tmp(i,k,1) = qrs(i,k,1)
 843             ncr_tmp(i,k) = ncr(i,k,3)
 844           enddo
 845         enddo
 846         call slope_rain(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
 847                      rslope3,work1,workn,its,ite,kts,kte)
 848         do k = kte, kts, -1
 849           do i = its, ite
 850             work1(i,k,1) = work1(i,k,1)/delz(i,k)
 851             workn(i,k) = workn(i,k)/delz(i,k)
 852           enddo
 853         enddo
 854       enddo
 855 !
 856 ! for semi
 857 !
 858       do k = kte, kts, -1
 859         do i = its, ite
 860           qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
 861           if(qsum(i,k) .gt. 1.e-15 ) then
 862             worka(i,k) = (work1(i,k,2)*qrs(i,k,2) + work1(i,k,3)*qrs(i,k,3)) &
 863                       /qsum(i,k)
 864           else
 865             worka(i,k) = 0.
 866           endif
 867           denqrs2(i,k) = den(i,k)*qrs(i,k,2)
 868           denqrs3(i,k) = den(i,k)*qrs(i,k,3)
 869         enddo
 870       enddo
 871       call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka,         &
 872                            denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1)
 873       do k = kts, kte
 874         do i = its, ite
 875           qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
 876           qrs(i,k,3) = max(denqrs3(i,k)/den(i,k),0.)
 877           fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k)
 878           fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k)
 879         enddo
 880       enddo
 881       do i = its, ite
 882         fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
 883         fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld
 884       enddo
 885       do k = kts, kte
 886         do i = its, ite
 887           qrs_tmp(i,k,1) = qrs(i,k,1)
 888           qrs_tmp(i,k,2) = qrs(i,k,2)
 889           qrs_tmp(i,k,3) = qrs(i,k,3)
 890           ncr_tmp(i,k) = ncr(i,k,3)
 891         enddo
 892       enddo
 893       call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
 894                      rslope3,work1,workn,its,ite,kts,kte)
 895 !
 896       do k = kte, kts, -1
 897         do i = its, ite
 898           supcol = t0c-t(i,k)
 899           n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
 900           if(t(i,k).gt.t0c) then
 901 !---------------------------------------------------------------
 902 ! psmlt: melting of snow [HL A33] [RH83 A25]
 903 !       (T>T0: QS->QR)
 904 !---------------------------------------------------------------
 905             xlf = xlf0
 906             work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
 907             if(qrs(i,k,2).gt.0.) then
 908               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
 909               psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2.       &
 910                          *n0sfac(i,k)*(precs1*rslope2(i,k,2)                 &
 911                          +precs2*work2(i,k)*coeres)/den(i,k)
 912               psmlt(i,k) = min(max(psmlt(i,k)*dtcld,-qrs(i,k,2)),0.)
 913 !-------------------------------------------------------------------
 914 ! nsmlt: melting of snow [LH A27]
 915 !       (T>T0: ->NR)
 916 !-------------------------------------------------------------------
 917               if(qrs(i,k,2).gt.qcrmin) then
 918                 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
 919                 ncr(i,k,3) = ncr(i,k,3) - sfac*psmlt(i,k)
 920               endif
 921 ! error correction based on Lei et al., (JGR, 2020)
 922               qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
 923               qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
 924               t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
 925             endif
 926 !---------------------------------------------------------------
 927 ! pgmlt: melting of graupel [HL A23]  [LFO 47]
 928 !       (T>T0: QG->QR)
 929 !---------------------------------------------------------------
 930             if(qrs(i,k,3).gt.0.) then
 931               coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
 932               pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*(precg1     &
 933                            *rslope2(i,k,3) + precg2*work2(i,k)*coeres)       &
 934                            /den(i,k)
 935               pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld,-qrs(i,k,3)),0.)
 936 !-------------------------------------------------------------------
 937 ! ngmlt: melting of graupel [LH A28]
 938 !       (T>T0: ->NR)
 939 !-------------------------------------------------------------------
 940               if(qrs(i,k,3).gt.qcrmin) then
 941                 gfac = rslope(i,k,3)*n0g/qrs(i,k,3)
 942                 ncr(i,k,3) = ncr(i,k,3) - gfac*pgmlt(i,k)
 943               endif
 944 ! error correction based on Lei et al., (JGR, 2020)
 945               qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
 946               qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
 947               t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
 948             endif
 949           endif
 950         enddo
 951       enddo
 952 !---------------------------------------------------------------
 953 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
 954 !---------------------------------------------------------------
 955       do k = kte, kts, -1
 956         do i = its, ite
 957           if(qci(i,k,2).le.0.) then
 958             work1c(i,k) = 0.
 959           else
 960             xmi = den(i,k)*qci(i,k,2)/xni(i,k)
 961             diameter  = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
 962             work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
 963           endif
 964         enddo
 965       enddo
 966 !
 967 !  forward semi-laglangian scheme (JH), PCM (piecewise constant),  (linear)
 968 !
 969       do k = kte, kts, -1
 970         do i = its, ite
 971           denqci(i,k) = den(i,k)*qci(i,k,2)
 972         enddo
 973       enddo
 974       call nislfv_rain_plmr(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci,denqci,  &
 975                            delqi,dtcld,1,0,0)
 976       do k = kts, kte
 977         do i = its, ite
 978           qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
 979         enddo
 980       enddo
 981       do i = its, ite
 982         fallc(i,1) = delqi(i)/delz(i,1)/dtcld
 983       enddo
 984 !
 985 !----------------------------------------------------------------
 986 !      rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
 987 !
 988       do i = its, ite
 989         fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
 990         fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
 991         fallsum_qg = fall(i,kts,3)
 992         if(fallsum.gt.0.) then
 993           rainncv(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rainncv(i)
 994           rain(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rain(i)
 995         endif
 996         If(fallsum_qsi.gt.0.) then
 997           tstepsnow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + tstepsnow(i)
 998           IF( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
 999             snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snowncv(i)
1000             snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
1001           ENDIF
1002         ENDIF
1003         IF(fallsum_qg.gt.0.) then
1004           tstepgraup(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000.              &
1005                             + tstepgraup(i)
1006           IF( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN
1007             graupelncv(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000.            &
1008                             + graupelncv(i)
1009             graupel(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. + graupel(i)
1010           ENDIF
1011         ENDIF
1012         IF ( PRESENT (snowncv)) THEN
1013           if(fallsum.gt.0.)sr(i)=(snowncv(i) + graupelncv(i))/(rainncv(i)+1.e-12)
1014         ELSE
1015           if(fallsum.gt.0.)sr(i)=(tstepsnow(i) + tstepgraup(i))/(rainncv(i)+1.e-12)
1016         ENDIF
1017       enddo
1018 !
1019 !---------------------------------------------------------------
1020 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
1021 !       (T>T0: QI->QC)
1022 !---------------------------------------------------------------
1023       do k = kts, kte
1024         do i = its, ite
1025           supcol = t0c-t(i,k)
1026           xlf = xls-xl(i,k)
1027           if(supcol.lt.0.) xlf = xlf0
1028           if(supcol.lt.0 .and. qci(i,k,2).gt.0.) then
1029             qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
1030 !---------------------------------------------------------------
1031 ! nimlt: instantaneous melting of cloud ice  [LH A18]
1032 !        (T>T0: ->NC)
1033 !--------------------------------------------------------------
1034             ncr(i,k,2) = ncr(i,k,2) + xni(i,k)
1035             t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
1036             qci(i,k,2) = 0.
1037           endif
1038 !---------------------------------------------------------------
1039 ! pihmf: homogeneous  of cloud water below -40c [HL A45]
1040 !        (T<-40C: QC->QI)
1041 !---------------------------------------------------------------
1042           if(supcol.gt.40. .and. qci(i,k,1).gt.0.) then
1043             qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
1044 !---------------------------------------------------------------
1045 ! nihmf: homogeneous  of cloud water below -40c [LH A17]
1046 !        (T<-40C: NC->)
1047 !---------------------------------------------------------------
1048             if(ncr(i,k,2).gt.0.) ncr(i,k,2) = 0.
1049             t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
1050             qci(i,k,1) = 0.
1051           endif
1052 !---------------------------------------------------------------
1053 ! pihtf: heterogeneous  of cloud water [HL A44]
1054 !        (T0>T>-40C: QC->QI)
1055 !---------------------------------------------------------------
1056           if(supcol.gt.0. .and. qci(i,k,1).gt.qmin) then
1057             supcolt=min(supcol,70.)
1058             pfrzdtc = min(pi*pi*pfrz1*(exp(pfrz2*supcolt)-1.)*denr/den(i,k)    &
1059                      *ncr(i,k,2)*rslopec3(i,k)*rslopec3(i,k)/18.*dtcld         &
1060                      ,qci(i,k,1))
1061 !---------------------------------------------------------------
1062 ! nihtf: heterogeneous  of cloud water [LH A16]
1063 !         (T0>T>-40C: NC->)
1064 !---------------------------------------------------------------
1065             if(ncr(i,k,2).gt.ncmin) then
1066               nfrzdtc = min(pi*pfrz1*(exp(pfrz2*supcolt)-1.)*ncr(i,k,2)        &
1067                       *rslopec3(i,k)/6.*dtcld,ncr(i,k,2))
1068               ncr(i,k,2) = ncr(i,k,2) - nfrzdtc
1069             endif
1070             qci(i,k,2) = qci(i,k,2) + pfrzdtc
1071             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
1072             qci(i,k,1) = qci(i,k,1)-pfrzdtc
1073           endif
1074 !---------------------------------------------------------------
1075 ! pgfrz:  freezing of rain water [HL A20] [LFO 45]
1076 !        (T<T0, QR->QG)
1077 !---------------------------------------------------------------
1078           if(supcol.gt.0. .and. qrs(i,k,1).gt.0.) then
1079             supcolt=min(supcol,70.)
1080             pfrzdtr = min(140.*(pi*pi)*pfrz1*ncr(i,k,3)*denr/den(i,k)          &
1081                   *(exp(pfrz2*supcolt)-1.)*rslope3(i,k,1)*rslope3(i,k,1)       &
1082                   *dtcld,qrs(i,k,1))
1083 !---------------------------------------------------------------
1084 ! ngfrz: freezing of rain water [LH A26]
1085 !        (T<T0, NR-> )
1086 !---------------------------------------------------------------
1087             if(ncr(i,k,3).gt.nrmin) then
1088               nfrzdtr = min(4.*pi*pfrz1*ncr(i,k,3)*(exp(pfrz2*supcolt)-1.)     &
1089                        *rslope3(i,k,1)*dtcld, ncr(i,k,3))
1090               ncr(i,k,3) = ncr(i,k,3) - nfrzdtr
1091             endif
1092             qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
1093             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
1094             qrs(i,k,1) = qrs(i,k,1) - pfrzdtr
1095           endif
1096         enddo
1097       enddo
1098 !
1099       do k = kts, kte
1100         do i = its, ite
1101           ncr(i,k,2) = max(ncr(i,k,2),0.0)
1102           ncr(i,k,3) = max(ncr(i,k,3),0.0)
1103         enddo
1104       enddo
1105 !
1106 !----------------------------------------------------------------
1107 !     update the slope parameters for microphysics computation
1108 !
1109       do k = kts, kte
1110         do i = its, ite
1111           qrs_tmp(i,k,1) = qrs(i,k,1)
1112           qrs_tmp(i,k,2) = qrs(i,k,2)
1113           qrs_tmp(i,k,3) = qrs(i,k,3)
1114           ncr_tmp(i,k) = ncr(i,k,3)
1115         enddo
1116       enddo
1117       call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
1118                      rslope3,work1,workn,its,ite,kts,kte)
1119       do k = kts, kte
1120         do i = its, ite
1121 !-----------------------------------------------------------------
1122 ! compute the mean-volume drop diameter                  [LH A10]
1123 ! for raindrop distribution
1124 !-----------------------------------------------------------------
1125           avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
1126 !
1127           if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin) then
1128             rslopec(i,k) = rslopecmax
1129             rslopec2(i,k) = rslopec2max
1130             rslopec3(i,k) = rslopec3max
1131           else
1132             rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
1133             rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
1134             rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
1135           endif
1136 !--------------------------------------------------------------------
1137 ! compute the mean-volume drop diameter                   [LH A7]
1138 ! for cloud-droplet distribution
1139 !--------------------------------------------------------------------
1140           avedia(i,k,1) = rslopec(i,k)
1141         enddo
1142       enddo
1143 !
1144       do k = kts, kte
1145         do i = its, ite
1146           work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
1147           work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
1148           work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
1149         enddo
1150       enddo
1151 !
1152 !===============================================================
1153 !
1154 ! warm rain processes
1155 !
1156 ! - follows the double-moment processes in Lim and Hong
1157 !
1158 !===============================================================
1159 !
1160       do k = kts, kte
1161         do i = its, ite
1162           supsat = max(q(i,k),qmin)-qs(i,k,1)
1163           satdt = supsat/dtcld
1164 !---------------------------------------------------------------
1165 ! praut: auto conversion rate from cloud to rain  [LH 9] [CP 17]
1166 !        (QC->QR)
1167 !--------------------------------------------------------------
1168          lencon  = 2.7e-2*den(i,k)*qci(i,k,1)*(1.e20/16.*rslopec2(i,k)         &
1169                   *rslopec2(i,k)-0.4)
1170          lenconcr = max(1.2*lencon, qcrmin)
1171          if(qci(i,k,1).gt.qcr(i,k).and.ncr(i,k,2).gt.ncmin) then
1172            praut(i,k) = qck1*qci(i,k,1)**(7./3.)*ncr(i,k,2)**(-1./3.)
1173            praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
1174 !----------------------------------------------------------------------
1175 ! nraut: auto conversion rate from cloud to rain [LH A6] [CP 18 & 19]
1176 !        (NC->NR)
1177 !------------------------------------------------------------------------
1178            nraut(i,k) = 3.5e9*den(i,k)*praut(i,k)
1179            if(qrs(i,k,1).gt.lenconcr)                                           &
1180            nraut(i,k) = ncr(i,k,3)/qrs(i,k,1)*praut(i,k)
1181            nraut(i,k) = min(nraut(i,k),ncr(i,k,2)/dtcld)
1182          endif
1183 !---------------------------------------------------------------
1184 ! pracw: accretion of cloud water by rain     [LH 10] [CP 22 & 23]
1185 !        (QC->QR)
1186 ! nracw: accretion of cloud water by rain     [LH A9]
1187 !        (NC->)
1188 !---------------------------------------------------------------
1189           if(qrs(i,k,1).ge.lenconcr) then
1190             if(avedia(i,k,2).ge.di100) then
1191               nracw(i,k) = min(ncrk1*ncr(i,k,2)*ncr(i,k,3)*(rslopec3(i,k)      &
1192                          + 24.*rslope3(i,k,1)),ncr(i,k,2)/dtcld)
1193               pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk1*ncr(i,k,2)          &
1194                          *ncr(i,k,3)*rslopec3(i,k)*(2.*rslopec3(i,k)           &
1195                          + 24.*rslope3(i,k,1)),qci(i,k,1)/dtcld)
1196             else
1197               nracw(i,k) = min(ncrk2*ncr(i,k,2)*ncr(i,k,3)*(2.*rslopec3(i,k)   &
1198                          *rslopec3(i,k)+5040.*rslope3(i,k,1)                   &
1199                          *rslope3(i,k,1)),ncr(i,k,2)/dtcld)
1200               pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk2*ncr(i,k,2)          &
1201                          *ncr(i,k,3)*rslopec3(i,k)*(6.*rslopec3(i,k)           &
1202                          *rslopec3(i,k)+5040.*rslope3(i,k,1)*rslope3(i,k,1))   &
1203                          ,qci(i,k,1)/dtcld)
1204             endif
1205           endif
1206 !----------------------------------------------------------------
1207 ! nccol: self collection of cloud water             [LH A8] [CP 24 & 25]
1208 !        (NC->)
1209 !----------------------------------------------------------------
1210           if(avedia(i,k,1).ge.di100) then
1211             nccol(i,k) = ncrk1*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)
1212           else
1213             nccol(i,k) = 2.*ncrk2*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)        &
1214                          *rslopec3(i,k)
1215           endif
1216 !----------------------------------------------------------------
1217 ! nrcol: self collection of rain-drops and break-up [LH A21] [CP 24 & 25]
1218 !        (NR->)
1219 !----------------------------------------------------------------
1220           if(qrs(i,k,1).ge.lenconcr) then
1221             if(avedia(i,k,2).lt.di100) then
1222               nrcol(i,k) = 5040.*ncrk2*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1)    &
1223                           *rslope3(i,k,1)
1224             elseif(avedia(i,k,2).ge.di100 .and. avedia(i,k,2).lt.di600) then
1225               nrcol(i,k) = 24.*ncrk1*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1)
1226             elseif(avedia(i,k,2).ge.di600 .and. avedia(i,k,2).lt.di2000) then
1227               coecol = -2.5e3*(avedia(i,k,2)-di600)
1228               nrcol(i,k) = 24.*exp(coecol)*ncrk1*ncr(i,k,3)*ncr(i,k,3)         &
1229                          *rslope3(i,k,1)
1230             else
1231               nrcol(i,k) = 0.
1232             endif
1233           endif
1234 !---------------------------------------------------------------
1235 ! prevp: evaporation/condensation rate of rain   [HL A41]
1236 !        (QV->QR or QR->QV)
1237 !---------------------------------------------------------------
1238           if(qrs(i,k,1).gt.0.) then
1239             coeres = rslope(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
1240             prevp(i,k) = (rh(i,k,1)-1.)*ncr(i,k,3)*(precr1*rslope(i,k,1)       &
1241                          + precr2*work2(i,k)*coeres)/work1(i,k,1)
1242             if(prevp(i,k).lt.0.) then
1243               prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
1244               prevp(i,k) = max(prevp(i,k),satdt/2)
1245 !----------------------------------------------------------------
1246 ! Nrevp: evaporation/condensation rate of rain   [LH A14]
1247 !        (NR->NCCN)
1248 !----------------------------------------------------------------
1249               if(prevp(i,k).eq.-qrs(i,k,1)/dtcld) then
1250                 ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,3)
1251                 ncr(i,k,3) = 0.
1252               endif
1253             else
1254 !
1255               prevp(i,k) = min(prevp(i,k),satdt/2)
1256             endif
1257           endif
1258         enddo
1259       enddo
1260 !
1261 !===============================================================
1262 !
1263 ! cold rain processes
1264 !
1265 ! - follows the revised ice microphysics processes in HDC
1266 ! - the processes same as in RH83 and RH84  and LFO behave
1267 !   following ice crystal hapits defined in HDC, inclduing
1268 !   intercept parameter for snow (n0s), ice crystal number
1269 !   concentration (ni), ice nuclei number concentration
1270 !   (n0i), ice diameter (d)
1271 !
1272 !===============================================================
1273 !
1274       do k = kts, kte
1275         do i = its, ite
1276           supcol = t0c-t(i,k)
1277           n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1278           supsat = max(q(i,k),qmin)-qs(i,k,2)
1279           satdt = supsat/dtcld
1280           ifsat = 0
1281 !-------------------------------------------------------------
1282 ! Ni: ice crystal number concentraiton   [HDC 5c]
1283 !-------------------------------------------------------------
1284 !         xni(i,k) = min(max(5.38e7*(den(i,k)                                  &
1285 !                      *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
1286           temp = (den(i,k)*max(qci(i,k,2),qmin))
1287           temp = sqrt(sqrt(temp*temp*temp))
1288           xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1289           eacrs = exp(0.07*(-supcol))
1290 !
1291           xmi = den(i,k)*qci(i,k,2)/xni(i,k)
1292           diameter  = min(dicon * sqrt(xmi),dimax)
1293           vt2i = 1.49e4*diameter**1.31
1294           vt2r=pvtr*rslopeb(i,k,1)*denfac(i,k)
1295           vt2s=pvts*rslopeb(i,k,2)*denfac(i,k)
1296           vt2g=pvtg*rslopeb(i,k,3)*denfac(i,k)
1297           qsum(i,k) = max((qrs(i,k,2)+qrs(i,k,3)),1.e-15)
1298           if(qsum(i,k) .gt. 1.e-15) then
1299             vt2ave=(vt2s*qrs(i,k,2)+vt2g*qrs(i,k,3))/(qsum(i,k))
1300           else
1301             vt2ave=0.
1302           endif
1303           if(supcol.gt.0. .and. qci(i,k,2).gt.qmin) then
1304             if(qrs(i,k,1).gt.qcrmin) then
1305 !-------------------------------------------------------------
1306 ! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
1307 !        (T<T0: QI->QR)
1308 !-------------------------------------------------------------
1309               acrfac = 6.*rslope2(i,k,1)+4.*diameter*rslope(i,k,1) + diameter**2
1310               praci(i,k) = pi*qci(i,k,2)*ncr(i,k,3)*abs(vt2r-vt2i)*acrfac/4.
1311               ! reduce collection efficiency (suggested by B. Wilt)
1312               praci(i,k) = praci(i,k)*min(max(0.0,qrs(i,k,1)/qci(i,k,2)),1.)**2
1313               praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
1314 !-------------------------------------------------------------
1315 ! piacr: Accretion of rain by cloud ice [HL A19] [LFO 26]
1316 !        (T<T0: QR->QS or QR->QG)
1317 !-------------------------------------------------------------
1318               piacr(i,k) = pi*pi*avtr*ncr(i,k,3)*denr*xni(i,k)*denfac(i,k)     &
1319                           *g7pbr*rslope3(i,k,1)*rslope2(i,k,1)*rslopeb(i,k,1)  &
1320                           /24./den(i,k)
1321               ! reduce collection efficiency (suggested by B. Wilt)
1322               piacr(i,k) = piacr(i,k)*min(max(0.0,qci(i,k,2)/qrs(i,k,1)),1.)**2
1323               piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
1324             endif
1325 !-------------------------------------------------------------
1326 ! niacr: Accretion of rain by cloud ice  [LH A25]
1327 !        (T<T0: NR->)
1328 !-------------------------------------------------------------
1329             if(ncr(i,k,3).gt.nrmin) then
1330               niacr(i,k) = pi*avtr*ncr(i,k,3)*xni(i,k)*denfac(i,k)*g4pbr       &
1331                           *rslope2(i,k,1)*rslopeb(i,k,1)/4.
1332               ! reduce collection efficiency (suggested by B. Wilt)
1333               niacr(i,k) = niacr(i,k)*min(max(0.0,qci(i,k,2)/qrs(i,k,1)),1.)**2
1334               niacr(i,k) = min(niacr(i,k),ncr(i,k,3)/dtcld)
1335             endif
1336 !-------------------------------------------------------------
1337 ! psaci: Accretion of cloud ice by snow [HDC 10]
1338 !        (T<T0: QI->QS)
1339 !-------------------------------------------------------------
1340             if(qrs(i,k,2).gt.qcrmin) then
1341               acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2)            &
1342                       + diameter**2*rslope(i,k,2)
1343               psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k)                 &
1344                           *abs(vt2ave-vt2i)*acrfac/4.
1345               psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
1346             endif
1347 !-------------------------------------------------------------
1348 ! pgaci: Accretion of cloud ice by graupel [HL A17] [LFO 41]
1349 !        (T<T0: QI->QG)
1350 !-------------------------------------------------------------
1351             if(qrs(i,k,3).gt.qcrmin) then
1352               egi = exp(0.07*(-supcol))
1353               acrfac = 2.*rslope3(i,k,3)+2.*diameter*rslope2(i,k,3)            &
1354                       + diameter**2*rslope(i,k,3)
1355               pgaci(i,k) = pi*egi*qci(i,k,2)*n0g*abs(vt2ave-vt2i)*acrfac/4.
1356               pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
1357             endif
1358           endif
1359 !-------------------------------------------------------------
1360 ! psacw: Accretion of cloud water by snow  [HL A7] [LFO 24]
1361 !        (T<T0: QC->QS, and T>=T0: QC->QR)
1362 !-------------------------------------------------------------
1363           if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,1).gt.qmin) then
1364             psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2)   &
1365               ! reduce collection efficiency (suggested by B. Wilt)
1366                         *min(max(0.0,qrs(i,k,2)/qci(i,k,1)),1.)**2             &
1367                         *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
1368           endif
1369 !-------------------------------------------------------------
1370 ! nsacw: Accretion of cloud water by snow  [LH A12]
1371 !        (NC ->)
1372 !-------------------------------------------------------------
1373          if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
1374            nsacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2)    &
1375               ! reduce collection efficiency (suggested by B. Wilt)
1376                        *min(max(0.0,qrs(i,k,2)/qci(i,k,1)),1.)**2              &
1377                        *ncr(i,k,2)*denfac(i,k),ncr(i,k,2)/dtcld)
1378          endif
1379 !-------------------------------------------------------------
1380 ! pgacw: Accretion of cloud water by graupel [HL A6] [LFO 40]
1381 !        (T<T0: QC->QG, and T>=T0: QC->QR)
1382 !-------------------------------------------------------------
1383           if(qrs(i,k,3).gt.qcrmin .and. qci(i,k,1).gt.qmin) then
1384             pgacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3)*qci(i,k,1)    &
1385               ! reduce collection efficiency (suggested by B. Wilt)
1386                         *min(max(0.0,qrs(i,k,3)/qci(i,k,1)),1.)**2             &
1387                         *denfac(i,k),qci(i,k,1)/dtcld)
1388           endif
1389 !-------------------------------------------------------------
1390 ! ngacw: Accretion of cloud water by graupel [LH A13]
1391 !        (NC->
1392 !-------------------------------------------------------------
1393           if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
1394             ngacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3)*ncr(i,k,2)    &
1395               ! reduce collection efficiency (suggested by B. Wilt)
1396                         *min(max(0.0,qrs(i,k,3)/qci(i,k,1)),1.)**2             &
1397                         *denfac(i,k),ncr(i,k,2)/dtcld)
1398           endif
1399 !-------------------------------------------------------------
1400 ! paacw: Accretion of cloud water by averaged snow/graupel
1401 !        (T<T0: QC->QG or QS, and T>=T0: QC->QR)
1402 !-------------------------------------------------------------
1403           if(qsum(i,k) .gt. 1.e-15 ) then
1404             paacw(i,k) = (qrs(i,k,2)*psacw(i,k)+qrs(i,k,3)*pgacw(i,k))/(qsum(i,k))
1405 !-------------------------------------------------------------
1406 ! naacw: Accretion of cloud water by averaged snow/graupel
1407 !        (Nc->)
1408 !-------------------------------------------------------------
1409             naacw(i,k) = (qrs(i,k,2)*nsacw(i,k)+qrs(i,k,3)*ngacw(i,k))/(qsum(i,k))
1410           endif
1411 !-------------------------------------------------------------
1412 ! pracs: Accretion of snow by rain [HL A11] [LFO 27]
1413 !         (T<T0: QS->QG)
1414 !-------------------------------------------------------------
1415           if(qrs(i,k,2).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin) then
1416             if(supcol.gt.0) then
1417               acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2)                        &
1418                       + 4.*rslope3(i,k,2)*rslope2(i,k,2)*rslope(i,k,1)         &
1419                       + 1.5*rslope2(i,k,2)*rslope2(i,k,2)*rslope2(i,k,1)
1420               pracs(i,k) = pi*pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2r-vt2ave)   &
1421                           *(dens/den(i,k))*acrfac
1422               ! reduce collection efficiency (suggested by B. Wilt)
1423               pracs(i,k) = pracs(i,k)*min(max(0.0,qrs(i,k,1)/qrs(i,k,2)),1.)**2
1424               pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
1425             endif
1426 !-------------------------------------------------------------
1427 ! psacr: Accretion of rain by snow [HL A10] [LFO 28]
1428 !         (T<T0:QR->QS or QR->QG) (T>=T0: enhance melting of snow)
1429 !-------------------------------------------------------------
1430             acrfac = 30.*rslope3(i,k,1)*rslope2(i,k,1)*rslope(i,k,2)           &
1431                      +10.*rslope2(i,k,1)*rslope2(i,k,1)*rslope2(i,k,2)         &
1432                      + 2.*rslope3(i,k,1)*rslope3(i,k,2)
1433             psacr(i,k) = pi*pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2ave-vt2r)     &
1434                         *(denr/den(i,k))*acrfac
1435               ! reduce collection efficiency (suggested by B. Wilt)
1436             psacr(i,k) = psacr(i,k)*min(max(0.0,qrs(i,k,2)/qrs(i,k,1)),1.)**2
1437             psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
1438           endif
1439           if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,3).gt.nrmin) then
1440 !-------------------------------------------------------------
1441 ! nsacr: Accretion of rain by snow  [LH A23]
1442 !       (T<T0:NR->)
1443 !-------------------------------------------------------------
1444             acrfac = 1.5*rslope2(i,k,1)*rslope(i,k,2)                          &
1445                     + 1.0*rslope(i,k,1)*rslope2(i,k,2)+.5*rslope3(i,k,2)
1446             nsacr(i,k) = pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2ave-vt2r)        &
1447                         *acrfac
1448               ! reduce collection efficiency (suggested by B. Wilt)
1449             nsacr(i,k) = nsacr(i,k)*min(max(0.0,qrs(i,k,2)/qrs(i,k,1)),1.)**2
1450             nsacr(i,k) = min(nsacr(i,k),ncr(i,k,3)/dtcld)
1451           endif
1452 !-------------------------------------------------------------
1453 ! pgacr: Accretion of rain by graupel [HL A12] [LFO 42]
1454 !         (T<T0: QR->QG) (T>=T0: enhance melting of graupel)
1455 !-------------------------------------------------------------
1456           if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin) then
1457             acrfac = 30.*rslope3(i,k,1)*rslope2(i,k,1)*rslope(i,k,3)           &
1458                     +10.*rslope2(i,k,1)*rslope2(i,k,1)*rslope2(i,k,3)          &
1459                     + 2.*rslope3(i,k,1)*rslope3(i,k,3)
1460             pgacr(i,k) = pi*pi*ncr(i,k,3)*n0g*abs(vt2ave-vt2r)*(denr/den(i,k)) &
1461                         *acrfac
1462               ! reduce collection efficiency (suggested by B. Wilt)
1463             pgacr(i,k) = pgacr(i,k)*min(max(0.0,qrs(i,k,3)/qrs(i,k,1)),1.)**2
1464             pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
1465           endif
1466 !-------------------------------------------------------------
1467 ! ngacr: Accretion of rain by graupel  [LH A24]
1468 !        (T<T0: NR->)
1469 !-------------------------------------------------------------
1470           if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,3).gt.nrmin) then
1471             acrfac = 1.5*rslope2(i,k,1)*rslope(i,k,3)                          &
1472                     + 1.0*rslope(i,k,1)*rslope2(i,k,3) + .5*rslope3(i,k,3)
1473             ngacr(i,k) = pi*ncr(i,k,3)*n0g*abs(vt2ave-vt2r)*acrfac
1474               ! reduce collection efficiency (suggested by B. Wilt)
1475             ngacr(i,k) = ngacr(i,k)*min(max(0.0,qrs(i,k,3)/qrs(i,k,1)),1.)**2
1476             ngacr(i,k) = min(ngacr(i,k),ncr(i,k,3)/dtcld)
1477           endif
1478 !
1479 !-------------------------------------------------------------
1480 ! pgacs: Accretion of snow by graupel [HL A13] [LFO 29]
1481 !        (QS->QG) : This process is eliminated in V3.0 with the
1482 !        new combined snow/graupel fall speeds
1483 !-------------------------------------------------------------
1484           if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,2).gt.qcrmin) then
1485             pgacs(i,k) = 0.
1486           endif
1487           if(supcol.le.0) then
1488             xlf = xlf0
1489 !-------------------------------------------------------------
1490 ! pseml: Enhanced melting of snow by accretion of water [HL A34]
1491 !        (T>=T0: QS->QR)
1492 !-------------------------------------------------------------
1493             if(qrs(i,k,2).gt.0.)                                               &
1494               pseml(i,k) = min(max(cliq*supcol*(paacw(i,k)+psacr(i,k))         &
1495                           /xlf,-qrs(i,k,2)/dtcld),0.)
1496 !--------------------------------------------------------------
1497 ! nseml: Enhanced melting of snow by accretion of water    [LH A29]
1498 !        (T>=T0: ->NR)
1499 !--------------------------------------------------------------
1500               if  (qrs(i,k,2).gt.qcrmin) then
1501                 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
1502                 nseml(i,k) = -sfac*pseml(i,k)
1503               endif
1504 !-------------------------------------------------------------
1505 ! pgeml: Enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22]
1506 !        (T>=T0: QG->QR)
1507 !-------------------------------------------------------------
1508             if(qrs(i,k,3).gt.0.)                                               &
1509               pgeml(i,k) = min(max(cliq*supcol*(paacw(i,k)+pgacr(i,k))/xlf     &
1510                           ,-qrs(i,k,3)/dtcld),0.)
1511 !--------------------------------------------------------------
1512 ! ngeml: Enhanced melting of graupel by accretion of water [LH A30]
1513 !         (T>=T0: -> NR)
1514 !--------------------------------------------------------------
1515               if (qrs(i,k,3).gt.qcrmin) then
1516                 gfac = rslope(i,k,3)*n0g/qrs(i,k,3)
1517                 ngeml(i,k) = -gfac*pgeml(i,k)
1518               endif
1519           endif
1520           if(supcol.gt.0) then
1521 !-------------------------------------------------------------
1522 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
1523 !       (T<T0: QV->QI or QI->QV)
1524 !-------------------------------------------------------------
1525             if(qci(i,k,2).gt.0. .and. ifsat.ne.1) then
1526               pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1527               supice = satdt-prevp(i,k)
1528               if(pidep(i,k).lt.0.) then
1529                 pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1530                 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1531               else
1532                 pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1533               endif
1534               if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1535             endif
1536 !-------------------------------------------------------------
1537 ! psdep: deposition/sublimation rate of snow [HDC 14]
1538 !        (T<T0: QV->QS or QS->QV)
1539 !-------------------------------------------------------------
1540             if(qrs(i,k,2).gt.0. .and. ifsat.ne.1) then
1541               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1542               psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2)   &
1543                            + precs2*work2(i,k)*coeres)/work1(i,k,2)
1544               supice = satdt-prevp(i,k)-pidep(i,k)
1545               if(psdep(i,k).lt.0.) then
1546                 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1547                 psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1548               else
1549                 psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1550               endif
1551               if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) ifsat = 1
1552             endif
1553 !-------------------------------------------------------------
1554 ! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46]
1555 !        (T<T0: QV->QG or QG->QV)
1556 !-------------------------------------------------------------
1557             if(qrs(i,k,3).gt.0. .and. ifsat.ne.1) then
1558               coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1559               pgdep(i,k) = (rh(i,k,2)-1.)*(precg1*rslope2(i,k,3)               &
1560                           + precg2*work2(i,k)*coeres)/work1(i,k,2)
1561               supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1562               if(pgdep(i,k).lt.0.) then
1563                 pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
1564                 pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
1565               else
1566                 pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
1567               endif
1568               if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge.          &
1569                 abs(satdt)) ifsat = 1
1570             endif
1571 !-------------------------------------------------------------
1572 ! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8]
1573 !       (T<T0: QV->QI)
1574 !-------------------------------------------------------------
1575             if(supsat.gt.0. .and. ifsat.ne.1) then
1576               supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
1577               xni0 = 1.e3*exp(0.1*supcol)
1578               roqi0 = 4.92e-11*xni0**1.33
1579               pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))/dtcld)
1580               pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1581             endif
1582 !
1583 !-------------------------------------------------------------
1584 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1585 !        (T<T0: QI->QS)
1586 !-------------------------------------------------------------
1587             if(qci(i,k,2).gt.0.) then
1588               qimax = roqimax/den(i,k)
1589               psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1590             endif
1591 !
1592 !-------------------------------------------------------------
1593 ! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37]
1594 !        (T<T0: QS->QG)
1595 !-------------------------------------------------------------
1596             if(qrs(i,k,2).gt.0.) then
1597               alpha2 = 1.e-3*exp(0.09*(-supcol))
1598               pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)),qrs(i,k,2)/dtcld)
1599             endif
1600           endif
1601 !
1602 !-------------------------------------------------------------
1603 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1604 !       (T>=T0: QS->QV)
1605 !-------------------------------------------------------------
1606           if(supcol.lt.0.) then
1607             if(qrs(i,k,2).gt.0. .and. rh(i,k,1).lt.1.) then
1608               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1609               psevp(i,k) = (rh(i,k,1)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2)   &
1610                            +precs2*work2(i,k)*coeres)/work1(i,k,1)
1611               psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1612             endif
1613 !-------------------------------------------------------------
1614 ! pgevp: Evaporation of melting graupel [HL A25] [RH84 A19]
1615 !       (T>=T0: QG->QV)
1616 !-------------------------------------------------------------
1617             if(qrs(i,k,3).gt.0. .and. rh(i,k,1).lt.1.) then
1618               coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1619               pgevp(i,k) = (rh(i,k,1)-1.)*(precg1*rslope2(i,k,3)               &
1620                          + precg2*work2(i,k)*coeres)/work1(i,k,1)
1621               pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.)
1622             endif
1623           endif
1624         enddo
1625       enddo
1626 !
1627 !
1628 !----------------------------------------------------------------
1629 !     check mass conservation of generation terms and feedback to the
1630 !     large scale
1631 !
1632       do k = kts, kte
1633         do i = its, ite
1634 !
1635           delta2=0.
1636           delta3=0.
1637           if(qrs(i,k,1).lt.1.e-4 .and. qrs(i,k,2).lt.1.e-4) delta2=1.
1638           if(qrs(i,k,1).lt.1.e-4) delta3=1.
1639           if(t(i,k).le.t0c) then
1640 !
1641 !     cloud water
1642 !
1643             value = max(qmin,qci(i,k,1))
1644             source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))             &
1645                     *dtcld
1646             if (source.gt.value) then
1647               factor = value/source
1648               praut(i,k) = praut(i,k)*factor
1649               pracw(i,k) = pracw(i,k)*factor
1650               paacw(i,k) = paacw(i,k)*factor
1651             endif
1652 !
1653 !     cloud ice
1654 !
1655             value = max(qmin,qci(i,k,2))
1656             source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k)+psaci(i,k)   &
1657                     +pgaci(i,k))*dtcld
1658             if (source.gt.value) then
1659               factor = value/source
1660               psaut(i,k) = psaut(i,k)*factor
1661               pigen(i,k) = pigen(i,k)*factor
1662               pidep(i,k) = pidep(i,k)*factor
1663               praci(i,k) = praci(i,k)*factor
1664               psaci(i,k) = psaci(i,k)*factor
1665               pgaci(i,k) = pgaci(i,k)*factor
1666             endif
1667 !
1668 !     rain
1669 !
1670             value = max(qmin,qrs(i,k,1))
1671             source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k)             &
1672                     +psacr(i,k)+pgacr(i,k))*dtcld
1673             if (source.gt.value) then
1674               factor = value/source
1675               praut(i,k) = praut(i,k)*factor
1676               prevp(i,k) = prevp(i,k)*factor
1677               pracw(i,k) = pracw(i,k)*factor
1678               piacr(i,k) = piacr(i,k)*factor
1679               psacr(i,k) = psacr(i,k)*factor
1680               pgacr(i,k) = pgacr(i,k)*factor
1681             endif
1682 !
1683 !     snow
1684 !
1685             value = max(qmin,qrs(i,k,2))
1686             source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k)             &
1687                      +piacr(i,k)*delta3+praci(i,k)*delta3                      &
1688                      -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2                 &
1689                      +psaci(i,k)-pgacs(i,k) )*dtcld
1690             if (source.gt.value) then
1691               factor = value/source
1692               psdep(i,k) = psdep(i,k)*factor
1693               psaut(i,k) = psaut(i,k)*factor
1694               pgaut(i,k) = pgaut(i,k)*factor
1695               paacw(i,k) = paacw(i,k)*factor
1696               piacr(i,k) = piacr(i,k)*factor
1697               praci(i,k) = praci(i,k)*factor
1698               psaci(i,k) = psaci(i,k)*factor
1699               pracs(i,k) = pracs(i,k)*factor
1700               psacr(i,k) = psacr(i,k)*factor
1701               pgacs(i,k) = pgacs(i,k)*factor
1702             endif
1703 !
1704 !     graupel
1705 !
1706             value = max(qmin,qrs(i,k,3))
1707             source = -(pgdep(i,k)+pgaut(i,k)                                   &
1708                      +piacr(i,k)*(1.-delta3)+praci(i,k)*(1.-delta3)            &
1709                      +psacr(i,k)*(1.-delta2)+pracs(i,k)*(1.-delta2)            &
1710                      +pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
1711             if (source.gt.value) then
1712               factor = value/source
1713               pgdep(i,k) = pgdep(i,k)*factor
1714               pgaut(i,k) = pgaut(i,k)*factor
1715               piacr(i,k) = piacr(i,k)*factor
1716               praci(i,k) = praci(i,k)*factor
1717               psacr(i,k) = psacr(i,k)*factor
1718               pracs(i,k) = pracs(i,k)*factor
1719               paacw(i,k) = paacw(i,k)*factor
1720               pgaci(i,k) = pgaci(i,k)*factor
1721               pgacr(i,k) = pgacr(i,k)*factor
1722               pgacs(i,k) = pgacs(i,k)*factor
1723             endif
1724 !
1725 !     cloud
1726 !
1727             value = max(ncmin,ncr(i,k,2))
1728             source = (nraut(i,k)+nccol(i,k)+nracw(i,k)                         &
1729                     +naacw(i,k)+naacw(i,k))*dtcld
1730             if (source.gt.value) then
1731               factor = value/source
1732               nraut(i,k) = nraut(i,k)*factor
1733               nccol(i,k) = nccol(i,k)*factor
1734               nracw(i,k) = nracw(i,k)*factor
1735               naacw(i,k) = naacw(i,k)*factor
1736             endif
1737 !
1738 !     rain
1739 !
1740             value = max(nrmin,ncr(i,k,3))
1741             source = (-nraut(i,k)+nrcol(i,k)+niacr(i,k)+nsacr(i,k)+ngacr(i,k)  &
1742                      )*dtcld
1743             if (source.gt.value) then
1744               factor = value/source
1745               nraut(i,k) = nraut(i,k)*factor
1746               nrcol(i,k) = nrcol(i,k)*factor
1747               niacr(i,k) = niacr(i,k)*factor
1748               nsacr(i,k) = nsacr(i,k)*factor
1749               ngacr(i,k) = ngacr(i,k)*factor
1750             endif
1751 !
1752             work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k))
1753 !     update
1754             q(i,k) = q(i,k)+work2(i,k)*dtcld
1755             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
1756                            +paacw(i,k)+paacw(i,k))*dtcld,0.)
1757             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
1758                            +prevp(i,k)-piacr(i,k)-pgacr(i,k)                   &
1759                            -psacr(i,k))*dtcld,0.)
1760             qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k)                 &
1761                            +psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k))       &
1762                            *dtcld,0.)
1763             qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k)      &
1764                            -pgaut(i,k)+piacr(i,k)*delta3                       &
1765                            +praci(i,k)*delta3+psaci(i,k)-pgacs(i,k)            &
1766                            -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2)          &
1767                            *dtcld,0.)
1768             qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k)                 &
1769                            +piacr(i,k)*(1.-delta3)                             &
1770                            +praci(i,k)*(1.-delta3)+psacr(i,k)*(1.-delta2)      &
1771                            +pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k)       &
1772                            +pgacr(i,k)+pgacs(i,k))*dtcld,0.)
1773             ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k)     &
1774                            -naacw(i,k)-naacw(i,k))*dtcld,0.)
1775             ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)-niacr(i,k)      &
1776                            -nsacr(i,k)-ngacr(i,k))*dtcld,0.)
1777             xlf = xls-xl(i,k)
1778             xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k))       &
1779                       -xl(i,k)*prevp(i,k)-xlf*(piacr(i,k)+paacw(i,k)           &
1780                       +paacw(i,k)+pgacr(i,k)+psacr(i,k))
1781             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1782           else
1783 !
1784 !     cloud water
1785 !
1786             value = max(qmin,qci(i,k,1))
1787             source= (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))              &
1788                    *dtcld
1789             if (source.gt.value) then
1790               factor = value/source
1791               praut(i,k) = praut(i,k)*factor
1792               pracw(i,k) = pracw(i,k)*factor
1793               paacw(i,k) = paacw(i,k)*factor
1794             endif
1795 !
1796 !     rain
1797 !
1798             value = max(qmin,qrs(i,k,1))
1799             source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k)             &
1800                      -pracw(i,k)-paacw(i,k)-prevp(i,k))*dtcld
1801             if (source.gt.value) then
1802               factor = value/source
1803               praut(i,k) = praut(i,k)*factor
1804               prevp(i,k) = prevp(i,k)*factor
1805               pracw(i,k) = pracw(i,k)*factor
1806               paacw(i,k) = paacw(i,k)*factor
1807               pseml(i,k) = pseml(i,k)*factor
1808               pgeml(i,k) = pgeml(i,k)*factor
1809             endif
1810 !
1811 !     snow
1812 !
1813             value = max(qcrmin,qrs(i,k,2))
1814             source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
1815             if (source.gt.value) then
1816               factor = value/source
1817               pgacs(i,k) = pgacs(i,k)*factor
1818               psevp(i,k) = psevp(i,k)*factor
1819               pseml(i,k) = pseml(i,k)*factor
1820             endif
1821 !
1822 !     graupel
1823 !
1824             value = max(qcrmin,qrs(i,k,3))
1825             source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
1826             if (source.gt.value) then
1827               factor = value/source
1828               pgacs(i,k) = pgacs(i,k)*factor
1829               pgevp(i,k) = pgevp(i,k)*factor
1830               pgeml(i,k) = pgeml(i,k)*factor
1831             endif
1832 !
1833 !     cloud
1834 !
1835             value = max(ncmin,ncr(i,k,2))
1836             source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+naacw(i,k)             &
1837                      +naacw(i,k))*dtcld
1838             if (source.gt.value) then
1839               factor = value/source
1840               nraut(i,k) = nraut(i,k)*factor
1841               nccol(i,k) = nccol(i,k)*factor
1842               nracw(i,k) = nracw(i,k)*factor
1843               naacw(i,k) = naacw(i,k)*factor
1844             endif
1845 !
1846 !     rain
1847 !
1848             value = max(nrmin,ncr(i,k,3))
1849             source = (-nraut(i,k)+nrcol(i,k)-nseml(i,k)-ngeml(i,k)             &
1850                       )*dtcld
1851             if (source.gt.value) then
1852               factor = value/source
1853               nraut(i,k) = nraut(i,k)*factor
1854               nrcol(i,k) = nrcol(i,k)*factor
1855               nseml(i,k) = nseml(i,k)*factor
1856               ngeml(i,k) = ngeml(i,k)*factor
1857             endif
1858 !
1859             work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
1860 !     update
1861             q(i,k) = q(i,k)+work2(i,k)*dtcld
1862             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
1863                     +paacw(i,k)+paacw(i,k))*dtcld,0.)
1864             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
1865                     +prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k)               &
1866                     -pgeml(i,k))*dtcld,0.)
1867             qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k)                 &
1868                     +pseml(i,k))*dtcld,0.)
1869             qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k)                 &
1870                     +pgeml(i,k))*dtcld,0.)
1871             ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k)     &
1872                    -naacw(i,k)-naacw(i,k))*dtcld,0.)
1873             ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)+nseml(i,k)      &
1874                            +ngeml(i,k))*dtcld,0.)
1875             xlf = xls-xl(i,k)
1876             xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k))              &
1877                       -xlf*(pseml(i,k)+pgeml(i,k))
1878             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1879           endif
1880         enddo
1881       enddo
1882 !
1883 ! Inline expansion for fpvs
1884 !         qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1885 !         qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1886       hsub = xls
1887       hvap = xlv0
1888       cvap = cpv
1889       ttp=t0c+0.01
1890       dldt=cvap-cliq
1891       xa=-dldt/rv
1892       xb=xa+hvap/(rv*ttp)
1893       dldti=cvap-cice
1894       xai=-dldti/rv
1895       xbi=xai+hsub/(rv*ttp)
1896       do k = kts, kte
1897         do i = its, ite
1898           tr=ttp/t(i,k)
1899           qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1900           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
1901           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1902           qs(i,k,1) = max(qs(i,k,1),qmin)
1903           tr=ttp/t(i,k)
1904           if(t(i,k).lt.ttp) then
1905             qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
1906           else
1907             qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1908           endif
1909           qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
1910           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
1911           qs(i,k,2) = max(qs(i,k,2),qmin)
1912           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
1913         enddo
1914       enddo
1915 !
1916      do k = kts,kte_in
1917        do i = its,ite
1918          qrs_tmp(i,k,1) = qrs(i,k,1)
1919          qrs_tmp(i,k,2) = qrs(i,k,2)
1920          qrs_tmp(i,k,3) = qrs(i,k,3)
1921          ncr_tmp(i,k)   = ncr(i,k,3)
1922        enddo
1923      enddo
1924 !
1925       call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
1926                      rslope3,work1,workn,its,ite,kts,kte)
1927       do k = kts, kte
1928         do i = its, ite
1929 !-----------------------------------------------------------------
1930 ! re-compute the mean-volume drop diameter       [LH A10]
1931 ! for raindrop distribution
1932 !-----------------------------------------------------------------
1933           avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
1934 !----------------------------------------------------------------
1935 ! Nrevp_s: evaporation/condensation rate of rain   [LH A14]
1936 !        (NR->NC)
1937 !----------------------------------------------------------------
1938           if(avedia(i,k,2).le.di82) then
1939             ncr(i,k,2) = ncr(i,k,2)+ncr(i,k,3)
1940             ncr(i,k,3) = 0.
1941 !----------------------------------------------------------------
1942 ! Prevp_s: evaporation/condensation rate of rain [LH A15] [KK 23]
1943 !        (QR->QC)
1944 !----------------------------------------------------------------
1945             qci(i,k,1) = qci(i,k,1)+qrs(i,k,1)
1946             qrs(i,k,1) = 0.
1947           endif
1948         enddo
1949       enddo
1950 !
1951       do k = kts, kte
1952         do i = its, ite
1953 !---------------------------------------------------------------
1954 ! rate of change of cloud drop concentration due to CCN activation
1955 ! pcact: QV -> QC                                 [LH 8]  [KK 14]
1956 ! ncact: NCCN -> NC                               [LH A2] [KK 12]
1957 !---------------------------------------------------------------
1958           if(rh(i,k,1).gt.1.) then
1959             ncact(i,k) = max(0.,((ncr(i,k,1)+ncr(i,k,2))                       &
1960                        *min(1.,(rh(i,k,1)/satmax)**actk) - ncr(i,k,2)))/dtcld
1961             ncact(i,k) =min(ncact(i,k),max(ncr(i,k,1),0.)/dtcld)
1962             pcact(i,k) = min(4.*pi*denr*(actr*1.E-6)**3*ncact(i,k)/            &
1963                          (3.*den(i,k)),max(q(i,k),0.)/dtcld)
1964             q(i,k) = max(q(i,k)-pcact(i,k)*dtcld,0.)
1965             qci(i,k,1) = max(qci(i,k,1)+pcact(i,k)*dtcld,0.)
1966             ncr(i,k,1) = max(ncr(i,k,1)-ncact(i,k)*dtcld,0.)
1967             ncr(i,k,2) = max(ncr(i,k,2)+ncact(i,k)*dtcld,0.)
1968             t(i,k) = t(i,k)+pcact(i,k)*xl(i,k)/cpm(i,k)*dtcld
1969           endif
1970 !---------------------------------------------------------------
1971 !  pcond:condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1972 !     if there exists additional water vapor condensated/if
1973 !     evaporation of cloud water is not enough to remove subsaturation
1974 !  (QV->QC or QC->QV)
1975 !---------------------------------------------------------------
1976           tr=ttp/t(i,k)
1977           qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1978           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
1979           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1980           qs(i,k,1) = max(qs(i,k,1),qmin)
1981           work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
1982           work2(i,k) = qci(i,k,1)+work1(i,k,1)
1983           pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1984           if(qci(i,k,1).gt.0. .and. work1(i,k,1).lt.0.)                        &
1985             pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1986 !----------------------------------------------------------------
1987 ! ncevp: evpration of Cloud number concentration  [LH A3]
1988 ! (NC->NCCN)
1989 !----------------------------------------------------------------
1990           if(pcond(i,k).eq.-qci(i,k,1)/dtcld) then
1991             ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,2)
1992 ! error correction based on Lei et al. (JGR, 2020)
1993             ncr(i,k,2) = 0.
1994           endif
1995 !
1996           q(i,k) = q(i,k)-pcond(i,k)*dtcld
1997           qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1998           t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1999         enddo
2000       enddo
2001 !
2002 !----------------------------------------------------------------
2003 !     padding for small values
2004 !
2005       do k = kts, kte
2006         do i = its, ite
2007           if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
2008           if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
2009           if(qrs(i,k,1).ge.qcrmin .and. ncr(i,k,3) .ge. nrmin) then
2010             lamdr_tmp(i,k) = exp(log(((pidnr*ncr(i,k,3))                       &
2011                          /(den(i,k)*qrs(i,k,1))))*((.33333333)))
2012             if(lamdr_tmp(i,k) .le. lamdarmin) then
2013               lamdr_tmp(i,k) = lamdarmin
2014               ncr(i,k,3) = den(i,k)*qrs(i,k,1)*lamdr_tmp(i,k)**3/pidnr
2015             elseif(lamdr_tmp(i,k) .ge. lamdarmax) then
2016               lamdr_tmp(i,k) = lamdarmax
2017               ncr(i,k,3) = den(i,k)*qrs(i,k,1)*lamdr_tmp(i,k)**3/pidnr
2018             endif
2019           endif
2020           if(qci(i,k,1).ge.qmin .and. ncr(i,k,2) .ge. ncmin ) then
2021             lamdc_tmp(i,k) = exp(log(((pidnc*ncr(i,k,2))                       &
2022                          /(den(i,k)*qci(i,k,1))))*((.33333333)))
2023             if(lamdc_tmp(i,k) .le. lamdacmin) then
2024               lamdc_tmp(i,k) = lamdacmin
2025               ncr(i,k,2) = den(i,k)*qci(i,k,1)*lamdc_tmp(i,k)**3/pidnc
2026             elseif(lamdc_tmp(i,k) .ge. lamdacmax) then
2027               lamdc_tmp(i,k) = lamdacmax
2028               ncr(i,k,2) = den(i,k)*qci(i,k,1)*lamdc_tmp(i,k)**3/pidnc
2029             endif
2030           endif
2031         enddo
2032       enddo
2033    enddo                  ! big loops
2034 !
2035    end subroutine wdm62d
2036 !------------------------------------------------------------------------------
2037    real function rgmma(x)
2038 !------------------------------------------------------------------------------
2039    implicit none
2040 !-----------------------------------------------------------------------------
2041 !  rgmma function:  use infinite product form
2042 !
2043    real, parameter :: euler = 0.577215664901532
2044    real    :: x, y
2045    integer :: i
2046 !------------------------------------------------------------------------------
2047    if(x.eq.1.) then
2048      rgmma=0.
2049    else
2050      rgmma = x*exp(euler*x)
2051      do i = 1,10000
2052        y = float(i)
2053        rgmma = rgmma*(1.000+x/y)*exp(-x/y)
2054      enddo
2055      rgmma = 1./rgmma
2056    endif
2057 !
2058    end function rgmma
2059 !
2060 !--------------------------------------------------------------------------
2061    real function fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
2062 !--------------------------------------------------------------------------
2063    implicit none
2064 !--------------------------------------------------------------------------
2065       REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti,         &
2066            xai,xbi,ttp,tr
2067       INTEGER ice
2068 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2069       ttp=t0c+0.01
2070       dldt=cvap-cliq
2071       xa=-dldt/rv
2072       xb=xa+hvap/(rv*ttp)
2073       dldti=cvap-cice
2074       xai=-dldti/rv
2075       xbi=xai+hsub/(rv*ttp)
2076       tr=ttp/t
2077       if(t.lt.ttp .and. ice.eq.1) then
2078         fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2079       else
2080         fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2081       endif
2082 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2083       END FUNCTION fpvs
2084 !-------------------------------------------------------------------
2085   SUBROUTINE wdm6init(den0,denr,dens,cl,cpv, ccn0, hail_opt, allowed_to_read) ! RAS
2086 !-------------------------------------------------------------------
2087   IMPLICIT NONE
2088 !-------------------------------------------------------------------
2089 !.... constants which may not be tunable
2090    REAL, INTENT(IN) :: den0,denr,dens,cl,cpv,ccn0
2091    INTEGER, INTENT(IN) :: hail_opt  ! RAS
2092    LOGICAL, INTENT(IN) :: allowed_to_read
2093
2094 ! RAS13.1 define graupel parameters as graupel-like or hail-like,
2095 !         depending on namelist option
2096       IF (hail_opt .eq. 1) THEN !Hail!
2097          n0g       = 4.e4
2098          deng      = 700.
2099          avtg      = 285.0
2100          bvtg      = 0.8
2101          lamdagmax = 2.e4
2102       ELSE !Graupel!
2103          n0g       = 4.e6
2104          deng      = 500
2105          avtg      = 330.0
2106          bvtg      = 0.8
2107          lamdagmax = 6.e4
2108       ENDIF
2109 !
2110    pi = 4.*atan(1.)
2111    xlv1 = cl-cpv
2112 !
2113    qc0  = 4./3.*pi*denr*r0**3.*xncr0/den0
2114    qc1  = 4./3.*pi*denr*r0**3.*xncr1/den0
2115    qck1 = .104*9.8*peaut/(denr)**(1./3.)/xmyu*den0**(4./3.) ! 4706.08203
2116    pidnc = pi*denr/6.
2117 !
2118    bvtr1 = 1.+bvtr
2119    bvtr2 = 2.+bvtr
2120    bvtr3 = 3.+bvtr
2121    bvtr4 = 4.+bvtr
2122    bvtr5 = 5.+bvtr
2123    bvtr6 = 6.+bvtr
2124    bvtr7 = 7.+bvtr
2125    bvtr2o5 = 2.5+.5*bvtr
2126    bvtr3o5 = 3.5+.5*bvtr
2127    g1pbr = rgmma(bvtr1)
2128    g2pbr = rgmma(bvtr2)
2129    g3pbr = rgmma(bvtr3)
2130    g4pbr = rgmma(bvtr4)            ! 17.837825
2131    g5pbr = rgmma(bvtr5)
2132    g6pbr = rgmma(bvtr6)
2133    g7pbr = rgmma(bvtr7)
2134    g5pbro2 = rgmma(bvtr2o5)
2135    g7pbro2 = rgmma(bvtr3o5)
2136    pvtr = avtr*g5pbr/24.
2137    pvtrn = avtr*g2pbr
2138    eacrr = 1.0
2139    pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
2140    precr1 = 2.*pi*1.56
2141    precr2 = 2.*pi*.31*avtr**.5*g7pbro2
2142    pidn0r =  pi*denr*n0r
2143    pidnr = 4.*pi*denr
2144 !
2145    xmmax = (dimax/dicon)**2
2146    roqimax = 2.08e22*dimax**8
2147 !
2148    bvts1 = 1.+bvts
2149    bvts2 = 2.5+.5*bvts
2150    bvts3 = 3.+bvts
2151    bvts4 = 4.+bvts
2152    g1pbs = rgmma(bvts1)    !.8875
2153    g3pbs = rgmma(bvts3)
2154    g4pbs = rgmma(bvts4)    ! 12.0786
2155    g5pbso2 = rgmma(bvts2)
2156    pvts = avts*g4pbs/6.
2157    pacrs = pi*n0s*avts*g3pbs*.25
2158    precs1 = 4.*n0s*.65
2159    precs2 = 4.*n0s*.44*avts**.5*g5pbso2
2160    pidn0s =  pi*dens*n0s
2161 !
2162    pacrc = pi*n0s*avts*g3pbs*.25*eacrc
2163 !
2164    bvtg1 = 1.+bvtg
2165    bvtg2 = 2.5+.5*bvtg
2166    bvtg3 = 3.+bvtg
2167    bvtg4 = 4.+bvtg
2168    g1pbg = rgmma(bvtg1)
2169    g3pbg = rgmma(bvtg3)
2170    g4pbg = rgmma(bvtg4)
2171    g5pbgo2 = rgmma(bvtg2)
2172    pacrg = pi*n0g*avtg*g3pbg*.25
2173    pvtg = avtg*g4pbg/6.
2174    precg1 = 2.*pi*n0g*.78
2175    precg2 = 2.*pi*n0g*.31*avtg**.5*g5pbgo2
2176    pidn0g =  pi*deng*n0g
2177 !
2178    rslopecmax = 1./lamdacmax
2179    rslopermax = 1./lamdarmax
2180    rslopesmax = 1./lamdasmax
2181    rslopegmax = 1./lamdagmax
2182    rsloperbmax = rslopermax ** bvtr
2183    rslopesbmax = rslopesmax ** bvts
2184    rslopegbmax = rslopegmax ** bvtg
2185    rslopec2max = rslopecmax * rslopecmax
2186    rsloper2max = rslopermax * rslopermax
2187    rslopes2max = rslopesmax * rslopesmax
2188    rslopeg2max = rslopegmax * rslopegmax
2189    rslopec3max = rslopec2max * rslopecmax
2190    rsloper3max = rsloper2max * rslopermax
2191    rslopes3max = rslopes2max * rslopesmax
2192    rslopeg3max = rslopeg2max * rslopegmax
2193
2194 !+---+-----------------------------------------------------------------+
2195 !..Set these variables needed for computing radar reflectivity.  These
2196 !.. get used within radar_init to create other variables used in the
2197 !.. radar module.
2198    xam_r = PI*denr/6.
2199    xbm_r = 3.
2200    xmu_r = 1.
2201    xam_s = PI*dens/6.
2202    xbm_s = 3.
2203    xmu_s = 0.
2204    xam_g = PI*deng/6.
2205    xbm_g = 3.
2206    xmu_g = 0.
2207
2208    call radar_init
2209 !+---+-----------------------------------------------------------------+
2210 !
2211   END SUBROUTINE wdm6init
2212 !------------------------------------------------------------------------------
2213       subroutine slope_wdm6(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
2214                             vt,vtn,its,ite,kts,kte)
2215   IMPLICIT NONE
2216   INTEGER       ::               its,ite, jts,jte, kts,kte
2217   REAL, DIMENSION( its:ite , kts:kte,3) ::                                     &
2218                                                                           qrs, &
2219                                                                        rslope, &
2220                                                                       rslopeb, &
2221                                                                       rslope2, &
2222                                                                       rslope3, &
2223                                                                            vt
2224   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
2225                                                                           ncr, &
2226                                                                           vtn, &
2227                                                                           den, &
2228                                                                        denfac, &
2229                                                                             t
2230   REAL, PARAMETER  :: t0c = 273.15
2231   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
2232                                                                        n0sfac
2233   REAL       ::  lamdar, lamdas, lamdag, x, y, z, supcol
2234   integer :: i, j, k
2235 !----------------------------------------------------------------
2236 !     size distributions: (x=mixing ratio, y=air density):
2237 !     valid for mixing ratio > 1.e-9 kg/kg.
2238 !
2239 !  Optimizatin : A**B => exp(log(A)*(B))
2240       lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
2241       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
2242       lamdag(x,y)=   sqrt(sqrt(pidn0g/(x*y)))      ! (pidn0g/(x*y))**.25
2243 !
2244       do k = kts, kte
2245         do i = its, ite
2246           supcol = t0c-t(i,k)
2247 !---------------------------------------------------------------
2248 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
2249 !---------------------------------------------------------------
2250           n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
2251           if(qrs(i,k,1).le.qcrmin .or. ncr(i,k).le.nrmin ) then
2252             rslope(i,k,1) = rslopermax
2253             rslopeb(i,k,1) = rsloperbmax
2254             rslope2(i,k,1) = rsloper2max
2255             rslope3(i,k,1) = rsloper3max
2256           else
2257             rslope(i,k,1) = min(1./lamdar(qrs(i,k,1),den(i,k),ncr(i,k)),1.e-3)
2258             rslopeb(i,k,1) = rslope(i,k,1)**bvtr
2259             rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
2260             rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
2261           endif
2262           if(qrs(i,k,2).le.qcrmin) then
2263             rslope(i,k,2) = rslopesmax
2264             rslopeb(i,k,2) = rslopesbmax
2265             rslope2(i,k,2) = rslopes2max
2266             rslope3(i,k,2) = rslopes3max
2267           else
2268             rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
2269             rslopeb(i,k,2) = rslope(i,k,2)**bvts
2270             rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
2271             rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
2272           endif
2273           if(qrs(i,k,3).le.qcrmin) then
2274             rslope(i,k,3) = rslopegmax
2275             rslopeb(i,k,3) = rslopegbmax
2276             rslope2(i,k,3) = rslopeg2max
2277             rslope3(i,k,3) = rslopeg3max
2278           else
2279             rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
2280             rslopeb(i,k,3) = rslope(i,k,3)**bvtg
2281             rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
2282             rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
2283           endif
2284           vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
2285           vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
2286           vt(i,k,3) = pvtg*rslopeb(i,k,3)*denfac(i,k)
2287           vtn(i,k) = pvtrn*rslopeb(i,k,1)*denfac(i,k)
2288           if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
2289           if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
2290           if(qrs(i,k,3).le.0.0) vt(i,k,3) = 0.0
2291           if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
2292         enddo
2293       enddo
2294   END subroutine slope_wdm6
2295 !-----------------------------------------------------------------------------
2296       subroutine slope_rain(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
2297                             vt,vtn,its,ite,kts,kte)
2298   IMPLICIT NONE
2299   INTEGER       ::               its,ite, jts,jte, kts,kte
2300   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
2301                                                                           qrs, &
2302                                                                           ncr, &
2303                                                                        rslope, &
2304                                                                       rslopeb, &
2305                                                                       rslope2, &
2306                                                                       rslope3, &
2307                                                                            vt, &
2308                                                                           vtn, &
2309                                                                           den, &
2310                                                                        denfac, &
2311                                                                             t
2312   REAL, PARAMETER  :: t0c = 273.15
2313   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
2314                                                                        n0sfac
2315   REAL       ::  lamdar, x, y, z, supcol
2316   integer :: i, j, k
2317 !----------------------------------------------------------------
2318 !     size distributions: (x=mixing ratio, y=air density):
2319 !     valid for mixing ratio > 1.e-9 kg/kg.
2320       lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
2321 !
2322       do k = kts, kte
2323         do i = its, ite
2324           if(qrs(i,k).le.qcrmin .or. ncr(i,k).le.nrmin) then
2325             rslope(i,k) = rslopermax
2326             rslopeb(i,k) = rsloperbmax
2327             rslope2(i,k) = rsloper2max
2328             rslope3(i,k) = rsloper3max
2329           else
2330             rslope(i,k) = min(1./lamdar(qrs(i,k),den(i,k),ncr(i,k)),1.e-3)
2331             rslopeb(i,k) = rslope(i,k)**bvtr
2332             rslope2(i,k) = rslope(i,k)*rslope(i,k)
2333             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
2334           endif
2335           vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
2336           vtn(i,k) = pvtrn*rslopeb(i,k)*denfac(i,k)
2337           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
2338           if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
2339         enddo
2340       enddo
2341   END subroutine slope_rain
2342 !------------------------------------------------------------------------------
2343       subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
2344                             vt,its,ite,kts,kte)
2345   IMPLICIT NONE
2346   INTEGER       ::               its,ite, jts,jte, kts,kte
2347   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
2348                                                                           qrs, &
2349                                                                        rslope, &
2350                                                                       rslopeb, &
2351                                                                       rslope2, &
2352                                                                       rslope3, &
2353                                                                            vt, &
2354                                                                           den, &
2355                                                                        denfac, &
2356                                                                             t
2357   REAL, PARAMETER  :: t0c = 273.15
2358   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
2359                                                                        n0sfac
2360   REAL       ::  lamdas, x, y, z, supcol
2361   integer :: i, j, k
2362 !----------------------------------------------------------------
2363 !     size distributions: (x=mixing ratio, y=air density):
2364 !     valid for mixing ratio > 1.e-9 kg/kg.
2365       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
2366 !
2367       do k = kts, kte
2368         do i = its, ite
2369           supcol = t0c-t(i,k)
2370 !---------------------------------------------------------------
2371 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
2372 !---------------------------------------------------------------
2373           n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
2374           if(qrs(i,k).le.qcrmin)then
2375             rslope(i,k) = rslopesmax
2376             rslopeb(i,k) = rslopesbmax
2377             rslope2(i,k) = rslopes2max
2378             rslope3(i,k) = rslopes3max
2379           else
2380             rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
2381             rslopeb(i,k) = rslope(i,k)**bvts
2382             rslope2(i,k) = rslope(i,k)*rslope(i,k)
2383             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
2384           endif
2385           vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
2386           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
2387         enddo
2388       enddo
2389   END subroutine slope_snow
2390 !----------------------------------------------------------------------------------
2391       subroutine slope_graup(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
2392                             vt,its,ite,kts,kte)
2393   IMPLICIT NONE
2394   INTEGER       ::               its,ite, jts,jte, kts,kte
2395   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
2396                                                                           qrs, &
2397                                                                        rslope, &
2398                                                                       rslopeb, &
2399                                                                       rslope2, &
2400                                                                       rslope3, &
2401                                                                            vt, &
2402                                                                           den, &
2403                                                                        denfac, &
2404                                                                             t
2405   REAL, PARAMETER  :: t0c = 273.15
2406   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
2407                                                                        n0sfac
2408   REAL       ::  lamdag, x, y, z, supcol
2409   integer :: i, j, k
2410 !----------------------------------------------------------------
2411 !     size distributions: (x=mixing ratio, y=air density):
2412 !     valid for mixing ratio > 1.e-9 kg/kg.
2413       lamdag(x,y)=   sqrt(sqrt(pidn0g/(x*y)))      ! (pidn0g/(x*y))**.25
2414 !
2415       do k = kts, kte
2416         do i = its, ite
2417 !---------------------------------------------------------------
2418 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
2419 !---------------------------------------------------------------
2420           if(qrs(i,k).le.qcrmin)then
2421             rslope(i,k) = rslopegmax
2422             rslopeb(i,k) = rslopegbmax
2423             rslope2(i,k) = rslopeg2max
2424             rslope3(i,k) = rslopeg3max
2425           else
2426             rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k))
2427             rslopeb(i,k) = rslope(i,k)**bvtg
2428             rslope2(i,k) = rslope(i,k)*rslope(i,k)
2429             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
2430           endif
2431           vt(i,k) = pvtg*rslopeb(i,k)*denfac(i,k)
2432           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
2433         enddo
2434       enddo
2435   END subroutine slope_graup
2436 !---------------------------------------------------------------------------------
2437 !-------------------------------------------------------------------
2438       SUBROUTINE nislfv_rain_plmr(im,km,denl,denfacl,tkl,dzl,wwl,rql,rnl,precip,dt,id,iter,rid)
2439 !-------------------------------------------------------------------
2440 !
2441 ! for non-iteration semi-Lagrangain forward advection for cloud
2442 ! with mass conservation and positive definite advection
2443 ! 2nd order interpolation with monotonic piecewise linear method
2444 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
2445 !
2446 ! dzl    depth of model layer in meter
2447 ! wwl    terminal velocity at model layer m/s
2448 ! rql    cloud density*mixing ration
2449 ! precip precipitation
2450 ! dt     time step
2451 ! id     kind of precip: 0 test case; 1 raindrop
2452 ! iter   how many time to guess mean terminal velocity: 0 pure forward.
2453 !        0 : use departure wind for advection
2454 !        1 : use mean wind for advection
2455 !        > 1 : use mean wind after iter-1 iterations
2456 !        rid : 1 for number 0 for mixing ratio
2457 !
2458 ! author: hann-ming henry juang <henry.juang@noaa.gov>
2459 !         implemented by song-you hong
2460 !
2461       implicit none
2462       integer  im,km,id
2463       real  dt
2464       real  dzl(im,km),wwl(im,km),rql(im,km),rnl(im,km),precip(im)
2465       real  denl(im,km),denfacl(im,km),tkl(im,km)
2466 !
2467       integer  i,k,n,m,kk,kb,kt,iter,rid
2468       real  tl,tl2,qql,dql,qqd
2469       real  th,th2,qqh,dqh
2470       real  zsum,qsum,dim,dip,c1,con1,fa1,fa2
2471       real  allold, allnew, zz, dzamin, cflmax, decfl
2472       real  dz(km), ww(km), qq(km), nr(km), wd(km), wa(km), wa2(km), was(km)
2473       real  den(km), denfac(km), tk(km)
2474       real  wi(km+1), zi(km+1), za(km+1)
2475       real  qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
2476       real  dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
2477 !
2478       precip(:) = 0.0
2479 !
2480       i_loop : do i=1,im
2481 ! -----------------------------------
2482       dz(:) = dzl(i,:)
2483       qq(:) = rql(i,:)
2484       nr(:) = rnl(i,:)
2485       if(rid .eq. 1) nr(:) = rnl(i,:)/denl(i,:)
2486       ww(:) = wwl(i,:)
2487       den(:) = denl(i,:)
2488       denfac(:) = denfacl(i,:)
2489       tk(:) = tkl(i,:)
2490 ! skip for no precipitation for all layers
2491       allold = 0.0
2492       do k=1,km
2493         allold = allold + qq(k)
2494       enddo
2495       if(allold.le.0.0) then
2496         cycle i_loop
2497       endif
2498 !
2499 ! compute interface values
2500       zi(1)=0.0
2501       do k=1,km
2502         zi(k+1) = zi(k)+dz(k)
2503       enddo
2504 !
2505 ! save departure wind
2506       wd(:) = ww(:)
2507       n=1
2508  100  continue
2509 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
2510 ! 2nd order interpolation to get wi
2511       wi(1) = ww(1)
2512       wi(km+1) = ww(km)
2513       do k=2,km
2514         wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
2515       enddo
2516 ! 3rd order interpolation to get wi
2517       fa1 = 9./16.
2518       fa2 = 1./16.
2519       wi(1) = ww(1)
2520       wi(2) = 0.5*(ww(2)+ww(1))
2521       do k=3,km-1
2522         wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
2523       enddo
2524       wi(km) = 0.5*(ww(km)+ww(km-1))
2525       wi(km+1) = ww(km)
2526 !
2527 ! terminate of top of raingroup
2528       do k=2,km
2529         if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
2530       enddo
2531 !
2532 ! diffusivity of wi
2533       con1 = 0.05
2534       do k=km,1,-1
2535         decfl = (wi(k+1)-wi(k))*dt/dz(k)
2536         if( decfl .gt. con1 ) then
2537           wi(k) = wi(k+1) - con1*dz(k)/dt
2538         endif
2539       enddo
2540 ! compute arrival point
2541       do k=1,km+1
2542         za(k) = zi(k) - wi(k)*dt
2543       enddo
2544 !
2545       do k=1,km
2546         dza(k) = za(k+1)-za(k)
2547       enddo
2548       dza(km+1) = zi(km+1) - za(km+1)
2549 !
2550 ! computer deformation at arrival point
2551       do k=1,km
2552         qa(k) = qq(k)*dz(k)/dza(k)
2553         qr(k) = qa(k)/den(k)
2554         if(rid .eq. 1) qr(k) = qa(K)
2555       enddo
2556       qa(km+1) = 0.0
2557 !     call maxmin(km,1,qa,' arrival points ')
2558 !
2559 ! compute arrival terminal velocity, and estimate mean terminal velocity
2560 ! then back to use mean terminal velocity
2561       if( n.le.iter ) then
2562         if(rid.eq.1) then
2563         call slope_rain(nr,qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,wa2,1,1,1,km)
2564         else
2565         call slope_rain(qr,nr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,wa2,1,1,1,km)
2566         endif
2567         if(rid.eq.1) wa(:) = wa2(:)
2568         if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
2569         do k=1,km
2570 !#ifdef DEBUG
2571 !        print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
2572 !#endif
2573 ! mean wind is average of departure and new arrival winds
2574           ww(k) = 0.5* ( wd(k)+wa(k) )
2575         enddo
2576         was(:) = wa(:)
2577         n=n+1
2578         go to 100
2579       endif
2580 !
2581 ! estimate values at arrival cell interface with monotone
2582       do k=2,km
2583         dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
2584         dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
2585         if( dip*dim.le.0.0 ) then
2586           qmi(k)=qa(k)
2587           qpi(k)=qa(k)
2588         else
2589           qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
2590           qmi(k)=2.0*qa(k)-qpi(k)
2591           if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
2592             qpi(k) = qa(k)
2593             qmi(k) = qa(k)
2594           endif
2595         endif
2596       enddo
2597       qpi(1)=qa(1)
2598       qmi(1)=qa(1)
2599       qmi(km+1)=qa(km+1)
2600       qpi(km+1)=qa(km+1)
2601 !
2602 ! interpolation to regular point
2603       qn = 0.0
2604       kb=1
2605       kt=1
2606       intp : do k=1,km
2607              kb=max(kb-1,1)
2608              kt=max(kt-1,1)
2609 ! find kb and kt
2610              if( zi(k).ge.za(km+1) ) then
2611                exit intp
2612              else
2613                find_kb : do kk=kb,km
2614                          if( zi(k).le.za(kk+1) ) then
2615                            kb = kk
2616                            exit find_kb
2617                          else
2618                            cycle find_kb
2619                          endif
2620                enddo find_kb
2621                find_kt : do kk=kt,km
2622                          if( zi(k+1).le.za(kk) ) then
2623                            kt = kk
2624                            exit find_kt
2625                          else
2626                            cycle find_kt
2627                          endif
2628                enddo find_kt
2629                kt = kt - 1
2630 ! compute q with piecewise constant method
2631                if( kt.eq.kb ) then
2632                  tl=(zi(k)-za(kb))/dza(kb)
2633                  th=(zi(k+1)-za(kb))/dza(kb)
2634                  tl2=tl*tl
2635                  th2=th*th
2636                  qqd=0.5*(qpi(kb)-qmi(kb))
2637                  qqh=qqd*th2+qmi(kb)*th
2638                  qql=qqd*tl2+qmi(kb)*tl
2639                  qn(k) = (qqh-qql)/(th-tl)
2640                else if( kt.gt.kb ) then
2641                  tl=(zi(k)-za(kb))/dza(kb)
2642                  tl2=tl*tl
2643                  qqd=0.5*(qpi(kb)-qmi(kb))
2644                  qql=qqd*tl2+qmi(kb)*tl
2645                  dql = qa(kb)-qql
2646                  zsum  = (1.-tl)*dza(kb)
2647                  qsum  = dql*dza(kb)
2648                  if( kt-kb.gt.1 ) then
2649                  do m=kb+1,kt-1
2650                    zsum = zsum + dza(m)
2651                    qsum = qsum + qa(m) * dza(m)
2652                  enddo
2653                  endif
2654                  th=(zi(k+1)-za(kt))/dza(kt)
2655                  th2=th*th
2656                  qqd=0.5*(qpi(kt)-qmi(kt))
2657                  dqh=qqd*th2+qmi(kt)*th
2658                  zsum  = zsum + th*dza(kt)
2659                  qsum  = qsum + dqh*dza(kt)
2660                  qn(k) = qsum/zsum
2661                endif
2662                cycle intp
2663              endif
2664 !
2665        enddo intp
2666 !
2667 ! rain out
2668       sum_precip: do k=1,km
2669                     if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
2670                       precip(i) = precip(i) + qa(k)*dza(k)
2671                       cycle sum_precip
2672                     else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
2673                       precip(i) = precip(i) + qa(k)*(0.0-za(k))
2674                       exit sum_precip
2675                     endif
2676                     exit sum_precip
2677       enddo sum_precip
2678 !
2679 ! replace the new values
2680       rql(i,:) = qn(:)
2681 !
2682 ! ----------------------------------
2683       enddo i_loop
2684 !
2685   END SUBROUTINE nislfv_rain_plmr
2686 !-------------------------------------------------------------------
2687       SUBROUTINE nislfv_rain_plm6(im,km,denl,denfacl,tkl,dzl,wwl,rql,rql2, precip1, precip2,dt,id,iter)
2688 !-------------------------------------------------------------------
2689 !
2690 ! for non-iteration semi-Lagrangain forward advection for cloud
2691 ! with mass conservation and positive definite advection
2692 ! 2nd order interpolation with monotonic piecewise linear method
2693 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
2694 !
2695 ! dzl    depth of model layer in meter
2696 ! wwl    terminal velocity at model layer m/s
2697 ! rql    cloud density*mixing ration
2698 ! precip precipitation
2699 ! dt     time step
2700 ! id     kind of precip: 0 test case; 1 raindrop
2701 ! iter   how many time to guess mean terminal velocity: 0 pure forward.
2702 !        0 : use departure wind for advection
2703 !        1 : use mean wind for advection
2704 !        > 1 : use mean wind after iter-1 iterations
2705 !
2706 ! author: hann-ming henry juang <henry.juang@noaa.gov>
2707 !         implemented by song-you hong
2708 !
2709       implicit none
2710       integer  im,km,id
2711       real  dt
2712       real  dzl(im,km),wwl(im,km),rql(im,km),rql2(im,km),precip(im),precip1(im),precip2(im)
2713       real  denl(im,km),denfacl(im,km),tkl(im,km)
2714 !
2715       integer  i,k,n,m,kk,kb,kt,iter,ist
2716       real  tl,tl2,qql,dql,qqd
2717       real  th,th2,qqh,dqh
2718       real  zsum,qsum,dim,dip,c1,con1,fa1,fa2
2719       real  allold, allnew, zz, dzamin, cflmax, decfl
2720       real  dz(km), ww(km), qq(km), qq2(km), wd(km), wa(km), wa2(km), was(km)
2721       real  den(km), denfac(km), tk(km)
2722       real  wi(km+1), zi(km+1), za(km+1)
2723       real  qn(km), qr(km),qr2(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
2724       real  dza(km+1), qa(km+1), qa2(km+1),qmi(km+1), qpi(km+1)
2725 !
2726       precip(:) = 0.0
2727       precip1(:) = 0.0
2728       precip2(:) = 0.0
2729 !
2730       i_loop : do i=1,im
2731 ! -----------------------------------
2732       dz(:) = dzl(i,:)
2733       qq(:) = rql(i,:)
2734       qq2(:) = rql2(i,:)
2735       ww(:) = wwl(i,:)
2736       den(:) = denl(i,:)
2737       denfac(:) = denfacl(i,:)
2738       tk(:) = tkl(i,:)
2739 ! skip for no precipitation for all layers
2740       allold = 0.0
2741       do k=1,km
2742         allold = allold + qq(k) + qq2(k)
2743       enddo
2744       if(allold.le.0.0) then
2745         cycle i_loop
2746       endif
2747 !
2748 ! compute interface values
2749       zi(1)=0.0
2750       do k=1,km
2751         zi(k+1) = zi(k)+dz(k)
2752       enddo
2753 !
2754 ! save departure wind
2755       wd(:) = ww(:)
2756       n=1
2757  100  continue
2758 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
2759 ! 2nd order interpolation to get wi
2760       wi(1) = ww(1)
2761       wi(km+1) = ww(km)
2762       do k=2,km
2763         wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
2764       enddo
2765 ! 3rd order interpolation to get wi
2766       fa1 = 9./16.
2767       fa2 = 1./16.
2768       wi(1) = ww(1)
2769       wi(2) = 0.5*(ww(2)+ww(1))
2770       do k=3,km-1
2771         wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
2772       enddo
2773       wi(km) = 0.5*(ww(km)+ww(km-1))
2774       wi(km+1) = ww(km)
2775 !
2776 ! terminate of top of raingroup
2777       do k=2,km
2778         if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
2779       enddo
2780 !
2781 ! diffusivity of wi
2782       con1 = 0.05
2783       do k=km,1,-1
2784         decfl = (wi(k+1)-wi(k))*dt/dz(k)
2785         if( decfl .gt. con1 ) then
2786           wi(k) = wi(k+1) - con1*dz(k)/dt
2787         endif
2788       enddo
2789 ! compute arrival point
2790       do k=1,km+1
2791         za(k) = zi(k) - wi(k)*dt
2792       enddo
2793 !
2794       do k=1,km
2795         dza(k) = za(k+1)-za(k)
2796       enddo
2797       dza(km+1) = zi(km+1) - za(km+1)
2798 !
2799 ! computer deformation at arrival point
2800       do k=1,km
2801         qa(k) = qq(k)*dz(k)/dza(k)
2802         qa2(k) = qq2(k)*dz(k)/dza(k)
2803         qr(k) = qa(k)/den(k)
2804         qr2(k) = qa2(k)/den(k)
2805       enddo
2806       qa(km+1) = 0.0
2807       qa2(km+1) = 0.0
2808 !     call maxmin(km,1,qa,' arrival points ')
2809 !
2810 ! compute arrival terminal velocity, and estimate mean terminal velocity
2811 ! then back to use mean terminal velocity
2812       if( n.le.iter ) then
2813         call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
2814         call slope_graup(qr2,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa2,1,1,1,km)
2815         do k = 1, km
2816           tmp(k) = max((qr(k)+qr2(k)), 1.E-15)
2817           IF ( tmp(k) .gt. 1.e-15 ) THEN
2818             wa(k) = (wa(k)*qr(k) + wa2(k)*qr2(k))/tmp(k)
2819           ELSE
2820             wa(k) = 0.
2821           ENDIF
2822         enddo
2823         if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
2824         do k=1,km
2825 !#ifdef DEBUG
2826 !        print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k), &
2827 !           ww(k),wa(k)
2828 !#endif
2829 ! mean wind is average of departure and new arrival winds
2830           ww(k) = 0.5* ( wd(k)+wa(k) )
2831         enddo
2832         was(:) = wa(:)
2833         n=n+1
2834         go to 100
2835       endif
2836       ist_loop : do ist = 1, 2
2837       if (ist.eq.2) then
2838        qa(:) = qa2(:)
2839       endif
2840 !
2841       precip(i) = 0.
2842 !
2843 ! estimate values at arrival cell interface with monotone
2844       do k=2,km
2845         dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
2846         dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
2847         if( dip*dim.le.0.0 ) then
2848           qmi(k)=qa(k)
2849           qpi(k)=qa(k)
2850         else
2851           qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
2852           qmi(k)=2.0*qa(k)-qpi(k)
2853           if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
2854             qpi(k) = qa(k)
2855             qmi(k) = qa(k)
2856           endif
2857         endif
2858       enddo
2859       qpi(1)=qa(1)
2860       qmi(1)=qa(1)
2861       qmi(km+1)=qa(km+1)
2862       qpi(km+1)=qa(km+1)
2863 !
2864 ! interpolation to regular point
2865       qn = 0.0
2866       kb=1
2867       kt=1
2868       intp : do k=1,km
2869              kb=max(kb-1,1)
2870              kt=max(kt-1,1)
2871 ! find kb and kt
2872              if( zi(k).ge.za(km+1) ) then
2873                exit intp
2874              else
2875                find_kb : do kk=kb,km
2876                          if( zi(k).le.za(kk+1) ) then
2877                            kb = kk
2878                            exit find_kb
2879                          else
2880                            cycle find_kb
2881                          endif
2882                enddo find_kb
2883                find_kt : do kk=kt,km
2884                          if( zi(k+1).le.za(kk) ) then
2885                            kt = kk
2886                            exit find_kt
2887                          else
2888                            cycle find_kt
2889                          endif
2890                enddo find_kt
2891                kt = kt - 1
2892 ! compute q with piecewise constant method
2893                if( kt.eq.kb ) then
2894                  tl=(zi(k)-za(kb))/dza(kb)
2895                  th=(zi(k+1)-za(kb))/dza(kb)
2896                  tl2=tl*tl
2897                  th2=th*th
2898                  qqd=0.5*(qpi(kb)-qmi(kb))
2899                  qqh=qqd*th2+qmi(kb)*th
2900                  qql=qqd*tl2+qmi(kb)*tl
2901                  qn(k) = (qqh-qql)/(th-tl)
2902                else if( kt.gt.kb ) then
2903                  tl=(zi(k)-za(kb))/dza(kb)
2904                  tl2=tl*tl
2905                  qqd=0.5*(qpi(kb)-qmi(kb))
2906                  qql=qqd*tl2+qmi(kb)*tl
2907                  dql = qa(kb)-qql
2908                  zsum  = (1.-tl)*dza(kb)
2909                  qsum  = dql*dza(kb)
2910                  if( kt-kb.gt.1 ) then
2911                  do m=kb+1,kt-1
2912                    zsum = zsum + dza(m)
2913                    qsum = qsum + qa(m) * dza(m)
2914                  enddo
2915                  endif
2916                  th=(zi(k+1)-za(kt))/dza(kt)
2917                  th2=th*th
2918                  qqd=0.5*(qpi(kt)-qmi(kt))
2919                  dqh=qqd*th2+qmi(kt)*th
2920                  zsum  = zsum + th*dza(kt)
2921                  qsum  = qsum + dqh*dza(kt)
2922                  qn(k) = qsum/zsum
2923                endif
2924                cycle intp
2925              endif
2926 !
2927        enddo intp
2928 !
2929 ! rain out
2930       sum_precip: do k=1,km
2931                     if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
2932                       precip(i) = precip(i) + qa(k)*dza(k)
2933                       cycle sum_precip
2934                     else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
2935                       precip(i) = precip(i) + qa(k)*(0.0-za(k))
2936                       exit sum_precip
2937                     endif
2938                     exit sum_precip
2939       enddo sum_precip
2940 !
2941 ! replace the new values
2942       if(ist.eq.1) then
2943         rql(i,:) = qn(:)
2944         precip1(i) = precip(i)
2945       else
2946         rql2(i,:) = qn(:)
2947         precip2(i) = precip(i)
2948       endif
2949       enddo ist_loop
2950 !
2951 ! ----------------------------------
2952       enddo i_loop
2953 !
2954   END SUBROUTINE nislfv_rain_plm6
2955
2956 !+---+-----------------------------------------------------------------+
2957       subroutine refl10cm_wdm6 (qv1d, qr1d, nr1d, qs1d, qg1d,           &
2958                        t1d, p1d, dBZ, kts, kte, ii, jj)
2959
2960       IMPLICIT NONE
2961
2962 !..Sub arguments
2963       INTEGER, INTENT(IN):: kts, kte, ii, jj
2964       REAL, DIMENSION(kts:kte), INTENT(IN)::                            &
2965                       qv1d, qr1d, nr1d, qs1d, qg1d, t1d, p1d
2966       REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
2967
2968 !..Local variables
2969       REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
2970       REAL, DIMENSION(kts:kte):: rr, nr, rs, rg
2971       REAL:: temp_C
2972
2973       DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg
2974       DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g
2975       DOUBLE PRECISION:: lamr, lams, lamg
2976       LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
2977
2978       REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
2979       DOUBLE PRECISION:: fmelt_s, fmelt_g
2980
2981       INTEGER:: i, k, k_0, kbot, n
2982       LOGICAL:: melti
2983
2984       DOUBLE PRECISION:: cback, x, eta, f_d
2985       REAL, PARAMETER:: R=287.
2986
2987 !+---+
2988
2989       do k = kts, kte
2990          dBZ(k) = -35.0
2991       enddo
2992
2993 !+---+-----------------------------------------------------------------+
2994 !..Put column of data into local arrays.
2995 !+---+-----------------------------------------------------------------+
2996       do k = kts, kte
2997          temp(k) = t1d(k)
2998          temp_C = min(-0.001, temp(K)-273.15)
2999          qv(k) = MAX(1.E-10, qv1d(k))
3000          pres(k) = p1d(k)
3001          rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
3002
3003          if (qr1d(k) .gt. 1.E-9) then
3004             rr(k) = qr1d(k)*rho(k)
3005             nr(k) = nr1d(k)*rho(k)
3006             lamr = (xam_r*xcrg(3)*xorg2*nr(k)/rr(k))**xobmr
3007             ilamr(k) = 1./lamr
3008             N0_r(k) = nr(k)*xorg2*lamr**xcre(2)
3009             L_qr(k) = .true.
3010          else
3011             rr(k) = 1.E-12
3012             nr(k) = 1.E-12
3013             L_qr(k) = .false.
3014          endif
3015
3016          if (qs1d(k) .gt. 1.E-9) then
3017             rs(k) = qs1d(k)*rho(k)
3018             N0_s(k) = min(n0smax, n0s*exp(-alpha*temp_C))
3019             lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1))
3020             ilams(k) = 1./lams
3021             L_qs(k) = .true.
3022          else
3023             rs(k) = 1.E-12
3024             L_qs(k) = .false.
3025          endif
3026
3027          if (qg1d(k) .gt. 1.E-9) then
3028             rg(k) = qg1d(k)*rho(k)
3029             N0_g(k) = n0g
3030             lamg = (xam_g*xcgg(3)*N0_g(k)/rg(k))**(1./xcge(1))
3031             ilamg(k) = 1./lamg
3032             L_qg(k) = .true.
3033          else
3034             rg(k) = 1.E-12
3035             L_qg(k) = .false.
3036          endif
3037       enddo
3038
3039 !+---+-----------------------------------------------------------------+
3040 !..Locate K-level of start of melting (k_0 is level above).
3041 !+---+-----------------------------------------------------------------+
3042       melti = .false.
3043       k_0 = kts
3044       do k = kte-1, kts, -1
3045          if ( (temp(k).gt.273.15) .and. L_qr(k)                         &
3046                                   .and. (L_qs(k+1).or.L_qg(k+1)) ) then
3047             k_0 = MAX(k+1, k_0)
3048             melti=.true.
3049             goto 195
3050          endif
3051       enddo
3052  195  continue
3053
3054 !+---+-----------------------------------------------------------------+
3055 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
3056 !.. and non-water-coated snow and graupel when below freezing are
3057 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
3058 !+---+-----------------------------------------------------------------+
3059
3060       do k = kts, kte
3061          ze_rain(k) = 1.e-22
3062          ze_snow(k) = 1.e-22
3063          ze_graupel(k) = 1.e-22
3064          if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
3065          if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI)     &
3066                                  * (xam_s/900.0)*(xam_s/900.0)          &
3067                                  * N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
3068          if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI)  &
3069                                     * (xam_g/900.0)*(xam_g/900.0)       &
3070                                     * N0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
3071       enddo
3072
3073
3074 !+---+-----------------------------------------------------------------+
3075 !..Special case of melting ice (snow/graupel) particles.  Assume the
3076 !.. ice is surrounded by the liquid water.  Fraction of meltwater is
3077 !.. extremely simple based on amount found above the melting level.
3078 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
3079 !.. routines).
3080 !+---+-----------------------------------------------------------------+
3081
3082       if (melti .and. k_0.ge.kts+1) then
3083        do k = k_0-1, kts, -1
3084
3085 !..Reflectivity contributed by melting snow
3086           if (L_qs(k) .and. L_qs(k_0) ) then
3087            fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
3088            eta = 0.d0
3089            lams = 1./ilams(k)
3090            do n = 1, nrbins
3091               x = xam_s * xxDs(n)**xbm_s
3092               call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
3093                     fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
3094                     CBACK, mixingrulestring_s, matrixstring_s,          &
3095                     inclusionstring_s, hoststring_s,                    &
3096                     hostmatrixstring_s, hostinclusionstring_s)
3097               f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
3098               eta = eta + f_d * CBACK * simpson(n) * xdts(n)
3099            enddo
3100            ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
3101           endif
3102
3103
3104 !..Reflectivity contributed by melting graupel
3105
3106           if (L_qg(k) .and. L_qg(k_0) ) then
3107            fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
3108            eta = 0.d0
3109            lamg = 1./ilamg(k)
3110            do n = 1, nrbins
3111               x = xam_g * xxDg(n)**xbm_g
3112               call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), &
3113                     fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
3114                     CBACK, mixingrulestring_g, matrixstring_g,          &
3115                     inclusionstring_g, hoststring_g,                    &
3116                     hostmatrixstring_g, hostinclusionstring_g)
3117               f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n))
3118               eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
3119            enddo
3120            ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
3121           endif
3122
3123        enddo
3124       endif
3125
3126       do k = kte, kts, -1
3127          dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
3128       enddo
3129
3130
3131       end subroutine refl10cm_wdm6
3132 !+---+-----------------------------------------------------------------+
3133
3134 !-----------------------------------------------------------------------
3135      subroutine effectRad_wdm6 (t, qc, nc, qi, qs, rho, qmin, t0c,        &
3136                                 re_qc, re_qi, re_qs, kts, kte, ii, jj)
3137
3138 !-----------------------------------------------------------------------
3139 !  Compute radiation effective radii of cloud water, ice, and snow for
3140 !  double-moment microphysics..
3141 !  These are entirely consistent with microphysics assumptions, not
3142 !  constant or otherwise ad hoc as is internal to most radiation
3143 !  schemes.
3144 !  Coded and implemented by Soo Ya Bae, KIAPS, January 2015.
3145 !-----------------------------------------------------------------------
3146
3147       implicit none
3148
3149 !..Sub arguments
3150       integer, intent(in) :: kts, kte, ii, jj
3151       real, intent(in) :: qmin
3152       real, intent(in) :: t0c
3153       real, dimension( kts:kte ), intent(in)::  t
3154       real, dimension( kts:kte ), intent(in)::  qc
3155       real, dimension( kts:kte ), intent(in)::  nc
3156       real, dimension( kts:kte ), intent(in)::  qi
3157       real, dimension( kts:kte ), intent(in)::  qs
3158       real, dimension( kts:kte ), intent(in)::  rho
3159       real, dimension( kts:kte ), intent(inout):: re_qc
3160       real, dimension( kts:kte ), intent(inout):: re_qi
3161       real, dimension( kts:kte ), intent(inout):: re_qs
3162 !..Local variables
3163       integer:: i,k
3164       integer :: inu_c
3165       real, dimension( kts:kte ):: ni
3166       real, dimension( kts:kte ):: rqc
3167       real, dimension( kts:kte ):: rnc
3168       real, dimension( kts:kte ):: rqi
3169       real, dimension( kts:kte ):: rni
3170       real, dimension( kts:kte ):: rqs
3171       real :: cdm2
3172       real :: temp
3173       real :: supcol, n0sfac, lamdas
3174       real :: diai      ! diameter of ice in m
3175       double precision :: lamc
3176       logical :: has_qc, has_qi, has_qs
3177 !..Minimum microphys values
3178       real, parameter :: R1 = 1.E-12
3179       real, parameter :: R2 = 1.E-6
3180       real, parameter :: pi = 3.1415926536
3181       real, parameter :: bm_r = 3.0
3182       real, parameter :: obmr = 1.0/bm_r
3183       real, parameter :: cdm  = 5./3.
3184 !-----------------------------------------------------------------------
3185       has_qc = .false.
3186       has_qi = .false.
3187       has_qs = .false.
3188
3189       cdm2 = rgmma(cdm)
3190
3191       do k=kts,kte
3192         ! for cloud
3193         rqc(k) = max(R1, qc(k)*rho(k))
3194         rnc(k) = max(R2, nc(k)*rho(k))
3195         if (rqc(k).gt.R1 .and. rnc(k).gt.R2) has_qc = .true.
3196         ! for ice
3197         rqi(k) = max(R1, qi(k)*rho(k))
3198         temp = (rho(k)*max(qi(k),qmin))
3199         temp = sqrt(sqrt(temp*temp*temp))
3200         ni(k) = min(max(5.38e7*temp,1.e3),1.e6)
3201         rni(k)= max(R2, ni(k)*rho(k))
3202         if (rqi(k).gt.R1 .and. rni(k).gt.R2) has_qi = .true.
3203         ! for snow
3204         rqs(k) = max(R1, qs(k)*rho(k))
3205         if (rqs(k).gt.R1) has_qs = .true.
3206       enddo
3207
3208       if (has_qc) then
3209         do k=kts,kte
3210           if (rqc(k).le.R1 .or. rnc(k).le.R2) CYCLE
3211           lamc = 2.*cdm2*(pidnc*nc(k)/rqc(k))**obmr
3212           re_qc(k) =  max(2.51e-6,min(sngl(1.0d0/lamc),50.e-6))
3213         enddo
3214       endif
3215
3216       if (has_qi) then
3217         do k=kts,kte
3218           if (rqi(k).le.R1 .or. rni(k).le.R2) CYCLE
3219           diai = 11.9*sqrt(rqi(k)/ni(k))
3220           re_qi(k) = max(10.01e-6,min(0.75*0.163*diai,125.e-6))
3221         enddo
3222       endif
3223
3224       if (has_qs) then
3225         do k=kts,kte
3226           if (rqs(k).le.R1) CYCLE
3227           supcol = t0c-t(k)
3228           n0sfac = max(min(exp(alpha*supcol),n0smax/n0s),1.)
3229           lamdas = sqrt(sqrt(pidn0s*n0sfac/rqs(k)))
3230           re_qs(k) = max(25.e-6,min(0.5*(1./lamdas),999.e-6))
3231         enddo
3232       endif
3233
3234       end subroutine effectRad_wdm6
3235 !-----------------------------------------------------------------------
3236
3237 END MODULE module_mp_wdm6