phys/module_mp_wdm5.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 !Including inline expansion statistical function
  10 MODULE module_mp_wdm5
  11 !
  12    USE module_mp_radar
  13    USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG
  14 !
  15    REAL, PARAMETER, PRIVATE :: dtcldcr     = 120. ! maximum time step for minor loops
  16    REAL, PARAMETER, PRIVATE :: n0r = 8.e6         ! intercept parameter rain
  17    REAL, PARAMETER, PRIVATE :: avtr = 841.9       ! a constant for terminal velocity of rain
  18    REAL, PARAMETER, PRIVATE :: bvtr = 0.8         ! a constant for terminal velocity of rain
  19    REAL, PARAMETER, PRIVATE :: r0 = .8e-5         ! 8 microm  in contrast to 10 micro m
  20    REAL, PARAMETER, PRIVATE :: peaut = .55        ! collection efficiency
  21    REAL, PARAMETER, PRIVATE :: xncr = 3.e8        ! maritime cloud in contrast to 3.e8 in tc80
  22    REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5    ! the dynamic viscosity kgm-1s-1
  23    REAL, PARAMETER, PRIVATE :: avts = 11.72       ! a constant for terminal velocity of snow
  24    REAL, PARAMETER, PRIVATE :: bvts = .41         ! a constant for terminal velocity of snow
  25    REAL, PARAMETER, PRIVATE :: n0smax =  1.e11    ! maximum n0s (t=-90C unlimited)
  26    REAL, PARAMETER, PRIVATE :: lamdacmax = 5.0e5  ! limited maximum value for slope parameter of cloud water
  27    REAL, PARAMETER, PRIVATE :: lamdacmin = 2.0e4  ! limited minimum value for slope parameter of cloud water
  28    REAL, PARAMETER, PRIVATE :: lamdarmax = 5.0e4  ! limited maximum value for slope parameter of rain
  29    REAL, PARAMETER, PRIVATE :: lamdarmin = 2.0e3  ! limited minimum value for slope parameter of rain
  30    REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5   ! limited maximum value for slope parameter of snow
  31    REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4   ! limited maximum value for slope parameter of graupel
  32    REAL, PARAMETER, PRIVATE :: dicon = 11.9       ! constant for the cloud-ice diamter
  33    REAL, PARAMETER, PRIVATE :: dimax = 500.e-6    ! limited maximum value for the cloud-ice diamter
  34    REAL, PARAMETER, PRIVATE :: n0s = 2.e6         ! temperature dependent intercept parameter snow
  35    REAL, PARAMETER, PRIVATE :: alpha = .12        ! .122 exponen factor for n0s
  36    REAL, PARAMETER, PRIVATE :: pfrz1 = 100.       ! constant in Biggs freezing
  37    REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66       ! constant in Biggs freezing
  38    REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9     ! minimun values for qr, qs, and qg
  39    REAL, PARAMETER, PRIVATE :: ncmin = 1.e1       ! minimum value for Nc
  40    REAL, PARAMETER, PRIVATE :: nrmin = 1.e-2      ! minimum value for Nr
  41    REAL, PARAMETER, PRIVATE :: eacrc = 1.0        ! Snow/cloud-water collection efficiency
  42 !
  43    REAL, PARAMETER, PRIVATE :: satmax = 1.0048    ! maximum saturation value for CCN activation
  44                                                   ! 1.008 for maritime air mass /1.0048 for conti
  45    REAL, PARAMETER, PRIVATE :: actk = 0.6         ! parameter for the CCN activation
  46    REAL, PARAMETER, PRIVATE :: actr = 1.5         ! radius of activated CCN drops
  47    REAL, PARAMETER, PRIVATE :: ncrk1 = 3.03e3     ! Long's collection kernel coefficient
  48    REAL, PARAMETER, PRIVATE :: ncrk2 = 2.59e15    ! Long's collection kernel coefficient
  49    REAL, PARAMETER, PRIVATE :: di100 = 1.e-4      ! parameter related with accretion and collection of cloud drops
  50    REAL, PARAMETER, PRIVATE :: di600 = 6.e-4      ! parameter related with accretion and collection of cloud drops
  51    REAL, PARAMETER, PRIVATE :: di2000 = 20.e-4    ! parameter related with accretion and collection of cloud drops
  52    REAL, PARAMETER, PRIVATE :: di82 = 82.e-6      ! dimater related with raindrops evaporation
  53    REAL, PARAMETER, PRIVATE :: di15 = 15.e-6      ! auto conversion takes place beyond this diameter
  54    REAL, SAVE ::                                      &
  55              qc0, qck1,pidnc,bvtr1,bvtr2,bvtr3,bvtr4, &
  56              bvtr5,bvtr7,bvtr2o5,bvtr3o5,g1pbr,g2pbr, &
  57              g3pbr,g4pbr,g5pbr,g7pbr,g5pbro2,g7pbro2, &
  58              pvtr,pvtrn,eacrr,pacrr, pi,              &
  59              precr1,precr2,xmmax,roqimax,bvts1,       &
  60              bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs,     &
  61              g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
  62              pidn0s,pidnr,xlv1,pacrc,                 &
  63              rslopecmax,rslopec2max,rslopec3max,      &
  64              rslopermax,rslopesmax,rslopegmax,        &
  65              rsloperbmax,rslopesbmax,rslopegbmax,     &
  66              rsloper2max,rslopes2max,rslopeg2max,     &
  67              rsloper3max,rslopes3max,rslopeg3max
  68 !
  69 ! Specifies code-inlining of fpvs function in WDM52D below. JM 20040507
  70 !
  71 CONTAINS
  72 !===================================================================
  73 !
  74   SUBROUTINE wdm5(th, q, qc, qr, qi, qs                            &
  75                  ,nn, nc, nr                                       &
  76                  ,den, pii, p, delz                                &
  77                  ,delt,g, cpd, cpv, ccn0, rd, rv, t0c              &
  78                  ,ep1, ep2, qmin                                   &
  79                  ,XLS, XLV0, XLF0, den0, denr                      &
  80                  ,cliq,cice,psat                                   &
  81                  ,rain, rainncv                                    &
  82                  ,snow, snowncv                                    &
  83                  ,sr                                               &
  84                  ,has_reqc, has_reqi, has_reqs                     &  ! for radiation
  85                  ,re_cloud, re_ice,   re_snow                      &  ! for radiation
  86                  ,refl_10cm, diagflag, do_radar_ref                &
  87                  ,itimestep                                        &
  88                  ,ids,ide, jds,jde, kds,kde                        &
  89                  ,ims,ime, jms,jme, kms,kme                        &
  90                  ,its,ite, jts,jte, kts,kte                        &
  91                                                                    )
  92 !-------------------------------------------------------------------
  93   IMPLICIT NONE
  94 !-------------------------------------------------------------------
  95 !
  96 !  This code is a WRF double-moment 5-class  mixed ice
  97 !  microphyiscs scheme (WDM5). The WDM microphysics scheme predicts
  98 !  number concentrations for warm rain species including clouds and
  99 !  rain. cloud condensation nuclei (CCN) is also predicted.
 100 !  The cold rain species including ice, snow, graupel follow the
 101 !  WRF single-moment 5-class microphysics (WSM5)
 102 !  in which theoretical background for WSM ice phase microphysics is
 103 !  based on Hong et al. (2004).
 104 !  The WDM scheme is described in Lim and Hong (2010).
 105 !  All units are in m.k.s. and source/sink terms in kgkg-1s-1.
 106 !
 107 !  WDM5 cloud scheme
 108 !
 109 !  Coded by Kyo-Sun Lim and Song-You Hong (Yonsei Univ.) Fall 2008
 110 !
 111 !  Implemented by Kyo-Sun Lim and Jimy Dudhia (NCAR) Winter 2008
 112 !
 113 !  further modifications :
 114 !        semi-lagrangian sedimentation (JH,2010),hong, aug 2009
 115 !        ==> higher accuracy and efficient at lower resolutions
 116 !        reflectivity computation from greg thompson, lim, jun 2011
 117 !        ==> only diagnostic, but with removal of too large drops
 118 !        effective radius of hydrometeors, bae from kiaps, jan 2015
 119 !        ==> consistency in solar insolation of rrtmg radiation
 120 !        bug fix in melting terms, bae from kiaps, nov 2015
 121 !        ==> density of air is divided, which has not been
 122 !
 123 !  Reference) Lim and Hong (LH, 2010) Mon. Wea. Rev.
 124 !             Juang and Hong (JH, 2010) Mon. Wea. Rev.
 125 !             Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
 126 !             Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
 127 !             Cohard and Pinty (CP, 2000) Quart. J. Roy. Meteor. Soc.
 128 !             Khairoutdinov and Kogan (KK, 2000) Mon. Wea. Rev.
 129 !             Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
 130 !
 131 !             Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
 132 !             Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
 133 !             Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
 134 !
 135   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
 136                                       ims,ime, jms,jme, kms,kme , &
 137                                       its,ite, jts,jte, kts,kte
 138   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
 139         INTENT(INOUT) ::                                          &
 140                                                              th,  &
 141                                                               q,  &
 142                                                               qc, &
 143                                                               qi, &
 144                                                               qr, &
 145                                                               qs, &
 146                                                               nn, &
 147                                                               nc, &
 148                                                               nr
 149   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
 150         INTENT(IN   ) ::                                          &
 151                                                              den, &
 152                                                              pii, &
 153                                                                p, &
 154                                                             delz
 155   REAL, INTENT(IN   ) ::                                    delt, &
 156                                                                g, &
 157                                                               rd, &
 158                                                               rv, &
 159                                                              t0c, &
 160                                                             den0, &
 161                                                              cpd, &
 162                                                              cpv, &
 163                                                             ccn0, &
 164                                                              ep1, &
 165                                                              ep2, &
 166                                                             qmin, &
 167                                                              XLS, &
 168                                                             XLV0, &
 169                                                             XLF0, &
 170                                                             cliq, &
 171                                                             cice, &
 172                                                             psat, &
 173                                                             denr
 174   INTEGER, INTENT(IN   ) ::                            itimestep
 175   REAL, DIMENSION( ims:ime , jms:jme ),                           &
 176         INTENT(INOUT) ::                                    rain, &
 177                                                          rainncv, &
 178                                                               sr
 179 ! for radiation connecting
 180   INTEGER, INTENT(IN)::                                           &
 181                                                         has_reqc, &
 182                                                         has_reqi, &
 183                                                         has_reqs
 184   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),                     &
 185         INTENT(INOUT)::                                           &
 186                                                         re_cloud, &
 187                                                           re_ice, &
 188                                                          re_snow
 189 !+---+-----------------------------------------------------------------+
 190   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT)::     &  ! GT
 191                                                        refl_10cm
 192 !+---+-----------------------------------------------------------------+
 193
 194   REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                 &
 195         INTENT(INOUT) ::                                    snow, &
 196                                                          snowncv
 197 ! LOCAL VAR
 198   REAL, DIMENSION( its:ite , kts:kte ) ::   t
 199   REAL, DIMENSION( its:ite , kts:kte, 2 ) ::   qci, qrs
 200   REAL, DIMENSION( its:ite , kts:kte, 3 ) ::   ncr
 201   CHARACTER*256 :: emess
 202   INTEGER :: mkx_test
 203   INTEGER ::               i,j,k
 204
 205 !+---+-----------------------------------------------------------------+
 206       REAL, DIMENSION(kts:kte):: qv1d, t1d, p1d, qr1d, nr1d, qs1d, dBZ
 207       LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
 208       INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
 209 !+---+-----------------------------------------------------------------+
 210 ! to calculate effective radius for radiation
 211   REAL, DIMENSION( kts:kte ) :: qc1d, nc1d, den1d
 212   REAL, DIMENSION( kts:kte ) :: qi1d
 213   REAL, DIMENSION( kts:kte ) :: re_qc, re_qi, re_qs
 214
 215 #ifndef RUN_ON_GPU
 216       IF (itimestep .eq. 1) THEN
 217         DO j=jms,jme
 218            DO k=kms,kme
 219            DO i=ims,ime
 220               nn(i,k,j) = ccn0
 221            ENDDO
 222            ENDDO
 223         ENDDO
 224       ENDIF
 225 !
 226       DO j=jts,jte
 227          DO k=kts,kte
 228          DO i=its,ite
 229             t(i,k)=th(i,k,j)*pii(i,k,j)
 230             qci(i,k,1) = qc(i,k,j)
 231             qci(i,k,2) = qi(i,k,j)
 232             qrs(i,k,1) = qr(i,k,j)
 233             qrs(i,k,2) = qs(i,k,j)
 234             ncr(i,k,1) = nn(i,k,j)
 235             ncr(i,k,2) = nc(i,k,j)
 236             ncr(i,k,3) = nr(i,k,j)
 237          ENDDO
 238          ENDDO
 239          !  Sending array starting locations of optional variables may cause
 240          !  troubles, so we explicitly change the call.
 241          CALL wdm52D(t, q(ims,kms,j), qci, qrs, ncr               &
 242                     ,den(ims,kms,j)                               &
 243                     ,p(ims,kms,j), delz(ims,kms,j)                &
 244                     ,delt,g, cpd, cpv, ccn0, rd, rv, t0c          &
 245                     ,ep1, ep2, qmin                               &
 246                     ,XLS, XLV0, XLF0, den0, denr                  &
 247                     ,cliq,cice,psat                               &
 248                     ,j                                            &
 249                     ,rain(ims,j),rainncv(ims,j)                   &
 250                     ,sr(ims,j)                                    &
 251                     ,ids,ide, jds,jde, kds,kde                    &
 252                     ,ims,ime, jms,jme, kms,kme                    &
 253                     ,its,ite, jts,jte, kts,kte                    &
 254                     ,snow(ims,j),snowncv(ims,j)                   &
 255                                                                   )
 256          DO K=kts,kte
 257          DO I=its,ite
 258             th(i,k,j)=t(i,k)/pii(i,k,j)
 259             qc(i,k,j) = qci(i,k,1)
 260             qi(i,k,j) = qci(i,k,2)
 261             qr(i,k,j) = qrs(i,k,1)
 262             qs(i,k,j) = qrs(i,k,2)
 263             nn(i,k,j) = ncr(i,k,1)
 264             nc(i,k,j) = ncr(i,k,2)
 265             nr(i,k,j) = ncr(i,k,3)
 266          ENDDO
 267          ENDDO
 268
 269 !+---+-----------------------------------------------------------------+
 270          IF ( PRESENT (diagflag) ) THEN
 271          if (diagflag .and. do_radar_ref == 1) then
 272             DO I=its,ite
 273                DO K=kts,kte
 274                   t1d(k)=th(i,k,j)*pii(i,k,j)
 275                   p1d(k)=p(i,k,j)
 276                   qv1d(k)=q(i,k,j)
 277                   qr1d(k)=qr(i,k,j)
 278                   nr1d(k)=nr(i,k,j)
 279                   qs1d(k)=qs(i,k,j)
 280                ENDDO
 281                call refl10cm_wdm5 (qv1d, qr1d, nr1d, qs1d,              &
 282                        t1d, p1d, dBZ, kts, kte, i, j)
 283                do k = kts, kte
 284                   refl_10cm(i,k,j) = MAX(-35., dBZ(k))
 285                enddo
 286             ENDDO
 287          endif
 288          ENDIF
 289
 290         if (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) then
 291           do i = its, ite
 292             do k = kts, kte
 293               re_qc(k) = RE_QC_BG
 294               re_qi(k) = RE_QI_BG
 295               re_qs(k) = RE_QS_BG
 296
 297               t1d(k)  = th(i,k,j)*pii(i,k,j)
 298               den1d(k)= den(i,k,j)
 299               qc1d(k) = qc(i,k,j)
 300               qi1d(k) = qi(i,k,j)
 301               qs1d(k) = qs(i,k,j)
 302               nc1d(k) = nc(i,k,j)
 303             enddo
 304             call effectRad_wdm5 (t1d, qc1d, nc1d, qi1d, qs1d, den1d,     &
 305                                 qmin, t0c, re_qc, re_qi, re_qs,         &
 306                                 kts, kte, i, j)
 307             do k = kts, kte
 308               re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc(k),  50.E-6))
 309               re_ice(i,k,j)   = MAX(RE_QI_BG, MIN(re_qi(k), 125.E-6))
 310               re_snow(i,k,j)  = MAX(RE_QS_BG, MIN(re_qs(k), 999.E-6))
 311             enddo ! k loop
 312           enddo   ! i loop
 313         endif     ! has_reqc, etc...
 314 !+---+-----------------------------------------------------------------+
 315
 316       ENDDO
 317 #else
 318       CALL get_wsm5_gpu_levels ( mkx_test )
 319       IF ( mkx_test .LT. kte ) THEN
 320         WRITE(emess,*)'Number of levels compiled for GPU WSM5 too small. ',    &
 321                       mkx_test,' < ',kte
 322         CALL wrf_error_fatal(emess)
 323       ENDIF
 324       CALL wsm5_host (                                                         &
 325                     th(its:ite,kts:kte,jts:jte), pii(its:ite,kts:kte,jts:jte)  &
 326                    ,q(its:ite,kts:kte,jts:jte), qc(its:ite,kts:kte,jts:jte)    &
 327                    ,qi(its:ite,kts:kte,jts:jte), qr(its:ite,kts:kte,jts:jte)   &
 328                    ,qs(its:ite,kts:kte,jts:jte), den(its:ite,kts:kte,jts:jte)  &
 329                    ,p(its:ite,kts:kte,jts:jte), delz(its:ite,kts:kte,jts:jte)  &
 330                    ,delt                                                       &
 331                    ,rain(its:ite,jts:jte),rainncv(its:ite,jts:jte)             &
 332                    ,snow(its:ite,jts:jte),snowncv(its:ite,jts:jte)             &
 333                    ,sr(its:ite,jts:jte)                                        &
 334                    ,its, ite,  jts, jte,  kts, kte                             &
 335                    ,its, ite,  jts, jte,  kts, kte                             &
 336                    ,its, ite,  jts, jte,  kts, kte                             &
 337           )
 338 #endif
 339   END SUBROUTINE wdm5
 340 !===================================================================
 341 !
 342   SUBROUTINE wdm52D(t, q, qci, qrs, ncr, den, p, delz             &
 343                    ,delt,g, cpd, cpv, ccn0, rd, rv, t0c           &
 344                    ,ep1, ep2, qmin                                &
 345                    ,XLS, XLV0, XLF0, den0, denr                   &
 346                    ,cliq,cice,psat                                &
 347                    ,lat                                           &
 348                    ,rain,rainncv                                  &
 349                    ,sr                                            &
 350                    ,ids,ide, jds,jde, kds,kde                     &
 351                    ,ims,ime, jms,jme, kms,kme                     &
 352                    ,its,ite, jts,jte, kts,kte                     &
 353                    ,snow,snowncv                                  &
 354                                                                   )
 355 !-------------------------------------------------------------------
 356   IMPLICIT NONE
 357 !-------------------------------------------------------------------
 358   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
 359                                       ims,ime, jms,jme, kms,kme , &
 360                                       its,ite, jts,jte, kts,kte,  &
 361                                       lat
 362   REAL, DIMENSION( its:ite , kts:kte ),                           &
 363         INTENT(INOUT) ::                                          &
 364                                                                t
 365   REAL, DIMENSION( its:ite , kts:kte, 2 ),                        &
 366         INTENT(INOUT) ::                                          &
 367                                                              qci, &
 368                                                              qrs
 369   REAL, DIMENSION( its:ite , kts:kte, 3 ),                        &
 370         INTENT(INOUT) ::                                          &
 371                                                              ncr
 372   REAL, DIMENSION( ims:ime , kms:kme ),                           &
 373         INTENT(INOUT) ::                                          &
 374                                                                q
 375   REAL, DIMENSION( ims:ime , kms:kme ),                           &
 376         INTENT(IN   ) ::                                          &
 377                                                              den, &
 378                                                                p, &
 379                                                             delz
 380   REAL, INTENT(IN   ) ::                                    delt, &
 381                                                                g, &
 382                                                              cpd, &
 383                                                              cpv, &
 384                                                             ccn0, &
 385                                                              t0c, &
 386                                                             den0, &
 387                                                               rd, &
 388                                                               rv, &
 389                                                              ep1, &
 390                                                              ep2, &
 391                                                             qmin, &
 392                                                              XLS, &
 393                                                             XLV0, &
 394                                                             XLF0, &
 395                                                             cliq, &
 396                                                             cice, &
 397                                                             psat, &
 398                                                             denr
 399   REAL, DIMENSION( ims:ime ),                                     &
 400         INTENT(INOUT) ::                                    rain, &
 401                                                          rainncv, &
 402                                                               sr
 403   REAL, DIMENSION( ims:ime ),     OPTIONAL,                       &
 404         INTENT(INOUT) ::                                    snow, &
 405                                                          snowncv
 406 ! LOCAL VAR
 407   REAL, DIMENSION( its:ite , kts:kte , 2) ::                      &
 408         rh, qs, rslope, rslope2, rslope3, rslopeb,                &
 409         falk, fall, work1, qrs_tmp
 410   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
 411         rslopec, rslopec2,rslopec3
 412   REAL, DIMENSION( its:ite , kts:kte,  2) ::                      &
 413         avedia
 414   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
 415         workn,falln,falkn
 416   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
 417         works
 418   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
 419         den_tmp, delz_tmp, ncr_tmp
 420   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
 421         lamdr_tmp
 422   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
 423         lamdc_tmp
 424   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
 425               falkc, work1c, work2c, fallc
 426   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
 427         pcact, praut, psaut, prevp, psdep, pracw, psaci, psacw,   &
 428         pigen, pidep, pcond,                                      &
 429         xl, cpm, work2, psmlt, psevp, denfac, xni,                &
 430         n0sfac, denqrs2, denqci
 431   REAL, DIMENSION( its:ite ) ::                                   &
 432         delqrs2, delqi
 433   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
 434         nraut, nracw, ncevp, nccol, nrcol,                        &
 435         nsacw, nseml, ncact
 436   REAL :: ifac, sfac
 437   REAL, DIMENSION(its:ite) :: tstepsnow
 438 !
 439 #define WSM_NO_CONDITIONAL_IN_VECTOR
 440 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
 441   REAL, DIMENSION(its:ite) :: xal, xbl
 442 #endif
 443 ! variables for optimization
 444   REAL, DIMENSION( its:ite )           :: tvec1
 445   INTEGER, DIMENSION( its:ite ) :: mnstep, numndt
 446   INTEGER, DIMENSION( its:ite ) :: mstep, numdt
 447   REAL, DIMENSION(its:ite) :: rmstep
 448   REAL dtcldden, rdelz, rdtcld
 449   LOGICAL, DIMENSION( its:ite ) :: flgcld
 450   REAL  ::                                                        &
 451             cpmcal, xlcal, lamdac, diffus,                        &
 452             viscos, xka, venfac, conden, diffac,                  &
 453             x, y, z, a, b, c, d, e,                               &
 454             ndt, qdt, holdrr, holdrs, supcol, supcolt, pvt,       &
 455             coeres, supsat, dtcld, xmi, eacrs, satdt,             &
 456             vt2i,vt2s,acrfac, coecol,                             &
 457             nfrzdtr, nfrzdtc,                                     &
 458             taucon, lencon, lenconcr,                             &
 459             qimax, diameter, xni0, roqi0,                         &
 460             fallsum, fallsum_qsi, xlwork2, factor, source,        &
 461             value, xlf, pfrzdtc, pfrzdtr, supice
 462   REAL :: temp
 463   REAL  :: holdc, holdci
 464   INTEGER :: i, j, k, mstepmax,                                                &
 465             iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
 466 ! Temporaries used for inlining fpvs function
 467   REAL  :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
 468   REAL  :: logtr
 469 !
 470 !=================================================================
 471 !   compute internal functions
 472 !
 473       cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
 474       xlcal(x) = xlv0-xlv1*(x-t0c)
 475 !----------------------------------------------------------------
 476 !     size distributions: (x=mixing ratio, y=air density):
 477 !     valid for mixing ratio > 1.e-9 kg/kg.
 478 !
 479 ! Optimizatin : A**B => exp(log(A)*(B))
 480       lamdac(x,y,z)= exp(log(((pidnc*z)/(x*y)))*((.33333333)))
 481 !
 482 !----------------------------------------------------------------
 483 !     diffus: diffusion coefficient of the water vapor
 484 !     viscos: kinematic viscosity(m2s-1)
 485 !     diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y
 486 !     viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y
 487 !     xka(x,y) = 1.414e3*viscos(x,y)*y
 488 !     diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
 489 !     venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333)))         &
 490 !                    /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
 491 !     conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
 492 !
 493 !
 494       idim = ite-its+1
 495       kdim = kte-kts+1
 496 !
 497 !----------------------------------------------------------------
 498 !     paddint 0 for negative values generated by dynamics
 499 !
 500       do k = kts, kte
 501         do i = its, ite
 502           qci(i,k,1) = max(qci(i,k,1),0.0)
 503           qrs(i,k,1) = max(qrs(i,k,1),0.0)
 504           qci(i,k,2) = max(qci(i,k,2),0.0)
 505           qrs(i,k,2) = max(qrs(i,k,2),0.0)
 506           ncr(i,k,1) = max(ncr(i,k,1),0.)
 507           ncr(i,k,2) = max(ncr(i,k,2),0.)
 508           ncr(i,k,3) = max(ncr(i,k,3),0.)
 509         enddo
 510       enddo
 511 !
 512 !     latent heat for phase changes and heat capacity. neglect the
 513 !     changes during microphysical process calculation
 514 !     emanuel(1994)
 515 !
 516       do k = kts, kte
 517         do i = its, ite
 518           cpm(i,k) = cpmcal(q(i,k))
 519           xl(i,k) = xlcal(t(i,k))
 520           delz_tmp(i,k) = delz(i,k)
 521           den_tmp(i,k) = den(i,k)
 522         enddo
 523       enddo
 524 !
 525 !----------------------------------------------------------------
 526 !    initialize the surface rain, snow
 527 !
 528       do i = its, ite
 529         rainncv(i) = 0.
 530         if(PRESENT (snowncv) .AND. PRESENT (snow)) snowncv(i) = 0.
 531         sr(i) = 0.
 532 ! new local array to catch step snow
 533         tstepsnow(i) = 0.
 534       enddo
 535 !
 536 !----------------------------------------------------------------
 537 !     compute the minor time steps.
 538 !
 539       loops = max(nint(delt/dtcldcr),1)
 540       dtcld = delt/loops
 541       if(delt.le.dtcldcr) dtcld = delt
 542 !
 543       do loop = 1,loops
 544 !
 545 !----------------------------------------------------------------
 546 !     initialize the large scale variables
 547 !
 548       do i = its, ite
 549         mstep(i) = 1
 550         mnstep(i) = 1
 551         flgcld(i) = .true.
 552       enddo
 553 !
 554 !     do k = kts, kte
 555 !       do i = its, ite
 556 !         denfac(i,k) = sqrt(den0/den(i,k))
 557 !       enddo
 558 !     enddo
 559       do k = kts, kte
 560         CALL VREC( tvec1(its), den(its,k), ite-its+1)
 561         do i = its, ite
 562           tvec1(i) = tvec1(i)*den0
 563         enddo
 564         CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
 565       enddo
 566 !
 567 ! Inline expansion for fpvs
 568 !         qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
 569 !         qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
 570       hsub = xls
 571       hvap = xlv0
 572       cvap = cpv
 573       ttp=t0c+0.01
 574       dldt=cvap-cliq
 575       xa=-dldt/rv
 576       xb=xa+hvap/(rv*ttp)
 577       dldti=cvap-cice
 578       xai=-dldti/rv
 579       xbi=xai+hsub/(rv*ttp)
 580 ! this is for compilers where the conditional inhibits vectorization
 581 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
 582       do k = kts, kte
 583         do i = its, ite
 584           if(t(i,k).lt.ttp) then
 585             xal(i) = xai
 586             xbl(i) = xbi
 587           else
 588             xal(i) = xa
 589             xbl(i) = xb
 590           endif
 591         enddo
 592         do i = its, ite
 593           tr=ttp/t(i,k)
 594           logtr=log(tr)
 595           qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
 596           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
 597           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
 598           qs(i,k,1) = max(qs(i,k,1),qmin)
 599           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
 600           qs(i,k,2)=psat*exp(logtr*(xal(i))+xbl(i)*(1.-tr))
 601           qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
 602           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
 603           qs(i,k,2) = max(qs(i,k,2),qmin)
 604           rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
 605         enddo
 606       enddo
 607 #else
 608       do k = kts, kte
 609         do i = its, ite
 610           tr=ttp/t(i,k)
 611           logtr=log(tr)
 612           qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
 613           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
 614           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
 615           qs(i,k,1) = max(qs(i,k,1),qmin)
 616           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
 617           if(t(i,k).lt.ttp) then
 618             qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
 619           else
 620             qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
 621           endif
 622           qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
 623           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
 624           qs(i,k,2) = max(qs(i,k,2),qmin)
 625           rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
 626         enddo
 627       enddo
 628 #endif
 629 !
 630 !----------------------------------------------------------------
 631 !     initialize the variables for microphysical physics
 632 !
 633 !
 634       do k = kts, kte
 635         do i = its, ite
 636           prevp(i,k) = 0.
 637           psdep(i,k) = 0.
 638           praut(i,k) = 0.
 639           psaut(i,k) = 0.
 640           pracw(i,k) = 0.
 641           psaci(i,k) = 0.
 642           psacw(i,k) = 0.
 643           pigen(i,k) = 0.
 644           pidep(i,k) = 0.
 645           pcond(i,k) = 0.
 646           psmlt(i,k) = 0.
 647           psevp(i,k) = 0.
 648           pcact(i,k) = 0.
 649           falk(i,k,1) = 0.
 650           falk(i,k,2) = 0.
 651           fall(i,k,1) = 0.
 652           fall(i,k,2) = 0.
 653           fallc(i,k) = 0.
 654           falkc(i,k) = 0.
 655           falln(i,k) = 0.
 656           falkn(i,k) = 0.
 657           xni(i,k) = 1.e3
 658           nsacw(i,k) = 0.
 659           nseml(i,k) = 0.
 660           nracw(i,k) = 0.
 661           nccol(i,k) = 0.
 662           nrcol(i,k) = 0.
 663           ncact(i,k) = 0.
 664           nraut(i,k) = 0.
 665           ncevp(i,k) = 0.
 666         enddo
 667       enddo
 668 !
 669 !----------------------------------------------------------------
 670 !     compute the fallout term:
 671 !     first, vertical terminal velosity for minor loops
 672 !
 673       do k = kts, kte
 674         do i = its, ite
 675           if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin)then
 676             rslopec(i,k) = rslopecmax
 677             rslopec2(i,k) = rslopec2max
 678             rslopec3(i,k) = rslopec3max
 679           else
 680             rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
 681             rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
 682             rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
 683           endif
 684 !-------------------------------------------------------------
 685 ! Ni: ice crystal number concentraiton   [HDC 5c]
 686 !-------------------------------------------------------------
 687 !         xni(i,k) = min(max(5.38e7*(den(i,k)                                  &
 688 !                   *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
 689           temp = (den(i,k)*max(qci(i,k,2),qmin))
 690           temp = sqrt(sqrt(temp*temp*temp))
 691           xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
 692         enddo
 693       enddo
 694       do k = kts, kte
 695         do i = its, ite
 696           qrs_tmp(i,k,1) = qrs(i,k,1)
 697           qrs_tmp(i,k,2) = qrs(i,k,2)
 698           ncr_tmp(i,k) = ncr(i,k,3)
 699         enddo
 700       enddo
 701       call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
 702                      rslope3,work1,workn,its,ite,kts,kte)
 703 !----------------------------------------------------------------
 704 !     compute the fallout term:
 705 !     first, vertical terminal velosity for minor loops
 706 !----------------------------------------------------------------
 707 !
 708 ! vt update for qr and nr
 709       mstepmax = 1
 710       numdt = 1
 711       do k = kte, kts, -1
 712         do i = its, ite
 713           work1(i,k,1) = work1(i,k,1)/delz(i,k)
 714           workn(i,k) = workn(i,k)/delz(i,k)
 715           numdt(i) = max(nint(max(work1(i,k,1),workn(i,k))*dtcld+.5),1)
 716           if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
 717         enddo
 718       enddo
 719       do i = its, ite
 720         if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
 721       enddo
 722 !
 723       do n = 1, mstepmax
 724         k = kte
 725         do i = its, ite
 726           if(n.le.mstep(i)) then
 727             falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
 728             falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
 729             fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
 730             falln(i,k) = falln(i,k)+falkn(i,k)
 731             qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
 732             ncr(i,k,3) = max(ncr(i,k,3)-falkn(i,k)*dtcld,0.)
 733           endif
 734         enddo
 735         do k = kte-1, kts, -1
 736           do i = its, ite
 737             if(n.le.mstep(i)) then
 738               falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
 739               falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
 740               fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
 741               falln(i,k) = falln(i,k)+falkn(i,k)
 742               qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1)           &
 743                           *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
 744               ncr(i,k,3) = max(ncr(i,k,3)-(falkn(i,k)-falkn(i,k+1)*delz(i,k+1) &
 745                           /delz(i,k))*dtcld,0.)
 746             endif
 747           enddo
 748         enddo
 749         do k = kts, kte
 750           do i = its, ite
 751             qrs_tmp(i,k,1) = qrs(i,k,1)
 752             ncr_tmp(i,k) = ncr(i,k,3)
 753           enddo
 754         enddo
 755         call slope_rain(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
 756                      rslope3,work1,workn,its,ite,kts,kte)
 757         do k = kte, kts, -1
 758           do i = its, ite
 759             work1(i,k,1) = work1(i,k,1)/delz(i,k)
 760             workn(i,k) = workn(i,k)/delz(i,k)
 761           enddo
 762         enddo
 763       enddo
 764 ! for semi
 765       do k = kte, kts, -1
 766         do i = its, ite
 767           works(i,k) = work1(i,k,2)
 768           denqrs2(i,k) = den(i,k)*qrs(i,k,2)
 769           if(qrs(i,k,2).le.0.0) works(i,k) = 0.0
 770         enddo
 771       enddo
 772       call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,works,denqrs2,  &
 773                            delqrs2,dtcld,2,1)
 774       do k = kts, kte
 775         do i = its, ite
 776           qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
 777           fall(i,k,2) = denqrs2(i,k)*works(i,k)/delz(i,k)
 778         enddo
 779       enddo
 780       do i = its, ite
 781         fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
 782       enddo
 783       do k = kts, kte
 784         do i = its, ite
 785           qrs_tmp(i,k,1) = qrs(i,k,1)
 786           qrs_tmp(i,k,2) = qrs(i,k,2)
 787           ncr_tmp(i,k) = ncr(i,k,3)
 788         enddo
 789       enddo
 790       call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
 791                      rslope3,work1,workn,its,ite,kts,kte)
 792 !
 793       do k = kte, kts, -1
 794         do i = its, ite
 795           supcol = t0c-t(i,k)
 796           n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
 797           if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
 798 !----------------------------------------------------------------
 799 ! psmlt: melting of snow [HL A33] [RH83 A25]
 800 !       (T>T0: QS->QR)
 801 !----------------------------------------------------------------
 802             xlf = xlf0
 803 !           work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
 804             work2(i,k)= (exp(log(((1.496e-6*((t(i,k))*sqrt(t(i,k)))        &
 805                         /((t(i,k))+120.)/(den(i,k)))/(8.794e-5             &
 806                         *exp(log(t(i,k))*(1.81))/p(i,k))))                 &
 807                         *((.3333333)))/sqrt((1.496e-6*((t(i,k))            &
 808                         *sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k))))        &
 809                         *sqrt(sqrt(den0/(den(i,k)))))
 810             coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
 811 !           psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2.       &
 812 !                       *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
 813 !                       *work2(i,k)*coeres)
 814             psmlt(i,k) = (1.414e3*(1.496e-6 * ((t(i,k))*sqrt(t(i,k)))      &
 815                          /((t(i,k))+120.)/(den(i,k)))*(den(i,k)))/xlf      &
 816                         *(t0c-t(i,k))*pi/2.*n0sfac(i,k)                    &
 817                         *(precs1*rslope2(i,k,2)+precs2*work2(i,k)*coeres)  &
 818                         /den(i,k)
 819             psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),-qrs(i,k,2)     &
 820                           /mstep(i)),0.)
 821 !-------------------------------------------------------------------
 822 ! nsmlt: melgin of snow  [LH A27]
 823 !       (T>T0: ->NR)
 824 !-------------------------------------------------------------------
 825             if(qrs(i,k,2).gt.qcrmin) then
 826               sfac = rslope(i,k,2)*n0s*n0sfac(i,k)*mstep(i)/qrs(i,k,2)
 827               ncr(i,k,3) = ncr(i,k,3) - sfac*psmlt(i,k)
 828             endif
 829             qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
 830             qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
 831             t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
 832           endif
 833         enddo
 834       enddo
 835 !---------------------------------------------------------------
 836 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
 837 !---------------------------------------------------------------
 838       do k = kte, kts, -1
 839         do i = its, ite
 840           if(qci(i,k,2).le.0.) then
 841             work1c(i,k) = 0.
 842           else
 843             xmi = den(i,k)*qci(i,k,2)/xni(i,k)
 844             diameter  = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
 845             work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
 846           endif
 847         enddo
 848       enddo
 849 !
 850 !  forward semi-laglangian scheme (JH), PCM (piecewise constant),  (linear)
 851 !
 852       do k = kte, kts, -1
 853         do i = its, ite
 854           denqci(i,k) = den(i,k)*qci(i,k,2)
 855         enddo
 856       enddo
 857       call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci,  &
 858                            delqi,dtcld,1,0)
 859       do k = kts, kte
 860         do i = its, ite
 861           qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
 862         enddo
 863       enddo
 864       do i = its, ite
 865         fallc(i,1) = delqi(i)/delz(i,1)/dtcld
 866       enddo
 867 !
 868 !----------------------------------------------------------------
 869 !      rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
 870 !
 871       do i = its, ite
 872         fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
 873         fallsum_qsi = fall(i,1,2)+fallc(i,1)
 874         if(fallsum.gt.0.) then
 875           rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rainncv(i)
 876           rain(i) = fallsum*delz(i,1)/denr*dtcld*1000.+rain(i)
 877         endif
 878        if(fallsum_qsi.gt.0.) then
 879           tstepsnow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + tstepsnow(i)
 880           if (PRESENT (snowncv) .and. PRESENT (snow)) then
 881             snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snowncv(i)
 882             snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.+snow(i)
 883           endif
 884         endif
 885         IF ( PRESENT (snowncv) ) THEN
 886           if(fallsum.gt.0.) sr(i) = snowncv(i)/(rainncv(i)+1.e-12)
 887         ELSE
 888           if(fallsum.gt.0.)sr(i)= tstepsnow(i)/(rainncv(i)+1.e-12)
 889         ENDIF
 890       enddo
 891 !
 892 !---------------------------------------------------------------
 893 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
 894 !       (T>T0: QI->QC)
 895 !---------------------------------------------------------------
 896       do k = kts, kte
 897         do i = its, ite
 898           supcol = t0c-t(i,k)
 899           xlf = xls-xl(i,k)
 900           if(supcol.lt.0.) xlf = xlf0
 901           if(supcol.lt.0 .and. qci(i,k,2).gt.0.) then
 902             qci(i,k,1) = qci(i,k,1)+qci(i,k,2)
 903 !---------------------------------------------------------------
 904 ! nimlt: instantaneous melting of cloud ice  [LH A18]
 905 !        (T>T0: ->NC)
 906 !--------------------------------------------------------------
 907             ncr(i,k,2) = ncr(i,k,2) + xni(i,k)
 908             t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
 909             qci(i,k,2) = 0.
 910           endif
 911 !---------------------------------------------------------------
 912 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
 913 !        (T<-40C: QC->QI)
 914 !---------------------------------------------------------------
 915           if(supcol.gt.40. .and. qci(i,k,1).gt.0.) then
 916             qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
 917 !---------------------------------------------------------------
 918 ! nihmf: homogeneous  of cloud water below -40c [LH A17]
 919 !        (T<-40C: NC->)
 920 !---------------------------------------------------------------
 921             if(ncr(i,k,2).gt.0.) ncr(i,k,2) = 0.
 922             t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
 923             qci(i,k,1) = 0.
 924           endif
 925 !---------------------------------------------------------------
 926 ! pihtf: heterogeneous freezing of cloud water [HL A44]
 927 !        (T0>T>-40C: QC->QI)
 928 !---------------------------------------------------------------
 929           if(supcol.gt.0. .and. qci(i,k,1).gt.0.) then
 930             supcolt=min(supcol,70.)
 931             pfrzdtc = min(pi*pi*pfrz1*(exp(pfrz2*supcolt)-1.)*denr/den(i,k)    &
 932                    *ncr(i,k,2)*rslopec3(i,k)*rslopec3(i,k)/18.*dtcld,qci(i,k,1))
 933 !---------------------------------------------------------------
 934 ! nihtf: heterogeneous  of cloud water  [LH A16]
 935 !         (T0>T>-40C: NC->)
 936 !---------------------------------------------------------------
 937             if(ncr(i,k,2).gt.ncmin) then
 938               nfrzdtc = min(pi*pfrz1*(exp(pfrz2*supcolt)-1.)*ncr(i,k,2)        &
 939                       *rslopec3(i,k)/6.*dtcld,ncr(i,k,2))
 940               ncr(i,k,2) = ncr(i,k,2) - nfrzdtc
 941             endif
 942             qci(i,k,2) = qci(i,k,2) + pfrzdtc
 943             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
 944             qci(i,k,1) = qci(i,k,1)-pfrzdtc
 945            endif
 946 !---------------------------------------------------------------
 947 ! psfrz: freezing of rain water [HL A20] [LFO 45]
 948 !        (T<T0, QR->QS)
 949 !---------------------------------------------------------------
 950           if(supcol.gt.0. .and. qrs(i,k,1).gt.0.) then
 951             supcolt=min(supcol,70.)
 952             pfrzdtr = min(140.*(pi*pi)*pfrz1*ncr(i,k,3)*denr/den(i,k)          &
 953                   *(exp(pfrz2*supcolt)-1.)*rslope3(i,k,1)*rslope3(i,k,1)       &
 954                   *dtcld,qrs(i,k,1))
 955 !---------------------------------------------------------------
 956 ! nsfrz: freezing of rain water [LH A26]
 957 !        (T<T0, NR-> )
 958 !---------------------------------------------------------------
 959             if(ncr(i,k,3).gt.nrmin) then
 960               nfrzdtr = min(4.*pi*pfrz1*ncr(i,k,3)*(exp(pfrz2*supcolt)-1.)     &
 961                        *rslope3(i,k,1)*dtcld,ncr(i,k,3))
 962               ncr(i,k,3) = ncr(i,k,3)-nfrzdtr
 963             endif
 964             qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
 965             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
 966             qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
 967           endif
 968         enddo
 969       enddo
 970 !
 971       do k = kts, kte
 972         do i = its, ite
 973           ncr(i,k,2) = max(ncr(i,k,2),0.0)
 974           ncr(i,k,3) = max(ncr(i,k,3),0.0)
 975         enddo
 976       enddo
 977 !----------------------------------------------------------------
 978 !     update the slope parameters for microphysics computation
 979 !
 980       do k = kts, kte
 981         do i = its, ite
 982           qrs_tmp(i,k,1) = qrs(i,k,1)
 983           qrs_tmp(i,k,2) = qrs(i,k,2)
 984           ncr_tmp(i,k) = ncr(i,k,3)
 985         enddo
 986       enddo
 987       call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
 988                      rslope3,work1,workn,its,ite,kts,kte)
 989       do k = kts, kte
 990         do i = its, ite
 991 !-----------------------------------------------------------------
 992 ! compute the mean-volume drop diameter                  [LH A10]
 993 ! for raindrop distribution
 994 !-----------------------------------------------------------------
 995           avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
 996           if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin) then
 997             rslopec(i,k) = rslopecmax
 998             rslopec2(i,k) = rslopec2max
 999             rslopec3(i,k) = rslopec3max
