phys/module_mp_lin.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2 !
   3
   4 MODULE module_mp_lin
   5
   6    USE     module_wrf_error
   7    USE module_mp_radar
   8 !
   9    REAL    , PARAMETER, PRIVATE ::       RH = 1.0
  10 !  REAL    , PARAMETER, PRIVATE ::   episp0 = 0.622*611.21
  11    REAL    , PARAMETER, PRIVATE ::     xnor = 8.0e6
  12    REAL    , PARAMETER, PRIVATE ::     xnos = 3.0e6
  13
  14 !  Lin
  15 !  REAL    , PARAMETER, PRIVATE ::     xnog = 4.0e4
  16 !  REAL    , PARAMETER, PRIVATE ::     rhograul = 917.
  17
  18 ! Hobbs
  19   REAL     , PARAMETER, PRIVATE ::     xnog = 4.0e6
  20   REAL     , PARAMETER, PRIVATE ::     rhograul = 400.
  21
  22 !
  23   REAL     , PARAMETER, PRIVATE ::                              &
  24              qi0 = 1.0e-3, ql0 = 7.0e-4, qs0 = 6.0E-4,          &
  25              xmi50 = 4.8e-10, xmi40 = 2.46e-10,                 &
  26              constb = 0.8, constd = 0.25,                       &
  27              o6 = 1./6.,  cdrag = 0.6,                          &
  28              avisc = 1.49628e-6, adiffwv = 8.7602e-5,           &
  29              axka = 1.4132e3, di50 = 1.0e-4, xmi = 4.19e-13,    &
  30              cw = 4.187e3, vf1s = 0.78, vf2s = 0.31,            &
  31              xni0 = 1.0e-2, xmnin = 1.05e-18, bni = 0.5,        &
  32              ci = 2.093e3
  33 CONTAINS
  34
  35 !-------------------------------------------------------------------
  36 !  Lin et al., 1983, JAM, 1065-1092, and
  37 !  Rutledge and Hobbs, 1984, JAS, 2949-2972
  38 !  May 2009 - Changes and corrections from P. Blossey (U. Washington)
  39 !-------------------------------------------------------------------
  40   SUBROUTINE lin_et_al(th                                          &
  41                       ,qv, ql, qr                                  &
  42                       ,qi, qs                                      &
  43                       ,rho, pii, p                                 &
  44                       ,dt_in                                       &
  45                       ,z,ht, dz8w                                  &
  46                       ,grav, cp, Rair, rvapor                      &
  47                       ,XLS, XLV, XLF, rhowater, rhosnow            &
  48                       ,EP2,SVP1,SVP2,SVP3,SVPT0                    &
  49                       , RAINNC, RAINNCV                            &
  50                       , SNOWNC, SNOWNCV                            &
  51                       , GRAUPELNC, GRAUPELNCV, SR                  &
  52                       ,refl_10cm, diagflag, do_radar_ref           &
  53                       ,ids,ide, jds,jde, kds,kde                   &
  54                       ,ims,ime, jms,jme, kms,kme                   &
  55                       ,its,ite, jts,jte, kts,kte                   &
  56                   ! Optional
  57                       ,qlsink, precr, preci, precs, precg          &
  58                       , F_QG,F_QNDROP                              &
  59                       , qg, qndrop                                 &
  60                                                                    )
  61 !-------------------------------------------------------------------
  62   IMPLICIT NONE
  63 !-------------------------------------------------------------------
  64 !
  65 ! Shuhua 12/17/99
  66 !
  67   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
  68                                       ims,ime, jms,jme, kms,kme , &
  69                                       its,ite, jts,jte, kts,kte
  70
  71   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
  72         INTENT(INOUT) ::                                          &
  73                                                               th, &
  74                                                               qv, &
  75                                                               ql, &
  76                                                               qr
  77
  78 !
  79   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
  80         INTENT(IN   ) ::                                          &
  81                                                              rho, &
  82                                                              pii, &
  83                                                                p, &
  84                                                             dz8w
  85
  86   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
  87         INTENT(IN   ) ::                                       z
  88
  89
  90
  91   REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) ::       ht
  92
  93   REAL, INTENT(IN   ) ::                                   dt_in, &
  94                                                             grav, &
  95                                                             Rair, &
  96                                                           rvapor, &
  97                                                               cp, &
  98                                                              XLS, &
  99                                                              XLV, &
 100                                                              XLF, &
 101                                                         rhowater, &
 102                                                          rhosnow
 103
 104   REAL, INTENT(IN   ) :: EP2,SVP1,SVP2,SVP3,SVPT0
 105
 106   REAL, DIMENSION( ims:ime , jms:jme ),                           &
 107         INTENT(INOUT) ::                                  RAINNC, &
 108                                                          RAINNCV, &
 109                                                               SR
 110
 111 !+---+-----------------------------------------------------------------+
 112   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT)::     &  ! GT
 113                                                        refl_10cm
 114 !+---+-----------------------------------------------------------------+
 115
 116 ! Optional
 117   REAL, DIMENSION( ims:ime , jms:jme ),                           &
 118         OPTIONAL,                                                 &
 119         INTENT(INOUT) ::                                  SNOWNC, &
 120                                                          SNOWNCV, &
 121                                                        GRAUPELNC, &
 122                                                       GRAUPELNCV
 123
 124   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
 125         OPTIONAL,                                                 &
 126         INTENT(INOUT) ::                                          &
 127                                                               qi, &
 128                                                               qs, &
 129                                                               qg, &
 130                                                           qndrop
 131
 132   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
 133         OPTIONAL, INTENT(OUT   ) ::                               &
 134         qlsink, & ! cloud water conversion to rain (/s)
 135         precr,  & ! rain precipitation rate at all levels (kg/m2/s)
 136         preci,  & ! ice precipitation rate at all levels (kg/m2/s)
 137         precs,  & ! snow precipitation rate at all levels (kg/m2/s)
 138         precg     ! graupel precipitation rate at all levels (kg/m2/s)
 139
 140   LOGICAL, INTENT(IN), OPTIONAL ::                F_QG, F_QNDROP
 141
 142 ! LOCAL VAR
 143
 144   INTEGER             ::                            min_q, max_q
 145
 146   REAL, DIMENSION( its:ite , jts:jte )                            &
 147                                ::        rain, snow, graupel,ice
 148
 149   REAL, DIMENSION( kts:kte )   ::                  qvz, qlz, qrz, &
 150                                                    qiz, qsz, qgz, &
 151                                                              thz, &
 152                                                      tothz, rhoz, &
 153                                                    orhoz, sqrhoz, &
 154                                                         prez, zz, &
 155                                   precrz, preciz, precsz, precgz, &
 156                                                          qndropz, &
 157                                                      dzw, preclw
 158
 159   LOGICAL :: flag_qg, flag_qndrop
 160 !
 161   REAL    ::         dt, pptrain, pptsnow, pptgraul, rhoe_s,      &
 162                      gindex, pptice
 163   real :: qndropconst
 164
 165   INTEGER ::               i,j,k
 166
 167 !+---+-----------------------------------------------------------------+
 168       REAL, DIMENSION(kts:kte):: qv1d, t1d, p1d, qr1d, qs1d, qg1d, dBZ
 169       LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
 170       INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
 171 !+---+-----------------------------------------------------------------+
 172
 173    flag_qg     = .false.
 174    flag_qndrop = .false.
 175    IF ( PRESENT ( f_qg     ) ) flag_qg     = f_qg
 176    IF ( PRESENT ( f_qndrop ) ) flag_qndrop = f_qndrop
 177 !
 178    dt=dt_in
 179    rhoe_s=1.29
 180    qndropconst=100.e6  !sg
 181    gindex=1.0
 182
 183    IF (.not.flag_qg) gindex=0.
 184
 185    j_loop:  DO j = jts, jte
 186    i_loop:  DO i = its, ite
 187 !
 188 !- write data from 3-D to 1-D
 189 !
 190    DO k = kts, kte
 191       qvz(k)=qv(i,k,j)
 192       qlz(k)=ql(i,k,j)
 193       qrz(k)=qr(i,k,j)
 194       thz(k)=th(i,k,j)
 195 !
 196       rhoz(k)=rho(i,k,j)
 197       orhoz(k)=1./rhoz(k)
 198       prez(k)=p(i,k,j)
 199       sqrhoz(k)=sqrt(rhoe_s*orhoz(k))
 200       tothz(k)=pii(i,k,j)
 201       zz(k)=z(i,k,j)
 202       dzw(k)=dz8w(i,k,j)
 203    END DO
 204
 205    IF (flag_qndrop .AND. PRESENT( qndrop )) THEN
 206      DO k = kts, kte
 207          qndropz(k)=qndrop(i,k,j)
 208       ENDDO
 209    ELSE
 210       DO k = kts, kte
 211          qndropz(k)=qndropconst
 212       ENDDO
 213    ENDIF
 214
 215    DO k = kts, kte
 216       qiz(k)=qi(i,k,j)
 217       qsz(k)=qs(i,k,j)
 218    ENDDO
 219
 220    IF ( flag_qg .AND. PRESENT( qg ) ) THEN
 221       DO k = kts, kte
 222          qgz(k)=qg(i,k,j)
 223       ENDDO
 224    ELSE
 225       DO k = kts, kte
 226          qgz(k)=0.
 227       ENDDO
 228    ENDIF
 229 !
 230    pptrain=0.
 231    pptsnow=0.
 232    pptgraul=0.
 233    pptice=0.
 234    CALL clphy1d(    dt, qvz, qlz, qrz, qiz, qsz, qgz,         &
 235                     qndropz,flag_qndrop,                      &
 236                     thz, tothz, rhoz, orhoz, sqrhoz,          &
 237                     prez, zz, dzw, ht(I,J), preclw,           &
 238                     precrz, preciz, precsz, precgz,           &
 239                     pptrain, pptsnow, pptgraul, pptice,       &
 240                     grav, cp, Rair, rvapor, gindex,           &
 241                     XLS, XLV, XLF, rhowater, rhosnow,         &
 242                     EP2,SVP1,SVP2,SVP3,SVPT0,                 &
 243                     kts, kte, i, j                            )
 244
 245 !
 246 ! Precipitation from cloud microphysics -- only for one time step
 247 !
 248 ! unit is transferred from m to mm
 249
 250 !
 251    rain(i,j)=pptrain
 252    snow(i,j)=pptsnow
 253    graupel(i,j)=pptgraul
 254    ice(i,j)=pptice
 255    sr(i,j)=(pptice+pptsnow+pptgraul)/(pptice+pptsnow+pptgraul+pptrain+1.e-12)
 256 !
 257    RAINNCV(i,j)= pptrain + pptsnow + pptgraul + pptice
 258    RAINNC(i,j)=RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
 259    IF(PRESENT(SNOWNCV))SNOWNCV(i,j)= pptsnow + pptice
 260    IF(PRESENT(SNOWNC))SNOWNC(i,j)=SNOWNC(i,j) + pptsnow + pptice
 261    IF(PRESENT(GRAUPELNCV))GRAUPELNCV(i,j)= pptgraul
 262    IF(PRESENT(GRAUPELNC))GRAUPELNC(i,j)=GRAUPELNC(i,j) + pptgraul
 263
 264 !
 265 !- update data from 1-D back to 3-D
 266 !
 267 !
 268    IF ( present(qlsink) .and. present(precr) ) THEN !sg beg
 269       DO k = kts, kte
 270          if(ql(i,k,j)>1.e-20) then
 271             qlsink(i,k,j)=-preclw(k)/ql(i,k,j)
 272          else
 273             qlsink(i,k,j)=0.
 274          endif
 275          precr(i,k,j)=precrz(k)
 276       END DO
 277    END IF                                          !sg end
 278
 279    DO k = kts, kte
 280       qv(i,k,j)=qvz(k)
 281       ql(i,k,j)=qlz(k)
 282       qr(i,k,j)=qrz(k)
 283       th(i,k,j)=thz(k)
 284    END DO
 285 !
 286    IF ( flag_qndrop .AND. PRESENT( qndrop ) ) THEN !sg beg
 287       DO k = kts, kte
 288          qndrop(i,k,j)=qndropz(k)
 289       ENDDO
 290    END IF                                          !sg end
 291
 292    DO k = kts, kte
 293       qi(i,k,j)=qiz(k)
 294       qs(i,k,j)=qsz(k)
 295    ENDDO
 296
 297    IF ( present(preci) ) THEN     !sg beg
 298       DO k = kts, kte
 299          preci(i,k,j)=preciz(k)
 300       ENDDO
 301    END IF
 302
 303    IF ( present(precs) ) THEN
 304       DO k = kts, kte
 305          precs(i,k,j)=precsz(k)
 306       ENDDO
 307    END IF                         !sg end
 308
 309    IF ( flag_qg .AND. PRESENT( qg ) ) THEN
 310       DO k = kts, kte
 311          qg(i,k,j)=qgz(k)
 312       ENDDO
 313
 314       IF ( present(precg) ) THEN  !sg beg
 315          DO k = kts, kte
 316             precg(i,k,j)=precgz(k)
 317          ENDDO                    !sg end
 318       END IF
 319    ELSE                           !sg beg
 320       IF ( present(precg) ) precg(i,:,j)=0.  !sg end
 321    ENDIF
 322
 323 !+---+-----------------------------------------------------------------+
 324          IF ( PRESENT (diagflag) ) THEN
 325          if (diagflag .and. do_radar_ref == 1) then
 326                DO K=kts,kte
 327                   t1d(k)=th(i,k,j)*pii(i,k,j)
 328                   p1d(k)=p(i,k,j)
 329                   qv1d(k)=qv(i,k,j)
 330                   qr1d(k)=qr(i,k,j)
 331                   qs1d(k)=qs(i,k,j)
 332                   qg1d(k)=qg(i,k,j)
 333                ENDDO
 334                call refl10cm_lin (qv1d, qr1d, qs1d, qg1d,              &
 335                        t1d, p1d, dBZ, kts, kte, i, j)
 336                do k = kts, kte
 337                   refl_10cm(i,k,j) = MAX(-35., dBZ(k))
 338                enddo
 339          endif
 340          ENDIF
 341 !+---+-----------------------------------------------------------------+
 342
 343 !
 344    ENDDO i_loop
 345    ENDDO j_loop
 346
 347    END SUBROUTINE lin_et_al
 348
 349
 350 !-----------------------------------------------------------------------
 351    SUBROUTINE clphy1d(dt, qvz, qlz, qrz, qiz, qsz, qgz,                &
 352                       qndropz,flag_qndrop,                             &
 353                       thz, tothz, rho, orho, sqrho,                    &
 354                       prez, zz, dzw, zsfc, preclw,                     &
 355                       precrz, preciz, precsz, precgz,                  &
 356                       pptrain, pptsnow, pptgraul,                      &
 357                       pptice, grav, cp, Rair, rvapor, gindex,          &
 358                       XLS, XLV, XLF, rhowater, rhosnow,                &
 359                       EP2,SVP1,SVP2,SVP3,SVPT0,                        &
 360                       kts, kte, i, j                                   )
 361 !-----------------------------------------------------------------------
 362     IMPLICIT NONE
 363 !-----------------------------------------------------------------------
 364 !  This program handles the vertical 1-D cloud micphysics
 365 !-----------------------------------------------------------------------
 366 ! avisc: constant in empirical formular for dynamic viscosity of air
 367 !         =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
 368 ! adiffwv: constant in empirical formular for diffusivity of water
 369 !          vapor in air
 370 !          = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
 371 ! axka: constant in empirical formular for thermal conductivity of air
 372 !       = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
 373 ! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
 374 ! xmi50: mass of a 50 micron ice crystal
 375 !        = 4.8e-10 [kg] =4.8e-7 [g]
 376 ! xmi40: mass of a 40 micron ice crystal
 377 !        = 2.46e-10 [kg] = 2.46e-7 [g]
 378 ! di50: diameter of a 50 micro (radius) ice crystal
 379 !       =1.0e-4 [m]
 380 ! xmi: mass of one cloud ice crystal
 381 !      =4.19e-13 [kg] = 4.19e-10 [g]
 382 ! oxmi=1.0/xmi
 383 !
 384 ! xni0=1.0e-2 [m-3] The value given in Lin et al. is wrong.(see
 385 !                   Hsie et al.(1980) and Rutledge and Hobbs(1983) )
 386 ! bni=0.5 [K-1]
 387 ! xmnin: mass of a natural ice nucleus
 388 !    = 1.05e-18 [kg] = 1.05e-15 [g] This values is suggested by
 389 !                    Hsie et al. (1980)
 390 !    = 1.0e-12 [kg] suggested by Rutlegde and Hobbs (1983)
 391 ! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3
 392 ! consta: constant in empirical formular for terminal
 393 !         velocity of raindrop
 394 !         =2115.0 [cm**(1-b)/s] = 2115.0*0.01**(1-b) [m**(1-b)/s]
 395 ! constb: constant in empirical formular for terminal
 396 !         velocity of raindrop
 397 !         =0.8
 398 ! xnor: intercept parameter of the raindrop size distribution
 399 !       = 0.08 cm-4 = 8.0e6 m-4
 400 ! araut: time sacle for autoconversion of cloud water to raindrops
 401 !       =1.0e-3 [s-1]
 402 ! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg]
 403 ! vf1r: ventilation factors for rain =0.78
 404 ! vf2r: ventilation factors for rain =0.31
 405 ! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3
 406 ! constc: constant in empirical formular for terminal
 407 !         velocity of snow
 408 !         =152.93 [cm**(1-d)/s] = 152.93*0.01**(1-d) [m**(1-d)/s]
 409 ! constd: constant in empirical formular for terminal
 410 !         velocity of snow
 411 !         =0.25
 412 ! xnos: intercept parameter of the snowflake size distribution
 413 ! vf1s: ventilation factors for snow =0.78
 414 ! vf2s: ventilation factors for snow =0.31
 415 !
 416 !----------------------------------------------------------------------
 417
 418   INTEGER, INTENT(IN   )               :: kts, kte, i, j
 419
 420   REAL,    DIMENSION( kts:kte ),                                      &
 421            INTENT(INOUT)               :: qvz, qlz, qrz, qiz, qsz,    &
 422                                           qndropz,                    &
 423                                           qgz, thz
 424
 425   REAL,    DIMENSION( kts:kte ),                                      &
 426            INTENT(IN   )               :: tothz, rho, orho, sqrho,    &
 427                                           prez, zz, dzw
 428
 429   REAL,    INTENT(IN   )               :: dt, grav, cp, Rair, rvapor, &
 430                                           XLS, XLV, XLF, rhowater,    &
 431                                           rhosnow,EP2,SVP1,SVP2,SVP3,SVPT0
 432
 433   REAL,    DIMENSION( kts:kte ), INTENT(OUT)               :: preclw, &
 434                     precrz, preciz, precsz, precgz
 435
 436   REAL,    INTENT(INOUT)               :: pptrain, pptsnow, pptgraul, pptice
 437
 438   REAL,    INTENT(IN   )               :: zsfc
 439   logical, intent(in)                  :: flag_qndrop !sg
 440
 441 ! local vars
 442
 443    REAL                                :: obp4, bp3, bp5, bp6, odp4,  &
 444                                           dp3, dp5, dp5o2
 445
 446
 447 ! temperary vars
 448
 449    REAL                                :: tmp, tmp0, tmp1, tmp2,tmp3,  &
 450                                           tmp4,delta2,delta3, delta4,  &
 451                                           tmpa,tmpb,tmpc,tmpd,alpha1,  &
 452                                           qic, abi,abr, abg, odtberg,  &
 453                                           vti50,eiw,eri,esi,esr, esw,  &
 454                                           erw,delrs,term0,term1,araut, &
 455                                           constg2, vf1r, vf2r,alpha2,  &
 456                                           Ap, Bp, egw, egi, egs, egr,  &
 457                                           constg, gdelta4, g1sdelt4,   &
 458                                           factor, tmp_r, tmp_s,tmp_g,  &
 459                                           qlpqi, rsat, a1, a2, xnin
 460
 461   INTEGER                              :: k
 462 !
 463   REAL, DIMENSION( kts:kte )    ::  oprez, tem, temcc, theiz, qswz,    &
 464                                     qsiz, qvoqswz, qvoqsiz, qvzodt,    &
 465                                     qlzodt, qizodt, qszodt, qrzodt,    &
 466                                     qgzodt
 467
 468   REAL, DIMENSION( kts:kte )    :: psnow, psaut, psfw,  psfi,  praci,  &
 469                                    piacr, psaci, psacw, psdep, pssub,  &
 470                                    pracs, psacr, psmlt, psmltevp,      &
 471                                    prain, praut, pracw, prevp, pvapor, &
 472                                    pclw,  pladj, pcli,  pimlt, pihom,  &
 473                                    pidw,  piadj, pgraupel, pgaut,      &
 474                                    pgfr,  pgacw, pgaci, pgacr, pgacs,  &
 475                                    pgacip,pgacrp,pgacsp,pgwet, pdry,   &
 476                                    pgsub, pgdep, pgmlt, pgmltevp,      &
 477                                    qschg, qgchg
 478 !
 479
 480   REAL, DIMENSION( kts:kte )    :: qvsbar, rs0, viscmu, visc, diffwv,  &
 481                                    schmidt, xka
 482
 483   REAL, DIMENSION( kts:kte )    :: vtr, vts, vtg,                      &
 484                                    vtrold, vtsold, vtgold, vtiold,     &
 485                                    xlambdar, xlambdas, xlambdag,       &
 486                                    olambdar, olambdas, olambdag
 487
 488   REAL                          :: episp0k, dtb, odtb, pi, pio4,       &
 489                                    pio6, oxLf, xLvocp, xLfocp, consta, &
 490                                    constc, ocdrag, gambp4, gamdp4,     &
 491                                    gam4pt5, Cpor, oxmi, gambp3, gamdp3,&
 492                                    gambp6, gam3pt5, gam2pt75, gambp5o2,&
 493                                    gamdp5o2, cwoxlf, ocp, xni50, es
 494 !
 495   REAL                          :: qvmin=1.e-20
 496   REAL                          :: gindex
 497   REAL                          :: temc1,save1,save2,xni50mx
 498
 499 ! for terminal velocity flux
 500
 501   INTEGER                       :: min_q, max_q
 502   REAL                          :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
 503   LOGICAL                       :: notlast
 504 !
 505   REAL                          :: tmp_tem, tmp_temcc !bloss
 506 !
 507   real :: liqconc, dis, beta, kappa, p0, xc, capn,rhocgs
 508
 509 !+---+-----------------------------------------------------------------+
 510       INTEGER:: NCALL=0
 511
 512       IF (NCALL .EQ. 0) THEN
 513 !..Set these variables needed for computing radar reflectivity.  These
 514 !.. get used within radar_init to create other variables used in the
 515 !.. radar module.
 516          xam_r = 3.14159*rhowater/6.
 517          xbm_r = 3.
 518          xmu_r = 0.
 519          xam_s = 3.14159*rhosnow/6.
 520          xbm_s = 3.
 521          xmu_s = 0.
 522          xam_g = 3.14159*rhograul/6.
 523          xbm_g = 3.
 524          xmu_g = 0.
 525
 526          call radar_init
 527          NCALL = 1
 528       ENDIF
 529 !+---+-----------------------------------------------------------------+
 530
 531 !sg: begin
 532 ! liqconc = liquid water content (g cm^-3)
 533 ! capn = droplet number concentration (# cm^-3)
 534 ! dis = relative dispersion (dimensionless) between 0.2 and 1.
 535 ! Written by Yangang Liu based on Liu et al., GRL 32, 2005.
 536 ! Autoconversion rate p = p0 * (threshold function)
 537 ! p0 = "base" autoconversion rate (g cm^-3 s^-1)
 538 ! kappa = constant in Long kernel = [kappa2 * (3/(4*pi*rhow))^3] in Liu papers
 539 ! beta = Condensation rate constant = (beta6)^6 in Liu papers
 540 ! xc = Normalized critical mass
 541 ! ***********************************************************
 542   if(flag_qndrop)then
 543      dis = 0.5 ! droplet dispersion, set to 0.5 per SG 8-Nov-2006
 544 !    Give  empirical constants
 545      kappa=1.1d10
 546 !    Calculate Condensation rate constant
 547      beta = (1.0d0+3.0d0*dis**2)*(1.0d0+4.0d0*dis**2)*    &
 548          (1.0d0+5.0d0*dis**2)/((1.0d0+dis**2)*(1.0d0+2.0d0*dis**2))
 549   endif
 550 !sg: end
 551
 552    dtb=dt
 553    odtb=1./dtb
 554    pi=acos(-1.)
 555    pio4=acos(-1.)/4.
 556    pio6=acos(-1.)/6.
 557    ocp=1./cp
 558    oxLf=1./xLf
 559    xLvocp=xLv/cp
 560    xLfocp=xLf/cp
 561    consta=2115.0*0.01**(1-constb)
 562    constc=152.93*0.01**(1-constd)
 563    ocdrag=1./Cdrag
 564 !  episp0k=RH*episp0
 565    episp0k=RH*ep2*1000.*svp1
 566 !
 567    gambp4=ggamma(constb+4.)
 568    gamdp4=ggamma(constd+4.)
 569    gam4pt5=ggamma(4.5)
 570    Cpor=cp/Rair
 571    oxmi=1.0/xmi
 572    gambp3=ggamma(constb+3.)
 573    gamdp3=ggamma(constd+3.)
 574    gambp6=ggamma(constb+6)
 575    gam3pt5=ggamma(3.5)
 576    gam2pt75=ggamma(2.75)
 577    gambp5o2=ggamma((constb+5.)/2.)
 578    gamdp5o2=ggamma((constd+5.)/2.)
 579    cwoxlf=cw/xlf
 580    delta2=0.
 581    delta3=0.
 582    delta4=0.
 583 !
 584 !-----------------------------------------------------------------------
 585 !     oprez       1./prez ( prez : pressure)
 586 !     qsw         saturated mixing ratio on water surface
 587 !     qsi         saturated mixing ratio on ice surface
 588 !     episp0k     RH*e*saturated pressure at 273.15 K
 589 !     qvoqsw      qv/qsw
 590 !     qvoqsi      qv/qsi
 591 !     qvzodt      qv/dt
 592 !     qlzodt      ql/dt
 593 !     qizodt      qi/dt
 594 !     qszodt      qs/dt
 595 !     qrzodt      qr/dt
 596 !     qgzodt      qg/dt
 597 !
 598 !     temcc       temperature in dregee C
 599 !
 600
 601       obp4=1.0/(constb+4.0)
 602       bp3=constb+3.0
 603       bp5=constb+5.0
 604       bp6=constb+6.0
 605       odp4=1.0/(constd+4.0)
 606       dp3=constd+3.0
 607       dp5=constd+5.0
 608       dp5o2=0.5*(constd+5.0)
 609 !
 610       do k=kts,kte
 611          oprez(k)=1./prez(k)
 612       enddo
 613
 614       do k=kts,kte
 615          qlz(k)=amax1( 0.0,qlz(k) )
 616          qiz(k)=amax1( 0.0,qiz(k) )
 617          qvz(k)=amax1( qvmin,qvz(k) )
 618          qsz(k)=amax1( 0.0,qsz(k) )
 619          qrz(k)=amax1( 0.0,qrz(k) )
 620          qgz(k)=amax1( 0.0,qgz(k) )
 621          qndropz(k)=amax1( 0.0,qndropz(k) )     !sg
 622 !
 623          tem(k)=thz(k)*tothz(k)
 624          temcc(k)=tem(k)-273.15
 625 !
 626 !        qswz(k)=episp0k*oprez(k)* &
 627 !               exp( svp2*temcc(k)/(tem(k)-svp3) )
 628          es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
 629          qswz(k)=ep2*es/(prez(k)-es)
 630          if (tem(k) .lt. 273.15 ) then
 631 !           qsiz(k)=episp0k*oprez(k)* &
 632 !                  exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
 633             es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
 634             qsiz(k)=ep2*es/(prez(k)-es)
 635             if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
 636          else
 637             qsiz(k)=qswz(k)
 638          endif
 639 !
 640          qvoqswz(k)=qvz(k)/qswz(k)
 641          qvoqsiz(k)=qvz(k)/qsiz(k)
 642
 643          theiz(k)=thz(k)+(xlvocp*qvz(k)-xlfocp*qiz(k))/tothz(k)
 644       enddo
 645
 646
 647 !
 648 !
 649 !-----------------------------------------------------------------------
 650 ! In this simple stable cloud parameterization scheme, only five
 651 ! forms of water substance (water vapor, cloud water, cloud ice,
 652 ! rain and snow are considered. The prognostic variables are total
 653 ! water (qp),cloud water (ql), and cloud ice (qi). Rain and snow are
 654 ! diagnosed following Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7).
 655 ! the micro physics are based on (1) Hsie et al.,1980, JAM, 950-977 ;
 656 ! (2) Lin et al., 1983, JAM, 1065-1092 ; (3) Rutledge and Hobbs, 1983,
 657 ! JAS, 1185-1206 ; (4) Rutledge and Hobbs, 1984, JAS, 2949-2972.
 658 !-----------------------------------------------------------------------
 659 !
 660 ! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3
 661 ! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3
 662 ! xnor: intercept parameter of the raindrop size distribution
 663 !       = 0.08 cm-4 = 8.0e6 m-4
 664 ! xnos: intercept parameter of the snowflake size distribution
 665 !       = 0.03 cm-4 = 3.0e6 m-4
 666 ! xnog: intercept parameter of the graupel size distribution
 667 !       = 4.0e-4 cm-4 = 4.0e4 m-4
 668 ! consta: constant in empirical formular for terminal
 669 !         velocity of raindrop
 670 !         =2115.0 [cm**(1-b)/s] = 2115.0*0.01**(1-b) [m**(1-b)/s]
 671 ! constb: constant in empirical formular for terminal
 672 !         velocity of raindrop
 673 !         =0.8
 674 ! constc: constant in empirical formular for terminal
 675 !         velocity of snow
 676 !         =152.93 [cm**(1-d)/s] = 152.93*0.01**(1-d) [m**(1-d)/s]
 677 ! constd: constant in empirical formular for terminal
 678 !         velocity of snow
 679 !         =0.25
 680 ! avisc: constant in empirical formular for dynamic viscosity of air
 681 !         =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
 682 ! adiffwv: constant in empirical formular for diffusivity of water
 683 !          vapor in air
 684 !          = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
 685 ! axka: constant in empirical formular for thermal conductivity of air
 686 !       = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
 687 ! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
 688 !      = 1.0e-3 g/g = 1.0e-3 kg/gk
 689 ! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg]
 690 !      = 2.0e-3 g/g = 2.0e-3 kg/gk
 691 ! qs0: mixing ratio threshold for snow aggregation
 692 !      = 6.0e-4 g/g = 6.0e-4 kg/gk
 693 ! xmi50: mass of a 50 micron ice crystal
 694 !        = 4.8e-10 [kg] =4.8e-7 [g]
 695 ! xmi40: mass of a 40 micron ice crystal
 696 !        = 2.46e-10 [kg] = 2.46e-7 [g]
 697 ! di50: diameter of a 50 micro (radius) ice crystal
 698 !       =1.0e-4 [m]
 699 ! xmi: mass of one cloud ice crystal
 700 !      =4.19e-13 [kg] = 4.19e-10 [g]
 701 ! oxmi=1.0/xmi
 702 !
 703
 704
 705 ! if gindex=1.0 include graupel
 706 ! if gindex=0. no graupel
 707 !
 708 !
 709       do k=kts,kte
 710          psnow(k)=0.0
 711          psaut(k)=0.0
 712          psfw(k)=0.0
 713          psfi(k)=0.0
 714          praci(k)=0.0
 715          piacr(k)=0.0
 716          psaci(k)=0.0
 717          psacw(k)=0.0
 718          psdep(k)=0.0
 719          pssub(k)=0.0
 720          pracs(k)=0.0
 721          psacr(k)=0.0
 722          psmlt(k)=0.0
 723          psmltevp(k)=0.0
 724 !
 725          prain(k)=0.0
 726          praut(k)=0.0
 727          pracw(k)=0.0
 728          prevp(k)=0.0
 729 !
 730          pvapor(k)=0.0
 731 !
 732          pclw(k)=0.0
 733          preclw(k)=0.0       !sg
 734          pladj(k)=0.0
 735 !
 736          pcli(k)=0.0
 737          pimlt(k)=0.0
 738          pihom(k)=0.0
 739          pidw(k)=0.0
 740          piadj(k)=0.0
 741       enddo
 742
 743 !
 744 !!! graupel
 745 !
 746       do k=kts,kte
 747          pgraupel(k)=0.0
 748          pgaut(k)=0.0
 749          pgfr(k)=0.0
 750          pgacw(k)=0.0
 751          pgaci(k)=0.0
 752          pgacr(k)=0.0
 753          pgacs(k)=0.0
 754          pgacip(k)=0.0
 755          pgacrP(k)=0.0
 756          pgacsp(k)=0.0
 757          pgwet(k)=0.0
 758          pdry(k)=0.0
 759          pgsub(k)=0.0
 760          pgdep(k)=0.0
 761          pgmlt(k)=0.0
 762          pgmltevp(k)=0.0
 763          qschg(k)=0.
 764          qgchg(k)=0.
 765       enddo
 766 !
 767 !
 768 ! Set rs0=episp0*oprez(k)
 769 ! episp0=e*saturated pressure at 273.15 K
 770 ! e     = 0.622
 771 !
 772       DO k=kts,kte
 773          rs0(k)=ep2*1000.*svp1/(prez(k)-1000.*svp1)
 774       END DO
 775 !
 776 !***********************************************************************
 777 ! Calculate precipitation fluxes due to terminal velocities.
 778 !***********************************************************************
 779 !
 780 !- Calculate termianl velocity (vt?)  of precipitation q?z
 781 !- Find maximum vt? to determine the small delta t
 782 !
 783 !-- rain
 784 !
 785     t_del_tv=0.
 786     del_tv=dtb
 787     notlast=.true.
 788     DO while (notlast)
 789 !
 790       min_q=kte
 791       max_q=kts-1
 792 !
 793       do k=kts,kte-1
 794          if (qrz(k) .gt. 1.0e-8) then
 795             min_q=min0(min_q,k)
 796             max_q=max0(max_q,k)
 797             tmp1=sqrt(pi*rhowater*xnor/rho(k)/qrz(k))
 798             tmp1=sqrt(tmp1)
 799             vtrold(k)=o6*consta*gambp4*sqrho(k)/tmp1**constb
 800             if (k .eq. 1) then
 801                del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtrold(k))
 802             else
 803                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtrold(k))
 804             endif
 805          else
 806             vtrold(k)=0.
 807          endif
 808       enddo
 809
 810       if (max_q .ge. min_q) then
 811 !
 812 !- Check if the summation of the small delta t >=  big delta t
 813 !             (t_del_tv)          (del_tv)             (dtb)
 814
 815          t_del_tv=t_del_tv+del_tv
 816 !
 817          if ( t_del_tv .ge. dtb ) then
 818               notlast=.false.
 819               del_tv=dtb+del_tv-t_del_tv
 820          endif
 821 !
 822          fluxin=0.
 823          do k=max_q,min_q,-1
 824             fluxout=rho(k)*vtrold(k)*qrz(k)
 825             flux=(fluxin-fluxout)/rho(k)/dzw(k)
 826             tmpqrz=qrz(k)
 827             qrz(k)=qrz(k)+del_tv*flux
 828             fluxin=fluxout
 829          enddo
 830          if (min_q .eq. 1) then
 831             pptrain=pptrain+fluxin*del_tv
 832          else
 833             qrz(min_q-1)=qrz(min_q-1)+del_tv*  &
 834                           fluxin/rho(min_q-1)/dzw(min_q-1)
 835          endif
 836 !
 837       else
 838          notlast=.false.
 839       endif
 840     ENDDO
 841
 842 !
 843 !-- snow
 844 !
 845     t_del_tv=0.
 846     del_tv=dtb
 847     notlast=.true.
 848
 849     DO while (notlast)
 850 !
 851       min_q=kte
 852       max_q=kts-1
 853 !
 854       do k=kts,kte-1
 855          if (qsz(k) .gt. 1.0e-8) then
 856             min_q=min0(min_q,k)
 857             max_q=max0(max_q,k)
 858             tmp1=sqrt(pi*rhosnow*xnos/rho(k)/qsz(k))
 859             tmp1=sqrt(tmp1)
 860             vtsold(k)=o6*constc*gamdp4*sqrho(k)/tmp1**constd
 861             if (k .eq. 1) then
 862                del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtsold(k))
 863             else
 864                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtsold(k))
 865             endif
 866          else
 867             vtsold(k)=0.
 868          endif
 869       enddo
 870
 871       if (max_q .ge. min_q) then
 872 !
 873 !
 874 !- Check if the summation of the small delta t >=  big delta t
 875 !             (t_del_tv)          (del_tv)             (dtb)
 876
 877          t_del_tv=t_del_tv+del_tv
 878
 879          if ( t_del_tv .ge. dtb ) then
 880               notlast=.false.
 881               del_tv=dtb+del_tv-t_del_tv
 882          endif
 883 !
 884          fluxin=0.
 885          do k=max_q,min_q,-1
 886             fluxout=rho(k)*vtsold(k)*qsz(k)
 887             flux=(fluxin-fluxout)/rho(k)/dzw(k)
 888             qsz(k)=qsz(k)+del_tv*flux
 889             qsz(k)=amax1(0.,qsz(k))
 890             fluxin=fluxout
 891          enddo
 892          if (min_q .eq. 1) then
 893             pptsnow=pptsnow+fluxin*del_tv
 894          else
 895             qsz(min_q-1)=qsz(min_q-1)+del_tv*  &
 896                          fluxin/rho(min_q-1)/dzw(min_q-1)
 897          endif
 898 !
 899       else
 900          notlast=.false.
 901       endif
 902
 903     ENDDO
 904 !
 905 !-- grauupel
 906 !
 907     t_del_tv=0.
 908     del_tv=dtb
 909     notlast=.true.
 910 !
 911     DO while (notlast)
 912 !
 913       min_q=kte
 914       max_q=kts-1
 915 !
 916       do k=kts,kte-1
 917          if (qgz(k) .gt. 1.0e-8) then
 918             min_q=min0(min_q,k)
 919             max_q=max0(max_q,k)
 920             tmp1=sqrt(pi*rhograul*xnog/rho(k)/qgz(k))
 921             tmp1=sqrt(tmp1)
 922             term0=sqrt(4.*grav*rhograul*0.33334/rho(k)/cdrag)
 923             vtgold(k)=o6*gam4pt5*term0*sqrt(1./tmp1)
 924             if (k .eq. 1) then
 925                del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtgold(k))
 926             else
 927                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtgold(k))
 928             endif
 929          else
 930             vtgold(k)=0.
 931          endif
 932       enddo
 933
 934       if (max_q .ge. min_q) then
 935 !
 936 !
 937 !- Check if the summation of the small delta t >=  big delta t
 938 !             (t_del_tv)          (del_tv)             (dtb)
 939
 940          t_del_tv=t_del_tv+del_tv
 941
 942          if ( t_del_tv .ge. dtb ) then
 943               notlast=.false.
 944               del_tv=dtb+del_tv-t_del_tv
 945          endif
 946
 947 !
 948          fluxin=0.
 949          do k=max_q,min_q,-1
 950             fluxout=rho(k)*vtgold(k)*qgz(k)
 951             flux=(fluxin-fluxout)/rho(k)/dzw(k)
 952             qgz(k)=qgz(k)+del_tv*flux
 953             qgz(k)=amax1(0.,qgz(k))
 954             fluxin=fluxout
 955          enddo
 956          if (min_q .eq. 1) then
 957             pptgraul=pptgraul+fluxin*del_tv
 958          else
 959             qgz(min_q-1)=qgz(min_q-1)+del_tv*  &
 960                          fluxin/rho(min_q-1)/dzw(min_q-1)
 961          endif
 962 !
 963       else
 964          notlast=.false.
 965       endif
 966 !
 967    ENDDO
 968
 969 !
 970 !-- cloud ice  (03/21/02) follow Vaughan T.J. Phillips at GFDL
 971 !
 972     t_del_tv=0.
 973     del_tv=dtb
 974     notlast=.true.
 975 !
 976     DO while (notlast)
 977 !
 978       min_q=kte
 979       max_q=kts-1
 980 !
 981       do k=kts,kte-1
 982          if (qiz(k) .gt. 1.0e-8) then
 983             min_q=min0(min_q,k)
 984             max_q=max0(max_q,k)
 985             vtiold(k)= 3.29 * (rho(k)* qiz(k))** 0.16  ! Heymsfield and Donner
 986             if (k .eq. 1) then
 987                del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtiold(k))
 988             else
 989                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtiold(k))
 990             endif
 991          else
 992             vtiold(k)=0.
 993          endif
 994       enddo
 995
 996       if (max_q .ge. min_q) then
 997 !
 998 !
 999 !- Check if the summation of the small delta t >=  big delta t
