phys/module_mp_wsm6r.F

   1 MODULE module_mp_wsm6r
   2 implicit none
   3 !
   4 ! parameters
   5    REAL, PARAMETER :: dtcldcr   = 120.
   6    REAL, PARAMETER :: n0r       = 8.e6
   7    REAL, PARAMETER :: n0g       = 4.e6
   8    REAL, PARAMETER :: avtr      = 841.9
   9    REAL, PARAMETER :: bvtr      = 0.8
  10    REAL, PARAMETER :: r0        = .8e-5     ! 8 microm  in contrast to 10 micro m
  11    REAL, PARAMETER :: peaut     = .55       ! collection efficiency
  12    REAL, PARAMETER :: xncr      = 3.e8      ! maritime cloud in contrast to 3.e8 in tc80
  13    REAL, PARAMETER :: xmyu      = 1.718e-5  ! the dynamic viscosity kgm-1s-1
  14    REAL, PARAMETER :: avts      = 11.72
  15    REAL, PARAMETER :: bvts      = .41
  16    REAL, PARAMETER :: avtg      = 330.
  17    REAL, PARAMETER :: bvtg      = 0.8
  18    REAL, PARAMETER :: deng      = 500.
  19    REAL, PARAMETER :: n0smax    = 1.e11     ! t                                                       = -90C unlimited
  20 !   REAL, PARAMETER :: betai     = .6
  21 !   REAL, PARAMETER :: xn0       = 1.e-2
  22 !   REAL, PARAMETER :: dicon     = 11.9
  23 !   REAL, PARAMETER :: di0       = 12.9e-6
  24    REAL, PARAMETER :: dimax     = 500.e-6
  25    REAL, PARAMETER :: n0s       = 2.e6      ! temperature dependent n0s
  26    REAL, PARAMETER :: alpha     = .12       ! .122 exponen factor for n0s
  27    REAL, PARAMETER :: pfrz1     = 100.
  28    REAL, PARAMETER :: pfrz2     = 0.66
  29    REAL, PARAMETER :: t40c      = 233.16
  30    REAL, PARAMETER :: eacrc     = 1.0       ! Esc
  31    REAL, PARAMETER :: eacrr     = 1.0       ! Erc
  32    REAL, PARAMETER :: dens      = 100.0
  33    REAL, PARAMETER :: qs0       = 6.e-4     ! pgaut
  34    REAL, PARAMETER :: g         = 9.81      !                                                         = 9.81
  35    REAL, PARAMETER :: rd        = 287.      ! gas constant for dry air (J/kg/K)                       = 287
  36    REAL, PARAMETER :: rv        = 461.6     ! gas constant for water vapor (J/kg/K)
  37    REAL, PARAMETER :: t0c       = 273.15    !                                                         = 273.15
  38    REAL, PARAMETER :: den0      = 1.28      ! density of 0 degree air      (kg/m^3)
  39    REAL, PARAMETER :: cpd       = 1004.5    ! heat capacity at constant pressure for dry air (J/kg/K) = 7.*rd/2.
  40    REAL, PARAMETER :: cpv       = 1846.4    ! heat capacity at constant pressure for vapor  (J/kg/K)  = 4.*r_v
  41 !   REAL, PARAMETER :: ep1       = 0.6083624 !                                                         = rv/rd-1.
  42    REAL, PARAMETER :: ep2       = 0.6217504 !                                                         = rd/rv
  43    REAL, PARAMETER :: qcrmin    = 1.e-9
  44    REAL, PARAMETER :: qmin      = 1.E-15    ! epsilon                                                 = 1.E-15
  45    REAL, PARAMETER :: xls       = 2.85E6    ! latent heat of sublimation   (J/kg)                     = 2.85E6
  46    REAL, PARAMETER :: xlv0      = 2.5E6     ! latent heat of vaporization  (J/kg)                     = 3.15E6
  47    REAL, PARAMETER :: xlf0      = 3.50E5    ! latent heat of melting       (J/kg)                     = 3.50E5
  48    REAL, PARAMETER :: cliq      = 4190.     !                                                         = 4190.
  49    REAL, PARAMETER :: cice      = 2106.     !                                                         = 2106
  50    REAL, PARAMETER :: psat      = 610.78    !                                                         = 610.78
  51    REAL, PARAMETER :: denr      = 1000.     ! water density                                           = 1000 (kg/m^3)
  52
  53
  54    REAL, SAVE :: pi , &
  55                  qc0     , qck1       , &
  56                  bvtr1   , bvtr2      , bvtr3    , bvtr4   , bvtr6   , &
  57                  g1pbr   , g3pbr      , g4pbr    , g5pbro2 , g6pbr   , pvtr    , &
  58                  bvts1   , bvts2      , bvts3    , bvts4   , &
  59                  g1pbs   , g3pbs      , g4pbs    , g5pbso2 , pvts    , &
  60                  bvtg1   , bvtg2      , bvtg3    , bvtg4   , &
  61                  g1pbg   , g3pbg      , g4pbg    , g5pbgo2 , pvtg    , &
  62                  roqimax , pidn0r     , pidn0s   , pidn0g  , xlv1    , &
  63                  vt2i    , vt2r       , vt2s     , vt2g    , egs     , egi     , &
  64                  vt2r_a  , vt2s_a     , vt2g_a   , vt2i_a  , &
  65                  fallr_a , falls_a    , fallg_a  , falli_a , &
  66                  pgfrz_a  , diffac_a   , diffac_b , pidep_a  , &
  67                  pgacs_a , pgacs_b    , pgacs_c  , pgacs_d , &
  68                  pgacr_a , pgacr_b    , pgacr_c  , pgacr_d , &
  69                  psacr_a , psacr_b    , psacr_c  , psacr_d , &
  70                  pracs_a , pracs_b    , pracs_c  , pracs_d , &
  71                  pgaci_a , pgaci_b    , pgaci_c  , pgaci_d , &
  72                  psevp_a , psevp_b    , pgevp_a  , pgevp_b , &
  73                  psmlt_a , psmlt_b    , pgmlt_a  , pgmlt_b , &
  74                  prevp_a , prevp_b    , psdep_a  , psdep_b , &
  75                  pgdep_a , pgdep_b    , &
  76                  praci_a , praci_b    , praci_c  , praci_d , &
  77                  psaci_a , psaci_b    , psaci_c  , psaci_d , &
  78                  pracw_a , piacr_a    , psacw_a  , pgacw_a
  79
  80 CONTAINS
  81
  82 !=======================================================================
  83 !
  84 !=======================================================================
  85    SUBROUTINE wsm6r(th, q, qc, qr, qi, qs, qg                        &
  86                   ,den, pii, p, delz                                &
  87                   ,delt                                             &
  88                   ,rain, rainncv                                    &
  89                   ,ids,ide, jds,jde, kds,kde                        &
  90                   ,ims,ime, jms,jme, kms,kme                        &
  91                   ,its,ite, jts,jte, kts,kte                        &
  92                                                                    )
  93 !-------------------------------------------------------------------
  94    IMPLICIT NONE
  95 !-------------------------------------------------------------------
  96    INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde,  &
  97                                        ims,ime, jms,jme, kms,kme,  &
  98                                        its,ite, jts,jte, kts,kte
  99    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
 100          INTENT(INOUT) ::                                          &
 101                                                                th, &
 102                                                                q , &
 103                                                                qc, &
 104                                                                qi, &
 105                                                                qr, &
 106                                                                qs, &
 107                                                                qg
 108    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
 109          INTENT(IN   ) ::                                          &
 110                                                               den, &
 111                                                               pii, &
 112                                                                 p, &
 113                                                              delz
 114    REAL, INTENT(IN   ) ::                                    delt
 115
 116    REAL, DIMENSION( ims:ime , jms:jme ),                           &
 117          INTENT(INOUT) ::                                    rain, &
 118                                                           rainncv
 119 !  LOCAL VAR
 120    REAL, DIMENSION( its:ite , kts:kte    ) ::  t
 121    REAL, DIMENSION( its:ite , kts:kte, 2 ) ::  qci
 122    REAL, DIMENSION( its:ite , kts:kte, 3 ) ::  qrs
 123    REAL, DIMENSION( ims:ime , kms:kme    ) ::  q2d, den2d, p2d, delz2d
 124    REAL, DIMENSION( ims:ime              ) ::  r1d, rcv1d
 125    REAL    ::   delt1
 126    INTEGER ::               i,j,k,ierr
 127
 128    delt1=delt
 129    DO j=jts,jte
 130       DO i=its,ite
 131          r1d  (i) = rain(i,j)
 132          rcv1d(i) = rainncv(i,j)
 133          DO k=kts,kte
 134             t      ( i,k   ) = th   ( i,k,j ) *pii ( i,k,j )
 135             qci    ( i,k,1 ) = qc   ( i,k,j )
 136             qci    ( i,k,2 ) = qi   ( i,k,j )
 137             qrs    ( i,k,1 ) = qr   ( i,k,j )
 138             qrs    ( i,k,2 ) = qs   ( i,k,j )
 139             qrs    ( i,k,3 ) = qg   ( i,k,j )
 140             q2d    ( i,k   ) = q    ( i,k,j )
 141             den2d  ( i,k   ) = den  ( i,k,j )
 142             p2d    ( i,k   ) = p    ( i,k,j )
 143             delz2d ( i,k   ) = delz ( i,k,j )
 144          ENDDO
 145       ENDDO
 146       !  Sending array starting locations of optional variables may cause
 147       !  troubles, so we explicitly change the call.
 148       CALL wsm62D(t, q2d, qci, qrs         &
 149                  ,den2d                    &
 150                  ,p2d, delz2d              &
 151                  ,delt1                    &
 152                  ,r1d,rcv1d                &
 153                  ,ims,ime, kms,kme         &
 154                  ,its,ite, kts,kte         )
 155
 156       DO I=its,ite
 157          rain(i,j)   =r1d(i)
 158          rainncv(i,j)=rcv1d(i)
 159          DO K=kts,kte
 160             th(i,k,j) = t(i,k)/pii(i,k,j)
 161             qc(i,k,j) = qci(i,k,1)
 162             qi(i,k,j) = qci(i,k,2)
 163             qr(i,k,j) = qrs(i,k,1)
 164             qs(i,k,j) = qrs(i,k,2)
 165             qg(i,k,j) = qrs(i,k,3)
 166             q (i,k,j) = q2d(i,k)
 167          ENDDO
 168       ENDDO
 169    ENDDO
 170
 171    END SUBROUTINE
 172
 173 !=======================================================================
 174 !
 175 !=======================================================================
 176    SUBROUTINE wsm62D(t, q, qci, qrs, den, p, delz ,delt ,rain,rainncv      &
 177                    ,ims,ime, kms,kme ,its,ite, kts,kte )
 178 !-------------------------------------------------------------------
 179    IMPLICIT NONE
 180 !-------------------------------------------------------------------
 181    integer ::   ims,ime, kms,kme,its,ite, kts,kte
 182
 183    REAL, DIMENSION( its:ite , kts:kte ),                           &
 184          INTENT(INOUT) ::                                          &
 185                                                                 t
 186    REAL, DIMENSION( its:ite , kts:kte, 2 ),                        &
 187          INTENT(INOUT) ::                                          &
 188                                                               qci
 189    REAL, DIMENSION( its:ite , kts:kte, 3 ),                        &
 190          INTENT(INOUT) ::                                          &
 191                                                               qrs
 192    REAL, DIMENSION( ims:ime , kms:kme ),                           &
 193          INTENT(INOUT) ::                                          &
 194                                                                 q
 195    REAL, DIMENSION( ims:ime , kms:kme ),                           &
 196          INTENT(IN   ) ::                                     den, &
 197                                                                 p, &
 198                                                              delz
 199    REAL, INTENT(IN   ) ::                                    delt
 200
 201    REAL, DIMENSION( ims:ime ),                                     &
 202          INTENT(INOUT) ::                                    rain, &
 203                                                           rainncv
 204 !  LOCAL VAR
 205    REAL, DIMENSION( its:ite , kts:kte , 3) ::                      &
 206          rh, qs, rslope, rslope2, rslope3, rslopeb,                &
 207          falk, fall ,work1
 208    REAL, DIMENSION( its:ite , kts:kte) ::                          &
 209          pracw, psacw, pgacw, pgacr, pgacs, psaci,  praci,    &
 210          piacr, pracs, psacr, pgaci, pseml, pgeml, fallc, &
 211          praut, psaut, pgaut, prevp, psdep, pgdep
 212    REAL, DIMENSION( its:ite , kts:kte ) ::                         &
 213          pigen, pidep, pcond, xl, cpm,  psevp,             &
 214          xni, pgevp,n0sfac,work2
 215 !   LOGICAL, DIMENSION( its:ite ) :: flgcld
 216    REAL  ::  dtcld,temp,temp0,supcol,supsat,satdt,eacrs,xmi,diameter,delta2,delta3
 217    INTEGER :: i,  k,  loop, loops
 218    real ::hsub, hvap, cvap, ttp, dldt, xa, xb, dldti, &
 219            xai ,xbi, tr, qs10 , qs11, qs20, qs21
 220    real :: fq, fqc, fqi, fqr, fqs, fqg, fallsum
 221
 222 !=================================================================
 223 !
 224    call wsm6rinit
 225 !----------------------------------------------------------------
 226 !  paddint 0 for negative values generated by dynamics
 227 !
 228    do k = kts, kte
 229       do i = its, ite
 230          q  (i,k  )=max(q  (i,k  ),0.)
 231          qci(i,k,1)=max(qci(i,k,1),0.)
 232          qrs(i,k,1)=max(qrs(i,k,1),0.)
 233          qci(i,k,2)=max(qci(i,k,2),0.)
 234          qrs(i,k,2)=max(qrs(i,k,2),0.)
 235          qrs(i,k,3)=max(qrs(i,k,3),0.)
 236       enddo
 237    enddo
 238 !
 239 !----------------------------------------------------------------
 240 !  compute the minor time steps.
 241 !
 242    loops = max(nint(delt/dtcldcr),1)
 243    dtcld = delt/loops
 244    if(delt.le.dtcldcr) dtcld = delt
 245 !
 246
 247    do loop = 1,loops
 248 !
 249 !----------------------------------------------------------------
 250 !     initialize the variables for microphysical physics
 251       call inimp(prevp,psdep,pgdep,praut,psaut,pgaut,pracw,praci,piacr,psaci,                 &
 252                   psacw,pracs,psacr,pgacw,pgaci,pgacr,pgacs,pigen,    &
 253                   pidep,pcond,pseml,pgeml,psevp,pgevp,falk,fall,fallc, xni,     &
 254                   kts, kte,its, ite)
 255 !----------------------------------------------------------------
 256 !     compute the fallout term:
 257 !     first, vertical terminal velosity for minor loops
 258
 259       call fallk(cpm,t,p,q,den,qrs,     &
 260                  delz,dtcld,falk,fall,                       &
 261                  kte, kts,its, ite,kme, kms,ims, ime)
 262
 263       call fallkc(qci,fallc,den,delz,dtcld,           &
 264                   kte,kts, its,ite, kme,kms, ims,ime)
 265
 266       call rainsc(fall,fallc,xl,t,q,qci,cpm,den,qrs, &
 267                   delz,rain,rainncv,dtcld,                  &
 268                   kte, kts,its, ite,kme, kms,ims, ime)
 269
 270       call warmr(t,q,qci,qrs,den,p,dtcld,xl,rh,qs,praut,pracw,&
 271                  prevp,ims,ime,kms,kme,its,ite,kts,kte)
 272 !
 273 ! cold rain processes
 274 !
 275 !
 276       call accret1(qci,den,qrs,t,q,     &
 277                    dtcld,praci,piacr,psaci,pgaci,psacw,pgacw,       &
 278                    ims,ime, kms,kme,its,ite,kts,kte)
 279
 280       call accret2(qrs,t,q,den,dtcld,  &
 281                    psacw,pgacw,pracs,psacr,pgacr,     &
 282                    pgacs,pseml,pgeml,                               &
 283                    ims,ime, kms,kme, its,ite, kts,kte)
 284
 285       call accret3(qrs,qci,rh,t,p,den,dtcld,   &
 286                    q,qs,psdep,pgdep,pigen,psaut,pgaut,psevp,     &
 287                    pgevp,pidep,                                       &
 288                    ims,ime, kms,kme,its,ite,kts,kte)
 289
 290       call pconadd(t,p,q,qci,qs,xl,cpm,dtcld,                       &
 291                       kte, kts,its, ite,kme, kms,ims, ime)
 292
 293    enddo                  ! big loops
 294 !
 295    END SUBROUTINE wsm62d
 296
 297 !=======================================================================
 298 !
 299 !=======================================================================
 300    subroutine calcrh(t,p,q,rh,qs)
 301
 302    REAL, INTENT(IN) :: t, q, p
 303    REAL, DIMENSION(3),  INTENT(OUT  ) :: rh, qs
 304
 305    real :: tr, qs10 , qs11, qs20, qs21
 306    real, parameter :: hsub = xls
 307    real, parameter :: hvap = xlv0
 308    real, parameter :: cvap = cpv
 309    real, parameter :: ttp  = t0c+0.01
 310    real, parameter :: dldt = cvap-cliq
 311    real, parameter :: xa   = -dldt/rv
 312    real, parameter :: xb   = xa+hvap/(rv*ttp)
 313    real, parameter :: dldti= cvap-cice
 314    real, parameter :: xai  = -dldti/rv
 315    real, parameter :: xbi  = xai+hsub/(rv*ttp)
 316    tr=ttp/t
 317    qs10  = psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
 318    qs11  = ep2 * qs10 / (p - qs10)
 319    qs(1) = qs11
 320    rh(1) = q / max(qs(1),qmin)
 321
 322    qs20 = psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
 323    qs21 = ep2 * qs20 / (p - qs20)
 324    qs(2) = qs21
 325    rh(2) = q / max(qs(2),qmin)
 326
 327    end subroutine
 328 !
 329 !=======================================================================
 330 !
 331 !=======================================================================
 332    SUBROUTINE wsm6rinit
 333 !-------------------------------------------------------------------
 334    IMPLICIT NONE
 335 !-------------------------------------------------------------------
 336 !.... constants which may not be tunable
 337    pi          = 4.*atan(1.)
 338    xlv1        = cliq-cpv
 339    qc0         = 4./3.*pi*denr*r0**3*xncr/den0  ! 0.419e-3 -- .61e-3
 340    qck1        = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
 341    bvtr1       = 1.+bvtr
 342    bvtr2       = 2.5+.5*bvtr
 343    bvtr3       = 3.+bvtr
 344    bvtr4       = 4.+bvtr
 345    bvtr6       = 6.+bvtr
 346    g1pbr       = rgmma(bvtr1)
 347    g3pbr       = rgmma(bvtr3)
 348    g4pbr       = rgmma(bvtr4)          ! 17.837825
 349    g6pbr       = rgmma(bvtr6)
 350    g5pbro2     = rgmma(bvtr2)          ! 1.8273
 351    pvtr        = avtr*g4pbr/6.
 352    roqimax     = 2.08e22*dimax**8
 353 !
 354    bvts1       = 1.+bvts
 355    bvts2       = 2.5+.5*bvts
 356    bvts3       = 3.+bvts
 357    bvts4       = 4.+bvts
 358    g1pbs       = rgmma(bvts1)    !.8875
 359    g3pbs       = rgmma(bvts3)
 360    g4pbs       = rgmma(bvts4)    ! 12.0786
 361    g5pbso2     = rgmma(bvts2)
 362    pvts        = avts*g4pbs/6.
 363    pidn0r      = pi*denr*n0r
 364    pidn0s      = pi*dens*n0s
 365
 366    bvtg1       = 1.+bvtg
 367    bvtg2       = 2.5+.5*bvtg
 368    bvtg3       = 3.+bvtg
 369    bvtg4       = 4.+bvtg
 370    g1pbg       = rgmma(bvtg1)
 371    g3pbg       = rgmma(bvtg3)
 372    g4pbg       = rgmma(bvtg4)
 373    g5pbgo2     = rgmma(bvtg2)
 374    pvtg        = avtg*g4pbg/6.
 375    pidn0g      = pi*deng*n0g
 376
 377    vt2r_a=pvtr*((pidn0r)**(-bvtr/4.))*sqrt(den0)
 378    vt2s_a=pvts*((pidn0s)**(-bvts/4.))*sqrt(den0)
 379    vt2g_a=pvtg*((pidn0g)**(-bvtg/4.))*sqrt(den0)
 380    vt2i_a=3.3
 381    fallr_a=vt2r_a
 382    falls_a=vt2s_a
 383    fallg_a=vt2g_a
 384    falli_a=vt2i_a
 385
 386    prevp_a=1.56*pi*n0r/sqrt(pidn0r)
 387    prevp_b=130.37*pi*sqrt(avtr)*n0r*((pidn0r)**(-(5.+bvtr)/8.))*sqrt(sqrt(den0))*g5pbro2
 388
 389    psdep_a=2.6*n0s/sqrt(pidn0s)
 390    psdep_b=370.08*sqrt(avts)*n0s*((pidn0s)**(-(5.+bvts)/8.))*sqrt(sqrt(den0))*g5pbso2
 391    psevp_a=psdep_a
 392    psevp_b=psdep_b
 393
 394    pgdep_a=1.56*pi*n0g/sqrt(pidn0g)
 395    pgdep_b=130.37*pi*sqrt(avtg)*n0g*((pidn0g)**(-(5.+bvtg)/8.))*sqrt(sqrt(den0))*g5pbgo2
 396    pgevp_a=pgdep_a
 397    pgevp_b=pgdep_b
 398
 399    psmlt_a=2.75e-3*pi*n0s/sqrt(pidn0s)/xlf0
 400    psmlt_b=0.391*pi*n0s*sqrt(sqrt(den0))*sqrt(avts)*((pidn0s)**(-(5.+bvts)/8.))*g5pbso2/xlf0
 401
 402    pgmlt_a= 3.3e-3*pi*n0g/sqrt(pidn0g)/xlf0
 403    pgmlt_b=0.276*pi*n0g*sqrt(sqrt(den0))*sqrt(avtg)*((pidn0g)**(-(5.+bvtg)/8.))*g5pbgo2/xlf0
 404
 405    praci_a=pi*n0r/4.
 406    praci_b=2./((pidn0r)**(3./4.))
 407    praci_c=3.245e-3/sqrt(pidn0r)
 408    praci_d=2.633e-6/sqrt(sqrt(pidn0r))
 409
 410    psaci_a=pi*n0s/4.
 411    psaci_b=2./((pidn0s)**(3./4.))
 412    psaci_c=3.245e-3/sqrt(pidn0s)
 413    psaci_d=2.633e-6/sqrt(sqrt(pidn0s))
 414
 415    pgaci_a=pi*n0g/4.
 416    pgaci_b=2./((pidn0g)**(3./4.))
 417    pgaci_c=3.245e-3/sqrt(pidn0g)
 418    pgaci_d=2.633e-6/sqrt(sqrt(pidn0g))
 419
 420    pracs_a=pi*n0r*pidn0s
 421    pracs_b=5./((pidn0s)**(3./2.))/sqrt(sqrt(pidn0r))
 422    pracs_c=2./((pidn0s)**(5./4.))/sqrt(pidn0r)
 423    pracs_d=.5/(pidn0s)/((pidn0r)**(3./4.))
 424
 425    psacr_a=pi*n0s*pidn0r
 426    psacr_b=5./((pidn0r)**(3./2.))/sqrt(sqrt(pidn0s))
 427    psacr_c=2./((pidn0r)**(5./4.))/sqrt(pidn0s)
 428    psacr_d=.5/(pidn0r)/((pidn0s)**(3./4.))
 429
 430    pgacr_a=pi*n0g*pidn0r
 431    pgacr_b=5./((pidn0r)**(3./2.))/sqrt(sqrt(pidn0g))
 432    pgacr_c=2./((pidn0r)**(5./4.))/sqrt(pidn0g)
 433    pgacr_d=.5/(pidn0r)/((pidn0g)**(3./4.))
 434
 435    pgacs_a=pi*n0g*pidn0s
 436    pgacs_b=5./((pidn0s)**(3./2.))/sqrt(sqrt(pidn0g))
 437    pgacs_c=2./((pidn0s)**(5./4.))/sqrt(pidn0g)
 438    pgacs_d=.5/(pidn0s)/((pidn0g)**(3./4.))
 439
 440    pidep_a=3.4927e5
 441
 442    diffac_a=4.7274e2
 443    diffac_b=1.1371e4
 444
 445    pgfrz_a=20.*pi*pfrz1/((pidn0r)**(3./4.))
 446
 447    piacr_a=5.38e7*pi*avtr*pidn0r*g6pbr*sqrt(den0)*((pidn0r)**(-(6.+bvtr)/4.))/24.
 448    pracw_a=.25*pi*avtr*n0r*g3pbr*sqrt(den0)*((pidn0r)**(-(3.+bvtr)/4.))
 449    psacw_a=.25*pi*avts*n0s*g3pbs*sqrt(den0)*((pidn0s)**(-(3.+bvts)/4.))
 450    pgacw_a=.25*pi*avtg*n0g*g3pbg*sqrt(den0)*((pidn0g)**(-(3.+bvtg)/4.))
 451    END SUBROUTINE wsm6rinit
 452
 453 !=======================================================================
 454 !
 455 !=======================================================================
 456    subroutine inimp(prevp,psdep,pgdep,praut,psaut,pgaut,pracw,praci,piacr,psaci,            &
 457           psacw,pracs,psacr,pgacw,pgaci,pgacr,pgacs,pigen,pidep,   &
 458           pcond,pseml,pgeml,psevp,pgevp,falk,fall,fallc,  &
 459           xni,kts, kte,its, ite)
 460
 461 !     initialize the variables for microphysical physics
 462    IMPLICIT NONE
 463
 464    integer :: kts, kte,its, ite, k, i
 465    REAL, DIMENSION( its:ite , kts:kte , 3) :: falk, fall
 466    REAL, DIMENSION( its:ite , kts:kte ) ::                         &
 467         xni, pgevp ,  pigen, pidep, pcond, fallc,             &
 468         pracw, psacw, pgacw, pgacr, pgacs, psaci, praci,    &
 469         piacr, pracs, psacr, pgaci, pseml, pgeml , psevp, &
 470         praut, psaut, pgaut, prevp, psdep, pgdep
 471
 472    do k = kts, kte
 473       do i = its, ite
 474          prevp(i,k) = 0.
 475          psdep(i,k) = 0.
 476          pgdep(i,k) = 0.
 477          praut(i,k) = 0.
 478          psaut(i,k) = 0.
 479          pgaut(i,k) = 0.
 480          pracw(i,k) = 0.
 481          praci(i,k) = 0.
 482          piacr(i,k) = 0.
 483          psaci(i,k) = 0.
 484          psacw(i,k) = 0.
 485          pracs(i,k) = 0.
 486          psacr(i,k) = 0.
 487          pgacw(i,k) = 0.
 488          pgaci(i,k) = 0.
 489          pgacr(i,k) = 0.
 490          pgacs(i,k) = 0.
 491          pigen(i,k) = 0.
 492          pidep(i,k) = 0.
 493          pcond(i,k) = 0.
 494          pseml(i,k) = 0.
 495          pgeml(i,k) = 0.
 496          psevp(i,k) = 0.
 497          pgevp(i,k) = 0.
 498          falk(i,k,1) = 0.
 499          falk(i,k,2) = 0.
 500          falk(i,k,3) = 0.
 501          fall(i,k,1) = 0.
 502          fall(i,k,2) = 0.
 503          fall(i,k,3) = 0.
 504          fallc(i,k) = 0.
 505          xni(i,k) = 1.e3
 506       enddo
 507    enddo
 508    end subroutine inimp
 509
 510 !=======================================================================
 511 !
 512 !=======================================================================
 513    subroutine fallk(cpm,t,p,q,den,qrs,delz,dtcld, &
 514                     falk,fall,kte, kts,its, ite,kme, kms,ims, ime)
 515    implicit none
 516    integer :: kte, kts,its, ite,kme, kms,ims, ime
 517    REAL, DIMENSION( its:ite , kts:kte , 3) :: qrs , falk, fall, work1
 518    REAL, DIMENSION( ims:ime , kms:kme )  :: delz, den ,p,q
 519    REAL, DIMENSION( its:ite , kts:kte )  :: psmlt,pgmlt,t,work2,cpm
 520    INTEGER, DIMENSION( its:ite ) :: mstep, numdt
 521    real ::  dtcld,coeres1,coeres2,coeresi,coeresh,xlf,psmlt0,pgmlt0,help_i,help_h,w1
 522    real ::  tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8
 523    integer :: mstepmax ,k, i,n,nw,jj
 524    real :: fqs, fqg, supcol, a, b, c, d
 525
 526    mstep=1
 527    mstepmax = 1
 528    numdt = 1
 529    do k = kte, kts, -1
 530       do i = its, ite
 531          work1(i,k,1) = vt2r_a*(den(i,k)**((bvtr-2.)/4.))*(max(qcrmin,qrs(i,k,1))**(bvtr/4.))/delz(i,k)
 532          work1(i,k,2) = vt2s_a*(den(i,k)**((bvts-2.)/4.))*(max(qcrmin,qrs(i,k,2))**(bvts/4.))/delz(i,k) &
 533                        *exp(-bvts*alpha*max(0.,min(90.,(t0c-t(i,k))))/4.)
 534          work1(i,k,3) = vt2g_a*(den(i,k)**((bvtg-2.)/4.))*(max(qcrmin,qrs(i,k,3))**(bvtg/4.))/delz(i,k)
 535          if(work1(i,k,1) .ge. work1(i,k,2) .and. work1(i,k,1) .ge. work1(i,k,3)) then
 536             w1=work1(i,k,1)
 537          elseif(work1(i,k,2) .ge. work1(i,k,1) .and. work1(i,k,2) .ge. work1(i,k,3)) then
 538             w1=work1(i,k,2)
 539          else
 540             w1=work1(i,k,3)
 541          endif
 542          nw=nint(w1*dtcld+.5)
 543          if(nw.gt.1) then
 544             numdt(i) = nw
 545          else
 546             numdt(i) = 1
 547          endif
 548
 549          if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
 550       enddo
 551    enddo
 552    do i = its, ite
 553       if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
 554    enddo
 555
 556    do n = 1, mstepmax
 557       do i = its, ite
 558          if(n.le.mstep(i)) then
 559             k=kte
 560             falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
 561             falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
 562             falk(i,k,3) = den(i,k)*qrs(i,k,3)*work1(i,k,3)/mstep(i)
 563
 564             fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
 565             fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
 566             fall(i,k,3) = fall(i,k,3)+falk(i,k,3)
 567
 568             do jj=1,3
 569                tmp1=min(falk(i,k,jj)*dtcld/den(i,k),qrs(i,k,jj))
 570                if(abs(tmp1)<qmin)then
 571                   tmp1=0.
 572                endif
 573                qrs(i,k,jj)=qrs(i,k,jj)-tmp1
 574             enddo
 575          endif
 576       enddo
 577       do k = kte-1, kts, -1
 578          do i = its, ite
 579             if(n.le.mstep(i)) then
 580                falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
 581                falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
 582                falk(i,k,3) = den(i,k)*qrs(i,k,3)*work1(i,k,3)/mstep(i)
 583
 584                fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
 585                fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
 586                fall(i,k,3) = fall(i,k,3)+falk(i,k,3)
 587                do jj=1,3
 588                   tmp2=min((falk(i,k,jj)-falk(i,k+1,jj)*delz(i,k+1)/delz(i,k))*dtcld/den(i,k),qrs(i,k,jj))
 589                    if(abs(tmp2)<qmin)then
 590                      tmp2=0.
 591                    endif
 592                    qrs(i,k,jj) = qrs(i,k,jj) - tmp2
 593
 594                enddo
 595
 596             endif
 597          enddo
 598       enddo
 599
 600
 601       do k = kte, kts, -1
 602          do i = its, ite
 603             if(n.le.mstep(i)) then !
 604    !---------------------------------------------------------------
 605    ! psmlt: melting of snow [RH83 A25]
 606    !       (T>T0: S->R) psmlt<0: min=-qrs(i,k,2), max=0
 607    !---------------------------------------------------------------
 608                !update xl, cpm
 609                cpm(i,k)=cpmcal(q(i,k))
 610                xlf = xlf0
 611                a=exp(alpha*max(0.,min(90.,(t0c-t(i,k))))/2.)
 612                b=exp(alpha*max(0.,min(90.,(t0c-t(i,k))))*(3-bvts)/8.)
 613                c=(t(i,k)**1.5      )*(t0c-t(i,k))/ (t(i,k)+120.)
 614                d=(t(i,k)**(3.88/6.))*(t0c-t(i,k))/((t(i,k)+120.)**(5./6.))
 615                psmlt0 = psmlt_a*a*c*sqrt(den(i,k)*max(qrs(i,k,2),qcrmin)) &
 616                       +psmlt_b*b*d*(p(i,k)**(1./3.))*(den(i,k)**((13.+3*bvts)/24.))*(max(qrs(i,k,2),qcrmin)**((5.+bvts)/8.))
 617                tmp3=psmlt0*dtcld/mstep(i)
 618                tmp4=-qrs(i,k,2)/mstep(i)
 619                if(tmp3 .gt.tmp4) then
 620                   tmp5=tmp3
 621                else
 622                   tmp5=tmp4
 623                endif
 624                if(tmp5 .lt. 0.) then
 625                   psmlt(i,k) = tmp5
 626                else
 627                   psmlt(i,k) = 0.
 628                endif
 629                if(abs(psmlt(i,k))<qmin)then
 630                   psmlt(i,k)=0.
 631                endif
 632                !update qs, qr, t
 633                qrs(i,k,2) = max(qrs(i,k,2) + psmlt(i,k),0.)
 634                qrs(i,k,1) = max(qrs(i,k,1) - psmlt(i,k),0.)
 635                t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
 636
 637             endif
 638          enddo
 639       enddo
 640    !---------------------------------------------------------------
 641    ! pgmlt: melting of graupel [LFO 47]
 642    !       (T>T0: G->R) pgmlt<0: min=-qrs(i,k,3), max=0
 643    !---------------------------------------------------------------
 644       do k = kte, kts, -1
 645          do i = its, ite
 646             if(n.le.mstep(i)) then
 647                !update xl, cpm
 648                xlf = xlf0
 649 !               cpm(i,k)=cpmcal(q(i,k)) ! not change
 650                c=(t(i,k)**1.5      )*(t0c-t(i,k))/ (t(i,k)+120.)
 651                d=(t(i,k)**(3.88/6.))*(t0c-t(i,k))/((t(i,k)+120.)**(5./6.))
 652                pgmlt0 = pgmlt_a*c*sqrt(den(i,k)*max(qrs(i,k,3),qcrmin)) &
 653                       +pgmlt_b*d*(p(i,k)**(1./3.))*(den(i,k)**((13.+3*bvtg)/24.))*(max(qrs(i,k,3),qcrmin)**((5.+bvtg)/8.))
 654                tmp6=pgmlt0*dtcld/mstep(i)
 655                tmp7=-qrs(i,k,3)/mstep(i)
 656                if(tmp6 .gt.tmp7) then
 657                   tmp8=tmp6
 658                else
 659                   tmp8=tmp7
 660                endif
 661                if(tmp8 .lt. 0.) then
 662                   pgmlt(i,k) = tmp8
 663                else
 664                   pgmlt(i,k) = 0.
 665                endif
 666                if(abs(pgmlt(i,k))<qmin)then
 667                   pgmlt(i,k)=0.
 668                endif
 669                !update qg, qr, t
 670                qrs(i,k,3) = max(qrs(i,k,3) + pgmlt(i,k),0.)
 671                qrs(i,k,1) = max(qrs(i,k,1) - pgmlt(i,k),0.)
 672                t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
 673             endif
 674          enddo
 675       enddo
 676    enddo
 677
 678    end subroutine fallk
 679
 680 !=======================================================================
 681 !
 682 !=======================================================================
 683    subroutine fallkc(qci,fallc,den,delz,dtcld, kte, kts,its, ite,kme, kms,ims, ime)
 684    implicit none
 685
 686    integer :: kte, kts,its, ite,kme, kms,ims, ime
 687    REAL, DIMENSION( its:ite , kts:kte, 2 ) ::  qci
 688    REAL, DIMENSION( ims:ime , kms:kme )  ::    delz, den
 689    REAL, DIMENSION( its:ite , kts:kte )  ::    falkc,work1c,work2c,xni,fallc
 690
 691    INTEGER, DIMENSION( its:ite ) :: mstep, numdt
 692    real :: dtcld ,xmi,diameter,temp1,temp2,temp3,temp4,temp5, temp0
 693    integer :: mstepmax ,k, i, n
 694
 695    mstepmax = 1
 696    mstep = 1
 697    numdt = 1
 698    do k = kte, kts, -1
 699       do i = its, ite
 700          work1c(i,k) = vt2i_a * ((den(i,k)*qci(i,k,2))**(1.31/8.))
 701          work2c(i,k) = work1c(i,k)/delz(i,k)
 702          numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
 703          if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
 704       enddo
 705    enddo
 706    do i = its, ite
 707       if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
 708    enddo
 709
 710    do n = 1, mstepmax
 711       k = kte
 712       do i = its, ite
 713          if(n.le.mstep(i)) then
 714             falkc(i,k) = falli_a*((den(i,k)*qci(i,k,2))**(9.31/8.))/delz(i,k)/mstep(i)
 715             fallc(i,k) = fallc(i,k)+falkc(i,k)
 716
 717             temp3=min(falkc(i,k)*dtcld/den(i,k),qci(i,k,2))
 718             if(abs(temp3)<qmin)then
 719                temp3=0.
 720             endif
 721             qci(i,k,2)=qci(i,k,2)-temp3
 722          endif
 723       enddo
 724       do k = kte-1, kts, -1
 725          do i = its, ite
 726             if(n.le.mstep(i)) then
 727                falkc(i,k) = falli_a*((den(i,k)*qci(i,k,2))**(9.31/8.))/delz(i,k)/mstep(i)
 728                fallc(i,k) = fallc(i,k)+falkc(i,k)
 729
 730                temp4=min((falkc(i,k)-falkc(i,k+1)*delz(i,k+1)/delz(i,k))*dtcld/den(i,k),qci(i,k,2))
 731                if(abs(temp4)<qmin)then
 732                   temp4=0.
 733                endif
 734                qci(i,k,2)=qci(i,k,2)-temp4
 735             endif
 736          enddo
 737       enddo
 738    enddo
 739    end subroutine fallkc
 740
 741 !=======================================================================
 742 !
 743 !=======================================================================
 744    subroutine rainsc(fall,fallc,xl,t,q,qci,cpm,den,qrs,delz,rain,rainncv,  &
 745                      dtcld,kte, kts,its, ite,kme, kms,ims, ime)
 746    implicit none
 747    integer :: kte, kts,its, ite,kme, kms,ims, ime
 748    REAL, DIMENSION( its:ite , kts:kte , 3) :: qrs , fall
 749    REAL, DIMENSION( its:ite , kts:kte , 2) :: qci
 750    REAL, DIMENSION( ims:ime , kms:kme )  :: delz,den ,q
 751    REAL, DIMENSION( its:ite , kts:kte )  :: xl,t,cpm, fallc
 752    real, DIMENSION( ims:ime ) :: rain,rainncv
 753    integer :: k, i
 754    real :: dtcld ,fallsum,supcol,xlf,temp,temp0,pfrzdtr,pfrzdtc
 755    real :: ft0, ft40, fsupcol, fqc, fqi, fqr, qtmp
 756
 757    do i = its, ite
 758       fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
 759       if(fallsum.gt.qmin) then
 760          !rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
 761          rainncv(i) = fallsum*delz(i,kts)/denr*dtcld*1000.
 762          rain   (i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rain(i)
 763       endif
 764    enddo
 765
 766    do k = kts, kte
 767       do i = its, ite
 768    !---------------------------------------------------------------
 769    ! pimlt: instantaneous melting of cloud ice [RH83 A28]
 770    !       (T>T0: I->C) pimlt=qci(i,k,2) t-
 771    !---------------------------------------------------------------
 772          !update xl, cpm
 773          xl (i,k)=xlcal (t(i,k))
 774 !         cpm(i,k)=cpmcal(q(i,k)) !not change
 775          xlf = xls-xl(i,k)
 776          supcol = t0c-t(i,k)
 777          if(supcol.lt.0.) xlf = xlf0
 778          call smoothif(t(i,k)    ,t0c,ft0,'t0')
 779          qtmp=ft0*max(qci(i,k,2),0.)
 780          if(abs(qtmp)<qmin)then
 781             qtmp=0.
 782          endif
 783          !update qc, qi, t
 784          qci(i,k,1) = max(qci(i,k,1) + qtmp,0.)
 785          qci(i,k,2) = max(qci(i,k,2) - qtmp,0.)
 786          t  (i,k  ) = t(i,k) - xlf/cpm(i,k)*qtmp
 787    !---------------------------------------------------------------
 788    ! pihmf: homogeneous freezing of cloud water below -40c
 789    !        (T<-40C: C->I) min=0,pihmf=qci(i,k,1) t+
 790    !---------------------------------------------------------------
 791          !update xl, cpm
 792          xl (i,k)=xlcal (t(i,k))
 793 !         cpm(i,k)=cpmcal(q(i,k)) !not change
 794          xlf = xls-xl(i,k)
 795          supcol = t0c-t(i,k)
 796          if(supcol.lt.0.) xlf = xlf0
 797          call smoothif(supcol   ,40.,ft40,'t0')
 798          qtmp=max(ft40*qci(i,k,1),0.)
 799          if(abs(qtmp)<qmin)then
 800             qtmp=0.
 801          endif            !update qc, qi, t
 802          qci(i,k,2) = max(qci(i,k,2) + qtmp,0.)
 803          qci(i,k,1) = max(qci(i,k,1) - qtmp,0.)
 804          t  (i,k  ) = t  (i,k  ) + xlf/cpm(i,k)*qtmp
 805    !---------------------------------------------------------------
 806    ! pihtf: heterogeneous freezing of cloud water
 807    !        (T0>T>-40C: C->I) max=qci(i,k,1),min=0. t+
 808    !---------------------------------------------------------------
 809          !update xl, cpm
 810          xl (i,k)=xlcal (t(i,k))
 811 !         cpm(i,k)=cpmcal(q(i,k)) !not change
 812          xlf = xls-xl(i,k)
 813          supcol = t0c-t(i,k)
 814          if(supcol.lt.0.) xlf = xlf0
 815          !t>-40C=t0c-40,t0c-t<40, supcol<40,-supcol>-40
 816          call smoothif(-supcol   ,-40.,ft40,'t0')
 817          ! t<0C, exp(pfrz2*supcol)-1.<0, pfrzdtc<0,
 818          pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.)*den(i,k)/denr &
 819                   /xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
 820          qtmp=max(ft40*pfrzdtc,0.)
 821          if(abs(qtmp)<qmin)then
 822             qtmp=0.
 823          endif            !update qc, qi, t
 824          qci(i,k,2) = max(qci(i,k,2) + qtmp,0.)
 825          qci(i,k,1) = max(qci(i,k,1) - qtmp,0.)
 826          t  (i,k  ) = t  (i,k  ) + xlf/cpm(i,k)*qtmp
 827    !---------------------------------------------------------------
 828    ! pgfrz: freezing of rain water [LFO 45]
 829    !        (T<T0, R->G) max=qrs(i,k,1),min=0. t+
 830    !---------------------------------------------------------------
 831          !update xl, cpm
 832          xl (i,k)=xlcal (t(i,k))
 833 !         cpm(i,k)=cpmcal(q(i,k))!not change
 834          xlf = xls-xl(i,k)
 835          supcol = t0c-t(i,k)
 836          if(supcol.lt.0.) xlf = xlf0
 837          if(qrs(i,k,1)>0.)then
 838             temp=pgfrz_a*(exp(pfrz2*supcol)-1.)*(den(i,k)**(3./4.))*(qrs(i,k,1)**(7./4.))
 839          else
 840             temp=0.
 841          endif
 842          pfrzdtr = min(temp*dtcld,qrs(i,k,1))
 843          qtmp=max(pfrzdtr,0.)
 844          if(abs(qtmp)<qmin)then
 845             qtmp=0.
 846          endif
 847          !update qg, qr, t
 848          qrs(i,k,3) = max(qrs(i,k,3) + qtmp,0.)
 849          qrs(i,k,1) = max(qrs(i,k,1) - qtmp,0.)
 850          t  (i,k  ) = t  (i,k  ) + xlf/cpm(i,k)*qtmp
 851       enddo
 852    enddo
 853    end subroutine rainsc
 854 !=======================================================================
 855 !
 856 !=======================================================================
 857    subroutine warmr(t, q, qci, qrs, den, p ,dtcld,xl,rh,qs,               &
 858                    praut,pracw,prevp,  &
 859                    ims,ime, kms,kme, its,ite, kts,kte)
 860 !-------------------------------------------------------------------
 861    IMPLICIT NONE
 862    integer ::  ims,ime, kms,kme,its,ite, kts,kte
 863 !------------------------------------------------------------------
 864    REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci
 865    REAL, DIMENSION( ims:ime , kms:kme ) :: q, den, p
 866    REAL, DIMENSION( its:ite , kts:kte , 3) ::                      &
 867          rh, qs, qrs,           &
 868           work1
 869    REAL, DIMENSION( its:ite , kts:kte) :: praut, prevp, pracw,xl,denfac,t,cpm
 870    REAL  ::  coeres,supsat,satdt,dtcld,praut1
 871    INTEGER :: i, k
 872    real :: fqv, fqc, fqr, fqc0, fprevp,prevp0,prevp1,temp, a, b, c, d, e
 873
 874    do k = kts, kte
 875       do i = its, ite
 876    !---------------------------------------------------------------
 877    ! praut: auto conversion rate from cloud to rain [HDC 16]
 878    !        (C->R) praut>0 max=qci(i,k,1)/dtcld, min=0.
 879    !---------------------------------------------------------------
 880          call smoothif(qci(i,k,1),qc0,fqc0,'q0') !qc0=5.03e-4
 881          if(qci(i,k,1)>0.)then ! x**a need x>0
 882             praut1 = fqc0*qck1*exp(log(qci(i,k,1))*(7./3.)) !(qci(i,k,1)**(7./3.))
 883          else
 884             praut1=0.
 885          endif
 886          praut(i,k) = min(praut1,qci(i,k,1)/dtcld)
 887          if(abs(praut(i,k))<qmin/dtcld)then
 888             praut(i,k)=0.
 889          endif
 890
 891          ! update qc,qr
 892          qci(i,k,1)=max(qci(i,k,1)-praut(i,k)*dtcld,0.)
 893          qrs(i,k,1)=max(qrs(i,k,1)+praut(i,k)*dtcld,0.)
 894          praut(i,k)=0.
 895
 896    !---------------------------------------------------------------
 897    ! pracw: accretion of cloud water by rain [LFO 51]
 898    !        (C->R) max=qci(i,k,1)/dtcld, min=0.
 899    !---------------------------------------------------------------
 900          !call smoothif(qrs(i,k,1),qcrmin,fqr,'q0')
 901          !call smoothif(qci(i,k,1),qmin  ,fqc,'q0')
 902          if(qrs(i,k,1)>0.and.qci(i,k,1)>0.)then
 903             pracw(i,k) = pracw_a*(den(i,k)**((1.+bvtr)/4.))*(qrs(i,k,1)**((3.+bvtr)/4.))*qci(i,k,1)
 904          else
 905             pracw(i,k) = 0.
 906          endif
 907          pracw(i,k) =max(min(pracw(i,k),qci(i,k,1)/dtcld),0.)
 908          if(abs(pracw(i,k))<qmin/dtcld)then
 909             pracw(i,k)=0.
 910          endif
 911          ! update qc,qr
 912          qci(i,k,1)=max(qci(i,k,1)-pracw(i,k)*dtcld,0.)
 913          qrs(i,k,1)=max(qrs(i,k,1)+pracw(i,k)*dtcld,0.)
 914          pracw(i,k)=0.
 915    !
 916    !---------------------------------------------------------------
 917    ! prevp: evaporation/condensation rate of rain [HDC 14]
 918    !        (V->R or R->V) rh(i,k,1)>1., prevp>0, V->R, min=0.,                max=satdt ;
 919    !                       rh(i,k,1)<1., prevp<0, R->V, min=-qrs(i,k,1)/dtcld, max=0.
 920    !---------------------------------------------------------------
 921          !update rh
 922          call calcrh(t(i,k),p(i,k),q(i,k),rh(i,k,:),qs(i,k,:))
 923          !update xl, cpm
 924          xl (i,k)=xlcal (t(i,k))
 925          cpm(i,k)=cpmcal(q(i,k))
 926
 927          supsat = q(i,k)-qs(i,k,1)
 928          satdt = supsat/dtcld
 929
 930          a=sqrt(den(i,k)*max(qrs(i,k,1),qcrmin))
 931          b=((t(i,k)+120.)**(1./6.))/(t(i,k)**(5.12/6.))*(p(i,k)**(1./3.))&
 932             *(den(i,k)**((13.+3.*bvtr)/24.))*(max(qrs(i,k,1),qcrmin)**((5.+bvtr)/8.))
 933          c=diffac_a*den(i,k)*xl(i,k)*xl(i,k)*(t(i,k)+120.)/rv/(t(i,k)**3.5)
 934          d=diffac_b*p(i,k)/(t(i,k)**1.81)/qs(i,k,1)
 935          e=(rh(i,k,1)-1.)/(c+d)
 936          prevp(i,k)=(prevp_a*a+prevp_b*b)*e
 937          if(prevp(i,k)<0.)then
 938             prevp(i,k) = min(max(prevp(i,k),-qrs(i,k,1)/dtcld),0.) ! R->V
 939          else
 940             prevp(i,k) = max(min(prevp(i,k),satdt),0.)             ! V->R
 941          endif
 942          if(abs(prevp(i,k))<qmin/dtcld)then
 943             prevp(i,k)=0.
 944          endif
 945          ! update q,qr,t
 946          q  (i,k  )=max(q  (i,k  )-prevp(i,k)*dtcld,0.)
 947          qrs(i,k,1)=max(qrs(i,k,1)+prevp(i,k)*dtcld,0.)
 948          t  (i,k  )=t  (i,k  )+prevp(i,k)*dtcld*xl(i,k)/cpm(i,k)
 949          prevp(i,k) = 0.
 950       enddo
 951    enddo
 952    end subroutine warmr
 953
 954
 955 !===================================================================
 956    subroutine accret1(qci,den,qrs,t,q,dtcld,  &
 957                       praci,piacr,psaci,pgaci,psacw,pgacw,                &
 958                       ims,ime, kms,kme,its,ite,kts,kte)
 959    IMPLICIT NONE
 960    integer ::         ims,ime, kms,kme,its,ite, kts,kte
 961 !-------------------------------------------------------------------
 962    REAL, DIMENSION( its:ite , kts:kte, 2 ) ::   qci
 963    REAL, DIMENSION( its:ite , kts:kte, 3 ) ::   qrs
 964    REAL, DIMENSION( ims:ime , kms:kme ) ::  den,q
 965    REAL, DIMENSION( its:ite , kts:kte) :: praci,piacr,psaci,pgaci,psacw,pgacw,t,xl,cpm
 966    REAL  ::  supcol,dtcld,eacrs,egi,praci1,piacr1,psaci1,pgaci1,temp,temp0
 967    INTEGER :: i, k
 968    real :: fsupcol, fqc, fqi, fqr, fqs, fqg, delta3, xlf, a, b, c, d, e
 969
 970    do k = kts, kte
 971      do i = its, ite
 972    !-------------------------------------------------------------
 973    ! praci: Accretion of cloud ice by rain [LFO 25]
 974    !        (T<T0: I->S or I->G) praci: min=0,max=qci(i,k,2)/dtcld
 975    !-------------------------------------------------------------
 976          supcol = t0c-t(i,k)
 977          call smoothif(supcol,    0.,fsupcol,'t0')
 978
 979          vt2r=vt2r_a*(den(i,k)**((bvtr-2.)/4.))*(max(qrs(i,k,1),qcrmin)**(bvtr/4.))
 980          vt2i=vt2i_a*((den(i,k)*max(qci(i,k,2),qmin))**(1.31/8.))
 981          b=((den(i,k)*max(qrs(i,k,1),qcrmin))**(3./4.))*max(qci(i,k,2),qmin)
 982          c=(den(i,k)**(5./8.))*sqrt(max(qrs(i,k,1),qcrmin)) *(max(qci(i,k,2),qmin)**(9./8.))
 983          d=sqrt(den(i,k))*sqrt(sqrt(max(qrs(i,k,1),qcrmin)))*(max(qci(i,k,2),qmin)**(5./4.))
 984          praci1=praci_a*abs(vt2r-vt2i)*(praci_b*b+praci_c*c+praci_d*d)
 985
 986          praci(i,k) = min(praci1,qci(i,k,2)/dtcld)
 987          praci(i,k)=fsupcol*praci(i,k)
 988          !update qi, qs, qg
 989          if(qrs(i,k,1).lt.1.e-4)then
 990             delta3=1.
 991          else
 992             delta3=0.
 993          endif
 994          if(abs(praci(i,k))<qmin/dtcld)then
 995             praci(i,k)=0.
 996          endif
 997          qci(i,k,2)=max(qci(i,k,2)-praci(i,k)           *dtcld,0.)
 998          qrs(i,k,2)=max(qrs(i,k,2)+praci(i,k)    *delta3*dtcld,0.)
 999          qrs(i,k,3)=max(qrs(i,k,3)+praci(i,k)*(1-delta3)*dtcld,0.)
