PageRenderTime 81ms CodeModel.GetById 13ms 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 
    

  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(i,k,2)
    1033. 

  1033.               supice = satdt-prevp(i,k)-pidep(i,k)
    1034. 

  1034.               if(psdep(i,k).lt.0.) then
    1035. 

  1035.                 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
    1036. 

  1036.                 psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
    1037. 

  1037.               else
    1038. 

  1038.                 psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
    1039. 

  1039.               endif
    1040. 

  1040.               if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt))          &
    1041. 

  1041.                 ifsat = 1
    1042. 

  1042.             endif
    1043. 

  1043. !-------------------------------------------------------------
    1044. 

  1044. ! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46]
    1045. 

  1045. !        (T<T0: V->G or G->V)
    1046. 

  1046. !-------------------------------------------------------------
    1047. 

  1047.             if(qrs(i,k,3).gt.0..and.ifsat.ne.1) then
    1048. 

  1048.               coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
    1049. 

  1049.               pgdep(i,k) = (rh(i,k,2)-1.)*(precg1*rslope2(i,k,3)               &
    1050. 

  1050.                               +precg2*work2(i,k)*coeres)/work1(i,k,2)
    1051. 

  1051.               supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
    1052. 

  1052.               if(pgdep(i,k).lt.0.) then
    1053. 

  1053.                 pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
    1054. 

  1054.                 pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
    1055. 

  1055.               else
    1056. 

  1056.                 pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
    1057. 

  1057.               endif
    1058. 

  1058.               if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge.          &
    1059. 

  1059.                 abs(satdt)) ifsat = 1
    1060. 

  1060.             endif
    1061. 

  1061. !-------------------------------------------------------------
    1062. 

  1062. ! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8]
    1063. 

  1063. !       (T<T0: V->I)
    1064. 

  1064. !-------------------------------------------------------------
    1065. 

  1065.             if(supsat.gt.0.and.ifsat.ne.1) then
    1066. 

  1066.               supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
    1067. 

  1067.               xni0 = 1.e3*exp(0.1*supcol)
    1068. 

  1068.               roqi0 = 4.92e-11*xni0**1.33
    1069. 

  1069.               pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))/dtcld)
    1070. 

  1070.               pigen(i,k) = min(min(pigen(i,k),satdt),supice)
    1071. 

  1071.             endif
    1072. 

  1072. !
    1073. 

  1073. !-------------------------------------------------------------
    1074. 

  1074. ! psaut: conversion(aggregation) of ice to snow [HDC 12]
    1075. 

  1075. !        (T<T0: I->S)
    1076. 

  1076. !-------------------------------------------------------------
    1077. 

  1077.             if(qci(i,k,2).gt.0.) then
    1078. 

  1078.               qimax = roqimax/den(i,k)
    1079. 

  1079.               psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
    1080. 

  1080.             endif
    1081. 

  1081. !
    1082. 

  1082. !-------------------------------------------------------------
    1083. 

  1083. ! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37]
    1084. 

  1084. !        (T<T0: S->G)
    1085. 

  1085. !-------------------------------------------------------------
    1086. 

  1086.             if(qrs(i,k,2).gt.0.) then
    1087. 

  1087.               alpha2 = 1.e-3*exp(0.09*(-supcol))
    1088. 

  1088.               pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)),qrs(i,k,2)/dtcld)
    1089. 

  1089.             endif
    1090. 

  1090.           endif
    1091. 

  1091. !
    1092. 

  1092. !-------------------------------------------------------------
    1093. 

  1093. ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
    1094. 

  1094. !       (T>=T0: S->V)
    1095. 

  1095. !-------------------------------------------------------------
    1096. 

  1096.           if(supcol.lt.0.) then
    1097. 

  1097.             if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) then
    1098. 

  1098.               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    1099. 

  1099.               psevp(i,k) = (rh(i,k,1)-1.)*n0sfac(i,k)*(precs1                  &
    1100. 

  1100.                            *rslope2(i,k,2)+precs2*work2(i,k)                   &
    1101. 

  1101.                            *coeres)/work1(i,k,1)
    1102. 

  1102.               psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
    1103. 

  1103.             endif
    1104. 

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

  1105. ! pgevp: Evaporation of melting graupel [HL A25] [RH84 A19]
    1106. 

  1106. !       (T>=T0: G->V)
    1107. 

  1107. !-------------------------------------------------------------
    1108. 

  1108.             if(qrs(i,k,3).gt.0..and.rh(i,k,1).lt.1.) then
    1109. 

  1109.               coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
    1110. 

  1110.               pgevp(i,k) = (rh(i,k,1)-1.)*(precg1*rslope2(i,k,3)               &
    1111. 

  1111.                          +precg2*work2(i,k)*coeres)/work1(i,k,1)
    1112. 

  1112.               pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.)
    1113. 

  1113.             endif
    1114. 

  1114.           endif
    1115. 

  1115.         enddo
    1116. 

  1116.       enddo
    1117. 

  1117. !
    1118. 

  1118. !
    1119. 

  1119. !----------------------------------------------------------------
    1120. 

  1120. !     check mass conservation of generation terms and feedback to the
    1121. 

  1121. !     large scale
    1122. 

  1122. !
    1123. 

  1123.       do k = kts, kte
    1124. 

  1124.         do i = its, ite
    1125. 

  1125. !
    1126. 

  1126.           delta2=0.
    1127. 

  1127.           delta3=0.
    1128. 

  1128.           if(qrs(i,k,1).lt.1.e-4.and.qrs(i,k,2).lt.1.e-4) delta2=1.
    1129. 

  1129.           if(qrs(i,k,1).lt.1.e-4) delta3=1.
    1130. 

  1130.           if(t(i,k).le.t0c) then
    1131. 

  1131. !
    1132. 

  1132. !     cloud water
    1133. 

  1133. !
    1134. 

  1134.             value = max(qmin,qci(i,k,1))
    1135. 

  1135.             source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
    1136. 

  1136.             if (source.gt.value) then
    1137. 

  1137.               factor = value/source
    1138. 

  1138.               praut(i,k) = praut(i,k)*factor
    1139. 

  1139.               pracw(i,k) = pracw(i,k)*factor
    1140. 

  1140.               paacw(i,k) = paacw(i,k)*factor
    1141. 

  1141.             endif
    1142. 

  1142. !
    1143. 

  1143. !     cloud ice
    1144. 

  1144. !
    1145. 

  1145.             value = max(qmin,qci(i,k,2))
    1146. 

  1146.             source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k)+psaci(i,k)   &     
    1147. 

  1147.                     +pgaci(i,k))*dtcld
    1148. 

  1148.             if (source.gt.value) then
    1149. 

  1149.               factor = value/source
    1150. 

  1150.               psaut(i,k) = psaut(i,k)*factor
    1151. 

  1151.               pigen(i,k) = pigen(i,k)*factor
    1152. 

  1152.               pidep(i,k) = pidep(i,k)*factor
    1153. 

  1153.               praci(i,k) = praci(i,k)*factor
    1154. 

  1154.               psaci(i,k) = psaci(i,k)*factor
    1155. 

  1155.               pgaci(i,k) = pgaci(i,k)*factor
    1156. 

  1156.             endif
    1157. 

  1157. !
    1158. 

  1158. !     rain
    1159. 

  1159. !
    1160. 

  1160.             value = max(qmin,qrs(i,k,1))
    1161. 

  1161.             source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k)+psacr(i,k)  &    
    1162. 

  1162.                      +pgacr(i,k))*dtcld
    1163. 

  1163.             if (source.gt.value) then
    1164. 

  1164.               factor = value/source
    1165. 

  1165.               praut(i,k) = praut(i,k)*factor
    1166. 

  1166.               prevp(i,k) = prevp(i,k)*factor
    1167. 

  1167.               pracw(i,k) = pracw(i,k)*factor
    1168. 

  1168.               piacr(i,k) = piacr(i,k)*factor
    1169. 

  1169.               psacr(i,k) = psacr(i,k)*factor
    1170. 

  1170.               pgacr(i,k) = pgacr(i,k)*factor
    1171. 

  1171.             endif
    1172. 

  1172. !
    1173. 

  1173. !     snow
    1174. 

  1174. !
    1175. 

  1175.             value = max(qmin,qrs(i,k,2))
    1176. 

  1176.             source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k)+piacr(i,k)  &        
    1177. 

  1177.                      *delta3+praci(i,k)*delta3-pracs(i,k)*(1.-delta2)          &
    1178. 

  1178.                      +psacr(i,k)*delta2+psaci(i,k)-pgacs(i,k) )*dtcld
    1179. 

  1179.             if (source.gt.value) then
    1180. 

  1180.               factor = value/source
    1181. 

  1181.               psdep(i,k) = psdep(i,k)*factor
    1182. 

  1182.               psaut(i,k) = psaut(i,k)*factor
    1183. 

  1183.               pgaut(i,k) = pgaut(i,k)*factor
    1184. 

  1184.               paacw(i,k) = paacw(i,k)*factor
    1185. 

  1185.               piacr(i,k) = piacr(i,k)*factor
    1186. 

  1186.               praci(i,k) = praci(i,k)*factor
    1187. 

  1187.               psaci(i,k) = psaci(i,k)*factor
    1188. 

  1188.               pracs(i,k) = pracs(i,k)*factor
    1189. 

  1189.               psacr(i,k) = psacr(i,k)*factor
    1190. 

  1190.               pgacs(i,k) = pgacs(i,k)*factor
    1191. 

  1191.             endif
    1192. 

  1192. !
    1193. 

  1193. !     graupel
    1194. 

  1194. !
    1195. 

  1195.             value = max(qmin,qrs(i,k,3))
    1196. 

  1196.             source = -(pgdep(i,k)+pgaut(i,k)                                   &
    1197. 

  1197.                      +piacr(i,k)*(1.-delta3)+praci(i,k)*(1.-delta3)            &
    1198. 

  1198.                      +psacr(i,k)*(1.-delta2)+pracs(i,k)*(1.-delta2)            &
    1199. 

  1199.                      +pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
    1200. 

  1200.             if (source.gt.value) then
    1201. 

  1201.               factor = value/source
    1202. 

  1202.               pgdep(i,k) = pgdep(i,k)*factor
    1203. 

  1203.               pgaut(i,k) = pgaut(i,k)*factor
    1204. 

  1204.               piacr(i,k) = piacr(i,k)*factor
    1205. 

  1205.               praci(i,k) = praci(i,k)*factor
    1206. 

  1206.               psacr(i,k) = psacr(i,k)*factor
    1207. 

  1207.               pracs(i,k) = pracs(i,k)*factor
    1208. 

  1208.               paacw(i,k) = paacw(i,k)*factor
    1209. 

  1209.               pgaci(i,k) = pgaci(i,k)*factor
    1210. 

  1210.               pgacr(i,k) = pgacr(i,k)*factor
    1211. 

  1211.               pgacs(i,k) = pgacs(i,k)*factor
    1212. 

  1212.             endif
    1213. 

  1213. !
    1214. 

  1214.             work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k))
    1215. 

  1215. !     update
    1216. 

  1216.             q(i,k) = q(i,k)+work2(i,k)*dtcld
    1217. 

  1217.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
    1218. 

  1218.                            +paacw(i,k)+paacw(i,k))*dtcld,0.)
    1219. 

  1219.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
    1220. 

  1220.                            +prevp(i,k)-piacr(i,k)-pgacr(i,k)                   &
    1221. 

  1221.                            -psacr(i,k))*dtcld,0.)
    1222. 

  1222.             qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k)                 &
    1223. 

  1223.                            +psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k))       &
    1224. 

  1224.                            *dtcld,0.)
    1225. 

  1225.             qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k)      &
    1226. 

  1226.                            -pgaut(i,k)+piacr(i,k)*delta3                       &
    1227. 

  1227.                            +praci(i,k)*delta3+psaci(i,k)-pgacs(i,k)            &
    1228. 

  1228.                            -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2)          &
    1229. 

  1229.                            *dtcld,0.)
    1230. 

  1230.             qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k)                 &
    1231. 

  1231.                            +piacr(i,k)*(1.-delta3)                             &
    1232. 

  1232.                            +praci(i,k)*(1.-delta3)+psacr(i,k)*(1.-delta2)      &
    1233. 

  1233.                            +pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k)       &
    1234. 

  1234.                            +pgacr(i,k)+pgacs(i,k))*dtcld,0.)
    1235. 

  1235.             xlf = xls-xl(i,k)
    1236. 

  1236.             xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k))       &
    1237. 

  1237.                       -xl(i,k)*prevp(i,k)-xlf*(piacr(i,k)+paacw(i,k)           &
    1238. 

  1238.                       +paacw(i,k)+pgacr(i,k)+psacr(i,k))
    1239. 

  1239.             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
    1240. 

  1240.           else
    1241. 

  1241. !
    1242. 

  1242. !     cloud water
    1243. 

  1243. !
    1244. 

  1244.             value = max(qmin,qci(i,k,1))
    1245. 

  1245.             source=(praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
    1246. 

  1246.             if (source.gt.value) then
    1247. 

  1247.               factor = value/source
    1248. 

  1248.               praut(i,k) = praut(i,k)*factor
    1249. 

  1249.               pracw(i,k) = pracw(i,k)*factor
    1250. 

  1250.               paacw(i,k) = paacw(i,k)*factor
    1251. 

  1251.             endif
    1252. 

  1252. !
    1253. 

  1253. !     rain
    1254. 

  1254. !
    1255. 

  1255.             value = max(qmin,qrs(i,k,1))
    1256. 

  1256.             source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k)-pracw(i,k)  &  
    1257. 

  1257.                      -paacw(i,k)-prevp(i,k))*dtcld
    1258. 

  1258.             if (source.gt.value) then
    1259. 

  1259.               factor = value/source
    1260. 

  1260.               praut(i,k) = praut(i,k)*factor
    1261. 

  1261.               prevp(i,k) = prevp(i,k)*factor
    1262. 

  1262.               pracw(i,k) = pracw(i,k)*factor
    1263. 

  1263.               paacw(i,k) = paacw(i,k)*factor
    1264. 

  1264.               pseml(i,k) = pseml(i,k)*factor
    1265. 

  1265.               pgeml(i,k) = pgeml(i,k)*factor
    1266. 

  1266.             endif
    1267. 

  1267. !
    1268. 

  1268. !     snow
    1269. 

  1269. !
    1270. 

  1270.             value = max(qcrmin,qrs(i,k,2))
    1271. 

  1271.             source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
    1272. 

  1272.             if (source.gt.value) then
    1273. 

  1273.               factor = value/source
    1274. 

  1274.               pgacs(i,k) = pgacs(i,k)*factor
    1275. 

  1275.               psevp(i,k) = psevp(i,k)*factor
    1276. 

  1276.               pseml(i,k) = pseml(i,k)*factor
    1277. 

  1277.             endif
    1278. 

  1278. !
    1279. 

  1279. !     graupel
    1280. 

  1280. !
    1281. 

  1281.             value = max(qcrmin,qrs(i,k,3))
    1282. 

  1282.             source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
    1283. 

  1283.             if (source.gt.value) then
    1284. 

  1284.               factor = value/source
    1285. 

  1285.               pgacs(i,k) = pgacs(i,k)*factor
    1286. 

  1286.               pgevp(i,k) = pgevp(i,k)*factor
    1287. 

  1287.               pgeml(i,k) = pgeml(i,k)*factor
    1288. 

  1288.             endif
    1289. 

  1289.             work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
    1290. 

  1290. !     update
    1291. 

  1291.             q(i,k) = q(i,k)+work2(i,k)*dtcld
    1292. 

  1292.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
    1293. 

  1293.                     +paacw(i,k)+paacw(i,k))*dtcld,0.)
    1294. 

  1294.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
    1295. 

  1295.                     +prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k)               &
    1296. 

  1296.                     -pgeml(i,k))*dtcld,0.)
    1297. 

  1297.             qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k)                 &
    1298. 

  1298.                     +pseml(i,k))*dtcld,0.)
    1299. 

  1299.             qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k)                 &
    1300. 

  1300.                     +pgeml(i,k))*dtcld,0.)
    1301. 

  1301.             xlf = xls-xl(i,k)
    1302. 

  1302.             xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k))              &
    1303. 

  1303.                       -xlf*(pseml(i,k)+pgeml(i,k))
    1304. 

  1304.             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
    1305. 

  1305.           endif
    1306. 

  1306.         enddo
    1307. 

  1307.       enddo
    1308. 

  1308. !
    1309. 

  1309. ! Inline expansion for fpvs
    1310. 

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

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

  1312.       hsub = xls
    1313. 

  1313.       hvap = xlv0
    1314. 

  1314.       cvap = cpv
    1315. 

  1315.       ttp=t0c+0.01
    1316. 

  1316.       dldt=cvap-cliq
    1317. 

  1317.       xa=-dldt/rv
    1318. 

  1318.       xb=xa+hvap/(rv*ttp)
    1319. 

  1319.       dldti=cvap-cice
    1320. 

  1320.       xai=-dldti/rv
    1321. 

  1321.       xbi=xai+hsub/(rv*ttp)
    1322. 

  1322.       do k = kts, kte
    1323. 

  1323.         do i = its, ite
    1324. 

  1324.           tr=ttp/t(i,k)
    1325. 

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

  1326.           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
    1327. 

  1327.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    1328. 

  1328.           qs(i,k,1) = max(qs(i,k,1),qmin)
    1329. 

  1329.           tr=ttp/t(i,k)
    1330. 

  1330.           if(t(i,k).lt.ttp) then
    1331. 

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

  1332.           else
    1333. 

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

  1334.           endif
    1335. 

  1335.           qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
    1336. 

  1336.           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
    1337. 

  1337.           qs(i,k,2) = max(qs(i,k,2),qmin)
    1338. 

  1338.         enddo
    1339. 

  1339.       enddo
    1340. 

  1340. !
    1341. 

  1341. !----------------------------------------------------------------
    1342. 

  1342. !  pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
    1343. 

  1343. !     if there exists additional water vapor condensated/if
    1344. 

  1344. !     evaporation of cloud water is not enough to remove subsaturation
    1345. 

  1345. !
    1346. 

  1346.       do k = kts, kte
    1347. 

  1347.         do i = its, ite
    1348. 

  1348.           work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
    1349. 

  1349.           work2(i,k) = qci(i,k,1)+work1(i,k,1)
    1350. 

  1350.           pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
    1351. 

  1351.           if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.)                          &
    1352. 

  1352.             pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
    1353. 

  1353.           q(i,k) = q(i,k)-pcond(i,k)*dtcld
    1354. 

  1354.           qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
    1355. 

  1355.           t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
    1356. 

  1356.         enddo
    1357. 

  1357.       enddo
    1358. 

  1358. !
    1359. 

  1359. !
    1360. 

  1360. !----------------------------------------------------------------
    1361. 

  1361. !     padding for small values
    1362. 

  1362. !
    1363. 

  1363.       do k = kts, kte
    1364. 

  1364.         do i = its, ite
    1365. 

  1365.           if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
    1366. 

  1366.           if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
    1367. 

  1367.         enddo
    1368. 

  1368.       enddo
    1369. 

  1369.       enddo                  ! big loops
    1370. 

  1370.   END SUBROUTINE wsm62d
    1371. 

  1371. ! ...................................................................
    1372. 

  1372.       REAL FUNCTION rgmma(x)
    1373. 

  1373. !-------------------------------------------------------------------
    1374. 

  1374.   IMPLICIT NONE
    1375. 

  1375. !-------------------------------------------------------------------
    1376. 

  1376. !     rgmma function:  use infinite product form
    1377. 

  1377.       REAL :: euler
    1378. 

  1378.       PARAMETER (euler=0.577215664901532)
    1379. 

  1379.       REAL :: x, y
    1380. 

  1380.       INTEGER :: i
    1381. 

  1381.       if(x.eq.1.)then
    1382. 

  1382.         rgmma=0.
    1383. 

  1383.           else
    1384. 

  1384.         rgmma=x*exp(euler*x)
    1385. 

  1385.         do i=1,10000
    1386. 

  1386.           y=float(i)
    1387. 

  1387.           rgmma=rgmma*(1.000+x/y)*exp(-x/y)
    1388. 

  1388.         enddo
    1389. 

  1389.         rgmma=1./rgmma
    1390. 

  1390.       endif
    1391. 

  1391.       END FUNCTION rgmma
    1392. 

  1392. !
    1393. 

  1393. !--------------------------------------------------------------------------
    1394. 

  1394.       REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
    1395. 

  1395. !--------------------------------------------------------------------------
    1396. 

  1396.       IMPLICIT NONE
    1397. 

  1397. !--------------------------------------------------------------------------
    1398. 

  1398.       REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti,         &
    1399. 

  1399.            xai,xbi,ttp,tr
    1400. 

  1400.       INTEGER ice
    1401. 

  1401. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1402. 

  1402.       ttp=t0c+0.01
    1403. 

  1403.       dldt=cvap-cliq
    1404. 

  1404.       xa=-dldt/rv
    1405. 

  1405.       xb=xa+hvap/(rv*ttp)
    1406. 

  1406.       dldti=cvap-cice
    1407. 

  1407.       xai=-dldti/rv
    1408. 

  1408.       xbi=xai+hsub/(rv*ttp)
    1409. 

  1409.       tr=ttp/t
    1410. 

  1410.       if(t.lt.ttp.and.ice.eq.1) then
    1411. 

  1411.         fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
    1412. 

  1412.       else
    1413. 

  1413.         fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
    1414. 

  1414.       endif
    1415. 

  1415. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1416. 

  1416.       END FUNCTION fpvs
    1417. 

  1417. !-------------------------------------------------------------------
    1418. 

  1418.   SUBROUTINE wsm6init(den0,denr,dens,cl,cpv,allowed_to_read)
    1419. 

  1419. !-------------------------------------------------------------------
    1420. 

  1420.   IMPLICIT NONE
    1421. 

  1421. !-------------------------------------------------------------------
    1422. 

  1422. !.... constants which may not be tunable
    1423. 

  1423.    REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
    1424. 

  1424.    LOGICAL, INTENT(IN) :: allowed_to_read
    1425. 

  1425. !
    1426. 

  1426.    pi = 4.*atan(1.)
    1427. 

  1427.    xlv1 = cl-cpv
    1428. 

  1428. !
    1429. 

  1429.    qc0  = 4./3.*pi*denr*r0**3*xncr/den0  ! 0.419e-3 -- .61e-3
    1430. 

  1430.    qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
    1431. 

  1431. !
    1432. 

  1432.    bvtr1 = 1.+bvtr
    1433. 

  1433.    bvtr2 = 2.5+.5*bvtr
    1434. 

  1434.    bvtr3 = 3.+bvtr
    1435. 

  1435.    bvtr4 = 4.+bvtr
    1436. 

  1436.    bvtr6 = 6.+bvtr
    1437. 

  1437.    g1pbr = rgmma(bvtr1)
    1438. 

  1438.    g3pbr = rgmma(bvtr3)
    1439. 

  1439.    g4pbr = rgmma(bvtr4)            ! 17.837825
    1440. 

  1440.    g6pbr = rgmma(bvtr6)
    1441. 

  1441.    g5pbro2 = rgmma(bvtr2)          ! 1.8273
    1442. 

  1442.    pvtr = avtr*g4pbr/6.
    1443. 

  1443.    eacrr = 1.0
    1444. 

  1444.    pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
    1445. 

  1445.    precr1 = 2.*pi*n0r*.78
    1446. 

  1446.    precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
    1447. 

  1447.    roqimax = 2.08e22*dimax**8
    1448. 

  1448. !
    1449. 

  1449.    bvts1 = 1.+bvts
    1450. 

  1450.    bvts2 = 2.5+.5*bvts
    1451. 

  1451.    bvts3 = 3.+bvts
    1452. 

  1452.    bvts4 = 4.+bvts
    1453. 

  1453.    g1pbs = rgmma(bvts1)    !.8875
    1454. 

  1454.    g3pbs = rgmma(bvts3)
    1455. 

  1455.    g4pbs = rgmma(bvts4)    ! 12.0786
    1456. 

  1456.    g5pbso2 = rgmma(bvts2)
    1457. 

  1457.    pvts = avts*g4pbs/6.
    1458. 

  1458.    pacrs = pi*n0s*avts*g3pbs*.25
    1459. 

  1459.    precs1 = 4.*n0s*.65
    1460. 

  1460.    precs2 = 4.*n0s*.44*avts**.5*g5pbso2
    1461. 

  1461.    pidn0r =  pi*denr*n0r
    1462. 

  1462.    pidn0s =  pi*dens*n0s
    1463. 

  1463. !
    1464. 

  1464.    pacrc = pi*n0s*avts*g3pbs*.25*eacrc
    1465. 

  1465. !
    1466. 

  1466.    bvtg1 = 1.+bvtg
    1467. 

  1467.    bvtg2 = 2.5+.5*bvtg
    1468. 

  1468.    bvtg3 = 3.+bvtg
    1469. 

  1469.    bvtg4 = 4.+bvtg
    1470. 

  1470.    g1pbg = rgmma(bvtg1)
    1471. 

  1471.    g3pbg = rgmma(bvtg3)
    1472. 

  1472.    g4pbg = rgmma(bvtg4)
    1473. 

  1473.    pacrg = pi*n0g*avtg*g3pbg*.25
    1474. 

  1474.    g5pbgo2 = rgmma(bvtg2)
    1475. 

  1475.    pvtg = avtg*g4pbg/6.
    1476. 

  1476.    precg1 = 2.*pi*n0g*.78
    1477. 

  1477.    precg2 = 2.*pi*n0g*.31*avtg**.5*g5pbgo2
    1478. 

  1478.    pidn0g =  pi*deng*n0g
    1479. 

  1479. !
    1480. 

  1480.    rslopermax = 1./lamdarmax
    1481. 

  1481.    rslopesmax = 1./lamdasmax
    1482. 

  1482.    rslopegmax = 1./lamdagmax
    1483. 

  1483.    rsloperbmax = rslopermax ** bvtr
    1484. 

  1484.    rslopesbmax = rslopesmax ** bvts
    1485. 

  1485.    rslopegbmax = rslopegmax ** bvtg
    1486. 

  1486.    rsloper2max = rslopermax * rslopermax
    1487. 

  1487.    rslopes2max = rslopesmax * rslopesmax
    1488. 

  1488.    rslopeg2max = rslopegmax * rslopegmax
    1489. 

  1489.    rsloper3max = rsloper2max * rslopermax
    1490. 

  1490.    rslopes3max = rslopes2max * rslopesmax
    1491. 

  1491.    rslopeg3max = rslopeg2max * rslopegmax
    1492. 

  1492. !
    1493. 

  1493.   END SUBROUTINE wsm6init
    1494. 

  1494. !------------------------------------------------------------------------------
    1495. 

  1495.       subroutine slope_wsm6(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
    1496. 

  1496.                             vt,its,ite,kts,kte)
    1497. 

  1497.   IMPLICIT NONE
    1498. 

  1498.   INTEGER       ::               its,ite, jts,jte, kts,kte
    1499. 

  1499.   REAL, DIMENSION( its:ite , kts:kte,3) ::                                     &
    1500. 

  1500.                                                                           qrs, &
    1501. 

  1501.                                                                        rslope, &
    1502. 

  1502.                                                                       rslopeb, &                                                 
    1503. 

  1503.                                                                       rslope2, &                                                 
    1504. 

  1504.                                                                       rslope3, &                                                 
    1505. 

  1505.                                                                            vt
    1506. 

  1506.   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
    1507. 

  1507.                                                                           den, &
    1508. 

  1508.                                                                        denfac, &
    1509. 

  1509.                                                                             t
    1510. 

  1510.   REAL, PARAMETER  :: t0c = 273.15
    1511. 

  1511.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    1512. 

  1512.                                                                        n0sfac
    1513. 

  1513.   REAL       ::  lamdar, lamdas, lamdag, x, y, z, supcol
    1514. 

  1514.   integer :: i, j, k
    1515. 

  1515. !----------------------------------------------------------------
    1516. 

  1516. !     size distributions: (x=mixing ratio, y=air density):
    1517. 

  1517. !     valid for mixing ratio > 1.e-9 kg/kg.
    1518. 

  1518.       lamdar(x,y)=   sqrt(sqrt(pidn0r/(x*y)))      ! (pidn0r/(x*y))**.25
    1519. 

  1519.       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
    1520. 

  1520.       lamdag(x,y)=   sqrt(sqrt(pidn0g/(x*y)))      ! (pidn0g/(x*y))**.25
    1521. 

  1521. !
    1522. 

  1522.       do k = kts, kte
    1523. 

  1523.         do i = its, ite
    1524. 

  1524.           supcol = t0c-t(i,k)
    1525. 

  1525. !---------------------------------------------------------------
    1526. 

  1526. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
    1527. 

  1527. !---------------------------------------------------------------
    1528. 

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

  1529.           if(qrs(i,k,1).le.qcrmin)then
    1530. 

  1530.             rslope(i,k,1) = rslopermax
    1531. 

  1531.             rslopeb(i,k,1) = rsloperbmax
    1532. 

  1532.             rslope2(i,k,1) = rsloper2max
    1533. 

  1533.             rslope3(i,k,1) = rsloper3max
    1534. 

  1534.           else
    1535. 

  1535.             rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
    1536. 

  1536.             rslopeb(i,k,1) = rslope(i,k,1)**bvtr
    1537. 

  1537.             rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
    1538. 

  1538.             rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
    1539. 

  1539.           endif
    1540. 

  1540.           if(qrs(i,k,2).le.qcrmin)then
    1541. 

  1541.             rslope(i,k,2) = rslopesmax
    1542. 

  1542.             rslopeb(i,k,2) = rslopesbmax
    1543. 

  1543.             rslope2(i,k,2) = rslopes2max
    1544. 

  1544.             rslope3(i,k,2) = rslopes3max
    1545. 

  1545.           else
    1546. 

  1546.             rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
    1547. 

  1547.             rslopeb(i,k,2) = rslope(i,k,2)**bvts
    1548. 

  1548.             rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
    1549. 

  1549.             rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
    1550. 

  1550.           endif
    1551. 

  1551.           if(qrs(i,k,3).le.qcrmin)then
    1552. 

  1552.             rslope(i,k,3) = rslopegmax
    1553. 

  1553.             rslopeb(i,k,3) = rslopegbmax
    1554. 

  1554.             rslope2(i,k,3) = rslopeg2max
    1555. 

  1555.             rslope3(i,k,3) = rslopeg3max
    1556. 

  1556.           else
    1557. 

  1557.             rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
    1558. 

  1558.             rslopeb(i,k,3) = rslope(i,k,3)**bvtg
    1559. 

  1559.             rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
    1560. 

  1560.             rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
    1561. 

  1561.           endif
    1562. 

  1562.           vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
    1563. 

  1563.           vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
    1564. 

  1564.           vt(i,k,3) = pvtg*rslopeb(i,k,3)*denfac(i,k)
    1565. 

  1565.           if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
    1566. 

  1566.           if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
    1567. 

  1567.           if(qrs(i,k,3).le.0.0) vt(i,k,3) = 0.0
    1568. 

  1568.         enddo
    1569. 

  1569.       enddo
    1570. 

  1570.   END subroutine slope_wsm6
    1571. 

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

  1572.       subroutine slope_rain(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   & 
    1573. 

  1573.                             vt,its,ite,kts,kte)
    1574. 

  1574.   IMPLICIT NONE
    1575. 

  1575.   INTEGER       ::               its,ite, jts,jte, kts,kte
    1576. 

  1576.   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
    1577. 

  1577.                                                                           qrs, &
    1578. 

  1578.                                                                        rslope, &
    1579. 

  1579.                                                                       rslopeb, &
    1580. 

  1580.                                                                       rslope2, &
    1581. 

  1581.                                                                       rslope3, &
    1582. 

  1582.                                                                            vt, &      
    1583. 

  1583.                                                                           den, &
    1584. 

  1584.                                                                        denfac, &
    1585. 

  1585.                                                                             t
    1586. 

  1586.   REAL, PARAMETER  :: t0c = 273.15
    1587. 

  1587.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    1588. 

  1588.                                                                        n0sfac
    1589. 

  1589.   REAL       ::  lamdar, x, y, z, supcol
    1590. 

  1590.   integer :: i, j, k
    1591. 

  1591. !----------------------------------------------------------------
    1592. 

  1592. !     size distributions: (x=mixing ratio, y=air density):
    1593. 

  1593. !     valid for mixing ratio > 1.e-9 kg/kg.
    1594. 

  1594.       lamdar(x,y)=   sqrt(sqrt(pidn0r/(x*y)))      ! (pidn0r/(x*y))**.25
    1595. 

  1595. !
    1596. 

  1596.       do k = kts, kte
    1597. 

  1597.         do i = its, ite
    1598. 

  1598.           if(qrs(i,k).le.qcrmin)then
    1599. 

  1599.             rslope(i,k) = rslopermax
    1600. 

  1600.             rslopeb(i,k) = rsloperbmax
    1601. 

  1601.             rslope2(i,k) = rsloper2max
    1602. 

  1602.             rslope3(i,k) = rsloper3max
    1603. 

  1603.           else
    1604. 

  1604.             rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
    1605. 

  1605.             rslopeb(i,k) = rslope(i,k)**bvtr
    1606. 

  1606.             rslope2(i,k) = rslope(i,k)*rslope(i,k)
    1607. 

  1607.             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
    1608. 

  1608.           endif
    1609. 

  1609.           vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
    1610. 

  1610.           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
    1611. 

  1611.         enddo
    1612. 

  1612.       enddo
    1613. 

  1613.   END subroutine slope_rain
    1614. 

  1614. !------------------------------------------------------------------------------
    1615. 

  1615.       subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
    1616. 

  1616.                             vt,its,ite,kts,kte)
    1617. 

  1617.   IMPLICIT NONE
    1618. 

  1618.   INTEGER       ::               its,ite, jts,jte, kts,kte
    1619. 

  1619.   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
    1620. 

  1620.                                                                           qrs, &
    1621. 

  1621.                                                                        rslope, &
    1622. 

  1622.                                                                       rslopeb, &
    1623. 

  1623.                                                                       rslope2, &
    1624. 

  1624.                                                                       rslope3, &
    1625. 

  1625.                                                                            vt, &  
    1626. 

  1626.                                                                           den, &
    1627. 

  1627.                                                                        denfac, &
    1628. 

  1628.                                                                             t
    1629. 

  1629.   REAL, PARAMETER  :: t0c = 273.15
    1630. 

  1630.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    1631. 

  1631.                                                                        n0sfac
    1632. 

  1632.   REAL       ::  lamdas, x, y, z, supcol
    1633. 

  1633.   integer :: i, j, k
    1634. 

  1634. !----------------------------------------------------------------
    1635. 

  1635. !     size distributions: (x=mixing ratio, y=air density):
    1636. 

  1636. !     valid for mixing ratio > 1.e-9 kg/kg.
    1637. 

  1637.       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
    1638. 

  1638. !
    1639. 

  1639.       do k = kts, kte
    1640. 

  1640.         do i = its, ite
    1641. 

  1641.           supcol = t0c-t(i,k)
    1642. 

  1642. !---------------------------------------------------------------
    1643. 

  1643. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
    1644. 

  1644. !---------------------------------------------------------------
    1645. 

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

  1646.           if(qrs(i,k).le.qcrmin)then
    1647. 

  1647.             rslope(i,k) = rslopesmax
    1648. 

  1648.             rslopeb(i,k) = rslopesbmax
    1649. 

  1649.             rslope2(i,k) = rslopes2max
    1650. 

  1650.             rslope3(i,k) = rslopes3max
    1651. 

  1651.           else
    1652. 

  1652.             rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
    1653. 

  1653.             rslopeb(i,k) = rslope(i,k)**bvts
    1654. 

  1654.             rslope2(i,k) = rslope(i,k)*rslope(i,k)
    1655. 

  1655.             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
    1656. 

  1656.           endif
    1657. 

  1657.           vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
    1658. 

  1658.           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
    1659. 

  1659.         enddo
    1660. 

  1660.       enddo
    1661. 

  1661.   END subroutine slope_snow
    1662. 

  1662. !----------------------------------------------------------------------------------
    1663. 

  1663.       subroutine slope_graup(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
    1664. 

  1664.                             vt,its,ite,kts,kte)
    1665. 

  1665.   IMPLICIT NONE
    1666. 

  1666.   INTEGER       ::               its,ite, jts,jte, kts,kte
    1667. 

  1667.   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
    1668. 

  1668.                                                                           qrs, &
    1669. 

  1669.                                                                        rslope, &
    1670. 

  1670.                                                                       rslopeb, &
    1671. 

  1671.                                                                       rslope2, &
    1672. 

  1672.                                                                       rslope3, &
    1673. 

  1673.                                                                            vt, &  
    1674. 

  1674.                                                                           den, &
    1675. 

  1675.                                                                        denfac, &
    1676. 

  1676.                                                                             t
    1677. 

  1677.   REAL, PARAMETER  :: t0c = 273.15
    1678. 

  1678.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    1679. 

  1679.                                                                        n0sfac
    1680. 

  1680.   REAL       ::  lamdag, x, y, z, supcol
    1681. 

  1681.   integer :: i, j, k
    1682. 

  1682. !----------------------------------------------------------------
    1683. 

  1683. !     size distributions: (x=mixing ratio, y=air density):
    1684. 

  1684. !     valid for mixing ratio > 1.e-9 kg/kg.
    1685. 

  1685.       lamdag(x,y)=   sqrt(sqrt(pidn0g/(x*y)))      ! (pidn0g/(x*y))**.25
    1686. 

  1686. !
    1687. 

  1687.       do k = kts, kte
    1688. 

  1688.         do i = its, ite
    1689. 

  1689. !---------------------------------------------------------------
    1690. 

  1690. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
    1691. 

  1691. !---------------------------------------------------------------
    1692. 

  1692.           if(qrs(i,k).le.qcrmin)then
    1693. 

  1693.             rslope(i,k) = rslopegmax
    1694. 

  1694.             rslopeb(i,k) = rslopegbmax
    1695. 

  1695.             rslope2(i,k) = rslopeg2max
    1696. 

  1696.             rslope3(i,k) = rslopeg3max
    1697. 

  1697.           else
    1698. 

  1698.             rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k))
    1699. 

  1699.             rslopeb(i,k) = rslope(i,k)**bvtg
    1700. 

  1700.             rslope2(i,k) = rslope(i,k)*rslope(i,k)
    1701. 

  1701.             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
    1702. 

  1702.           endif
    1703. 

  1703.           vt(i,k) = pvtg*rslopeb(i,k)*denfac(i,k)
    1704. 

  1704.           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
    1705. 

  1705.         enddo
    1706. 

  1706.       enddo
    1707. 

  1707.   END subroutine slope_graup
    1708. 

  1708. !---------------------------------------------------------------------------------
    1709. 

  1709. !-------------------------------------------------------------------
    1710. 

  1710.       SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
    1711. 

  1711. !-------------------------------------------------------------------
    1712. 

  1712. !
    1713. 

  1713. ! for non-iteration semi-Lagrangain forward advection for cloud
    1714. 

  1714. ! with mass conservation and positive definite advection
    1715. 

  1715. ! 2nd order interpolation with monotonic piecewise linear method
    1716. 

  1716. ! this routine is under assumption of decfl < 1 for semi_Lagrangian
    1717. 

  1717. !
    1718. 

  1718. ! dzl    depth of model layer in meter
    1719. 

  1719. ! wwl    terminal velocity at model layer m/s
    1720. 

  1720. ! rql    cloud density*mixing ration
    1721. 

  1721. ! precip precipitation
    1722. 

  1722. ! dt     time step
    1723. 

  1723. ! id     kind of precip: 0 test case; 1 raindrop
    1724. 

  1724. ! iter   how many time to guess mean terminal velocity: 0 pure forward.
    1725. 

  1725. !        0 : use departure wind for advection
    1726. 

  1726. !        1 : use mean wind for advection
    1727. 

  1727. !        > 1 : use mean wind after iter-1 iterations
    1728. 

  1728. !
    1729. 

  1729. ! author: hann-ming henry juang <henry.juang@noaa.gov>
    1730. 

  1730. !         implemented by song-you hong
    1731. 

  1731. !
    1732. 

  1732.       implicit none
    1733. 

  1733.       integer  im,km,id
    1734. 

  1734.       real  dt
    1735. 

  1735.       real  dzl(im,km),wwl(im,km),rql(im,km),precip(im)
    1736. 

  1736.       real  denl(im,km),denfacl(im,km),tkl(im,km)
    1737. 

  1737. !
    1738. 

  1738.       integer  i,k,n,m,kk,kb,kt,iter
    1739. 

  1739.       real  tl,tl2,qql,dql,qqd
    1740. 

  1740.       real  th,th2,qqh,dqh
    1741. 

  1741.       real  zsum,qsum,dim,dip,c1,con1,fa1,fa2
    1742. 

  1742.       real  allold, allnew, zz, dzamin, cflmax, decfl
    1743. 

  1743.       real  dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
    1744. 

  1744.       real  den(km), denfac(km), tk(km)
    1745. 

  1745.       real  wi(km+1), zi(km+1), za(km+1)
    1746. 

  1746.       real  qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
    1747. 

  1747.       real  dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
    1748. 

  1748. !
    1749. 

  1749.       precip(:) = 0.0
    1750. 

  1750. !
    1751. 

  1751.       i_loop : do i=1,im
    1752. 

  1752. ! -----------------------------------
    1753. 

  1753.       dz(:) = dzl(i,:)
    1754. 

  1754.       qq(:) = rql(i,:)
    1755. 

  1755.       ww(:) = wwl(i,:)
    1756. 

  1756.       den(:) = denl(i,:)
    1757. 

  1757.       denfac(:) = denfacl(i,:)
    1758. 

  1758.       tk(:) = tkl(i,:)
    1759. 

  1759. ! skip for no precipitation for all layers
    1760. 

  1760.       allold = 0.0
    1761. 

  1761.       do k=1,km
    1762. 

  1762.         allold = allold + qq(k)
    1763. 

  1763.       enddo
    1764. 

  1764.       if(allold.le.0.0) then
    1765. 

  1765.         cycle i_loop
    1766. 

  1766.       endif
    1767. 

  1767. !
    1768. 

  1768. ! compute interface values
    1769. 

  1769.       zi(1)=0.0
    1770. 

  1770.       do k=1,km
    1771. 

  1771.         zi(k+1) = zi(k)+dz(k)
    1772. 

  1772.       enddo
    1773. 

  1773. !
    1774. 

  1774. ! save departure wind
    1775. 

  1775.       wd(:) = ww(:)
    1776. 

  1776.       n=1
    1777. 

  1777.  100  continue
    1778. 

  1778. ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
    1779. 

  1779. ! 2nd order interpolation to get wi
    1780. 

  1780.       wi(1) = ww(1)
    1781. 

  1781.       wi(km+1) = ww(km)
    1782. 

  1782.       do k=2,km
    1783. 

  1783.         wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
    1784. 

  1784.       enddo
    1785. 

  1785. ! 3rd order interpolation to get wi
    1786. 

  1786.       fa1 = 9./16.
    1787. 

  1787.       fa2 = 1./16.
    1788. 

  1788.       wi(1) = ww(1)
    1789. 

  1789.       wi(2) = 0.5*(ww(2)+ww(1))
    1790. 

  1790.       do k=3,km-1
    1791. 

  1791.         wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
    1792. 

  1792.       enddo
    1793. 

  1793.       wi(km) = 0.5*(ww(km)+ww(km-1))
    1794. 

  1794.       wi(km+1) = ww(km)
    1795. 

  1795. !
    1796. 

  1796. ! terminate of top of raingroup
    1797. 

  1797.       do k=2,km
    1798. 

  1798.         if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
    1799. 

  1799.       enddo
    1800. 

  1800. !
    1801. 

  1801. ! diffusivity of wi
    1802. 

  1802.       con1 = 0.05
    1803. 

  1803.       do k=km,1,-1
    1804. 

  1804.         decfl = (wi(k+1)-wi(k))*dt/dz(k)
    1805. 

  1805.         if( decfl .gt. con1 ) then
    1806. 

  1806.           wi(k) = wi(k+1) - con1*dz(k)/dt
    1807. 

  1807.         endif
    1808. 

  1808.       enddo
    1809. 

  1809. ! compute arrival point
    1810. 

  1810.       do k=1,km+1
    1811. 

  1811.         za(k) = zi(k) - wi(k)*dt
    1812. 

  1812.       enddo
    1813. 

  1813. !
    1814. 

  1814.       do k=1,km
    1815. 

  1815.         dza(k) = za(k+1)-za(k)
    1816. 

  1816.       enddo
    1817. 

  1817.       dza(km+1) = zi(km+1) - za(km+1)
    1818. 

  1818. !
    1819. 

  1819. ! computer deformation at arrival point
    1820. 

  1820.       do k=1,km
    1821. 

  1821.         qa(k) = qq(k)*dz(k)/dza(k)
    1822. 

  1822.         qr(k) = qa(k)/den(k)
    1823. 

  1823.       enddo
    1824. 

  1824.       qa(km+1) = 0.0
    1825. 

  1825. !     call maxmin(km,1,qa,' arrival points ')
    1826. 

  1826. !
    1827. 

  1827. ! compute arrival terminal velocity, and estimate mean terminal velocity
    1828. 

  1828. ! then back to use mean terminal velocity
    1829. 

  1829.       if( n.le.iter ) then
    1830. 

  1830.         call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
    1831. 

  1831.         if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
    1832. 

  1832.         do k=1,km
    1833. 

  1833. !#ifdef DEBUG
    1834. 

  1834. !        print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
    1835. 

  1835. !#endif
    1836. 

  1836. ! mean wind is average of departure and new arrival winds
    1837. 

  1837.           ww(k) = 0.5* ( wd(k)+wa(k) )
    1838. 

  1838.         enddo
    1839. 

  1839.         was(:) = wa(:)
    1840. 

  1840.         n=n+1
    1841. 

  1841.         go to 100
    1842. 

  1842.       endif
    1843. 

  1843. !
    1844. 

  1844. ! estimate values at arrival cell interface with monotone
    1845. 

  1845.       do k=2,km
    1846. 

  1846.         dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
    1847. 

  1847.         dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
    1848. 

  1848.         if( dip*dim.le.0.0 ) then
    1849. 

  1849.           qmi(k)=qa(k)
    1850. 

  1850.           qpi(k)=qa(k)
    1851. 

  1851.         else
    1852. 

  1852.           qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
    1853. 

  1853.           qmi(k)=2.0*qa(k)-qpi(k)
    1854. 

  1854.           if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
    1855. 

  1855.             qpi(k) = qa(k)
    1856. 

  1856.             qmi(k) = qa(k)
    1857. 

  1857.           endif
    1858. 

  1858.         endif
    1859. 

  1859.       enddo
    1860. 

  1860.       qpi(1)=qa(1)
    1861. 

  1861.       qmi(1)=qa(1)
    1862. 

  1862.       qmi(km+1)=qa(km+1)
    1863. 

  1863.       qpi(km+1)=qa(km+1)
    1864. 

  1864. !
    1865. 

  1865. ! interpolation to regular point
    1866. 

  1866.       qn = 0.0
    1867. 

  1867.       kb=1
    1868. 

  1868.       kt=1
    1869. 

  1869.       intp : do k=1,km
    1870. 

  1870.              kb=max(kb-1,1)
    1871. 

  1871.              kt=max(kt-1,1)
    1872. 

  1872. ! find kb and kt
    1873. 

  1873.              if( zi(k).ge.za(km+1) ) then
    1874. 

  1874.                exit intp
    1875. 

  1875.              else
    1876. 

  1876.                find_kb : do kk=kb,km
    1877. 

  1877.                          if( zi(k).le.za(kk+1) ) then
    1878. 

  1878.                            kb = kk
    1879. 

  1879.                            exit find_kb
    1880. 

  1880.                          else
    1881. 

  1881.                            cycle find_kb
    1882. 

  1882.                          endif
    1883. 

  1883.                enddo find_kb
    1884. 

  1884.                find_kt : do kk=kt,km
    1885. 

  1885.                          if( zi(k+1).le.za(kk) ) then
    1886. 

  1886.                            kt = kk
    1887. 

  1887.                            exit find_kt
    1888. 

  1888.                          else
    1889. 

  1889.                            cycle find_kt
    1890. 

  1890.                          endif
    1891. 

  1891.                enddo find_kt
    1892. 

  1892.                kt = kt - 1
    1893. 

  1893. ! compute q with piecewise constant method
    1894. 

  1894.                if( kt.eq.kb ) then
    1895. 

  1895.                  tl=(zi(k)-za(kb))/dza(kb)
    1896. 

  1896.                  th=(zi(k+1)-za(kb))/dza(kb)
    1897. 

  1897.                  tl2=tl*tl
    1898. 

  1898.                  th2=th*th
    1899. 

  1899.                  qqd=0.5*(qpi(kb)-qmi(kb))
    1900. 

  1900.                  qqh=qqd*th2+qmi(kb)*th
    1901. 

  1901.                  qql=qqd*tl2+qmi(kb)*tl
    1902. 

  1902.                  qn(k) = (qqh-qql)/(th-tl)
    1903. 

  1903.                else if( kt.gt.kb ) then
    1904. 

  1904.                  tl=(zi(k)-za(kb))/dza(kb)
    1905. 

  1905.                  tl2=tl*tl
    1906. 

  1906.                  qqd=0.5*(qpi(kb)-qmi(kb))
    1907. 

  1907.                  qql=qqd*tl2+qmi(kb)*tl
    1908. 

  1908.                  dql = qa(kb)-qql
    1909. 

  1909.                  zsum  = (1.-tl)*dza(kb)
    1910. 

  1910.                  qsum  = dql*dza(kb)
    1911. 

  1911.                  if( kt-kb.gt.1 ) then
    1912. 

  1912.                  do m=kb+1,kt-1
    1913. 

  1913.                    zsum = zsum + dza(m)
    1914. 

  1914.                    qsum = qsum + qa(m) * dza(m)
    1915. 

  1915.                  enddo
    1916. 

  1916.                  endif
    1917. 

  1917.                  th=(zi(k+1)-za(kt))/dza(kt)
    1918. 

  1918.                  th2=th*th
    1919. 

  1919.                  qqd=0.5*(qpi(kt)-qmi(kt))
    1920. 

  1920.                  dqh=qqd*th2+qmi(kt)*th
    1921. 

  1921.                  zsum  = zsum + th*dza(kt)
    1922. 

  1922.                  qsum  = qsum + dqh*dza(kt)
    1923. 

  1923.                  qn(k) = qsum/zsum
    1924. 

  1924.                endif
    1925. 

  1925.                cycle intp
    1926. 

  1926.              endif
    1927. 

  1927. !
    1928. 

  1928.        enddo intp
    1929. 

  1929. !
    1930. 

  1930. ! rain out
    1931. 

  1931.       sum_precip: do k=1,km
    1932. 

  1932.                     if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
    1933. 

  1933.                       precip(i) = precip(i) + qa(k)*dza(k)
    1934. 

  1934.                       cycle sum_precip
    1935. 

  1935.                     else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
    1936. 

  1936.                       precip(i) = precip(i) + qa(k)*(0.0-za(k))
    1937. 

  1937.                       exit sum_precip
    1938. 

  1938.                     endif
    1939. 

  1939.                     exit sum_precip
    1940. 

  1940.       enddo sum_precip
    1941. 

  1941. !
    1942. 

  1942. ! replace the new values
    1943. 

  1943.       rql(i,:) = qn(:)
    1944. 

  1944. !
    1945. 

  1945. ! ----------------------------------
    1946. 

  1946.       enddo i_loop
    1947. 

  1947. !
    1948. 

  1948.   END SUBROUTINE nislfv_rain_plm
    1949. 

  1949. !-------------------------------------------------------------------
    1950. 

  1950.       SUBROUTINE nislfv_rain_plm6(im,km,denl,denfacl,tkl,dzl,wwl,rql,rql2, precip1, precip2,dt,id,iter)
    1951. 

  1951. !-------------------------------------------------------------------
    1952. 

  1952. !
    1953. 

  1953. ! for non-iteration semi-Lagrangain forward advection for cloud
    1954. 

  1954. ! with mass conservation and positive definite advection
    1955. 

  1955. ! 2nd order interpolation with monotonic piecewise linear method
    1956. 

  1956. ! this routine is under assumption of decfl < 1 for semi_Lagrangian
    1957. 

  1957. !
    1958. 

  1958. ! dzl    depth of model layer in meter
    1959. 

  1959. ! wwl    terminal velocity at model layer m/s
    1960. 

  1960. ! rql    cloud density*mixing ration
    1961. 

  1961. ! precip precipitation
    1962. 

  1962. ! dt     time step
    1963. 

  1963. ! id     kind of precip: 0 test case; 1 raindrop
    1964. 

  1964. ! iter   how many time to guess mean terminal velocity: 0 pure forward.
    1965. 

  1965. !        0 : use departure wind for advection
    1966. 

  1966. !        1 : use mean wind for advection
    1967. 

  1967. !        > 1 : use mean wind after iter-1 iterations
    1968. 

  1968. !
    1969. 

  1969. ! author: hann-ming henry juang <henry.juang@noaa.gov>
    1970. 

  1970. !         implemented by song-you hong
    1971. 

  1971. !
    1972. 

  1972.       implicit none
    1973. 

  1973.       integer  im,km,id
    1974. 

  1974.       real  dt
    1975. 

  1975.       real  dzl(im,km),wwl(im,km),rql(im,km),rql2(im,km),precip(im),precip1(im),precip2(im)
    1976. 

  1976.       real  denl(im,km),denfacl(im,km),tkl(im,km)
    1977. 

  1977. !
    1978. 

  1978.       integer  i,k,n,m,kk,kb,kt,iter,ist
    1979. 

  1979.       real  tl,tl2,qql,dql,qqd
    1980. 

  1980.       real  th,th2,qqh,dqh
    1981. 

  1981.       real  zsum,qsum,dim,dip,c1,con1,fa1,fa2
    1982. 

  1982.       real  allold, allnew, zz, dzamin, cflmax, decfl
    1983. 

  1983.       real  dz(km), ww(km), qq(km), qq2(km), wd(km), wa(km), wa2(km), was(km)
    1984. 

  1984.       real  den(km), denfac(km), tk(km)
    1985. 

  1985.       real  wi(km+1), zi(km+1), za(km+1)
    1986. 

  1986.       real  qn(km), qr(km),qr2(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
    1987. 

  1987.       real  dza(km+1), qa(km+1), qa2(km+1),qmi(km+1), qpi(km+1)
    1988. 

  1988. !
    1989. 

  1989.       precip(:) = 0.0
    1990. 

  1990.       precip1(:) = 0.0
    1991. 

  1991.       precip2(:) = 0.0
    1992. 

  1992. !
    1993. 

  1993.       i_loop : do i=1,im
    1994. 

  1994. ! -----------------------------------
    1995. 

  1995.       dz(:) = dzl(i,:)
    1996. 

  1996.       qq(:) = rql(i,:)
    1997. 

  1997.       qq2(:) = rql2(i,:)
    1998. 

  1998.       ww(:) = wwl(i,:)
    1999. 

  1999.       den(:) = denl(i,:)
    2000. 

  2000.       denfac(:) = denfacl(i,:)
    2001. 

  2001.       tk(:) = tkl(i,:)
    2002. 

  2002. ! skip for no precipitation for all layers
    2003. 

  2003.       allold = 0.0
    2004. 

  2004.       do k=1,km
    2005. 

  2005.         allold = allold + qq(k)
    2006. 

  2006.       enddo
    2007. 

  2007.       if(allold.le.0.0) then
    2008. 

  2008.         cycle i_loop
    2009. 

  2009.       endif
    2010. 

  2010. !
    2011. 

  2011. ! compute interface values
    2012. 

  2012.       zi(1)=0.0
    2013. 

  2013.       do k=1,km
    2014. 

  2014.         zi(k+1) = zi(k)+dz(k)
    2015. 

  2015.       enddo
    2016. 

  2016. !
    2017. 

  2017. ! save departure wind
    2018. 

  2018.       wd(:) = ww(:)
    2019. 

  2019.       n=1
    2020. 

  2020.  100  continue
    2021. 

  2021. ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
    2022. 

  2022. ! 2nd order interpolation to get wi
    2023. 

  2023.       wi(1) = ww(1)
    2024. 

  2024.       wi(km+1) = ww(km)
    2025. 

  2025.       do k=2,km
    2026. 

  2026.         wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
    2027. 

  2027.       enddo
    2028. 

  2028. ! 3rd order interpolation to get wi
    2029. 

  2029.       fa1 = 9./16.
    2030. 

  2030.       fa2 = 1./16.
    2031. 

  2031.       wi(1) = ww(1)
    2032. 

  2032.       wi(2) = 0.5*(ww(2)+ww(1))
    2033. 

  2033.       do k=3,km-1
    2034. 

  2034.         wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
    2035. 

  2035.       enddo
    2036. 

  2036.       wi(km) = 0.5*(ww(km)+ww(km-1))
    2037. 

  2037.       wi(km+1) = ww(km)
    2038. 

  2038. !
    2039. 

  2039. ! terminate of top of raingroup
    2040. 

  2040.       do k=2,km
    2041. 

  2041.         if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
    2042. 

  2042.       enddo
    2043. 

  2043. !
    2044. 

  2044. ! diffusivity of wi
    2045. 

  2045.       con1 = 0.05
    2046. 

  2046.       do k=km,1,-1
    2047. 

  2047.         decfl = (wi(k+1)-wi(k))*dt/dz(k)
    2048. 

  2048.         if( decfl .gt. con1 ) then
    2049. 

  2049.           wi(k) = wi(k+1) - con1*dz(k)/dt
    2050. 

  2050.         endif
    2051. 

  2051.       enddo
    2052. 

  2052. ! compute arrival point
    2053. 

  2053.       do k=1,km+1
    2054. 

  2054.         za(k) = zi(k) - wi(k)*dt
    2055. 

  2055.       enddo
    2056. 

  2056. !
    2057. 

  2057.       do k=1,km
    2058. 

  2058.         dza(k) = za(k+1)-za(k)
    2059. 

  2059.       enddo
    2060. 

  2060.       dza(km+1) = zi(km+1) - za(km+1)
    2061. 

  2061. !
    2062. 

  2062. ! computer deformation at arrival point
    2063. 

  2063.       do k=1,km
    2064. 

  2064.         qa(k) = qq(k)*dz(k)/dza(k)
    2065. 

  2065.         qa2(k) = qq2(k)*dz(k)/dza(k)
    2066. 

  2066.         qr(k) = qa(k)/den(k)
    2067. 

  2067.         qr2(k) = qa2(k)/den(k)
    2068. 

  2068.       enddo
    2069. 

  2069.       qa(km+1) = 0.0
    2070. 

  2070.       qa2(km+1) = 0.0
    2071. 

  2071. !     call maxmin(km,1,qa,' arrival points ')
    2072. 

  2072. !
    2073. 

  2073. ! compute arrival terminal velocity, and estimate mean terminal velocity
    2074. 

  2074. ! then back to use mean terminal velocity
    2075. 

  2075.       if( n.le.iter ) then
    2076. 

  2076.         call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
    2077. 

  2077.         call slope_graup(qr2,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa2,1,1,1,km)
    2078. 

  2078.         do k = 1, km
    2079. 

  2079.           tmp(k) = max((qr(k)+qr2(k)), 1.E-15)
    2080. 

  2080.           IF ( tmp(k) .gt. 1.e-15 ) THEN
    2081. 

  2081.             wa(k) = (wa(k)*qr(k) + wa2(k)*qr2(k))/tmp(k)
    2082. 

  2082.           ELSE
    2083. 

  2083.             wa(k) = 0.
    2084. 

  2084.           ENDIF
    2085. 

  2085.         enddo
    2086. 

  2086.         if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
    2087. 

  2087.         do k=1,km
    2088. 

  2088. !#ifdef DEBUG
    2089. 

  2089. !        print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k), &
    2090. 

  2090. !           ww(k),wa(k)
    2091. 

  2091. !#endif
    2092. 

  2092. ! mean wind is average of departure and new arrival winds
    2093. 

  2093.           ww(k) = 0.5* ( wd(k)+wa(k) )
    2094. 

  2094.         enddo
    2095. 

  2095.         was(:) = wa(:)
    2096. 

  2096.         n=n+1
    2097. 

  2097.         go to 100
    2098. 

  2098.       endif
    2099. 

  2099.       ist_loop : do ist = 1, 2
    2100. 

  2100.       if (ist.eq.2) then
    2101. 

  2101.        qa(:) = qa2(:)
    2102. 

  2102.       endif
    2103. 

  2103. !
    2104. 

  2104.       precip(i) = 0.
    2105. 

  2105. !
    2106. 

  2106. ! estimate values at arrival cell interface with monotone
    2107. 

  2107.       do k=2,km
    2108. 

  2108.         dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
    2109. 

  2109.         dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
    2110. 

  2110.         if( dip*dim.le.0.0 ) then
    2111. 

  2111.           qmi(k)=qa(k)
    2112. 

  2112.           qpi(k)=qa(k)
    2113. 

  2113.         else
    2114. 

  2114.           qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
    2115. 

  2115.           qmi(k)=2.0*qa(k)-qpi(k)
    2116. 

  2116.           if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
    2117. 

  2117.             qpi(k) = qa(k)
    2118. 

  2118.             qmi(k) = qa(k)
    2119. 

  2119.           endif
    2120. 

  2120.         endif
    2121. 

  2121.       enddo
    2122. 

  2122.       qpi(1)=qa(1)
    2123. 

  2123.       qmi(1)=qa(1)
    2124. 

  2124.       qmi(km+1)=qa(km+1)
    2125. 

  2125.       qpi(km+1)=qa(km+1)
    2126. 

  2126. !
    2127. 

  2127. ! interpolation to regular point
    2128. 

  2128.       qn = 0.0
    2129. 

  2129.       kb=1
    2130. 

  2130.       kt=1
    2131. 

  2131.       intp : do k=1,km
    2132. 

  2132.              kb=max(kb-1,1)
    2133. 

  2133.              kt=max(kt-1,1)
    2134. 

  2134. ! find kb and kt
    2135. 

  2135.              if( zi(k).ge.za(km+1) ) then
    2136. 

  2136.                exit intp
    2137. 

  2137.              else
    2138. 

  2138.                find_kb : do kk=kb,km
    2139. 

  2139.                          if( zi(k).le.za(kk+1) ) then
    2140. 

  2140.                            kb = kk
    2141. 

  2141.                            exit find_kb
    2142. 

  2142.                          else
    2143. 

  2143.                            cycle find_kb
    2144. 

  2144.                          endif
    2145. 

  2145.                enddo find_kb
    2146. 

  2146.                find_kt : do kk=kt,km
    2147. 

  2147.                          if( zi(k+1).le.za(kk) ) then
    2148. 

  2148.                            kt = kk
    2149. 

  2149.                            exit find_kt
    2150. 

  2150.                          else
    2151. 

  2151.                            cycle find_kt
    2152. 

  2152.                          endif
    2153. 

  2153.                enddo find_kt
    2154. 

  2154.                kt = kt - 1
    2155. 

  2155. ! compute q with piecewise constant method
    2156. 

  2156.                if( kt.eq.kb ) then
    2157. 

  2157.                  tl=(zi(k)-za(kb))/dza(kb)
    2158. 

  2158.                  th=(zi(k+1)-za(kb))/dza(kb)
    2159. 

  2159.                  tl2=tl*tl
    2160. 

  2160.                  th2=th*th
    2161. 

  2161.                  qqd=0.5*(qpi(kb)-qmi(kb))
    2162. 

  2162.                  qqh=qqd*th2+qmi(kb)*th
    2163. 

  2163.                  qql=qqd*tl2+qmi(kb)*tl
    2164. 

  2164.                  qn(k) = (qqh-qql)/(th-tl)
    2165. 

  2165.                else if( kt.gt.kb ) then
    2166. 

  2166.                  tl=(zi(k)-za(kb))/dza(kb)
    2167. 

  2167.                  tl2=tl*tl
    2168. 

  2168.                  qqd=0.5*(qpi(kb)-qmi(kb))
    2169. 

  2169.                  qql=qqd*tl2+qmi(kb)*tl
    2170. 

  2170.                  dql = qa(kb)-qql
    2171. 

  2171.                  zsum  = (1.-tl)*dza(kb)
    2172. 

  2172.                  qsum  = dql*dza(kb)
    2173. 

  2173.                  if( kt-kb.gt.1 ) then
    2174. 

  2174.                  do m=kb+1,kt-1
    2175. 

  2175.                    zsum = zsum + dza(m)
    2176. 

  2176.                    qsum = qsum + qa(m) * dza(m)
    2177. 

  2177.                  enddo
    2178. 

  2178.                  endif
    2179. 

  2179.                  th=(zi(k+1)-za(kt))/dza(kt)
    2180. 

  2180.                  th2=th*th
    2181. 

  2181.                  qqd=0.5*(qpi(kt)-qmi(kt))
    2182. 

  2182.                  dqh=qqd*th2+qmi(kt)*th
    2183. 

  2183.                  zsum  = zsum + th*dza(kt)
    2184. 

  2184.                  qsum  = qsum + dqh*dza(kt)
    2185. 

  2185.                  qn(k) = qsum/zsum
    2186. 

  2186.                endif
    2187. 

  2187.                cycle intp
    2188. 

  2188.              endif
    2189. 

  2189. !
    2190. 

  2190.        enddo intp
    2191. 

  2191. !
    2192. 

  2192. ! rain out
    2193. 

  2193.       sum_precip: do k=1,km
    2194. 

  2194.                     if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
    2195. 

  2195.                       precip(i) = precip(i) + qa(k)*dza(k)
    2196. 

  2196.                       cycle sum_precip
    2197. 

  2197.                     else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
    2198. 

  2198.                       precip(i) = precip(i) + qa(k)*(0.0-za(k))
    2199. 

  2199.                       exit sum_precip
    2200. 

  2200.                     endif
    2201. 

  2201.                     exit sum_precip
    2202. 

  2202.       enddo sum_precip
    2203. 

  2203. !
    2204. 

  2204. ! replace the new values
    2205. 

  2205.       if(ist.eq.1) then
    2206. 

  2206.         rql(i,:) = qn(:)
    2207. 

  2207.         precip1(i) = precip(i)
    2208. 

  2208.       else
    2209. 

  2209.         rql2(i,:) = qn(:)
    2210. 

  2210.         precip2(i) = precip(i)
    2211. 

  2211.       endif
    2212. 

  2212.       enddo ist_loop
    2213. 

  2213. !
    2214. 

  2214. ! ----------------------------------
    2215. 

  2215.       enddo i_loop
    2216. 

  2216. !
    2217. 

  2217.   END SUBROUTINE nislfv_rain_plm6
    2218. 

  2218. END MODULE module_mp_wsm6
    2219. 

  2219.