1000           else
1001             rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
1002             rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
1003             rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
1004           endif
1005 !--------------------------------------------------------------------
1006 ! compute the mean-volume drop diameter                   [LH A7]
1007 ! for cloud-droplet distribution
1008 !--------------------------------------------------------------------
1009           avedia(i,k,1) = rslopec(i,k)
1010         enddo
1011       enddo
1012 !----------------------------------------------------------------
1013 !     work1:  the thermodynamic term in the denominator associated with
1014 !             heat conduction and vapor diffusion
1015 !             (ry88, y93, h85)
1016 !     work2: parameter associated with the ventilation effects(y93)
1017 !
1018       do k = kts, kte
1019         do i = its, ite
1020 !         work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
1021           work1(i,k,1) = ((((den(i,k))*(xl(i,k))*(xl(i,k)))*((t(i,k))+120.)    &
1022                        *(den(i,k)))/(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))&
1023                        *(den(i,k))*(rv*(t(i,k))*(t(i,k)))))                    &
1024                       +  p(i,k)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k))*(1.81))))
1025 !         work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
1026           work1(i,k,2) = ((((den(i,k))*(xls)*(xls))*((t(i,k))+120.)*(den(i,k)))&
1027                        /(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))*(den(i,k)) &
1028                        *(rv*(t(i,k))*(t(i,k))))                                &
1029                       + p(i,k)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k))*(1.81)))))
1030 !         work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
1031           work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) &
1032                      *p(i,k))/(((t(i,k))+120.)*den(i,k)*(8.794e-5              &
1033                      *exp(log(t(i,k))*(1.81))))))*sqrt(sqrt(den0/(den(i,k))))) &
1034                       /sqrt((1.496e-6*((t(i,k))*sqrt(t(i,k))))                 &
1035                       /(((t(i,k))+120.)*den(i,k)))
1036         enddo
1037       enddo
1038 !
1039 !===============================================================
1040 !
1041 ! warm rain processes
1042 !
1043 ! - follows the processes in RH83 and LFO except for autoconcersion
1044 !
1045 !===============================================================
1046 !
1047       do k = kts, kte
1048         do i = its, ite
1049           supsat = max(q(i,k),qmin)-qs(i,k,1)
1050           satdt = supsat/dtcld
1051 !---------------------------------------------------------------
1052 ! praut: auto conversion rate from cloud to rain [LH 9] [CP 17]
1053 !        (QC->QR)
1054 !---------------------------------------------------------------
1055           lencon  = 2.7e-2*den(i,k)*qci(i,k,1)*(1.e20/16.*rslopec2(i,k)        &
1056                    *rslopec2(i,k)-0.4)
1057           lenconcr = max(1.2*lencon,qcrmin)
1058           if(avedia(i,k,1).gt.di15) then
1059             taucon = 3.7/den(i,k)/qci(i,k,1)/(0.5e6*rslopec(i,k)-7.5)
1060             taucon = max(taucon, 1.e-12)
1061             praut(i,k) = lencon/(taucon*den(i,k))
1062             praut(i,k) = min(max(praut(i,k),0.),qci(i,k,1)/dtcld)
1063 !---------------------------------------------------------------
1064 ! nraut: auto conversion rate from cloud to rain [LH A6][CP 18 & 19]
1065 !        (NC->NR)
1066 !---------------------------------------------------------------
1067             nraut(i,k) = 3.5e9*den(i,k)*praut(i,k)
1068             if(qrs(i,k,1).gt.lenconcr)                                         &
1069             nraut(i,k) = ncr(i,k,3)/qrs(i,k,1)*praut(i,k)
1070             nraut(i,k) = min(nraut(i,k),ncr(i,k,2)/dtcld)
1071           endif
1072 !---------------------------------------------------------------
1073 ! pracw: accretion of cloud water by rain   [LH 10][CP 22 & 23]
1074 !        (QC->QR)
1075 ! nracw: accretion of cloud water by rain   [LH A9]
1076 !        (NC->)
1077 !---------------------------------------------------------------
1078           if(qrs(i,k,1).ge.lenconcr) then
1079             if(avedia(i,k,2).ge.di100) then
1080               nracw(i,k) = min(ncrk1*ncr(i,k,2)*ncr(i,k,3)*(rslopec3(i,k)      &
1081                          + 24.*rslope3(i,k,1)),ncr(i,k,2)/dtcld)
1082               pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk1*ncr(i,k,2)          &
1083                          *ncr(i,k,3)*rslopec3(i,k)*(2.*rslopec3(i,k)           &
1084                          + 24.*rslope3(i,k,1)),qci(i,k,1)/dtcld)
1085             else
1086               nracw(i,k) = min(ncrk2*ncr(i,k,2)*ncr(i,k,3)*(2.*rslopec3(i,k)   &
1087                          *rslopec3(i,k)+5040.*rslope3(i,k,1)                   &
1088                          *rslope3(i,k,1)),ncr(i,k,2)/dtcld)
1089               pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk2*ncr(i,k,2)          &
1090                          *ncr(i,k,3)*rslopec3(i,k)*(6.*rslopec3(i,k)           &
1091                          *rslopec3(i,k)+5040.*rslope3(i,k,1)                   &
1092                          *rslope3(i,k,1)),qci(i,k,1)/dtcld)
1093             endif
1094           endif
1095 !----------------------------------------------------------------
1096 ! nccol: self collection of cloud water         [LH A8][CP 24 & 25]
1097 !        (NC->)
1098 !----------------------------------------------------------------
1099           if(avedia(i,k,1).ge.di100) then
1100             nccol(i,k) = ncrk1*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)
1101           else
1102             nccol(i,k) = 2.*ncrk2*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)        &
1103                         *rslopec3(i,k)
1104           endif
1105 !----------------------------------------------------------------
1106 ! nrcol: self collection of rain-drops and break-up [LH A21][CP 24 & 25]
1107 !        (NR->)
1108 !----------------------------------------------------------------
1109           if(qrs(i,k,1).ge.lenconcr) then
1110             if(avedia(i,k,2).lt.di100) then
1111               nrcol(i,k) = 5040.*ncrk2*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1)    &
1112                           *rslope3(i,k,1)
1113             elseif(avedia(i,k,2).ge.di100 .and. avedia(i,k,2).lt.di600) then
1114               nrcol(i,k) = 24.*ncrk1*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1)
1115             elseif(avedia(i,k,2).ge.di600 .and. avedia(i,k,2).lt.di2000) then
1116               coecol = -2.5e3*(avedia(i,k,2)-di600)
1117               nrcol(i,k) = 24.*exp(coecol)*ncrk1*ncr(i,k,3)*ncr(i,k,3)         &
1118                          *rslope3(i,k,1)
1119             else
1120               nrcol(i,k) = 0.
1121             endif
1122           endif
1123 !---------------------------------------------------------------
1124 ! prevp: evaporation/condensation rate of rain  [HL A41]
1125 !        (QV->QR or QR->QV)
1126 !---------------------------------------------------------------
1127           if(qrs(i,k,1).gt.0.) then
1128             coeres = rslope(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
1129             prevp(i,k) = (rh(i,k,1)-1.)*ncr(i,k,3)*(precr1*rslope(i,k,1)       &
1130                          +precr2*work2(i,k)*coeres)/work1(i,k,1)
1131             if(prevp(i,k).lt.0.) then
1132               prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
1133               prevp(i,k) = max(prevp(i,k),satdt/2)
1134 !----------------------------------------------------------------
1135 ! Nrevp: evaporation/condensation rate of rain  [LH A14]
1136 !        (NR->NCCN)
1137 !----------------------------------------------------------------
1138               if(prevp(i,k).eq.-qrs(i,k,1)/dtcld) then
1139                  ncr(i,k,1) = ncr(i,k,1) + ncr(i,k,3)
1140                  ncr(i,k,3) = 0.
1141               endif
1142             else
1143 !
1144               prevp(i,k) = min(prevp(i,k),satdt/2)
1145             endif
1146           endif
1147         enddo
1148       enddo
1149 !
1150 !===============================================================
1151 !
1152 ! cold rain processes
1153 !
1154 ! - follows the revised ice microphysics processes in HDC
1155 ! - the processes same as in RH83 and RH84  and LFO behave
1156 !   following ice crystal hapits defined in HDC, inclduing
1157 !   intercept parameter for snow (n0s), ice crystal number
1158 !   concentration (ni), ice nuclei number concentration
1159 !   (n0i), ice diameter (d)
1160 !
1161 !===============================================================
1162 !
1163       rdtcld = 1./dtcld
1164       do k = kts, kte
1165         do i = its, ite
1166           supcol = t0c-t(i,k)
1167           n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1168           supsat = max(q(i,k),qmin)-qs(i,k,2)
1169           satdt = supsat/dtcld
1170           ifsat = 0
1171 !-------------------------------------------------------------
1172 ! Ni: ice crystal number concentraiton   [HDC 5c]
1173 !-------------------------------------------------------------
1174 !         xni(i,k) = min(max(5.38e7*(den(i,k)                                  &
1175 !                      *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
1176           temp = (den(i,k)*max(qci(i,k,2),qmin))
1177           temp = sqrt(sqrt(temp*temp*temp))
1178           xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1179           eacrs = exp(0.07*(-supcol))
1180 !
1181           if(supcol.gt.0) then
1182             if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,2).gt.qmin) then
1183               xmi = den(i,k)*qci(i,k,2)/xni(i,k)
1184               diameter  = min(dicon * sqrt(xmi),dimax)
1185               vt2i = 1.49e4*diameter**1.31
1186               vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
1187 !-------------------------------------------------------------
1188 ! psaci: Accretion of cloud ice by rain [HDC 10]
1189 !        (T<T0: QI->QS)
1190 !-------------------------------------------------------------
1191               acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2)            &
1192                       + diameter**2*rslope(i,k,2)
1193               psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k)*abs(vt2s-vt2i)  &
1194                          *acrfac/4.
1195             endif
1196           endif
1197 !-------------------------------------------------------------
1198 ! psacw: Accretion of cloud water by snow  [HL A7] [LFO 24]
1199 !        (T<T0: QC->QS, and T>=T0: QC->QR)
1200 !-------------------------------------------------------------
1201           if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,1).gt.qmin) then
1202             psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2)   &
1203                         *qci(i,k,1)*denfac(i,k),qci(i,k,1)*rdtcld)
1204           endif
1205 !-------------------------------------------------------------
1206 ! nsacw: Accretion of cloud water by snow  [LH A12]
1207 !        (NC ->)
1208 !-------------------------------------------------------------
1209          if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
1210            nsacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2)    &
1211                        *ncr(i,k,2)*denfac(i,k),ncr(i,k,2)/dtcld)
1212          endif
1213          if(supcol.le.0) then
1214            xlf = xlf0
1215 !--------------------------------------------------------------
1216 ! nseml: Enhanced melting of snow by accretion of water  [LH A29]
1217 !        (T>=T0: ->NR)
1218 !--------------------------------------------------------------
1219               if  (qrs(i,k,2).gt.qcrmin) then
1220                 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
1221                 nseml(i,k) = -sfac*min(max(cliq*supcol*(psacw(i,k))/xlf        &
1222                             ,-qrs(i,k,2)/dtcld),0.)
1223               endif
1224          endif
1225          if(supcol.gt.0) then
1226 !-------------------------------------------------------------
1227 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
1228 !       (T<T0: QV->QI or QI->QV)
1229 !-------------------------------------------------------------
1230             if(qci(i,k,2).gt.0 .and. ifsat.ne.1) then
1231               xmi = den(i,k)*qci(i,k,2)/xni(i,k)
1232               diameter = dicon * sqrt(xmi)
1233               pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1234               supice = satdt-prevp(i,k)
1235               if(pidep(i,k).lt.0.) then
1236 !               pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1237 !               pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1238                 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
1239                 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
1240               else
1241 !               pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1242                 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
1243               endif
1244               if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1245             endif
1246 !-------------------------------------------------------------
1247 ! psdep: deposition/sublimation rate of snow [HDC 14]
1248 !        (QV->QS or QS->QV)
1249 !-------------------------------------------------------------
1250             if(qrs(i,k,2).gt.0. .and. ifsat.ne.1) then
1251               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1252               psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2)   &
1253                            + precs2*work2(i,k)*coeres)/work1(i,k,2)
1254               supice = satdt-prevp(i,k)-pidep(i,k)
1255               if(psdep(i,k).lt.0.) then
1256 !               psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1257 !               psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1258                 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
1259                 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
1260               else
1261 !             psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1262                 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
1263               endif
1264               if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) ifsat = 1
1265             endif
1266 !-------------------------------------------------------------
1267 ! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
1268 !       (T<T0: QV->QI)
1269 !-------------------------------------------------------------
1270             if(supsat.gt.0 .and. ifsat.ne.1) then
1271               supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1272               xni0 = 1.e3*exp(0.1*supcol)
1273               roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
1274               pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))*rdtcld)
1275               pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1276             endif
1277 !
1278 !-------------------------------------------------------------
1279 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1280 !       (T<T0: QI->QS)
1281 !-------------------------------------------------------------
1282             if(qci(i,k,2).gt.0.) then
1283               qimax = roqimax/den(i,k)
1284 !             psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1285               psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
1286             endif
1287           endif
1288 !-------------------------------------------------------------
1289 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1290 !       (T>T0: QS->QV)
1291 !-------------------------------------------------------------
1292           if(supcol.lt.0.) then
1293             if(qrs(i,k,2).gt.0. .and. rh(i,k,1).lt.1.)                         &
1294               psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
1295 !              psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1296               psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
1297           endif
1298         enddo
1299       enddo
1300 !
1301 !
1302 !----------------------------------------------------------------
1303 !     check mass conservation of generation terms and feedback to the
1304 !     large scale
1305 !
1306       do k = kts, kte
1307         do i = its, ite
1308           if(t(i,k).le.t0c) then
1309 !
1310 !     Q_cloud water
1311 !
1312             value = max(qmin,qci(i,k,1))
1313             source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1314             if (source.gt.value) then
1315               factor = value/source
1316               praut(i,k) = praut(i,k)*factor
1317               pracw(i,k) = pracw(i,k)*factor
1318               psacw(i,k) = psacw(i,k)*factor
1319             endif
1320 !
1321 !     Q_cloud ice
1322 !
1323             value = max(qmin,qci(i,k,2))
1324             source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
1325             if (source.gt.value) then
1326               factor = value/source
1327               psaut(i,k) = psaut(i,k)*factor
1328               psaci(i,k) = psaci(i,k)*factor
1329               pigen(i,k) = pigen(i,k)*factor
1330               pidep(i,k) = pidep(i,k)*factor
1331             endif
1332 !
1333 !     Q_rain
1334 !
1335 !
1336             value = max(qmin,qrs(i,k,1))
1337             source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
1338             if (source.gt.value) then
1339               factor = value/source
1340               praut(i,k) = praut(i,k)*factor
1341               pracw(i,k) = pracw(i,k)*factor
1342               prevp(i,k) = prevp(i,k)*factor
1343             endif
1344 !
1345 !    Q_snow
1346 !
1347             value = max(qmin,qrs(i,k,2))
1348             source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
1349             if (source.gt.value) then
1350               factor = value/source
1351               psdep(i,k) = psdep(i,k)*factor
1352               psaut(i,k) = psaut(i,k)*factor
1353               psaci(i,k) = psaci(i,k)*factor
1354               psacw(i,k) = psacw(i,k)*factor
1355             endif
1356 !
1357 !     N_cloud
1358 !
1359             value = max(ncmin,ncr(i,k,2))
1360             source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+nsacw(i,k))*dtcld
1361             if (source.gt.value) then
1362               factor = value/source
1363               nraut(i,k) = nraut(i,k)*factor
1364               nccol(i,k) = nccol(i,k)*factor
1365               nracw(i,k) = nracw(i,k)*factor
1366               nsacw(i,k) = nsacw(i,k)*factor
1367             endif
1368 !
1369 !     N_rain
1370 !
1371             value = max(nrmin,ncr(i,k,3))
1372             source = (-nraut(i,k)+nrcol(i,k))*dtcld
1373             if (source.gt.value) then
1374               factor = value/source
1375               nraut(i,k) = nraut(i,k)*factor
1376               nrcol(i,k) = nrcol(i,k)*factor
1377             endif
1378 !
1379             work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
1380 !     update
1381             q(i,k) = q(i,k)+work2(i,k)*dtcld
1382             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)+psacw(i,k)      &
1383                         )*dtcld,0.)
1384             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)+prevp(i,k)      &
1385                         )*dtcld,0.)
1386             qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k)-pigen(i,k)      &
1387                         -pidep(i,k))*dtcld,0.)
1388             qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+psaci(i,k)      &
1389                         +psacw(i,k))*dtcld,0.)
1390             ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k)     &
1391                          -nsacw(i,k))*dtcld,0.)
1392             ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k))                &
1393                         *dtcld,0.)
1394             xlf = xls-xl(i,k)
1395             xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k))                  &
1396                       -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
1397             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1398           else
1399 !
1400 !     Q_cloud water
1401 !
1402             value = max(qmin,qci(i,k,1))
1403             source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1404             if (source.gt.value) then
1405               factor = value/source
1406               praut(i,k) = praut(i,k)*factor
1407               pracw(i,k) = pracw(i,k)*factor
1408               psacw(i,k) = psacw(i,k)*factor
1409             endif
1410 !
1411 !     Q_rain
1412 !
1413             value = max(qmin,qrs(i,k,1))
1414             source = (-praut(i,k)-pracw(i,k)-prevp(i,k)                        &
1415                     -psacw(i,k))*dtcld
1416             if (source.gt.value) then
1417               factor = value/source
1418               praut(i,k) = praut(i,k)*factor
1419               pracw(i,k) = pracw(i,k)*factor
1420               prevp(i,k) = prevp(i,k)*factor
1421               psacw(i,k) = psacw(i,k)*factor
1422             endif
1423 !
1424 !     Q_snow
1425 !
1426             value = max(qcrmin,qrs(i,k,2))
1427             source=(-psevp(i,k))*dtcld
1428             if (source.gt.value) then
1429               factor = value/source
1430               psevp(i,k) = psevp(i,k)*factor
1431             endif
1432 !
1433 !     N_cloud
1434 !
1435             value = max(ncmin,ncr(i,k,2))
1436             source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+nsacw(i,k))*dtcld
1437             if (source.gt.value) then
1438               factor = value/source
1439               nraut(i,k) = nraut(i,k)*factor
1440               nccol(i,k) = nccol(i,k)*factor
1441               nracw(i,k) = nracw(i,k)*factor
1442               nsacw(i,k) = nsacw(i,k)*factor
1443             endif
1444 !
1445 !     N_rain
1446 !
1447             value = max(nrmin,ncr(i,k,3))
1448             source = (-nraut(i,k)-nseml(i,k)+nrcol(i,k))*dtcld
1449             if (source.gt.value) then
1450               factor = value/source
1451               nraut(i,k) = nraut(i,k)*factor
1452               nseml(i,k) = nseml(i,k)*factor
1453               nrcol(i,k) = nrcol(i,k)*factor
1454             endif
1455             work2(i,k)=-(prevp(i,k)+psevp(i,k))
1456 !     update
1457             q(i,k) = q(i,k)+work2(i,k)*dtcld
1458             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)+psacw(i,k)      &
1459                         )*dtcld,0.)
1460             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)+prevp(i,k)      &
1461                         +psacw(i,k))*dtcld,0.)
1462             qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1463             ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k)     &
1464                         -nsacw(i,k))*dtcld,0.)
1465             ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)+nseml(i,k)-nrcol(i,k)      &
1466                         )*dtcld,0.)
1467             xlf = xls-xl(i,k)
1468             xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1469             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1470           endif
1471         enddo
1472       enddo
1473 !
1474 ! Inline expansion for fpvs
1475 !         qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1476 !         qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1477       hsub = xls
1478       hvap = xlv0
1479       cvap = cpv
1480       ttp=t0c+0.01
1481       dldt=cvap-cliq
1482       xa=-dldt/rv
1483       xb=xa+hvap/(rv*ttp)
1484       dldti=cvap-cice
1485       xai=-dldti/rv
1486       xbi=xai+hsub/(rv*ttp)
1487       do k = kts, kte
1488       do i = its, ite
1489           tr=ttp/t(i,k)
1490           logtr = log(tr)
1491           qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1492           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
1493           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1494           qs(i,k,1) = max(qs(i,k,1),qmin)
1495           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
1496         enddo
1497       enddo
1498 !
1499       call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
1500                      rslope3,work1,workn,its,ite,kts,kte)
1501       do k = kts, kte
1502         do i = its, ite
1503 !-----------------------------------------------------------------
1504 ! re-compute the mean-volume drop diameter            [LH A10]
1505 ! for raindrop distribution
1506 !-----------------------------------------------------------------
1507           avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
1508 !----------------------------------------------------------------
1509 ! Nrevp_s: evaporation/condensation rate of rain   [LH A14]
1510 !        (NR->NC)
1511 !----------------------------------------------------------------
1512           if(avedia(i,k,2).le.di82) then
1513             ncr(i,k,2) = ncr(i,k,2)+ncr(i,k,3)
1514             ncr(i,k,3) = 0.
1515 !----------------------------------------------------------------
1516 ! Prevp_s: evaporation/condensation rate of rain [LH A15] [KK 23]
1517 !        (QR->QC)
1518 !----------------------------------------------------------------
1519             qci(i,k,1) = qci(i,k,1)+qrs(i,k,1)
1520             qrs(i,k,1) = 0.
1521           endif
1522         enddo
1523       enddo
1524 !
1525       do k = kts, kte
1526         do i = its, ite
1527 !-------------------------------------------------------------------
1528 ! rate of change of cloud drop concentration due to CCN activation
1529 ! pcact: QV -> QC                                 [LH 8]  [KK 14]
1530 ! ncact: NCCN -> NC                               [LH A2] [KK 12]
1531 !-------------------------------------------------------------------
1532           if(rh(i,k,1).gt.1.) then
1533             ncact(i,k) = max(0.,((ncr(i,k,1)+ncr(i,k,2))                       &
1534                        *min(1.,(rh(i,k,1)/satmax)**actk) - ncr(i,k,2)))/dtcld
1535             ncact(i,k) =min(ncact(i,k),max(ncr(i,k,1),0.)/dtcld)
1536             pcact(i,k) = min(4.*pi*denr*(actr*1.E-6)**3*ncact(i,k)/            &
1537                          (3.*den(i,k)),max(q(i,k),0.)/dtcld)
1538             q(i,k) = max(q(i,k)-pcact(i,k)*dtcld,0.)
1539             qci(i,k,1) = max(qci(i,k,1)+pcact(i,k)*dtcld,0.)
1540             ncr(i,k,1) = max(ncr(i,k,1)-ncact(i,k)*dtcld,0.)
1541             ncr(i,k,2) = max(ncr(i,k,2)+ncact(i,k)*dtcld,0.)
1542             t(i,k) = t(i,k)+pcact(i,k)*xl(i,k)/cpm(i,k)*dtcld
1543           endif
1544 !---------------------------------------------------------------------
1545 !  pcond:condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1546 !     if there exists additional water vapor condensated/if
1547 !     evaporation of cloud water is not enough to remove subsaturation
1548 !  (QV->QC or QC->QV)
1549 !---------------------------------------------------------------------
1550           tr=ttp/t(i,k)
1551           qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1552           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
1553           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1554           qs(i,k,1) = max(qs(i,k,1),qmin)
1555           work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/(1.+(xl(i,k))         &
1556                        *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))/((t(i,k))        &
1557                        *(t(i,k)))))
1558           work2(i,k) = qci(i,k,1)+work1(i,k,1)
1559           pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1560           if(qci(i,k,1).gt.0. .and. work1(i,k,1).lt.0.)                        &
1561             pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1562 !----------------------------------------------------------------------
1563 ! ncevp: evpration of Cloud number concentration  [LH A3]
1564 ! (NC->NCCN)
1565 !----------------------------------------------------------------------
1566           if(pcond(i,k).eq.-qci(i,k,1)/dtcld) then
1567             ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,2)
1568             ncr(i,k,2) = 0.
1569           endif
1570 !
1571           q(i,k) = q(i,k)-pcond(i,k)*dtcld
1572           qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1573           t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1574         enddo
1575       enddo
1576 !
1577 !----------------------------------------------------------------
1578 !     padding for small values
1579 !
1580       do k = kts, kte
1581         do i = its, ite
1582           if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1583           if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1584           if(qrs(i,k,1).ge.qcrmin .and. ncr(i,k,3) .ge. nrmin) then
1585             lamdr_tmp(i,k) = exp(log(((pidnr*ncr(i,k,3))                       &
1586                          /(den(i,k)*qrs(i,k,1))))*((.33333333)))
1587             if(lamdr_tmp(i,k) .le. lamdarmin) then
1588               lamdr_tmp(i,k) = lamdarmin
1589               ncr(i,k,3) = den(i,k)*qrs(i,k,1)*lamdr_tmp(i,k)**3/pidnr
1590             elseif(lamdr_tmp(i,k) .ge. lamdarmax) then
1591               lamdr_tmp(i,k) = lamdarmax
1592               ncr(i,k,3) = den(i,k)*qrs(i,k,1)*lamdr_tmp(i,k)**3/pidnr
1593             endif
1594           endif
1595           if(qci(i,k,1).ge.qmin .and. ncr(i,k,2) .ge. ncmin ) then
1596             lamdc_tmp(i,k) = exp(log(((pidnc*ncr(i,k,2))                       &
1597                          /(den(i,k)*qci(i,k,1))))*((.33333333)))
1598             if(lamdc_tmp(i,k) .le. lamdacmin) then
1599               lamdc_tmp(i,k) = lamdacmin
1600               ncr(i,k,2) = den(i,k)*qci(i,k,1)*lamdc_tmp(i,k)**3/pidnc
1601             elseif(lamdc_tmp(i,k) .ge. lamdacmax) then
1602               lamdc_tmp(i,k) = lamdacmax
1603               ncr(i,k,2) = den(i,k)*qci(i,k,1)*lamdc_tmp(i,k)**3/pidnc
1604             endif
1605           endif
1606         enddo
1607       enddo
1608       enddo                  ! big loops
1609   END SUBROUTINE wdm52d
1610 ! ...................................................................
1611       REAL FUNCTION rgmma(x)
1612 !-------------------------------------------------------------------
1613   IMPLICIT NONE
1614 !-------------------------------------------------------------------
1615 !     rgmma function:  use infinite product form
1616       REAL :: euler
1617       PARAMETER (euler=0.577215664901532)
1618       REAL :: x, y
1619       INTEGER :: i
1620       if(x.eq.1.)then
1621         rgmma=0.
1622           else
1623         rgmma=x*exp(euler*x)
1624         do i=1,10000
1625           y=float(i)
1626           rgmma=rgmma*(1.000+x/y)*exp(-x/y)
1627         enddo
1628         rgmma=1./rgmma
1629       endif
1630       END FUNCTION rgmma
1631 !
1632 !--------------------------------------------------------------------------
1633       REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
1634 !--------------------------------------------------------------------------
1635       IMPLICIT NONE
1636 !--------------------------------------------------------------------------
1637       REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti,         &
1638            xai,xbi,ttp,tr
1639       INTEGER ice
1640 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1641       ttp=t0c+0.01
1642       dldt=cvap-cliq
1643       xa=-dldt/rv
1644       xb=xa+hvap/(rv*ttp)
1645       dldti=cvap-cice
1646       xai=-dldti/rv
1647       xbi=xai+hsub/(rv*ttp)
1648       tr=ttp/t
1649       if(t.lt.ttp .and. ice.eq.1) then
1650         fpvs=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
1651       else
1652         fpvs=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1653       endif
1654 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1655       END FUNCTION fpvs
1656 !-------------------------------------------------------------------
1657   SUBROUTINE wdm5init(den0,denr,dens,cl,cpv,ccn0,allowed_to_read)
1658 !-------------------------------------------------------------------
1659   IMPLICIT NONE
1660 !-------------------------------------------------------------------
1661 !.... constants which may not be tunable
1662    REAL, INTENT(IN) :: den0,denr,dens,cl,cpv,ccn0
1663    LOGICAL, INTENT(IN) :: allowed_to_read
1664 !
1665    pi = 4.*atan(1.)
1666    xlv1 = cl-cpv
1667 !
1668    qc0  = 4./3.*pi*denr*r0**3*xncr/den0  ! 0.419e-3 -- .61e-3
1669    qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
1670    pidnc = pi*denr/6.
1671 !
1672    bvtr1 = 1.+bvtr
1673    bvtr2 = 2.+bvtr
1674    bvtr3 = 3.+bvtr
1675    bvtr4 = 4.+bvtr
1676    bvtr5 = 5.+bvtr
1677    bvtr7 = 7.+bvtr
1678    bvtr2o5 = 2.5+.5*bvtr
1679    bvtr3o5 = 3.5+.5*bvtr
1680    g1pbr = rgmma(bvtr1)
1681    g2pbr = rgmma(bvtr2)
1682    g3pbr = rgmma(bvtr3)
1683    g4pbr = rgmma(bvtr4)            ! 17.837825
1684    g5pbr = rgmma(bvtr5)
1685    g7pbr = rgmma(bvtr7)
1686    g5pbro2 = rgmma(bvtr2o5)
1687    g7pbro2 = rgmma(bvtr3o5)
1688    pvtr = avtr*g5pbr/24.
1689    pvtrn = avtr*g2pbr
1690    eacrr = 1.0
1691    pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
1692    precr1 = 2.*pi*1.56
1693    precr2 = 2.*pi*.31*avtr**.5*g7pbro2
1694    pidn0r =  pi*denr*n0r
1695    pidnr = 4.*pi*denr
1696    xmmax = (dimax/dicon)**2
1697    roqimax = 2.08e22*dimax**8
1698 !
1699    bvts1 = 1.+bvts
1700    bvts2 = 2.5+.5*bvts
1701    bvts3 = 3.+bvts
1702    bvts4 = 4.+bvts
1703    g1pbs = rgmma(bvts1)    !.8875
1704    g3pbs = rgmma(bvts3)
1705    g4pbs = rgmma(bvts4)    ! 12.0786
1706    g5pbso2 = rgmma(bvts2)
1707    pvts = avts*g4pbs/6.
1708    pacrs = pi*n0s*avts*g3pbs*.25
1709    precs1 = 4.*n0s*.65
1710    precs2 = 4.*n0s*.44*avts**.5*g5pbso2
1711    pidn0s =  pi*dens*n0s
1712    pacrc = pi*n0s*avts*g3pbs*.25*eacrc
1713 !
1714    rslopecmax = 1./lamdacmax
1715    rslopermax = 1./lamdarmax
1716    rslopesmax = 1./lamdasmax
1717    rsloperbmax = rslopermax ** bvtr
1718    rslopesbmax = rslopesmax ** bvts
1719    rslopec2max = rslopecmax * rslopecmax
1720    rsloper2max = rslopermax * rslopermax
1721    rslopes2max = rslopesmax * rslopesmax
1722    rslopec3max = rslopec2max * rslopecmax
1723    rsloper3max = rsloper2max * rslopermax
1724    rslopes3max = rslopes2max * rslopesmax
1725 !
1726 !+---+-----------------------------------------------------------------+
1727 !..Set these variables needed for computing radar reflectivity.  These
1728 !.. get used within radar_init to create other variables used in the
1729 !.. radar module.
1730    xam_r = PI*denr/6.
1731    xbm_r = 3.
1732    xmu_r = 0.
1733    xam_s = PI*dens/6.
1734    xbm_s = 3.
1735    xmu_s = 0.
1736    xam_g = PI*dens/6.      !  These 3 variables for graupel are set but unused.
1737    xbm_g = 3.
1738    xmu_g = 0.
1739
1740    call radar_init
1741 !+---+-----------------------------------------------------------------+
1742
1743   END SUBROUTINE wdm5init
1744 !------------------------------------------------------------------------------
1745       subroutine slope_wdm5(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1746                             vt,vtn,its,ite,kts,kte)
1747   IMPLICIT NONE
1748   INTEGER       ::               its,ite, jts,jte, kts,kte
1749   REAL, DIMENSION( its:ite , kts:kte,2) ::                                     &
1750                                                                           qrs, &
1751                                                                        rslope, &
1752                                                                       rslopeb, &
1753                                                                       rslope2, &
1754                                                                       rslope3, &
1755                                                                            vt
1756   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
1757                                                                           ncr, &
1758                                                                           vtn, &
1759                                                                           den, &
1760                                                                        denfac, &
1761                                                                             t
1762   REAL, PARAMETER  :: t0c = 273.15
1763   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
1764                                                                        n0sfac
1765   REAL       ::  lamdar, lamdas, x, y, z, supcol
1766   integer :: i, j, k
1767 !----------------------------------------------------------------
1768 !     size distributions: (x=mixing ratio, y=air density):
1769 !     valid for mixing ratio > 1.e-9 kg/kg.
1770 !
1771 !  Optimizatin : A**B => exp(log(A)*(B))
1772       lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
1773       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
1774 !
1775       do k = kts, kte
1776         do i = its, ite
1777           supcol = t0c-t(i,k)
1778 !---------------------------------------------------------------
1779 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1780 !---------------------------------------------------------------
1781           n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1782           if(qrs(i,k,1).le.qcrmin .or. ncr(i,k).le.nrmin ) then
1783             rslope(i,k,1) = rslopermax
1784             rslopeb(i,k,1) = rsloperbmax
1785             rslope2(i,k,1) = rsloper2max
1786             rslope3(i,k,1) = rsloper3max
1787           else
1788             rslope(i,k,1) = min(1./lamdar(qrs(i,k,1),den(i,k),ncr(i,k)),1.e-3)
1789             rslopeb(i,k,1) = rslope(i,k,1)**bvtr
1790             rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1791             rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1792           endif
1793           if(qrs(i,k,2).le.qcrmin) then
1794             rslope(i,k,2) = rslopesmax
1795             rslopeb(i,k,2) = rslopesbmax
1796             rslope2(i,k,2) = rslopes2max
1797             rslope3(i,k,2) = rslopes3max
1798           else
1799             rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
1800             rslopeb(i,k,2) = rslope(i,k,2)**bvts
1801             rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1802             rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1803           endif
1804           vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
1805           vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
1806           vtn(i,k) = pvtrn*rslopeb(i,k,1)*denfac(i,k)
1807           if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
1808           if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
1809           if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
1810         enddo
1811       enddo
1812   END subroutine slope_wdm5
1813 !-----------------------------------------------------------------------------
1814       subroutine slope_rain(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
1815                             vt,vtn,its,ite,kts,kte)
1816   IMPLICIT NONE
1817   INTEGER       ::               its,ite, jts,jte, kts,kte
1818   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
1819                                                                           qrs, &
1820                                                                           ncr, &
1821                                                                        rslope, &
1822                                                                       rslopeb, &
1823                                                                       rslope2, &
1824                                                                       rslope3, &
1825                                                                            vt, &
1826                                                                           vtn, &
1827                                                                           den, &
1828                                                                        denfac, &
1829                                                                             t
1830   REAL, PARAMETER  :: t0c = 273.15
1831   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
1832                                                                        n0sfac
1833   REAL       ::  lamdar, x, y, z, supcol
1834   integer :: i, j, k
1835 !----------------------------------------------------------------
1836 !     size distributions: (x=mixing ratio, y=air density):
1837 !     valid for mixing ratio > 1.e-9 kg/kg.
1838       lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
1839 !
1840       do k = kts, kte
1841         do i = its, ite
1842           if(qrs(i,k).le.qcrmin .or. ncr(i,k).le.nrmin) then
1843             rslope(i,k) = rslopermax
1844             rslopeb(i,k) = rsloperbmax
1845             rslope2(i,k) = rsloper2max
1846             rslope3(i,k) = rsloper3max
1847           else
1848             rslope(i,k) = min(1./lamdar(qrs(i,k),den(i,k),ncr(i,k)),1.e-3)
1849             rslopeb(i,k) = rslope(i,k)**bvtr
1850             rslope2(i,k) = rslope(i,k)*rslope(i,k)
1851             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1852           endif
1853           vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
1854           vtn(i,k) = pvtrn*rslopeb(i,k)*denfac(i,k)
1855           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1856           if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
1857         enddo
1858       enddo
1859   END subroutine slope_rain
1860 !------------------------------------------------------------------------------
1861       subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
1862                             vt,its,ite,kts,kte)
1863   IMPLICIT NONE
1864   INTEGER       ::               its,ite, jts,jte, kts,kte
1865   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
1866                                                                           qrs, &
1867                                                                        rslope, &
1868                                                                       rslopeb, &
1869                                                                       rslope2, &
1870                                                                       rslope3, &
1871                                                                            vt, &
1872                                                                           den, &
1873                                                                        denfac, &
1874                                                                             t
1875   REAL, PARAMETER  :: t0c = 273.15
1876   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
1877                                                                        n0sfac
1878   REAL       ::  lamdas, x, y, z, supcol
1879   integer :: i, j, k
1880 !----------------------------------------------------------------
1881 !     size distributions: (x=mixing ratio, y=air density):
1882 !     valid for mixing ratio > 1.e-9 kg/kg.
1883       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
1884 !
1885       do k = kts, kte
1886         do i = its, ite
1887           supcol = t0c-t(i,k)
1888 !---------------------------------------------------------------
1889 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1890 !---------------------------------------------------------------
1891           n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1892           if(qrs(i,k).le.qcrmin)then
1893             rslope(i,k) = rslopesmax
1894             rslopeb(i,k) = rslopesbmax
1895             rslope2(i,k) = rslopes2max
1896             rslope3(i,k) = rslopes3max
1897           else
1898             rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
1899             rslopeb(i,k) = rslope(i,k)**bvts
1900             rslope2(i,k) = rslope(i,k)*rslope(i,k)
1901             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1902           endif
1903           vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
1904           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1905         enddo
1906       enddo
1907   END subroutine slope_snow
1908 !-------------------------------------------------------------------
1909       SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
1910 !-------------------------------------------------------------------
1911 !
1912 ! for non-iteration semi-Lagrangain forward advection for cloud
1913 ! with mass conservation and positive definite advection
1914 ! 2nd order interpolation with monotonic piecewise linear method
1915 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
1916 !
1917 ! dzl    depth of model layer in meter
1918 ! wwl    terminal velocity at model layer m/s
1919 ! rql    cloud density*mixing ration
1920 ! precip precipitation
1921 ! dt     time step
1922 ! id     kind of precip: 0 test case; 1 raindrop  2: snow
1923 ! iter   how many time to guess mean terminal velocity: 0 pure forward.
1924 !        0 : use departure wind for advection
1925 !        1 : use mean wind for advection
1926 !        > 1 : use mean wind after iter-1 iterations
1927 !
1928 ! author: hann-ming henry juang <henry.juang@noaa.gov>
1929 !         implemented by song-you hong
1930 !
1931       implicit none
1932       integer  im,km,id
1933       real  dt
1934       real  dzl(im,km),wwl(im,km),rql(im,km),precip(im)
1935       real  denl(im,km),denfacl(im,km),tkl(im,km)
1936 !
1937       integer  i,k,n,m,kk,kb,kt,iter
1938       real  tl,tl2,qql,dql,qqd
1939       real  th,th2,qqh,dqh
1940       real  zsum,qsum,dim,dip,c1,con1,fa1,fa2
1941       real  allold, allnew, zz, dzamin, cflmax, decfl
1942       real  dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
1943       real  den(km), denfac(km), tk(km)
1944       real  wi(km+1), zi(km+1), za(km+1)
1945       real  qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
1946       real  dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
1947 !
1948       precip(:) = 0.0
1949 !
1950       i_loop : do i=1,im
1951 ! -----------------------------------
1952       dz(:) = dzl(i,:)
1953       qq(:) = rql(i,:)
1954       ww(:) = wwl(i,:)
1955       den(:) = denl(i,:)
1956       denfac(:) = denfacl(i,:)
1957       tk(:) = tkl(i,:)
1958 ! skip for no precipitation for all layers
1959       allold = 0.0
1960       do k=1,km
1961         allold = allold + qq(k)
1962       enddo
1963       if(allold.le.0.0) then
1964         cycle i_loop
1965       endif
1966 !
1967 ! compute interface values
1968       zi(1)=0.0
1969       do k=1,km
1970         zi(k+1) = zi(k)+dz(k)
1971       enddo
1972 !
1973 ! save departure wind
1974       wd(:) = ww(:)
1975       n=1
1976  100  continue
1977 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
1978 ! 2nd order interpolation to get wi
1979       wi(1) = ww(1)
1980       wi(km+1) = ww(km)
1981       do k=2,km
1982         wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
1983       enddo
1984 ! 3rd order interpolation to get wi
1985       fa1 = 9./16.
1986       fa2 = 1./16.
1987       wi(1) = ww(1)
1988       wi(2) = 0.5*(ww(2)+ww(1))
1989       do k=3,km-1
1990         wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
1991       enddo
1992       wi(km) = 0.5*(ww(km)+ww(km-1))
1993       wi(km+1) = ww(km)
1994 !
1995 ! terminate of top of raingroup
1996       do k=2,km
1997         if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1998       enddo
1999 !
2000 ! diffusivity of wi
2001       con1 = 0.05
2002       do k=km,1,-1
2003         decfl = (wi(k+1)-wi(k))*dt/dz(k)
2004         if( decfl .gt. con1 ) then
2005           wi(k) = wi(k+1) - con1*dz(k)/dt
2006         endif
2007       enddo
2008 ! compute arrival point
2009       do k=1,km+1
2010         za(k) = zi(k) - wi(k)*dt
2011       enddo
2012 !
2013       do k=1,km
2014         dza(k) = za(k+1)-za(k)
2015       enddo
2016       dza(km+1) = zi(km+1) - za(km+1)
2017 !
2018 ! computer deformation at arrival point
2019       do k=1,km
2020         qa(k) = qq(k)*dz(k)/dza(k)
2021         qr(k) = qa(k)/den(k)
2022       enddo
2023       qa(km+1) = 0.0
2024 !     call maxmin(km,1,qa,' arrival points ')
2025 !
2026 ! compute arrival terminal velocity, and estimate mean terminal velocity
2027 ! then back to use mean terminal velocity
2028       if( n.le.iter ) then
2029 !        if (id.eq.1) then
2030 !
2031 !          call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
2032 !        else
2033           call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
2034 !        endif
2035         if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
2036         do k=1,km
2037 !#ifdef DEBUG
2038 !        print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
2039 !#endif
2040 ! mean wind is average of departure and new arrival winds
2041           ww(k) = 0.5* ( wd(k)+wa(k) )
2042         enddo
2043         was(:) = wa(:)
2044         n=n+1
2045         go to 100
2046       endif
2047 !
2048 ! estimate values at arrival cell interface with monotone
2049       do k=2,km
2050         dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
2051         dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
2052         if( dip*dim.le.0.0 ) then
2053           qmi(k)=qa(k)
2054           qpi(k)=qa(k)
2055         else
2056           qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
2057           qmi(k)=2.0*qa(k)-qpi(k)
2058           if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
2059             qpi(k) = qa(k)
2060             qmi(k) = qa(k)
2061           endif
2062         endif
2063       enddo
2064       qpi(1)=qa(1)
2065       qmi(1)=qa(1)
2066       qmi(km+1)=qa(km+1)
2067       qpi(km+1)=qa(km+1)
2068 !
2069 ! interpolation to regular point
2070       qn = 0.0
2071       kb=1
2072       kt=1
2073       intp : do k=1,km
2074              kb=max(kb-1,1)
2075              kt=max(kt-1,1)
2076 ! find kb and kt
2077              if( zi(k).ge.za(km+1) ) then
2078                exit intp
2079              else
2080                find_kb : do kk=kb,km
2081                          if( zi(k).le.za(kk+1) ) then
2082                            kb = kk
2083                            exit find_kb
2084                          else
2085                            cycle find_kb
2086                          endif
2087                enddo find_kb
2088                find_kt : do kk=kt,km
2089                          if( zi(k+1).le.za(kk) ) then
2090                            kt = kk
2091                            exit find_kt
2092                          else
2093                            cycle find_kt
2094                          endif
2095                enddo find_kt
2096                kt = kt - 1
2097 ! compute q with piecewise constant method
2098                if( kt.eq.kb ) then
2099                  tl=(zi(k)-za(kb))/dza(kb)
2100                  th=(zi(k+1)-za(kb))/dza(kb)
2101                  tl2=tl*tl
2102                  th2=th*th
2103                  qqd=0.5*(qpi(kb)-qmi(kb))
2104                  qqh=qqd*th2+qmi(kb)*th
2105                  qql=qqd*tl2+qmi(kb)*tl
2106                  qn(k) = (qqh-qql)/(th-tl)
2107                else if( kt.gt.kb ) then
2108                  tl=(zi(k)-za(kb))/dza(kb)
2109                  tl2=tl*tl
2110                  qqd=0.5*(qpi(kb)-qmi(kb))
2111                  qql=qqd*tl2+qmi(kb)*tl
2112                  dql = qa(kb)-qql
2113                  zsum  = (1.-tl)*dza(kb)
2114                  qsum  = dql*dza(kb)
2115                  if( kt-kb.gt.1 ) then
2116                  do m=kb+1,kt-1
2117                    zsum = zsum + dza(m)
2118                    qsum = qsum + qa(m) * dza(m)
2119                  enddo
2120                  endif
2121                  th=(zi(k+1)-za(kt))/dza(kt)
2122                  th2=th*th
2123                  qqd=0.5*(qpi(kt)-qmi(kt))
2124                  dqh=qqd*th2+qmi(kt)*th
2125                  zsum  = zsum + th*dza(kt)
2126                  qsum  = qsum + dqh*dza(kt)
2127                  qn(k) = qsum/zsum
2128                endif
2129                cycle intp
2130              endif
2131 !
2132        enddo intp
2133 !
2134 ! rain out
2135       sum_precip: do k=1,km
2136                     if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
2137                       precip(i) = precip(i) + qa(k)*dza(k)
2138                       cycle sum_precip
2139                     else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
2140                       precip(i) = precip(i) + qa(k)*(0.0-za(k))
2141                       exit sum_precip
2142                     endif
2143                     exit sum_precip
2144       enddo sum_precip
2145 !
2146 ! replace the new values
2147       rql(i,:) = qn(:)
2148 !
2149 ! ----------------------------------
2150       enddo i_loop
2151 !
2152   END SUBROUTINE nislfv_rain_plm
2153
2154
2155 !+---+-----------------------------------------------------------------+
2156       subroutine refl10cm_wdm5 (qv1d, qr1d, nr1d, qs1d,                 &
2157                        t1d, p1d, dBZ, kts, kte, ii, jj)
2158
2159       IMPLICIT NONE
2160
2161 !..Sub arguments
2162       INTEGER, INTENT(IN):: kts, kte, ii, jj
2163       REAL, DIMENSION(kts:kte), INTENT(IN)::                            &
2164                       qv1d, qr1d, nr1d, qs1d, t1d, p1d
2165       REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
2166
2167 !..Local variables
2168       REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
2169       REAL, DIMENSION(kts:kte):: rr, nr, rs
2170       REAL:: temp_C
2171
2172       DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams
2173       DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s
2174       DOUBLE PRECISION:: lamr, lams
2175       LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs
2176
2177       REAL, DIMENSION(kts:kte):: ze_rain, ze_snow
2178       DOUBLE PRECISION:: fmelt_s
2179
2180       INTEGER:: i, k, k_0, kbot, n
2181       LOGICAL:: melti
2182
2183       DOUBLE PRECISION:: cback, x, eta, f_d
2184       REAL, PARAMETER:: R=287.
2185
2186 !+---+
2187
2188       do k = kts, kte
2189          dBZ(k) = -35.0
2190       enddo
2191
2192 !+---+-----------------------------------------------------------------+
2193 !..Put column of data into local arrays.
2194 !+---+-----------------------------------------------------------------+
2195       do k = kts, kte
2196          temp(k) = t1d(k)
2197          temp_C = min(-0.001, temp(K)-273.15)
2198          qv(k) = MAX(1.E-10, qv1d(k))
2199          pres(k) = p1d(k)
2200          rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
2201
2202          if (qr1d(k) .gt. 1.E-9) then
2203             rr(k) = qr1d(k)*rho(k)
2204             nr(k) = nr1d(k)*rho(k)
2205             lamr = (xam_r*xcrg(3)*xorg2*nr(k)/rr(k))**xobmr
2206             ilamr(k) = 1./lamr
2207             N0_r(k) = nr(k)*xorg2*lamr**xcre(2)
2208             L_qr(k) = .true.
2209          else
2210             rr(k) = 1.E-12
2211             nr(k) = 1.E-12
2212             L_qr(k) = .false.
2213          endif
2214
2215          if (qs1d(k) .gt. 1.E-9) then
2216             rs(k) = qs1d(k)*rho(k)
2217             N0_s(k) = min(n0smax, n0s*exp(-alpha*temp_C))
2218             lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1))
2219             ilams(k) = 1./lams
2220             L_qs(k) = .true.
2221          else
2222             rs(k) = 1.E-12
2223             L_qs(k) = .false.
2224          endif
2225
2226       enddo
2227
2228 !+---+-----------------------------------------------------------------+
2229 !..Locate K-level of start of melting (k_0 is level above).
2230 !+---+-----------------------------------------------------------------+
2231       melti = .false.
2232       k_0 = kts
2233       do k = kte-1, kts, -1
2234          if ( (temp(k).gt.273.15) .and. L_qr(k) .and. L_qs(k+1) ) then
2235             k_0 = MAX(k+1, k_0)
2236             melti=.true.
2237             goto 195
2238          endif
2239       enddo
2240  195  continue
2241
2242 !+---+-----------------------------------------------------------------+
2243 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
2244 !.. and non-water-coated snow and graupel when below freezing are
2245 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
2246 !+---+-----------------------------------------------------------------+
2247
2248       do k = kts, kte
2249          ze_rain(k) = 1.e-22
2250          ze_snow(k) = 1.e-22
2251          if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
2252          if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI)     &
2253                                  * (xam_s/900.0)*(xam_s/900.0)          &
2254                                  * N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
2255       enddo
2256
2257
2258 !+---+-----------------------------------------------------------------+
2259 !..Special case of melting ice (snow/graupel) particles.  Assume the
2260 !.. ice is surrounded by the liquid water.  Fraction of meltwater is
2261 !.. extremely simple based on amount found above the melting level.
2262 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
2263 !.. routines).
2264 !+---+-----------------------------------------------------------------+
2265
2266       if (melti .and. k_0.ge.kts+1) then
2267        do k = k_0-1, kts, -1
2268
2269 !..Reflectivity contributed by melting snow
2270           if (L_qs(k) .and. L_qs(k_0) ) then
2271            fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
2272            eta = 0.d0
2273            lams = 1./ilams(k)
2274            do n = 1, nrbins
2275               x = xam_s * xxDs(n)**xbm_s
2276               call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
2277                     fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
2278                     CBACK, mixingrulestring_s, matrixstring_s,          &
2279                     inclusionstring_s, hoststring_s,                    &
2280                     hostmatrixstring_s, hostinclusionstring_s)
2281               f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
2282               eta = eta + f_d * CBACK * simpson(n) * xdts(n)
2283            enddo
2284            ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
2285           endif
2286        enddo
2287       endif
2288
2289       do k = kte, kts, -1
2290          dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k))*1.d18)
2291       enddo
2292
2293
2294       end subroutine refl10cm_wdm5
2295 !+---+-----------------------------------------------------------------+
2296
2297 !-----------------------------------------------------------------------
2298      subroutine effectRad_wdm5 (t, qc, nc, qi, qs, rho, qmin, t0c,        &
2299                                 re_qc, re_qi, re_qs, kts, kte, ii, jj)
2300
2301 !-----------------------------------------------------------------------
2302 !  Compute radiation effective radii of cloud water, ice, and snow for
2303 !  double-moment microphysics..
2304 !  These are entirely consistent with microphysics assumptions, not
2305 !  constant or otherwise ad hoc as is internal to most radiation
2306 !  schemes.
2307 !  Coded and implemented by Soo Ya Bae, KIAPS, January 2015.
2308 !-----------------------------------------------------------------------
2309      implicit none
2310
2311 !..Sub arguments
2312       integer, intent(in) :: kts, kte, ii, jj
2313       real, intent(in) :: qmin
2314       real, intent(in) :: t0c
2315       real, dimension( kts:kte ), intent(in)::  t
2316       real, dimension( kts:kte ), intent(in)::  qc
2317       real, dimension( kts:kte ), intent(in)::  nc
2318       real, dimension( kts:kte ), intent(in)::  qi
2319       real, dimension( kts:kte ), intent(in)::  qs
2320       real, dimension( kts:kte ), intent(in)::  rho
2321       real, dimension( kts:kte ), intent(inout):: re_qc
2322       real, dimension( kts:kte ), intent(inout):: re_qi
2323       real, dimension( kts:kte ), intent(inout):: re_qs
2324 !..Local variables
2325       integer:: i,k
2326       integer :: inu_c
2327       real, dimension( kts:kte ):: ni
2328       real, dimension( kts:kte ):: rqc
2329       real, dimension( kts:kte ):: rnc
2330       real, dimension( kts:kte ):: rqi
2331       real, dimension( kts:kte ):: rni
2332       real, dimension( kts:kte ):: rqs
2333       real :: cdm2
2334       real :: temp
2335       real :: supcol, n0sfac, lamdas
2336       real :: diai      ! diameter of ice in m
2337       double precision :: lamc
2338       logical :: has_qc, has_qi, has_qs
2339 !..Minimum microphys values
2340       real, parameter :: R1 = 1.E-12
2341       real, parameter :: R2 = 1.E-6
2342       real, parameter :: pi = 3.1415926536
2343       real, parameter :: bm_r = 3.0
2344       real, parameter :: obmr = 1.0/bm_r
2345       real, parameter :: cdm  = 5./3.
2346
2347       has_qc = .false.
2348       has_qi = .false.
2349       has_qs = .false.
2350
2351       cdm2 = rgmma(cdm)
2352
2353       do k=kts,kte
2354         ! for cloud
2355         rqc(k) = max(R1, qc(k)*rho(k))
2356         rnc(k) = max(R2, nc(k)*rho(k))
2357         if (rqc(k).gt.R1 .and. rnc(k).gt.R2) has_qc = .true.
2358         ! for ice
2359         rqi(k) = max(R1, qi(k)*rho(k))
2360         temp = (rho(k)*max(qi(k),qmin))
2361         temp = sqrt(sqrt(temp*temp*temp))
2362         ni(k) = min(max(5.38e7*temp,1.e3),1.e6)
2363         rni(k)= max(R2, ni(k)*rho(k))
2364         if (rqi(k).gt.R1 .and. rni(k).gt.R2) has_qi = .true.
2365         ! for snow
2366         rqs(k) = max(R1, qs(k)*rho(k))
2367         if (rqs(k).gt.R1) has_qs = .true.
2368       enddo
2369
2370       if (has_qc) then
2371         do k=kts,kte
2372           if (rqc(k).le.R1 .or. rnc(k).le.R2) CYCLE
2373           lamc = 2.*cdm2*(pidnc*nc(k)/rqc(k))**obmr
2374           re_qc(k) =  max(2.51e-6,min(sngl(1.0d0/lamc),50.e-6))
2375         enddo
2376       endif
2377
2378       if (has_qi) then
2379         do k=kts,kte
2380           if (rqi(k).le.R1 .or. rni(k).le.R2) CYCLE
2381           diai = 11.9*sqrt(rqi(k)/ni(k))
2382           re_qi(k) = max(10.01e-6,min(0.75*0.163*diai,125.e-6))
2383         enddo
2384       endif
2385
2386       if (has_qs) then
2387         do k=kts,kte
2388           if (rqs(k).le.R1) CYCLE
2389           supcol = t0c-t(k)
2390           n0sfac = max(min(exp(alpha*supcol),n0smax/n0s),1.)
2391           lamdas = sqrt(sqrt(pidn0s*n0sfac/rqs(k)))
2392           re_qs(k) = max(25.e-6,min(0.5*(1./lamdas),999.e-6))
2393         enddo
2394       endif
2395
2396       end subroutine effectRad_wdm5
2397 !-----------------------------------------------------------------------
2398
2399 END MODULE module_mp_wdm5