1000 !             (t_del_tv)          (del_tv)             (dtb)
1001
1002          t_del_tv=t_del_tv+del_tv
1003
1004          if ( t_del_tv .ge. dtb ) then
1005               notlast=.false.
1006               del_tv=dtb+del_tv-t_del_tv
1007          endif
1008
1009          fluxin=0.
1010          do k=max_q,min_q,-1
1011             fluxout=rho(k)*vtiold(k)*qiz(k)
1012             flux=(fluxin-fluxout)/rho(k)/dzw(k)
1013             qiz(k)=qiz(k)+del_tv*flux
1014             qiz(k)=amax1(0.,qiz(k))
1015             fluxin=fluxout
1016          enddo
1017          if (min_q .eq. 1) then
1018             pptice=pptice+fluxin*del_tv
1019          else
1020             qiz(min_q-1)=qiz(min_q-1)+del_tv*  &
1021                          fluxin/rho(min_q-1)/dzw(min_q-1)
1022          endif
1023 !
1024       else
1025          notlast=.false.
1026       endif
1027 !
1028    ENDDO
1029    do k=kts,kte-1                         !sg beg
1030       precrz(k)=rho(k)*vtrold(k)*qrz(k)
1031       preciz(k)=rho(k)*vtiold(k)*qiz(k)
1032       precsz(k)=rho(k)*vtsold(k)*qsz(k)
1033       precgz(k)=rho(k)*vtgold(k)*qgz(k)
1034    enddo                                  !sg end
1035    precrz(kte)=0. !wig - top level never set for vtXold vars
1036    preciz(kte)=0. !wig
1037    precsz(kte)=0. !wig
1038    precgz(kte)=0. !wig
1039
1040
1041 ! Microphysics processes
1042 !
1043       DO 2000 k=kts,kte
1044 !
1045          qvzodt(k)=amax1( 0.0,odtb*qvz(k) )
1046          qlzodt(k)=amax1( 0.0,odtb*qlz(k) )
1047          qizodt(k)=amax1( 0.0,odtb*qiz(k) )
1048          qszodt(k)=amax1( 0.0,odtb*qsz(k) )
1049          qrzodt(k)=amax1( 0.0,odtb*qrz(k) )
1050          qgzodt(k)=amax1( 0.0,odtb*qgz(k) )
1051 !***********************************************************************
1052 !*****   diagnose mixing ratios (qrz,qsz), terminal                *****
1053 !*****   velocities (vtr,vts), and slope parameters in size        *****
1054 !*****   distribution(xlambdar,xlambdas) of rain and snow          *****
1055 !*****   follows Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7)   *****
1056 !***********************************************************************
1057 !
1058 !**** assuming no cloud water can exist in the top two levels due to
1059 !**** radiation consideration
1060 !
1061 !!  if
1062 !!     unsaturated,
1063 !!     no cloud water, rain, ice, snow and graupel
1064 !!  then
1065 !!     skip these processes and jump to line 2000
1066 !
1067 !
1068         tmp=qiz(k)+qlz(k)+qsz(k)+qrz(k)+qgz(k)*gindex
1069         if( qvz(k)+qlz(k)+qiz(k) .lt. qsiz(k)  &
1070             .and. tmp .eq. 0.0 ) go to 2000
1071
1072 !! calculate terminal velocity of rain
1073 !
1074         if (qrz(k) .gt. 1.0e-8) then
1075             tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k))
1076             xlambdar(k)=sqrt(tmp1)
1077             olambdar(k)=1.0/xlambdar(k)
1078             vtrold(k)=o6*consta*gambp4*sqrho(k)*olambdar(k)**constb
1079         else
1080             vtrold(k)=0.
1081             olambdar(k)=0.
1082         endif
1083 !
1084 !       if (qrz(k) .gt. 1.0e-12) then
1085         if (qrz(k) .gt. 1.0e-8) then
1086             tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k))
1087             xlambdar(k)=sqrt(tmp1)
1088             olambdar(k)=1.0/xlambdar(k)
1089             vtr(k)=o6*consta*gambp4*sqrho(k)*olambdar(k)**constb
1090         else
1091             vtr(k)=0.
1092             olambdar(k)=0.
1093         endif
1094 !
1095 !! calculate terminal velocity of snow
1096 !
1097         if (qsz(k) .gt. 1.0e-8) then
1098             tmp1=sqrt(pi*rhosnow*xnos*orho(k)/qsz(k))
1099             xlambdas(k)=sqrt(tmp1)
1100             olambdas(k)=1.0/xlambdas(k)
1101             vtsold(k)=o6*constc*gamdp4*sqrho(k)*olambdas(k)**constd
1102         else
1103             vtsold(k)=0.
1104             olambdas(k)=0.
1105         endif
1106 !
1107 !       if (qsz(k) .gt. 1.0e-12) then
1108         if (qsz(k) .gt. 1.0e-8) then
1109             tmp1=sqrt(pi*rhosnow*xnos*orho(k)/qsz(k))
1110             xlambdas(k)=sqrt(tmp1)
1111             olambdas(k)=1.0/xlambdas(k)
1112             vts(k)=o6*constc*gamdp4*sqrho(k)*olambdas(k)**constd
1113         else
1114             vts(k)=0.
1115             olambdas(k)=0.
1116         endif
1117 !
1118 !! calculate terminal velocity of graupel
1119 !
1120         if (qgz(k) .gt. 1.0e-8) then
1121             tmp1=sqrt( pi*rhograul*xnog*orho(k)/qgz(k))
1122             xlambdag(k)=sqrt(tmp1)
1123             olambdag(k)=1.0/xlambdag(k)
1124             term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag)
1125             vtgold(k)=o6*gam4pt5*term0*sqrt(olambdag(k))
1126         else
1127             vtgold(k)=0.
1128             olambdag(k)=0.
1129         endif
1130 !
1131 !       if (qgz(k) .gt. 1.0e-12) then
1132         if (qgz(k) .gt. 1.0e-8) then
1133             tmp1=sqrt( pi*rhograul*xnog*orho(k)/qgz(k))
1134             xlambdag(k)=sqrt(tmp1)
1135             olambdag(k)=1.0/xlambdag(k)
1136             term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag)
1137             vtg(k)=o6*gam4pt5*term0*sqrt(olambdag(k))
1138         else
1139             vtg(k)=0.
1140             olambdag(k)=0.
1141         endif
1142 !
1143 !***********************************************************************
1144 !*****  compute viscosity,difusivity,thermal conductivity, and    ******
1145 !*****  Schmidt number                                            ******
1146 !***********************************************************************
1147 !c------------------------------------------------------------------
1148 !c      viscmu: dynamic viscosity of air kg/m/s
1149 !c      visc: kinematic viscosity of air = viscmu/rho (m2/s)
1150 !c      avisc=1.49628e-6 kg/m/s=1.49628e-5 g/cm/s
1151 !c      viscmu=1.718e-5 kg/m/s in RH
1152 !c      diffwv: Diffusivity of water vapor in air
1153 !c      adiffwv = 8.7602e-5 (8.794e-5 in MM5) kgm/s3
1154 !c              = 8.7602 (8.794 in MM5) gcm/s3
1155 !c      diffwv(k)=2.26e-5 m2/s
1156 !c      schmidt: Schmidt number=visc/diffwv
1157 !c      xka: thermal conductivity of air J/m/s/K (Kgm/s3/K)
1158 !c      xka(k)=2.43e-2 J/m/s/K in RH
1159 !c      axka=1.4132e3 (1.414e3 in MM5) m2/s2/k = 1.4132e7 cm2/s2/k
1160 !c------------------------------------------------------------------
1161
1162         viscmu(k)=avisc*tem(k)**1.5/(tem(k)+120.0)
1163         visc(k)=viscmu(k)*orho(k)
1164         diffwv(k)=adiffwv*tem(k)**1.81*oprez(k)
1165         schmidt(k)=visc(k)/diffwv(k)
1166         xka(k)=axka*viscmu(k)
1167
1168         if (tem(k) .lt. 273.15) then
1169
1170 !
1171 !***********************************************************************
1172 !*********        snow production processes for T < 0 C       **********
1173 !***********************************************************************
1174 !c
1175 !c (1) ICE CRYSTAL AGGREGATION TO SNOW (Psaut): Lin (21)
1176 !c!    psaut=alpha1*(qi-qi0)
1177 !c!    alpha1=1.0e-3*exp(0.025*(T-T0))
1178 !c
1179 !          alpha1=1.0e-3*exp( 0.025*temcc(k) )
1180
1181            alpha1=1.0e-3*exp( 0.025*temcc(k) )
1182 !
1183            if(temcc(k) .lt. -20.0) then
1184              tmp1=-7.6+4.0*exp( -0.2443e-3*(abs(temcc(k))-20)**2.455 )
1185              qic=1.0e-3*exp(tmp1)*orho(k)
1186            else
1187              qic=qi0
1188            end if
1189 !testing
1190 !          tmp1=amax1( 0.0,alpha1*(qiz(k)-qic) )
1191 !          psaut(k)=amin1( tmp1,qizodt(k) )
1192
1193            tmp1=odtb*(qiz(k)-qic)*(1.0-exp(-alpha1*dtb))
1194            psaut(k)=amax1( 0.0,tmp1 )
1195
1196 !c
1197 !c (2) BERGERON PROCESS TRANSFER OF CLOUD WATER TO SNOW (Psfw)
1198 !c     this process only considered when -31 C < T < 0 C
1199 !c     Lin (33) and Hsie (17)
1200 !c
1201 !c!
1202 !c!    parama1 and parama2 functions must be user supplied
1203 !c!
1204
1205 ! testing
1206           if( qlz(k) .gt. 1.0e-10 ) then
1207             temc1=amax1(-30.99,temcc(k))
1208 !           print*,'temc1',temc1,qlz(k)
1209             a1=parama1( temc1 )
1210             a2=parama2( temc1 )
1211             tmp1=1.0-a2
1212 !! change unit from cgs to mks
1213             a1=a1*0.001**tmp1
1214 !c!   dtberg is the time needed for a crystal to grow from 40 to 50 um
1215 !c !  odtberg=1.0/dtberg
1216             odtberg=(a1*tmp1)/(xmi50**tmp1-xmi40**tmp1)
1217 !
1218 !c!   compute terminal velocity of a 50 micron ice cystal
1219 !
1220             vti50=constc*di50**constd*sqrho(k)
1221 !
1222             eiw=1.0
1223             save1=a1*xmi50**a2
1224             save2=0.25*pi*eiw*rho(k)*di50*di50*vti50
1225 !
1226             tmp2=( save1 + save2*qlz(k) )
1227 !
1228 !! maximum number of 50 micron crystals limited by the amount
1229 !!  of supercool water
1230 !
1231             xni50mx=qlzodt(k)/tmp2
1232 !
1233 !!   number of 50 micron crystals produced
1234 !
1235 !
1236             xni50=qiz(k)*( 1.0-exp(-dtb*odtberg) )/xmi50
1237             xni50=amin1(xni50,xni50mx)
1238 !
1239             tmp3=odtb*tmp2/save2*( 1.0-exp(-save2*xni50*dtb) )
1240             psfw(k)=amin1( tmp3,qlzodt(k) )
1241 !testing
1242 !           psfw(k)=0.
1243
1244 !0915     if( temcc(k).gt.-30.99 ) then
1245 !0915        a1=parama1( temcc(k) )
1246 !0915        a2=parama2( temcc(k) )
1247 !0915        tmp1=1.0-a2
1248 !!     change unit from cgs to mks
1249 !0915        a1=a1*0.001**tmp1
1250
1251 !c!    dtberg is the time needed for a crystal to grow from 40 to 50 um
1252 !c!    odtberg=1.0/dtberg
1253 !0915        odtberg=(a1*tmp1)/(xmi50**tmp1-xmi40**tmp1)
1254
1255 !c!    number of 50 micron crystals produced
1256 !0915        xni50=qiz(k)*dtb*odtberg/xmi50
1257
1258 !c!    need to calculate the terminal velocity of a 50 micron
1259 !c!    ice cystal
1260 !0915        vti50=constc*di50**constd*sqrho(k)
1261 !0915        eiw=1.0
1262 !0915        tmp2=xni50*( a1*xmi50**a2 + &
1263 !0915             0.25*qlz(k)*pi*eiw*rho(k)*di50*di50*vti50 )
1264 !0915        psfw(k)=amin1( tmp2,qlzodt(k) )
1265 !0915        psfw(k)=0.
1266 !c
1267 !c (3) REDUCTION OF CLOUD ICE BY BERGERON PROCESS (Psfi): Lin (34)
1268 !c     this process only considered when -31 C < T < 0 C
1269 !c
1270             tmp1=xni50*xmi50-psfw(k)
1271             psfi(k)=amin1(tmp1,qizodt(k))
1272 ! testing
1273 !           psfi(k)=0.
1274           end if
1275 !
1276
1277 !0915        tmp1=qiz(k)*odtberg
1278 !0915        psfi(k)=amin1(tmp1,qizodt(k))
1279 ! testing
1280 !0915        psfi(k)=0.
1281 !0915     end if
1282 !
1283           if(qrz(k) .le. 0.0) go to 1000
1284 !
1285 ! Processes (4) and (5) only need when qrz > 0.0
1286 !
1287 !c
1288 !c (4) CLOUD ICE ACCRETION BY RAIN (Praci): Lin (25)
1289 !c     may produce snow or graupel
1290 !c
1291           eri=1.0
1292 !0915     tmp1=qiz(k)*pio4*eri*xnor*consta*sqrho(k)
1293 !0915     tmp2=tmp1*gambp3*olambdar(k)**bp3
1294 !0915     praci(k)=amin1( tmp2,qizodt(k) )
1295
1296           save1=pio4*eri*xnor*consta*sqrho(k)
1297           tmp1=save1*gambp3*olambdar(k)**bp3
1298           praci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) )
1299
1300 !c
1301 !c (5) RAIN ACCRETION BY CLOUD ICE (Piacr): Lin (26)
1302 !c
1303 !0915     tmp2=tmp1*rho(k)*pio6*rhowater*gambp6*oxmi* &
1304 !0915              olambdar(k)**bp6
1305 !0915     piacr(k)=amin1( tmp2,qrzodt(k) )
1306
1307           tmp2=qiz(k)*save1*rho(k)*pio6*rhowater*gambp6*oxmi* &
1308                    olambdar(k)**bp6
1309           piacr(k)=amin1( tmp2,qrzodt(k) )
1310
1311 !
1312 1000      continue
1313 !
1314           if(qsz(k) .le. 0.0) go to 1200
1315 !
1316 ! Compute the following processes only when qsz > 0.0
1317 !
1318 !c
1319 !c (6) ICE CRYSTAL ACCRETION BY SNOW (Psaci): Lin (22)
1320 !c
1321           esi=exp( 0.025*temcc(k) )
1322           save1=pio4*xnos*constc*gamdp3*sqrho(k)* &
1323                olambdas(k)**dp3
1324           tmp1=esi*save1
1325           psaci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) )
1326
1327 !0915     tmp1=pio4*xnos*constc*gamdp3*sqrho(k)* &
1328 !0915          olambdas(k)**dp3
1329 !0915     tmp2=qiz(k)*esi*tmp1
1330 !0915     psaci(k)=amin1( tmp2,qizodt(k) )
1331 !c
1332 !c (7) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
1333 !c
1334           esw=1.0
1335           tmp1=esw*save1
1336           psacw(k)=qlzodt(K)*( 1.0-exp(-tmp1*dtb) )
1337
1338 !0915     tmp2=qlz(k)*esw*tmp1
1339 !0915     psacw(k)=amin1( tmp2,qlzodt(k) )
1340 !c
1341 !c (8) DEPOSITION/SUBLIMATION OF SNOW (Psdep/Pssub): Lin (31)
1342 !c     includes consideration of ventilation effect
1343 !c
1344 !c     abi=2*pi*(Si-1)/rho/(A"+B")
1345 !c
1346           tmpa=rvapor*xka(k)*tem(k)*tem(k)
1347           tmpb=xls*xls*rho(k)*qsiz(k)*diffwv(k)
1348           tmpc=tmpa*qsiz(k)*diffwv(k)
1349           abi=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
1350 !
1351 !c     vf1s,vf2s=ventilation factors for snow
1352 !c     vf1s=0.78,vf2s=0.31 in LIN
1353 !
1354           tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
1355           tmp2=abi*xnos*( vf1s*olambdas(k)*olambdas(k)+ &
1356                     vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
1357           tmp3=odtb*( qvz(k)-qsiz(k) )
1358 !
1359 !bloss: BEGIN
1360 !the old implementation would give the wrong results if olambdas(k) ==0
1361 !which can lead to positive pssub, i.e. tmp2=0 but tmp3> 0
1362 !          if( tmp2 .le. 0.0) then
1363           if( tmp3 .le. 0.0) then
1364             tmp2=amax1( tmp2,tmp3)
1365             pssub(k)=amin1(0.,amax1( tmp2,-qszodt(k) ))
1366 !bloss: END
1367             psdep(k)=0.0
1368           else
1369             psdep(k)=amin1( tmp2,tmp3 )
1370             pssub(k)=0.0
1371           end if
1372
1373 !0915     psdep(k)=amax1(0.0,tmp2)
1374 !0915     pssub(k)=amin1(0.0,tmp2)
1375 !0915     pssub(k)=amax1( pssub(k),-qszodt(k) )
1376 !
1377           if(qrz(k) .le. 0.0) go to 1200
1378 !
1379 ! Compute processes (9) and (10) only when qsz > 0.0 and qrz > 0.0
1380 !
1381 !c
1382 !c (9) ACCRETION OF SNOW BY RAIN (Pracs): Lin (27)
1383 !c
1384           esr=1.0
1385           tmpa=olambdar(k)*olambdar(k)
1386           tmpb=olambdas(k)*olambdas(k)
1387           tmpc=olambdar(k)*olambdas(k)
1388           tmp1=pi*pi*esr*xnor*xnos*abs( vtr(k)-vts(k) )*orho(k)
1389           tmp2=tmpb*tmpb*olambdar(k)*(5.0*tmpb+2.0*tmpc+0.5*tmpa)
1390           tmp3=tmp1*rhosnow*tmp2
1391           pracs(k)=amin1( tmp3,qszodt(k) )
1392 !c
1393 !c (10)  ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
1394 !c
1395           tmp3=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1396           tmp4=tmp1*rhowater*tmp3
1397           psacr(k)=amin1( tmp4,qrzodt(k) )
1398 !
1399 1200      continue
1400 !
1401         else
1402 !
1403 !***********************************************************************
1404 !*********        snow production processes for T > 0 C       **********
1405 !***********************************************************************
1406 !
1407          if (qsz(k) .le. 0.0) go to 1400
1408 !c
1409 !c (1) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
1410 !c
1411             esw=1.0
1412
1413             tmp1=esw*pio4*xnos*constc*gamdp3*sqrho(k)* &
1414                  olambdas(k)**dp3
1415             psacw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
1416
1417 !0915       tmp1=pio4*xnos*constc*gamdp3*sqrho(k)* &
1418 !0915            olambdas(k)**dp3
1419 !0915       tmp2=qlz(k)*esw*tmp1
1420 !0915       psacw(k)=amin1( tmp2,qlzodt(k) )
1421 !c
1422 !c (2) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
1423 !c
1424             esr=1.0
1425             tmpa=olambdar(k)*olambdar(k)
1426             tmpb=olambdas(k)*olambdas(k)
1427             tmpc=olambdar(k)*olambdas(k)
1428             tmp1=pi*pi*esr*xnor*xnos*abs( vtr(k)-vts(k) )*orho(k)
1429             tmp2=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1430             tmp3=tmp1*rhowater*tmp2
1431             psacr(k)=amin1( tmp3,qrzodt(k) )
1432 !c
1433 !c (3) MELTING OF SNOW (Psmlt): Lin (32)
1434 !c     Psmlt is negative value
1435 !
1436             delrs=rs0(k)-qvz(k)
1437             term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- &
1438                   xka(k)*temcc(k) )
1439             tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
1440             tmp2=xnos*( vf1s*olambdas(k)*olambdas(k)+  &
1441                  vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
1442             tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( psacw(k)+psacr(k) )
1443             tmp4=amin1(0.0,tmp3)
1444             psmlt(k)=amax1( tmp4,-qszodt(k) )
1445 !c
1446 !c (4) EVAPORATION OF MELTING SNOW (Psmltevp): HR (A27)
1447 !c     but use Lin et al. coefficience
1448 !c     Psmltevp is a negative value
1449 !c
1450             tmpa=rvapor*xka(k)*tem(k)*tem(k)
1451             tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
1452             tmpc=tmpa*qswz(k)*diffwv(k)
1453             tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb )
1454
1455 !      abr=2.0*pi*(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb)
1456
1457             abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1458 !
1459 !**** allow evaporation to occur when RH less than 90%
1460 !**** here not using 100% because the evaporation cooling
1461 !**** of temperature is not taking into account yet; hence,
1462 !**** the qsw value is a little bit larger. This will avoid
1463 !**** evaporation can generate cloud.
1464 !
1465 !c    vf1s,vf2s=ventilation factors for snow
1466 !c    vf1s=0.78,vf2s=0.31 in LIN
1467 !
1468             tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
1469             tmp2=abr*xnos*( vf1s*olambdas(k)*olambdas(k)+  &
1470                  vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
1471             tmp3=amin1(0.0,tmp2)
1472             tmp3=amax1( tmp3,tmpd )
1473             psmltevp(k)=amax1( tmp3,-qszodt(k) )
1474 1400     continue
1475 !
1476         end if
1477
1478 !***********************************************************************
1479 !*********           rain production processes                **********
1480 !***********************************************************************
1481 !
1482 !c
1483 !c (1) AUTOCONVERSION OF RAIN (Praut): RH
1484 !sg: begin
1485         if(flag_qndrop)then
1486            if( qndropz(k) >= 1. ) then
1487 !         Liu et al. autoconversion scheme
1488               rhocgs=rho(k)*1.e-3
1489               liqconc=rhocgs*qlz(k)   ! (kg/kg) to (g/cm3)
1490               capn=1.0e-3*rhocgs*qndropz(k)   ! (#/kg) to (#/cm3)
1491 !         rate function
1492               if(liqconc.gt.1.e-10)then
1493                  p0=(kappa*beta/capn)*(liqconc*liqconc*liqconc)
1494                  xc=9.7d-17*capn*sqrt(capn)/(liqconc*liqconc)
1495 !         Calculate autoconversion rate (g/g/s)
1496                  if(xc.lt.10.)then
1497                     praut(k)=(p0/rhocgs) * ( 0.5d0*(xc*xc+2*xc+2.0d0)* &
1498                          (1.0d0+xc)*exp(-2.0d0*xc) )
1499                  endif
1500               endif
1501            endif
1502         else
1503 !sg: end
1504 !c          araut=afa*rho
1505 !c          afa=0.001 Rate coefficient for autoconvergence
1506 !c
1507 !c          araut=1.0e-3
1508 !c
1509             araut=0.001
1510 !testing
1511 !           tmp1=amax1( 0.0,araut*(qlz(k)-ql0) )
1512 !           praut(k)=amin1( tmp1,qlzodt(k) )
1513             tmp1=odtb*(qlz(k)-ql0)*( 1.0-exp(-araut*dtb) )
1514             praut(k)=amax1( 0.0,tmp1 )
1515         endif !sg
1516
1517 !c
1518 !c (2) ACCRETION OF CLOUD WATER BY RAIN (Pracw): Lin (51)
1519 !c
1520          erw=1.0
1521 !        tmp1=qlz(k)*pio4*erw*xnor*consta*sqrho(k)
1522 !        tmp2=tmp1*gambp3*olambdar(k)**bp3
1523 !        pracw(k)=amin1( tmp2,qlzodt(k) )
1524
1525         tmp1=pio4*erw*xnor*consta*sqrho(k)* &
1526              gambp3*olambdar(k)**bp3
1527         pracw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
1528
1529 !c
1530 !c (3) EVAPORATION OF RAIN (Prevp): Lin (52)
1531 !c     Prevp is negative value
1532 !c
1533 !c     Sw=qvoqsw : saturation ratio
1534 !c
1535          tmpa=rvapor*xka(k)*tem(k)*tem(k)
1536          tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
1537          tmpc=tmpa*qswz(k)*diffwv(k)
1538 !bloss: BEGIN
1539 !         tmpd=amin1(0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb)
1540
1541 ! set max allowed evaporation to 90% of the amount that
1542 !   would induce saturation wrt liquid in a single step.
1543          tmpd = qswz(k)*xlv/(rvapor*tem(k)**2) ! d(qsat_liq)/dT
1544          tmpd = min( 0., 0.9*odtb*(qvz(k) + qlz(k) - qswz(k)) &
1545                           / (1. + xlvocp * tmpd) )
1546
1547          abr=2.0*pi*(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb)
1548
1549 !         abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1550 !bloss: END
1551 !
1552 !c     vf1r,vf2r=ventilation factors for rain
1553 !c     vf1r=0.78,vf2r=0.31 in RH, LIN  and MM5
1554 !
1555          vf1r=0.78
1556          vf2r=0.31
1557          tmp1=consta*sqrho(k)*olambdar(k)**bp5/visc(k)
1558          tmp2=abr*xnor*( vf1r*olambdar(k)*olambdar(k)+  &
1559               vf2r*schmidt(k)**0.33334*gambp5o2*sqrt(tmp1) )
1560          tmp3=amin1( 0.0,tmp2 )
1561          tmp3=amax1( tmp3,tmpd )
1562          prevp(k)=amax1( tmp3,-qrzodt(k) )
1563
1564 !
1565 !      if(iout .gt. 0) write(20,*)'tmp1,tmp2,tmp3=',tmp1,tmp2,tmp3
1566 !      if(iout .gt. 0) write(20,*)'qlz,qiz,qrz=',qlz(k),qiz(k),qrz(k)
1567 !      if(iout .gt. 0) write(20,*)'tem,qsz,qvz=',tem(k),qsz(k),qvz(k)
1568
1569
1570
1571 !     if (gindex .eq. 0.) goto 900
1572 !
1573       if (tem(k) .lt. 273.15) then
1574 !
1575 !
1576 !-- graupel
1577 !***********************************************************************
1578 !*********        graupel production processes for T < 0 C    **********
1579 !***********************************************************************
1580 !c
1581 !c (1) AUTOCONVERSION OF SNOW TO FORM GRAUPEL (Pgaut): Lin (37)
1582 !c     pgaut=alpha2*(qsz-qs0)
1583 !c     qs0=6.0E-4
1584 !c     alpha2=1.0e-3*exp(0.09*temcc(k))      Lin (38)
1585 !
1586             alpha2=1.0e-3*exp(0.09*temcc(k))
1587 !
1588
1589 ! testing
1590 !           tmp1=alpha2*(qsz(k)-qs0)
1591 !           tmp1=amax1(0.0,tmp1)
1592 !           pgaut(k)=amin1( tmp1,qszodt(k) )
1593
1594             tmp1=odtb*(qsz(k)-qs0)*(1.0-exp(-alpha2*dtb))
1595             pgaut(k)=amax1( 0.0,tmp1 )
1596
1597 !c
1598 !c (2) FREEZING OF RAIN TO FORM GRAUPEL  (Pgfr): Lin (45)
1599 !c     positive value
1600 !c     Constant in Bigg freezing Aplume=Ap=0.66 /k
1601 !c     Constant in raindrop freezing equ. Bplume=Bp=100./m/m/m/s
1602 !
1603
1604             if (qrz(k) .gt. 1.e-8 ) then
1605                Bp=100.
1606                Ap=0.66
1607                tmp1=olambdar(k)*olambdar(k)*olambdar(k)
1608                tmp2=20.*pi*pi*Bp*xnor*rhowater*orho(k)*  &
1609                     (exp(-Ap*temcc(k))-1.0)*tmp1*tmp1*olambdar(k)
1610                Pgfr(k)=amin1( tmp2,qrzodt(k) )
1611             else
1612                Pgfr(k)=0
1613             endif
1614
1615 !c
1616 !c       if (qgz(k) = 0.0) skip the other step below about graupel
1617 !c
1618          if (qgz(k) .eq. 0.0) goto 4000
1619
1620 !c
1621 !c       Comparing Pgwet(wet process) and Pdry(dry process),
1622 !c       we will pick up the small one.
1623 !c
1624
1625 !c       ---------------
1626 !c      | dry processes |
1627 !c       ---------------
1628 !c
1629 !c (3)   ACCRETION OF CLOUD WATER BY GRAUPEL  (Pgacw): Lin (40)
1630 !c       egw=1.0
1631 !c       Cdrag=0.6 drag coefficients for hairstone
1632 !c       constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1633 !c
1634          egw=1.0
1635          constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag)
1636          tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5
1637          tmp2=qlz(k)*egw*tmp1
1638          Pgacw(k)=amin1( tmp2,qlzodt(k) )
1639 !c
1640 !c (4)   ACCRETION OF ICE CRYSTAL BY GRAUPEL  (Pgaci): Lin (41)
1641 !c       egi=1.   for wet growth
1642 !c       egi=0.1  for dry growth
1643 !c
1644          egi=0.1
1645          tmp2=qiz(k)*egi*tmp1
1646          pgaci(k)=amin1( tmp2,qizodt(k) )
1647
1648
1649 !c
1650 !c (5)   ACCRETION OF SNOW BY GRAUPEL  (Pgacs) : Lin (29)
1651 !c       Compute processes (6) only when qsz > 0.0 and qgz > 0.0
1652 !c
1653          egs=exp(0.09*temcc(k))
1654          tmpa=olambdas(k)*olambdas(k)
1655          tmpb=olambdag(k)*olambdag(k)
1656          tmpc=olambdas(k)*olambdag(k)
1657          tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k)
1658          tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1659          tmp3=tmp1*egs*rhosnow*tmp2
1660          Pgacs(k)=amin1( tmp3,qszodt(k) )
1661
1662
1663 !c
1664 !c (6)   ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
1665 !c       Compute processes (6) only when qrz > 0.0 and qgz > 0.0
1666 !c       egr=1.
1667 !c
1668          egr=1.
1669          tmpa=olambdar(k)*olambdar(k)
1670          tmpb=olambdag(k)*olambdag(k)
1671          tmpc=olambdar(k)*olambdag(k)
1672          tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k)
1673          tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1674          tmp3=tmp1*egr*rhowater*tmp2
1675          pgacr(k)=amin1( tmp3,qrzodt(k) )
1676
1677 !c
1678 !c (7)   Calculate total dry process effect Pdry(k)
1679 !c
1680          Pdry(k)=Pgacw(k)+pgaci(k)+Pgacs(k)+pgacr(k)
1681
1682 !c       ---------------
1683 !c      | wet processes |
1684 !c       ---------------
1685 !c
1686 !c (3)   ACCRETION OF ICE CRYSTAL BY GRAUPEL  (Pgacip): Lin (41)
1687 !c       egi=1.   for wet growth
1688 !c       egi=0.1  for dry growth
1689 !c
1690          tmp2=10.*pgaci(k)
1691          pgacip(k)=amin1( tmp2,qizodt(k) )
1692
1693 !c
1694 !c (4)   ACCRETION OF SNOW BY GRAUPEL  ((Pgacsp) : Lin (29)
1695 !c       Compute processes (6) only when qsz > 0.0 and qgz > 0.0
1696 !c       egs=exp(0.09*(tem(k)-273.15))  when T < 273.15 k
1697 !c
1698          tmp3=Pgacs(k)*1.0/egs
1699          Pgacsp(k)=amin1( tmp3,qszodt(k) )
1700
1701 !c
1702 !c (5)   WET GROWTH OF GRAUPEL (Pgwet) : Lin (43)
1703 !c       may involve Pgacs or Pgaci and
1704 !c       must include PPgacw or Pgacr, or both.
1705 !c       ( The amount of Pgacw which is not able
1706 !c       to freeze is shed to rain. )
1707          IF(temcc(k).gt.-40.)THEN
1708
1709              term0=constg*olambdag(k)**5.5/visc(k)
1710
1711 !c
1712 !c    vf1s,vf2s=ventilation factors for graupel
1713 !c    vf1s=0.78,vf2s=0.31 in LIN
1714 !c    Cdrag=0.6  drag coefficient for hairstone
1715 !c    constg2=vf1s*olambdag(k)*olambdag(k)+
1716 !c            vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
1717
1718              delrs=rs0(k)-qvz(k)
1719              tmp0=1./(xlf+cw*temcc(k))
1720              tmp1=2.*pi*xnog*(rho(k)*xlv*diffwv(k)*delrs-xka(k)*  &
1721                   temcc(k))*orho(k)*tmp0
1722              constg2=vf1s*olambdag(k)*olambdag(k)+  &
1723                      vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
1724              tmp3=tmp1*constg2+(Pgacip(k)+Pgacsp(k))*  &
1725                   (1-Ci*temcc(k)*tmp0)
1726              tmp3=amax1(0.0,tmp3)
1727              Pgwet(k)=amin1(tmp3,qlzodt(k)+qszodt(k)+qizodt(k) ) !bloss
1728
1729 !c
1730 !c     Comparing Pgwet(wet process) and Pdry(dry process),
1731 !c     we will apply the small one.
1732 !c     if dry processes then delta4=1.0
1733 !c     if wet processes then delta4=0.0
1734 !
1735          if ( Pdry(k) .lt. Pgwet(k) ) then
1736             delta4=1.0
1737          else
1738             delta4=0.0
1739          endif
1740          ELSE
1741             delta4=1.0
1742          ENDIF
1743
1744 !c
1745 !c
1746 !c (6)   Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k)
1747 !c       if Pgacrp(k) > 0. then some of the rain is frozen to hail
1748 !c       if Pgacrp(k) < 0. then some of the cloud water collected
1749 !c                            by the hail is unable to freeze and is
1750 !c                            shed as rain.
1751 !c
1752             Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k)
1753
1754 !c
1755 !c (8)   DEPOSITION/SUBLIMATION OF GRAUPEL  (Pgdep/Pgsub): Lin (46)
1756 !c       includes ventilation effect
1757 !c       constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1758 !c       constg2=vf1s*olambdag(k)*olambdag(k)+
1759 !c             vf2s*schmidt(k)**0.33334*gam2pt75*constg
1760 !c
1761 !c       abg=2*pi*(Si-1)/rho/(A"+B")
1762 !c
1763             tmpa=rvapor*xka(k)*tem(k)*tem(k)
1764             tmpb=xls*xls*rho(k)*qsiz(k)*diffwv(k)
1765             tmpc=tmpa*qsiz(k)*diffwv(k)
1766             abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
1767 !c
1768 !c     vf1s,vf2s=ventilation factors for graupel
1769 !c     vf1s=0.78,vf2s=0.31 in LIN
1770 !c     Cdrag=0.6  drag coefficient for hairstone
1771 !c
1772             term0=constg*olambdag(k)**5.5/visc(k)
1773             constg2=vf1s*olambdag(k)*olambdag(k)+  &
1774                     vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
1775             tmp2=abg*xnog*constg2
1776             pgdep(k)=amax1(0.0,tmp2)
1777             pgsub(k)=amin1(0.0,tmp2)
1778             pgsub(k)=amax1( pgsub(k),-qgzodt(k) )
1779
1780  4000    continue
1781         else
1782 !
1783 !***********************************************************************
1784 !*********      graupel production processes for T > 0 C      **********
1785 !***********************************************************************
1786 !
1787 !c
1788 !c (1) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40)
1789 !c     egw=1.0
1790 !c     Cdrag=0.6 drag coefficients for hairstone
1791 !c     constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1792
1793             egw=1.0
1794             constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag)
1795             tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5
1796             tmp2=qlz(k)*egw*tmp1
1797             Pgacw(k)=amin1( tmp2,qlzodt(k) )
1798
1799 !c
1800 !c (2) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
1801 !c     Compute processes (5) only when qrz > 0.0 and qgz > 0.0
1802 !c     egr=1.
1803 !c
1804             egr=1.
1805             tmpa=olambdar(k)*olambdar(k)
1806             tmpb=olambdag(k)*olambdag(k)
1807             tmpc=olambdar(k)*olambdag(k)
1808             tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k)
1809             tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1810             tmp3=tmp1*egr*rhowater*tmp2
1811             pgacr(k)=amin1( tmp3,qrzodt(k) )
1812
1813
1814 !c
1815 !c (3) GRAUPEL MELTING TO FORM RAIN (Pgmlt): Lin (47)
1816 !c     Pgmlt is negative value
1817 !c     constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1818 !c     constg2=vf1s*olambdag(k)*olambdag(k)+
1819 !c             vf2s*schmidt(k)**0.33334*gam2pt75*constg
1820 !c     Cdrag=0.6  drag coefficients for hairstone
1821 !
1822             delrs=rs0(k)-qvz(k)
1823             term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- &
1824                   xka(k)*temcc(k) )
1825             term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag) &
1826                   *olambdag(k)**5.5/visc(k)
1827
1828             constg2=vf1s*olambdag(k)*olambdag(k)+ &
1829                     vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
1830             tmp2=xnog*constg2
1831             tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( pgacw(k)+pgacr(k) )
1832             tmp4=amin1(0.0,tmp3)
1833             pgmlt(k)=amax1( tmp4,-qgzodt(k) )
1834
1835
1836 !c
1837 !c (4) EVAPORATION OF MELTING GRAUPEL (Pgmltevp) : HR (A19)
1838 !c     but use Lin et al. coefficience
1839 !c     Pgmltevp is a negative value
1840 !c     abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
1841 !c
1842             tmpa=rvapor*xka(k)*tem(k)*tem(k)
1843             tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
1844             tmpc=tmpa*qswz(k)*diffwv(k)
1845             tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb )
1846
1847 !c
1848 !c     abg=2*pi*(Si-1)/rho/(A"+B")
1849 !c
1850             abg=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1851 !
1852 !**** allow evaporation to occur when RH less than 90%
1853 !**** here not using 100% because the evaporation cooling
1854 !**** of temperature is not taking into account yet; hence,
1855 !**** the qgw value is a little bit larger. This will avoid
1856 !**** evaporation can generate cloud.
1857 !
1858 !c    vf1s,vf2s=ventilation factors for snow
1859 !c    vf1s=0.78,vf2s=0.31 in LIN
1860 !c    constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1861 !c    constg2=vf1s*olambdag(k)*olambdag(k)+
1862 !c            vf2s*schmidt(k)**0.33334*gam2pt75*constg
1863 !
1864             tmp2=abg*xnog*constg2
1865             tmp3=amin1(0.0,tmp2)
1866             tmp3=amax1( tmp3,tmpd )
1867             pgmltevp(k)=amax1( tmp3,-qgzodt(k) )
1868
1869 !c
1870 !c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29)
1871 !c     Compute processes (3) only when qsz > 0.0 and qgz > 0.0
1872 !c     egs=1.0
1873 !c
1874            egs=1.
1875            tmpa=olambdas(k)*olambdas(k)
1876            tmpb=olambdag(k)*olambdag(k)
1877            tmpc=olambdas(k)*olambdag(k)
1878            tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k)
1879            tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1880            tmp3=tmp1*egs*rhosnow*tmp2
1881            Pgacs(k)=amin1( tmp3,qszodt(k) )
1882
1883         endif
1884
1885
1886 !
1887   900   continue
1888
1889 !cc
1890 !c
1891 !c**********************************************************************
1892 !c*****     combine all processes together and avoid negative      *****
1893 !c*****     water substances
1894 !***********************************************************************
1895 !c
1896       if ( temcc(k) .lt. 0.0) then
1897 !,delta4,1.-delta4
1898 !c
1899 !c  gdelta4=gindex*delta4
1900 !c  g1sdelt4=gindex*(1.-delta4)
1901 !c
1902            gdelta4=gindex*delta4
1903            g1sdelt4=gindex*(1.-delta4)
1904 !c
1905 !c  combined water vapor depletions
1906 !c
1907 !cc graupel
1908            tmp=psdep(k)+pgdep(k)*gindex
1909            if ( tmp .gt. qvzodt(k) ) then
1910               factor=qvzodt(k)/tmp
1911               psdep(k)=psdep(k)*factor
1912               pgdep(k)=pgdep(k)*factor*gindex
1913            end if
1914 !c
1915 !c  combined cloud water depletions
1916 !c
1917            tmp=praut(k)+psacw(k)+psfw(k)+pracw(k)+gindex*Pgacw(k)
1918            if ( tmp .gt. qlzodt(k) ) then
1919               factor=qlzodt(k)/tmp
1920               praut(k)=praut(k)*factor
1921               psacw(k)=psacw(k)*factor
1922               psfw(k)=psfw(k)*factor
1923               pracw(k)=pracw(k)*factor
1924 !cc graupel
1925               Pgacw(k)=Pgacw(k)*factor*gindex
1926            end if
1927 !c
1928 !c  combined cloud ice depletions
1929 !c
1930            tmp=psaut(k)+psaci(k)+praci(k)+psfi(k)+Pgaci(k)*gdelta4 &
1931                +Pgacip(k)*g1sdelt4
1932            if (tmp .gt. qizodt(k) ) then
1933               factor=qizodt(k)/tmp
1934               psaut(k)=psaut(k)*factor
1935               psaci(k)=psaci(k)*factor
1936               praci(k)=praci(k)*factor
1937               psfi(k)=psfi(k)*factor
1938 !cc graupel
1939               Pgaci(k)=Pgaci(k)*factor*gdelta4
1940               Pgacip(k)=Pgacip(k)*factor*g1sdelt4
1941            endif
1942 !c
1943 !c  combined all rain processes
1944 !c
1945           tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
1946                 +Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
1947                 +Pgacrp(k)*g1sdelt4
1948           if (tmp_r .gt. qrzodt(k) ) then
1949              factor=qrzodt(k)/tmp_r
1950              piacr(k)=piacr(k)*factor
1951              psacr(k)=psacr(k)*factor
1952              prevp(k)=prevp(k)*factor
1953 !cc graupel
1954              Pgfr(k)=Pgfr(k)*factor*gindex
1955              Pgacr(k)=Pgacr(k)*factor*gdelta4
1956              Pgacrp(k)=Pgacrp(k)*factor*g1sdelt4
1957           endif
1958
1959 !c
1960 !c     if qrz < 1.0E-4 and qsz < 1.0E-4  then delta2=1.
1961 !c                  (all Pracs and Psacr become to snow)
1962 !c     if qrz >= 1.0E-4 or qsz >= 1.0E-4 then delta2=0.
1963 !c                 (all Pracs and Psacr become to graupel)
1964 !c
1965           if (qrz(k) .lt. 1.0E-4 .and. qsz(k) .lt. 1.0E-4) then
1966              delta2=1.0
1967           else
1968              delta2=0.0
1969           endif
1970 !
1971 !cc graupel
1972
1973 !c
1974 !c  if qrz(k) < 1.0e-4 then delta3=1. means praci(k) -->  qs
1975 !c                                          piacr(k) -->  qs
1976 !c  if qrz(k) > 1.0e-4 then delta3=0. means praci(k) -->  qg
1977 !c                                          piacr(k) -->  qg : Lin (20)
1978
1979           if (qrz(k) .lt. 1.0e-4) then
1980              delta3=1.0
1981           else
1982              delta3=0.0
1983           endif
1984 !
1985 !c
1986 !c     if gindex = 0.(no graupel) then delta2=1.0
1987 !c                                     delta3=1.0
1988 !c
1989           if (gindex .eq. 0.) then
1990               delta2=1.0
1991               delta3=1.0
1992          endif
1993 !
1994 !c
1995 !c   combined all snow processes
1996 !c
1997           tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
1998                  psfi(k)+praci(k)*delta3+piacr(k)*delta3+ &
1999                  psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ &
2000                  Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- &
2001                  Psacr(k)*delta2
2002           if ( tmp_s .gt. qszodt(k) ) then
2003              factor=qszodt(k)/tmp_s
2004              pssub(k)=pssub(k)*factor
2005              Pracs(k)=Pracs(k)*factor
2006 !cc graupel
2007              Pgaut(k)=Pgaut(k)*factor*gindex
2008              Pgacs(k)=Pgacs(k)*factor*gdelta4
2009              Pgacsp(k)=Pgacsp(k)*factor*g1sdelt4
2010           endif
2011
2012 !cc graupel
2013 !
2014
2015 !          if (gindex .eq. 0.) goto 998
2016 !c
2017 !c  combined all graupel processes
2018 !c
2019            if(delta4.lt.0.5) then
2020              !Re-define pgwet to account for limiting of pgacrp,
2021              !         pgacip, pgacw and pgacsp above
2022              pgwet(k) = pgacrp(k) + pgacw(k) + pgacip(k) + pgacsp(k)
2023            end if
2024            tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4  &
2025                  -Pgacr(k)*delta4-Pgacs(k)*delta4  &
2026                  -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k)  &
2027                  -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2)  &
2028                  -praci(k)*(1-delta3)-piacr(k)*(1-delta3)
2029            if (tmp_g .gt. qgzodt(k)) then
2030                factor=qgzodt(k)/tmp_g
2031                pgsub(k)=pgsub(k)*factor
2032            endif
2033
2034   998      continue
2035 !c
2036 !c  calculate new water substances, thetae, tem, and qvsbar
2037 !c
2038
2039 !cc graupel
2040          pvapor(k)=-pssub(k)-psdep(k)-prevp(k)-pgsub(k)*gindex &
2041                    -pgdep(k)*gindex
2042          qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k) )
2043          pclw(k)=-praut(k)-pracw(k)-psacw(k)-psfw(k)-pgacw(k)*gindex
2044          if(flag_qndrop)then
2045            if( qlz(k) > 1e-20 ) &
2046               qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) )  !sg
2047          endif
2048          qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
2049          pcli(k)=-psaut(k)-psfi(k)-psaci(k)-praci(k)-pgaci(k)*gdelta4 &
2050                  -Pgacip(k)*g1sdelt4
2051          qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
2052          tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
2053                 +Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
2054                 +Pgacrp(k)*g1sdelt4
2055  232     format(i2,1x,6(f9.3,1x))
2056          prain(k)=-tmp_r
2057          qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
2058          tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+  &
2059                 psfi(k)+praci(k)*delta3+piacr(k)*delta3+  &
2060                 psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+  &
2061                 Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)-  &
2062                 Psacr(k)*delta2
2063          psnow(k)=-tmp_s
2064          qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
2065          qschg(k)=qschg(k)+psnow(k)
2066          qschg(k)=psnow(k)
2067 !cc graupel
2068          tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 &
2069                -Pgacr(k)*delta4-Pgacs(k)*delta4 &
2070                -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
2071                -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) &
2072                -praci(k)*(1-delta3)-piacr(k)*(1-delta3)
2073  252     format(i2,1x,6(f12.9,1x))
2074  262     format(i2,1x,7(f12.9,1x))
2075          pgraupel(k)=-tmp_g
2076          pgraupel(k)=pgraupel(k)*gindex
2077          qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
2078 !        qgchg(k)=qgchg(k)+pgraupel(k)
2079          qgchg(k)=pgraupel(k)
2080          qgz(k)=qgz(k)*gindex
2081
2082          tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
2083          theiz(k)=theiz(k)+dtb*tmp
2084          thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2085          tem(k)=thz(k)*tothz(k)
2086
2087          temcc(k)=tem(k)-273.15
2088
2089 !bloss: Moves computation of saturation mixing ratio below
2090
2091 !         if( temcc(k) .lt. -40.0 ) qswz(k)=qsiz(k)
2092 !         qlpqi=qlz(k)+qiz(k)
2093 !         if ( qlpqi .eq. 0.0 ) then
2094 !            qvsbar(k)=qsiz(k)
2095 !         else
2096 !            qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2097 !         endif
2098
2099 !
2100       else
2101 !c
2102 !c  combined cloud water depletions
2103 !c
2104           tmp=praut(k)+psacw(k)+pracw(k)+pgacw(k)*gindex
2105           if ( tmp .gt. qlzodt(k) ) then
2106              factor=qlzodt(k)/tmp
2107              praut(k)=praut(k)*factor
2108              psacw(k)=psacw(k)*factor
2109              pracw(k)=pracw(k)*factor
2110 !cc graupel
2111              pgacw(k)=pgacw(k)*factor*gindex
2112           end if
2113 !c
2114 !c  combined all snow processes
2115 !c
2116           tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
2117           if (tmp_s .gt. qszodt(k) ) then
2118              factor=qszodt(k)/tmp_s
2119              psmlt(k)=psmlt(k)*factor
2120              psmltevp(k)=psmltevp(k)*factor
2121 !cc graupel
2122              Pgacs(k)=Pgacs(k)*factor*gindex
2123           endif
2124
2125 !c
2126 !c
2127 !cc graupel
2128 !c
2129 !         if (gindex .eq. 0.) goto 997
2130
2131 !c
2132 !c  combined all graupel processes
2133 !c
2134             tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
2135             if (tmp_g .gt. qgzodt(k)) then
2136               factor=qgzodt(k)/tmp_g
2137               pgmltevp(k)=pgmltevp(k)*factor
2138               pgmlt(k)=pgmlt(k)*factor
2139             endif
2140 !c
2141   997     continue
2142
2143 !c
2144 !c  combined all rain processes
2145 !c
2146           tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
2147                 +pgmlt(k)*gindex-pgacw(k)*gindex
2148           if (tmp_r .gt. qrzodt(k) ) then
2149              factor=qrzodt(k)/tmp_r
2150              prevp(k)=prevp(k)*factor
2151           endif
2152 !c
2153 !c
2154 !c  calculate new water substances and thetae
2155 !c
2156
2157
2158           pvapor(k)=-psmltevp(k)-prevp(k)-pgmltevp(k)*gindex
2159           qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k))
2160           pclw(k)=-praut(k)-pracw(k)-psacw(k)-pgacw(k)*gindex
2161           if(flag_qndrop)then
2162              if( qlz(k) > 1e-20 ) &
2163                qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) )  !sg
2164           endif
2165           qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
2166           pcli(k)=0.0
2167           qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
2168           tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
2169                 +pgmlt(k)*gindex-pgacw(k)*gindex
2170  242      format(i2,1x,7(f9.6,1x))
2171           prain(k)=-tmp_r
2172           tmpqrz=qrz(k)
2173           qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
2174           tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
2175           psnow(k)=-tmp_s
2176           qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
2177 !         qschg(k)=qschg(k)+psnow(k)
2178           qschg(k)=psnow(k)
2179 !cc graupel
2180
2181           tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
2182 !         write(*,272)k,pgmlt(k),pgacs(k),pgmltevp(k),
2183  272      format(i2,1x,3(f12.9,1x))
2184           pgraupel(k)=-tmp_g*gindex
2185           qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
2186 !         qgchg(k)=qgchg(k)+pgraupel(k)
2187           qgchg(k)=pgraupel(k)
2188           qgz(k)=qgz(k)*gindex
2189 !
2190           tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
2191           theiz(k)=theiz(k)+dtb*tmp
2192           thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2193
2194           tem(k)=thz(k)*tothz(k)
2195           temcc(k)=tem(k)-273.15
2196 !         qswz(k)=episp0k*oprez(k)*  &
2197 !                exp( svp2*temcc(k)/(tem(k)-svp3) )
2198
2199 !bloss: Moved computation of saturation mixing ratio below
2200
2201 !          es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2202 !          qswz(k)=ep2*es/(prez(k)-es)
2203 !          qsiz(k)=qswz(k)
2204 !          qvsbar(k)=qswz(k)
2205 !
2206       end if
2207       preclw(k)=pclw(k)        !sg
2208
2209 !
2210 !***********************************************************************
2211 !**********              saturation adjustment                **********
2212 !***********************************************************************
2213 !bloss: BEGIN
2214 !
2215 !     compute saturation specific humidity at the temperature
2216 !     which would be experienced if all cloud liquid and ice
2217 !     were to become vapor.  This will make for a consistent
2218 !     check of saturation.  Previously, qvsbar was computed
2219 !     without accounting for the change in temperature due to
2220 !     evaporation/sublimation, and as a result would
2221 !     incorrectly identify some points as subsaturated.
2222
2223       tmp_tem = tem(k)          ! updated temperature from above.
2224
2225       if(qlz(k)+qiz(k).gt. 0.) then
2226          ! account for latent heat if all qlz and qiz were converted to vapor
2227          tmp_tem = tem(k) - xlvocp*qlz(k) - (xlvocp+xlfocp)*qiz(k)
2228       end if
2229
2230      ! compute temperature in celsius
2231       tmp_temcc = tmp_tem - 273.15
2232
2233       ! estimate lower bound on saturation vapor pressure
2234       !   (ice if <0C, liquid aboce)
2235       if (tmp_temcc .lt. 0.) then
2236          es=1000.*svp1*exp( 21.8745584*(tmp_tem-273.16)/(tmp_tem-7.66) ) !ice
2237       else
2238          es=1000.*svp1*exp( svp2*tmp_temcc/(tmp_tem-svp3) ) !liquid
2239       end if
2240
2241       ! compute lower bound on saturation mixing ratio.
2242       qvsbar(k)=ep2*es/(prez(k)-es)
2243 !
2244 !bloss: END
2245 !
2246 !    allow supersaturation exits linearly from 0% at 500 mb to 50%
2247 !    above 300 mb
2248 !    5.0e-5=1.0/(500mb-300mb)
2249 !
2250          rsat=1.0+0.5*(50000.0-prez(k))*5.0e-5
2251          rsat=amax1(1.0,rsat)
2252          rsat=amin1(1.5,rsat)
2253          rsat=1.0
2254          if( qvz(k)+qlz(k)+qiz(k) .lt. rsat*qvsbar(k) ) then
2255
2256 !c
2257 !c   unsaturated
2258 !c
2259           qvz(k)=qvz(k)+qlz(k)+qiz(k)
2260           qlz(k)=0.0
2261           qiz(k)=0.0
2262
2263           thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2264           tem(k)=thz(k)*tothz(k)
2265           temcc(k)=tem(k)-273.15
2266
2267           go to 1800
2268 !
2269         else
2270 !c
2271 !c   saturated
2272 !c
2273 !
2274           pladj(k)=qlz(k)
2275           piadj(k)=qiz(k)
2276 !
2277
2278           CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
2279                       k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0 )
2280
2281 !
2282           pladj(k)=odtb*(qlz(k)-pladj(k))
2283           piadj(k)=odtb*(qiz(k)-piadj(k))
2284 !
2285           pclw(k)=pclw(k)+pladj(k)
2286           pcli(k)=pcli(k)+piadj(k)
2287           pvapor(k)=pvapor(k)-( pladj(k)+piadj(k) )
2288 !
2289           thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2290           tem(k)=thz(k)*tothz(k)
2291
2292           temcc(k)=tem(k)-273.15
2293
2294 !         qswz(k)=episp0k*oprez(k)*  &
2295 !                 exp( svp2*temcc(k)/(tem(k)-svp3) )
2296           es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2297           qswz(k)=ep2*es/(prez(k)-es)
2298           if (tem(k) .lt. 273.15 ) then
2299 !            qsiz(k)=episp0k*oprez(k)* &
2300 !                   exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2301              es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2302              qsiz(k)=ep2*es/(prez(k)-es)
2303              if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
2304           else
2305              qsiz(k)=qswz(k)
2306           endif
2307           qlpqi=qlz(k)+qiz(k)
2308           if ( qlpqi .eq. 0.0 ) then
2309              qvsbar(k)=qsiz(k)
2310           else
2311              qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2312           endif
2313
2314         end if
2315
2316 !
2317 !***********************************************************************
2318 !*****     melting and freezing of cloud ice and cloud water       *****
2319 !***********************************************************************
2320         qlpqi=qlz(k)+qiz(k)
2321         if(qlpqi .le. 0.0) go to 1800
2322 !
2323 !c
2324 !c (1)  HOMOGENEOUS NUCLEATION WHEN T< -40 C (Pihom)
2325 !c
2326         if(temcc(k) .lt. -40.0) pihom(k)=qlz(k)*odtb
2327 !c
2328 !c (2)  MELTING OF ICE CRYSTAL WHEN T> 0 C (Pimlt)
2329 !c
2330         if(temcc(k) .gt. 0.0) pimlt(k)=qiz(k)*odtb
2331 !c
2332 !c (3) PRODUCTION OF CLOUD ICE BY BERGERON PROCESS (Pidw): Hsie (p957)
2333 !c     this process only considered when -31 C < T < 0 C
2334 !c
2335         if(temcc(k) .lt. 0.0 .and. temcc(k) .gt. -31.0) then
2336 !c!
2337 !c!   parama1 and parama2 functions must be user supplied
2338 !c!
2339           a1=parama1( temcc(k) )
2340           a2=parama2( temcc(k) )
2341 !! change unit from cgs to mks
2342           a1=a1*0.001**(1.0-a2)
2343           xnin=xni0*exp(-bni*temcc(k))
2344           pidw(k)=xnin*orho(k)*(a1*xmnin**a2)
2345         end if
2346 !
2347         pcli(k)=pcli(k)+pihom(k)-pimlt(k)+pidw(k)
2348         pclw(k)=pclw(k)-pihom(k)+pimlt(k)-pidw(k)
2349         qlz(k)=amax1( 0.0,qlz(k)+dtb*(-pihom(k)+pimlt(k)-pidw(k)) )
2350         qiz(k)=amax1( 0.0,qiz(k)+dtb*(pihom(k)-pimlt(k)+pidw(k)) )
2351
2352 !
2353         CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
2354                     k, xLvocp, xLfocp, episp0k ,EP2,SVP1,SVP2,SVP3,SVPT0)
2355
2356         thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2357         tem(k)=thz(k)*tothz(k)
2358
2359         temcc(k)=tem(k)-273.15
2360
2361 !       qswz(k)=episp0k*oprez(k)* &
2362 !              exp( svp2*temcc(k)/(tem(k)-svp3) )
2363         es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2364         qswz(k)=ep2*es/(prez(k)-es)
2365
2366         if (tem(k) .lt. 273.15 ) then
2367 !          qsiz(k)=episp0k*oprez(k)* &
2368 !                 exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2369            es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2370            qsiz(k)=ep2*es/(prez(k)-es)
2371            if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
2372         else
2373            qsiz(k)=qswz(k)
2374         endif
2375         qlpqi=qlz(k)+qiz(k)
2376         if ( qlpqi .eq. 0.0 ) then
2377            qvsbar(k)=qsiz(k)
2378         else
2379            qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2380         endif
2381
2382 1800  continue
2383 !
2384 !***********************************************************************
2385 !**********    integrate the productions of rain and snow     **********
2386 !***********************************************************************
2387 !c
2388
2389 2000  continue
2390
2391
2392 !---------------------------------------------------------------------
2393
2394 !
2395 !***********************************************************************
2396 !******      Write terms in cloud physics to time series dataset   *****
2397 !***********************************************************************
2398 !
2399 !     open(unit=24,form='formatted',status='new',
2400 !    &     file='cloud.dat')
2401
2402 !9030  format(10e12.6)
2403
2404 !      write(24,*)'tmp'
2405 !      write(24,9030) (tem(k),k=kts+1,kte)
2406 !      write(24,*)'qiz'
2407 !      write(24,9030) (qiz(k),k=kts+1,kte)
2408 !      write(24,*)'qsz'
2409 !      write(24,9030) (qsz(k),k=kts+1,kte)
2410 !      write(24,*)'qrz'
2411 !      write(24,9030) (qrz(k),k=kts+1,kte)
2412 !      write(24,*)'qgz'
2413 !      write(24,9030) (qgz(k),k=kts+1,kte)
2414 !      write(24,*)'qvoqsw'
2415 !      write(24,9030) (qvoqswz(k),k=kts+1,kte)
2416 !      write(24,*)'qvoqsi'
2417 !      write(24,9030) (qvoqsiz(k),k=kts+1,kte)
2418 !      write(24,*)'vtr'
2419 !      write(24,9030) (vtr(k),k=kts+1,kte)
2420 !      write(24,*)'vts'
2421 !      write(24,9030) (vts(k),k=kts+1,kte)
2422 !      write(24,*)'vtg'
2423 !      write(24,9030) (vtg(k),k=kts+1,kte)
2424 !      write(24,*)'pclw'
2425 !      write(24,9030) (pclw(k),k=kts+1,kte)
2426 !      write(24,*)'pvapor'
2427 !      write(24,9030) (pvapor(k),k=kts+1,kte)
2428 !      write(24,*)'pcli'
2429 !      write(24,9030) (pcli(k),k=kts+1,kte)
2430 !      write(24,*)'pimlt'
2431 !      write(24,9030) (pimlt(k),k=kts+1,kte)
2432 !      write(24,*)'pihom'
2433 !      write(24,9030) (pihom(k),k=kts+1,kte)
2434 !      write(24,*)'pidw'
2435 !      write(24,9030) (pidw(k),k=kts+1,kte)
2436 !      write(24,*)'prain'
2437 !      write(24,9030) (prain(k),k=kts+1,kte)
2438 !      write(24,*)'praut'
2439 !      write(24,9030) (praut(k),k=kts+1,kte)
2440 !      write(24,*)'pracw'
2441 !      write(24,9030) (pracw(k),k=kts+1,kte)
2442 !      write(24,*)'prevp'
2443 !      write(24,9030) (prevp(k),k=kts+1,kte)
2444 !      write(24,*)'psnow'
2445 !      write(24,9030) (psnow(k),k=kts+1,kte)
2446 !      write(24,*)'psaut'
2447 !      write(24,9030) (psaut(k),k=kts+1,kte)
2448 !      write(24,*)'psfw'
2449 !      write(24,9030) (psfw(k),k=kts+1,kte)
2450 !      write(24,*)'psfi'
2451 !      write(24,9030) (psfi(k),k=kts+1,kte)
2452 !      write(24,*)'praci'
2453 !      write(24,9030) (praci(k),k=kts+1,kte)
2454 !      write(24,*)'piacr'
2455 !      write(24,9030) (piacr(k),k=kts+1,kte)
2456 !      write(24,*)'psaci'
2457 !      write(24,9030) (psaci(k),k=kts+1,kte)
2458 !      write(24,*)'psacw'
2459 !      write(24,9030) (psacw(k),k=kts+1,kte)
2460 !      write(24,*)'psdep'
2461 !      write(24,9030) (psdep(k),k=kts+1,kte)
2462 !      write(24,*)'pssub'
2463 !      write(24,9030) (pssub(k),k=kts+1,kte)
2464 !      write(24,*)'pracs'
2465 !      write(24,9030) (pracs(k),k=kts+1,kte)
2466 !      write(24,*)'psacr'
2467 !      write(24,9030) (psacr(k),k=kts+1,kte)
2468 !      write(24,*)'psmlt'
2469 !      write(24,9030) (psmlt(k),k=kts+1,kte)
2470 !      write(24,*)'psmltevp'
2471 !      write(24,9030) (psmltevp(k),k=kts+1,kte)
2472 !      write(24,*)'pladj'
2473 !      write(24,9030) (pladj(k),k=kts+1,kte)
2474 !      write(24,*)'piadj'
2475 !      write(24,9030) (piadj(k),k=kts+1,kte)
2476 !      write(24,*)'pgraupel'
2477 !      write(24,9030) (pgraupel(k),k=kts+1,kte)
2478 !      write(24,*)'pgaut'
2479 !      write(24,9030) (pgaut(k),k=kts+1,kte)
2480 !      write(24,*)'pgfr'
2481 !      write(24,9030) (pgfr(k),k=kts+1,kte)
2482 !      write(24,*)'pgacw'
2483 !      write(24,9030) (pgacw(k),k=kts+1,kte)
2484 !      write(24,*)'pgaci'
2485 !      write(24,9030) (pgaci(k),k=kts+1,kte)
2486 !      write(24,*)'pgacr'
2487 !      write(24,9030) (pgacr(k),k=kts+1,kte)
2488 !      write(24,*)'pgacs'
2489 !      write(24,9030) (pgacs(k),k=kts+1,kte)
2490 !      write(24,*)'pgacip'
2491 !      write(24,9030) (pgacip(k),k=kts+1,kte)
2492 !      write(24,*)'pgacrP'
2493 !      write(24,9030) (pgacrP(k),k=kts+1,kte)
2494 !      write(24,*)'pgacsp'
2495 !      write(24,9030) (pgacsp(k),k=kts+1,kte)
2496 !      write(24,*)'pgwet'
2497 !      write(24,9030) (pgwet(k),k=kts+1,kte)
2498 !      write(24,*)'pdry'
2499 !      write(24,9030) (pdry(k),k=kts+1,kte)
2500 !      write(24,*)'pgsub'
2501 !      write(24,9030) (pgsub(k),k=kts+1,kte)
2502 !      write(24,*)'pgdep'
2503 !      write(24,9030) (pgdep(k),k=kts+1,kte)
2504 !      write(24,*)'pgmlt'
2505 !      write(24,9030) (pgmlt(k),k=kts+1,kte)
2506 !      write(24,*)'pgmltevp'
2507 !      write(24,9030) (pgmltevp(k),k=kts+1,kte)
2508
2509
2510
2511 !**** below if qv < qvmin then qv=qvmin, ql=0.0, and qi=0.0
2512 !
2513       do k=kts+1,kte
2514          if ( qvz(k) .lt. qvmin ) then
2515             qlz(k)=0.0
2516             qiz(k)=0.0
2517             qvz(k)=amax1( qvmin,qvz(k)+qlz(k)+qiz(k) )
2518          end if
2519       enddo
2520 !
2521   END SUBROUTINE clphy1d
2522
2523
2524 !---------------------------------------------------------------------
2525 !                         SATURATED ADJUSTMENT
2526 !---------------------------------------------------------------------
2527       SUBROUTINE satadj(qvz, qlz, qiz, prez, theiz, thz, tothz,      &
2528                         kts, kte, k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0)
2529 !---------------------------------------------------------------------
2530    IMPLICIT NONE
2531 !---------------------------------------------------------------------
2532 !  This program use Newton's method for finding saturated temperature
2533 !  and saturation mixing ratio.
2534 !
2535 ! In this saturation adjustment scheme we assume
2536 ! (1)  the saturation mixing ratio is the mass weighted average of
2537 !      saturation values over liquid water (qsw), and ice (qsi)
2538 !      following Lord et al., 1984 and Tao, 1989
2539 !
2540 ! (2) the percentage of cloud liquid and cloud ice will
2541 !      be fixed during the saturation calculation
2542 !---------------------------------------------------------------------
2543 !
2544
2545   INTEGER, INTENT(IN   )             :: kts, kte, k
2546
2547   REAL,      DIMENSION( kts:kte ),                                   &
2548                        INTENT(INOUT) :: qvz, qlz, qiz
2549 !
2550   REAL,      DIMENSION( kts:kte ),                                   &
2551                        INTENT(IN   ) :: prez, theiz, tothz
2552
2553   REAL,     INTENT(IN   )            :: xLvocp, xLfocp, episp0k
2554   REAL,     INTENT(IN   )            :: EP2,SVP1,SVP2,SVP3,SVPT0
2555
2556 ! LOCAL VARS
2557
2558   INTEGER                            :: n
2559
2560   REAL, DIMENSION( kts:kte )         :: thz, tem, temcc, qsiz,       &
2561                                         qswz, qvsbar
2562
2563   REAL ::   qsat, qlpqi, ratql, t0, t1, tmp1, ratqi, tsat, absft,    &
2564             denom1, denom2, dqvsbar, ftsat, dftsat, qpz,             &
2565             gindex, es
2566
2567   REAL ::   tem_noliqice, qsat_noliqice !bloss
2568 !
2569 !---------------------------------------------------------------------
2570
2571       thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2572
2573       tem(k)=tothz(k)*thz(k)
2574
2575 !bloss: BEGIN
2576       tem_noliqice = tem(k) - xlvocp*qlz(k) - (xLvocp+xLfocp)*qiz(k)
2577
2578       if (tem_noliqice .gt. 273.15) then
2579          es=1000.*svp1*exp( svp2*(tem_noliqice-svpt0)/(tem_noliqice-svp3) )
2580          qsat_noliqice=ep2*es/(prez(k)-es)
2581       else
2582         qsat_noliqice=episp0k/prez(k)*  &
2583              exp( 21.8745584*(tem_noliqice-273.15)/(tem_noliqice-7.66) )
2584       end if
2585 !bloss: END
2586       qpz=qvz(k)+qlz(k)+qiz(k)
2587       if (qpz .lt. qsat_noliqice) then
2588          qvz(k)=qpz
2589          qiz(k)=0.0
2590          qlz(k)=0.0
2591          go to 400
2592       end if
2593       qlpqi=qlz(k)+qiz(k)
2594       if( qlpqi .ge. 1.0e-5) then
2595         ratql=qlz(k)/qlpqi
2596         ratqi=qiz(k)/qlpqi
2597       else
2598         t0=273.15
2599 !       t1=233.15
2600         t1=248.15
2601         tmp1=( t0-tem(k) )/(t0-t1)
2602         tmp1=amin1(1.0,tmp1)
2603         tmp1=amax1(0.0,tmp1)
2604         ratqi=tmp1
2605         ratql=1.0-tmp1
2606       end if
2607 !
2608 !
2609 !--  saturation mixing ratios over water and ice
2610 !--  at the outset we will follow Bolton 1980 MWR for
2611 !--  the water and Murray JAS 1967 for the ice
2612 !
2613 !-- dqvsbar=d(qvsbar)/dT
2614 !-- ftsat=F(Tsat)
2615 !-- dftsat=d(F(T))/dT
2616 !
2617 !  First guess of tsat
2618
2619       tsat=tem(k)
2620       absft=1.0
2621 !
2622       do 200 n=1,20
2623          denom1=1.0/(tsat-svp3)
2624          denom2=1.0/(tsat-7.66)
2625 !        qswz(k)=episp0k/prez(k)*  &
2626 !                exp( svp2*denom1*(tsat-273.15) )
2627          es=1000.*svp1*exp( svp2*denom1*(tsat-svpt0) )
2628          qswz(k)=ep2*es/(prez(k)-es)
2629          if (tem(k) .lt. 273.15) then
2630 !           qsiz(k)=episp0k/prez(k)*  &
2631 !                   exp( 21.8745584*denom2*(tsat-273.15) )
2632             es=1000.*svp1*exp( 21.8745584*denom2*(tsat-273.15) )
2633             qsiz(k)=ep2*es/(prez(k)-es)
2634             if (tem(k) .lt. 233.15) qswz(k)=qsiz(k)
2635          else
2636             qsiz(k)=qswz(k)
2637          endif
2638          qvsbar(k)=ratql*qswz(k)+ratqi*qsiz(k)
2639 !
2640 !        if( absft .lt. 0.01 .and. n .gt. 3 ) go to 300
2641          if( absft .lt. 0.01 ) go to 300
2642 !
2643          dqvsbar=ratql*qswz(k)*svp2*243.5*denom1*denom1+  &
2644                  ratqi*qsiz(k)*21.8745584*265.5*denom2*denom2
2645          ftsat=tsat+(xlvocp+ratqi*xlfocp)*qvsbar(k)-  &
2646                tothz(k)*theiz(k)-xlfocp*ratqi*(qvz(k)+qlz(k)+qiz(k))
2647          dftsat=1.0+(xlvocp+ratqi*xlfocp)*dqvsbar
2648          tsat=tsat-ftsat/dftsat
2649          absft=abs(ftsat)
2650
2651 200   continue
2652 9020  format(1x,'point can not converge, absft,n=',e12.5,i5)
2653 !
2654 300   continue
2655       if( qpz .gt. qvsbar(k) ) then
2656         qvz(k)=qvsbar(k)
2657         qiz(k)=ratqi*( qpz-qvz(k) )
2658         qlz(k)=ratql*( qpz-qvz(k) )
2659       else
2660         qvz(k)=qpz
2661         qiz(k)=0.0
2662         qlz(k)=0.0
2663       end if
2664  400  continue
2665
2666    END SUBROUTINE satadj
2667
2668
2669 !----------------------------------------------------------------
2670    REAL FUNCTION parama1(temp)
2671 !----------------------------------------------------------------
2672    IMPLICIT NONE
2673 !----------------------------------------------------------------
2674 !  This program calculate the parameter for crystal growth rate
2675 !  in Bergeron process
2676 !----------------------------------------------------------------
2677
2678       REAL, INTENT (IN   )   :: temp
2679       REAL, DIMENSION(32)    :: a1
2680       INTEGER                :: i1, i1p1
2681       REAL                   :: ratio
2682
2683       data a1/0.100e-10,0.7939e-7,0.7841e-6,0.3369e-5,0.4336e-5, &
2684               0.5285e-5,0.3728e-5,0.1852e-5,0.2991e-6,0.4248e-6, &
2685               0.7434e-6,0.1812e-5,0.4394e-5,0.9145e-5,0.1725e-4, &
2686               0.3348e-4,0.1725e-4,0.9175e-5,0.4412e-5,0.2252e-5, &
2687               0.9115e-6,0.4876e-6,0.3473e-6,0.4758e-6,0.6306e-6, &
2688               0.8573e-6,0.7868e-6,0.7192e-6,0.6513e-6,0.5956e-6, &
2689               0.5333e-6,0.4834e-6/
2690
2691       i1=int(-temp)+1
2692       i1p1=i1+1
2693       ratio=-(temp)-float(i1-1)
2694       parama1=a1(i1)+ratio*( a1(i1p1)-a1(i1) )
2695
2696    END FUNCTION parama1
2697
2698 !----------------------------------------------------------------
2699    REAL FUNCTION parama2(temp)
2700 !----------------------------------------------------------------
2701    IMPLICIT NONE
2702 !----------------------------------------------------------------
2703 !  This program calculate the parameter for crystal growth rate
2704 !  in Bergeron process
2705 !----------------------------------------------------------------
2706
2707       REAL, INTENT (IN   )   :: temp
2708       REAL, DIMENSION(32)    :: a2
2709       INTEGER                :: i1, i1p1
2710       REAL                   :: ratio
2711
2712       data a2/0.0100,0.4006,0.4831,0.5320,0.5307,0.5319,0.5249, &
2713               0.4888,0.3849,0.4047,0.4318,0.4771,0.5183,0.5463, &
2714               0.5651,0.5813,0.5655,0.5478,0.5203,0.4906,0.4447, &
2715               0.4126,0.3960,0.4149,0.4320,0.4506,0.4483,0.4460, &
2716               0.4433,0.4413,0.4382,0.4361/
2717       i1=int(-temp)+1
2718       i1p1=i1+1
2719       ratio=-(temp)-float(i1-1)
2720       parama2=a2(i1)+ratio*( a2(i1p1)-a2(i1) )
2721
2722    END FUNCTION parama2
2723
2724 !----------------------------------------------------------------
2725    REAL FUNCTION ggamma(X)
2726 !----------------------------------------------------------------
2727    IMPLICIT NONE
2728 !----------------------------------------------------------------
2729       REAL, INTENT(IN   ) :: x
2730       REAL, DIMENSION(8)  :: B
2731       INTEGER             ::j, K1
2732       REAL                ::PF, G1TO2 ,TEMP
2733
2734       DATA B/-.577191652,.988205891,-.897056937,.918206857,  &
2735              -.756704078,.482199394,-.193527818,.035868343/
2736
2737       PF=1.
2738       TEMP=X
2739       DO 10 J=1,200
2740       IF (TEMP .LE. 2) GO TO 20
2741       TEMP=TEMP-1.
2742    10 PF=PF*TEMP
2743   100 FORMAT(//,5X,'module_mp_lin: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
2744       WRITE(wrf_err_message,100)X
2745       CALL wrf_error_fatal(wrf_err_message)
2746    20 G1TO2=1.
2747       TEMP=TEMP - 1.
2748       DO 30 K1=1,8
2749    30 G1TO2=G1TO2 + B(K1)*TEMP**K1
2750       ggamma=PF*G1TO2
2751
2752       END FUNCTION ggamma
2753
2754 !----------------------------------------------------------------
2755
2756 !+---+-----------------------------------------------------------------+
2757       subroutine refl10cm_lin (qv1d, qr1d, qs1d, qg1d,                  &
2758                        t1d, p1d, dBZ, kts, kte, ii, jj)
2759
2760       IMPLICIT NONE
2761
2762 !..Sub arguments
2763       INTEGER, INTENT(IN):: kts, kte, ii, jj
2764       REAL, DIMENSION(kts:kte), INTENT(IN)::                            &
2765                       qv1d, qr1d, qs1d, qg1d, t1d, p1d
2766       REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
2767
2768 !..Local variables
2769       REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
2770       REAL, DIMENSION(kts:kte):: rr, rs, rg
2771
2772       DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg
2773       DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g
2774       DOUBLE PRECISION:: lamr, lams, lamg
2775       LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
2776
2777       REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
2778       DOUBLE PRECISION:: fmelt_s, fmelt_g
2779
2780       INTEGER:: i, k, k_0, kbot, n
2781       LOGICAL:: melti
2782
2783       DOUBLE PRECISION:: cback, x, eta, f_d
2784       REAL, PARAMETER:: R=287.
2785       REAL, PARAMETER:: PIx=3.1415926536
2786
2787 !+---+
2788
2789       do k = kts, kte
2790          dBZ(k) = -35.0
2791       enddo
2792
2793 !+---+-----------------------------------------------------------------+
2794 !..Put column of data into local arrays.
2795 !+---+-----------------------------------------------------------------+
2796       do k = kts, kte
2797          temp(k) = t1d(k)
2798          qv(k) = MAX(1.E-10, qv1d(k))
2799          pres(k) = p1d(k)
2800          rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
2801
2802          if (qr1d(k) .gt. 1.E-9) then
2803             rr(k) = qr1d(k)*rho(k)
2804             N0_r(k) = xnor
2805             lamr = (xam_r*xcrg(3)*N0_r(k)/rr(k))**(1./xcre(1))
2806             ilamr(k) = 1./lamr
2807             L_qr(k) = .true.
2808          else
2809             rr(k) = 1.E-12
2810             L_qr(k) = .false.
2811          endif
2812
2813          if (qs1d(k) .gt. 1.E-9) then
2814             rs(k) = qs1d(k)*rho(k)
2815             N0_s(k) = xnos
2816             lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1))
2817             ilams(k) = 1./lams
2818             L_qs(k) = .true.
2819          else
2820             rs(k) = 1.E-12
2821             L_qs(k) = .false.
2822          endif
2823
2824          if (qg1d(k) .gt. 1.E-9) then
2825             rg(k) = qg1d(k)*rho(k)
2826             N0_g(k) = xnog
2827             lamg = (xam_g*xcgg(3)*N0_g(k)/rg(k))**(1./xcge(1))
2828             ilamg(k) = 1./lamg
2829             L_qg(k) = .true.
2830          else
2831             rg(k) = 1.E-12
2832             L_qg(k) = .false.
2833          endif
2834       enddo
2835
2836 !+---+-----------------------------------------------------------------+
2837 !..Locate K-level of start of melting (k_0 is level above).
2838 !+---+-----------------------------------------------------------------+
2839       melti = .false.
2840       k_0 = kts
2841       do k = kte-1, kts, -1
2842          if ( (temp(k).gt.273.15) .and. L_qr(k)                         &
2843                                   .and. (L_qs(k+1).or.L_qg(k+1)) ) then
2844             k_0 = MAX(k+1, k_0)
2845             melti=.true.
2846             goto 195
2847          endif
2848       enddo
2849  195  continue
2850
2851 !+---+-----------------------------------------------------------------+
2852 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
2853 !.. and non-water-coated snow and graupel when below freezing are
2854 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
2855 !+---+-----------------------------------------------------------------+
2856
2857       do k = kts, kte
2858          ze_rain(k) = 1.e-22
2859          ze_snow(k) = 1.e-22
2860          ze_graupel(k) = 1.e-22
2861          if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
2862          if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PIx)*(6.0/PIx)     &
2863                                  * (xam_s/900.0)*(xam_s/900.0)          &
2864                                  * N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
2865          if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PIx)*(6.0/PIx)  &
2866                                     * (xam_g/900.0)*(xam_g/900.0)       &
2867                                     * N0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
2868       enddo
2869
2870
2871 !+---+-----------------------------------------------------------------+
2872 !..Special case of melting ice (snow/graupel) particles.  Assume the
2873 !.. ice is surrounded by the liquid water.  Fraction of meltwater is
2874 !.. extremely simple based on amount found above the melting level.
2875 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
2876 !.. routines).
2877 !+---+-----------------------------------------------------------------+
2878
2879       if (melti .and. k_0.ge.kts+1) then
2880        do k = k_0-1, kts, -1
2881
2882 !..Reflectivity contributed by melting snow
2883           if (L_qs(k) .and. L_qs(k_0) ) then
2884            fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
2885            eta = 0.d0
2886            lams = 1./ilams(k)
2887            do n = 1, nrbins
2888               x = xam_s * xxDs(n)**xbm_s
2889               call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
2890                     fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
2891                     CBACK, mixingrulestring_s, matrixstring_s,          &
2892                     inclusionstring_s, hoststring_s,                    &
2893                     hostmatrixstring_s, hostinclusionstring_s)
2894               f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
2895               eta = eta + f_d * CBACK * simpson(n) * xdts(n)
2896            enddo
2897            ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
2898           endif
2899
2900
2901 !..Reflectivity contributed by melting graupel
2902
2903           if (L_qg(k) .and. L_qg(k_0) ) then
2904            fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
2905            eta = 0.d0
2906            lamg = 1./ilamg(k)
2907            do n = 1, nrbins
2908               x = xam_g * xxDg(n)**xbm_g
2909               call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), &
2910                     fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
2911                     CBACK, mixingrulestring_g, matrixstring_g,          &
2912                     inclusionstring_g, hoststring_g,                    &
2913                     hostmatrixstring_g, hostinclusionstring_g)
2914               f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n))
2915               eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
2916            enddo
2917            ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
2918           endif
2919
2920        enddo
2921       endif
2922
2923       do k = kte, kts, -1
2924          dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
2925       enddo
2926
2927       end subroutine refl10cm_lin
2928 !+---+-----------------------------------------------------------------+
2929
2930 END MODULE module_mp_lin