PageRenderTime 63ms CodeModel.GetById 15ms RepoModel.GetById 0ms app.codeStats 0ms

/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

Large files files are truncated, but you can click here to view the full file

    

  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)…
    1494.

Large files files are truncated, but you can click here to view the full file