PageRenderTime 62ms CodeModel.GetById 12ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_mp_lin.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2701 lines | 1365 code | 245 blank | 1091 comment | 91 complexity | 8ed673c9633037adad62d5d18307c8af MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. !WRF:MODEL_LAYER:PHYSICS
    2. 

  2. !
    3. 

  3. 

  4. MODULE module_mp_lin
    5. 

  5. 

  6.    USE     module_wrf_error
    7. 

  7. !
    8. 

  8.    REAL    , PARAMETER, PRIVATE ::       RH = 1.0
    9. 

  9. !  REAL    , PARAMETER, PRIVATE ::   episp0 = 0.622*611.21
    10. 

  10.    REAL    , PARAMETER, PRIVATE ::     xnor = 8.0e6
    11. 

  11.    REAL    , PARAMETER, PRIVATE ::     xnos = 3.0e6
    12. 

  12. 

  13. !  Lin
    14. 

  14. !  REAL    , PARAMETER, PRIVATE ::     xnog = 4.0e4
    15. 

  15. !  REAL    , PARAMETER, PRIVATE ::     rhograul = 917.
    16. 

  16. 

  17. ! Hobbs
    18. 

  18.   REAL     , PARAMETER, PRIVATE ::     xnog = 4.0e6
    19. 

  19.   REAL     , PARAMETER, PRIVATE ::     rhograul = 400.
    20. 

  20. 

  21. !
    22. 

  22.   REAL     , PARAMETER, PRIVATE ::                              &
    23. 

  23.              qi0 = 1.0e-3, ql0 = 7.0e-4, qs0 = 6.0E-4,          &
    24. 

  24.              xmi50 = 4.8e-10, xmi40 = 2.46e-10,                 &
    25. 

  25.              constb = 0.8, constd = 0.25,                       &
    26. 

  26.              o6 = 1./6.,  cdrag = 0.6,                          &
    27. 

  27.              avisc = 1.49628e-6, adiffwv = 8.7602e-5,           &
    28. 

  28.              axka = 1.4132e3, di50 = 1.0e-4, xmi = 4.19e-13,    &
    29. 

  29.              cw = 4.187e3, vf1s = 0.78, vf2s = 0.31,            &
    30. 

  30.              xni0 = 1.0e-2, xmnin = 1.05e-18, bni = 0.5,        &
    31. 

  31.              ci = 2.093e3
    32. 

  32. CONTAINS
    33. 

  33. 

  34. !-------------------------------------------------------------------
    35. 

  35. !  Lin et al., 1983, JAM, 1065-1092, and
    36. 

  36. !  Rutledge and Hobbs, 1984, JAS, 2949-2972
    37. 

  37. !  May 2009 - Changes and corrections from P. Blossey (U. Washington)
    38. 

  38. !-------------------------------------------------------------------
    39. 

  39.   SUBROUTINE lin_et_al(th                                          &
    40. 

  40.                       ,qv, ql, qr                                  &
    41. 

  41.                       ,qi, qs                                      &
    42. 

  42.                       ,rho, pii, p                                 &
    43. 

  43.                       ,dt_in                                       &
    44. 

  44.                       ,z,ht, dz8w                                  &
    45. 

  45.                       ,grav, cp, Rair, rvapor                      &
    46. 

  46.                       ,XLS, XLV, XLF, rhowater, rhosnow            &
    47. 

  47.                       ,EP2,SVP1,SVP2,SVP3,SVPT0                    &
    48. 

  48.                       , RAINNC, RAINNCV                            &
    49. 

  49.                       , SNOWNC, SNOWNCV                            &
    50. 

  50.                       , GRAUPELNC, GRAUPELNCV, SR                  &
    51. 

  51.                       ,ids,ide, jds,jde, kds,kde                   &
    52. 

  52.                       ,ims,ime, jms,jme, kms,kme                   &
    53. 

  53.                       ,its,ite, jts,jte, kts,kte                   &
    54. 

  54.                   ! Optional 
    55. 

  55.                       ,qlsink, precr, preci, precs, precg          &
    56. 

  56.                       , F_QG,F_QNDROP                              &
    57. 

  57.                       , qg, qndrop                                 &
    58. 

  58.                                                                    )
    59. 

  59. !-------------------------------------------------------------------
    60. 

  60.   IMPLICIT NONE
    61. 

  61. !-------------------------------------------------------------------
    62. 

  62. !
    63. 

  63. ! Shuhua 12/17/99
    64. 

  64. !
    65. 

  65.   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
    66. 

  66.                                       ims,ime, jms,jme, kms,kme , &
    67. 

  67.                                       its,ite, jts,jte, kts,kte
    68. 

  68. 

  69.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    70. 

  70.         INTENT(INOUT) ::                                          &
    71. 

  71.                                                               th, &
    72. 

  72.                                                               qv, &
    73. 

  73.                                                               ql, &
    74. 

  74.                                                               qr
    75. 

  75. 

  76. !
    77. 

  77.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    78. 

  78.         INTENT(IN   ) ::                                          &
    79. 

  79.                                                              rho, &
    80. 

  80.                                                              pii, &
    81. 

  81.                                                                p, &
    82. 

  82.                                                             dz8w
    83. 

  83. 

  84.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    85. 

  85.         INTENT(IN   ) ::                                       z
    86. 

  86. 

  87. 

  88. 

  89.   REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) ::       ht
    90. 

  90. 

  91.   REAL, INTENT(IN   ) ::                                   dt_in, &
    92. 

  92.                                                             grav, &
    93. 

  93.                                                             Rair, &
    94. 

  94.                                                           rvapor, &
    95. 

  95.                                                               cp, &
    96. 

  96.                                                              XLS, &
    97. 

  97.                                                              XLV, &
    98. 

  98.                                                              XLF, &
    99. 

  99.                                                         rhowater, &
    100. 

  100.                                                          rhosnow
    101. 

  101. 

  102.   REAL, INTENT(IN   ) :: EP2,SVP1,SVP2,SVP3,SVPT0
    103. 

  103. 

  104.   REAL, DIMENSION( ims:ime , jms:jme ),                           &
    105. 

  105.         INTENT(INOUT) ::                                  RAINNC, &
    106. 

  106.                                                          RAINNCV, &
    107. 

  107.                                                               SR
    108. 

  108. 

  109. ! Optional
    110. 

  110.   REAL, DIMENSION( ims:ime , jms:jme ),                           &
    111. 

  111.         OPTIONAL,                                                 &
    112. 

  112.         INTENT(INOUT) ::                                  SNOWNC, &
    113. 

  113.                                                          SNOWNCV, &
    114. 

  114.                                                        GRAUPELNC, &
    115. 

  115.                                                       GRAUPELNCV
    116. 

  116. 

  117.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    118. 

  118.         OPTIONAL,                                                 &
    119. 

  119.         INTENT(INOUT) ::                                          &
    120. 

  120.                                                               qi, &
    121. 

  121.                                                               qs, &
    122. 

  122.                                                               qg, &
    123. 

  123.                                                           qndrop
    124. 

  124. 

  125.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    126. 

  126.         OPTIONAL, INTENT(OUT   ) ::                               &
    127. 

  127.         qlsink, & ! cloud water conversion to rain (/s)
    128. 

  128.         precr,  & ! rain precipitation rate at all levels (kg/m2/s)
    129. 

  129.         preci,  & ! ice precipitation rate at all levels (kg/m2/s)
    130. 

  130.         precs,  & ! snow precipitation rate at all levels (kg/m2/s)
    131. 

  131.         precg     ! graupel precipitation rate at all levels (kg/m2/s)
    132. 

  132. 

  133.   LOGICAL, INTENT(IN), OPTIONAL ::                F_QG, F_QNDROP
    134. 

  134. 

  135. ! LOCAL VAR
    136. 

  136. 

  137.   INTEGER             ::                            min_q, max_q
    138. 

  138. 

  139.   REAL, DIMENSION( its:ite , jts:jte )                            &
    140. 

  140.                                ::        rain, snow, graupel,ice
    141. 

  141. 

  142.   REAL, DIMENSION( kts:kte )   ::                  qvz, qlz, qrz, &
    143. 

  143.                                                    qiz, qsz, qgz, &
    144. 

  144.                                                              thz, &
    145. 

  145.                                                      tothz, rhoz, &
    146. 

  146.                                                    orhoz, sqrhoz, &
    147. 

  147.                                                         prez, zz, &
    148. 

  148.                                   precrz, preciz, precsz, precgz, &
    149. 

  149.                                                          qndropz, &
    150. 

  150.                                                      dzw, preclw
    151. 

  151. 

  152.   LOGICAL :: flag_qg, flag_qndrop
    153. 

  153. !
    154. 

  154.   REAL    ::         dt, pptrain, pptsnow, pptgraul, rhoe_s,      &
    155. 

  155.                      gindex, pptice
    156. 

  156.   real :: qndropconst
    157. 

  157. 

  158.   INTEGER ::               i,j,k
    159. 

  159. !
    160. 

  160.    flag_qg     = .false.
    161. 

  161.    flag_qndrop = .false.
    162. 

  162.    IF ( PRESENT ( f_qg     ) ) flag_qg     = f_qg
    163. 

  163.    IF ( PRESENT ( f_qndrop ) ) flag_qndrop = f_qndrop
    164. 

  164. !
    165. 

  165.    dt=dt_in
    166. 

  166.    rhoe_s=1.29
    167. 

  167.    qndropconst=100.e6  !sg
    168. 

  168.    gindex=1.0
    169. 

  169. 

  170.    IF (.not.flag_qg) gindex=0.
    171. 

  171. 

  172.    j_loop:  DO j = jts, jte
    173. 

  173.    i_loop:  DO i = its, ite
    174. 

  174. !
    175. 

  175. !- write data from 3-D to 1-D
    176. 

  176. !
    177. 

  177.    DO k = kts, kte
    178. 

  178.       qvz(k)=qv(i,k,j)
    179. 

  179.       qlz(k)=ql(i,k,j)
    180. 

  180.       qrz(k)=qr(i,k,j)
    181. 

  181.       thz(k)=th(i,k,j)
    182. 

  182. !
    183. 

  183.       rhoz(k)=rho(i,k,j)
    184. 

  184.       orhoz(k)=1./rhoz(k)
    185. 

  185.       prez(k)=p(i,k,j)
    186. 

  186.       sqrhoz(k)=sqrt(rhoe_s*orhoz(k))
    187. 

  187.       tothz(k)=pii(i,k,j)
    188. 

  188.       zz(k)=z(i,k,j)
    189. 

  189.       dzw(k)=dz8w(i,k,j)
    190. 

  190.    END DO
    191. 

  191. 

  192.    IF (flag_qndrop .AND. PRESENT( qndrop )) THEN
    193. 

  193.      DO k = kts, kte
    194. 

  194.          qndropz(k)=qndrop(i,k,j)
    195. 

  195.       ENDDO
    196. 

  196.    ELSE
    197. 

  197.       DO k = kts, kte
    198. 

  198.          qndropz(k)=qndropconst
    199. 

  199.       ENDDO
    200. 

  200.    ENDIF
    201. 

  201. 

  202.    DO k = kts, kte
    203. 

  203.       qiz(k)=qi(i,k,j)
    204. 

  204.       qsz(k)=qs(i,k,j)
    205. 

  205.    ENDDO
    206. 

  206. 

  207.    IF ( flag_qg .AND. PRESENT( qg ) ) THEN
    208. 

  208.       DO k = kts, kte
    209. 

  209.          qgz(k)=qg(i,k,j)
    210. 

  210.       ENDDO
    211. 

  211.    ELSE
    212. 

  212.       DO k = kts, kte
    213. 

  213.          qgz(k)=0.
    214. 

  214.       ENDDO
    215. 

  215.    ENDIF
    216. 

  216. !
    217. 

  217.    pptrain=0.
    218. 

  218.    pptsnow=0.
    219. 

  219.    pptgraul=0.
    220. 

  220.    pptice=0.
    221. 

  221.    CALL clphy1d(    dt, qvz, qlz, qrz, qiz, qsz, qgz,         &
    222. 

  222.                     qndropz,flag_qndrop,                      &
    223. 

  223.                     thz, tothz, rhoz, orhoz, sqrhoz,          &
    224. 

  224.                     prez, zz, dzw, ht(I,J), preclw,           &
    225. 

  225.                     precrz, preciz, precsz, precgz,           &
    226. 

  226.                     pptrain, pptsnow, pptgraul, pptice,       &
    227. 

  227.                     grav, cp, Rair, rvapor, gindex,           &
    228. 

  228.                     XLS, XLV, XLF, rhowater, rhosnow,         &
    229. 

  229.                     EP2,SVP1,SVP2,SVP3,SVPT0,                 &
    230. 

  230.                     kts, kte, i, j                            )
    231. 

  231. 

  232. !
    233. 

  233. ! Precipitation from cloud microphysics -- only for one time step
    234. 

  234. !
    235. 

  235. ! unit is transferred from m to mm
    236. 

  236. 

  237. !
    238. 

  238.    rain(i,j)=pptrain
    239. 

  239.    snow(i,j)=pptsnow
    240. 

  240.    graupel(i,j)=pptgraul
    241. 

  241.    ice(i,j)=pptice
    242. 

  242.    sr(i,j)=(pptice+pptsnow+pptgraul)/(pptice+pptsnow+pptgraul+pptrain+1.e-12)
    243. 

  243. !
    244. 

  244.    RAINNCV(i,j)= pptrain + pptsnow + pptgraul + pptice
    245. 

  245.    RAINNC(i,j)=RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
    246. 

  246.    IF(PRESENT(SNOWNCV))SNOWNCV(i,j)= pptsnow + pptice
    247. 

  247.    IF(PRESENT(SNOWNC))SNOWNC(i,j)=SNOWNC(i,j) + pptsnow + pptice
    248. 

  248.    IF(PRESENT(GRAUPELNCV))GRAUPELNCV(i,j)= pptgraul
    249. 

  249.    IF(PRESENT(GRAUPELNC))GRAUPELNC(i,j)=GRAUPELNC(i,j) + pptgraul
    250. 

  250. 

  251. !
    252. 

  252. !- update data from 1-D back to 3-D
    253. 

  253. !
    254. 

  254. !
    255. 

  255.    IF ( present(qlsink) .and. present(precr) ) THEN !sg beg
    256. 

  256.       DO k = kts, kte
    257. 

  257.          if(ql(i,k,j)>1.e-20) then
    258. 

  258.             qlsink(i,k,j)=-preclw(k)/ql(i,k,j)
    259. 

  259.          else
    260. 

  260.             qlsink(i,k,j)=0.
    261. 

  261.          endif
    262. 

  262.          precr(i,k,j)=precrz(k)
    263. 

  263.       END DO
    264. 

  264.    END IF                                          !sg end
    265. 

  265. 

  266.    DO k = kts, kte
    267. 

  267.       qv(i,k,j)=qvz(k)
    268. 

  268.       ql(i,k,j)=qlz(k)
    269. 

  269.       qr(i,k,j)=qrz(k)
    270. 

  270.       th(i,k,j)=thz(k)
    271. 

  271.    END DO
    272. 

  272. !
    273. 

  273.    IF ( flag_qndrop .AND. PRESENT( qndrop ) ) THEN !sg beg
    274. 

  274.       DO k = kts, kte
    275. 

  275.          qndrop(i,k,j)=qndropz(k)
    276. 

  276.       ENDDO
    277. 

  277.    END IF                                          !sg end
    278. 

  278. 

  279.    DO k = kts, kte
    280. 

  280.       qi(i,k,j)=qiz(k)
    281. 

  281.       qs(i,k,j)=qsz(k)
    282. 

  282.    ENDDO
    283. 

  283. 

  284.    IF ( present(preci) ) THEN     !sg beg
    285. 

  285.       DO k = kts, kte
    286. 

  286.          preci(i,k,j)=preciz(k)
    287. 

  287.       ENDDO
    288. 

  288.    END IF
    289. 

  289.       
    290. 

  290.    IF ( present(precs) ) THEN
    291. 

  291.       DO k = kts, kte
    292. 

  292.          precs(i,k,j)=precsz(k)
    293. 

  293.       ENDDO
    294. 

  294.    END IF                         !sg end
    295. 

  295.       
    296. 

  296.    IF ( flag_qg .AND. PRESENT( qg ) ) THEN
    297. 

  297.       DO k = kts, kte
    298. 

  298.          qg(i,k,j)=qgz(k)
    299. 

  299.       ENDDO
    300. 

  300. 

  301.       IF ( present(precg) ) THEN  !sg beg
    302. 

  302.          DO k = kts, kte
    303. 

  303.             precg(i,k,j)=precgz(k)
    304. 

  304.          ENDDO                    !sg end
    305. 

  305.       END IF
    306. 

  306.    ELSE                           !sg beg
    307. 

  307.       IF ( present(precg) ) precg(i,:,j)=0.  !sg end
    308. 

  308.    ENDIF
    309. 

  309. !
    310. 

  310.    ENDDO i_loop
    311. 

  311.    ENDDO j_loop
    312. 

  312. 

  313.    END SUBROUTINE lin_et_al
    314. 

  314. 

  315. 

  316. !-----------------------------------------------------------------------
    317. 

  317.    SUBROUTINE clphy1d(dt, qvz, qlz, qrz, qiz, qsz, qgz,                &
    318. 

  318.                       qndropz,flag_qndrop,                             &
    319. 

  319.                       thz, tothz, rho, orho, sqrho,                    &
    320. 

  320.                       prez, zz, dzw, zsfc, preclw,                     &
    321. 

  321.                       precrz, preciz, precsz, precgz,                  &
    322. 

  322.                       pptrain, pptsnow, pptgraul,                      &
    323. 

  323.                       pptice, grav, cp, Rair, rvapor, gindex,          &
    324. 

  324.                       XLS, XLV, XLF, rhowater, rhosnow,                &
    325. 

  325.                       EP2,SVP1,SVP2,SVP3,SVPT0,                        &
    326. 

  326.                       kts, kte, i, j                                   )
    327. 

  327. !-----------------------------------------------------------------------
    328. 

  328.     IMPLICIT NONE
    329. 

  329. !-----------------------------------------------------------------------
    330. 

  330. !  This program handles the vertical 1-D cloud micphysics
    331. 

  331. !-----------------------------------------------------------------------
    332. 

  332. ! avisc: constant in empirical formular for dynamic viscosity of air
    333. 

  333. !         =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
    334. 

  334. ! adiffwv: constant in empirical formular for diffusivity of water
    335. 

  335. !          vapor in air
    336. 

  336. !          = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
    337. 

  337. ! axka: constant in empirical formular for thermal conductivity of air
    338. 

  338. !       = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
    339. 

  339. ! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
    340. 

  340. ! xmi50: mass of a 50 micron ice crystal
    341. 

  341. !        = 4.8e-10 [kg] =4.8e-7 [g]
    342. 

  342. ! xmi40: mass of a 40 micron ice crystal
    343. 

  343. !        = 2.46e-10 [kg] = 2.46e-7 [g]
    344. 

  344. ! di50: diameter of a 50 micro (radius) ice crystal
    345. 

  345. !       =1.0e-4 [m]
    346. 

  346. ! xmi: mass of one cloud ice crystal
    347. 

  347. !      =4.19e-13 [kg] = 4.19e-10 [g]
    348. 

  348. ! oxmi=1.0/xmi
    349. 

  349. !
    350. 

  350. ! xni0=1.0e-2 [m-3] The value given in Lin et al. is wrong.(see
    351. 

  351. !                   Hsie et al.(1980) and Rutledge and Hobbs(1983) )
    352. 

  352. ! bni=0.5 [K-1]
    353. 

  353. ! xmnin: mass of a natural ice nucleus
    354. 

  354. !    = 1.05e-18 [kg] = 1.05e-15 [g] This values is suggested by
    355. 

  355. !                    Hsie et al. (1980)
    356. 

  356. !    = 1.0e-12 [kg] suggested by Rutlegde and Hobbs (1983)
    357. 

  357. ! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3
    358. 

  358. ! consta: constant in empirical formular for terminal
    359. 

  359. !         velocity of raindrop
    360. 

  360. !         =2115.0 [cm**(1-b)/s] = 2115.0*0.01**(1-b) [m**(1-b)/s]
    361. 

  361. ! constb: constant in empirical formular for terminal
    362. 

  362. !         velocity of raindrop
    363. 

  363. !         =0.8
    364. 

  364. ! xnor: intercept parameter of the raindrop size distribution
    365. 

  365. !       = 0.08 cm-4 = 8.0e6 m-4
    366. 

  366. ! araut: time sacle for autoconversion of cloud water to raindrops
    367. 

  367. !       =1.0e-3 [s-1]
    368. 

  368. ! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg]
    369. 

  369. ! vf1r: ventilation factors for rain =0.78
    370. 

  370. ! vf2r: ventilation factors for rain =0.31
    371. 

  371. ! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3
    372. 

  372. ! constc: constant in empirical formular for terminal
    373. 

  373. !         velocity of snow
    374. 

  374. !         =152.93 [cm**(1-d)/s] = 152.93*0.01**(1-d) [m**(1-d)/s]
    375. 

  375. ! constd: constant in empirical formular for terminal
    376. 

  376. !         velocity of snow
    377. 

  377. !         =0.25
    378. 

  378. ! xnos: intercept parameter of the snowflake size distribution
    379. 

  379. ! vf1s: ventilation factors for snow =0.78
    380. 

  380. ! vf2s: ventilation factors for snow =0.31
    381. 

  381. !
    382. 

  382. !----------------------------------------------------------------------
    383. 

  383. 

  384.   INTEGER, INTENT(IN   )               :: kts, kte, i, j
    385. 

  385. 

  386.   REAL,    DIMENSION( kts:kte ),                                      &
    387. 

  387.            INTENT(INOUT)               :: qvz, qlz, qrz, qiz, qsz,    &
    388. 

  388.                                           qndropz,                    &
    389. 

  389.                                           qgz, thz
    390. 

  390. 

  391.   REAL,    DIMENSION( kts:kte ),                                      &
    392. 

  392.            INTENT(IN   )               :: tothz, rho, orho, sqrho,    &
    393. 

  393.                                           prez, zz, dzw
    394. 

  394. 

  395.   REAL,    INTENT(IN   )               :: dt, grav, cp, Rair, rvapor, &
    396. 

  396.                                           XLS, XLV, XLF, rhowater,    &
    397. 

  397.                                           rhosnow,EP2,SVP1,SVP2,SVP3,SVPT0
    398. 

  398. 

  399.   REAL,    DIMENSION( kts:kte ), INTENT(OUT)               :: preclw, &
    400. 

  400.   		    precrz, preciz, precsz, precgz
    401. 

  401. 

  402.   REAL,    INTENT(INOUT)               :: pptrain, pptsnow, pptgraul, pptice
    403. 

  403. 

  404.   REAL,    INTENT(IN   )               :: zsfc
    405. 

  405.   logical, intent(in)                  :: flag_qndrop !sg
    406. 

  406. 

  407. ! local vars
    408. 

  408. 

  409.    REAL                                :: obp4, bp3, bp5, bp6, odp4,  &
    410. 

  410.                                           dp3, dp5, dp5o2
    411. 

  411. 

  412. 

  413. ! temperary vars
    414. 

  414. 

  415.    REAL                                :: tmp, tmp0, tmp1, tmp2,tmp3,  &
    416. 

  416.                                           tmp4,delta2,delta3, delta4,  &
    417. 

  417.                                           tmpa,tmpb,tmpc,tmpd,alpha1,  &
    418. 

  418.                                           qic, abi,abr, abg, odtberg,  &
    419. 

  419.                                           vti50,eiw,eri,esi,esr, esw,  &
    420. 

  420.                                           erw,delrs,term0,term1,araut, &
    421. 

  421.                                           constg2, vf1r, vf2r,alpha2,  &
    422. 

  422.                                           Ap, Bp, egw, egi, egs, egr,  &
    423. 

  423.                                           constg, gdelta4, g1sdelt4,   &
    424. 

  424.                                           factor, tmp_r, tmp_s,tmp_g,  &
    425. 

  425.                                           qlpqi, rsat, a1, a2, xnin
    426. 

  426. 

  427.   INTEGER                              :: k
    428. 

  428. !
    429. 

  429.   REAL, DIMENSION( kts:kte )    ::  oprez, tem, temcc, theiz, qswz,    &
    430. 

  430.                                     qsiz, qvoqswz, qvoqsiz, qvzodt,    &
    431. 

  431.                                     qlzodt, qizodt, qszodt, qrzodt,    &
    432. 

  432.                                     qgzodt
    433. 

  433. 

  434.   REAL, DIMENSION( kts:kte )    :: psnow, psaut, psfw,  psfi,  praci,  &
    435. 

  435.                                    piacr, psaci, psacw, psdep, pssub,  &
    436. 

  436.                                    pracs, psacr, psmlt, psmltevp,      &
    437. 

  437.                                    prain, praut, pracw, prevp, pvapor, &
    438. 

  438.                                    pclw,  pladj, pcli,  pimlt, pihom,  &
    439. 

  439.                                    pidw,  piadj, pgraupel, pgaut,      &
    440. 

  440.                                    pgfr,  pgacw, pgaci, pgacr, pgacs,  &
    441. 

  441.                                    pgacip,pgacrp,pgacsp,pgwet, pdry,   &
    442. 

  442.                                    pgsub, pgdep, pgmlt, pgmltevp,      &
    443. 

  443.                                    qschg, qgchg
    444. 

  444. !
    445. 

  445. 

  446.   REAL, DIMENSION( kts:kte )    :: qvsbar, rs0, viscmu, visc, diffwv,  &
    447. 

  447.                                    schmidt, xka
    448. 

  448. 

  449.   REAL, DIMENSION( kts:kte )    :: vtr, vts, vtg,                      &
    450. 

  450.                                    vtrold, vtsold, vtgold, vtiold,     &
    451. 

  451.                                    xlambdar, xlambdas, xlambdag,       &
    452. 

  452.                                    olambdar, olambdas, olambdag
    453. 

  453. 

  454.   REAL                          :: episp0k, dtb, odtb, pi, pio4,       &
    455. 

  455.                                    pio6, oxLf, xLvocp, xLfocp, consta, &
    456. 

  456.                                    constc, ocdrag, gambp4, gamdp4,     &
    457. 

  457.                                    gam4pt5, Cpor, oxmi, gambp3, gamdp3,&
    458. 

  458.                                    gambp6, gam3pt5, gam2pt75, gambp5o2,&
    459. 

  459.                                    gamdp5o2, cwoxlf, ocp, xni50, es
    460. 

  460. !
    461. 

  461.   REAL                          :: qvmin=1.e-20
    462. 

  462.   REAL                          :: gindex
    463. 

  463.   REAL                          :: temc1,save1,save2,xni50mx
    464. 

  464. 

  465. ! for terminal velocity flux
    466. 

  466. 

  467.   INTEGER                       :: min_q, max_q
    468. 

  468.   REAL                          :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
    469. 

  469.   LOGICAL                       :: notlast
    470. 

  470. !
    471. 

  471.   REAL                          :: tmp_tem, tmp_temcc !bloss
    472. 

  472. !
    473. 

  473.   real :: liqconc, dis, beta, kappa, p0, xc, capn,rhocgs
    474. 

  474. 

  475. !sg: begin
    476. 

  476. ! liqconc = liquid water content (g cm^-3)
    477. 

  477. ! capn = droplet number concentration (# cm^-3)
    478. 

  478. ! dis = relative dispersion (dimensionless) between 0.2 and 1.
    479. 

  479. ! Written by Yangang Liu based on Liu et al., GRL 32, 2005.
    480. 

  480. ! Autoconversion rate p = p0 * (threshold function)
    481. 

  481. ! p0 = "base" autoconversion rate (g cm^-3 s^-1)
    482. 

  482. ! kappa = constant in Long kernel = [kappa2 * (3/(4*pi*rhow))^3] in Liu papers
    483. 

  483. ! beta = Condensation rate constant = (beta6)^6 in Liu papers
    484. 

  484. ! xc = Normalized critical mass
    485. 

  485. ! ***********************************************************
    486. 

  486.   if(flag_qndrop)then
    487. 

  487.      dis = 0.5 ! droplet dispersion, set to 0.5 per SG 8-Nov-2006
    488. 

  488. !    Give  empirical constants
    489. 

  489.      kappa=1.1d10
    490. 

  490. !    Calculate Condensation rate constant
    491. 

  491.      beta = (1.0d0+3.0d0*dis**2)*(1.0d0+4.0d0*dis**2)*    &
    492. 

  492.          (1.0d0+5.0d0*dis**2)/((1.0d0+dis**2)*(1.0d0+2.0d0*dis**2))
    493. 

  493.   endif
    494. 

  494. !sg: end
    495. 

  495. 

  496.    dtb=dt
    497. 

  497.    odtb=1./dtb
    498. 

  498.    pi=acos(-1.)
    499. 

  499.    pio4=acos(-1.)/4.
    500. 

  500.    pio6=acos(-1.)/6.
    501. 

  501.    ocp=1./cp
    502. 

  502.    oxLf=1./xLf
    503. 

  503.    xLvocp=xLv/cp
    504. 

  504.    xLfocp=xLf/cp
    505. 

  505.    consta=2115.0*0.01**(1-constb)
    506. 

  506.    constc=152.93*0.01**(1-constd)
    507. 

  507.    ocdrag=1./Cdrag
    508. 

  508. !  episp0k=RH*episp0
    509. 

  509.    episp0k=RH*ep2*1000.*svp1
    510. 

  510. !
    511. 

  511.    gambp4=ggamma(constb+4.)
    512. 

  512.    gamdp4=ggamma(constd+4.)
    513. 

  513.    gam4pt5=ggamma(4.5)
    514. 

  514.    Cpor=cp/Rair
    515. 

  515.    oxmi=1.0/xmi
    516. 

  516.    gambp3=ggamma(constb+3.)
    517. 

  517.    gamdp3=ggamma(constd+3.)
    518. 

  518.    gambp6=ggamma(constb+6)
    519. 

  519.    gam3pt5=ggamma(3.5)
    520. 

  520.    gam2pt75=ggamma(2.75)
    521. 

  521.    gambp5o2=ggamma((constb+5.)/2.)
    522. 

  522.    gamdp5o2=ggamma((constd+5.)/2.)
    523. 

  523.    cwoxlf=cw/xlf
    524. 

  524.    delta2=0.
    525. 

  525.    delta3=0.
    526. 

  526.    delta4=0.
    527. 

  527. !
    528. 

  528. !-----------------------------------------------------------------------
    529. 

  529. !     oprez       1./prez ( prez : pressure)
    530. 

  530. !     qsw         saturated mixing ratio on water surface
    531. 

  531. !     qsi         saturated mixing ratio on ice surface
    532. 

  532. !     episp0k     RH*e*saturated pressure at 273.15 K
    533. 

  533. !     qvoqsw      qv/qsw
    534. 

  534. !     qvoqsi      qv/qsi
    535. 

  535. !     qvzodt      qv/dt
    536. 

  536. !     qlzodt      ql/dt
    537. 

  537. !     qizodt      qi/dt
    538. 

  538. !     qszodt      qs/dt
    539. 

  539. !     qrzodt      qr/dt
    540. 

  540. !     qgzodt      qg/dt
    541. 

  541. !
    542. 

  542. !     temcc       temperature in dregee C
    543. 

  543. !
    544. 

  544. 

  545.       obp4=1.0/(constb+4.0)
    546. 

  546.       bp3=constb+3.0
    547. 

  547.       bp5=constb+5.0
    548. 

  548.       bp6=constb+6.0
    549. 

  549.       odp4=1.0/(constd+4.0)
    550. 

  550.       dp3=constd+3.0
    551. 

  551.       dp5=constd+5.0
    552. 

  552.       dp5o2=0.5*(constd+5.0)
    553. 

  553. !
    554. 

  554.       do k=kts,kte
    555. 

  555.          oprez(k)=1./prez(k)
    556. 

  556.       enddo
    557. 

  557. 

  558.       do k=kts,kte
    559. 

  559.          qlz(k)=amax1( 0.0,qlz(k) )
    560. 

  560.          qiz(k)=amax1( 0.0,qiz(k) )
    561. 

  561.          qvz(k)=amax1( qvmin,qvz(k) )
    562. 

  562.          qsz(k)=amax1( 0.0,qsz(k) )
    563. 

  563.          qrz(k)=amax1( 0.0,qrz(k) )
    564. 

  564.          qgz(k)=amax1( 0.0,qgz(k) )
    565. 

  565.          qndropz(k)=amax1( 0.0,qndropz(k) )     !sg
    566. 

  566. !
    567. 

  567.          tem(k)=thz(k)*tothz(k)
    568. 

  568.          temcc(k)=tem(k)-273.15
    569. 

  569. !
    570. 

  570. !        qswz(k)=episp0k*oprez(k)* &
    571. 

  571. !               exp( svp2*temcc(k)/(tem(k)-svp3) )
    572. 

  572.          es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
    573. 

  573.          qswz(k)=ep2*es/(prez(k)-es)
    574. 

  574.          if (tem(k) .lt. 273.15 ) then
    575. 

  575. !           qsiz(k)=episp0k*oprez(k)* &
    576. 

  576. !                  exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
    577. 

  577.             es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
    578. 

  578.             qsiz(k)=ep2*es/(prez(k)-es)
    579. 

  579.             if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
    580. 

  580.          else
    581. 

  581.             qsiz(k)=qswz(k)
    582. 

  582.          endif
    583. 

  583. !
    584. 

  584.          qvoqswz(k)=qvz(k)/qswz(k)
    585. 

  585.          qvoqsiz(k)=qvz(k)/qsiz(k)
    586. 

  586. 

  587.          theiz(k)=thz(k)+(xlvocp*qvz(k)-xlfocp*qiz(k))/tothz(k)
    588. 

  588.       enddo
    589. 

  589. 

  590. 

  591. !
    592. 

  592. !
    593. 

  593. !-----------------------------------------------------------------------
    594. 

  594. ! In this simple stable cloud parameterization scheme, only five
    595. 

  595. ! forms of water substance (water vapor, cloud water, cloud ice,
    596. 

  596. ! rain and snow are considered. The prognostic variables are total
    597. 

  597. ! water (qp),cloud water (ql), and cloud ice (qi). Rain and snow are
    598. 

  598. ! diagnosed following Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7).
    599. 

  599. ! the micro physics are based on (1) Hsie et al.,1980, JAM, 950-977 ;
    600. 

  600. ! (2) Lin et al., 1983, JAM, 1065-1092 ; (3) Rutledge and Hobbs, 1983,
    601. 

  601. ! JAS, 1185-1206 ; (4) Rutledge and Hobbs, 1984, JAS, 2949-2972.
    602. 

  602. !-----------------------------------------------------------------------
    603. 

  603. !
    604. 

  604. ! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3
    605. 

  605. ! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3
    606. 

  606. ! xnor: intercept parameter of the raindrop size distribution
    607. 

  607. !       = 0.08 cm-4 = 8.0e6 m-4
    608. 

  608. ! xnos: intercept parameter of the snowflake size distribution
    609. 

  609. !       = 0.03 cm-4 = 3.0e6 m-4
    610. 

  610. ! xnog: intercept parameter of the graupel size distribution
    611. 

  611. !       = 4.0e-4 cm-4 = 4.0e4 m-4
    612. 

  612. ! consta: constant in empirical formular for terminal
    613. 

  613. !         velocity of raindrop
    614. 

  614. !         =2115.0 [cm**(1-b)/s] = 2115.0*0.01**(1-b) [m**(1-b)/s]
    615. 

  615. ! constb: constant in empirical formular for terminal
    616. 

  616. !         velocity of raindrop
    617. 

  617. !         =0.8
    618. 

  618. ! constc: constant in empirical formular for terminal
    619. 

  619. !         velocity of snow
    620. 

  620. !         =152.93 [cm**(1-d)/s] = 152.93*0.01**(1-d) [m**(1-d)/s]
    621. 

  621. ! constd: constant in empirical formular for terminal
    622. 

  622. !         velocity of snow
    623. 

  623. !         =0.25
    624. 

  624. ! avisc: constant in empirical formular for dynamic viscosity of air
    625. 

  625. !         =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
    626. 

  626. ! adiffwv: constant in empirical formular for diffusivity of water
    627. 

  627. !          vapor in air
    628. 

  628. !          = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
    629. 

  629. ! axka: constant in empirical formular for thermal conductivity of air
    630. 

  630. !       = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
    631. 

  631. ! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
    632. 

  632. !      = 1.0e-3 g/g = 1.0e-3 kg/gk
    633. 

  633. ! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg]
    634. 

  634. !      = 2.0e-3 g/g = 2.0e-3 kg/gk
    635. 

  635. ! qs0: mixing ratio threshold for snow aggregation
    636. 

  636. !      = 6.0e-4 g/g = 6.0e-4 kg/gk
    637. 

  637. ! xmi50: mass of a 50 micron ice crystal
    638. 

  638. !        = 4.8e-10 [kg] =4.8e-7 [g]
    639. 

  639. ! xmi40: mass of a 40 micron ice crystal
    640. 

  640. !        = 2.46e-10 [kg] = 2.46e-7 [g]
    641. 

  641. ! di50: diameter of a 50 micro (radius) ice crystal
    642. 

  642. !       =1.0e-4 [m]
    643. 

  643. ! xmi: mass of one cloud ice crystal
    644. 

  644. !      =4.19e-13 [kg] = 4.19e-10 [g]
    645. 

  645. ! oxmi=1.0/xmi
    646. 

  646. !
    647. 

  647. 

  648. 

  649. ! if gindex=1.0 include graupel
    650. 

  650. ! if gindex=0. no graupel
    651. 

  651. !
    652. 

  652. !
    653. 

  653.       do k=kts,kte
    654. 

  654.          psnow(k)=0.0
    655. 

  655.          psaut(k)=0.0
    656. 

  656.          psfw(k)=0.0
    657. 

  657.          psfi(k)=0.0
    658. 

  658.          praci(k)=0.0
    659. 

  659.          piacr(k)=0.0
    660. 

  660.          psaci(k)=0.0
    661. 

  661.          psacw(k)=0.0
    662. 

  662.          psdep(k)=0.0
    663. 

  663.          pssub(k)=0.0
    664. 

  664.          pracs(k)=0.0
    665. 

  665.          psacr(k)=0.0
    666. 

  666.          psmlt(k)=0.0
    667. 

  667.          psmltevp(k)=0.0
    668. 

  668. !
    669. 

  669.          prain(k)=0.0
    670. 

  670.          praut(k)=0.0
    671. 

  671.          pracw(k)=0.0
    672. 

  672.          prevp(k)=0.0
    673. 

  673. !
    674. 

  674.          pvapor(k)=0.0
    675. 

  675. !
    676. 

  676.          pclw(k)=0.0
    677. 

  677.          preclw(k)=0.0       !sg
    678. 

  678.          pladj(k)=0.0
    679. 

  679. !
    680. 

  680.          pcli(k)=0.0
    681. 

  681.          pimlt(k)=0.0
    682. 

  682.          pihom(k)=0.0
    683. 

  683.          pidw(k)=0.0
    684. 

  684.          piadj(k)=0.0
    685. 

  685.       enddo
    686. 

  686. 

  687. !
    688. 

  688. !!! graupel
    689. 

  689. !
    690. 

  690.       do k=kts,kte
    691. 

  691.          pgraupel(k)=0.0
    692. 

  692.          pgaut(k)=0.0
    693. 

  693.          pgfr(k)=0.0
    694. 

  694.          pgacw(k)=0.0
    695. 

  695.          pgaci(k)=0.0
    696. 

  696.          pgacr(k)=0.0
    697. 

  697.          pgacs(k)=0.0
    698. 

  698.          pgacip(k)=0.0
    699. 

  699.          pgacrP(k)=0.0
    700. 

  700.          pgacsp(k)=0.0
    701. 

  701.          pgwet(k)=0.0
    702. 

  702.          pdry(k)=0.0
    703. 

  703.          pgsub(k)=0.0
    704. 

  704.          pgdep(k)=0.0
    705. 

  705.          pgmlt(k)=0.0
    706. 

  706.          pgmltevp(k)=0.0
    707. 

  707.          qschg(k)=0.
    708. 

  708.          qgchg(k)=0.
    709. 

  709.       enddo
    710. 

  710. !
    711. 

  711. !
    712. 

  712. ! Set rs0=episp0*oprez(k)
    713. 

  713. ! episp0=e*saturated pressure at 273.15 K
    714. 

  714. ! e     = 0.622
    715. 

  715. !
    716. 

  716.       DO k=kts,kte
    717. 

  717.          rs0(k)=ep2*1000.*svp1/(prez(k)-1000.*svp1)
    718. 

  718.       END DO
    719. 

  719. !
    720. 

  720. !***********************************************************************
    721. 

  721. ! Calculate precipitation fluxes due to terminal velocities.
    722. 

  722. !***********************************************************************
    723. 

  723. !
    724. 

  724. !- Calculate termianl velocity (vt?)  of precipitation q?z
    725. 

  725. !- Find maximum vt? to determine the small delta t
    726. 

  726. !
    727. 

  727. !-- rain
    728. 

  728. !
    729. 

  729.     t_del_tv=0.
    730. 

  730.     del_tv=dtb
    731. 

  731.     notlast=.true.
    732. 

  732.     DO while (notlast)
    733. 

  733. !
    734. 

  734.       min_q=kte
    735. 

  735.       max_q=kts-1
    736. 

  736. !
    737. 

  737.       do k=kts,kte-1
    738. 

  738.          if (qrz(k) .gt. 1.0e-8) then
    739. 

  739.             min_q=min0(min_q,k)
    740. 

  740.             max_q=max0(max_q,k)
    741. 

  741.             tmp1=sqrt(pi*rhowater*xnor/rho(k)/qrz(k))
    742. 

  742.             tmp1=sqrt(tmp1)
    743. 

  743.             vtrold(k)=o6*consta*gambp4*sqrho(k)/tmp1**constb
    744. 

  744.             if (k .eq. 1) then
    745. 

  745.                del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtrold(k))
    746. 

  746.             else
    747. 

  747.                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtrold(k))
    748. 

  748.             endif
    749. 

  749.          else
    750. 

  750.             vtrold(k)=0.
    751. 

  751.          endif
    752. 

  752.       enddo
    753. 

  753. 

  754.       if (max_q .ge. min_q) then
    755. 

  755. !
    756. 

  756. !- Check if the summation of the small delta t >=  big delta t
    757. 

  757. !             (t_del_tv)          (del_tv)             (dtb)
    758. 

  758. 

  759.          t_del_tv=t_del_tv+del_tv
    760. 

  760. !
    761. 

  761.          if ( t_del_tv .ge. dtb ) then
    762. 

  762.               notlast=.false.
    763. 

  763.               del_tv=dtb+del_tv-t_del_tv
    764. 

  764.          endif
    765. 

  765. !
    766. 

  766.          fluxin=0.
    767. 

  767.          do k=max_q,min_q,-1
    768. 

  768.             fluxout=rho(k)*vtrold(k)*qrz(k)
    769. 

  769.             flux=(fluxin-fluxout)/rho(k)/dzw(k)
    770. 

  770.             tmpqrz=qrz(k)
    771. 

  771.             qrz(k)=qrz(k)+del_tv*flux
    772. 

  772.             fluxin=fluxout
    773. 

  773.          enddo
    774. 

  774.          if (min_q .eq. 1) then
    775. 

  775.             pptrain=pptrain+fluxin*del_tv
    776. 

  776.          else
    777. 

  777.             qrz(min_q-1)=qrz(min_q-1)+del_tv*  &
    778. 

  778.                           fluxin/rho(min_q-1)/dzw(min_q-1)
    779. 

  779.          endif
    780. 

  780. !
    781. 

  781.       else
    782. 

  782.          notlast=.false.
    783. 

  783.       endif
    784. 

  784.     ENDDO
    785. 

  785. 

  786. !
    787. 

  787. !-- snow
    788. 

  788. !
    789. 

  789.     t_del_tv=0.
    790. 

  790.     del_tv=dtb
    791. 

  791.     notlast=.true.
    792. 

  792. 

  793.     DO while (notlast)
    794. 

  794. !
    795. 

  795.       min_q=kte
    796. 

  796.       max_q=kts-1
    797. 

  797. !
    798. 

  798.       do k=kts,kte-1
    799. 

  799.          if (qsz(k) .gt. 1.0e-8) then
    800. 

  800.             min_q=min0(min_q,k)
    801. 

  801.             max_q=max0(max_q,k)
    802. 

  802.             tmp1=sqrt(pi*rhosnow*xnos/rho(k)/qsz(k))
    803. 

  803.             tmp1=sqrt(tmp1)
    804. 

  804.             vtsold(k)=o6*constc*gamdp4*sqrho(k)/tmp1**constd
    805. 

  805.             if (k .eq. 1) then
    806. 

  806.                del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtsold(k))
    807. 

  807.             else
    808. 

  808.                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtsold(k))
    809. 

  809.             endif
    810. 

  810.          else
    811. 

  811.             vtsold(k)=0.
    812. 

  812.          endif
    813. 

  813.       enddo
    814. 

  814. 

  815.       if (max_q .ge. min_q) then
    816. 

  816. !
    817. 

  817. !
    818. 

  818. !- Check if the summation of the small delta t >=  big delta t
    819. 

  819. !             (t_del_tv)          (del_tv)             (dtb)
    820. 

  820. 

  821.          t_del_tv=t_del_tv+del_tv
    822. 

  822. 

  823.          if ( t_del_tv .ge. dtb ) then
    824. 

  824.               notlast=.false.
    825. 

  825.               del_tv=dtb+del_tv-t_del_tv
    826. 

  826.          endif
    827. 

  827. !
    828. 

  828.          fluxin=0.
    829. 

  829.          do k=max_q,min_q,-1
    830. 

  830.             fluxout=rho(k)*vtsold(k)*qsz(k)
    831. 

  831.             flux=(fluxin-fluxout)/rho(k)/dzw(k)
    832. 

  832.             qsz(k)=qsz(k)+del_tv*flux
    833. 

  833.             qsz(k)=amax1(0.,qsz(k))
    834. 

  834.             fluxin=fluxout
    835. 

  835.          enddo
    836. 

  836.          if (min_q .eq. 1) then
    837. 

  837.             pptsnow=pptsnow+fluxin*del_tv
    838. 

  838.          else
    839. 

  839.             qsz(min_q-1)=qsz(min_q-1)+del_tv*  &
    840. 

  840.                          fluxin/rho(min_q-1)/dzw(min_q-1)
    841. 

  841.          endif
    842. 

  842. !
    843. 

  843.       else
    844. 

  844.          notlast=.false.
    845. 

  845.       endif
    846. 

  846. 

  847.     ENDDO
    848. 

  848. !
    849. 

  849. !-- grauupel
    850. 

  850. !
    851. 

  851.     t_del_tv=0.
    852. 

  852.     del_tv=dtb
    853. 

  853.     notlast=.true.
    854. 

  854. !
    855. 

  855.     DO while (notlast)
    856. 

  856. !
    857. 

  857.       min_q=kte
    858. 

  858.       max_q=kts-1
    859. 

  859. !
    860. 

  860.       do k=kts,kte-1
    861. 

  861.          if (qgz(k) .gt. 1.0e-8) then
    862. 

  862.             min_q=min0(min_q,k)
    863. 

  863.             max_q=max0(max_q,k)
    864. 

  864.             tmp1=sqrt(pi*rhograul*xnog/rho(k)/qgz(k))
    865. 

  865.             tmp1=sqrt(tmp1)
    866. 

  866.             term0=sqrt(4.*grav*rhograul*0.33334/rho(k)/cdrag)
    867. 

  867.             vtgold(k)=o6*gam4pt5*term0*sqrt(1./tmp1)
    868. 

  868.             if (k .eq. 1) then
    869. 

  869.                del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtgold(k))
    870. 

  870.             else
    871. 

  871.                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtgold(k))
    872. 

  872.             endif
    873. 

  873.          else
    874. 

  874.             vtgold(k)=0.
    875. 

  875.          endif
    876. 

  876.       enddo
    877. 

  877. 

  878.       if (max_q .ge. min_q) then
    879. 

  879. !
    880. 

  880. !
    881. 

  881. !- Check if the summation of the small delta t >=  big delta t
    882. 

  882. !             (t_del_tv)          (del_tv)             (dtb)
    883. 

  883. 

  884.          t_del_tv=t_del_tv+del_tv
    885. 

  885. 

  886.          if ( t_del_tv .ge. dtb ) then
    887. 

  887.               notlast=.false.
    888. 

  888.               del_tv=dtb+del_tv-t_del_tv
    889. 

  889.          endif
    890. 

  890. 

  891. !
    892. 

  892.          fluxin=0.
    893. 

  893.          do k=max_q,min_q,-1
    894. 

  894.             fluxout=rho(k)*vtgold(k)*qgz(k)
    895. 

  895.             flux=(fluxin-fluxout)/rho(k)/dzw(k)
    896. 

  896.             qgz(k)=qgz(k)+del_tv*flux
    897. 

  897.             qgz(k)=amax1(0.,qgz(k))
    898. 

  898.             fluxin=fluxout
    899. 

  899.          enddo
    900. 

  900.          if (min_q .eq. 1) then
    901. 

  901.             pptgraul=pptgraul+fluxin*del_tv
    902. 

  902.          else
    903. 

  903.             qgz(min_q-1)=qgz(min_q-1)+del_tv*  &
    904. 

  904.                          fluxin/rho(min_q-1)/dzw(min_q-1)
    905. 

  905.          endif
    906. 

  906. !
    907. 

  907.       else
    908. 

  908.          notlast=.false.
    909. 

  909.       endif
    910. 

  910. !
    911. 

  911.    ENDDO
    912. 

  912. 

  913. !
    914. 

  914. !-- cloud ice  (03/21/02) follow Vaughan T.J. Phillips at GFDL
    915. 

  915. !
    916. 

  916.     t_del_tv=0.
    917. 

  917.     del_tv=dtb
    918. 

  918.     notlast=.true.
    919. 

  919. !
    920. 

  920.     DO while (notlast)
    921. 

  921. !
    922. 

  922.       min_q=kte
    923. 

  923.       max_q=kts-1
    924. 

  924. !
    925. 

  925.       do k=kts,kte-1
    926. 

  926.          if (qiz(k) .gt. 1.0e-8) then
    927. 

  927.             min_q=min0(min_q,k)
    928. 

  928.             max_q=max0(max_q,k)
    929. 

  929.             vtiold(k)= 3.29 * (rho(k)* qiz(k))** 0.16  ! Heymsfield and Donner
    930. 

  930.             if (k .eq. 1) then
    931. 

  931.                del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtiold(k))
    932. 

  932.             else
    933. 

  933.                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtiold(k))
    934. 

  934.             endif
    935. 

  935.          else
    936. 

  936.             vtiold(k)=0.
    937. 

  937.          endif
    938. 

  938.       enddo
    939. 

  939. 

  940.       if (max_q .ge. min_q) then
    941. 

  941. !
    942. 

  942. !
    943. 

  943. !- Check if the summation of the small delta t >=  big delta t
    944. 

  944. !             (t_del_tv)          (del_tv)             (dtb)
    945. 

  945. 

  946.          t_del_tv=t_del_tv+del_tv
    947. 

  947. 

  948.          if ( t_del_tv .ge. dtb ) then
    949. 

  949.               notlast=.false.
    950. 

  950.               del_tv=dtb+del_tv-t_del_tv
    951. 

  951.          endif
    952. 

  952. 

  953.          fluxin=0.
    954. 

  954.          do k=max_q,min_q,-1
    955. 

  955.             fluxout=rho(k)*vtiold(k)*qiz(k)
    956. 

  956.             flux=(fluxin-fluxout)/rho(k)/dzw(k)
    957. 

  957.             qiz(k)=qiz(k)+del_tv*flux
    958. 

  958.             qiz(k)=amax1(0.,qiz(k))
    959. 

  959.             fluxin=fluxout
    960. 

  960.          enddo
    961. 

  961.          if (min_q .eq. 1) then
    962. 

  962.             pptice=pptice+fluxin*del_tv
    963. 

  963.          else
    964. 

  964.             qiz(min_q-1)=qiz(min_q-1)+del_tv*  &
    965. 

  965.                          fluxin/rho(min_q-1)/dzw(min_q-1)
    966. 

  966.          endif
    967. 

  967. !
    968. 

  968.       else
    969. 

  969.          notlast=.false.
    970. 

  970.       endif
    971. 

  971. !
    972. 

  972.    ENDDO
    973. 

  973.    do k=kts,kte-1                         !sg beg
    974. 

  974.       precrz(k)=rho(k)*vtrold(k)*qrz(k)
    975. 

  975.       preciz(k)=rho(k)*vtiold(k)*qiz(k)
    976. 

  976.       precsz(k)=rho(k)*vtsold(k)*qsz(k)
    977. 

  977.       precgz(k)=rho(k)*vtgold(k)*qgz(k)
    978. 

  978.    enddo                                  !sg end
    979. 

  979.    precrz(kte)=0. !wig - top level never set for vtXold vars
    980. 

  980.    preciz(kte)=0. !wig
    981. 

  981.    precsz(kte)=0. !wig
    982. 

  982.    precgz(kte)=0. !wig
    983. 

  983.    
    984. 

  984. 

  985. ! Microphysics processes
    986. 

  986. !
    987. 

  987.       DO 2000 k=kts,kte
    988. 

  988. !
    989. 

  989.          qvzodt(k)=amax1( 0.0,odtb*qvz(k) )
    990. 

  990.          qlzodt(k)=amax1( 0.0,odtb*qlz(k) )
    991. 

  991.          qizodt(k)=amax1( 0.0,odtb*qiz(k) )
    992. 

  992.          qszodt(k)=amax1( 0.0,odtb*qsz(k) )
    993. 

  993.          qrzodt(k)=amax1( 0.0,odtb*qrz(k) )
    994. 

  994.          qgzodt(k)=amax1( 0.0,odtb*qgz(k) )
    995. 

  995. !***********************************************************************
    996. 

  996. !*****   diagnose mixing ratios (qrz,qsz), terminal                *****
    997. 

  997. !*****   velocities (vtr,vts), and slope parameters in size        *****
    998. 

  998. !*****   distribution(xlambdar,xlambdas) of rain and snow          *****
    999. 

  999. !*****   follows Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7)   *****
    1000. 

  1000. !***********************************************************************
    1001. 

  1001. !
    1002. 

  1002. !**** assuming no cloud water can exist in the top two levels due to
    1003. 

  1003. !**** radiation consideration
    1004. 

  1004. !
    1005. 

  1005. !!  if
    1006. 

  1006. !!     unsaturated,
    1007. 

  1007. !!     no cloud water, rain, ice, snow and graupel
    1008. 

  1008. !!  then
    1009. 

  1009. !!     skip these processes and jump to line 2000
    1010. 

  1010. !
    1011. 

  1011. !
    1012. 

  1012.         tmp=qiz(k)+qlz(k)+qsz(k)+qrz(k)+qgz(k)*gindex
    1013. 

  1013.         if( qvz(k)+qlz(k)+qiz(k) .lt. qsiz(k)  &
    1014. 

  1014.             .and. tmp .eq. 0.0 ) go to 2000
    1015. 

  1015. 

  1016. !! calculate terminal velocity of rain
    1017. 

  1017. !
    1018. 

  1018.         if (qrz(k) .gt. 1.0e-8) then
    1019. 

  1019.             tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k))
    1020. 

  1020.             xlambdar(k)=sqrt(tmp1)
    1021. 

  1021.             olambdar(k)=1.0/xlambdar(k)
    1022. 

  1022.             vtrold(k)=o6*consta*gambp4*sqrho(k)*olambdar(k)**constb
    1023. 

  1023.         else
    1024. 

  1024.             vtrold(k)=0.
    1025. 

  1025.             olambdar(k)=0.
    1026. 

  1026.         endif
    1027. 

  1027. !
    1028. 

  1028. !       if (qrz(k) .gt. 1.0e-12) then
    1029. 

  1029.         if (qrz(k) .gt. 1.0e-8) then
    1030. 

  1030.             tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k))
    1031. 

  1031.             xlambdar(k)=sqrt(tmp1)
    1032. 

  1032.             olambdar(k)=1.0/xlambdar(k)
    1033. 

  1033.             vtr(k)=o6*consta*gambp4*sqrho(k)*olambdar(k)**constb
    1034. 

  1034.         else
    1035. 

  1035.             vtr(k)=0.
    1036. 

  1036.             olambdar(k)=0.
    1037. 

  1037.         endif
    1038. 

  1038. !
    1039. 

  1039. !! calculate terminal velocity of snow
    1040. 

  1040. !
    1041. 

  1041.         if (qsz(k) .gt. 1.0e-8) then
    1042. 

  1042.             tmp1=sqrt(pi*rhosnow*xnos*orho(k)/qsz(k))
    1043. 

  1043.             xlambdas(k)=sqrt(tmp1)
    1044. 

  1044.             olambdas(k)=1.0/xlambdas(k)
    1045. 

  1045.             vtsold(k)=o6*constc*gamdp4*sqrho(k)*olambdas(k)**constd
    1046. 

  1046.         else
    1047. 

  1047.             vtsold(k)=0.
    1048. 

  1048.             olambdas(k)=0.
    1049. 

  1049.         endif
    1050. 

  1050. !
    1051. 

  1051. !       if (qsz(k) .gt. 1.0e-12) then
    1052. 

  1052.         if (qsz(k) .gt. 1.0e-8) then
    1053. 

  1053.             tmp1=sqrt(pi*rhosnow*xnos*orho(k)/qsz(k))
    1054. 

  1054.             xlambdas(k)=sqrt(tmp1)
    1055. 

  1055.             olambdas(k)=1.0/xlambdas(k)
    1056. 

  1056.             vts(k)=o6*constc*gamdp4*sqrho(k)*olambdas(k)**constd
    1057. 

  1057.         else
    1058. 

  1058.             vts(k)=0.
    1059. 

  1059.             olambdas(k)=0.
    1060. 

  1060.         endif
    1061. 

  1061. !
    1062. 

  1062. !! calculate terminal velocity of graupel
    1063. 

  1063. !
    1064. 

  1064.         if (qgz(k) .gt. 1.0e-8) then
    1065. 

  1065.             tmp1=sqrt( pi*rhograul*xnog*orho(k)/qgz(k))
    1066. 

  1066.             xlambdag(k)=sqrt(tmp1)
    1067. 

  1067.             olambdag(k)=1.0/xlambdag(k)
    1068. 

  1068.             term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag)
    1069. 

  1069.             vtgold(k)=o6*gam4pt5*term0*sqrt(olambdag(k))
    1070. 

  1070.         else
    1071. 

  1071.             vtgold(k)=0.
    1072. 

  1072.             olambdag(k)=0.
    1073. 

  1073.         endif
    1074. 

  1074. !
    1075. 

  1075. !       if (qgz(k) .gt. 1.0e-12) then
    1076. 

  1076.         if (qgz(k) .gt. 1.0e-8) then
    1077. 

  1077.             tmp1=sqrt( pi*rhograul*xnog*orho(k)/qgz(k))
    1078. 

  1078.             xlambdag(k)=sqrt(tmp1)
    1079. 

  1079.             olambdag(k)=1.0/xlambdag(k)
    1080. 

  1080.             term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag)
    1081. 

  1081.             vtg(k)=o6*gam4pt5*term0*sqrt(olambdag(k))
    1082. 

  1082.         else
    1083. 

  1083.             vtg(k)=0.
    1084. 

  1084.             olambdag(k)=0.
    1085. 

  1085.         endif
    1086. 

  1086. !
    1087. 

  1087. !***********************************************************************
    1088. 

  1088. !*****  compute viscosity,difusivity,thermal conductivity, and    ******
    1089. 

  1089. !*****  Schmidt number                                            ******
    1090. 

  1090. !***********************************************************************
    1091. 

  1091. !c------------------------------------------------------------------
    1092. 

  1092. !c      viscmu: dynamic viscosity of air kg/m/s
    1093. 

  1093. !c      visc: kinematic viscosity of air = viscmu/rho (m2/s)
    1094. 

  1094. !c      avisc=1.49628e-6 kg/m/s=1.49628e-5 g/cm/s
    1095. 

  1095. !c      viscmu=1.718e-5 kg/m/s in RH
    1096. 

  1096. !c      diffwv: Diffusivity of water vapor in air
    1097. 

  1097. !c      adiffwv = 8.7602e-5 (8.794e-5 in MM5) kgm/s3
    1098. 

  1098. !c              = 8.7602 (8.794 in MM5) gcm/s3
    1099. 

  1099. !c      diffwv(k)=2.26e-5 m2/s
    1100. 

  1100. !c      schmidt: Schmidt number=visc/diffwv
    1101. 

  1101. !c      xka: thermal conductivity of air J/m/s/K (Kgm/s3/K)
    1102. 

  1102. !c      xka(k)=2.43e-2 J/m/s/K in RH
    1103. 

  1103. !c      axka=1.4132e3 (1.414e3 in MM5) m2/s2/k = 1.4132e7 cm2/s2/k
    1104. 

  1104. !c------------------------------------------------------------------
    1105. 

  1105. 

  1106.         viscmu(k)=avisc*tem(k)**1.5/(tem(k)+120.0)
    1107. 

  1107.         visc(k)=viscmu(k)*orho(k)
    1108. 

  1108.         diffwv(k)=adiffwv*tem(k)**1.81*oprez(k)
    1109. 

  1109.         schmidt(k)=visc(k)/diffwv(k)
    1110. 

  1110.         xka(k)=axka*viscmu(k)
    1111. 

  1111. 

  1112.         if (tem(k) .lt. 273.15) then
    1113. 

  1113. 

  1114. !
    1115. 

  1115. !***********************************************************************
    1116. 

  1116. !*********        snow production processes for T < 0 C       **********
    1117. 

  1117. !***********************************************************************
    1118. 

  1118. !c
    1119. 

  1119. !c (1) ICE CRYSTAL AGGREGATION TO SNOW (Psaut): Lin (21)
    1120. 

  1120. !c!    psaut=alpha1*(qi-qi0)
    1121. 

  1121. !c!    alpha1=1.0e-3*exp(0.025*(T-T0))
    1122. 

  1122. !c
    1123. 

  1123. !          alpha1=1.0e-3*exp( 0.025*temcc(k) )
    1124. 

  1124. 

  1125.            alpha1=1.0e-3*exp( 0.025*temcc(k) )
    1126. 

  1126. !
    1127. 

  1127.            if(temcc(k) .lt. -20.0) then
    1128. 

  1128.              tmp1=-7.6+4.0*exp( -0.2443e-3*(abs(temcc(k))-20)**2.455 )
    1129. 

  1129.              qic=1.0e-3*exp(tmp1)*orho(k)
    1130. 

  1130.            else
    1131. 

  1131.              qic=qi0
    1132. 

  1132.            end if
    1133. 

  1133. !testing
    1134. 

  1134. !          tmp1=amax1( 0.0,alpha1*(qiz(k)-qic) )
    1135. 

  1135. !          psaut(k)=amin1( tmp1,qizodt(k) )
    1136. 

  1136. 

  1137.            tmp1=odtb*(qiz(k)-qic)*(1.0-exp(-alpha1*dtb))
    1138. 

  1138.            psaut(k)=amax1( 0.0,tmp1 )
    1139. 

  1139. 

  1140. !c
    1141. 

  1141. !c (2) BERGERON PROCESS TRANSFER OF CLOUD WATER TO SNOW (Psfw)
    1142. 

  1142. !c     this process only considered when -31 C < T < 0 C
    1143. 

  1143. !c     Lin (33) and Hsie (17)
    1144. 

  1144. !c
    1145. 

  1145. !c!
    1146. 

  1146. !c!    parama1 and parama2 functions must be user supplied
    1147. 

  1147. !c!
    1148. 

  1148. 

  1149. ! testing
    1150. 

  1150.           if( qlz(k) .gt. 1.0e-10 ) then
    1151. 

  1151.             temc1=amax1(-30.99,temcc(k))
    1152. 

  1152. !           print*,'temc1',temc1,qlz(k)
    1153. 

  1153.             a1=parama1( temc1 )
    1154. 

  1154.             a2=parama2( temc1 )
    1155. 

  1155.             tmp1=1.0-a2
    1156. 

  1156. !! change unit from cgs to mks
    1157. 

  1157.             a1=a1*0.001**tmp1
    1158. 

  1158. !c!   dtberg is the time needed for a crystal to grow from 40 to 50 um
    1159. 

  1159. !c !  odtberg=1.0/dtberg
    1160. 

  1160.             odtberg=(a1*tmp1)/(xmi50**tmp1-xmi40**tmp1)
    1161. 

  1161. !
    1162. 

  1162. !c!   compute terminal velocity of a 50 micron ice cystal
    1163. 

  1163. !
    1164. 

  1164.             vti50=constc*di50**constd*sqrho(k)
    1165. 

  1165. !
    1166. 

  1166.             eiw=1.0
    1167. 

  1167.             save1=a1*xmi50**a2
    1168. 

  1168.             save2=0.25*pi*eiw*rho(k)*di50*di50*vti50
    1169. 

  1169. !
    1170. 

  1170.             tmp2=( save1 + save2*qlz(k) )
    1171. 

  1171. !
    1172. 

  1172. !! maximum number of 50 micron crystals limited by the amount
    1173. 

  1173. !!  of supercool water
    1174. 

  1174. !
    1175. 

  1175.             xni50mx=qlzodt(k)/tmp2
    1176. 

  1176. !
    1177. 

  1177. !!   number of 50 micron crystals produced
    1178. 

  1178. !
    1179. 

  1179. !
    1180. 

  1180.             xni50=qiz(k)*( 1.0-exp(-dtb*odtberg) )/xmi50
    1181. 

  1181.             xni50=amin1(xni50,xni50mx)
    1182. 

  1182. !
    1183. 

  1183.             tmp3=odtb*tmp2/save2*( 1.0-exp(-save2*xni50*dtb) )
    1184. 

  1184.             psfw(k)=amin1( tmp3,qlzodt(k) )
    1185. 

  1185. !testing
    1186. 

  1186. !           psfw(k)=0.
    1187. 

  1187. 

  1188. !0915     if( temcc(k).gt.-30.99 ) then
    1189. 

  1189. !0915        a1=parama1( temcc(k) )
    1190. 

  1190. !0915        a2=parama2( temcc(k) )
    1191. 

  1191. !0915        tmp1=1.0-a2
    1192. 

  1192. !!     change unit from cgs to mks
    1193. 

  1193. !0915        a1=a1*0.001**tmp1
    1194. 

  1194. 

  1195. !c!    dtberg is the time needed for a crystal to grow from 40 to 50 um
    1196. 

  1196. !c!    odtberg=1.0/dtberg
    1197. 

  1197. !0915        odtberg=(a1*tmp1)/(xmi50**tmp1-xmi40**tmp1)
    1198. 

  1198. 

  1199. !c!    number of 50 micron crystals produced
    1200. 

  1200. !0915        xni50=qiz(k)*dtb*odtberg/xmi50
    1201. 

  1201. 

  1202. !c!    need to calculate the terminal velocity of a 50 micron
    1203. 

  1203. !c!    ice cystal
    1204. 

  1204. !0915        vti50=constc*di50**constd*sqrho(k)
    1205. 

  1205. !0915        eiw=1.0
    1206. 

  1206. !0915        tmp2=xni50*( a1*xmi50**a2 + &
    1207. 

  1207. !0915             0.25*qlz(k)*pi*eiw*rho(k)*di50*di50*vti50 )
    1208. 

  1208. !0915        psfw(k)=amin1( tmp2,qlzodt(k) )
    1209. 

  1209. !0915        psfw(k)=0.
    1210. 

  1210. !c
    1211. 

  1211. !c (3) REDUCTION OF CLOUD ICE BY BERGERON PROCESS (Psfi): Lin (34)
    1212. 

  1212. !c     this process only considered when -31 C < T < 0 C
    1213. 

  1213. !c
    1214. 

  1214.             tmp1=xni50*xmi50-psfw(k)
    1215. 

  1215.             psfi(k)=amin1(tmp1,qizodt(k))
    1216. 

  1216. ! testing
    1217. 

  1217. !           psfi(k)=0.
    1218. 

  1218.           end if
    1219. 

  1219. !
    1220. 

  1220. 

  1221. !0915        tmp1=qiz(k)*odtberg
    1222. 

  1222. !0915        psfi(k)=amin1(tmp1,qizodt(k))
    1223. 

  1223. ! testing
    1224. 

  1224. !0915        psfi(k)=0.
    1225. 

  1225. !0915     end if
    1226. 

  1226. !
    1227. 

  1227.           if(qrz(k) .le. 0.0) go to 1000
    1228. 

  1228. !
    1229. 

  1229. ! Processes (4) and (5) only need when qrz > 0.0
    1230. 

  1230. !
    1231. 

  1231. !c
    1232. 

  1232. !c (4) CLOUD ICE ACCRETION BY RAIN (Praci): Lin (25)
    1233. 

  1233. !c     may produce snow or graupel
    1234. 

  1234. !c
    1235. 

  1235.           eri=1.0
    1236. 

  1236. !0915     tmp1=qiz(k)*pio4*eri*xnor*consta*sqrho(k)
    1237. 

  1237. !0915     tmp2=tmp1*gambp3*olambdar(k)**bp3
    1238. 

  1238. !0915     praci(k)=amin1( tmp2,qizodt(k) )
    1239. 

  1239. 

  1240.           save1=pio4*eri*xnor*consta*sqrho(k)
    1241. 

  1241.           tmp1=save1*gambp3*olambdar(k)**bp3
    1242. 

  1242.           praci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) )
    1243. 

  1243. 

  1244. !c
    1245. 

  1245. !c (5) RAIN ACCRETION BY CLOUD ICE (Piacr): Lin (26)
    1246. 

  1246. !c
    1247. 

  1247. !0915     tmp2=tmp1*rho(k)*pio6*rhowater*gambp6*oxmi* &
    1248. 

  1248. !0915              olambdar(k)**bp6
    1249. 

  1249. !0915     piacr(k)=amin1( tmp2,qrzodt(k) )
    1250. 

  1250. 

  1251.           tmp2=qiz(k)*save1*rho(k)*pio6*rhowater*gambp6*oxmi* &
    1252. 

  1252.                    olambdar(k)**bp6
    1253. 

  1253.           piacr(k)=amin1( tmp2,qrzodt(k) )
    1254. 

  1254. 

  1255. !
    1256. 

  1256. 1000      continue
    1257. 

  1257. !
    1258. 

  1258.           if(qsz(k) .le. 0.0) go to 1200
    1259. 

  1259. !
    1260. 

  1260. ! Compute the following processes only when qsz > 0.0
    1261. 

  1261. !
    1262. 

  1262. !c
    1263. 

  1263. !c (6) ICE CRYSTAL ACCRETION BY SNOW (Psaci): Lin (22)
    1264. 

  1264. !c
    1265. 

  1265.           esi=exp( 0.025*temcc(k) )
    1266. 

  1266.           save1=pio4*xnos*constc*gamdp3*sqrho(k)* &
    1267. 

  1267.                olambdas(k)**dp3
    1268. 

  1268.           tmp1=esi*save1
    1269. 

  1269.           psaci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) )
    1270. 

  1270. 

  1271. !0915     tmp1=pio4*xnos*constc*gamdp3*sqrho(k)* &
    1272. 

  1272. !0915          olambdas(k)**dp3
    1273. 

  1273. !0915     tmp2=qiz(k)*esi*tmp1
    1274. 

  1274. !0915     psaci(k)=amin1( tmp2,qizodt(k) )
    1275. 

  1275. !c
    1276. 

  1276. !c (7) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
    1277. 

  1277. !c
    1278. 

  1278.           esw=1.0
    1279. 

  1279.           tmp1=esw*save1
    1280. 

  1280.           psacw(k)=qlzodt(K)*( 1.0-exp(-tmp1*dtb) )
    1281. 

  1281. 

  1282. !0915     tmp2=qlz(k)*esw*tmp1
    1283. 

  1283. !0915     psacw(k)=amin1( tmp2,qlzodt(k) )
    1284. 

  1284. !c
    1285. 

  1285. !c (8) DEPOSITION/SUBLIMATION OF SNOW (Psdep/Pssub): Lin (31)
    1286. 

  1286. !c     includes consideration of ventilation effect
    1287. 

  1287. !c
    1288. 

  1288. !c     abi=2*pi*(Si-1)/rho/(A"+B")
    1289. 

  1289. !c
    1290. 

  1290.           tmpa=rvapor*xka(k)*tem(k)*tem(k)
    1291. 

  1291.           tmpb=xls*xls*rho(k)*qsiz(k)*diffwv(k)
    1292. 

  1292.           tmpc=tmpa*qsiz(k)*diffwv(k)
    1293. 

  1293.           abi=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
    1294. 

  1294. !
    1295. 

  1295. !c     vf1s,vf2s=ventilation factors for snow
    1296. 

  1296. !c     vf1s=0.78,vf2s=0.31 in LIN
    1297. 

  1297. !
    1298. 

  1298.           tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
    1299. 

  1299.           tmp2=abi*xnos*( vf1s*olambdas(k)*olambdas(k)+ &
    1300. 

  1300.                     vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
    1301. 

  1301.           tmp3=odtb*( qvz(k)-qsiz(k) )
    1302. 

  1302. !
    1303. 

  1303. !bloss: BEGIN
    1304. 

  1304. !the old implementation would give the wrong results if olambdas(k) ==0
    1305. 

  1305. !which can lead to positive pssub, i.e. tmp2=0 but tmp3> 0
    1306. 

  1306. !          if( tmp2 .le. 0.0) then
    1307. 

  1307.           if( tmp3 .le. 0.0) then
    1308. 

  1308.             tmp2=amax1( tmp2,tmp3)
    1309. 

  1309.             pssub(k)=amin1(0.,amax1( tmp2,-qszodt(k) ))
    1310. 

  1310. !bloss: END
    1311. 

  1311.             psdep(k)=0.0
    1312. 

  1312.           else
    1313. 

  1313.             psdep(k)=amin1( tmp2,tmp3 )
    1314. 

  1314.             pssub(k)=0.0
    1315. 

  1315.           end if
    1316. 

  1316. 

  1317. !0915     psdep(k)=amax1(0.0,tmp2)
    1318. 

  1318. !0915     pssub(k)=amin1(0.0,tmp2)
    1319. 

  1319. !0915     pssub(k)=amax1( pssub(k),-qszodt(k) )
    1320. 

  1320. !
    1321. 

  1321.           if(qrz(k) .le. 0.0) go to 1200
    1322. 

  1322. !
    1323. 

  1323. ! Compute processes (9) and (10) only when qsz > 0.0 and qrz > 0.0
    1324. 

  1324. !
    1325. 

  1325. !c
    1326. 

  1326. !c (9) ACCRETION OF SNOW BY RAIN (Pracs): Lin (27)
    1327. 

  1327. !c
    1328. 

  1328.           esr=1.0
    1329. 

  1329.           tmpa=olambdar(k)*olambdar(k)
    1330. 

  1330.           tmpb=olambdas(k)*olambdas(k)
    1331. 

  1331.           tmpc=olambdar(k)*olambdas(k)
    1332. 

  1332.           tmp1=pi*pi*esr*xnor*xnos*abs( vtr(k)-vts(k) )*orho(k)
    1333. 

  1333.           tmp2=tmpb*tmpb*olambdar(k)*(5.0*tmpb+2.0*tmpc+0.5*tmpa)
    1334. 

  1334.           tmp3=tmp1*rhosnow*tmp2
    1335. 

  1335.           pracs(k)=amin1( tmp3,qszodt(k) )
    1336. 

  1336. !c
    1337. 

  1337. !c (10)  ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
    1338. 

  1338. !c
    1339. 

  1339.           tmp3=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
    1340. 

  1340.           tmp4=tmp1*rhowater*tmp3
    1341. 

  1341.           psacr(k)=amin1( tmp4,qrzodt(k) )
    1342. 

  1342. !
    1343. 

  1343. 1200      continue
    1344. 

  1344. !
    1345. 

  1345.         else
    1346. 

  1346. !
    1347. 

  1347. !***********************************************************************
    1348. 

  1348. !*********        snow production processes for T > 0 C       **********
    1349. 

  1349. !***********************************************************************
    1350. 

  1350. !
    1351. 

  1351.          if (qsz(k) .le. 0.0) go to 1400
    1352. 

  1352. !c
    1353. 

  1353. !c (1) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
    1354. 

  1354. !c
    1355. 

  1355.             esw=1.0
    1356. 

  1356. 

  1357.             tmp1=esw*pio4*xnos*constc*gamdp3*sqrho(k)* &
    1358. 

  1358.                  olambdas(k)**dp3
    1359. 

  1359.             psacw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
    1360. 

  1360. 

  1361. !0915       tmp1=pio4*xnos*constc*gamdp3*sqrho(k)* &
    1362. 

  1362. !0915            olambdas(k)**dp3
    1363. 

  1363. !0915       tmp2=qlz(k)*esw*tmp1
    1364. 

  1364. !0915       psacw(k)=amin1( tmp2,qlzodt(k) )
    1365. 

  1365. !c
    1366. 

  1366. !c (2) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
    1367. 

  1367. !c
    1368. 

  1368.             esr=1.0
    1369. 

  1369.             tmpa=olambdar(k)*olambdar(k)
    1370. 

  1370.             tmpb=olambdas(k)*olambdas(k)
    1371. 

  1371.             tmpc=olambdar(k)*olambdas(k)
    1372. 

  1372.             tmp1=pi*pi*esr*xnor*xnos*abs( vtr(k)-vts(k) )*orho(k)
    1373. 

  1373.             tmp2=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
    1374. 

  1374.             tmp3=tmp1*rhowater*tmp2
    1375. 

  1375.             psacr(k)=amin1( tmp3,qrzodt(k) )
    1376. 

  1376. !c
    1377. 

  1377. !c (3) MELTING OF SNOW (Psmlt): Lin (32)
    1378. 

  1378. !c     Psmlt is negative value
    1379. 

  1379. !
    1380. 

  1380.             delrs=rs0(k)-qvz(k)
    1381. 

  1381.             term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- &
    1382. 

  1382.                   xka(k)*temcc(k) )
    1383. 

  1383.             tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
    1384. 

  1384.             tmp2=xnos*( vf1s*olambdas(k)*olambdas(k)+  &
    1385. 

  1385.                  vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
    1386. 

  1386.             tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( psacw(k)+psacr(k) )
    1387. 

  1387.             tmp4=amin1(0.0,tmp3)
    1388. 

  1388.             psmlt(k)=amax1( tmp4,-qszodt(k) )
    1389. 

  1389. !c
    1390. 

  1390. !c (4) EVAPORATION OF MELTING SNOW (Psmltevp): HR (A27)
    1391. 

  1391. !c     but use Lin et al. coefficience
    1392. 

  1392. !c     Psmltevp is a negative value
    1393. 

  1393. !c
    1394. 

  1394.             tmpa=rvapor*xka(k)*tem(k)*tem(k)
    1395. 

  1395.             tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
    1396. 

  1396.             tmpc=tmpa*qswz(k)*diffwv(k)
    1397. 

  1397.             tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb )
    1398. 

  1398. 

  1399. !      abr=2.0*pi*(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb)
    1400. 

  1400. 

  1401.             abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
    1402. 

  1402. !
    1403. 

  1403. !**** allow evaporation to occur when RH less than 90%
    1404. 

  1404. !**** here not using 100% because the evaporation cooling
    1405. 

  1405. !**** of temperature is not taking into account yet; hence,
    1406. 

  1406. !**** the qsw value is a little bit larger. This will avoid
    1407. 

  1407. !**** evaporation can generate cloud.
    1408. 

  1408. !
    1409. 

  1409. !c    vf1s,vf2s=ventilation factors for snow
    1410. 

  1410. !c    vf1s=0.78,vf2s=0.31 in LIN
    1411. 

  1411. !
    1412. 

  1412.             tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
    1413. 

  1413.             tmp2=abr*xnos*( vf1s*olambdas(k)*olambdas(k)+  &
    1414. 

  1414.                  vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
    1415. 

  1415.             tmp3=amin1(0.0,tmp2)
    1416. 

  1416.             tmp3=amax1( tmp3,tmpd )
    1417. 

  1417.             psmltevp(k)=amax1( tmp3,-qszodt(k) )
    1418. 

  1418. 1400     continue
    1419. 

  1419. !
    1420. 

  1420.         end if
    1421. 

  1421. 

  1422. !***********************************************************************
    1423. 

  1423. !*********           rain production processes                **********
    1424. 

  1424. !***********************************************************************
    1425. 

  1425. !
    1426. 

  1426. !c
    1427. 

  1427. !c (1) AUTOCONVERSION OF RAIN (Praut): RH
    1428. 

  1428. !sg: begin
    1429. 

  1429.         if(flag_qndrop)then
    1430. 

  1430.            if( qndropz(k) >= 1. ) then
    1431. 

  1431. !         Liu et al. autoconversion scheme
    1432. 

  1432.               rhocgs=rho(k)*1.e-3
    1433. 

  1433.               liqconc=rhocgs*qlz(k)   ! (kg/kg) to (g/cm3)
    1434. 

  1434.               capn=1.0e-3*rhocgs*qndropz(k)   ! (#/kg) to (#/cm3)
    1435. 

  1435. !         rate function
    1436. 

  1436.               if(liqconc.gt.1.e-10)then
    1437. 

  1437.                  p0=(kappa*beta/capn)*(liqconc*liqconc*liqconc)
    1438. 

  1438.                  xc=9.7d-17*capn*sqrt(capn)/(liqconc*liqconc)
    1439. 

  1439. !         Calculate autoconversion rate (g/g/s)
    1440. 

  1440.                  if(xc.lt.10.)then
    1441. 

  1441.                     praut(k)=(p0/rhocgs) * ( 0.5d0*(xc*xc+2*xc+2.0d0)* &
    1442. 

  1442.                          (1.0d0+xc)*exp(-2.0d0*xc) )
    1443. 

  1443.                  endif
    1444. 

  1444.               endif
    1445. 

  1445.            endif
    1446. 

  1446.         else
    1447. 

  1447. !sg: end
    1448. 

  1448. !c          araut=afa*rho
    1449. 

  1449. !c          afa=0.001 Rate coefficient for autoconvergence
    1450. 

  1450. !c
    1451. 

  1451. !c          araut=1.0e-3
    1452. 

  1452. !c
    1453. 

  1453.             araut=0.001
    1454. 

  1454. !testing
    1455. 

  1455. !           tmp1=amax1( 0.0,araut*(qlz(k)-ql0) )
    1456. 

  1456. !           praut(k)=amin1( tmp1,qlzodt(k) )
    1457. 

  1457.             tmp1=odtb*(qlz(k)-ql0)*( 1.0-exp(-araut*dtb) )
    1458. 

  1458.             praut(k)=amax1( 0.0,tmp1 )
    1459. 

  1459.         endif !sg
    1460. 

  1460. 

  1461. !c
    1462. 

  1462. !c (2) ACCRETION OF CLOUD WATER BY RAIN (Pracw): Lin (51)
    1463. 

  1463. !c
    1464. 

  1464.          erw=1.0
    1465. 

  1465. !        tmp1=qlz(k)*pio4*erw*xnor*consta*sqrho(k)
    1466. 

  1466. !        tmp2=tmp1*gambp3*olambdar(k)**bp3
    1467. 

  1467. !        pracw(k)=amin1( tmp2,qlzodt(k) )
    1468. 

  1468. 

  1469.         tmp1=pio4*erw*xnor*consta*sqrho(k)* &
    1470. 

  1470.              gambp3*olambdar(k)**bp3
    1471. 

  1471.         pracw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
    1472. 

  1472. 

  1473. !c
    1474. 

  1474. !c (3) EVAPORATION OF RAIN (Prevp): Lin (52)
    1475. 

  1475. !c     Prevp is negative value
    1476. 

  1476. !c
    1477. 

  1477. !c     Sw=qvoqsw : saturation ratio
    1478. 

  1478. !c
    1479. 

  1479.          tmpa=rvapor*xka(k)*tem(k)*tem(k)
    1480. 

  1480.          tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
    1481. 

  1481.          tmpc=tmpa*qswz(k)*diffwv(k)
    1482. 

  1482. !bloss: BEGIN
    1483. 

  1483. !         tmpd=amin1(0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb)
    1484. 

  1484. 

  1485. ! set max allowed evaporation to 90% of the amount that
    1486. 

  1486. !   would induce saturation wrt liquid in a single step.
    1487. 

  1487.          tmpd = qswz(k)*xlv/(rvapor*tem(k)**2) ! d(qsat_liq)/dT 
    1488. 

  1488.          tmpd = min( 0., 0.9*odtb*(qvz(k) + qlz(k) - qswz(k)) &
    1489. 

  1489.                           / (1. + xlvocp * tmpd) )
    1490. 

  1490. 

  1491.          abr=2.0*pi*(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb)
    1492. 

  1492. 

  1493. !         abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
    1494. 

  1494. !bloss: END
    1495. 

  1495. !
    1496. 

  1496. !c     vf1r,vf2r=ventilation factors for rain
    1497. 

  1497. !c     vf1r=0.78,vf2r=0.31 in RH, LIN  and MM5
    1498. 

  1498. !
    1499. 

  1499.          vf1r=0.78
    1500. 

  1500.          vf2r=0.31
    1501. 

  1501.          tmp1=consta*sqrho(k)*olambdar(k)**bp5/visc(k)
    1502. 

  1502.          tmp2=abr*xnor*( vf1r*olambdar(k)*olambdar(k)+  &
    1503. 

  1503.               vf2r*schmidt(k)**0.33334*gambp5o2*sqrt(tmp1) )
    1504. 

  1504.          tmp3=amin1( 0.0,tmp2 )
    1505. 

  1505.          tmp3=amax1( tmp3,tmpd )
    1506. 

  1506.          prevp(k)=amax1( tmp3,-qrzodt(k) )
    1507. 

  1507. 

  1508. !
    1509. 

  1509. !      if(iout .gt. 0) write(20,*)'tmp1,tmp2,tmp3=',tmp1,tmp2,tmp3
    1510. 

  1510. !      if(iout .gt. 0) write(20,*)'qlz,qiz,qrz=',qlz(k),qiz(k),qrz(k)
    1511. 

  1511. !      if(iout .gt. 0) write(20,*)'tem,qsz,qvz=',tem(k),qsz(k),qvz(k)
    1512. 

  1512. 

  1513. 

  1514. 

  1515. !     if (gindex .eq. 0.) goto 900
    1516. 

  1516. !
    1517. 

  1517.       if (tem(k) .lt. 273.15) then
    1518. 

  1518. !
    1519. 

  1519. !
    1520. 

  1520. !-- graupel
    1521. 

  1521. !***********************************************************************
    1522. 

  1522. !*********        graupel production processes for T < 0 C    **********
    1523. 

  1523. !***********************************************************************
    1524. 

  1524. !c
    1525. 

  1525. !c (1) AUTOCONVERSION OF SNOW TO FORM GRAUPEL (Pgaut): Lin (37)
    1526. 

  1526. !c     pgaut=alpha2*(qsz-qs0)
    1527. 

  1527. !c     qs0=6.0E-4
    1528. 

  1528. !c     alpha2=1.0e-3*exp(0.09*temcc(k))      Lin (38)
    1529. 

  1529. !
    1530. 

  1530.             alpha2=1.0e-3*exp(0.09*temcc(k))
    1531. 

  1531. !
    1532. 

  1532. 

  1533. ! testing
    1534. 

  1534. !           tmp1=alpha2*(qsz(k)-qs0)
    1535. 

  1535. !           tmp1=amax1(0.0,tmp1)
    1536. 

  1536. !           pgaut(k)=amin1( tmp1,qszodt(k) )
    1537. 

  1537. 

  1538.             tmp1=odtb*(qsz(k)-qs0)*(1.0-exp(-alpha2*dtb))
    1539. 

  1539.             pgaut(k)=amax1( 0.0,tmp1 )
    1540. 

  1540. 

  1541. !c
    1542. 

  1542. !c (2) FREEZING OF RAIN TO FORM GRAUPEL  (Pgfr): Lin (45)
    1543. 

  1543. !c     positive value
    1544. 

  1544. !c     Constant in Bigg freezing Aplume=Ap=0.66 /k
    1545. 

  1545. !c     Constant in raindrop freezing equ. Bplume=Bp=100./m/m/m/s
    1546. 

  1546. !
    1547. 

  1547. 

  1548.             if (qrz(k) .gt. 1.e-8 ) then
    1549. 

  1549.                Bp=100.
    1550. 

  1550.                Ap=0.66
    1551. 

  1551.                tmp1=olambdar(k)*olambdar(k)*olambdar(k)
    1552. 

  1552.                tmp2=20.*pi*pi*Bp*xnor*rhowater*orho(k)*  &
    1553. 

  1553.                     (exp(-Ap*temcc(k))-1.0)*tmp1*tmp1*olambdar(k)
    1554. 

  1554.                Pgfr(k)=amin1( tmp2,qrzodt(k) )
    1555. 

  1555.             else
    1556. 

  1556.                Pgfr(k)=0
    1557. 

  1557.             endif
    1558. 

  1558. 

  1559. !c
    1560. 

  1560. !c       if (qgz(k) = 0.0) skip the other step below about graupel
    1561. 

  1561. !c
    1562. 

  1562.          if (qgz(k) .eq. 0.0) goto 4000
    1563. 

  1563. 

  1564. !c
    1565. 

  1565. !c       Comparing Pgwet(wet process) and Pdry(dry process),
    1566. 

  1566. !c       we will pick up the small one.
    1567. 

  1567. !c
    1568. 

  1568. 

  1569. !c       ---------------
    1570. 

  1570. !c      | dry processes |
    1571. 

  1571. !c       ---------------
    1572. 

  1572. !c
    1573. 

  1573. !c (3)   ACCRETION OF CLOUD WATER BY GRAUPEL  (Pgacw): Lin (40)
    1574. 

  1574. !c       egw=1.0
    1575. 

  1575. !c       Cdrag=0.6 drag coefficients for hairstone
    1576. 

  1576. !c       constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
    1577. 

  1577. !c
    1578. 

  1578.          egw=1.0
    1579. 

  1579.          constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag)
    1580. 

  1580.          tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5
    1581. 

  1581.          tmp2=qlz(k)*egw*tmp1
    1582. 

  1582.          Pgacw(k)=amin1( tmp2,qlzodt(k) )
    1583. 

  1583. !c
    1584. 

  1584. !c (4)   ACCRETION OF ICE CRYSTAL BY GRAUPEL  (Pgaci): Lin (41)
    1585. 

  1585. !c       egi=1.   for wet growth
    1586. 

  1586. !c       egi=0.1  for dry growth
    1587. 

  1587. !c
    1588. 

  1588.          egi=0.1
    1589. 

  1589.          tmp2=qiz(k)*egi*tmp1
    1590. 

  1590.          pgaci(k)=amin1( tmp2,qizodt(k) )
    1591. 

  1591. 

  1592. 

  1593. !c
    1594. 

  1594. !c (5)   ACCRETION OF SNOW BY GRAUPEL  (Pgacs) : Lin (29)
    1595. 

  1595. !c       Compute processes (6) only when qsz > 0.0 and qgz > 0.0
    1596. 

  1596. !c
    1597. 

  1597.          egs=exp(0.09*temcc(k))
    1598. 

  1598.          tmpa=olambdas(k)*olambdas(k)
    1599. 

  1599.          tmpb=olambdag(k)*olambdag(k)
    1600. 

  1600.          tmpc=olambdas(k)*olambdag(k)
    1601. 

  1601.          tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k)
    1602. 

  1602.          tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
    1603. 

  1603.          tmp3=tmp1*egs*rhosnow*tmp2
    1604. 

  1604.          Pgacs(k)=amin1( tmp3,qszodt(k) )
    1605. 

  1605. 

  1606. 

  1607. !c
    1608. 

  1608. !c (6)   ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
    1609. 

  1609. !c       Compute processes (6) only when qrz > 0.0 and qgz > 0.0
    1610. 

  1610. !c       egr=1.
    1611. 

  1611. !c
    1612. 

  1612.          egr=1.
    1613. 

  1613.          tmpa=olambdar(k)*olambdar(k)
    1614. 

  1614.          tmpb=olambdag(k)*olambdag(k)
    1615. 

  1615.          tmpc=olambdar(k)*olambdag(k)
    1616. 

  1616.          tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k)
    1617. 

  1617.          tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
    1618. 

  1618.          tmp3=tmp1*egr*rhowater*tmp2
    1619. 

  1619.          pgacr(k)=amin1( tmp3,qrzodt(k) )
    1620. 

  1620. 

  1621. !c
    1622. 

  1622. !c (7)   Calculate total dry process effect Pdry(k)
    1623. 

  1623. !c
    1624. 

  1624.          Pdry(k)=Pgacw(k)+pgaci(k)+Pgacs(k)+pgacr(k)
    1625. 

  1625. 

  1626. !c       ---------------
    1627. 

  1627. !c      | wet processes |
    1628. 

  1628. !c       ---------------
    1629. 

  1629. !c
    1630. 

  1630. !c (3)   ACCRETION OF ICE CRYSTAL BY GRAUPEL  (Pgacip): Lin (41)
    1631. 

  1631. !c       egi=1.   for wet growth
    1632. 

  1632. !c       egi=0.1  for dry growth
    1633. 

  1633. !c
    1634. 

  1634.          tmp2=10.*pgaci(k)
    1635. 

  1635.          pgacip(k)=amin1( tmp2,qizodt(k) )
    1636. 

  1636. 

  1637. !c
    1638. 

  1638. !c (4)   ACCRETION OF SNOW BY GRAUPEL  ((Pgacsp) : Lin (29)
    1639. 

  1639. !c       Compute processes (6) only when qsz > 0.0 and qgz > 0.0
    1640. 

  1640. !c       egs=exp(0.09*(tem(k)-273.15))  when T < 273.15 k
    1641. 

  1641. !c
    1642. 

  1642.          tmp3=Pgacs(k)*1.0/egs
    1643. 

  1643.          Pgacsp(k)=amin1( tmp3,qszodt(k) )
    1644. 

  1644. 

  1645. !c
    1646. 

  1646. !c (5)   WET GROWTH OF GRAUPEL (Pgwet) : Lin (43)
    1647. 

  1647. !c       may involve Pgacs or Pgaci and
    1648. 

  1648. !c       must include PPgacw or Pgacr, or both.
    1649. 

  1649. !c       ( The amount of Pgacw which is not able
    1650. 

  1650. !c       to freeze is shed to rain. )
    1651. 

  1651.          IF(temcc(k).gt.-40.)THEN
    1652. 

  1652. 

  1653.              term0=constg*olambdag(k)**5.5/visc(k)
    1654. 

  1654. 

  1655. !c
    1656. 

  1656. !c    vf1s,vf2s=ventilation factors for graupel
    1657. 

  1657. !c    vf1s=0.78,vf2s=0.31 in LIN
    1658. 

  1658. !c    Cdrag=0.6  drag coefficient for hairstone
    1659. 

  1659. !c    constg2=vf1s*olambdag(k)*olambdag(k)+
    1660. 

  1660. !c            vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
    1661. 

  1661. 

  1662.              delrs=rs0(k)-qvz(k)
    1663. 

  1663.              tmp0=1./(xlf+cw*temcc(k))
    1664. 

  1664.              tmp1=2.*pi*xnog*(rho(k)*xlv*diffwv(k)*delrs-xka(k)*  &
    1665. 

  1665.                   temcc(k))*orho(k)*tmp0
    1666. 

  1666.              constg2=vf1s*olambdag(k)*olambdag(k)+  &
    1667. 

  1667.                      vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
    1668. 

  1668.              tmp3=tmp1*constg2+(Pgacip(k)+Pgacsp(k))*  &
    1669. 

  1669.                   (1-Ci*temcc(k)*tmp0)
    1670. 

  1670.              tmp3=amax1(0.0,tmp3)
    1671. 

  1671.              Pgwet(k)=amin1(tmp3,qlzodt(k)+qszodt(k)+qizodt(k) ) !bloss
    1672. 

  1672. 

  1673. !c
    1674. 

  1674. !c     Comparing Pgwet(wet process) and Pdry(dry process),
    1675. 

  1675. !c     we will apply the small one.
    1676. 

  1676. !c     if dry processes then delta4=1.0
    1677. 

  1677. !c     if wet processes then delta4=0.0
    1678. 

  1678. !
    1679. 

  1679.          if ( Pdry(k) .lt. Pgwet(k) ) then
    1680. 

  1680.             delta4=1.0
    1681. 

  1681.          else
    1682. 

  1682.             delta4=0.0
    1683. 

  1683.          endif
    1684. 

  1684.          ELSE
    1685. 

  1685.             delta4=1.0
    1686. 

  1686.          ENDIF
    1687. 

  1687. 

  1688. !c
    1689. 

  1689. !c
    1690. 

  1690. !c (6)   Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k)
    1691. 

  1691. !c       if Pgacrp(k) > 0. then some of the rain is frozen to hail
    1692. 

  1692. !c       if Pgacrp(k) < 0. then some of the cloud water collected
    1693. 

  1693. !c                            by the hail is unable to freeze and is
    1694. 

  1694. !c                            shed as rain.
    1695. 

  1695. !c
    1696. 

  1696.             Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k)
    1697. 

  1697. 

  1698. !c
    1699. 

  1699. !c (8)   DEPOSITION/SUBLIMATION OF GRAUPEL  (Pgdep/Pgsub): Lin (46)
    1700. 

  1700. !c       includes ventilation effect
    1701. 

  1701. !c       constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
    1702. 

  1702. !c       constg2=vf1s*olambdag(k)*olambdag(k)+
    1703. 

  1703. !c             vf2s*schmidt(k)**0.33334*gam2pt75*constg
    1704. 

  1704. !c
    1705. 

  1705. !c       abg=2*pi*(Si-1)/rho/(A"+B")
    1706. 

  1706. !c
    1707. 

  1707.             tmpa=rvapor*xka(k)*tem(k)*tem(k)
    1708. 

  1708.             tmpb=xls*xls*rho(k)*qsiz(k)*diffwv(k)
    1709. 

  1709.             tmpc=tmpa*qsiz(k)*diffwv(k)
    1710. 

  1710.             abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
    1711. 

  1711. !c
    1712. 

  1712. !c     vf1s,vf2s=ventilation factors for graupel
    1713. 

  1713. !c     vf1s=0.78,vf2s=0.31 in LIN
    1714. 

  1714. !c     Cdrag=0.6  drag coefficient for hairstone
    1715. 

  1715. !c
    1716. 

  1716.             term0=constg*olambdag(k)**5.5/visc(k)
    1717. 

  1717.             constg2=vf1s*olambdag(k)*olambdag(k)+  &
    1718. 

  1718.                     vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
    1719. 

  1719.             tmp2=abg*xnog*constg2
    1720. 

  1720.             pgdep(k)=amax1(0.0,tmp2)
    1721. 

  1721.             pgsub(k)=amin1(0.0,tmp2)
    1722. 

  1722.             pgsub(k)=amax1( pgsub(k),-qgzodt(k) )
    1723. 

  1723. 

  1724.  4000    continue
    1725. 

  1725.         else
    1726. 

  1726. !
    1727. 

  1727. !***********************************************************************
    1728. 

  1728. !*********      graupel production processes for T > 0 C      **********
    1729. 

  1729. !***********************************************************************
    1730. 

  1730. !
    1731. 

  1731. !c
    1732. 

  1732. !c (1) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40)
    1733. 

  1733. !c     egw=1.0
    1734. 

  1734. !c     Cdrag=0.6 drag coefficients for hairstone
    1735. 

  1735. !c     constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
    1736. 

  1736. 

  1737.             egw=1.0
    1738. 

  1738.             constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag)
    1739. 

  1739.             tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5
    1740. 

  1740.             tmp2=qlz(k)*egw*tmp1
    1741. 

  1741.             Pgacw(k)=amin1( tmp2,qlzodt(k) )
    1742. 

  1742. 

  1743. !c
    1744. 

  1744. !c (2) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
    1745. 

  1745. !c     Compute processes (5) only when qrz > 0.0 and qgz > 0.0
    1746. 

  1746. !c     egr=1.
    1747. 

  1747. !c
    1748. 

  1748.             egr=1.
    1749. 

  1749.             tmpa=olambdar(k)*olambdar(k)
    1750. 

  1750.             tmpb=olambdag(k)*olambdag(k)
    1751. 

  1751.             tmpc=olambdar(k)*olambdag(k)
    1752. 

  1752.             tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k)
    1753. 

  1753.             tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
    1754. 

  1754.             tmp3=tmp1*egr*rhowater*tmp2
    1755. 

  1755.             pgacr(k)=amin1( tmp3,qrzodt(k) )
    1756. 

  1756. 

  1757. 

  1758. !c
    1759. 

  1759. !c (3) GRAUPEL MELTING TO FORM RAIN (Pgmlt): Lin (47)
    1760. 

  1760. !c     Pgmlt is negative value
    1761. 

  1761. !c     constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
    1762. 

  1762. !c     constg2=vf1s*olambdag(k)*olambdag(k)+
    1763. 

  1763. !c             vf2s*schmidt(k)**0.33334*gam2pt75*constg
    1764. 

  1764. !c     Cdrag=0.6  drag coefficients for hairstone
    1765. 

  1765. !
    1766. 

  1766.             delrs=rs0(k)-qvz(k)
    1767. 

  1767.             term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- &
    1768. 

  1768.                   xka(k)*temcc(k) )
    1769. 

  1769.             term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag) &
    1770. 

  1770.                   *olambdag(k)**5.5/visc(k)
    1771. 

  1771. 

  1772.             constg2=vf1s*olambdag(k)*olambdag(k)+ &
    1773. 

  1773.                     vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
    1774. 

  1774.             tmp2=xnog*constg2
    1775. 

  1775.             tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( pgacw(k)+pgacr(k) )
    1776. 

  1776.             tmp4=amin1(0.0,tmp3)
    1777. 

  1777.             pgmlt(k)=amax1( tmp4,-qgzodt(k) )
    1778. 

  1778. 

  1779. 

  1780. !c
    1781. 

  1781. !c (4) EVAPORATION OF MELTING GRAUPEL (Pgmltevp) : HR (A19)
    1782. 

  1782. !c     but use Lin et al. coefficience
    1783. 

  1783. !c     Pgmltevp is a negative value
    1784. 

  1784. !c     abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
    1785. 

  1785. !c
    1786. 

  1786.             tmpa=rvapor*xka(k)*tem(k)*tem(k)
    1787. 

  1787.             tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
    1788. 

  1788.             tmpc=tmpa*qswz(k)*diffwv(k)
    1789. 

  1789.             tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb )
    1790. 

  1790. 

  1791. !c
    1792. 

  1792. !c     abg=2*pi*(Si-1)/rho/(A"+B")
    1793. 

  1793. !c
    1794. 

  1794.             abg=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
    1795. 

  1795. !
    1796. 

  1796. !**** allow evaporation to occur when RH less than 90%
    1797. 

  1797. !**** here not using 100% because the evaporation cooling
    1798. 

  1798. !**** of temperature is not taking into account yet; hence,
    1799. 

  1799. !**** the qgw value is a little bit larger. This will avoid
    1800. 

  1800. !**** evaporation can generate cloud.
    1801. 

  1801. !
    1802. 

  1802. !c    vf1s,vf2s=ventilation factors for snow
    1803. 

  1803. !c    vf1s=0.78,vf2s=0.31 in LIN
    1804. 

  1804. !c    constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
    1805. 

  1805. !c    constg2=vf1s*olambdag(k)*olambdag(k)+
    1806. 

  1806. !c            vf2s*schmidt(k)**0.33334*gam2pt75*constg
    1807. 

  1807. !
    1808. 

  1808.             tmp2=abg*xnog*constg2
    1809. 

  1809.             tmp3=amin1(0.0,tmp2)
    1810. 

  1810.             tmp3=amax1( tmp3,tmpd )
    1811. 

  1811.             pgmltevp(k)=amax1( tmp3,-qgzodt(k) )
    1812. 

  1812. 

  1813. !c
    1814. 

  1814. !c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29)
    1815. 

  1815. !c     Compute processes (3) only when qsz > 0.0 and qgz > 0.0
    1816. 

  1816. !c     egs=1.0
    1817. 

  1817. !c
    1818. 

  1818.            egs=1.
    1819. 

  1819.            tmpa=olambdas(k)*olambdas(k)
    1820. 

  1820.            tmpb=olambdag(k)*olambdag(k)
    1821. 

  1821.            tmpc=olambdas(k)*olambdag(k)
    1822. 

  1822.            tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k)
    1823. 

  1823.            tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
    1824. 

  1824.            tmp3=tmp1*egs*rhosnow*tmp2
    1825. 

  1825.            Pgacs(k)=amin1( tmp3,qszodt(k) )
    1826. 

  1826. 

  1827.         endif
    1828. 

  1828. 

  1829. 

  1830. !
    1831. 

  1831.   900   continue
    1832. 

  1832. 

  1833. !cc
    1834. 

  1834. !c
    1835. 

  1835. !c**********************************************************************
    1836. 

  1836. !c*****     combine all processes together and avoid negative      *****
    1837. 

  1837. !c*****     water substances
    1838. 

  1838. !***********************************************************************
    1839. 

  1839. !c
    1840. 

  1840.       if ( temcc(k) .lt. 0.0) then
    1841. 

  1841. !,delta4,1.-delta4
    1842. 

  1842. !c
    1843. 

  1843. !c  gdelta4=gindex*delta4
    1844. 

  1844. !c  g1sdelt4=gindex*(1.-delta4)
    1845. 

  1845. !c
    1846. 

  1846.            gdelta4=gindex*delta4
    1847. 

  1847.            g1sdelt4=gindex*(1.-delta4)
    1848. 

  1848. !c
    1849. 

  1849. !c  combined water vapor depletions
    1850. 

  1850. !c
    1851. 

  1851. !cc graupel
    1852. 

  1852.            tmp=psdep(k)+pgdep(k)*gindex
    1853. 

  1853.            if ( tmp .gt. qvzodt(k) ) then
    1854. 

  1854.               factor=qvzodt(k)/tmp
    1855. 

  1855.               psdep(k)=psdep(k)*factor
    1856. 

  1856.               pgdep(k)=pgdep(k)*factor*gindex
    1857. 

  1857.            end if
    1858. 

  1858. !c
    1859. 

  1859. !c  combined cloud water depletions
    1860. 

  1860. !c
    1861. 

  1861.            tmp=praut(k)+psacw(k)+psfw(k)+pracw(k)+gindex*Pgacw(k)
    1862. 

  1862.            if ( tmp .gt. qlzodt(k) ) then
    1863. 

  1863.               factor=qlzodt(k)/tmp
    1864. 

  1864.               praut(k)=praut(k)*factor
    1865. 

  1865.               psacw(k)=psacw(k)*factor
    1866. 

  1866.               psfw(k)=psfw(k)*factor
    1867. 

  1867.               pracw(k)=pracw(k)*factor
    1868. 

  1868. !cc graupel
    1869. 

  1869.               Pgacw(k)=Pgacw(k)*factor*gindex
    1870. 

  1870.            end if
    1871. 

  1871. !c
    1872. 

  1872. !c  combined cloud ice depletions
    1873. 

  1873. !c
    1874. 

  1874.            tmp=psaut(k)+psaci(k)+praci(k)+psfi(k)+Pgaci(k)*gdelta4 &
    1875. 

  1875.                +Pgacip(k)*g1sdelt4
    1876. 

  1876.            if (tmp .gt. qizodt(k) ) then
    1877. 

  1877.               factor=qizodt(k)/tmp
    1878. 

  1878.               psaut(k)=psaut(k)*factor
    1879. 

  1879.               psaci(k)=psaci(k)*factor
    1880. 

  1880.               praci(k)=praci(k)*factor
    1881. 

  1881.               psfi(k)=psfi(k)*factor
    1882. 

  1882. !cc graupel
    1883. 

  1883.               Pgaci(k)=Pgaci(k)*factor*gdelta4
    1884. 

  1884.               Pgacip(k)=Pgacip(k)*factor*g1sdelt4
    1885. 

  1885.            endif
    1886. 

  1886. !c
    1887. 

  1887. !c  combined all rain processes
    1888. 

  1888. !c
    1889. 

  1889.           tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
    1890. 

  1890.                 +Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
    1891. 

  1891.                 +Pgacrp(k)*g1sdelt4
    1892. 

  1892.           if (tmp_r .gt. qrzodt(k) ) then
    1893. 

  1893.              factor=qrzodt(k)/tmp_r
    1894. 

  1894.              piacr(k)=piacr(k)*factor
    1895. 

  1895.              psacr(k)=psacr(k)*factor
    1896. 

  1896.              prevp(k)=prevp(k)*factor
    1897. 

  1897. !cc graupel
    1898. 

  1898.              Pgfr(k)=Pgfr(k)*factor*gindex
    1899. 

  1899.              Pgacr(k)=Pgacr(k)*factor*gdelta4
    1900. 

  1900.              Pgacrp(k)=Pgacrp(k)*factor*g1sdelt4
    1901. 

  1901.           endif
    1902. 

  1902. 

  1903. !c
    1904. 

  1904. !c     if qrz < 1.0E-4 and qsz < 1.0E-4  then delta2=1.
    1905. 

  1905. !c                  (all Pracs and Psacr become to snow)
    1906. 

  1906. !c     if qrz >= 1.0E-4 or qsz >= 1.0E-4 then delta2=0.
    1907. 

  1907. !c                 (all Pracs and Psacr become to graupel)
    1908. 

  1908. !c
    1909. 

  1909.           if (qrz(k) .lt. 1.0E-4 .and. qsz(k) .lt. 1.0E-4) then
    1910. 

  1910.              delta2=1.0
    1911. 

  1911.           else
    1912. 

  1912.              delta2=0.0
    1913. 

  1913.           endif
    1914. 

  1914. !
    1915. 

  1915. !cc graupel
    1916. 

  1916. 

  1917. !c
    1918. 

  1918. !c  if qrz(k) < 1.0e-4 then delta3=1. means praci(k) -->  qs
    1919. 

  1919. !c                                          piacr(k) -->  qs
    1920. 

  1920. !c  if qrz(k) > 1.0e-4 then delta3=0. means praci(k) -->  qg
    1921. 

  1921. !c                                          piacr(k) -->  qg : Lin (20)
    1922. 

  1922. 

  1923.           if (qrz(k) .lt. 1.0e-4) then
    1924. 

  1924.              delta3=1.0
    1925. 

  1925.           else
    1926. 

  1926.              delta3=0.0
    1927. 

  1927.           endif
    1928. 

  1928. !
    1929. 

  1929. !c
    1930. 

  1930. !c     if gindex = 0.(no graupel) then delta2=1.0
    1931. 

  1931. !c                                     delta3=1.0
    1932. 

  1932. !c
    1933. 

  1933.           if (gindex .eq. 0.) then
    1934. 

  1934.               delta2=1.0
    1935. 

  1935.               delta3=1.0
    1936. 

  1936.          endif
    1937. 

  1937. !
    1938. 

  1938. !c
    1939. 

  1939. !c   combined all snow processes
    1940. 

  1940. !c
    1941. 

  1941.           tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
    1942. 

  1942.                  psfi(k)+praci(k)*delta3+piacr(k)*delta3+ &
    1943. 

  1943.                  psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ &
    1944. 

  1944.                  Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- &
    1945. 

  1945.                  Psacr(k)*delta2
    1946. 

  1946.           if ( tmp_s .gt. qszodt(k) ) then
    1947. 

  1947.              factor=qszodt(k)/tmp_s
    1948. 

  1948.              pssub(k)=pssub(k)*factor
    1949. 

  1949.              Pracs(k)=Pracs(k)*factor
    1950. 

  1950. !cc graupel
    1951. 

  1951.              Pgaut(k)=Pgaut(k)*factor*gindex
    1952. 

  1952.              Pgacs(k)=Pgacs(k)*factor*gdelta4
    1953. 

  1953.              Pgacsp(k)=Pgacsp(k)*factor*g1sdelt4
    1954. 

  1954.           endif
    1955. 

  1955. 

  1956. !cc graupel
    1957. 

  1957. !
    1958. 

  1958. 

  1959. !          if (gindex .eq. 0.) goto 998
    1960. 

  1960. !c
    1961. 

  1961. !c  combined all graupel processes
    1962. 

  1962. !c
    1963. 

  1963.            if(delta4.lt.0.5) then
    1964. 

  1964.              !Re-define pgwet to account for limiting of pgacrp,
    1965. 

  1965.              !         pgacip, pgacw and pgacsp above
    1966. 

  1966.              pgwet(k) = pgacrp(k) + pgacw(k) + pgacip(k) + pgacsp(k)
    1967. 

  1967.            end if
    1968. 

  1968.            tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4  &
    1969. 

  1969.                  -Pgacr(k)*delta4-Pgacs(k)*delta4  &
    1970. 

  1970.                  -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k)  &
    1971. 

  1971.                  -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2)  &
    1972. 

  1972.                  -praci(k)*(1-delta3)-piacr(k)*(1-delta3)
    1973. 

  1973.            if (tmp_g .gt. qgzodt(k)) then
    1974. 

  1974.                factor=qgzodt(k)/tmp_g
    1975. 

  1975.                pgsub(k)=pgsub(k)*factor
    1976. 

  1976.            endif
    1977. 

  1977. 

  1978.   998      continue
    1979. 

  1979. !c
    1980. 

  1980. !c  calculate new water substances, thetae, tem, and qvsbar
    1981. 

  1981. !c
    1982. 

  1982. 

  1983. !cc graupel
    1984. 

  1984.          pvapor(k)=-pssub(k)-psdep(k)-prevp(k)-pgsub(k)*gindex &
    1985. 

  1985.                    -pgdep(k)*gindex
    1986. 

  1986.          qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k) )
    1987. 

  1987.          pclw(k)=-praut(k)-pracw(k)-psacw(k)-psfw(k)-pgacw(k)*gindex
    1988. 

  1988. 	 if(flag_qndrop)then
    1989. 

  1989.            if( qlz(k) > 1e-20 ) &
    1990. 

  1990.               qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) )  !sg
    1991. 

  1991. 	 endif
    1992. 

  1992.          qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
    1993. 

  1993.          pcli(k)=-psaut(k)-psfi(k)-psaci(k)-praci(k)-pgaci(k)*gdelta4 &
    1994. 

  1994.                  -Pgacip(k)*g1sdelt4
    1995. 

  1995.          qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
    1996. 

  1996.          tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
    1997. 

  1997.                 +Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
    1998. 

  1998.                 +Pgacrp(k)*g1sdelt4
    1999. 

  1999.  232     format(i2,1x,6(f9.3,1x))
    2000. 

  2000.          prain(k)=-tmp_r
    2001. 

  2001.          qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
    2002. 

  2002.          tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+  &
    2003. 

  2003.                 psfi(k)+praci(k)*delta3+piacr(k)*delta3+  &
    2004. 

  2004.                 psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+  &
    2005. 

  2005.                 Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)-  &
    2006. 

  2006.                 Psacr(k)*delta2
    2007. 

  2007.          psnow(k)=-tmp_s
    2008. 

  2008.          qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
    2009. 

  2009.          qschg(k)=qschg(k)+psnow(k)
    2010. 

  2010.          qschg(k)=psnow(k)
    2011. 

  2011. !cc graupel
    2012. 

  2012.          tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 &
    2013. 

  2013.                -Pgacr(k)*delta4-Pgacs(k)*delta4 &
    2014. 

  2014.                -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
    2015. 

  2015.                -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) &
    2016. 

  2016.                -praci(k)*(1-delta3)-piacr(k)*(1-delta3)
    2017. 

  2017.  252     format(i2,1x,6(f12.9,1x))
    2018. 

  2018.  262     format(i2,1x,7(f12.9,1x))
    2019. 

  2019.          pgraupel(k)=-tmp_g
    2020. 

  2020.          pgraupel(k)=pgraupel(k)*gindex
    2021. 

  2021.          qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
    2022. 

  2022. !        qgchg(k)=qgchg(k)+pgraupel(k)
    2023. 

  2023.          qgchg(k)=pgraupel(k)
    2024. 

  2024.          qgz(k)=qgz(k)*gindex
    2025. 

  2025. 

  2026.          tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
    2027. 

  2027.          theiz(k)=theiz(k)+dtb*tmp
    2028. 

  2028.          thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
    2029. 

  2029.          tem(k)=thz(k)*tothz(k)
    2030. 

  2030. 

  2031.          temcc(k)=tem(k)-273.15
    2032. 

  2032. 

  2033. !bloss: Moves computation of saturation mixing ratio below
    2034. 

  2034. 

  2035. !         if( temcc(k) .lt. -40.0 ) qswz(k)=qsiz(k)
    2036. 

  2036. !         qlpqi=qlz(k)+qiz(k)
    2037. 

  2037. !         if ( qlpqi .eq. 0.0 ) then
    2038. 

  2038. !            qvsbar(k)=qsiz(k)
    2039. 

  2039. !         else
    2040. 

  2040. !            qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
    2041. 

  2041. !         endif
    2042. 

  2042. 

  2043. !
    2044. 

  2044.       else
    2045. 

  2045. !c
    2046. 

  2046. !c  combined cloud water depletions
    2047. 

  2047. !c
    2048. 

  2048.           tmp=praut(k)+psacw(k)+pracw(k)+pgacw(k)*gindex
    2049. 

  2049.           if ( tmp .gt. qlzodt(k) ) then
    2050. 

  2050.              factor=qlzodt(k)/tmp
    2051. 

  2051.              praut(k)=praut(k)*factor
    2052. 

  2052.              psacw(k)=psacw(k)*factor
    2053. 

  2053.              pracw(k)=pracw(k)*factor
    2054. 

  2054. !cc graupel
    2055. 

  2055.              pgacw(k)=pgacw(k)*factor*gindex
    2056. 

  2056.           end if
    2057. 

  2057. !c
    2058. 

  2058. !c  combined all snow processes
    2059. 

  2059. !c
    2060. 

  2060.           tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
    2061. 

  2061.           if (tmp_s .gt. qszodt(k) ) then
    2062. 

  2062.              factor=qszodt(k)/tmp_s
    2063. 

  2063.              psmlt(k)=psmlt(k)*factor
    2064. 

  2064.              psmltevp(k)=psmltevp(k)*factor
    2065. 

  2065. !cc graupel
    2066. 

  2066.              Pgacs(k)=Pgacs(k)*factor*gindex
    2067. 

  2067.           endif
    2068. 

  2068. 

  2069. !c
    2070. 

  2070. !c
    2071. 

  2071. !cc graupel
    2072. 

  2072. !c
    2073. 

  2073. !         if (gindex .eq. 0.) goto 997
    2074. 

  2074. 

  2075. !c
    2076. 

  2076. !c  combined all graupel processes
    2077. 

  2077. !c
    2078. 

  2078.             tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
    2079. 

  2079.             if (tmp_g .gt. qgzodt(k)) then
    2080. 

  2080.               factor=qgzodt(k)/tmp_g
    2081. 

  2081.               pgmltevp(k)=pgmltevp(k)*factor
    2082. 

  2082.               pgmlt(k)=pgmlt(k)*factor
    2083. 

  2083.             endif
    2084. 

  2084. !c
    2085. 

  2085.   997     continue
    2086. 

  2086. 

  2087. !c
    2088. 

  2088. !c  combined all rain processes
    2089. 

  2089. !c
    2090. 

  2090.           tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
    2091. 

  2091.                 +pgmlt(k)*gindex-pgacw(k)*gindex
    2092. 

  2092.           if (tmp_r .gt. qrzodt(k) ) then
    2093. 

  2093.              factor=qrzodt(k)/tmp_r
    2094. 

  2094.              prevp(k)=prevp(k)*factor
    2095. 

  2095.           endif
    2096. 

  2096. !c
    2097. 

  2097. !c
    2098. 

  2098. !c  calculate new water substances and thetae
    2099. 

  2099. !c
    2100. 

  2100. 

  2101. 

  2102.           pvapor(k)=-psmltevp(k)-prevp(k)-pgmltevp(k)*gindex
    2103. 

  2103.           qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k))
    2104. 

  2104.           pclw(k)=-praut(k)-pracw(k)-psacw(k)-pgacw(k)*gindex
    2105. 

  2105.           if(flag_qndrop)then
    2106. 

  2106. 	     if( qlz(k) > 1e-20 ) &
    2107. 

  2107.                qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) )  !sg
    2108. 

  2108. 	  endif
    2109. 

  2109.           qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
    2110. 

  2110.           pcli(k)=0.0
    2111. 

  2111.           qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
    2112. 

  2112.           tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
    2113. 

  2113.                 +pgmlt(k)*gindex-pgacw(k)*gindex
    2114. 

  2114.  242      format(i2,1x,7(f9.6,1x))
    2115. 

  2115.           prain(k)=-tmp_r
    2116. 

  2116.           tmpqrz=qrz(k)
    2117. 

  2117.           qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
    2118. 

  2118.           tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
    2119. 

  2119.           psnow(k)=-tmp_s
    2120. 

  2120.           qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
    2121. 

  2121. !         qschg(k)=qschg(k)+psnow(k)
    2122. 

  2122.           qschg(k)=psnow(k)
    2123. 

  2123. !cc graupel
    2124. 

  2124. 

  2125.           tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
    2126. 

  2126. !         write(*,272)k,pgmlt(k),pgacs(k),pgmltevp(k),
    2127. 

  2127.  272      format(i2,1x,3(f12.9,1x))
    2128. 

  2128.           pgraupel(k)=-tmp_g*gindex
    2129. 

  2129.           qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
    2130. 

  2130. !         qgchg(k)=qgchg(k)+pgraupel(k)
    2131. 

  2131.           qgchg(k)=pgraupel(k)
    2132. 

  2132.           qgz(k)=qgz(k)*gindex
    2133. 

  2133. !
    2134. 

  2134.           tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
    2135. 

  2135.           theiz(k)=theiz(k)+dtb*tmp
    2136. 

  2136.           thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
    2137. 

  2137. 

  2138.           tem(k)=thz(k)*tothz(k)
    2139. 

  2139.           temcc(k)=tem(k)-273.15
    2140. 

  2140. !         qswz(k)=episp0k*oprez(k)*  &
    2141. 

  2141. !                exp( svp2*temcc(k)/(tem(k)-svp3) )
    2142. 

  2142. 

  2143. !bloss: Moved computation of saturation mixing ratio below
    2144. 

  2144. 

  2145. !          es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
    2146. 

  2146. !          qswz(k)=ep2*es/(prez(k)-es)
    2147. 

  2147. !          qsiz(k)=qswz(k)
    2148. 

  2148. !          qvsbar(k)=qswz(k)
    2149. 

  2149. !
    2150. 

  2150.       end if
    2151. 

  2151.       preclw(k)=pclw(k)        !sg
    2152. 

  2152. 

  2153. !
    2154. 

  2154. !***********************************************************************
    2155. 

  2155. !**********              saturation adjustment                **********
    2156. 

  2156. !***********************************************************************
    2157. 

  2157. !bloss: BEGIN
    2158. 

  2158. !
    2159. 

  2159. !     compute saturation specific humidity at the temperature
    2160. 

  2160. !     which would be experienced if all cloud liquid and ice
    2161. 

  2161. !     were to become vapor.  This will make for a consistent
    2162. 

  2162. !     check of saturation.  Previously, qvsbar was computed
    2163. 

  2163. !     without accounting for the change in temperature due to
    2164. 

  2164. !     evaporation/sublimation, and as a result would
    2165. 

  2165. !     incorrectly identify some points as subsaturated.
    2166. 

  2166. 

  2167.       tmp_tem = tem(k)          ! updated temperature from above.
    2168. 

  2168. 

  2169.       if(qlz(k)+qiz(k).gt. 0.) then
    2170. 

  2170.          ! account for latent heat if all qlz and qiz were converted to vapor
    2171. 

  2171.          tmp_tem = tem(k) - xlvocp*qlz(k) - (xlvocp+xlfocp)*qiz(k)
    2172. 

  2172.       end if
    2173. 

  2173. 

  2174.      ! compute temperature in celsius
    2175. 

  2175.       tmp_temcc = tmp_tem - 273.15
    2176. 

  2176. 

  2177.       ! estimate lower bound on saturation vapor pressure 
    2178. 

  2178.       !   (ice if <0C, liquid aboce)
    2179. 

  2179.       if (tmp_temcc .lt. 0.) then
    2180. 

  2180.          es=1000.*svp1*exp( 21.8745584*(tmp_tem-273.16)/(tmp_tem-7.66) ) !ice
    2181. 

  2181.       else
    2182. 

  2182.          es=1000.*svp1*exp( svp2*tmp_temcc/(tmp_tem-svp3) ) !liquid
    2183. 

  2183.       end if
    2184. 

  2184. 

  2185.       ! compute lower bound on saturation mixing ratio.
    2186. 

  2186.       qvsbar(k)=ep2*es/(prez(k)-es)
    2187. 

  2187. !
    2188. 

  2188. !bloss: END
    2189. 

  2189. !
    2190. 

  2190. !    allow supersaturation exits linearly from 0% at 500 mb to 50%
    2191. 

  2191. !    above 300 mb
    2192. 

  2192. !    5.0e-5=1.0/(500mb-300mb)
    2193. 

  2193. !
    2194. 

  2194.          rsat=1.0+0.5*(50000.0-prez(k))*5.0e-5
    2195. 

  2195.          rsat=amax1(1.0,rsat)
    2196. 

  2196.          rsat=amin1(1.5,rsat)
    2197. 

  2197.          rsat=1.0
    2198. 

  2198.          if( qvz(k)+qlz(k)+qiz(k) .lt. rsat*qvsbar(k) ) then
    2199. 

  2199. 

  2200. !c
    2201. 

  2201. !c   unsaturated
    2202. 

  2202. !c
    2203. 

  2203.           qvz(k)=qvz(k)+qlz(k)+qiz(k)
    2204. 

  2204.           qlz(k)=0.0
    2205. 

  2205.           qiz(k)=0.0
    2206. 

  2206. 

  2207.           thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
    2208. 

  2208.           tem(k)=thz(k)*tothz(k)
    2209. 

  2209.           temcc(k)=tem(k)-273.15
    2210. 

  2210. 

  2211.           go to 1800
    2212. 

  2212. !
    2213. 

  2213.         else
    2214. 

  2214. !c
    2215. 

  2215. !c   saturated
    2216. 

  2216. !c
    2217. 

  2217. !
    2218. 

  2218.           pladj(k)=qlz(k)
    2219. 

  2219.           piadj(k)=qiz(k)
    2220. 

  2220. !
    2221. 

  2221. 

  2222.           CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
    2223. 

  2223.                       k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0 )
    2224. 

  2224. 

  2225. !
    2226. 

  2226.           pladj(k)=odtb*(qlz(k)-pladj(k))
    2227. 

  2227.           piadj(k)=odtb*(qiz(k)-piadj(k))
    2228. 

  2228. !
    2229. 

  2229.           pclw(k)=pclw(k)+pladj(k)
    2230. 

  2230.           pcli(k)=pcli(k)+piadj(k)
    2231. 

  2231.           pvapor(k)=pvapor(k)-( pladj(k)+piadj(k) )
    2232. 

  2232. !
    2233. 

  2233.           thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
    2234. 

  2234.           tem(k)=thz(k)*tothz(k)
    2235. 

  2235. 

  2236.           temcc(k)=tem(k)-273.15
    2237. 

  2237. 

  2238. !         qswz(k)=episp0k*oprez(k)*  &
    2239. 

  2239. !                 exp( svp2*temcc(k)/(tem(k)-svp3) )
    2240. 

  2240.           es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
    2241. 

  2241.           qswz(k)=ep2*es/(prez(k)-es)
    2242. 

  2242.           if (tem(k) .lt. 273.15 ) then
    2243. 

  2243. !            qsiz(k)=episp0k*oprez(k)* &
    2244. 

  2244. !                   exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
    2245. 

  2245.              es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
    2246. 

  2246.              qsiz(k)=ep2*es/(prez(k)-es)
    2247. 

  2247.              if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
    2248. 

  2248.           else
    2249. 

  2249.              qsiz(k)=qswz(k)
    2250. 

  2250.           endif
    2251. 

  2251.           qlpqi=qlz(k)+qiz(k)
    2252. 

  2252.           if ( qlpqi .eq. 0.0 ) then
    2253. 

  2253.              qvsbar(k)=qsiz(k)
    2254. 

  2254.           else
    2255. 

  2255.              qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
    2256. 

  2256.           endif
    2257. 

  2257. 

  2258.         end if
    2259. 

  2259. 

  2260. !
    2261. 

  2261. !***********************************************************************
    2262. 

  2262. !*****     melting and freezing of cloud ice and cloud water       *****
    2263. 

  2263. !***********************************************************************
    2264. 

  2264.         qlpqi=qlz(k)+qiz(k)
    2265. 

  2265.         if(qlpqi .le. 0.0) go to 1800
    2266. 

  2266. !
    2267. 

  2267. !c
    2268. 

  2268. !c (1)  HOMOGENEOUS NUCLEATION WHEN T< -40 C (Pihom)
    2269. 

  2269. !c
    2270. 

  2270.         if(temcc(k) .lt. -40.0) pihom(k)=qlz(k)*odtb
    2271. 

  2271. !c
    2272. 

  2272. !c (2)  MELTING OF ICE CRYSTAL WHEN T> 0 C (Pimlt)
    2273. 

  2273. !c
    2274. 

  2274.         if(temcc(k) .gt. 0.0) pimlt(k)=qiz(k)*odtb
    2275. 

  2275. !c
    2276. 

  2276. !c (3) PRODUCTION OF CLOUD ICE BY BERGERON PROCESS (Pidw): Hsie (p957)
    2277. 

  2277. !c     this process only considered when -31 C < T < 0 C
    2278. 

  2278. !c
    2279. 

  2279.         if(temcc(k) .lt. 0.0 .and. temcc(k) .gt. -31.0) then
    2280. 

  2280. !c!
    2281. 

  2281. !c!   parama1 and parama2 functions must be user supplied
    2282. 

  2282. !c!
    2283. 

  2283.           a1=parama1( temcc(k) )
    2284. 

  2284.           a2=parama2( temcc(k) )
    2285. 

  2285. !! change unit from cgs to mks
    2286. 

  2286.           a1=a1*0.001**(1.0-a2)
    2287. 

  2287.           xnin=xni0*exp(-bni*temcc(k))
    2288. 

  2288.           pidw(k)=xnin*orho(k)*(a1*xmnin**a2)
    2289. 

  2289.         end if
    2290. 

  2290. !
    2291. 

  2291.         pcli(k)=pcli(k)+pihom(k)-pimlt(k)+pidw(k)
    2292. 

  2292.         pclw(k)=pclw(k)-pihom(k)+pimlt(k)-pidw(k)
    2293. 

  2293.         qlz(k)=amax1( 0.0,qlz(k)+dtb*(-pihom(k)+pimlt(k)-pidw(k)) )
    2294. 

  2294.         qiz(k)=amax1( 0.0,qiz(k)+dtb*(pihom(k)-pimlt(k)+pidw(k)) )
    2295. 

  2295. 

  2296. !
    2297. 

  2297.         CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
    2298. 

  2298.                     k, xLvocp, xLfocp, episp0k ,EP2,SVP1,SVP2,SVP3,SVPT0)
    2299. 

  2299. 

  2300.         thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
    2301. 

  2301.         tem(k)=thz(k)*tothz(k)
    2302. 

  2302. 

  2303.         temcc(k)=tem(k)-273.15
    2304. 

  2304. 

  2305. !       qswz(k)=episp0k*oprez(k)* &
    2306. 

  2306. !              exp( svp2*temcc(k)/(tem(k)-svp3) )
    2307. 

  2307.         es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
    2308. 

  2308.         qswz(k)=ep2*es/(prez(k)-es)
    2309. 

  2309. 

  2310.         if (tem(k) .lt. 273.15 ) then
    2311. 

  2311. !          qsiz(k)=episp0k*oprez(k)* &
    2312. 

  2312. !                 exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
    2313. 

  2313.            es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
    2314. 

  2314.            qsiz(k)=ep2*es/(prez(k)-es)
    2315. 

  2315.            if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
    2316. 

  2316.         else
    2317. 

  2317.            qsiz(k)=qswz(k)
    2318. 

  2318.         endif
    2319. 

  2319.         qlpqi=qlz(k)+qiz(k)
    2320. 

  2320.         if ( qlpqi .eq. 0.0 ) then
    2321. 

  2321.            qvsbar(k)=qsiz(k)
    2322. 

  2322.         else
    2323. 

  2323.            qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
    2324. 

  2324.         endif
    2325. 

  2325. 

  2326. 1800  continue
    2327. 

  2327. !
    2328. 

  2328. !***********************************************************************
    2329. 

  2329. !**********    integrate the productions of rain and snow     **********
    2330. 

  2330. !***********************************************************************
    2331. 

  2331. !c
    2332. 

  2332. 

  2333. 2000  continue
    2334. 

  2334. 

  2335. 

  2336. !---------------------------------------------------------------------
    2337. 

  2337. 

  2338. !
    2339. 

  2339. !***********************************************************************
    2340. 

  2340. !******      Write terms in cloud physics to time series dataset   *****
    2341. 

  2341. !***********************************************************************
    2342. 

  2342. !
    2343. 

  2343. !     open(unit=24,form='formatted',status='new',
    2344. 

  2344. !    &     file='cloud.dat')
    2345. 

  2345. 

  2346. !9030  format(10e12.6)
    2347. 

  2347. 

  2348. !      write(24,*)'tmp'
    2349. 

  2349. !      write(24,9030) (tem(k),k=kts+1,kte)
    2350. 

  2350. !      write(24,*)'qiz'
    2351. 

  2351. !      write(24,9030) (qiz(k),k=kts+1,kte)
    2352. 

  2352. !      write(24,*)'qsz'
    2353. 

  2353. !      write(24,9030) (qsz(k),k=kts+1,kte)
    2354. 

  2354. !      write(24,*)'qrz'
    2355. 

  2355. !      write(24,9030) (qrz(k),k=kts+1,kte)
    2356. 

  2356. !      write(24,*)'qgz'
    2357. 

  2357. !      write(24,9030) (qgz(k),k=kts+1,kte)
    2358. 

  2358. !      write(24,*)'qvoqsw'
    2359. 

  2359. !      write(24,9030) (qvoqswz(k),k=kts+1,kte)
    2360. 

  2360. !      write(24,*)'qvoqsi'
    2361. 

  2361. !      write(24,9030) (qvoqsiz(k),k=kts+1,kte)
    2362. 

  2362. !      write(24,*)'vtr'
    2363. 

  2363. !      write(24,9030) (vtr(k),k=kts+1,kte)
    2364. 

  2364. !      write(24,*)'vts'
    2365. 

  2365. !      write(24,9030) (vts(k),k=kts+1,kte)
    2366. 

  2366. !      write(24,*)'vtg'
    2367. 

  2367. !      write(24,9030) (vtg(k),k=kts+1,kte)
    2368. 

  2368. !      write(24,*)'pclw'
    2369. 

  2369. !      write(24,9030) (pclw(k),k=kts+1,kte)
    2370. 

  2370. !      write(24,*)'pvapor'
    2371. 

  2371. !      write(24,9030) (pvapor(k),k=kts+1,kte)
    2372. 

  2372. !      write(24,*)'pcli'
    2373. 

  2373. !      write(24,9030) (pcli(k),k=kts+1,kte)
    2374. 

  2374. !      write(24,*)'pimlt'
    2375. 

  2375. !      write(24,9030) (pimlt(k),k=kts+1,kte)
    2376. 

  2376. !      write(24,*)'pihom'
    2377. 

  2377. !      write(24,9030) (pihom(k),k=kts+1,kte)
    2378. 

  2378. !      write(24,*)'pidw'
    2379. 

  2379. !      write(24,9030) (pidw(k),k=kts+1,kte)
    2380. 

  2380. !      write(24,*)'prain'
    2381. 

  2381. !      write(24,9030) (prain(k),k=kts+1,kte)
    2382. 

  2382. !      write(24,*)'praut'
    2383. 

  2383. !      write(24,9030) (praut(k),k=kts+1,kte)
    2384. 

  2384. !      write(24,*)'pracw'
    2385. 

  2385. !      write(24,9030) (pracw(k),k=kts+1,kte)
    2386. 

  2386. !      write(24,*)'prevp'
    2387. 

  2387. !      write(24,9030) (prevp(k),k=kts+1,kte)
    2388. 

  2388. !      write(24,*)'psnow'
    2389. 

  2389. !      write(24,9030) (psnow(k),k=kts+1,kte)
    2390. 

  2390. !      write(24,*)'psaut'
    2391. 

  2391. !      write(24,9030) (psaut(k),k=kts+1,kte)
    2392. 

  2392. !      write(24,*)'psfw'
    2393. 

  2393. !      write(24,9030) (psfw(k),k=kts+1,kte)
    2394. 

  2394. !      write(24,*)'psfi'
    2395. 

  2395. !      write(24,9030) (psfi(k),k=kts+1,kte)
    2396. 

  2396. !      write(24,*)'praci'
    2397. 

  2397. !      write(24,9030) (praci(k),k=kts+1,kte)
    2398. 

  2398. !      write(24,*)'piacr'
    2399. 

  2399. !      write(24,9030) (piacr(k),k=kts+1,kte)
    2400. 

  2400. !      write(24,*)'psaci'
    2401. 

  2401. !      write(24,9030) (psaci(k),k=kts+1,kte)
    2402. 

  2402. !      write(24,*)'psacw'
    2403. 

  2403. !      write(24,9030) (psacw(k),k=kts+1,kte)
    2404. 

  2404. !      write(24,*)'psdep'
    2405. 

  2405. !      write(24,9030) (psdep(k),k=kts+1,kte)
    2406. 

  2406. !      write(24,*)'pssub'
    2407. 

  2407. !      write(24,9030) (pssub(k),k=kts+1,kte)
    2408. 

  2408. !      write(24,*)'pracs'
    2409. 

  2409. !      write(24,9030) (pracs(k),k=kts+1,kte)
    2410. 

  2410. !      write(24,*)'psacr'
    2411. 

  2411. !      write(24,9030) (psacr(k),k=kts+1,kte)
    2412. 

  2412. !      write(24,*)'psmlt'
    2413. 

  2413. !      write(24,9030) (psmlt(k),k=kts+1,kte)
    2414. 

  2414. !      write(24,*)'psmltevp'
    2415. 

  2415. !      write(24,9030) (psmltevp(k),k=kts+1,kte)
    2416. 

  2416. !      write(24,*)'pladj'
    2417. 

  2417. !      write(24,9030) (pladj(k),k=kts+1,kte)
    2418. 

  2418. !      write(24,*)'piadj'
    2419. 

  2419. !      write(24,9030) (piadj(k),k=kts+1,kte)
    2420. 

  2420. !      write(24,*)'pgraupel'
    2421. 

  2421. !      write(24,9030) (pgraupel(k),k=kts+1,kte)
    2422. 

  2422. !      write(24,*)'pgaut'
    2423. 

  2423. !      write(24,9030) (pgaut(k),k=kts+1,kte)
    2424. 

  2424. !      write(24,*)'pgfr'
    2425. 

  2425. !      write(24,9030) (pgfr(k),k=kts+1,kte)
    2426. 

  2426. !      write(24,*)'pgacw'
    2427. 

  2427. !      write(24,9030) (pgacw(k),k=kts+1,kte)
    2428. 

  2428. !      write(24,*)'pgaci'
    2429. 

  2429. !      write(24,9030) (pgaci(k),k=kts+1,kte)
    2430. 

  2430. !      write(24,*)'pgacr'
    2431. 

  2431. !      write(24,9030) (pgacr(k),k=kts+1,kte)
    2432. 

  2432. !      write(24,*)'pgacs'
    2433. 

  2433. !      write(24,9030) (pgacs(k),k=kts+1,kte)
    2434. 

  2434. !      write(24,*)'pgacip'
    2435. 

  2435. !      write(24,9030) (pgacip(k),k=kts+1,kte)
    2436. 

  2436. !      write(24,*)'pgacrP'
    2437. 

  2437. !      write(24,9030) (pgacrP(k),k=kts+1,kte)
    2438. 

  2438. !      write(24,*)'pgacsp'
    2439. 

  2439. !      write(24,9030) (pgacsp(k),k=kts+1,kte)
    2440. 

  2440. !      write(24,*)'pgwet'
    2441. 

  2441. !      write(24,9030) (pgwet(k),k=kts+1,kte)
    2442. 

  2442. !      write(24,*)'pdry'
    2443. 

  2443. !      write(24,9030) (pdry(k),k=kts+1,kte)
    2444. 

  2444. !      write(24,*)'pgsub'
    2445. 

  2445. !      write(24,9030) (pgsub(k),k=kts+1,kte)
    2446. 

  2446. !      write(24,*)'pgdep'
    2447. 

  2447. !      write(24,9030) (pgdep(k),k=kts+1,kte)
    2448. 

  2448. !      write(24,*)'pgmlt'
    2449. 

  2449. !      write(24,9030) (pgmlt(k),k=kts+1,kte)
    2450. 

  2450. !      write(24,*)'pgmltevp'
    2451. 

  2451. !      write(24,9030) (pgmltevp(k),k=kts+1,kte)
    2452. 

  2452. 

  2453. 

  2454. 

  2455. !**** below if qv < qvmin then qv=qvmin, ql=0.0, and qi=0.0
    2456. 

  2456. !
    2457. 

  2457.       do k=kts+1,kte
    2458. 

  2458.          if ( qvz(k) .lt. qvmin ) then
    2459. 

  2459.             qlz(k)=0.0
    2460. 

  2460.             qiz(k)=0.0
    2461. 

  2461.             qvz(k)=amax1( qvmin,qvz(k)+qlz(k)+qiz(k) )
    2462. 

  2462.          end if
    2463. 

  2463.       enddo
    2464. 

  2464. !
    2465. 

  2465.   END SUBROUTINE clphy1d
    2466. 

  2466. 

  2467. 

  2468. !---------------------------------------------------------------------
    2469. 

  2469. !                         SATURATED ADJUSTMENT
    2470. 

  2470. !---------------------------------------------------------------------
    2471. 

  2471.       SUBROUTINE satadj(qvz, qlz, qiz, prez, theiz, thz, tothz,      &
    2472. 

  2472.                         kts, kte, k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0)
    2473. 

  2473. !---------------------------------------------------------------------
    2474. 

  2474.    IMPLICIT NONE
    2475. 

  2475. !---------------------------------------------------------------------
    2476. 

  2476. !  This program use Newton's method for finding saturated temperature
    2477. 

  2477. !  and saturation mixing ratio.
    2478. 

  2478. !
    2479. 

  2479. ! In this saturation adjustment scheme we assume
    2480. 

  2480. ! (1)  the saturation mixing ratio is the mass weighted average of
    2481. 

  2481. !      saturation values over liquid water (qsw), and ice (qsi)
    2482. 

  2482. !      following Lord et al., 1984 and Tao, 1989
    2483. 

  2483. !
    2484. 

  2484. ! (2) the percentage of cloud liquid and cloud ice will
    2485. 

  2485. !      be fixed during the saturation calculation
    2486. 

  2486. !---------------------------------------------------------------------
    2487. 

  2487. !
    2488. 

  2488. 

  2489.   INTEGER, INTENT(IN   )             :: kts, kte, k
    2490. 

  2490. 

  2491.   REAL,      DIMENSION( kts:kte ),                                   &
    2492. 

  2492.                        INTENT(INOUT) :: qvz, qlz, qiz
    2493. 

  2493. !
    2494. 

  2494.   REAL,      DIMENSION( kts:kte ),                                   &
    2495. 

  2495.                        INTENT(IN   ) :: prez, theiz, tothz
    2496. 

  2496. 

  2497.   REAL,     INTENT(IN   )            :: xLvocp, xLfocp, episp0k
    2498. 

  2498.   REAL,     INTENT(IN   )            :: EP2,SVP1,SVP2,SVP3,SVPT0
    2499. 

  2499. 

  2500. ! LOCAL VARS
    2501. 

  2501. 

  2502.   INTEGER                            :: n
    2503. 

  2503. 

  2504.   REAL, DIMENSION( kts:kte )         :: thz, tem, temcc, qsiz,       &
    2505. 

  2505.                                         qswz, qvsbar
    2506. 

  2506. 

  2507.   REAL ::   qsat, qlpqi, ratql, t0, t1, tmp1, ratqi, tsat, absft,    &
    2508. 

  2508.             denom1, denom2, dqvsbar, ftsat, dftsat, qpz,             &
    2509. 

  2509.             gindex, es
    2510. 

  2510. 

  2511.   REAL ::   tem_noliqice, qsat_noliqice !bloss
    2512. 

  2512. !
    2513. 

  2513. !---------------------------------------------------------------------
    2514. 

  2514. 

  2515.       thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
    2516. 

  2516. 

  2517.       tem(k)=tothz(k)*thz(k)
    2518. 

  2518. 

  2519. !bloss: BEGIN
    2520. 

  2520.       tem_noliqice = tem(k) - xlvocp*qlz(k) - (xLvocp+xLfocp)*qiz(k)
    2521. 

  2521. 

  2522.       if (tem_noliqice .gt. 273.15) then
    2523. 

  2523.          es=1000.*svp1*exp( svp2*(tem_noliqice-svpt0)/(tem_noliqice-svp3) )
    2524. 

  2524.          qsat_noliqice=ep2*es/(prez(k)-es)
    2525. 

  2525.       else
    2526. 

  2526.         qsat_noliqice=episp0k/prez(k)*  &
    2527. 

  2527.              exp( 21.8745584*(tem_noliqice-273.15)/(tem_noliqice-7.66) )
    2528. 

  2528.       end if
    2529. 

  2529. !bloss: END
    2530. 

  2530.       qpz=qvz(k)+qlz(k)+qiz(k)
    2531. 

  2531.       if (qpz .lt. qsat_noliqice) then
    2532. 

  2532.          qvz(k)=qpz
    2533. 

  2533.          qiz(k)=0.0
    2534. 

  2534.          qlz(k)=0.0
    2535. 

  2535.          go to 400
    2536. 

  2536.       end if
    2537. 

  2537.       qlpqi=qlz(k)+qiz(k)
    2538. 

  2538.       if( qlpqi .ge. 1.0e-5) then
    2539. 

  2539.         ratql=qlz(k)/qlpqi
    2540. 

  2540.         ratqi=qiz(k)/qlpqi
    2541. 

  2541.       else
    2542. 

  2542.         t0=273.15
    2543. 

  2543. !       t1=233.15
    2544. 

  2544.         t1=248.15
    2545. 

  2545.         tmp1=( t0-tem(k) )/(t0-t1)
    2546. 

  2546.         tmp1=amin1(1.0,tmp1)
    2547. 

  2547.         tmp1=amax1(0.0,tmp1)
    2548. 

  2548.         ratqi=tmp1
    2549. 

  2549.         ratql=1.0-tmp1
    2550. 

  2550.       end if
    2551. 

  2551. !
    2552. 

  2552. !
    2553. 

  2553. !--  saturation mixing ratios over water and ice
    2554. 

  2554. !--  at the outset we will follow Bolton 1980 MWR for
    2555. 

  2555. !--  the water and Murray JAS 1967 for the ice
    2556. 

  2556. !
    2557. 

  2557. !-- dqvsbar=d(qvsbar)/dT
    2558. 

  2558. !-- ftsat=F(Tsat)
    2559. 

  2559. !-- dftsat=d(F(T))/dT
    2560. 

  2560. !
    2561. 

  2561. !  First guess of tsat
    2562. 

  2562. 

  2563.       tsat=tem(k)
    2564. 

  2564.       absft=1.0
    2565. 

  2565. !
    2566. 

  2566.       do 200 n=1,20
    2567. 

  2567.          denom1=1.0/(tsat-svp3)
    2568. 

  2568.          denom2=1.0/(tsat-7.66)
    2569. 

  2569. !        qswz(k)=episp0k/prez(k)*  &
    2570. 

  2570. !                exp( svp2*denom1*(tsat-273.15) )
    2571. 

  2571.          es=1000.*svp1*exp( svp2*denom1*(tsat-svpt0) )
    2572. 

  2572.          qswz(k)=ep2*es/(prez(k)-es)
    2573. 

  2573.          if (tem(k) .lt. 273.15) then
    2574. 

  2574. !           qsiz(k)=episp0k/prez(k)*  &
    2575. 

  2575. !                   exp( 21.8745584*denom2*(tsat-273.15) )
    2576. 

  2576.             es=1000.*svp1*exp( 21.8745584*denom2*(tsat-273.15) )
    2577. 

  2577.             qsiz(k)=ep2*es/(prez(k)-es)
    2578. 

  2578.             if (tem(k) .lt. 233.15) qswz(k)=qsiz(k)
    2579. 

  2579.          else
    2580. 

  2580.             qsiz(k)=qswz(k)
    2581. 

  2581.          endif
    2582. 

  2582.          qvsbar(k)=ratql*qswz(k)+ratqi*qsiz(k)
    2583. 

  2583. !
    2584. 

  2584. !        if( absft .lt. 0.01 .and. n .gt. 3 ) go to 300
    2585. 

  2585.          if( absft .lt. 0.01 ) go to 300
    2586. 

  2586. !
    2587. 

  2587.          dqvsbar=ratql*qswz(k)*svp2*243.5*denom1*denom1+  &
    2588. 

  2588.                  ratqi*qsiz(k)*21.8745584*265.5*denom2*denom2
    2589. 

  2589.          ftsat=tsat+(xlvocp+ratqi*xlfocp)*qvsbar(k)-  &
    2590. 

  2590.                tothz(k)*theiz(k)-xlfocp*ratqi*(qvz(k)+qlz(k)+qiz(k))
    2591. 

  2591.          dftsat=1.0+(xlvocp+ratqi*xlfocp)*dqvsbar
    2592. 

  2592.          tsat=tsat-ftsat/dftsat
    2593. 

  2593.          absft=abs(ftsat)
    2594. 

  2594. 

  2595. 200   continue
    2596. 

  2596. 9020  format(1x,'point can not converge, absft,n=',e12.5,i5)
    2597. 

  2597. !
    2598. 

  2598. 300   continue
    2599. 

  2599.       if( qpz .gt. qvsbar(k) ) then
    2600. 

  2600.         qvz(k)=qvsbar(k)
    2601. 

  2601.         qiz(k)=ratqi*( qpz-qvz(k) )
    2602. 

  2602.         qlz(k)=ratql*( qpz-qvz(k) )
    2603. 

  2603.       else
    2604. 

  2604.         qvz(k)=qpz
    2605. 

  2605.         qiz(k)=0.0
    2606. 

  2606.         qlz(k)=0.0
    2607. 

  2607.       end if
    2608. 

  2608.  400  continue
    2609. 

  2609. 

  2610.    END SUBROUTINE satadj
    2611. 

  2611. 

  2612. 

  2613. !----------------------------------------------------------------
    2614. 

  2614.    REAL FUNCTION parama1(temp)
    2615. 

  2615. !----------------------------------------------------------------
    2616. 

  2616.    IMPLICIT NONE
    2617. 

  2617. !----------------------------------------------------------------
    2618. 

  2618. !  This program calculate the parameter for crystal growth rate
    2619. 

  2619. !  in Bergeron process
    2620. 

  2620. !----------------------------------------------------------------
    2621. 

  2621. 

  2622.       REAL, INTENT (IN   )   :: temp
    2623. 

  2623.       REAL, DIMENSION(32)    :: a1
    2624. 

  2624.       INTEGER                :: i1, i1p1
    2625. 

  2625.       REAL                   :: ratio
    2626. 

  2626. 

  2627.       data a1/0.100e-10,0.7939e-7,0.7841e-6,0.3369e-5,0.4336e-5, &
    2628. 

  2628.               0.5285e-5,0.3728e-5,0.1852e-5,0.2991e-6,0.4248e-6, &
    2629. 

  2629.               0.7434e-6,0.1812e-5,0.4394e-5,0.9145e-5,0.1725e-4, &
    2630. 

  2630.               0.3348e-4,0.1725e-4,0.9175e-5,0.4412e-5,0.2252e-5, &
    2631. 

  2631.               0.9115e-6,0.4876e-6,0.3473e-6,0.4758e-6,0.6306e-6, &
    2632. 

  2632.               0.8573e-6,0.7868e-6,0.7192e-6,0.6513e-6,0.5956e-6, &
    2633. 

  2633.               0.5333e-6,0.4834e-6/
    2634. 

  2634. 

  2635.       i1=int(-temp)+1
    2636. 

  2636.       i1p1=i1+1
    2637. 

  2637.       ratio=-(temp)-float(i1-1)
    2638. 

  2638.       parama1=a1(i1)+ratio*( a1(i1p1)-a1(i1) )
    2639. 

  2639. 

  2640.    END FUNCTION parama1
    2641. 

  2641. 

  2642. !----------------------------------------------------------------
    2643. 

  2643.    REAL FUNCTION parama2(temp)
    2644. 

  2644. !----------------------------------------------------------------
    2645. 

  2645.    IMPLICIT NONE
    2646. 

  2646. !----------------------------------------------------------------
    2647. 

  2647. !  This program calculate the parameter for crystal growth rate
    2648. 

  2648. !  in Bergeron process
    2649. 

  2649. !----------------------------------------------------------------
    2650. 

  2650. 

  2651.       REAL, INTENT (IN   )   :: temp
    2652. 

  2652.       REAL, DIMENSION(32)    :: a2
    2653. 

  2653.       INTEGER                :: i1, i1p1
    2654. 

  2654.       REAL                   :: ratio
    2655. 

  2655. 

  2656.       data a2/0.0100,0.4006,0.4831,0.5320,0.5307,0.5319,0.5249, &
    2657. 

  2657.               0.4888,0.3849,0.4047,0.4318,0.4771,0.5183,0.5463, &
    2658. 

  2658.               0.5651,0.5813,0.5655,0.5478,0.5203,0.4906,0.4447, &
    2659. 

  2659.               0.4126,0.3960,0.4149,0.4320,0.4506,0.4483,0.4460, &
    2660. 

  2660.               0.4433,0.4413,0.4382,0.4361/
    2661. 

  2661.       i1=int(-temp)+1
    2662. 

  2662.       i1p1=i1+1
    2663. 

  2663.       ratio=-(temp)-float(i1-1)
    2664. 

  2664.       parama2=a2(i1)+ratio*( a2(i1p1)-a2(i1) )
    2665. 

  2665. 

  2666.    END FUNCTION parama2
    2667. 

  2667. 

  2668. !----------------------------------------------------------------
    2669. 

  2669.    REAL FUNCTION ggamma(X)
    2670. 

  2670. !----------------------------------------------------------------
    2671. 

  2671.    IMPLICIT NONE
    2672. 

  2672. !----------------------------------------------------------------
    2673. 

  2673.       REAL, INTENT(IN   ) :: x
    2674. 

  2674.       REAL, DIMENSION(8)  :: B
    2675. 

  2675.       INTEGER             ::j, K1
    2676. 

  2676.       REAL                ::PF, G1TO2 ,TEMP
    2677. 

  2677. 

  2678.       DATA B/-.577191652,.988205891,-.897056937,.918206857,  &
    2679. 

  2679.              -.756704078,.482199394,-.193527818,.035868343/
    2680. 

  2680. 

  2681.       PF=1.
    2682. 

  2682.       TEMP=X
    2683. 

  2683.       DO 10 J=1,200
    2684. 

  2684.       IF (TEMP .LE. 2) GO TO 20
    2685. 

  2685.       TEMP=TEMP-1.
    2686. 

  2686.    10 PF=PF*TEMP
    2687. 

  2687.   100 FORMAT(//,5X,'module_mp_lin: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
    2688. 

  2688.       WRITE(wrf_err_message,100)X
    2689. 

  2689.       CALL wrf_error_fatal(wrf_err_message)
    2690. 

  2690.    20 G1TO2=1.
    2691. 

  2691.       TEMP=TEMP - 1.
    2692. 

  2692.       DO 30 K1=1,8
    2693. 

  2693.    30 G1TO2=G1TO2 + B(K1)*TEMP**K1
    2694. 

  2694.       ggamma=PF*G1TO2
    2695. 

  2695. 

  2696.       END FUNCTION ggamma
    2697. 

  2697. 

  2698. !----------------------------------------------------------------
    2699. 

  2699. 

  2700. END MODULE module_mp_lin
    2701. 

  2701. 

  2702.