PageRenderTime 57ms CodeModel.GetById 20ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_mp_wsm6.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2218 lines | 1636 code | 1 blank | 581 comment | 120 complexity | c10535f792aaff054b814b969f7f2130 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. #if ( RWORDSIZE == 4 )
    2. 

  2. #  define VREC vsrec
    3. 

  3. #  define VSQRT vssqrt
    4. 

  4. #else
    5. 

  5. #  define VREC vrec
    6. 

  6. #  define VSQRT vsqrt
    7. 

  7. #endif
    8. 

  8. 

  9. MODULE module_mp_wsm6
    10. 

  10. !
    11. 

  11. !
    12. 

  12.    REAL, PARAMETER, PRIVATE :: dtcldcr     = 120. ! maximum time step for minor loops
    13. 

  13.    REAL, PARAMETER, PRIVATE :: n0r = 8.e6         ! intercept parameter rain
    14. 

  14.    REAL, PARAMETER, PRIVATE :: n0g = 4.e6         ! intercept parameter graupel
    15. 

  15.    REAL, PARAMETER, PRIVATE :: avtr = 841.9       ! a constant for terminal velocity of rain
    16. 

  16.    REAL, PARAMETER, PRIVATE :: bvtr = 0.8         ! a constant for terminal velocity of rain
    17. 

  17.    REAL, PARAMETER, PRIVATE :: r0 = .8e-5         ! 8 microm  in contrast to 10 micro m
    18. 

  18.    REAL, PARAMETER, PRIVATE :: peaut = .55        ! collection efficiency
    19. 

  19.    REAL, PARAMETER, PRIVATE :: xncr = 3.e8        ! maritime cloud in contrast to 3.e8 in tc80
    20. 

  20.    REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5    ! the dynamic viscosity kgm-1s-1
    21. 

  21.    REAL, PARAMETER, PRIVATE :: avts = 11.72       ! a constant for terminal velocity of snow
    22. 

  22.    REAL, PARAMETER, PRIVATE :: bvts = .41         ! a constant for terminal velocity of snow
    23. 

  23.    REAL, PARAMETER, PRIVATE :: avtg = 330.        ! a constant for terminal velocity of graupel
    24. 

  24.    REAL, PARAMETER, PRIVATE :: bvtg = 0.8         ! a constant for terminal velocity of graupel
    25. 

  25.    REAL, PARAMETER, PRIVATE :: deng = 500.        ! density of graupel
    26. 

  26.    REAL, PARAMETER, PRIVATE :: n0smax =  1.e11    ! maximum n0s (t=-90C unlimited)
    27. 

  27.    REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4   ! limited maximum value for slope parameter of rain
    28. 

  28.    REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5   ! limited maximum value for slope parameter of snow
    29. 

  29.    REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4   ! limited maximum value for slope parameter of graupel
    30. 

  30.    REAL, PARAMETER, PRIVATE :: dicon = 11.9       ! constant for the cloud-ice diamter
    31. 

  31.    REAL, PARAMETER, PRIVATE :: dimax = 500.e-6    ! limited maximum value for the cloud-ice diamter
    32. 

  32.    REAL, PARAMETER, PRIVATE :: n0s = 2.e6         ! temperature dependent intercept parameter snow
    33. 

  33.    REAL, PARAMETER, PRIVATE :: alpha = .12        ! .122 exponen factor for n0s
    34. 

  34.    REAL, PARAMETER, PRIVATE :: pfrz1 = 100.       ! constant in Biggs freezing
    35. 

  35.    REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66       ! constant in Biggs freezing
    36. 

  36.    REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9     ! minimun values for qr, qs, and qg
    37. 

  37.    REAL, PARAMETER, PRIVATE :: eacrc = 1.0        ! Snow/cloud-water collection efficiency
    38. 

  38.    REAL, PARAMETER, PRIVATE :: dens  =  100.0     ! Density of snow
    39. 

  39.    REAL, PARAMETER, PRIVATE :: qs0   =  6.e-4     ! threshold amount for aggretion to occur
    40. 

  40.    REAL, SAVE ::                                      &
    41. 

  41.              qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
    42. 

  42.              g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr,    &
    43. 

  43.              bvtr6,g6pbr,                             &
    44. 

  44.              precr1,precr2,roqimax,bvts1,             &
    45. 

  45.              bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs,     &
    46. 

  46.              g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
    47. 

  47.              pidn0s,xlv1,pacrc,pi,                    &
    48. 

  48.              bvtg1,bvtg2,bvtg3,bvtg4,g1pbg,           &
    49. 

  49.              g3pbg,g4pbg,g5pbgo2,pvtg,pacrg,          &
    50. 

  50.              precg1,precg2,pidn0g,                    &
    51. 

  51.              rslopermax,rslopesmax,rslopegmax,        &
    52. 

  52.              rsloperbmax,rslopesbmax,rslopegbmax,     &
    53. 

  53.              rsloper2max,rslopes2max,rslopeg2max,     &
    54. 

  54.              rsloper3max,rslopes3max,rslopeg3max
    55. 

  55. CONTAINS
    56. 

  56. !===================================================================
    57. 

  57. !
    58. 

  58.   SUBROUTINE wsm6(th, q, qc, qr, qi, qs, qg                        &
    59. 

  59.                  ,den, pii, p, delz                                &
    60. 

  60.                  ,delt,g, cpd, cpv, rd, rv, t0c                    &
    61. 

  61.                  ,ep1, ep2, qmin                                   &
    62. 

  62.                  ,XLS, XLV0, XLF0, den0, denr                      &
    63. 

  63.                  ,cliq,cice,psat                                   &
    64. 

  64.                  ,rain, rainncv                                    &
    65. 

  65.                  ,snow, snowncv                                    &
    66. 

  66.                  ,sr                                               &
    67. 

  67.                  ,graupel, graupelncv                              &
    68. 

  68.                  ,ids,ide, jds,jde, kds,kde                        &
    69. 

  69.                  ,ims,ime, jms,jme, kms,kme                        &
    70. 

  70.                  ,its,ite, jts,jte, kts,kte                        &
    71. 

  71.                                                                    )
    72. 

  72. !-------------------------------------------------------------------
    73. 

  73.   IMPLICIT NONE
    74. 

  74. !-------------------------------------------------------------------
    75. 

  75.   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
    76. 

  76.                                       ims,ime, jms,jme, kms,kme , &
    77. 

  77.                                       its,ite, jts,jte, kts,kte
    78. 

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

  79.         INTENT(INOUT) ::                                          &
    80. 

  80.                                                              th,  &
    81. 

  81.                                                               q,  &
    82. 

  82.                                                               qc, &
    83. 

  83.                                                               qi, &
    84. 

  84.                                                               qr, &
    85. 

  85.                                                               qs, &
    86. 

  86.                                                               qg
    87. 

  87.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    88. 

  88.         INTENT(IN   ) ::                                          &
    89. 

  89.                                                              den, &
    90. 

  90.                                                              pii, &
    91. 

  91.                                                                p, &
    92. 

  92.                                                             delz
    93. 

  93.   REAL, INTENT(IN   ) ::                                    delt, &
    94. 

  94.                                                                g, &
    95. 

  95.                                                               rd, &
    96. 

  96.                                                               rv, &
    97. 

  97.                                                              t0c, &
    98. 

  98.                                                             den0, &
    99. 

  99.                                                              cpd, &
    100. 

  100.                                                              cpv, &
    101. 

  101.                                                              ep1, &
    102. 

  102.                                                              ep2, &
    103. 

  103.                                                             qmin, &
    104. 

  104.                                                              XLS, &
    105. 

  105.                                                             XLV0, &
    106. 

  106.                                                             XLF0, &
    107. 

  107.                                                             cliq, &
    108. 

  108.                                                             cice, &
    109. 

  109.                                                             psat, &
    110. 

  110.                                                             denr
    111. 

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

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

  113.                                                          rainncv, &
    114. 

  114.                                                               sr
    115. 

  115.   REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                 &
    116. 

  116.         INTENT(INOUT) ::                                    snow, &
    117. 

  117.                                                          snowncv
    118. 

  118.   REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                 &
    119. 

  119.         INTENT(INOUT) ::                                 graupel, &
    120. 

  120.                                                       graupelncv
    121. 

  121. ! LOCAL VAR
    122. 

  122.   REAL, DIMENSION( its:ite , kts:kte ) ::   t
    123. 

  123.   REAL, DIMENSION( its:ite , kts:kte, 2 ) ::   qci
    124. 

  124.   REAL, DIMENSION( its:ite , kts:kte, 3 ) ::   qrs
    125. 

  125.   INTEGER ::               i,j,k
    126. 

  126. !-------------------------------------------------------------------
    127. 

  127.       DO j=jts,jte
    128. 

  128.          DO k=kts,kte
    129. 

  129.          DO i=its,ite
    130. 

  130.             t(i,k)=th(i,k,j)*pii(i,k,j)
    131. 

  131.             qci(i,k,1) = qc(i,k,j)
    132. 

  132.             qci(i,k,2) = qi(i,k,j)
    133. 

  133.             qrs(i,k,1) = qr(i,k,j)
    134. 

  134.             qrs(i,k,2) = qs(i,k,j)
    135. 

  135.             qrs(i,k,3) = qg(i,k,j)
    136. 

  136.          ENDDO
    137. 

  137.          ENDDO
    138. 

  138.          !  Sending array starting locations of optional variables may cause
    139. 

  139.          !  troubles, so we explicitly change the call.
    140. 

  140.          CALL wsm62D(t, q(ims,kms,j), qci, qrs                     &
    141. 

  141.                     ,den(ims,kms,j)                                &
    142. 

  142.                     ,p(ims,kms,j), delz(ims,kms,j)                 &
    143. 

  143.                     ,delt,g, cpd, cpv, rd, rv, t0c                 &
    144. 

  144.                     ,ep1, ep2, qmin                                &
    145. 

  145.                     ,XLS, XLV0, XLF0, den0, denr                   &
    146. 

  146.                     ,cliq,cice,psat                                &
    147. 

  147.                     ,j                                             &
    148. 

  148.                     ,rain(ims,j),rainncv(ims,j)                    &
    149. 

  149.                     ,sr(ims,j)                                     &
    150. 

  150.                     ,ids,ide, jds,jde, kds,kde                     &
    151. 

  151.                     ,ims,ime, jms,jme, kms,kme                     &
    152. 

  152.                     ,its,ite, jts,jte, kts,kte                     &
    153. 

  153.                     ,snow,snowncv                                  &
    154. 

  154.                     ,graupel,graupelncv                            &
    155. 

  155.                                                                    )
    156. 

  156.          DO K=kts,kte
    157. 

  157.          DO I=its,ite
    158. 

  158.             th(i,k,j)=t(i,k)/pii(i,k,j)
    159. 

  159.             qc(i,k,j) = qci(i,k,1)
    160. 

  160.             qi(i,k,j) = qci(i,k,2)
    161. 

  161.             qr(i,k,j) = qrs(i,k,1)
    162. 

  162.             qs(i,k,j) = qrs(i,k,2)
    163. 

  163.             qg(i,k,j) = qrs(i,k,3)
    164. 

  164.          ENDDO
    165. 

  165.          ENDDO
    166. 

  166.       ENDDO
    167. 

  167.   END SUBROUTINE wsm6
    168. 

  168. !===================================================================
    169. 

  169. !
    170. 

  170.   SUBROUTINE wsm62D(t, q                                          &   
    171. 

  171.                    ,qci, qrs, den, p, delz                        &
    172. 

  172.                    ,delt,g, cpd, cpv, rd, rv, t0c                 &
    173. 

  173.                    ,ep1, ep2, qmin                                &
    174. 

  174.                    ,XLS, XLV0, XLF0, den0, denr                   &
    175. 

  175.                    ,cliq,cice,psat                                &
    176. 

  176.                    ,lat                                           &
    177. 

  177.                    ,rain,rainncv                                  &
    178. 

  178.                    ,sr                                            &
    179. 

  179.                    ,ids,ide, jds,jde, kds,kde                     &
    180. 

  180.                    ,ims,ime, jms,jme, kms,kme                     &
    181. 

  181.                    ,its,ite, jts,jte, kts,kte                     &
    182. 

  182.                    ,snow,snowncv                                  &
    183. 

  183.                    ,graupel,graupelncv                            &
    184. 

  184.                                                                   )
    185. 

  185. !-------------------------------------------------------------------
    186. 

  186.   IMPLICIT NONE
    187. 

  187. !-------------------------------------------------------------------
    188. 

  188. !
    189. 

  189. !  This code is a 6-class GRAUPEL phase microphyiscs scheme (WSM6) of the 
    190. 

  190. !  Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
    191. 

  191. !  number concentration is a function of temperature, and seperate assumption
    192. 

  192. !  is developed, in which ice crystal number concentration is a function
    193. 

  193. !  of ice amount. A theoretical background of the ice-microphysics and related
    194. 

  194. !  processes in the WSMMPs are described in Hong et al. (2004).
    195. 

  195. !  All production terms in the WSM6 scheme are described in Hong and Lim (2006).
    196. 

  196. !  All units are in m.k.s. and source/sink terms in kgkg-1s-1.
    197. 

  197. !
    198. 

  198. !  WSM6 cloud scheme
    199. 

  199. !
    200. 

  200. !  Coded by Song-You Hong and Jeong-Ock Jade Lim (Yonsei Univ.)
    201. 

  201. !           Summer 2003
    202. 

  202. !
    203. 

  203. !  Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
    204. 

  204. !           Summer 2004
    205. 

  205. !
    206. 

  206. !  History :  semi-lagrangian scheme sedimentation(JH), and clean up
    207. 

  207. !             Hong, August 2009
    208. 

  208. !
    209. 

  209. !  Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
    210. 

  210. !             Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
    211. 

  211. !             Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
    212. 

  212. !             Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
    213. 

  213. !             Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
    214. 

  214. !             Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
    215. 

  215. !             Juang and Hong (JH, 2010) Mon. Wea. Rev.
    216. 

  216. !
    217. 

  217.   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
    218. 

  218.                                       ims,ime, jms,jme, kms,kme , &
    219. 

  219.                                       its,ite, jts,jte, kts,kte,  &
    220. 

  220.                                       lat
    221. 

  221.   REAL, DIMENSION( its:ite , kts:kte ),                           &
    222. 

  222.         INTENT(INOUT) ::                                          &
    223. 

  223.                                                                t
    224. 

  224.   REAL, DIMENSION( its:ite , kts:kte, 2 ),                        &
    225. 

  225.         INTENT(INOUT) ::                                          &
    226. 

  226.                                                              qci
    227. 

  227.   REAL, DIMENSION( its:ite , kts:kte, 3 ),                        &
    228. 

  228.         INTENT(INOUT) ::                                          &
    229. 

  229.                                                              qrs
    230. 

  230.   REAL, DIMENSION( ims:ime , kms:kme ),                           &
    231. 

  231.         INTENT(INOUT) ::                                          &
    232. 

  232.                                                                q
    233. 

  233.   REAL, DIMENSION( ims:ime , kms:kme ),                           &
    234. 

  234.         INTENT(IN   ) ::                                          &
    235. 

  235.                                                              den, &
    236. 

  236.                                                                p, &
    237. 

  237.                                                             delz
    238. 

  238.   REAL, INTENT(IN   ) ::                                    delt, &
    239. 

  239.                                                                g, &
    240. 

  240.                                                              cpd, &
    241. 

  241.                                                              cpv, &
    242. 

  242.                                                              t0c, &
    243. 

  243.                                                             den0, &
    244. 

  244.                                                               rd, &
    245. 

  245.                                                               rv, &
    246. 

  246.                                                              ep1, &
    247. 

  247.                                                              ep2, &
    248. 

  248.                                                             qmin, &
    249. 

  249.                                                              XLS, &
    250. 

  250.                                                             XLV0, &
    251. 

  251.                                                             XLF0, &
    252. 

  252.                                                             cliq, &
    253. 

  253.                                                             cice, &
    254. 

  254.                                                             psat, &
    255. 

  255.                                                             denr
    256. 

  256.   REAL, DIMENSION( ims:ime ),                                     &
    257. 

  257.         INTENT(INOUT) ::                                    rain, &
    258. 

  258.                                                          rainncv, &
    259. 

  259.                                                               sr
    260. 

  260.   REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL,                  &
    261. 

  261.         INTENT(INOUT) ::                                    snow, &
    262. 

  262.                                                          snowncv
    263. 

  263.   REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL,                  &
    264. 

  264.         INTENT(INOUT) ::                                 graupel, &
    265. 

  265.                                                       graupelncv
    266. 

  266. ! LOCAL VAR
    267. 

  267.   REAL, DIMENSION( its:ite , kts:kte , 3) ::                      &
    268. 

  268.                                                               rh, &
    269. 

  269.                                                               qs, &
    270. 

  270.                                                           rslope, &
    271. 

  271.                                                          rslope2, &
    272. 

  272.                                                          rslope3, &
    273. 

  273.                                                          rslopeb, &
    274. 

  274.                                                          qrs_tmp, & 
    275. 

  275.                                                             falk, &
    276. 

  276.                                                             fall, &
    277. 

  277.                                                            work1
    278. 

  278.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    279. 

  279.                                                            fallc, &
    280. 

  280.                                                            falkc, &
    281. 

  281.                                                           work1c, &
    282. 

  282.                                                           work2c, &
    283. 

  283.                                                            workr, &
    284. 

  284.                                                            worka
    285. 

  285.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    286. 

  286.                                                          den_tmp, &
    287. 

  287.                                                         delz_tmp
    288. 

  288.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    289. 

  289.                                                            pigen, &
    290. 

  290.                                                            pidep, &
    291. 

  291.                                                            pcond, &
    292. 

  292.                                                            prevp, &
    293. 

  293.                                                            psevp, &
    294. 

  294.                                                            pgevp, &
    295. 

  295.                                                            psdep, &
    296. 

  296.                                                            pgdep, &
    297. 

  297.                                                            praut, &
    298. 

  298.                                                            psaut, &
    299. 

  299.                                                            pgaut, &
    300. 

  300.                                                            piacr, &
    301. 

  301.                                                            pracw, &
    302. 

  302.                                                            praci, &
    303. 

  303.                                                            pracs, &
    304. 

  304.                                                            psacw, &
    305. 

  305.                                                            psaci, &
    306. 

  306.                                                            psacr, &
    307. 

  307.                                                            pgacw, &
    308. 

  308.                                                            pgaci, &
    309. 

  309.                                                            pgacr, &
    310. 

  310.                                                            pgacs, &
    311. 

  311.                                                            paacw, &
    312. 

  312.                                                            psmlt, &
    313. 

  313.                                                            pgmlt, &
    314. 

  314.                                                            pseml, &
    315. 

  315.                                                            pgeml
    316. 

  316.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    317. 

  317.                                                             qsum, &
    318. 

  318.                                                               xl, &
    319. 

  319.                                                              cpm, &
    320. 

  320.                                                            work2, &
    321. 

  321.                                                           denfac, &
    322. 

  322.                                                              xni, &
    323. 

  323.                                                          denqrs1, &
    324. 

  324.                                                          denqrs2, &
    325. 

  325.                                                          denqrs3, &
    326. 

  326.                                                           denqci, & 
    327. 

  327.                                                           n0sfac
    328. 

  328.   REAL, DIMENSION( its:ite ) ::                          delqrs1, &
    329. 

  329.                                                          delqrs2, &
    330. 

  330.                                                          delqrs3, &
    331. 

  331.                                                            delqi  
    332. 

  332.   REAL, DIMENSION( its:ite ) ::                        tstepsnow, &
    333. 

  333.                                                       tstepgraup
    334. 

  334.   INTEGER, DIMENSION( its:ite ) ::                         mstep, &
    335. 

  335.                                                            numdt
    336. 

  336.   LOGICAL, DIMENSION( its:ite ) ::                        flgcld
    337. 

  337.   REAL  ::                                                        &
    338. 

  338.             cpmcal, xlcal, diffus,                                &
    339. 

  339.             viscos, xka, venfac, conden, diffac,                  &
    340. 

  340.             x, y, z, a, b, c, d, e,                               &
    341. 

  341.             qdt, holdrr, holdrs, holdrg, supcol, supcolt, pvt,    &
    342. 

  342.             coeres, supsat, dtcld, xmi, eacrs, satdt,             &
    343. 

  343.             qimax, diameter, xni0, roqi0,                         &
    344. 

  344.             fallsum, fallsum_qsi, fallsum_qg,                     &
    345. 

  345.             vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi,                   &
    346. 

  346.             xlwork2, factor, source, value,                       &
    347. 

  347.             xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3  
    348. 

  348.   REAL  :: vt2ave
    349. 

  349.   REAL  :: holdc, holdci
    350. 

  350.   INTEGER :: i, j, k, mstepmax,                                   &
    351. 

  351.             iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
    352. 

  352. ! Temporaries used for inlining fpvs function
    353. 

  353.   REAL  :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
    354. 

  354. ! variables for optimization
    355. 

  355.   REAL, DIMENSION( its:ite ) ::                             tvec1
    356. 

  356.   REAL                       ::                              temp
    357. 

  357. !
    358. 

  358. !=================================================================
    359. 

  359. !   compute internal functions
    360. 

  360. !
    361. 

  361.       cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
    362. 

  362.       xlcal(x) = xlv0-xlv1*(x-t0c)
    363. 

  363. !----------------------------------------------------------------
    364. 

  364. !     diffus: diffusion coefficient of the water vapor
    365. 

  365. !     viscos: kinematic viscosity(m2s-1)
    366. 

  366. !     Optimizatin : A**B => exp(log(A)*(B))
    367. 

  367. !
    368. 

  368.       diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y        ! 8.794e-5*x**1.81/y
    369. 

  369.       viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y  ! 1.496e-6*x**1.5/(x+120.)/y
    370. 

  370.       xka(x,y) = 1.414e3*viscos(x,y)*y
    371. 

  371.       diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
    372. 

  372.       venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333)))         &
    373. 

  373.                      /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
    374. 

  374.       conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
    375. 

  375. !
    376. 

  376. !
    377. 

  377.       idim = ite-its+1
    378. 

  378.       kdim = kte-kts+1
    379. 

  379. !
    380. 

  380. !----------------------------------------------------------------
    381. 

  381. !     paddint 0 for negative values generated by dynamics
    382. 

  382. !
    383. 

  383.       do k = kts, kte
    384. 

  384.         do i = its, ite
    385. 

  385.           qci(i,k,1) = max(qci(i,k,1),0.0)
    386. 

  386.           qrs(i,k,1) = max(qrs(i,k,1),0.0)
    387. 

  387.           qci(i,k,2) = max(qci(i,k,2),0.0)
    388. 

  388.           qrs(i,k,2) = max(qrs(i,k,2),0.0)
    389. 

  389.           qrs(i,k,3) = max(qrs(i,k,3),0.0)
    390. 

  390.         enddo
    391. 

  391.       enddo
    392. 

  392. !
    393. 

  393. !----------------------------------------------------------------
    394. 

  394. !     latent heat for phase changes and heat capacity. neglect the
    395. 

  395. !     changes during microphysical process calculation
    396. 

  396. !     emanuel(1994)
    397. 

  397. !
    398. 

  398.       do k = kts, kte
    399. 

  399.         do i = its, ite
    400. 

  400.           cpm(i,k) = cpmcal(q(i,k))
    401. 

  401.           xl(i,k) = xlcal(t(i,k))
    402. 

  402.         enddo
    403. 

  403.       enddo
    404. 

  404.       do k = kts, kte
    405. 

  405.         do i = its, ite
    406. 

  406.           delz_tmp(i,k) = delz(i,k)
    407. 

  407.           den_tmp(i,k) = den(i,k)
    408. 

  408.         enddo
    409. 

  409.       enddo
    410. 

  410. !
    411. 

  411. !----------------------------------------------------------------
    412. 

  412. !    initialize the surface rain, snow, graupel
    413. 

  413. !
    414. 

  414.       do i = its, ite
    415. 

  415.         rainncv(i) = 0.
    416. 

  416.         if(PRESENT (snowncv) .AND. PRESENT (snow)) snowncv(i,lat) = 0.
    417. 

  417.         if(PRESENT (graupelncv) .AND. PRESENT (graupel)) graupelncv(i,lat) = 0.
    418. 

  418.         sr(i) = 0.
    419. 

  419. ! new local array to catch step snow and graupel
    420. 

  420.         tstepsnow(i) = 0.
    421. 

  421.         tstepgraup(i) = 0.
    422. 

  422.       enddo
    423. 

  423. !
    424. 

  424. !----------------------------------------------------------------
    425. 

  425. !     compute the minor time steps.
    426. 

  426. !
    427. 

  427.       loops = max(nint(delt/dtcldcr),1)
    428. 

  428.       dtcld = delt/loops
    429. 

  429.       if(delt.le.dtcldcr) dtcld = delt
    430. 

  430. !
    431. 

  431.       do loop = 1,loops
    432. 

  432. !
    433. 

  433. !----------------------------------------------------------------
    434. 

  434. !     initialize the large scale variables
    435. 

  435. !
    436. 

  436.       do i = its, ite
    437. 

  437.         mstep(i) = 1
    438. 

  438.         flgcld(i) = .true.
    439. 

  439.       enddo
    440. 

  440. !
    441. 

  441. !     do k = kts, kte
    442. 

  442. !       do i = its, ite
    443. 

  443. !         denfac(i,k) = sqrt(den0/den(i,k))
    444. 

  444. !       enddo
    445. 

  445. !     enddo
    446. 

  446.       do k = kts, kte
    447. 

  447.         CALL VREC( tvec1(its), den(its,k), ite-its+1)
    448. 

  448.         do i = its, ite
    449. 

  449.           tvec1(i) = tvec1(i)*den0
    450. 

  450.         enddo
    451. 

  451.         CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
    452. 

  452.       enddo
    453. 

  453. !
    454. 

  454. ! Inline expansion for fpvs
    455. 

  455. !         qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
    456. 

  456. !         qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
    457. 

  457.       hsub = xls
    458. 

  458.       hvap = xlv0
    459. 

  459.       cvap = cpv
    460. 

  460.       ttp=t0c+0.01
    461. 

  461.       dldt=cvap-cliq
    462. 

  462.       xa=-dldt/rv
    463. 

  463.       xb=xa+hvap/(rv*ttp)
    464. 

  464.       dldti=cvap-cice
    465. 

  465.       xai=-dldti/rv
    466. 

  466.       xbi=xai+hsub/(rv*ttp)
    467. 

  467.       do k = kts, kte
    468. 

  468.         do i = its, ite
    469. 

  469.           tr=ttp/t(i,k)
    470. 

  470.           qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
    471. 

  471.           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
    472. 

  472.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    473. 

  473.           qs(i,k,1) = max(qs(i,k,1),qmin)
    474. 

  474.           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
    475. 

  475.           tr=ttp/t(i,k)
    476. 

  476.           if(t(i,k).lt.ttp) then
    477. 

  477.             qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
    478. 

  478.           else
    479. 

  479.             qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
    480. 

  480.           endif
    481. 

  481.           qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
    482. 

  482.           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
    483. 

  483.           qs(i,k,2) = max(qs(i,k,2),qmin)
    484. 

  484.           rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
    485. 

  485.         enddo
    486. 

  486.       enddo
    487. 

  487. !
    488. 

  488. !----------------------------------------------------------------
    489. 

  489. !     initialize the variables for microphysical physics
    490. 

  490. !
    491. 

  491. !
    492. 

  492.       do k = kts, kte
    493. 

  493.         do i = its, ite
    494. 

  494.           prevp(i,k) = 0.
    495. 

  495.           psdep(i,k) = 0.
    496. 

  496.           pgdep(i,k) = 0.
    497. 

  497.           praut(i,k) = 0.
    498. 

  498.           psaut(i,k) = 0.
    499. 

  499.           pgaut(i,k) = 0.
    500. 

  500.           pracw(i,k) = 0.
    501. 

  501.           praci(i,k) = 0.
    502. 

  502.           piacr(i,k) = 0.
    503. 

  503.           psaci(i,k) = 0.
    504. 

  504.           psacw(i,k) = 0.
    505. 

  505.           pracs(i,k) = 0.
    506. 

  506.           psacr(i,k) = 0.
    507. 

  507.           pgacw(i,k) = 0.
    508. 

  508.           paacw(i,k) = 0.
    509. 

  509.           pgaci(i,k) = 0.
    510. 

  510.           pgacr(i,k) = 0.
    511. 

  511.           pgacs(i,k) = 0.
    512. 

  512.           pigen(i,k) = 0.
    513. 

  513.           pidep(i,k) = 0.
    514. 

  514.           pcond(i,k) = 0.
    515. 

  515.           psmlt(i,k) = 0.
    516. 

  516.           pgmlt(i,k) = 0.
    517. 

  517.           pseml(i,k) = 0.
    518. 

  518.           pgeml(i,k) = 0.
    519. 

  519.           psevp(i,k) = 0.
    520. 

  520.           pgevp(i,k) = 0.
    521. 

  521.           falk(i,k,1) = 0.
    522. 

  522.           falk(i,k,2) = 0.
    523. 

  523.           falk(i,k,3) = 0.
    524. 

  524.           fall(i,k,1) = 0.
    525. 

  525.           fall(i,k,2) = 0.
    526. 

  526.           fall(i,k,3) = 0.
    527. 

  527.           fallc(i,k) = 0.
    528. 

  528.           falkc(i,k) = 0.
    529. 

  529.           xni(i,k) = 1.e3
    530. 

  530.         enddo
    531. 

  531.       enddo
    532. 

  532. !-------------------------------------------------------------
    533. 

  533. ! Ni: ice crystal number concentraiton   [HDC 5c]
    534. 

  534. !-------------------------------------------------------------
    535. 

  535.       do k = kts, kte
    536. 

  536.         do i = its, ite
    537. 

  537.           temp = (den(i,k)*max(qci(i,k,2),qmin))
    538. 

  538.           temp = sqrt(sqrt(temp*temp*temp))
    539. 

  539.           xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
    540. 

  540.         enddo
    541. 

  541.       enddo
    542. 

  542. !
    543. 

  543. !----------------------------------------------------------------
    544. 

  544. !     compute the fallout term:
    545. 

  545. !     first, vertical terminal velosity for minor loops
    546. 

  546. !----------------------------------------------------------------
    547. 

  547.       do k = kts, kte
    548. 

  548.         do i = its, ite
    549. 

  549.           qrs_tmp(i,k,1) = qrs(i,k,1)
    550. 

  550.           qrs_tmp(i,k,2) = qrs(i,k,2)
    551. 

  551.           qrs_tmp(i,k,3) = qrs(i,k,3)
    552. 

  552.         enddo
    553. 

  553.       enddo
    554. 

  554.       call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, & 
    555. 

  555.                      work1,its,ite,kts,kte)
    556. 

  556. !
    557. 

  557.       do k = kte, kts, -1
    558. 

  558.         do i = its, ite
    559. 

  559.           workr(i,k) = work1(i,k,1)
    560. 

  560.           qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
    561. 

  561.           IF ( qsum(i,k) .gt. 1.e-15 ) THEN
    562. 

  562.             worka(i,k) = (work1(i,k,2)*qrs(i,k,2) + work1(i,k,3)*qrs(i,k,3)) &
    563. 

  563.                       /qsum(i,k)
    564. 

  564.           ELSE
    565. 

  565.             worka(i,k) = 0.
    566. 

  566.           ENDIF
    567. 

  567.           denqrs1(i,k) = den(i,k)*qrs(i,k,1)
    568. 

  568.           denqrs2(i,k) = den(i,k)*qrs(i,k,2)
    569. 

  569.           denqrs3(i,k) = den(i,k)*qrs(i,k,3)
    570. 

  570.           if(qrs(i,k,1).le.0.0) workr(i,k) = 0.0
    571. 

  571.         enddo
    572. 

  572.       enddo
    573. 

  573.       call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,workr,denqrs1,  &
    574. 

  574.                            delqrs1,dtcld,1,1)
    575. 

  575.       call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka,         & 
    576. 

  576.                            denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1)
    577. 

  577.       do k = kts, kte
    578. 

  578.         do i = its, ite
    579. 

  579.           qrs(i,k,1) = max(denqrs1(i,k)/den(i,k),0.)
    580. 

  580.           qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
    581. 

  581.           qrs(i,k,3) = max(denqrs3(i,k)/den(i,k),0.)
    582. 

  582.           fall(i,k,1) = denqrs1(i,k)*workr(i,k)/delz(i,k)
    583. 

  583.           fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k)
    584. 

  584.           fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k)
    585. 

  585.         enddo
    586. 

  586.       enddo
    587. 

  587.       do i = its, ite
    588. 

  588.         fall(i,1,1) = delqrs1(i)/delz(i,1)/dtcld
    589. 

  589.         fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
    590. 

  590.         fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld
    591. 

  591.       enddo
    592. 

  592.       do k = kts, kte
    593. 

  593.         do i = its, ite
    594. 

  594.           qrs_tmp(i,k,1) = qrs(i,k,1)
    595. 

  595.           qrs_tmp(i,k,2) = qrs(i,k,2)
    596. 

  596.           qrs_tmp(i,k,3) = qrs(i,k,3)
    597. 

  597.         enddo
    598. 

  598.       enddo
    599. 

  599.       call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
    600. 

  600.                      work1,its,ite,kts,kte)
    601. 

  601. !
    602. 

  602.       do k = kte, kts, -1 
    603. 

  603.         do i = its, ite
    604. 

  604.           supcol = t0c-t(i,k)
    605. 

  605.           n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
    606. 

  606.           if(t(i,k).gt.t0c) then
    607. 

  607. !---------------------------------------------------------------
    608. 

  608. ! psmlt: melting of snow [HL A33] [RH83 A25]
    609. 

  609. !       (T>T0: S->R)
    610. 

  610. !---------------------------------------------------------------
    611. 

  611.             xlf = xlf0
    612. 

  612.             work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
    613. 

  613.             if(qrs(i,k,2).gt.0.) then
    614. 

  614.               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    615. 

  615.               psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2.       &
    616. 

  616.                          *n0sfac(i,k)*(precs1*rslope2(i,k,2)                 &
    617. 

  617.                          +precs2*work2(i,k)*coeres)
    618. 

  618.               psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),                &
    619. 

  619.                           -qrs(i,k,2)/mstep(i)),0.)
    620. 

  620.               qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
    621. 

  621.               qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
    622. 

  622.               t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
    623. 

  623.             endif
    624. 

  624. !---------------------------------------------------------------
    625. 

  625. ! pgmlt: melting of graupel [HL A23]  [LFO 47]
    626. 

  626. !       (T>T0: G->R)
    627. 

  627. !---------------------------------------------------------------
    628. 

  628.             if(qrs(i,k,3).gt.0.) then
    629. 

  629.               coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
    630. 

  630.               pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf                          &
    631. 

  631.                          *(t0c-t(i,k))*(precg1*rslope2(i,k,3)                &
    632. 

  632.                          +precg2*work2(i,k)*coeres)
    633. 

  633.               pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i),                &
    634. 

  634.                           -qrs(i,k,3)/mstep(i)),0.)                          
    635. 

  635.               qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
    636. 

  636.               qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
    637. 

  637.               t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
    638. 

  638.             endif
    639. 

  639.           endif
    640. 

  640.         enddo
    641. 

  641.       enddo
    642. 

  642. !---------------------------------------------------------------
    643. 

  643. ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
    644. 

  644. !---------------------------------------------------------------
    645. 

  645.       do k = kte, kts, -1
    646. 

  646.         do i = its, ite
    647. 

  647.           if(qci(i,k,2).le.0.) then
    648. 

  648.             work1c(i,k) = 0.
    649. 

  649.           else
    650. 

  650.             xmi = den(i,k)*qci(i,k,2)/xni(i,k)
    651. 

  651.             diameter  = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
    652. 

  652.             work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
    653. 

  653.           endif
    654. 

  654.         enddo
    655. 

  655.       enddo
    656. 

  656. !
    657. 

  657. !  forward semi-laglangian scheme (JH), PCM (piecewise constant),  (linear)
    658. 

  658. !
    659. 

  659.       do k = kte, kts, -1
    660. 

  660.         do i = its, ite
    661. 

  661.           denqci(i,k) = den(i,k)*qci(i,k,2)
    662. 

  662.         enddo
    663. 

  663.       enddo
    664. 

  664.       call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci,  &
    665. 

  665.                            delqi,dtcld,1,0)
    666. 

  666.       do k = kts, kte
    667. 

  667.         do i = its, ite
    668. 

  668.           qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
    669. 

  669.         enddo
    670. 

  670.       enddo
    671. 

  671.       do i = its, ite
    672. 

  672.         fallc(i,1) = delqi(i)/delz(i,1)/dtcld
    673. 

  673.       enddo
    674. 

  674. !
    675. 

  675. !----------------------------------------------------------------
    676. 

  676. !      rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
    677. 

  677. !
    678. 

  678.       do i = its, ite
    679. 

  679.         fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
    680. 

  680.         fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
    681. 

  681.         fallsum_qg = fall(i,kts,3)
    682. 

  682.         if(fallsum.gt.0.) then
    683. 

  683.           rainncv(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rainncv(i)
    684. 

  684.           rain(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rain(i)
    685. 

  685.         endif
    686. 

  686.         if(fallsum_qsi.gt.0.) then
    687. 

  687.           tstepsnow(i)   = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.            &
    688. 

  688.                            +tstepsnow(i)
    689. 

  689.         IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
    690. 

  690.           snowncv(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.            & 
    691. 

  691.                            +snowncv(i,lat)
    692. 

  692.           snow(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i,lat)
    693. 

  693.         ENDIF
    694. 

  694.         endif
    695. 

  695.         if(fallsum_qg.gt.0.) then
    696. 

  696.           tstepgraup(i)  = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.            &
    697. 

  697.                            +tstepgraup(i)
    698. 

  698.         IF ( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN
    699. 

  699.           graupelncv(i,lat) = fallsum_qg*delz(i,kts)/denr*dtcld*1000.          &   
    700. 

  700.                               + graupelncv(i,lat)
    701. 

  701.           graupel(i,lat) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. + graupel(i,lat)
    702. 

  702.         ENDIF
    703. 

  703.         endif
    704. 

  704. !       if(fallsum.gt.0.)sr(i)=(snowncv(i,lat) + graupelncv(i,lat))/(rainncv(i)+1.e-12)
    705. 

  705.         if(fallsum.gt.0.)sr(i)=(tstepsnow(i) + tstepgraup(i))/(rainncv(i)+1.e-12)
    706. 

  706.       enddo
    707. 

  707. !
    708. 

  708. !---------------------------------------------------------------
    709. 

  709. ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
    710. 

  710. !       (T>T0: I->C)
    711. 

  711. !---------------------------------------------------------------
    712. 

  712.       do k = kts, kte
    713. 

  713.         do i = its, ite
    714. 

  714.           supcol = t0c-t(i,k)
    715. 

  715.           xlf = xls-xl(i,k)
    716. 

  716.           if(supcol.lt.0.) xlf = xlf0
    717. 

  717.           if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
    718. 

  718.             qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
    719. 

  719.             t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
    720. 

  720.             qci(i,k,2) = 0.
    721. 

  721.           endif
    722. 

  722. !---------------------------------------------------------------
    723. 

  723. ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
    724. 

  724. !        (T<-40C: C->I)
    725. 

  725. !---------------------------------------------------------------
    726. 

  726.           if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
    727. 

  727.             qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
    728. 

  728.             t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
    729. 

  729.             qci(i,k,1) = 0.
    730. 

  730.           endif
    731. 

  731. !---------------------------------------------------------------
    732. 

  732. ! pihtf: heterogeneous freezing of cloud water [HL A44]
    733. 

  733. !        (T0>T>-40C: C->I)
    734. 

  734. !---------------------------------------------------------------
    735. 

  735.           if(supcol.gt.0..and.qci(i,k,1).gt.qmin) then
    736. 

  736. !           pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.)                         &
    737. 

  737. !              *den(i,k)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
    738. 

  738.             supcolt=min(supcol,50.)
    739. 

  739.             pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.)                        &
    740. 

  740.             *den(i,k)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
    741. 

  741.             qci(i,k,2) = qci(i,k,2) + pfrzdtc
    742. 

  742.             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
    743. 

  743.             qci(i,k,1) = qci(i,k,1)-pfrzdtc
    744. 

  744.           endif
    745. 

  745. !---------------------------------------------------------------
    746. 

  746. ! pgfrz: freezing of rain water [HL A20] [LFO 45]
    747. 

  747. !        (T<T0, R->G)
    748. 

  748. !---------------------------------------------------------------
    749. 

  749.           if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
    750. 

  750. !           pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k)                    &
    751. 

  751. !                 *(exp(pfrz2*supcol)-1.)*rslope3(i,k,1)**2                    &
    752. 

  752. !                 *rslope(i,k,1)*dtcld,qrs(i,k,1))
    753. 

  753.             temp = rslope3(i,k,1)
    754. 

  754.             temp = temp*temp*rslope(i,k,1)
    755. 

  755.             supcolt=min(supcol,50.)
    756. 

  756.             pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k)                  &
    757. 

  757.                   *(exp(pfrz2*supcolt)-1.)*temp*dtcld,                         &
    758. 

  758.                   qrs(i,k,1))
    759. 

  759.             qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
    760. 

  760.             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
    761. 

  761.             qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
    762. 

  762.           endif
    763. 

  763.         enddo
    764. 

  764.       enddo
    765. 

  765. !
    766. 

  766. !
    767. 

  767. !----------------------------------------------------------------
    768. 

  768. !     update the slope parameters for microphysics computation
    769. 

  769. !
    770. 

  770.       do k = kts, kte
    771. 

  771.         do i = its, ite
    772. 

  772.           qrs_tmp(i,k,1) = qrs(i,k,1)
    773. 

  773.           qrs_tmp(i,k,2) = qrs(i,k,2)
    774. 

  774.           qrs_tmp(i,k,3) = qrs(i,k,3)
    775. 

  775.         enddo
    776. 

  776.       enddo
    777. 

  777.       call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
    778. 

  778.                      work1,its,ite,kts,kte)
    779. 

  779. !------------------------------------------------------------------
    780. 

  780. !     work1:  the thermodynamic term in the denominator associated with
    781. 

  781. !             heat conduction and vapor diffusion
    782. 

  782. !             (ry88, y93, h85)
    783. 

  783. !     work2: parameter associated with the ventilation effects(y93)
    784. 

  784. !
    785. 

  785.       do k = kts, kte
    786. 

  786.         do i = its, ite
    787. 

  787.           work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
    788. 

  788.           work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
    789. 

  789.           work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
    790. 

  790.         enddo
    791. 

  791.       enddo
    792. 

  792. !
    793. 

  793. !===============================================================
    794. 

  794. !
    795. 

  795. ! warm rain processes
    796. 

  796. !
    797. 

  797. ! - follows the processes in RH83 and LFO except for autoconcersion
    798. 

  798. !
    799. 

  799. !===============================================================
    800. 

  800. !
    801. 

  801.       do k = kts, kte
    802. 

  802.         do i = its, ite
    803. 

  803.           supsat = max(q(i,k),qmin)-qs(i,k,1)
    804. 

  804.           satdt = supsat/dtcld
    805. 

  805. !---------------------------------------------------------------
    806. 

  806. ! praut: auto conversion rate from cloud to rain [HDC 16]
    807. 

  807. !        (C->R)
    808. 

  808. !---------------------------------------------------------------
    809. 

  809.           if(qci(i,k,1).gt.qc0) then
    810. 

  810.             praut(i,k) = qck1*qci(i,k,1)**(7./3.)
    811. 

  811.             praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
    812. 

  812.           endif
    813. 

  813. !---------------------------------------------------------------
    814. 

  814. ! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
    815. 

  815. !        (C->R)
    816. 

  816. !---------------------------------------------------------------
    817. 

  817.           if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
    818. 

  818.             pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1)               &
    819. 

  819.                         *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
    820. 

  820.           endif
    821. 

  821. !---------------------------------------------------------------
    822. 

  822. ! prevp: evaporation/condensation rate of rain [HDC 14]
    823. 

  823. !        (V->R or R->V)
    824. 

  824. !---------------------------------------------------------------
    825. 

  825.           if(qrs(i,k,1).gt.0.) then
    826. 

  826.             coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
    827. 

  827.             prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1)                 &
    828. 

  828.                          +precr2*work2(i,k)*coeres)/work1(i,k,1)
    829. 

  829.             if(prevp(i,k).lt.0.) then
    830. 

  830.               prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
    831. 

  831.               prevp(i,k) = max(prevp(i,k),satdt/2)
    832. 

  832.             else
    833. 

  833.               prevp(i,k) = min(prevp(i,k),satdt/2)
    834. 

  834.             endif
    835. 

  835.           endif
    836. 

  836.         enddo
    837. 

  837.       enddo
    838. 

  838. !
    839. 

  839. !===============================================================
    840. 

  840. !
    841. 

  841. ! cold rain processes
    842. 

  842. !
    843. 

  843. ! - follows the revised ice microphysics processes in HDC
    844. 

  844. ! - the processes same as in RH83 and RH84  and LFO behave
    845. 

  845. !   following ice crystal hapits defined in HDC, inclduing
    846. 

  846. !   intercept parameter for snow (n0s), ice crystal number
    847. 

  847. !   concentration (ni), ice nuclei number concentration
    848. 

  848. !   (n0i), ice diameter (d)
    849. 

  849. !
    850. 

  850. !===============================================================
    851. 

  851. !
    852. 

  852.       do k = kts, kte
    853. 

  853.         do i = its, ite
    854. 

  854.           supcol = t0c-t(i,k)
    855. 

  855.           n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
    856. 

  856.           supsat = max(q(i,k),qmin)-qs(i,k,2)
    857. 

  857.           satdt = supsat/dtcld
    858. 

  858.           ifsat = 0
    859. 

  859. !-------------------------------------------------------------
    860. 

  860. ! Ni: ice crystal number concentraiton   [HDC 5c]
    861. 

  861. !-------------------------------------------------------------
    862. 

  862. !         xni(i,k) = min(max(5.38e7*(den(i,k)                                  &
    863. 

  863. !                      *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
    864. 

  864.           temp = (den(i,k)*max(qci(i,k,2),qmin))
    865. 

  865.           temp = sqrt(sqrt(temp*temp*temp))
    866. 

  866.           xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
    867. 

  867.           eacrs = exp(0.07*(-supcol))
    868. 

  868. !
    869. 

  869.           xmi = den(i,k)*qci(i,k,2)/xni(i,k)
    870. 

  870.           diameter  = min(dicon * sqrt(xmi),dimax)
    871. 

  871.           vt2i = 1.49e4*diameter**1.31
    872. 

  872.           vt2r=pvtr*rslopeb(i,k,1)*denfac(i,k)
    873. 

  873.           vt2s=pvts*rslopeb(i,k,2)*denfac(i,k)
    874. 

  874.           vt2g=pvtg*rslopeb(i,k,3)*denfac(i,k)
    875. 

  875.           qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
    876. 

  876.           if(qsum(i,k) .gt. 1.e-15) then
    877. 

  877.           vt2ave=(vt2s*qrs(i,k,2)+vt2g*qrs(i,k,3))/(qsum(i,k))
    878. 

  878.           else
    879. 

  879.           vt2ave=0.
    880. 

  880.           endif
    881. 

  881.           if(supcol.gt.0.and.qci(i,k,2).gt.qmin) then
    882. 

  882.             if(qrs(i,k,1).gt.qcrmin) then
    883. 

  883. !-------------------------------------------------------------
    884. 

  884. ! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
    885. 

  885. !        (T<T0: I->R)
    886. 

  886. !-------------------------------------------------------------
    887. 

  887.               acrfac = 2.*rslope3(i,k,1)+2.*diameter*rslope2(i,k,1)            &
    888. 

  888.                       +diameter**2*rslope(i,k,1)
    889. 

  889.               praci(i,k) = pi*qci(i,k,2)*n0r*abs(vt2r-vt2i)*acrfac/4.
    890. 

  890.               praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
    891. 

  891. !-------------------------------------------------------------
    892. 

  892. ! piacr: Accretion of rain by cloud ice [HL A19] [LFO 26]
    893. 

  893. !        (T<T0: R->S or R->G)
    894. 

  894. !-------------------------------------------------------------
    895. 

  895.               piacr(i,k) = pi**2*avtr*n0r*denr*xni(i,k)*denfac(i,k)            &
    896. 

  896.                           *g6pbr*rslope3(i,k,1)*rslope3(i,k,1)                 &
    897. 

  897.                           *rslopeb(i,k,1)/24./den(i,k)
    898. 

  898.               piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
    899. 

  899.             endif
    900. 

  900. !-------------------------------------------------------------
    901. 

  901. ! psaci: Accretion of cloud ice by snow [HDC 10]
    902. 

  902. !        (T<T0: I->S)
    903. 

  903. !-------------------------------------------------------------
    904. 

  904.             if(qrs(i,k,2).gt.qcrmin) then
    905. 

  905.               acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2)            &
    906. 

  906.                       +diameter**2*rslope(i,k,2)
    907. 

  907.               psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k)                 &
    908. 

  908.                           *abs(vt2ave-vt2i)*acrfac/4.
    909. 

  909.               psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
    910. 

  910.             endif
    911. 

  911. !-------------------------------------------------------------
    912. 

  912. ! pgaci: Accretion of cloud ice by graupel [HL A17] [LFO 41]
    913. 

  913. !        (T<T0: I->G)
    914. 

  914. !-------------------------------------------------------------
    915. 

  915.             if(qrs(i,k,3).gt.qcrmin) then
    916. 

  916.               egi = exp(0.07*(-supcol))
    917. 

  917.               acrfac = 2.*rslope3(i,k,3)+2.*diameter*rslope2(i,k,3)            &
    918. 

  918.                       +diameter**2*rslope(i,k,3)
    919. 

  919.               pgaci(i,k) = pi*egi*qci(i,k,2)*n0g*abs(vt2ave-vt2i)*acrfac/4.
    920. 

  920.               pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
    921. 

  921.             endif
    922. 

  922.           endif
    923. 

  923. !-------------------------------------------------------------
    924. 

  924. ! psacw: Accretion of cloud water by snow  [HL A7] [LFO 24]
    925. 

  925. !        (T<T0: C->S, and T>=T0: C->R)
    926. 

  926. !-------------------------------------------------------------
    927. 

  927.           if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
    928. 

  928.             psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2)   &    
    929. 

  929.                         *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
    930. 

  930.           endif
    931. 

  931. !-------------------------------------------------------------
    932. 

  932. ! pgacw: Accretion of cloud water by graupel [HL A6] [LFO 40]
    933. 

  933. !        (T<T0: C->G, and T>=T0: C->R)
    934. 

  934. !-------------------------------------------------------------
    935. 

  935.           if(qrs(i,k,3).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
    936. 

  936.             pgacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3)               &
    937. 

  937.                         *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
    938. 

  938.           endif
    939. 

  939. !-------------------------------------------------------------
    940. 

  940. ! paacw: Accretion of cloud water by averaged snow/graupel 
    941. 

  941. !        (T<T0: C->G or S, and T>=T0: C->R) 
    942. 

  942. !-------------------------------------------------------------
    943. 

  943.           if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,3).gt.qcrmin) then
    944. 

  944.             paacw(i,k) = (qrs(i,k,2)*psacw(i,k)+qrs(i,k,3)*pgacw(i,k))         & 
    945. 

  945.                         /(qsum(i,k))
    946. 

  946.            endif      
    947. 

  947. !-------------------------------------------------------------
    948. 

  948. ! pracs: Accretion of snow by rain [HL A11] [LFO 27]
    949. 

  949. !         (T<T0: S->G)
    950. 

  950. !-------------------------------------------------------------
    951. 

  951.           if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
    952. 

  952.             if(supcol.gt.0) then
    953. 

  953.               acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2)*rslope(i,k,1)          &
    954. 

  954.                       +2.*rslope3(i,k,2)*rslope2(i,k,2)*rslope2(i,k,1)         &
    955. 

  955.                       +.5*rslope2(i,k,2)*rslope2(i,k,2)*rslope3(i,k,1)
    956. 

  956.               pracs(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2r-vt2ave)          &
    957. 

  957.                           *(dens/den(i,k))*acrfac
    958. 

  958.               pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
    959. 

  959.             endif
    960. 

  960. !-------------------------------------------------------------
    961. 

  961. ! psacr: Accretion of rain by snow [HL A10] [LFO 28]
    962. 

  962. !         (T<T0:R->S or R->G) (T>=T0: enhance melting of snow)
    963. 

  963. !-------------------------------------------------------------
    964. 

  964.             acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,2)            &
    965. 

  965.                     +2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,2)           &
    966. 

  966.                     +.5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,2)
    967. 

  967.             psacr(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2ave-vt2r)            &
    968. 

  968.                         *(denr/den(i,k))*acrfac
    969. 

  969.             psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
    970. 

  970.           endif
    971. 

  971. !-------------------------------------------------------------
    972. 

  972. ! pgacr: Accretion of rain by graupel [HL A12] [LFO 42]
    973. 

  973. !         (T<T0: R->G) (T>=T0: enhance melting of graupel)
    974. 

  974. !-------------------------------------------------------------
    975. 

  975.           if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
    976. 

  976.             acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,3)            &
    977. 

  977.                     +2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,3)           &
    978. 

  978.                     +.5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,3)
    979. 

  979.             pgacr(i,k) = pi**2*n0r*n0g*abs(vt2ave-vt2r)*(denr/den(i,k))        &
    980. 

  980.                         *acrfac
    981. 

  981.             pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
    982. 

  982.           endif
    983. 

  983. !
    984. 

  984. !-------------------------------------------------------------
    985. 

  985. ! pgacs: Accretion of snow by graupel [HL A13] [LFO 29]
    986. 

  986. !        (S->G): This process is eliminated in V3.0 with the 
    987. 

  987. !        new combined snow/graupel fall speeds
    988. 

  988. !-------------------------------------------------------------
    989. 

  989.           if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,2).gt.qcrmin) then
    990. 

  990.             pgacs(i,k) = 0.
    991. 

  991.           endif
    992. 

  992.           if(supcol.le.0) then
    993. 

  993.             xlf = xlf0
    994. 

  994. !-------------------------------------------------------------
    995. 

  995. ! pseml: Enhanced melting of snow by accretion of water [HL A34]
    996. 

  996. !        (T>=T0: S->R)
    997. 

  997. !-------------------------------------------------------------
    998. 

  998.             if(qrs(i,k,2).gt.0.)                                               &
    999. 

  999.               pseml(i,k) = min(max(cliq*supcol*(paacw(i,k)+psacr(i,k))         &
    1000. 

  1000.                           /xlf,-qrs(i,k,2)/dtcld),0.)
    1001. 

  1001. !-------------------------------------------------------------
    1002. 

  1002. ! pgeml: Enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22]
    1003. 

  1003. !        (T>=T0: G->R)
    1004. 

  1004. !-------------------------------------------------------------
    1005. 

  1005.             if(qrs(i,k,3).gt.0.)                                               &
    1006. 

  1006.               pgeml(i,k) = min(max(cliq*supcol*(paacw(i,k)+pgacr(i,k))         &
    1007. 

  1007.                           /xlf,-qrs(i,k,3)/dtcld),0.)
    1008. 

  1008.           endif
    1009. 

  1009.           if(supcol.gt.0) then
    1010. 

  1010. !-------------------------------------------------------------
    1011. 

  1011. ! pidep: Deposition/Sublimation rate of ice [HDC 9]
    1012. 

  1012. !       (T<T0: V->I or I->V)
    1013. 

  1013. !-------------------------------------------------------------
    1014. 

  1014.             if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
    1015. 

  1015.               pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
    1016. 

  1016.               supice = satdt-prevp(i,k)
    1017. 

  1017.               if(pidep(i,k).lt.0.) then
    1018. 

  1018.                 pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
    1019. 

  1019.                 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
    1020. 

  1020.               else
    1021. 

  1021.                 pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
    1022. 

  1022.               endif
    1023. 

  1023.               if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
    1024. 

  1024.             endif
    1025. 

  1025. !-------------------------------------------------------------
    1026. 

  1026. ! psdep: deposition/sublimation rate of snow [HDC 14]
    1027. 

  1027. !        (T<T0: V->S or S->V)
    1028. 

  1028. !-------------------------------------------------------------
    1029. 

  1029.             if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
    1030. 

  1030.               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    1031. 

  1031.               psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2)   &    
    1032. 

  1032.                            + precs2*work2(i,k)*coeres)/work1…
    1033.

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