1000          praci(i,k)=0.
1001
1002    !-------------------------------------------------------------
1003    ! piacr: Accretion of rain by cloud ice [LFO 26]
1004    !        (T<T0: R->S or R->G) piacr: min=0,max=qrs(i,k,1)/dtcld
1005    !-------------------------------------------------------------
1006 !         supcol = t0c-t(i,k) !not change
1007          call smoothif(supcol,    0.,fsupcol,'t0')
1008          !call smoothif(qci(i,k,2),qmin  ,fqi,'q0')
1009          !call smoothif(qrs(i,k,1),qcrmin,fqr,'q0')
1010          !update cpm
1011          cpm(i,k)=cpmcal(q(i,k))
1012          xl(i,k)=xlcal(t(i,k))
1013          xlf = xls-xl(i,k)
1014          if(supcol.lt.0.) xlf = xlf0
1015
1016          if(qci(i,k,2)>0..and.qrs(i,k,1)>0.)then
1017             !piacr_a=1.75e5
1018             piacr1=piacr_a*(den(i,k)**((3.+bvtr)/4.))*(qci(i,k,2)**0.75) &
1019                       *(qrs(i,k,1)**((6.+bvtr)/4.))
1020          else
1021             piacr1=0.
1022          endif
1023
1024          piacr(i,k) = min(piacr1,qrs(i,k,1)/dtcld)
1025          piacr(i,k)=fsupcol*piacr(i,k)
1026        ! update qr,qs,qg,t
1027          if(qrs(i,k,1).lt.1.e-4)then
1028             delta3=1.
1029          else
1030             delta3=0.
1031          endif
1032          if(abs(piacr(i,k))<qmin/dtcld)then
1033             piacr(i,k)=0.
1034          endif
1035          qrs(i,k,1)=max(qrs(i,k,1)-piacr(i,k)           *dtcld,0.)
1036          qrs(i,k,2)=max(qrs(i,k,2)+piacr(i,k)    *delta3*dtcld,0.)
1037          qrs(i,k,3)=max(qrs(i,k,3)+piacr(i,k)*(1-delta3)*dtcld,0.)
1038          t(i,k)=t(i,k)+piacr(i,k)*dtcld*xlf/cpm(i,k)
1039          piacr(i,k)=0.
1040    !-------------------------------------------------------------
1041    ! psaci: Accretion of cloud ice by snow [HDC 10]
1042    !        (T<T0: I->S) psaci: min=0, max=qci(i,k,2)/dtcld
1043    !-------------------------------------------------------------
1044          supcol = t0c-t(i,k)
1045          call smoothif(supcol,    0.,fsupcol,'t0')
1046          !call smoothif(qci(i,k,2),qmin  ,fqi,'q0')
1047          !call smoothif(qrs(i,k,2),qcrmin,fqs,'q0')
1048
1049          eacrs = min(exp(0.07*(-supcol)),1.)
1050
1051          vt2s=vt2s_a*(den(i,k)**((bvts-2.)/4.))*(max(qrs(i,k,2),qcrmin)**(bvts/4.))&
1052                *exp(-alpha*bvts*max(0.,min(90.,t0c-t(i,k)))/4.)
1053          vt2i=vt2i_a*((den(i,k)*max(qci(i,k,2),qmin))**(1.31/8.))
1054          a=exp(    alpha*max(0.,min(90.,t0c-t(i,k))))
1055          b=exp(-3.*alpha*max(0.,min(90.,t0c-t(i,k)))/4.)&
1056             *((den(i,k)*max(qrs(i,k,2),qcrmin))**(3./4.))*max(qci(i,k,2),qmin)
1057          c=exp(-   alpha*max(0.,min(90.,t0c-t(i,k)))/2.)&
1058             *(den(i,k)**(5./8.))*sqrt(max(qrs(i,k,2),qcrmin)) *(max(qci(i,k,2),qmin)**(9./8.))
1059          d=exp(-   alpha*max(0.,min(90.,t0c-t(i,k)))/4.)&
1060             *sqrt(den(i,k))*sqrt(sqrt(max(qrs(i,k,2),qcrmin)))*(max(qci(i,k,2),qmin)**(5./4.))
1061          psaci1=psaci_a*a*abs(vt2s-vt2i)*(psaci_b*b+psaci_c*c+psaci_d*d)*eacrs
1062
1063          psaci(i,k) = min(psaci1,qci(i,k,2)/dtcld)
1064          psaci(i,k)=fsupcol*psaci(i,k)
1065          if(abs(psaci(i,k))<qmin/dtcld)then
1066             psaci(i,k)=0.
1067          endif
1068          !update qi, qs
1069          qci(i,k,2)=max(qci(i,k,2)-psaci(i,k)*dtcld,0.)
1070          qrs(i,k,2)=max(qrs(i,k,2)+psaci(i,k)*dtcld,0.)
1071          psaci(i,k)=0.
1072
1073    !-------------------------------------------------------------
1074    ! pgaci: Accretion of cloud ice by graupel [LFO 41]
1075    !        (T<T0: I->G) pgaci:min=0,max=qci(i,k,2)/dtcld
1076    !-------------------------------------------------------------
1077 !         supcol = t0c-t(i,k) !not change
1078 !         call smoothif(supcol,    0.,fsupcol,'t0')
1079          !call smoothif(qci(i,k,2),qmin  ,fqi,'q0')
1080          !call smoothif(qrs(i,k,3),qcrmin,fqg,'q0')
1081          egi = eacrs !min(exp(0.07*(-supcol)),1.)
1082
1083          vt2g=vt2g_a*(den(i,k)**((bvtg-2.)/4.))*(max(qrs(i,k,3),qcrmin)**(bvtg/4.))
1084          vt2i=vt2i_a*((den(i,k)*max(qci(i,k,2),qmin))**(1.31/8.))
1085          b=((den(i,k)*max(qrs(i,k,3),qcrmin))**(3./4.))*max(qci(i,k,2),qmin)
1086          c=(den(i,k)**(5./8.))*sqrt(max(qrs(i,k,3),qcrmin) )*(max(qci(i,k,2),qmin)**(9./8.))
1087          d=sqrt(den(i,k))*sqrt(sqrt(max(qrs(i,k,3),qcrmin)))*(max(qci(i,k,2),qmin)**(5./4.))
1088          pgaci1=pgaci_a*abs(vt2g-vt2i)*(pgaci_b*b+pgaci_c*c+pgaci_d*d)*egi
1089
1090          pgaci(i,k) = min(pgaci1,qci(i,k,2)/dtcld)
1091          pgaci(i,k)=fsupcol*pgaci(i,k)
1092          if(abs(pgaci(i,k))<qmin/dtcld)then
1093             pgaci(i,k)=0.
1094          endif
1095
1096          ! update qi,qg
1097          qci(i,k,2)=max(qci(i,k,2)-pgaci(i,k)*dtcld,0.)
1098          qrs(i,k,3)=max(qrs(i,k,3)+pgaci(i,k)*dtcld,0.)
1099          pgaci(i,k)=0.
1100
1101    !-------------------------------------------------------------
1102    ! psacw: Accretion of cloud water by snow  [LFO 24]
1103    !        (T<T0: C->G, and T>=T0: C->R) psacw:min=0,max=qci(i,k,1)/dtcld
1104    !-------------------------------------------------------------
1105 !         supcol = t0c-t(i,k) !not change
1106 !         call smoothif(supcol,    0.,fsupcol,'t0')
1107          !call smoothif(qrs(i,k,2),qcrmin,fqs,'q0')
1108          !call smoothif(qci(i,k,1),qmin  ,fqc,'q0')
1109          !update cpm
1110 !         cpm(i,k)=cpmcal(q(i,k)) !not change
1111          xl(i,k)=xlcal(t(i,k))
1112          xlf = xls-xl(i,k)
1113          if(supcol.lt.0.) xlf = xlf0
1114
1115          if(qrs(i,k,2)>0..and.qci(i,k,1)>0.)then
1116             a=exp((1.-bvts)*alpha*max(0.,min(90.,t0c-t(i,k)))/4.)
1117             psacw(i,k)=psacw_a*a*(den(i,k)**((1.+bvts)/4.))*(qrs(i,k,2)**((3.+bvts)/4.))*qci(i,k,1)
1118          else
1119             psacw(i,k)=0.
1120          endif
1121          psacw(i,k) =max(min(psacw(i,k),qci(i,k,1)/dtcld),0.)
1122          psacw(i,k)=fsupcol*psacw(i,k)
1123          if(abs(psacw(i,k))<qmin/dtcld)then
1124             psacw(i,k)=0.
1125          endif
1126          ! update qc, qr, qg, t
1127          qci(i,k,1)=max(qci(i,k,1)-             psacw(i,k)*dtcld,0.)
1128          qrs(i,k,1)=max(qrs(i,k,1)+(1.-fsupcol)*psacw(i,k)*dtcld,0.) !t>=t0
1129          qrs(i,k,3)=max(qrs(i,k,3)+fsupcol     *psacw(i,k)*dtcld,0.) !t<t0
1130          t  (i,k  )=t  (i,k  )+fsupcol*psacw(i,k)*dtcld*xlf/cpm(i,k)
1131          !t>=t0 pseml
1132          psacw(i,k)=(1-fsupcol)*psacw(i,k)
1133    !-------------------------------------------------------------
1134    ! pgacw: Accretion of cloud water by graupel [LFO 40]
1135    !        (T<T0: C->G, and T>=T0: C->R) pgacw:min=0.,max=qci(i,k,1)/dtcld
1136    !-------------------------------------------------------------
1137          supcol = t0c-t(i,k)
1138          call smoothif(supcol,    0.,fsupcol,'t0')
1139          !call smoothif(qrs(i,k,3),qcrmin,fqg,'q0')
1140          !call smoothif(qci(i,k,1),qmin  ,fqc,'q0')
1141 !         cpm(i,k)=cpmcal(q(i,k)) !not change
1142          xl(i,k)=xlcal(t(i,k))
1143          xlf = xls-xl(i,k)
1144          if(supcol.lt.0.) xlf = xlf0
1145          if(qrs(i,k,3)>0..and.qci(i,k,1)>0.)then
1146             pgacw(i,k)=pgacw_a*(den(i,k)**((1.+bvtg)/4.))*(qrs(i,k,3)**((3.+bvtg)/4.))*qci(i,k,1)
1147          else
1148             pgacw(i,k)=0.
1149          endif
1150          pgacw(i,k) =max(min(pgacw(i,k),qci(i,k,1)/dtcld),0.)
1151          !pgacw(i,k)=fqg*fqc*pgacw(i,k)
1152          if(abs(pgacw(i,k))<qmin/dtcld)then
1153             pgacw(i,k)=0.
1154          endif
1155          ! update qc, qr, qg, t
1156          qci(i,k,1)=max(qci(i,k,1)-             pgacw(i,k)*dtcld,0.)
1157          qrs(i,k,1)=max(qrs(i,k,1)+(1.-fsupcol)*pgacw(i,k)*dtcld,0.) !t>=t0
1158          qrs(i,k,3)=max(qrs(i,k,3)+    fsupcol *pgacw(i,k)*dtcld,0.) !t<t0
1159          t  (i,k  )=t  (i,k  )+fsupcol*pgacw(i,k)*dtcld*xlf/cpm(i,k)
1160          ! t>=t0 pgeml
1161          pgacw(i,k)=(1-fsupcol)*pgacw(i,k)
1162       enddo
1163    enddo
1164
1165    end  subroutine accret1
1166
1167 !=======================================================================
1168 !
1169 !=======================================================================
1170    subroutine accret2(qrs,t,q, den,dtcld,   &
1171                      psacw,pgacw,pracs,psacr,pgacr,pgacs,pseml,pgeml,  &
1172                      ims,ime, kms,kme, its,ite, kts,kte)
1173    IMPLICIT NONE
1174    integer ::         ims,ime, kms,kme ,its,ite, kts,kte
1175 !-------------------------------------------------------------------
1176    REAL, DIMENSION( its:ite , kts:kte, 3 ) :: qrs
1177    REAL, DIMENSION( ims:ime , kms:kme ) :: den,q
1178    REAL, DIMENSION( its:ite , kts:kte) :: psacw,pgacw,pracs,psacr,pgacr,pgacs,pseml,pgeml,  &
1179                                           t, xl, cpm
1180    REAL  ::  supcol,vt2r,vt2s,vt2g,dtcld,xlf,egs
1181    REAL  ::  acrfac1,acrfac2,acrfac3,acrfac4,pracs1,psacr1,pgacr1,pgacs1
1182    INTEGER :: i,  k
1183    real :: fsupcol, ft0, fqr, fqs, fqg, temp1, delta2, a, b, c, d
1184
1185    do k = kts, kte
1186       do i = its, ite
1187
1188    !-------------------------------------------------------------
1189    ! pracs: Accretion of snow by rain [LFO 27]
1190    !         (T<T0: S->G) pracs: min=0., max=qrs(i,k,2)/dtcld
1191    !-------------------------------------------------------------
1192          supcol = t0c-t(i,k)
1193          call smoothif(supcol,0.,fsupcol,'t0')
1194          !call smoothif(qrs(i,k,1),1.e-4,fqr,'q0')
1195          !call smoothif(qrs(i,k,2),1.e-4,fqs,'q0')
1196
1197          vt2r=vt2r_a*(den(i,k)**((bvtr-2.)/4.))*(max(qrs(i,k,1),qcrmin)**(bvtr/4.))
1198          vt2s=vt2s_a*(den(i,k)**((bvts-2.)/4.))*(max(qrs(i,k,2),qcrmin)**(bvts/4.))&
1199                *exp(-alpha*bvts*max(0.,min(90.,t0c-t(i,k)))/4.)
1200          a=exp(    alpha*max(0.,min(90.,t0c-t(i,k)))   )
1201          b=exp(-3.*alpha*max(0.,min(90.,t0c-t(i,k)))/2.)&
1202             *(den(i,k)**(3./4.))*(max(qrs(i,k,2),qcrmin)**(3./2.))*sqrt(sqrt(max(qrs(i,k,1),qcrmin)))
1203          c=exp(-5.*alpha*max(0.,min(90.,t0c-t(i,k)))/4.)&
1204             *(den(i,k)**(3./4.))*(max(qrs(i,k,2),qcrmin)**(5./4.))*sqrt     (max(qrs(i,k,1),qcrmin))
1205          d=exp(-   alpha*max(0.,min(90.,t0c-t(i,k)))   )&
1206             *(den(i,k)**(3./4.))*(max(qrs(i,k,2),qcrmin))         *         (max(qrs(i,k,1),qcrmin)**(3./4.))
1207          pracs1=pracs_a*a*abs(vt2r-vt2s)*(pracs_b*b+pracs_c*c+pracs_d*d)
1208
1209          pracs(i,k) = min(pracs1,qrs(i,k,2)/dtcld)
1210          pracs(i,k)=fsupcol*pracs(i,k)
1211          if(abs(pracs(i,k))<qmin/dtcld)then
1212             pracs(i,k)=0.
1213          endif
1214          !update qs qg
1215          qrs(i,k,2)=max(qrs(i,k,2)-pracs(i,k)*dtcld,0.)
1216          qrs(i,k,3)=max(qrs(i,k,3)+pracs(i,k)*dtcld,0.)
1217          pracs(i,k)=0.
1218    !-------------------------------------------------------------
1219    ! psacr: Accretion of rain by snow [LFO 28]
1220    !         (T< T0: R->S or R->G)                 min=0.,max=qrs(i,k,1)/dtcld
1221    !         (T>=T0: S->R enhance melting of snow) min=0.,max=qrs(i,k,2)/dtcld
1222    !-------------------------------------------------------------
1223 !         supcol = t0c-t(i,k) !not change
1224 !         call smoothif(supcol,0.,fsupcol,'t0')
1225          !call smoothif(qrs(i,k,1),qcrmin,fqr,'q0')
1226          !call smoothif(qrs(i,k,2),qcrmin,fqs,'q0')
1227          !update cpm
1228          cpm(i,k)=cpmcal(q(i,k))
1229          xl(i,k)=xlcal(t(i,k))
1230          xlf = xls-xl(i,k)
1231          if(supcol.lt.0.) xlf = xlf0
1232
1233
1234          vt2r=vt2r_a*(den(i,k)**((bvtr-2.)/4.))*(max(qrs(i,k,1),qcrmin)**(bvtr/4.))
1235          vt2s=vt2s_a*(den(i,k)**((bvts-2.)/4.))*(max(qrs(i,k,2),qcrmin)**(bvts/4.))&
1236                *exp(-alpha*bvts*max(0.,min(90.,t0c-t(i,k)))/4.)
1237          a=exp(    alpha*max(0.,min(90.,t0c-t(i,k)))   )
1238          b=exp(-   alpha*max(0.,min(90.,t0c-t(i,k)))/4.)&
1239             *(den(i,k)**(3./4.))*(max(qrs(i,k,1),qcrmin)**(3./2.))*sqrt(sqrt(max(qrs(i,k,2),qcrmin)))
1240          c=exp(-   alpha*max(0.,min(90.,t0c-t(i,k)))/2.)&
1241             *(den(i,k)**(3./4.))*(max(qrs(i,k,1),qcrmin)**(5./4.))*     sqrt(max(qrs(i,k,2),qcrmin))
1242          d=exp(-3.*alpha*max(0.,min(90.,t0c-t(i,k)))/4.)&
1243             *(den(i,k)**(3./4.))*(max(qrs(i,k,1),qcrmin))         *         (max(qrs(i,k,2),qcrmin)**(3./4.))
1244          psacr1=psacr_a*a*abs(vt2r-vt2s)*(psacr_b*b+psacr_c*c+psacr_d*d)
1245
1246          if(supcol>0.)then ! T<T0
1247             psacr(i,k) = min(psacr1,qrs(i,k,1)/dtcld)
1248          else ! T>=T0
1249             psacr(i,k) = min(psacr1,qrs(i,k,2)/dtcld)
1250          endif
1251          !psacr(i,k)=fqr*fqs*psacr(i,k)
1252          if(abs(psacr(i,k))<qmin/dtcld)then
1253             psacr(i,k)=0.
1254          endif
1255          !update qr qs qg
1256          if(qrs(i,k,1).lt.1.e-4.and.qrs(i,k,2).lt.1.e-4)then
1257             delta2=1.
1258          else
1259             delta2=0.
1260          endif
1261          qrs(i,k,1)=max(qrs(i,k,1)-fsupcol           *psacr(i,k)*dtcld,0.)
1262          qrs(i,k,2)=max(qrs(i,k,2)+fsupcol*   delta2 *psacr(i,k)*dtcld,0.)
1263          qrs(i,k,3)=max(qrs(i,k,3)+fsupcol*(1-delta2)*psacr(i,k)*dtcld,0.)
1264          t  (i,k  )=t  (i,k  )+fsupcol*psacr(i,k)*dtcld*xlf/cpm(i,k)
1265          ! t>=t0 pseml
1266          psacr(i,k)=(1-fsupcol)*psacr(i,k)
1267
1268    !-------------------------------------------------------------
1269    ! pgacr: Accretion of rain by graupel [LFO 42]
1270    !         (T< T0: R->G)                            min=0.,max=qrs(i,k,1)/dtcld
1271    !         (T>=T0: G->R enhance melting of graupel) min=0.,max=qrs(i,k,3)/dtcld
1272    !-------------------------------------------------------------
1273          supcol = t0c-t(i,k)
1274          call smoothif(supcol,0.,fsupcol,'t0')
1275          !call smoothif(qrs(i,k,3),qcrmin,fqg,'q0')
1276          !call smoothif(qrs(i,k,1),qcrmin,fqr,'q0')
1277          !update cpm
1278 !         cpm(i,k)=cpmcal(q(i,k)) !not change
1279          xl(i,k)=xlcal(t(i,k))
1280          xlf = xls-xl(i,k)
1281          if(supcol.lt.0.) xlf = xlf0
1282
1283          vt2r=vt2r_a*(den(i,k)**((bvtr-2.)/4.))*(max(qrs(i,k,1),qcrmin)**(bvtr/4.))
1284          vt2g=vt2g_a*(den(i,k)**((bvtg-2.)/4.))*(max(qrs(i,k,3),qcrmin)**(bvtg/4.))
1285          b=(den(i,k)**(3./4.))*(max(qrs(i,k,1),qcrmin)**(3./2.))*sqrt(sqrt(max(qrs(i,k,3),qcrmin)))
1286          c=(den(i,k)**(3./4.))*(max(qrs(i,k,1),qcrmin)**(5./4.))*     sqrt(max(qrs(i,k,3),qcrmin))
1287          d=(den(i,k)**(3./4.))*(max(qrs(i,k,1),qcrmin))         *         (max(qrs(i,k,3),qcrmin)**(3./4.))
1288          pgacr1=pgacr_a*abs(vt2r-vt2g)*(pgacr_b*b+pgacr_c*c+pgacr_d*d)
1289
1290          if(supcol>0.)then ! T<T0
1291             pgacr(i,k) = min(pgacr1,qrs(i,k,1)/dtcld)
1292          else
1293             pgacr(i,k) = min(pgacr1,qrs(i,k,3)/dtcld)
1294          endif
1295          !pgacr(i,k)=fqg*fqr*pgacr(i,k)
1296          if(abs(pgacr(i,k))<qmin/dtcld)then
1297             pgacr(i,k)=0.
1298          endif
1299          ! t<t0 update qr,qg,t
1300          qrs(i,k,1)=max(qrs(i,k,1)-fsupcol*pgacr(i,k)*dtcld,0.)
1301          qrs(i,k,3)=max(qrs(i,k,3)+fsupcol*pgacr(i,k)*dtcld,0.)
1302          t  (i,k  )=t  (i,k  )+fsupcol*pgacr(i,k)*dtcld*xlf/cpm(i,k)
1303          ! t>=t0 pgeml
1304          pgacr(i,k)=(1-fsupcol)*pgacr(i,k)
1305    !-------------------------------------------------------------
1306    ! pgacs: Accretion of snow by graupel [LFO 29]
1307    !        (S->G) min=0,max=qrs(i,k,2)/dtcld
1308    !-------------------------------------------------------------
1309          supcol = t0c-t(i,k)
1310          call smoothif(supcol,0.,fsupcol,'t0')
1311          !call smoothif(qrs(i,k,3),qcrmin,fqg,'q0')
1312          !call smoothif(qrs(i,k,2),qcrmin,fqs,'q0')
1313          egs = min(exp(-0.09*supcol),1.)
1314
1315          vt2g=vt2g_a*(den(i,k)**((bvtg-2.)/4.))*(max(qrs(i,k,3),qcrmin)**(bvtg/4.))
1316          vt2s=vt2s_a*(den(i,k)**((bvts-2.)/4.))*(max(qrs(i,k,2),qcrmin)**(bvts/4.))&
1317                *exp(-alpha*bvts*max(0.,min(90.,t0c-t(i,k)))/4.)
1318          a=exp(    alpha*max(0.,min(90.,t0c-t(i,k)))   )
1319          b=exp(-3.*alpha*max(0.,min(90.,t0c-t(i,k)))/2.)*(den(i,k)**(3./4.))&
1320             *(max(qrs(i,k,2),qcrmin)**(3./2.))*sqrt(sqrt(max(qrs(i,k,3),qcrmin)))
1321          c=exp(-5.*alpha*max(0.,min(90.,t0c-t(i,k)))/4.)*(den(i,k)**(3./4.))&
1322             *(max(qrs(i,k,2),qcrmin)**(5./4.))*     sqrt(max(qrs(i,k,3),qcrmin))
1323          d=exp(-   alpha*max(0.,min(90.,t0c-t(i,k)))   )*(den(i,k)**(3./4.))&
1324             *(max(qrs(i,k,2),qcrmin))         *         (max(qrs(i,k,3),qcrmin)**(3./4.))
1325          pgacs1=pgacs_a*a*abs(vt2g-vt2s)*(pgacs_b*b+pgacs_c*c+pgacs_d*d)*egs
1326
1327          pgacs(i,k) = min(pgacs1,qrs(i,k,2)/dtcld)
1328          !pgacs(i,k)=fqg*fqs*pgacs(i,k)
1329          if(abs(pgacs(i,k))<qmin/dtcld)then
1330             pgacs(i,k)=0.
1331          endif
1332          !update qs, qg
1333          qrs(i,k,2)=max(qrs(i,k,2)-pgacs(i,k)*dtcld,0.)
1334          qrs(i,k,3)=max(qrs(i,k,3)+pgacs(i,k)*dtcld,0.)
1335          pgacs(i,k)=0.
1336
1337    !-------------------------------------------------------------
1338    ! pseml: Enhanced melting of snow by accretion of water
1339    !        (T>=T0: S->R) pseml<0 max=0,min=-qrs(i,k,2)/dtcld
1340    !-------------------------------------------------------------
1341 !         supcol = t0c-t(i,k) ! not change
1342          !update cpm
1343 !         cpm(i,k)=cpmcal(q(i,k)) ! not change
1344          xl(i,k)=xlcal(t(i,k))
1345          xlf = xls-xl(i,k)
1346          if(supcol.lt.0.) xlf = xlf0
1347          call smoothif(t(i,k),t0c,ft0,'t0')
1348          call smoothif(qrs(i,k,2),0.,fqs,'q+')
1349          pseml(i,k) = min(max(cliq*supcol*(psacw(i,k)+psacr(i,k))  &
1350                        /xlf,-qrs(i,k,2)/dtcld),0.)
1351          pseml(i,k)=ft0*fqs*pseml(i,k)
1352          if(abs(pseml(i,k))<qmin/dtcld)then
1353             pseml(i,k)=0.
1354          endif
1355          !update qs, qr, t
1356          qrs(i,k,1)=max(qrs(i,k,1)-pseml(i,k)*dtcld,0.)
1357          qrs(i,k,2)=max(qrs(i,k,2)+pseml(i,k)*dtcld,0.)
1358          t  (i,k  )=t  (i,k  )+pseml(i,k)*dtcld*xlf/cpm(i,k)
1359          pseml(i,k)=0.
1360          psacw(i,k)=0.
1361          psacr(i,k)=0.
1362    !-------------------------------------------------------------
1363    ! pgeml: Enhanced melting of graupel by accretion of water [RH84 A21-A22]
1364    !        (T>=T0: G->R) pgeml<0 max=0,min=-qrs(i,k,3)/dtcld
1365    !-------------------------------------------------------------
1366          supcol = t0c-t(i,k)
1367          !update cpm
1368 !         cpm(i,k)=cpmcal(q(i,k)) ! not change
1369          xl(i,k)=xlcal(t(i,k))
1370          xlf = xls-xl(i,k)
1371          if(supcol.lt.0.) xlf = xlf0
1372          call smoothif(t(i,k),t0c,ft0,'t0')
1373          call smoothif(qrs(i,k,3),0.,fqg,'q+')
1374          pgeml(i,k) = min(max(cliq*supcol*(pgacw(i,k)+pgacr(i,k))  &
1375                        /xlf,-qrs(i,k,3)/dtcld),0.)
1376          pgeml(i,k)=ft0*fqg*pgeml(i,k)
1377          if(abs(pgeml(i,k))<qmin/dtcld)then
1378             pgeml(i,k)=0.
1379          endif
1380          !update qg, qr, t
1381          qrs(i,k,1)=max(qrs(i,k,1)-pgeml(i,k)*dtcld,0.)
1382          qrs(i,k,3)=max(qrs(i,k,3)+pgeml(i,k)*dtcld,0.)
1383          t  (i,k  )=t  (i,k  )+pgeml(i,k)*dtcld*xlf/cpm(i,k)
1384          pgeml(i,k)=0.
1385          pgacw(i,k)=0.
1386          pgacr(i,k)=0.
1387       enddo
1388    enddo
1389    end subroutine accret2
1390
1391 !=======================================================================
1392 !
1393 !=======================================================================
1394    subroutine accret3(qrs,qci,rh,t,p,den,dtcld,  &
1395                      q,qs,psdep,pgdep,pigen,psaut,pgaut,psevp,pgevp,pidep,       &
1396                      ims,ime, kms,kme,its,ite,kts,kte)
1397    IMPLICIT NONE
1398    integer ::   ims,ime, kms,kme,its,ite, kts,kte
1399 !-------------------------------------------------------------------
1400    REAL, DIMENSION( its:ite , kts:kte, 2 ) ::       qci
1401    REAL, DIMENSION( ims:ime , kms:kme ) ::          den,q,p
1402    REAL, DIMENSION( its:ite , kts:kte , 3) :: qrs,rh,qs
1403    REAL, DIMENSION( its:ite , kts:kte) :: pigen,psevp,pgevp,pidep,t,xl,cpm,&
1404                                            psdep,pgdep,psaut,pgaut
1405    REAL :: supcol, dtcld,satdt,supsat,qimax, diameter,xni0,roqi0,supice1,supice2,supice3,supice4,alpha2
1406    real :: pidep0,pidep1,psdep0, pgdep3,pigen0,psevp0,pgevp0,coeres1,coeres2,coeres3,coeres4
1407    real :: temp0,temp,xmi
1408    INTEGER :: i,  k
1409    real :: fqi, fqr, fqv, fqs, fqg, frh, ft0, fpidep, fpsdep, fpgdep, fsupcol, fsupsat, pidep2
1410    real :: value01, factor01, source01, vice, a, b, c, d, e, f, g
1411
1412    do k = kts, kte
1413       do i = its, ite
1414    !
1415    !-------------------------------------------------------------
1416    ! pidep: Deposition/Sublimation rate of ice [HDC 9]
1417    !       (T<T0: V->I or I->V)
1418    !       rh(i,k,2)>1.,pidep>0: V->I, min=0,        max=satdt
1419    !       rh(i,k,2)<1.,pidep<0: I->V, min=-qi/dtcld,max=0,
1420    !-------------------------------------------------------------
1421          !update supcol
1422          supcol = t0c-t(i,k)
1423          !update rh qs
1424          call calcrh(t(i,k),p(i,k),q(i,k),rh(i,k,:),qs(i,k,:))
1425          !update satdt
1426          supsat = q(i,k)-qs(i,k,2)
1427          satdt = supsat/dtcld
1428          !update xl, cpm
1429          xl (i,k)=xlcal (t(i,k))
1430          cpm(i,k)=cpmcal(q(i,k))
1431          call smoothif(supcol,0.,fsupcol,'t+')
1432
1433          if(qci(i,k,2)>0.)then
1434             b=diffac_a*den(i,k)*xls*xls*(t(i,k)+120.)/rv/(t(i,k)**3.5)
1435             c=diffac_b*p(i,k)/(t(i,k)**1.81)/qs(i,k,2)
1436             a=(rh(i,k,2)-1.)/(b+c)
1437             pidep0 = pidep_a*a*((den(i,k)*qci(i,k,2))**(7./8.))
1438          else
1439             pidep0 = 0.
1440          endif
1441
1442          if(pidep0<0.)then
1443             pidep(i,k) = min(max(pidep0,-qci(i,k,2)/dtcld),0.) !qci(i,k,2)/dtcld*dtcld /= qci(i,k,2)
1444          else
1445             pidep(i,k) = max(min(pidep0,satdt),0.)
1446          endif
1447          pidep(i,k)=fsupcol*pidep(i,k)
1448          if(abs(pidep(i,k))<qmin/dtcld)then
1449             pidep(i,k)=0.
1450          endif
1451          !update qv, qi, t
1452          q  (i,k  )=max(q  (i,k  )-pidep(i,k)*dtcld,0.)
1453          qci(i,k,2)=max(qci(i,k,2)+pidep(i,k)*dtcld,0.)
1454          t  (i,k  )=t  (i,k  )+pidep(i,k)*dtcld*xls/cpm(i,k)
1455          pidep(i,k)=0.
1456    !
1457    !-------------------------------------------------------------
1458    ! psdep: deposition/sublimation rate of snow [HDC 14]
1459    !        (T<T0: V->S or S->V)
1460    !       rh(i,k,2)>1.,psdep>0: V->S, min=0,        max=satdt
1461    !       rh(i,k,2)<1.,psdep<0: S->V, min=-qs/dtcld,max=0,
1462    !-------------------------------------------------------------
1463          !update supcol
1464          supcol = t0c-t(i,k)
1465          !update rh qs
1466          call calcrh(t(i,k),p(i,k),q(i,k),rh(i,k,:),qs(i,k,:))
1467          !update satdt
1468          supsat = q(i,k)-qs(i,k,2)
1469          satdt = supsat/dtcld
1470          !update xl, cpm
1471          xl (i,k)=xlcal (t(i,k))
1472          cpm(i,k)=cpmcal(q(i,k))
1473
1474          call smoothif(supcol,0.,fsupcol,'t+')
1475 !         call smoothif(qrs(i,k,2),0.,fqs,'q+')
1476 !         call smoothif(q  (i,k  ),0.,fqv,'q+')
1477
1478          a=exp(alpha*max(0.,min(90.,t0c-t(i,k)))/2.)*sqrt(den(i,k)*max(qrs(i,k,2),qcrmin))
1479          b=exp((3.-bvts)*alpha*max(0.,min(90.,t0c-t(i,k)))/8.)*((t(i,k)+120.)**(1./6.))/(t(i,k)**(5.12/6.)) &
1480             *(p(i,k)**(1./3.))*(den(i,k)**((13.+3.*bvts)/24.))*(max(qrs(i,k,2),qcrmin)**((5.+bvts)/8.))
1481          c=diffac_a*den(i,k)*xls*xls*(t(i,k)+120.)/rv/(t(i,k)**3.5)
1482          d=diffac_b*p(i,k)/(t(i,k)**1.81)/qs(i,k,2)
1483          e=(rh(i,k,2)-1.)/(c+d)
1484          psdep0 = e*(psdep_a*a+psdep_b*b)
1485
1486          if(psdep0<0.)then
1487             psdep(i,k) = min(max(psdep0,-qrs(i,k,2)/dtcld),0.)
1488          else
1489             psdep(i,k) = max(min(psdep0,satdt),0.)
1490          endif
1491          psdep(i,k)=fsupcol*psdep(i,k)
1492          if(abs(psdep(i,k))<qmin/dtcld)then
1493             psdep(i,k)=0.
1494          endif
1495          ! update q, qs, t
1496          q  (i,k  )=max(q  (i,k  )-psdep(i,k)*dtcld,0.)
1497          qrs(i,k,2)=max(qrs(i,k,2)+psdep(i,k)*dtcld,0.)
1498          t  (i,k  )=t  (i,k  )+psdep(i,k)*dtcld*xls/cpm(i,k)
1499          psdep(i,k)=0.
1500    !
1501    !------------------------------------------------------------
1502    ! pgdep: deposition/sublimation rate of graupel [LFO 46]
1503    !        (T<T0: V->G or G->V)
1504    !       rh(i,k,2)>1.,pgdep>0: V->G, min=0,        max=satdt
1505    !       rh(i,k,2)<1.,pgdep<0: G->V, min=-qg/dtcld,max=0,
1506    !------------------------------------------------------------
1507          !update supcol
1508          supcol = t0c-t(i,k)
1509          !update rh qs
1510          call calcrh(t(i,k),p(i,k),q(i,k),rh(i,k,:),qs(i,k,:))
1511          !update satdt
1512          supsat = q(i,k)-qs(i,k,2)
1513          satdt = supsat/dtcld
1514          !update xl, cpm
1515          xl (i,k)=xlcal (t(i,k))
1516          cpm(i,k)=cpmcal(q(i,k))
1517
1518          call smoothif(supcol,0.,fsupcol,'t+')
1519 !         call smoothif(qrs(i,k,3),0.,fqg,'q+')
1520 !         call smoothif(q  (i,k  ),0.,fqv,'q+')
1521
1522          a=sqrt(den(i,k)*max(qrs(i,k,3),qcrmin))
1523          b=((t(i,k)+120.)**(1./6.))/(t(i,k)**(5.12/6.))*(p(i,k)**(1./3.)) &
1524             *(den(i,k)**((13.+3.*bvtg)/24.))*(max(qrs(i,k,3),qcrmin)**((5.+bvtg)/8.))
1525          c=diffac_a*den(i,k)*xls*xls*(t(i,k)+120.)/rv/(t(i,k)**3.5)
1526          d=diffac_b*p(i,k)/(t(i,k)**1.81)/qs(i,k,2)
1527          e=(rh(i,k,2)-1.)/(c+d)
1528          pgdep3 = e*(pgdep_a*a+pgdep_b*b)
1529
1530          if(pgdep3<0.)then
1531             pgdep(i,k) = min(max(pgdep3,-qrs(i,k,3)/dtcld),0.)
1532          else
1533             pgdep(i,k) = max(min(pgdep3,satdt),0.) ! V->G !
1534          endif
1535          pgdep(i,k)=fsupcol*pgdep(i,k)
1536          if(abs(pgdep(i,k))<qmin/dtcld)then
1537             pgdep(i,k)=0.
1538          endif
1539          !update q, qg, t
1540          q  (i,k  )=max(q  (i,k  )-pgdep(i,k)*dtcld,0.)
1541          qrs(i,k,3)=max(qrs(i,k,3)+pgdep(i,k)*dtcld,0.)
1542          t  (i,k  )=t  (i,k  )+pgdep(i,k)*dtcld*xls/cpm(i,k)
1543          pgdep(i,k)=0.
1544    !-------------------------------------------------------------
1545    ! pigen: generation(nucleation) of ice from vapor [HDC 7-8]
1546    !       (T<T0: V->I) min=0,max=min(roqi0/den(i,k)-max(qci(i,k,2),0.))/dtcld,satdt)
1547    !-------------------------------------------------------------
1548          !update supcol
1549          supcol = t0c-t(i,k)
1550          cpm(i,k)=cpmcal(q(i,k))
1551          !update rh qs
1552          call calcrh(t(i,k),p(i,k),q(i,k),rh(i,k,:),qs(i,k,:))
1553          !update satdt
1554          supsat = q(i,k)-qs(i,k,2)
1555          satdt = supsat/dtcld
1556          call smoothif(supsat,0.,fsupsat,'q+')
1557          call smoothif(supcol,0.,fsupcol,'t+')
1558
1559          xni0 = 1.e3*exp(0.1*supcol)
1560          roqi0 = 4.92e-11*xni0**1.33
1561          pigen0 = min((roqi0/den(i,k)-max(qci(i,k,2),0.))/dtcld,satdt)
1562          pigen(i,k)=max(pigen0,0.) !
1563          pigen(i,k)=fsupcol*fsupsat*pigen(i,k)
1564          if(abs(pigen(i,k))<qmin/dtcld)then
1565             pigen(i,k)=0.
1566          endif
1567          !update q, qi, t
1568          q  (i,k  )=max(q  (i,k  )-pigen(i,k)*dtcld,0.)
1569          qci(i,k,2)=max(qci(i,k,2)+pigen(i,k)*dtcld,0.)
1570          t  (i,k  )=t  (i,k  )+pigen(i,k)*dtcld*xls/cpm(i,k)
1571          pigen(i,k )=0.
1572    !
1573    !------------------------------------------------------------
1574    ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1575    !        (T<T0: I->S) psaut>0, min=0,max=(qci(i,k,2)-qimax)/dtcld
1576    !-------------------------------------------------------------
1577          !update supcol
1578          supcol = t0c-t(i,k)
1579          call smoothif(supcol,0.,fsupcol,'t+')
1580 !         call smoothif(qci(i,k,2),0.,fqi,'q+')
1581          qimax = roqimax/den(i,k)
1582          psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1583          psaut(i,k) = fsupcol*psaut(i,k)
1584          if(abs(psaut(i,k))<qmin/dtcld)then
1585             psaut(i,k)=0.
1586          endif
1587          !update qi, qs
1588          qci(i,k,2)=max(qci(i,k,2)-psaut(i,k)*dtcld,0.)
1589          qrs(i,k,2)=max(qrs(i,k,2)+psaut(i,k)*dtcld,0.)
1590          psaut(i,k)=0.
1591    !
1592    !-------------------------------------------------------------
1593    ! pgaut: conversion(aggregation) of snow to graupel [LFO 37]
1594    !        (T<T0: S->G) pgaut>0 min=0.,max=qrs(i,k,2)/dtcld
1595    !-------------------------------------------------------------
1596          !update supcol
1597 !         supcol = t0c-t(i,k) ! not change
1598 !         call smoothif(supcol,0.,fsupcol,'t0')
1599 !         call smoothif(qrs(i,k,2),0.,fqs,'q+')
1600          alpha2 = 1.e-3*exp(0.09*(-supcol))
1601          pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0))         &
1602                       ,qrs(i,k,2)/dtcld)
1603          pgaut(i,k)=fsupcol*pgaut(i,k)
1604          if(abs(pgaut(i,k))<qmin/dtcld)then
1605             pgaut(i,k)=0.
1606          endif
1607          !update qs, qg
1608          qrs(i,k,2)=max(qrs(i,k,2)-pgaut(i,k)*dtcld,0.)
1609          qrs(i,k,3)=max(qrs(i,k,3)+pgaut(i,k)*dtcld,0.)
1610          pgaut(i,k)=0.
1611    !
1612    !-------------------------------------------------------------
1613    ! psevp: Evaporation of melting snow [RH83 A27]
1614    !       (T>=T0: S->V) rh<1., psevp<0, min=-qrs(i,k,2)/dtcld, max=0.
1615    !-------------------------------------------------------------
1616 !         supcol = t0c-t(i,k) ! not change
1617          !update rh qs
1618          call calcrh(t(i,k),p(i,k),q(i,k),rh(i,k,:),qs(i,k,:))
1619          !update xl, cpm
1620          xl (i,k)=xlcal (t(i,k))
1621          cpm(i,k)=cpmcal(q(i,k))
1622          call smoothif(t(i,k),    t0c,ft0,'t+')
1623
1624          a=exp(alpha*max(0.,min(90.,t0c-t(i,k)))/2.)*sqrt(den(i,k)*max(qrs(i,k,2),qcrmin))
1625          b=exp((3.-bvts)*alpha*max(0.,min(90.,t0c-t(i,k)))/8.)*((t(i,k)+120.)**(1./6.)) &
1626             /(t(i,k)**(5.12/6.))*(p(i,k)**(1./3.))*(den(i,k)**((13.+3.*bvts)/24.))*(max(qrs(i,k,2),qcrmin)**((5.+bvts)/8.))
1627          c=diffac_a*den(i,k)*xl(i,k)*xl(i,k)*(t(i,k)+120.)/rv/(t(i,k)**3.5)
1628          d=diffac_b*p(i,k)/(t(i,k)**1.81)/qs(i,k,1)
1629          e=(rh(i,k,1)-1.)/(c+d)
1630          psevp0 = e*(psevp_a*a+psevp_b*b)
1631
1632          psevp(i,k) = min(max(psevp0,-qrs(i,k,2)/dtcld),0.)
1633          psevp(i,k) = ft0*psevp(i,k)
1634          if(abs(psevp(i,k))<qmin/dtcld)then
1635             psevp(i,k)=0.
1636          endif
1637          !update q, qs, t
1638          q  (i,k  )=max(q  (i,k  )-psevp(i,k)*dtcld,0.)
1639          qrs(i,k,2)=max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1640          t  (i,k  )=t  (i,k  )+psevp(i,k)*dtcld*xls/cpm(i,k)
1641          psevp(i,k )=0.
1642    !
1643    !-------------------------------------------------------------
1644    ! pgevp: Evaporation of melting graupel [RH84 A19]
1645    !       (T>=T0: G->V) rh<1., pgevp<0, min=-qrs(i,k,3)/dtcld, max=0.
1646    !-------------------------------------------------------------
1647          supcol = t0c-t(i,k)
1648          !update rh qs
1649          call calcrh(t(i,k),p(i,k),q(i,k),rh(i,k,:),qs(i,k,:))
1650          !update xl, cpm
1651          xl (i,k)=xlcal (t(i,k))
1652          cpm(i,k)=cpmcal(q(i,k))
1653
1654          call smoothif(t(i,k),    t0c,ft0,'t0')
1655
1656          a=sqrt(den(i,k)*max(qrs(i,k,3),qcrmin))
1657          b=((t(i,k)+120.)**(1./6.))/(t(i,k)**(5.12/6.))*(p(i,k)**(1./3.)) &
1658             *(den(i,k)**((13.+3.*bvtg)/24.))*(max(qrs(i,k,3),qcrmin)**((5.+bvtg)/8.))
1659          c=diffac_a*den(i,k)*xl(i,k)*xl(i,k)*(t(i,k)+120.)/rv/(t(i,k)**3.5)
1660          d=diffac_b*p(i,k)/(t(i,k)**1.81)/qs(i,k,1)
1661          e=(rh(i,k,1)-1.)/(c+d)
1662          pgevp0 = e*(pgevp_a*a+pgevp_b*b)
1663
1664          pgevp(i,k) = min(max(pgevp0,-qrs(i,k,3)/dtcld),0.)
1665          pgevp(i,k) = ft0*pgevp(i,k)
1666          if(abs(pgevp(i,k))<qmin/dtcld)then
1667             pgevp(i,k)=0.
1668          endif
1669          !update q, qg, t
1670          q  (i,k  )=max(q  (i,k  )-pgevp(i,k)*dtcld,0.)
1671          qrs(i,k,3)=max(qrs(i,k,3)+pgevp(i,k)*dtcld,0.)
1672          t  (i,k  )=t  (i,k  )+pgevp(i,k)*dtcld*xls/cpm(i,k)
1673          pgevp(i,k )=0.
1674       enddo
1675    enddo
1676    end subroutine accret3
1677
1678
1679 !=======================================================================
1680 !
1681 !=======================================================================
1682    subroutine pconadd(t,p,q,qci,qs,xl,cpm,dtcld,kte, kts,its, ite,kme, kms,ims, ime)
1683    IMPLICIT NONE
1684    integer :: ims,ime, kms,kme ,its,ite, kts,kte
1685    REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci
1686    REAL, DIMENSION( its:ite , kts:kte ) :: t,xl,pcond,work2,cpm
1687    REAL, DIMENSION( its:ite , kts:kte, 3 ) :: qs,work1,rh
1688    REAL, DIMENSION( ims:ime , kms:kme ) ::q, p
1689    integer :: k,i
1690    real :: hsub,hvap,cvap,ttp,dldt,xa,xb,dldti,xai,xbi,tr,dtcld,qs1,qs2,qs3,qs4,w1,q1
1691    real :: tmp1, tmp2, f1, f2,qs0
1692
1693    do k = kts, kte
1694       do i = its, ite
1695          !update qs
1696          call calcrh(t(i,k),p(i,k),q(i,k),rh(i,k,:),qs(i,k,:))
1697          !update xl, cpm
1698          xl (i,k)=xlcal (t(i,k))
1699          cpm(i,k)=cpmcal(q(i,k))
1700    !----------------------------------------------------------------
1701    !  pcond: condensational/evaporational rate of cloud water [RH83 A6]
1702    !     if there exists additional water vapor condensated/if
1703    !     evaporation of cloud water is not enough to remove subsaturation
1704    !q>qs, work1>0, pcond>0 V->C min=0,                 max=q(i,k)/dtcld
1705    !q<qs, work1<0, pcond<0 C->V min=-qci(i,k,1)/dtcld, max=0,
1706          work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
1707
1708          if(work1(i,k,1)>0.)then
1709             pcond(i,k)=min(work1(i,k,1),max(q(i,k),0.))/dtcld
1710          else
1711             pcond(i,k)=max(work1(i,k,1),-qci(i,k,1))/dtcld
1712          endif
1713
1714          if(abs(pcond(i,k))<qmin/dtcld)then
1715             pcond(i,k)=0.
1716          endif
1717          q(i,k)    = max(q(i,k)    -pcond(i,k)*dtcld,0.)
1718          qci(i,k,1)= max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1719          t(i,k)    = t(i,k)    +pcond(i,k)*dtcld*xl(i,k)/cpm(i,k)
1720          pcond(i,k)=0.
1721       enddo
1722    enddo
1723
1724    end subroutine pconadd
1725
1726 !=======================================================================
1727 !
1728 !=======================================================================
1729    subroutine smoothif(x,a,f,opt)
1730    implicit none
1731    real, intent(in)  :: x, a
1732    character(len=2), intent(in) :: opt
1733    real, intent(out) :: f
1734
1735    real(kind=8) :: k1, a1, x1, c1, f1, k, b
1736
1737    x1=x
1738    a1=a
1739    if(opt(1:1)=='q')then
1740      c1=1.E-15
1741    else
1742      c1=1.E-9
1743    endif
1744    k1=747./c1
1745    if(opt(2:2)=="+")then
1746      b=a1+710./k1
1747    else
1748      b=a1
1749    endif
1750    k=-k1*(x1-b)
1751    f1=1./(1.+exp(k))
1752    f=f1
1753    end subroutine
1754
1755 !
1756 !
1757 !=======================================================================
1758 !
1759 !=======================================================================
1760    REAL FUNCTION rgmma(x)
1761 !-------------------------------------------------------------------
1762    IMPLICIT NONE
1763 !-------------------------------------------------------------------
1764 !  rgmma function:  use infinite product form
1765    REAL :: euler
1766    PARAMETER (euler=0.577215664901532)
1767    REAL :: x, y
1768    INTEGER :: i
1769    if(x.eq.1.)then
1770       rgmma=0.
1771    else
1772       rgmma=x*exp(euler*x)
1773       do i=1,10000
1774          y=float(i)
1775          rgmma=rgmma*(1.000+x/y)*exp(-x/y)
1776       enddo
1777       rgmma=1./rgmma
1778    endif
1779    END FUNCTION rgmma
1780 !
1781 !
1782 !=======================================================================
1783 !
1784 !=======================================================================
1785 !   compute internal functions
1786    function cpmcal(x)
1787    implicit none
1788    real:: cpmcal,x
1789    cpmcal= cpd+x*(cpv-cpd)
1790    end function
1791 !
1792 !=======================================================================
1793 !
1794 !=======================================================================
1795    function xlcal(x)
1796    implicit none
1797    real:: xlcal,x
1798    xlcal= xlv0-xlv1*(x-t0c)
1799    end function
1800
1801
1802 !=======================================================================
1803 ! a:t, b:q, c:qs, d:xl, e:cpm
1804 !=======================================================================
1805    function conden(a,b,c,d,e)
1806    implicit none
1807    real :: conden,a,b,c,d,e
1808    real :: f
1809    conden = (b-c)/(1.+d*d/(rv*e)*c/(a*a))
1810    end function
1811 END MODULE module_mp_wsm6r
1812 !