PageRenderTime 81ms CodeModel.GetById 18ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_mp_wdm6.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2657 lines | 1935 code | 1 blank | 721 comment | 156 complexity | 487eded20e93a0325d996a56a299f8ab 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_wdm6
    10. 

  10. !
    11. 

  11. !
    12. 

  12. !
    13. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  28.    REAL, PARAMETER, PRIVATE :: lamdacmax = 1.e10  ! limited maximum value for slope parameter of cloud water 
    29. 

  29.    REAL, PARAMETER, PRIVATE :: lamdarmax = 1.e8   ! limited maximum value for slope parameter of rain
    30. 

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

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

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

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

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

  35.    REAL, PARAMETER, PRIVATE :: alpha = .12        ! .122 exponen factor for n0s
    36. 

  36.    REAL, PARAMETER, PRIVATE :: pfrz1 = 100.       ! constant in Biggs freezing
    37. 

  37.    REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66       ! constant in Biggs freezing
    38. 

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

  39.    REAL, PARAMETER, PRIVATE :: ncmin = 1.e1       ! minimum value for Nc 
    40. 

  40.    REAL, PARAMETER, PRIVATE :: nrmin = 1.e-2      ! minimum value for Nr
    41. 

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

  42.    REAL, PARAMETER, PRIVATE :: dens  =  100.0     ! Density of snow
    43. 

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

  44. !
    45. 

  45.    REAL, PARAMETER, PRIVATE :: satmax = 1.0048    ! maximum saturation value for CCN activation 
    46. 

  46.                                                   ! 1.008 for maritime /1.0048 for conti 
    47. 

  47.    REAL, PARAMETER, PRIVATE :: actk = 0.6         ! parameter for the CCN activation
    48. 

  48.    REAL, PARAMETER, PRIVATE :: actr = 1.5         ! radius of activated CCN drops
    49. 

  49.    REAL, PARAMETER, PRIVATE :: ncrk1 = 3.03e3     ! Long's collection kernel coefficient
    50. 

  50.    REAL, PARAMETER, PRIVATE :: ncrk2 = 2.59e15    ! Long's collection kernel coefficient
    51. 

  51.    REAL, PARAMETER, PRIVATE :: di100 = 1.e-4      ! parameter related with accretion and collection of cloud drops
    52. 

  52.    REAL, PARAMETER, PRIVATE :: di600 = 6.e-4      ! parameter related with accretion and collection of cloud drops
    53. 

  53.    REAL, PARAMETER, PRIVATE :: di2000 = 2000.e-6  ! parameter related with accretion and collection of cloud drops 
    54. 

  54.    REAL, PARAMETER, PRIVATE :: di82    = 82.e-6   ! dimater related with raindrops evaporation
    55. 

  55.    REAL, PARAMETER, PRIVATE :: di15    = 15.e-6   ! auto conversion takes place beyond this diameter 
    56. 

  56. !
    57. 

  57.    REAL, SAVE ::                                           &
    58. 

  58.              qc0,qck1,pidnc,bvtr1,bvtr2,bvtr3,bvtr4,bvtr5, &
    59. 

  59.              bvtr6,bvtr7, bvtr2o5,bvtr3o5,                 &
    60. 

  60.              g1pbr,g2pbr,g3pbr,g4pbr,g5pbr,g6pbr,g7pbr,    &
    61. 

  61.              g5pbro2,g7pbro2,pi,                           &
    62. 

  62.              pvtr,pvtrn,eacrr,pacrr,pidn0r,pidnr,          &
    63. 

  63.              precr1,precr2,xmmax,roqimax,bvts1,bvts2,      &
    64. 

  64.              bvts3,bvts4,g1pbs,g3pbs,g4pbs,g5pbso2,        &
    65. 

  65.              pvts,pacrs,precs1,precs2,pidn0s,xlv1,pacrc,   &
    66. 

  66.              bvtg1,bvtg2,bvtg3,bvtg4,g1pbg,g3pbg,g4pbg,    &
    67. 

  67.              g5pbgo2,pvtg,pacrg,precg1,precg2,pidn0g,      &
    68. 

  68.              rslopecmax,rslopec2max,rslopec3max,           &
    69. 

  69.              rslopermax,rslopesmax,rslopegmax,             &
    70. 

  70.              rsloperbmax,rslopesbmax,rslopegbmax,          &
    71. 

  71.              rsloper2max,rslopes2max,rslopeg2max,          &
    72. 

  72.              rsloper3max,rslopes3max,rslopeg3max
    73. 

  73. CONTAINS
    74. 

  74. !===================================================================
    75. 

  75. !
    76. 

  76.   SUBROUTINE wdm6(th, q, qc, qr, qi, qs, qg,               &
    77. 

  77.                     nn, nc, nr,                            &
    78. 

  78.                     den, pii, p, delz,                     &
    79. 

  79.                     delt,g, cpd, cpv, ccn0, rd, rv, t0c,   &
    80. 

  80.                     ep1, ep2, qmin,                        &
    81. 

  81.                     XLS, XLV0, XLF0, den0, denr,           &
    82. 

  82.                     cliq,cice,psat,                        &
    83. 

  83.                     rain, rainncv,                         &
    84. 

  84.                     snow, snowncv,                         &
    85. 

  85.                     sr,                                    &
    86. 

  86.                     graupel, graupelncv,                   &
    87. 

  87.                     itimestep,                             &
    88. 

  88.                     ids,ide, jds,jde, kds,kde,             &
    89. 

  89.                     ims,ime, jms,jme, kms,kme,             &
    90. 

  90.                     its,ite, jts,jte, kts,kte              &
    91. 

  91.                                                            )
    92. 

  92. !-------------------------------------------------------------------
    93. 

  93.   IMPLICIT NONE
    94. 

  94. !-------------------------------------------------------------------
    95. 

  95. !
    96. 

  96. !  This code is a WRF double-moment 6-class GRAUPEL phase 
    97. 

  97. !  microphyiscs scheme (WDM6). The WDM microphysics scheme predicts 
    98. 

  98. !  number concentrations for warm rain species including clouds and 
    99. 

  99. !  rain. cloud condensation nuclei (CCN) is also predicted.
    100. 

  100. !  The cold rain species including ice, snow, graupel follow the 
    101. 

  101. !  WRF single-moment 6-class microphysics (WSM6, Hong and Lim 2006)
    102. 

  102. !  in which theoretical background for WSM ice phase microphysics is 
    103. 

  103. !  based on Hong et al. (2004). A new mixed-phase terminal velocity
    104. 

  104. !  for precipitating ice is introduced in WSM6 (Dudhia et al. 2008). 
    105. 

  105. !  The WDM scheme is described in Lim and Hong (2010).
    106. 

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

  107. !
    108. 

  108. !  WDM6 cloud scheme
    109. 

  109. !
    110. 

  110. !  Coded by Kyo-Sun Lim and Song-You Hong (Yonsei Univ.) Fall 2008
    111. 

  111. !
    112. 

  112. !  Implemented by Kyo-Sun Lim and Jimy Dudhia (NCAR) Winter 2008
    113. 

  113. !
    114. 

  114. !  Reference) Lim and Hong (LH, 2010) Mon. Wea. Rev. 
    115. 

  115. !             Juang and Hong (JH, 2010) Mon. Wea. Rev.
    116. 

  116. !             Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev. 
    117. 

  117. !             Hong and Lim (HL, 2006) J. Korean Meteor. Soc. 
    118. 

  118. !             Cohard and Pinty (CP, 2000) Quart. J. Roy. Meteor. Soc.
    119. 

  119. !             Khairoutdinov and Kogan (KK, 2000) Mon. Wea. Rev.
    120. 

  120. !             Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
    121. 

  121. !             
    122. 

  122. !             Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
    123. 

  123. !             Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
    124. 

  124. !             Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
    125. 

  125. !
    126. 

  126.   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
    127. 

  127.                                       ims,ime, jms,jme, kms,kme , &
    128. 

  128.                                       its,ite, jts,jte, kts,kte
    129. 

  129.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    130. 

  130.         INTENT(INOUT) ::                                          &
    131. 

  131.                                                               th, &
    132. 

  132.                                                                q, &
    133. 

  133.                                                               qc, &
    134. 

  134.                                                               qi, &
    135. 

  135.                                                               qr, &
    136. 

  136.                                                               qs, &
    137. 

  137.                                                               qg, &
    138. 

  138.                                                               nn, & 
    139. 

  139.                                                               nc, &
    140. 

  140.                                                               nr
    141. 

  141.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    142. 

  142.         INTENT(IN   ) ::                                          &
    143. 

  143.                                                              den, &
    144. 

  144.                                                              pii, &
    145. 

  145.                                                                p, &
    146. 

  146.                                                             delz
    147. 

  147.   REAL, INTENT(IN   ) ::                                    delt, &
    148. 

  148.                                                                g, &
    149. 

  149.                                                               rd, &
    150. 

  150.                                                               rv, &
    151. 

  151.                                                              t0c, &
    152. 

  152.                                                             den0, &
    153. 

  153.                                                              cpd, &
    154. 

  154.                                                              cpv, &
    155. 

  155.                                                             ccn0, &
    156. 

  156.                                                              ep1, &
    157. 

  157.                                                              ep2, &
    158. 

  158.                                                             qmin, &
    159. 

  159.                                                              XLS, &
    160. 

  160.                                                             XLV0, &
    161. 

  161.                                                             XLF0, &
    162. 

  162.                                                             cliq, &
    163. 

  163.                                                             cice, &
    164. 

  164.                                                             psat, &
    165. 

  165.                                                             denr
    166. 

  166.   INTEGER, INTENT(IN   ) ::                            itimestep
    167. 

  167.   REAL, DIMENSION( ims:ime , jms:jme ),                           &
    168. 

  168.         INTENT(INOUT) ::                                    rain, &
    169. 

  169.                                                          rainncv, &
    170. 

  170.                                                               sr
    171. 

  171.   REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                 &
    172. 

  172.         INTENT(INOUT) ::                                    snow, &
    173. 

  173.                                                          snowncv
    174. 

  174.   REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                 &
    175. 

  175.         INTENT(INOUT) ::                                 graupel, &
    176. 

  176.                                                         graupelncv
    177. 

  177. ! LOCAL VAR
    178. 

  178.   REAL, DIMENSION( its:ite , kts:kte ) ::   t
    179. 

  179.   REAL, DIMENSION( its:ite , kts:kte, 2 ) ::   qci
    180. 

  180.   REAL, DIMENSION( its:ite , kts:kte, 3 ) ::   qrs, ncr
    181. 

  181.   INTEGER ::               i,j,k
    182. 

  182. !-------------------------------------------------------------------
    183. 

  183.       IF (itimestep .eq. 1) THEN 
    184. 

  184.         DO j=jms,jme
    185. 

  185.            DO k=kms,kme    
    186. 

  186.            DO i=ims,ime
    187. 

  187.               nn(i,k,j) = ccn0   
    188. 

  188.            ENDDO
    189. 

  189.            ENDDO
    190. 

  190.         ENDDO
    191. 

  191.       ENDIF
    192. 

  192. !
    193. 

  193.       DO j=jts,jte
    194. 

  194.          DO k=kts,kte
    195. 

  195.          DO i=its,ite
    196. 

  196.             t(i,k)=th(i,k,j)*pii(i,k,j)
    197. 

  197.             qci(i,k,1) = qc(i,k,j)
    198. 

  198.             qci(i,k,2) = qi(i,k,j)
    199. 

  199.             qrs(i,k,1) = qr(i,k,j)
    200. 

  200.             qrs(i,k,2) = qs(i,k,j)
    201. 

  201.             qrs(i,k,3) = qg(i,k,j)
    202. 

  202.             ncr(i,k,1) = nn(i,k,j)
    203. 

  203.             ncr(i,k,2) = nc(i,k,j)
    204. 

  204.             ncr(i,k,3) = nr(i,k,j)     
    205. 

  205.          ENDDO
    206. 

  206.          ENDDO
    207. 

  207.          !  Sending array starting locations of optional variables may cause
    208. 

  208.          !  troubles, so we explicitly change the call.
    209. 

  209.          CALL wdm62D(t, q(ims,kms,j), qci, qrs, ncr               &
    210. 

  210.                     ,den(ims,kms,j)                               &
    211. 

  211.                     ,p(ims,kms,j), delz(ims,kms,j)                &
    212. 

  212.                     ,delt,g, cpd, cpv, ccn0, rd, rv, t0c          &
    213. 

  213.                     ,ep1, ep2, qmin                               &
    214. 

  214.                     ,XLS, XLV0, XLF0, den0, denr                  &
    215. 

  215.                     ,cliq,cice,psat                               &
    216. 

  216.                     ,j                                            &
    217. 

  217.                     ,rain(ims,j),rainncv(ims,j)                   &
    218. 

  218.                     ,sr(ims,j)                                    &
    219. 

  219.                     ,ids,ide, jds,jde, kds,kde                    &
    220. 

  220.                     ,ims,ime, jms,jme, kms,kme                    &
    221. 

  221.                     ,its,ite, jts,jte, kts,kte                    &
    222. 

  222.                     ,snow(ims,j),snowncv(ims,j)                   &
    223. 

  223.                     ,graupel(ims,j),graupelncv(ims,j)             & 
    224. 

  224.                                                                    )
    225. 

  225.          DO K=kts,kte
    226. 

  226.          DO I=its,ite
    227. 

  227.             th(i,k,j)=t(i,k)/pii(i,k,j)
    228. 

  228.             qc(i,k,j) = qci(i,k,1)
    229. 

  229.             qi(i,k,j) = qci(i,k,2)
    230. 

  230.             qr(i,k,j) = qrs(i,k,1)
    231. 

  231.             qs(i,k,j) = qrs(i,k,2)
    232. 

  232.             qg(i,k,j) = qrs(i,k,3)
    233. 

  233.             nn(i,k,j) = ncr(i,k,1)
    234. 

  234.             nc(i,k,j) = ncr(i,k,2)
    235. 

  235.             nr(i,k,j) = ncr(i,k,3)   
    236. 

  236.          ENDDO
    237. 

  237.          ENDDO
    238. 

  238.       ENDDO
    239. 

  239.   END SUBROUTINE wdm6
    240. 

  240. !===================================================================
    241. 

  241. !
    242. 

  242.   SUBROUTINE wdm62D(t, q, qci, qrs, ncr, den, p, delz             &
    243. 

  243.                    ,delt,g, cpd, cpv, ccn0, rd, rv, t0c           &
    244. 

  244.                    ,ep1, ep2, qmin                                &
    245. 

  245.                    ,XLS, XLV0, XLF0, den0, denr                   &
    246. 

  246.                    ,cliq,cice,psat                                &
    247. 

  247.                    ,lat                                           &
    248. 

  248.                    ,rain,rainncv                                  &
    249. 

  249.                    ,sr                                            &
    250. 

  250.                    ,ids,ide, jds,jde, kds,kde                     &
    251. 

  251.                    ,ims,ime, jms,jme, kms,kme                     &
    252. 

  252.                    ,its,ite, jts,jte, kts,kte                     &
    253. 

  253.                    ,snow,snowncv                                  &
    254. 

  254.                    ,graupel,graupelncv                            &
    255. 

  255.                                                                   )
    256. 

  256. !-------------------------------------------------------------------
    257. 

  257.   IMPLICIT NONE
    258. 

  258. !-------------------------------------------------------------------
    259. 

  259.   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
    260. 

  260.                                       ims,ime, jms,jme, kms,kme , &
    261. 

  261.                                       its,ite, jts,jte, kts,kte , &
    262. 

  262.                                       lat
    263. 

  263.   REAL, DIMENSION( its:ite , kts:kte ),                           &
    264. 

  264.         INTENT(INOUT) ::                                          &
    265. 

  265.                                                                t
    266. 

  266.   REAL, DIMENSION( its:ite , kts:kte, 2 ),                        &
    267. 

  267.         INTENT(INOUT) ::                                          &
    268. 

  268.                                                              qci
    269. 

  269.   REAL, DIMENSION( its:ite , kts:kte, 3 ),                        &
    270. 

  270.         INTENT(INOUT) ::                                          &
    271. 

  271.                                                              qrs, &
    272. 

  272.                                                              ncr 
    273. 

  273.   REAL, DIMENSION( ims:ime , kms:kme ),                           &
    274. 

  274.         INTENT(INOUT) ::                                          &
    275. 

  275.                                                                q
    276. 

  276.   REAL, DIMENSION( ims:ime , kms:kme ),                           &
    277. 

  277.         INTENT(IN   ) ::                                          &
    278. 

  278.                                                              den, &
    279. 

  279.                                                                p, &
    280. 

  280.                                                             delz
    281. 

  281.   REAL, INTENT(IN   ) ::                                    delt, &
    282. 

  282.                                                                g, &
    283. 

  283.                                                              cpd, &
    284. 

  284.                                                              cpv, &
    285. 

  285.                                                             ccn0, &
    286. 

  286.                                                              t0c, &
    287. 

  287.                                                             den0, &
    288. 

  288.                                                               rd, &
    289. 

  289.                                                               rv, &
    290. 

  290.                                                              ep1, &
    291. 

  291.                                                              ep2, &
    292. 

  292.                                                             qmin, &
    293. 

  293.                                                              XLS, &
    294. 

  294.                                                             XLV0, &
    295. 

  295.                                                             XLF0, &
    296. 

  296.                                                             cliq, &
    297. 

  297.                                                             cice, &
    298. 

  298.                                                             psat, &
    299. 

  299.                                                             denr
    300. 

  300.   REAL, DIMENSION( ims:ime ),                                     &
    301. 

  301.         INTENT(INOUT) ::                                    rain, &
    302. 

  302.                                                          rainncv, &
    303. 

  303.                                                               sr
    304. 

  304.   REAL, DIMENSION( ims:ime ),     OPTIONAL,                       &
    305. 

  305.         INTENT(INOUT) ::                                    snow, &
    306. 

  306.                                                          snowncv
    307. 

  307.   REAL, DIMENSION( ims:ime ),     OPTIONAL,                       &
    308. 

  308.         INTENT(INOUT) ::                                 graupel, &
    309. 

  309.                                                       graupelncv
    310. 

  310. ! LOCAL VAR
    311. 

  311.   REAL, DIMENSION( its:ite , kts:kte , 3) ::                      &
    312. 

  312.         rh, qs, rslope, rslope2, rslope3, rslopeb,                &
    313. 

  313.         falk, fall, work1, qrs_tmp
    314. 

  314.   REAL, DIMENSION( its:ite , kts:kte ) ::                         & 
    315. 

  315.         rslopec, rslopec2,rslopec3 
    316. 

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

  317.         avedia 
    318. 

  318.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    319. 

  319.         workn,falln,falkn
    320. 

  320.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    321. 

  321.         worka,workr
    322. 

  322.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    323. 

  323.         den_tmp, delz_tmp, ncr_tmp 
    324. 

  324.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    325. 

  325.         falkc, work1c, work2c, fallc
    326. 

  326.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    327. 

  327.         pcact, prevp, psdep, pgdep, praut, psaut, pgaut,          &
    328. 

  328.         pracw, psacw, pgacw, pgacr, pgacs, psaci, pgmlt, praci,   &
    329. 

  329.         piacr, pracs, psacr, pgaci, pseml, pgeml
    330. 

  330.   REAL, DIMENSION( its:ite , kts:kte ) :: paacw
    331. 

  331.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    332. 

  332.         nraut, nracw, ncevp, nccol, nrcol,                        &
    333. 

  333.         nsacw, ngacw, niacr, nsacr, ngacr, naacw,                 &
    334. 

  334.         nseml, ngeml, ncact 
    335. 

  335.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    336. 

  336.         pigen, pidep, pcond, xl, cpm, work2, psmlt, psevp,        &
    337. 

  337.         denfac, xni, pgevp,n0sfac, qsum,                          &
    338. 

  338.         denqrs1, denqr1, denqrs2, denqrs3, denncr3, denqci  
    339. 

  339.   REAL, DIMENSION( its:ite ) ::                                   &
    340. 

  340.         delqrs1, delqrs2, delqrs3, delncr3, delqi
    341. 

  341.   REAL, DIMENSION( its:ite ) :: tstepsnow, tstepgraup
    342. 

  342.   REAL :: gfac, sfac
    343. 

  343. ! variables for optimization
    344. 

  344.   REAL, DIMENSION( its:ite )           :: tvec1
    345. 

  345.   REAL :: temp
    346. 

  346.   INTEGER, DIMENSION( its:ite ) :: mnstep, numndt
    347. 

  347.   INTEGER, DIMENSION( its:ite ) :: mstep, numdt
    348. 

  348.   LOGICAL, DIMENSION( its:ite ) :: flgcld
    349. 

  349.   REAL  ::                                                        &
    350. 

  350.             cpmcal, xlcal, lamdac,                                &
    351. 

  351.             diffus,                                               &
    352. 

  352.             viscos, xka, venfac, conden, diffac,                  &
    353. 

  353.             x, y, z, a, b, c, d, e,                               &
    354. 

  354.             ndt, qdt, holdrr, holdrs, holdrg, supcol, supcolt,    &
    355. 

  355.             pvt, coeres, supsat, dtcld, xmi, eacrs, satdt,        &
    356. 

  356.             qimax, diameter, xni0, roqi0,                         &
    357. 

  357.             fallsum, fallsum_qsi, fallsum_qg,                     &
    358. 

  358.             vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi,                   &
    359. 

  359.             xlwork2, factor, source, value, coecol,               &
    360. 

  360.             nfrzdtr, nfrzdtc,                                     &
    361. 

  361.             taucon, lencon, lenconcr,                       &
    362. 

  362.             xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3 
    363. 

  363.   REAL  :: vt2ave
    364. 

  364.   REAL  :: holdc, holdci
    365. 

  365. !
    366. 

  366.   INTEGER :: i, j, k, mstepmax,                                                &
    367. 

  367.             iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
    368. 

  368. ! Temporaries used for inlining fpvs function
    369. 

  369.   REAL  :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
    370. 

  370. !
    371. 

  371. !=================================================================
    372. 

  372. !   compute internal functions
    373. 

  373. !
    374. 

  374.       cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
    375. 

  375.       xlcal(x) = xlv0-xlv1*(x-t0c)
    376. 

  376. !----------------------------------------------------------------
    377. 

  377. !     size distributions: (x=mixing ratio, y=air density):
    378. 

  378. !     valid for mixing ratio > 1.e-9 kg/kg.
    379. 

  379. !
    380. 

  380. ! Optimizatin : A**B => exp(log(A)*(B)) 
    381. 

  381.       lamdac(x,y,z)= exp(log(((pidnc*z)/(x*y)))*((.33333333)))
    382. 

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

  383. !     diffus: diffusion coefficient of the water vapor
    384. 

  384. !     viscos: kinematic viscosity(m2s-1)
    385. 

  385. !
    386. 

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

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

  388.       xka(x,y) = 1.414e3*viscos(x,y)*y
    389. 

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

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

  391.                      /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
    392. 

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

  393. !
    394. 

  394.       idim = ite-its+1
    395. 

  395.       kdim = kte-kts+1
    396. 

  396. !
    397. 

  397. !----------------------------------------------------------------
    398. 

  398. !     paddint 0 for negative values generated by dynamics
    399. 

  399. !
    400. 

  400.       do k = kts, kte
    401. 

  401.         do i = its, ite
    402. 

  402.           qci(i,k,1) = max(qci(i,k,1),0.0)
    403. 

  403.           qrs(i,k,1) = max(qrs(i,k,1),0.0)
    404. 

  404.           qci(i,k,2) = max(qci(i,k,2),0.0)
    405. 

  405.           qrs(i,k,2) = max(qrs(i,k,2),0.0)
    406. 

  406.           qrs(i,k,3) = max(qrs(i,k,3),0.0)
    407. 

  407.           ncr(i,k,1) = max(ncr(i,k,1),0.0)
    408. 

  408.           ncr(i,k,2) = max(ncr(i,k,2),0.0)
    409. 

  409.           ncr(i,k,3) = max(ncr(i,k,3),0.0) 
    410. 

  410.         enddo
    411. 

  411.       enddo
    412. 

  412. !
    413. 

  413. !----------------------------------------------------------------
    414. 

  414. !     latent heat for phase changes and heat capacity. neglect the
    415. 

  415. !     changes during microphysical process calculation
    416. 

  416. !     emanuel(1994)
    417. 

  417. !
    418. 

  418.       do k = kts, kte
    419. 

  419.         do i = its, ite
    420. 

  420.           cpm(i,k) = cpmcal(q(i,k))
    421. 

  421.           xl(i,k) = xlcal(t(i,k))
    422. 

  422.         enddo
    423. 

  423.       enddo
    424. 

  424.       do k = kts, kte
    425. 

  425.         do i = its, ite
    426. 

  426.           delz_tmp(i,k) = delz(i,k)
    427. 

  427.           den_tmp(i,k) = den(i,k)
    428. 

  428.         enddo
    429. 

  429.       enddo
    430. 

  430. !
    431. 

  431. !----------------------------------------------------------------
    432. 

  432. !    initialize the surface rain, snow, graupel
    433. 

  433. !
    434. 

  434.       do i = its, ite
    435. 

  435.         rainncv(i) = 0.
    436. 

  436.         if(PRESENT (snowncv) .AND. PRESENT (snow)) snowncv(i) = 0.
    437. 

  437.         if(PRESENT (graupelncv) .AND. PRESENT (graupel)) graupelncv(i) = 0.
    438. 

  438.         sr(i) = 0.
    439. 

  439. ! new local array to catch step snow and graupel
    440. 

  440.         tstepsnow(i) = 0.
    441. 

  441.         tstepgraup(i) = 0.
    442. 

  442.       enddo
    443. 

  443. !
    444. 

  444. !----------------------------------------------------------------
    445. 

  445. !     compute the minor time steps.
    446. 

  446. !
    447. 

  447.       loops = max(nint(delt/dtcldcr),1)
    448. 

  448.       dtcld = delt/loops
    449. 

  449.       if(delt.le.dtcldcr) dtcld = delt
    450. 

  450. !
    451. 

  451.       do loop = 1,loops
    452. 

  452. !
    453. 

  453. !----------------------------------------------------------------
    454. 

  454. !     initialize the large scale variables
    455. 

  455. !
    456. 

  456.       do i = its, ite
    457. 

  457.         mstep(i) = 1
    458. 

  458.         mnstep(i) = 1
    459. 

  459.         flgcld(i) = .true.
    460. 

  460.       enddo
    461. 

  461. !
    462. 

  462.       do k = kts, kte
    463. 

  463.         CALL VREC( tvec1(its), den(its,k), ite-its+1)
    464. 

  464.         do i = its, ite
    465. 

  465.           tvec1(i) = tvec1(i)*den0
    466. 

  466.         enddo
    467. 

  467.         CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
    468. 

  468.       enddo
    469. 

  469. !
    470. 

  470. ! Inline expansion for fpvs
    471. 

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

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

  473.       hsub = xls
    474. 

  474.       hvap = xlv0
    475. 

  475.       cvap = cpv
    476. 

  476.       ttp=t0c+0.01
    477. 

  477.       dldt=cvap-cliq
    478. 

  478.       xa=-dldt/rv
    479. 

  479.       xb=xa+hvap/(rv*ttp)
    480. 

  480.       dldti=cvap-cice
    481. 

  481.       xai=-dldti/rv
    482. 

  482.       xbi=xai+hsub/(rv*ttp)
    483. 

  483.       do k = kts, kte
    484. 

  484.         do i = its, ite
    485. 

  485.           tr=ttp/t(i,k)
    486. 

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

  487.           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
    488. 

  488.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    489. 

  489.           qs(i,k,1) = max(qs(i,k,1),qmin)
    490. 

  490.           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
    491. 

  491.           tr=ttp/t(i,k)
    492. 

  492.           if(t(i,k).lt.ttp) then
    493. 

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

  494.           else
    495. 

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

  496.           endif
    497. 

  497.           qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
    498. 

  498.           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
    499. 

  499.           qs(i,k,2) = max(qs(i,k,2),qmin)
    500. 

  500.           rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
    501. 

  501.         enddo
    502. 

  502.       enddo
    503. 

  503. !
    504. 

  504. !----------------------------------------------------------------
    505. 

  505. !     initialize the variables for microphysical physics
    506. 

  506. !
    507. 

  507. !
    508. 

  508.       do k = kts, kte
    509. 

  509.         do i = its, ite
    510. 

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

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

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

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

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

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

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

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

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

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

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

  521.           pracs(i,k) = 0.
    522. 

  522.           psacr(i,k) = 0.
    523. 

  523.           pgacw(i,k) = 0.
    524. 

  524.           paacw(i,k) = 0.
    525. 

  525.           pgaci(i,k) = 0.
    526. 

  526.           pgacr(i,k) = 0.
    527. 

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

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

  529.           pidep(i,k) = 0.
    530. 

  530.           pcond(i,k) = 0.
    531. 

  531.           psmlt(i,k) = 0.
    532. 

  532.           pgmlt(i,k) = 0.
    533. 

  533.           pseml(i,k) = 0.
    534. 

  534.           pgeml(i,k) = 0.
    535. 

  535.           psevp(i,k) = 0.
    536. 

  536.           pgevp(i,k) = 0.
    537. 

  537.           pcact(i,k) = 0.
    538. 

  538.           falk(i,k,1) = 0.
    539. 

  539.           falk(i,k,2) = 0.
    540. 

  540.           falk(i,k,3) = 0.
    541. 

  541.           fall(i,k,1) = 0.
    542. 

  542.           fall(i,k,2) = 0.
    543. 

  543.           fall(i,k,3) = 0.
    544. 

  544.           fallc(i,k) = 0.
    545. 

  545.           falkc(i,k) = 0.
    546. 

  546.           falln(i,k) =0.
    547. 

  547.           falkn(i,k) =0.
    548. 

  548.           xni(i,k) = 1.e3
    549. 

  549.           nsacw(i,k) = 0.
    550. 

  550.           ngacw(i,k) = 0.
    551. 

  551.           naacw(i,k) = 0.
    552. 

  552.           niacr(i,k) = 0.
    553. 

  553.           nsacr(i,k) = 0.
    554. 

  554.           ngacr(i,k) = 0.
    555. 

  555.           nseml(i,k) = 0.
    556. 

  556.           ngeml(i,k) = 0.
    557. 

  557.           nracw(i,k) = 0.
    558. 

  558.           nccol(i,k) = 0.
    559. 

  559.           nrcol(i,k) = 0.
    560. 

  560.           ncact(i,k) = 0.
    561. 

  561.           nraut(i,k) = 0.
    562. 

  562.           ncevp(i,k) = 0.
    563. 

  563.         enddo
    564. 

  564.       enddo
    565. 

  565.       do k = kts, kte
    566. 

  566.         do i = its, ite
    567. 

  567.           if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin ) then
    568. 

  568.             rslopec(i,k) = rslopecmax
    569. 

  569.             rslopec2(i,k) = rslopec2max
    570. 

  570.             rslopec3(i,k) = rslopec3max
    571. 

  571.           else
    572. 

  572.             rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
    573. 

  573.             rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
    574. 

  574.             rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
    575. 

  575.           endif
    576. 

  576. !-------------------------------------------------------------
    577. 

  577. ! Ni: ice crystal number concentraiton   [HDC 5c]
    578. 

  578. !-------------------------------------------------------------
    579. 

  579.           temp = (den(i,k)*max(qci(i,k,2),qmin))
    580. 

  580.           temp = sqrt(sqrt(temp*temp*temp))
    581. 

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

  582.         enddo
    583. 

  583.       enddo
    584. 

  584. !----------------------------------------------------------------
    585. 

  585. !     compute the fallout term:
    586. 

  586. !     first, vertical terminal velosity for minor loops
    587. 

  587. !----------------------------------------------------------------
    588. 

  588.       do k = kts, kte
    589. 

  589.         do i = its, ite
    590. 

  590.           qrs_tmp(i,k,1) = qrs(i,k,1)
    591. 

  591.           qrs_tmp(i,k,2) = qrs(i,k,2)
    592. 

  592.           qrs_tmp(i,k,3) = qrs(i,k,3)
    593. 

  593.           ncr_tmp(i,k) = ncr(i,k,3)
    594. 

  594.         enddo
    595. 

  595.       enddo
    596. 

  596.       call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, & 
    597. 

  597.                      rslope3,work1,workn,its,ite,kts,kte)
    598. 

  598. !
    599. 

  599. ! vt update for qr and nr
    600. 

  600.       mstepmax = 1
    601. 

  601.       numdt = 1
    602. 

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

  603.         do i = its, ite
    604. 

  604.           work1(i,k,1) = work1(i,k,1)/delz(i,k)
    605. 

  605.           workn(i,k) = workn(i,k)/delz(i,k)
    606. 

  606.           numdt(i) = max(nint(max(work1(i,k,1),workn(i,k))*dtcld+.5),1)
    607. 

  607.           if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
    608. 

  608.         enddo
    609. 

  609.       enddo
    610. 

  610.       do i = its, ite
    611. 

  611.         if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
    612. 

  612.       enddo
    613. 

  613. !
    614. 

  614.       do n = 1, mstepmax
    615. 

  615.         k = kte
    616. 

  616.         do i = its, ite
    617. 

  617.           if(n.le.mstep(i)) then
    618. 

  618.             falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
    619. 

  619.             falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
    620. 

  620.             fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
    621. 

  621.             falln(i,k) = falln(i,k)+falkn(i,k)
    622. 

  622.             qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
    623. 

  623.             ncr(i,k,3) = max(ncr(i,k,3)-falkn(i,k)*dtcld,0.)
    624. 

  624.           endif
    625. 

  625.         enddo
    626. 

  626.         do k = kte-1, kts, -1
    627. 

  627.           do i = its, ite
    628. 

  628.             if(n.le.mstep(i)) then
    629. 

  629.               falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
    630. 

  630.               falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
    631. 

  631.               fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
    632. 

  632.               falln(i,k) = falln(i,k)+falkn(i,k)
    633. 

  633.               qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1)           &
    634. 

  634.                           *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
    635. 

  635.               ncr(i,k,3) = max(ncr(i,k,3)-(falkn(i,k)-falkn(i,k+1)*delz(i,k+1) &
    636. 

  636.                           /delz(i,k))*dtcld,0.)
    637. 

  637.             endif
    638. 

  638.           enddo
    639. 

  639.         enddo
    640. 

  640.         do k = kts, kte
    641. 

  641.           do i = its, ite
    642. 

  642.             qrs_tmp(i,k,1) = qrs(i,k,1)
    643. 

  643.             ncr_tmp(i,k) = ncr(i,k,3)
    644. 

  644.           enddo
    645. 

  645.         enddo
    646. 

  646.         call slope_rain(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
    647. 

  647.                      rslope3,work1,workn,its,ite,kts,kte)
    648. 

  648.         do k = kte, kts, -1
    649. 

  649.           do i = its, ite
    650. 

  650.             work1(i,k,1) = work1(i,k,1)/delz(i,k)
    651. 

  651.             workn(i,k) = workn(i,k)/delz(i,k)
    652. 

  652.           enddo
    653. 

  653.         enddo
    654. 

  654.       enddo
    655. 

  655. ! for semi
    656. 

  656.       do k = kte, kts, -1
    657. 

  657.         do i = its, ite
    658. 

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

  659.           if(qsum(i,k) .gt. 1.e-15 ) then
    660. 

  660.             worka(i,k) = (work1(i,k,2)*qrs(i,k,2) + work1(i,k,3)*qrs(i,k,3)) &
    661. 

  661.                       /qsum(i,k)
    662. 

  662.           else
    663. 

  663.             worka(i,k) = 0.
    664. 

  664.           endif
    665. 

  665.           denqrs2(i,k) = den(i,k)*qrs(i,k,2)
    666. 

  666.           denqrs3(i,k) = den(i,k)*qrs(i,k,3)
    667. 

  667.         enddo
    668. 

  668.       enddo
    669. 

  669.       call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka,         &
    670. 

  670.                            denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1)
    671. 

  671.       do k = kts, kte
    672. 

  672.         do i = its, ite
    673. 

  673.           qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
    674. 

  674.           qrs(i,k,3) = max(denqrs3(i,k)/den(i,k),0.)
    675. 

  675.           fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k)
    676. 

  676.           fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k)
    677. 

  677.         enddo
    678. 

  678.       enddo
    679. 

  679.       do i = its, ite
    680. 

  680.         fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
    681. 

  681.         fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld
    682. 

  682.       enddo
    683. 

  683.       do k = kts, kte
    684. 

  684.         do i = its, ite
    685. 

  685.           qrs_tmp(i,k,1) = qrs(i,k,1)
    686. 

  686.           qrs_tmp(i,k,2) = qrs(i,k,2)
    687. 

  687.           qrs_tmp(i,k,3) = qrs(i,k,3)
    688. 

  688.           ncr_tmp(i,k) = ncr(i,k,3)
    689. 

  689.         enddo
    690. 

  690.       enddo
    691. 

  691.       call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
    692. 

  692.                      rslope3,work1,workn,its,ite,kts,kte)
    693. 

  693. !
    694. 

  694.       do k = kte, kts, -1
    695. 

  695.         do i = its, ite
    696. 

  696.           supcol = t0c-t(i,k)
    697. 

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

  698.           if(t(i,k).gt.t0c) then
    699. 

  699. !---------------------------------------------------------------
    700. 

  700. ! psmlt: melting of snow [HL A33] [RH83 A25]
    701. 

  701. !       (T>T0: QS->QR)
    702. 

  702. !---------------------------------------------------------------
    703. 

  703.             xlf = xlf0
    704. 

  704.             work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
    705. 

  705.             if(qrs(i,k,2).gt.0.) then
    706. 

  706.               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    707. 

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

  708.                          *n0sfac(i,k)*(precs1*rslope2(i,k,2)                 &
    709. 

  709.                          +precs2*work2(i,k)*coeres)
    710. 

  710.               psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),-qrs(i,k,2)     &
    711. 

  711.                          /mstep(i)),0.)
    712. 

  712.               qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
    713. 

  713.               qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
    714. 

  714. !-------------------------------------------------------------------
    715. 

  715. ! nsmlt: melting of snow [LH A27]
    716. 

  716. !       (T>T0: ->NR)
    717. 

  717. !-------------------------------------------------------------------
    718. 

  718.               if(qrs(i,k,2).gt.qcrmin) then
    719. 

  719.                 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
    720. 

  720.                 ncr(i,k,3) = ncr(i,k,3) - sfac*psmlt(i,k)
    721. 

  721.               endif
    722. 

  722.               t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
    723. 

  723.             endif
    724. 

  724. !---------------------------------------------------------------
    725. 

  725. ! pgmlt: melting of graupel [HL A23]  [LFO 47]
    726. 

  726. !       (T>T0: QG->QR)
    727. 

  727. !---------------------------------------------------------------
    728. 

  728.             if(qrs(i,k,3).gt.0.) then
    729. 

  729.               coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
    730. 

  730.               pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*(precg1     &
    731. 

  731.                            *rslope2(i,k,3) + precg2*work2(i,k)*coeres)
    732. 

  732.               pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i),                &
    733. 

  733.                           -qrs(i,k,3)/mstep(i)),0.)
    734. 

  734.               qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
    735. 

  735.               qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
    736. 

  736. !-------------------------------------------------------------------
    737. 

  737. ! ngmlt: melting of graupel [LH A28]
    738. 

  738. !       (T>T0: ->NR)
    739. 

  739. !-------------------------------------------------------------------
    740. 

  740.               if(qrs(i,k,3).gt.qcrmin) then
    741. 

  741.                 gfac = rslope(i,k,3)*n0g/qrs(i,k,3)
    742. 

  742.                 ncr(i,k,3) = ncr(i,k,3) - gfac*pgmlt(i,k)
    743. 

  743.               endif
    744. 

  744.               t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
    745. 

  745.             endif
    746. 

  746.           endif
    747. 

  747.         enddo
    748. 

  748.       enddo
    749. 

  749. !---------------------------------------------------------------
    750. 

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

  751. !---------------------------------------------------------------
    752. 

  752.       do k = kte, kts, -1
    753. 

  753.         do i = its, ite
    754. 

  754.           if(qci(i,k,2).le.0.) then
    755. 

  755.             work1c(i,k) = 0.
    756. 

  756.           else
    757. 

  757.             xmi = den(i,k)*qci(i,k,2)/xni(i,k)
    758. 

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

  759.             work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
    760. 

  760.           endif
    761. 

  761.         enddo
    762. 

  762.       enddo
    763. 

  763. !
    764. 

  764. !  forward semi-laglangian scheme (JH), PCM (piecewise constant),  (linear)
    765. 

  765. !
    766. 

  766.       do k = kte, kts, -1
    767. 

  767.         do i = its, ite
    768. 

  768.           denqci(i,k) = den(i,k)*qci(i,k,2)
    769. 

  769.         enddo
    770. 

  770.       enddo
    771. 

  771.       call nislfv_rain_plmr(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci,denqci,  &
    772. 

  772.                            delqi,dtcld,1,0,0)
    773. 

  773.       do k = kts, kte
    774. 

  774.         do i = its, ite
    775. 

  775.           qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
    776. 

  776.         enddo
    777. 

  777.       enddo
    778. 

  778.       do i = its, ite
    779. 

  779.         fallc(i,1) = delqi(i)/delz(i,1)/dtcld
    780. 

  780.       enddo
    781. 

  781. !
    782. 

  782. !----------------------------------------------------------------
    783. 

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

  784. !
    785. 

  785.       do i = its, ite
    786. 

  786.         fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
    787. 

  787.         fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
    788. 

  788.         fallsum_qg = fall(i,kts,3)
    789. 

  789.         if(fallsum.gt.0.) then
    790. 

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

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

  792.         endif
    793. 

  793.         If(fallsum_qsi.gt.0.) then
    794. 

  794.           tstepsnow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + tstepsnow(i)
    795. 

  795.           IF( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
    796. 

  796.             snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snowncv(i)     
    797. 

  797.             snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
    798. 

  798.           ENDIF
    799. 

  799.         ENDIF
    800. 

  800.         IF(fallsum_qg.gt.0.) then
    801. 

  801.           tstepgraup(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000.              &
    802. 

  802.                             + tstepgraup(i)
    803. 

  803.           IF( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN
    804. 

  804.             graupelncv(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000.            &  
    805. 

  805.                             + graupelncv(i)
    806. 

  806.             graupel(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. + graupel(i)
    807. 

  807.           ENDIF
    808. 

  808.         ENDIF
    809. 

  809.         if(fallsum.gt.0.) sr(i) = (tstepsnow(i) + tstepgraup(i))               &
    810. 

  810.                                   /(rainncv(i)+1.e-12)
    811. 

  811.       enddo
    812. 

  812. !
    813. 

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

  814. ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
    815. 

  815. !       (T>T0: QI->QC)
    816. 

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

  817.       do k = kts, kte
    818. 

  818.         do i = its, ite
    819. 

  819.           supcol = t0c-t(i,k)
    820. 

  820.           xlf = xls-xl(i,k)
    821. 

  821.           if(supcol.lt.0.) xlf = xlf0
    822. 

  822.           if(supcol.lt.0 .and. qci(i,k,2).gt.0.) then
    823. 

  823.             qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
    824. 

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

  825. ! nimlt: instantaneous melting of cloud ice  [LH A18]
    826. 

  826. !        (T>T0: ->NC)
    827. 

  827. !--------------------------------------------------------------
    828. 

  828.             ncr(i,k,2) = ncr(i,k,2) + xni(i,k)
    829. 

  829.             t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
    830. 

  830.             qci(i,k,2) = 0.
    831. 

  831.           endif
    832. 

  832. !---------------------------------------------------------------
    833. 

  833. ! pihmf: homogeneous  of cloud water below -40c [HL A45]
    834. 

  834. !        (T<-40C: QC->QI)
    835. 

  835. !---------------------------------------------------------------
    836. 

  836.           if(supcol.gt.40. .and. qci(i,k,1).gt.0.) then
    837. 

  837.             qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
    838. 

  838. !---------------------------------------------------------------
    839. 

  839. ! nihmf: homogeneous  of cloud water below -40c [LH A17]
    840. 

  840. !        (T<-40C: NC->) 
    841. 

  841. !---------------------------------------------------------------
    842. 

  842.             if(ncr(i,k,2).gt.0.) ncr(i,k,2) = 0. 
    843. 

  843.             t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
    844. 

  844.             qci(i,k,1) = 0.
    845. 

  845.           endif
    846. 

  846. !---------------------------------------------------------------
    847. 

  847. ! pihtf: heterogeneous  of cloud water [HL A44]
    848. 

  848. !        (T0>T>-40C: QC->QI)
    849. 

  849. !---------------------------------------------------------------
    850. 

  850.           if(supcol.gt.0. .and. qci(i,k,1).gt.qmin) then
    851. 

  851.             supcolt=min(supcol,70.)
    852. 

  852.             pfrzdtc = min(pi*pi*pfrz1*(exp(pfrz2*supcolt)-1.)*denr/den(i,k)    & 
    853. 

  853.                      *ncr(i,k,2)*rslopec3(i,k)*rslopec3(i,k)/18.*dtcld         &
    854. 

  854.                      ,qci(i,k,1))
    855. 

  855. !---------------------------------------------------------------
    856. 

  856. ! nihtf: heterogeneous  of cloud water [LH A16]
    857. 

  857. !         (T0>T>-40C: NC->) 
    858. 

  858. !---------------------------------------------------------------
    859. 

  859.             if(ncr(i,k,2).gt.ncmin) then
    860. 

  860.               nfrzdtc = min(pi*pfrz1*(exp(pfrz2*supcolt)-1.)*ncr(i,k,2)        &
    861. 

  861.                       *rslopec3(i,k)/6.*dtcld,ncr(i,k,2))
    862. 

  862.               ncr(i,k,2) = ncr(i,k,2) - nfrzdtc
    863. 

  863.             endif
    864. 

  864.             qci(i,k,2) = qci(i,k,2) + pfrzdtc
    865. 

  865.             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
    866. 

  866.             qci(i,k,1) = qci(i,k,1)-pfrzdtc
    867. 

  867.           endif 
    868. 

  868. !---------------------------------------------------------------
    869. 

  869. ! pgfrz:  freezing of rain water [HL A20] [LFO 45]
    870. 

  870. !        (T<T0, QR->QG)
    871. 

  871. !---------------------------------------------------------------
    872. 

  872.           if(supcol.gt.0. .and. qrs(i,k,1).gt.0.) then
    873. 

  873.             supcolt=min(supcol,70.)
    874. 

  874.             pfrzdtr = min(140.*(pi*pi)*pfrz1*ncr(i,k,3)*denr/den(i,k)          &
    875. 

  875.                   *(exp(pfrz2*supcolt)-1.)*rslope3(i,k,1)*rslope3(i,k,1)       & 
    876. 

  876.                   *dtcld,qrs(i,k,1))        
    877. 

  877. !---------------------------------------------------------------
    878. 

  878. ! ngfrz: freezing of rain water [LH A26]
    879. 

  879. !        (T<T0, NR-> ) 
    880. 

  880. !---------------------------------------------------------------
    881. 

  881.             if(ncr(i,k,3).gt.nrmin) then
    882. 

  882.               nfrzdtr = min(4.*pi*pfrz1*ncr(i,k,3)*(exp(pfrz2*supcolt)-1.)     &
    883. 

  883.                        *rslope3(i,k,1)*dtcld, ncr(i,k,3)) 
    884. 

  884.               ncr(i,k,3) = ncr(i,k,3) - nfrzdtr
    885. 

  885.             endif
    886. 

  886.             qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
    887. 

  887.             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
    888. 

  888.             qrs(i,k,1) = qrs(i,k,1) - pfrzdtr
    889. 

  889.           endif
    890. 

  890.         enddo
    891. 

  891.       enddo
    892. 

  892. !
    893. 

  893.       do k = kts, kte
    894. 

  894.         do i = its, ite
    895. 

  895.           ncr(i,k,2) = max(ncr(i,k,2),0.0)
    896. 

  896.           ncr(i,k,3) = max(ncr(i,k,3),0.0)
    897. 

  897.         enddo
    898. 

  898.       enddo
    899. 

  899. !
    900. 

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

  901. !     update the slope parameters for microphysics computation
    902. 

  902. !
    903. 

  903.       do k = kts, kte
    904. 

  904.         do i = its, ite
    905. 

  905.           qrs_tmp(i,k,1) = qrs(i,k,1)
    906. 

  906.           qrs_tmp(i,k,2) = qrs(i,k,2)
    907. 

  907.           qrs_tmp(i,k,3) = qrs(i,k,3)
    908. 

  908.           ncr_tmp(i,k) = ncr(i,k,3)
    909. 

  909.         enddo
    910. 

  910.       enddo
    911. 

  911.       call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
    912. 

  912.                      rslope3,work1,workn,its,ite,kts,kte)
    913. 

  913.       do k = kts, kte
    914. 

  914.         do i = its, ite
    915. 

  915. !-----------------------------------------------------------------
    916. 

  916. ! compute the mean-volume drop diameter                  [LH A10] 
    917. 

  917. ! for raindrop distribution                          
    918. 

  918. !-----------------------------------------------------------------
    919. 

  919.           avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
    920. 

  920. !
    921. 

  921.           if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin) then
    922. 

  922.             rslopec(i,k) = rslopecmax
    923. 

  923.             rslopec2(i,k) = rslopec2max
    924. 

  924.             rslopec3(i,k) = rslopec3max
    925. 

  925.           else
    926. 

  926.             rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
    927. 

  927.             rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
    928. 

  928.             rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
    929. 

  929.           endif
    930. 

  930. !--------------------------------------------------------------------
    931. 

  931. ! compute the mean-volume drop diameter                   [LH A7]
    932. 

  932. ! for cloud-droplet distribution
    933. 

  933. !--------------------------------------------------------------------
    934. 

  934.           avedia(i,k,1) = rslopec(i,k)
    935. 

  935.         enddo
    936. 

  936.       enddo
    937. 

  937. !
    938. 

  938.       do k = kts, kte
    939. 

  939.         do i = its, ite
    940. 

  940.           work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
    941. 

  941.           work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
    942. 

  942.           work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
    943. 

  943.         enddo
    944. 

  944.       enddo
    945. 

  945. !
    946. 

  946. !===============================================================
    947. 

  947. !
    948. 

  948. ! warm rain processes
    949. 

  949. !
    950. 

  950. ! - follows the double-moment processes in Lim and Hong
    951. 

  951. !
    952. 

  952. !===============================================================
    953. 

  953. !
    954. 

  954.       do k = kts, kte
    955. 

  955.         do i = its, ite
    956. 

  956.           supsat = max(q(i,k),qmin)-qs(i,k,1)
    957. 

  957.           satdt = supsat/dtcld
    958. 

  958. !---------------------------------------------------------------
    959. 

  959. ! praut: auto conversion rate from cloud to rain  [LH 9] [CP 17]
    960. 

  960. !        (QC->QR)
    961. 

  961. !--------------------------------------------------------------
    962. 

  962.           lencon  = 2.7e-2*den(i,k)*qci(i,k,1)*(1.e20/16.*rslopec2(i,k)        &
    963. 

  963.                    *rslopec2(i,k)-0.4)
    964. 

  964.           lenconcr = max(1.2*lencon, qcrmin)
    965. 

  965.           if(avedia(i,k,1).gt.di15) then
    966. 

  966.             taucon = 3.7/den(i,k)/qci(i,k,1)/(0.5e6*rslopec(i,k)-7.5)
    967. 

  967.             taucon = max(taucon, 1.e-12)
    968. 

  968.             praut(i,k) = lencon/(taucon*den(i,k))
    969. 

  969.             praut(i,k) = min(max(praut(i,k),0.),qci(i,k,1)/dtcld)
    970. 

  970. !---------------------------------------------------------------
    971. 

  971. ! nraut: auto conversion rate from cloud to rain [LH A6] [CP 18 & 19]
    972. 

  972. !        (NC->NR)
    973. 

  973. !---------------------------------------------------------------
    974. 

  974.             nraut(i,k) = 3.5e9*den(i,k)*praut(i,k)
    975. 

  975.             if(qrs(i,k,1).gt.lenconcr)                                         &
    976. 

  976.             nraut(i,k) = ncr(i,k,3)/qrs(i,k,1)*praut(i,k)
    977. 

  977.             nraut(i,k) = min(nraut(i,k),ncr(i,k,2)/dtcld)
    978. 

  978.           endif
    979. 

  979. !---------------------------------------------------------------
    980. 

  980. ! pracw: accretion of cloud water by rain     [LH 10] [CP 22 & 23]
    981. 

  981. !        (QC->QR)
    982. 

  982. ! nracw: accretion of cloud water by rain     [LH A9]
    983. 

  983. !        (NC->)
    984. 

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

  985.           if(qrs(i,k,1).ge.lenconcr) then
    986. 

  986.             if(avedia(i,k,2).ge.di100) then
    987. 

  987.               nracw(i,k) = min(ncrk1*ncr(i,k,2)*ncr(i,k,3)*(rslopec3(i,k)      &
    988. 

  988.                          + 24.*rslope3(i,k,1)),ncr(i,k,2)/dtcld)
    989. 

  989.               pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk1*ncr(i,k,2)          &
    990. 

  990.                          *ncr(i,k,3)*rslopec3(i,k)*(2.*rslopec3(i,k)           &
    991. 

  991.                          + 24.*rslope3(i,k,1)),qci(i,k,1)/dtcld)   
    992. 

  992.             else
    993. 

  993.               nracw(i,k) = min(ncrk2*ncr(i,k,2)*ncr(i,k,3)*(2.*rslopec3(i,k)   &
    994. 

  994.                          *rslopec3(i,k)+5040.*rslope3(i,k,1)                   &
    995. 

  995.                          *rslope3(i,k,1)),ncr(i,k,2)/dtcld)
    996. 

  996.               pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk2*ncr(i,k,2)          &
    997. 

  997.                          *ncr(i,k,3)*rslopec3(i,k)*(6.*rslopec3(i,k)           &     
    998. 

  998.                          *rslopec3(i,k)+5040.*rslope3(i,k,1)*rslope3(i,k,1))   & 
    999. 

  999.                          ,qci(i,k,1)/dtcld)
    1000. 

  1000.             endif
    1001. 

  1001.           endif 
    1002. 

  1002. !----------------------------------------------------------------
    1003. 

  1003. ! nccol: self collection of cloud water             [LH A8] [CP 24 & 25]     
    1004. 

  1004. !        (NC->)
    1005. 

  1005. !----------------------------------------------------------------
    1006. 

  1006.           if(avedia(i,k,1).ge.di100) then
    1007. 

  1007.             nccol(i,k) = ncrk1*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)
    1008. 

  1008.           else
    1009. 

  1009.             nccol(i,k) = 2.*ncrk2*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)        &     
    1010. 

  1010.                          *rslopec3(i,k)
    1011. 

  1011.           endif
    1012. 

  1012. !----------------------------------------------------------------
    1013. 

  1013. ! nrcol: self collection of rain-drops and break-up [LH A21] [CP 24 & 25]
    1014. 

  1014. !        (NR->)
    1015. 

  1015. !----------------------------------------------------------------
    1016. 

  1016.           if(qrs(i,k,1).ge.lenconcr) then
    1017. 

  1017.             if(avedia(i,k,2).lt.di100) then 
    1018. 

  1018.               nrcol(i,k) = 5040.*ncrk2*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1)    &
    1019. 

  1019.                           *rslope3(i,k,1)
    1020. 

  1020.             elseif(avedia(i,k,2).ge.di100 .and. avedia(i,k,2).lt.di600) then
    1021. 

  1021.               nrcol(i,k) = 24.*ncrk1*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1)
    1022. 

  1022.             elseif(avedia(i,k,2).ge.di600 .and. avedia(i,k,2).lt.di2000) then
    1023. 

  1023.               coecol = -2.5e3*(avedia(i,k,2)-di600) 
    1024. 

  1024.               nrcol(i,k) = 24.*exp(coecol)*ncrk1*ncr(i,k,3)*ncr(i,k,3)         &
    1025. 

  1025.                          *rslope3(i,k,1)
    1026. 

  1026.             else
    1027. 

  1027.               nrcol(i,k) = 0.
    1028. 

  1028.             endif
    1029. 

  1029.           endif
    1030. 

  1030. !---------------------------------------------------------------
    1031. 

  1031. ! prevp: evaporation/condensation rate of rain   [HL A41]  
    1032. 

  1032. !        (QV->QR or QR->QV)
    1033. 

  1033. !---------------------------------------------------------------
    1034. 

  1034.           if(qrs(i,k,1).gt.0.) then
    1035. 

  1035.             coeres = rslope(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
    1036. 

  1036.             prevp(i,k) = (rh(i,k,1)-1.)*ncr(i,k,3)*(precr1*rslope(i,k,1)       &
    1037. 

  1037.                          + precr2*work2(i,k)*coeres)/work1(i,k,1)
    1038. 

  1038.             if(prevp(i,k).lt.0.) then
    1039. 

  1039.               prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
    1040. 

  1040.               prevp(i,k) = max(prevp(i,k),satdt/2)
    1041. 

  1041. !----------------------------------------------------------------
    1042. 

  1042. ! Nrevp: evaporation/condensation rate of rain   [LH A14] 
    1043. 

  1043. !        (NR->NCCN) 
    1044. 

  1044. !----------------------------------------------------------------
    1045. 

  1045.               if(prevp(i,k).eq.-qrs(i,k,1)/dtcld) then
    1046. 

  1046.                 ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,3)
    1047. 

  1047.                 ncr(i,k,3) = 0.
    1048. 

  1048.               endif
    1049. 

  1049.             else
    1050. 

  1050. !
    1051. 

  1051.               prevp(i,k) = min(prevp(i,k),satdt/2)
    1052. 

  1052.             endif
    1053. 

  1053.           endif
    1054. 

  1054.         enddo
    1055. 

  1055.       enddo
    1056. 

  1056. !
    1057. 

  1057. !===============================================================
    1058. 

  1058. !
    1059. 

  1059. ! cold rain processes
    1060. 

  1060. !
    1061. 

  1061. ! - follows the revised ice microphysics processes in HDC
    1062. 

  1062. ! - the processes same as in RH83 and RH84  and LFO behave
    1063. 

  1063. !   following ice crystal hapits defined in HDC, inclduing
    1064. 

  1064. !   intercept parameter for snow (n0s), ice crystal number
    1065. 

  1065. !   concentration (ni), ice nuclei number concentration
    1066. 

  1066. !   (n0i), ice diameter (d)
    1067. 

  1067. !
    1068. 

  1068. !===============================================================
    1069. 

  1069. !
    1070. 

  1070.       do k = kts, kte
    1071. 

  1071.         do i = its, ite
    1072. 

  1072.           supcol = t0c-t(i,k)
    1073. 

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

  1074.           supsat = max(q(i,k),qmin)-qs(i,k,2)
    1075. 

  1075.           satdt = supsat/dtcld
    1076. 

  1076.           ifsat = 0
    1077. 

  1077. !-------------------------------------------------------------
    1078. 

  1078. ! Ni: ice crystal number concentraiton   [HDC 5c]
    1079. 

  1079. !-------------------------------------------------------------
    1080. 

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

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

  1082.           temp = (den(i,k)*max(qci(i,k,2),qmin))
    1083. 

  1083.           temp = sqrt(sqrt(temp*temp*temp))
    1084. 

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

  1085.           eacrs = exp(0.07*(-supcol))
    1086. 

  1086. !
    1087. 

  1087.           xmi = den(i,k)*qci(i,k,2)/xni(i,k)
    1088. 

  1088.           diameter  = min(dicon * sqrt(xmi),dimax)
    1089. 

  1089.           vt2i = 1.49e4*diameter**1.31
    1090. 

  1090.           vt2r=pvtr*rslopeb(i,k,1)*denfac(i,k)
    1091. 

  1091.           vt2s=pvts*rslopeb(i,k,2)*denfac(i,k)
    1092. 

  1092.           vt2g=pvtg*rslopeb(i,k,3)*denfac(i,k)
    1093. 

  1093.           qsum(i,k) = max((qrs(i,k,2)+qrs(i,k,3)),1.e-15)
    1094. 

  1094.           if(qsum(i,k) .gt. 1.e-15) then
    1095. 

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

  1096.           else    
    1097. 

  1097.             vt2ave=0.
    1098. 

  1098.           endif
    1099. 

  1099.           if(supcol.gt.0. .and. qci(i,k,2).gt.qmin) then
    1100. 

  1100.             if(qrs(i,k,1).gt.qcrmin) then
    1101. 

  1101. !-------------------------------------------------------------
    1102. 

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

  1103. !        (T<T0: QI->QR)
    1104. 

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

  1105.               acrfac = 6.*rslope2(i,k,1)+4.*diameter*rslope(i,k,1) + diameter**2
    1106. 

  1106.               praci(i,k) = pi*qci(i,k,2)*ncr(i,k,3)*abs(vt2r-vt2i)*acrfac/4.
    1107. 

  1107.               praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
    1108. 

  1108. !-------------------------------------------------------------
    1109. 

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

  1110. !        (T<T0: QR->QS or QR->QG)
    1111. 

  1111. !-------------------------------------------------------------
    1112. 

  1112.               piacr(i,k) = pi*pi*avtr*ncr(i,k,3)*denr*xni(i,k)*denfac(i,k)     &
    1113. 

  1113.                           *g7pbr*rslope3(i,k,1)*rslope2(i,k,1)*rslopeb(i,k,1)  &
    1114. 

  1114.                           /24./den(i,k)
    1115. 

  1115.               piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
    1116. 

  1116.             endif
    1117. 

  1117. !-------------------------------------------------------------
    1118. 

  1118. ! niacr: Accretion of rain by cloud ice  [LH A25]
    1119. 

  1119. !        (T<T0: NR->) 
    1120. 

  1120. !-------------------------------------------------------------
    1121. 

  1121.             if(ncr(i,k,3).gt.nrmin) then
    1122. 

  1122.               niacr(i,k) = pi*avtr*ncr(i,k,3)*xni(i,k)*denfac(i,k)*g4pbr       &
    1123. 

  1123.                           *rslope2(i,k,1)*rslopeb(i,k,1)/4.
    1124. 

  1124.               niacr(i,k) = min(niacr(i,k),ncr(i,k,3)/dtcld)
    1125. 

  1125.             endif
    1126. 

  1126. !-------------------------------------------------------------
    1127. 

  1127. ! psaci: Accretion of cloud ice by snow [HDC 10]
    1128. 

  1128. !        (T<T0: QI->QS)
    1129. 

  1129. !-------------------------------------------------------------
    1130. 

  1130.             if(qrs(i,k,2).gt.qcrmin) then
    1131. 

  1131.               acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2)            &
    1132. 

  1132.                       + diameter**2*rslope(i,k,2)
    1133. 

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

  1134.                           *abs(vt2ave-vt2i)*acrfac/4.
    1135. 

  1135.               psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
    1136. 

  1136.             endif
    1137. 

  1137. !-------------------------------------------------------------
    1138. 

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

  1139. !        (T<T0: QI->QG)
    1140. 

  1140. !-------------------------------------------------------------
    1141. 

  1141.             if(qrs(i,k,3).gt.qcrmin) then
    1142. 

  1142.               egi = exp(0.07*(-supcol))
    1143. 

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

  1144.                       + diameter**2*rslope(i,k,3)
    1145. 

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

  1146.               pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
    1147. 

  1147.             endif
    1148. 

  1148.           endif
    1149. 

  1149. !-------------------------------------------------------------
    1150. 

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

  1151. !        (T<T0: QC->QS, and T>=T0: QC->QR)
    1152. 

  1152. !-------------------------------------------------------------
    1153. 

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

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

  1155.                         *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
    1156. 

  1156.           endif
    1157. 

  1157. !-------------------------------------------------------------
    1158. 

  1158. ! nsacw: Accretion of cloud water by snow  [LH A12]
    1159. 

  1159. !        (NC ->) 
    1160. 

  1160. !-------------------------------------------------------------
    1161. 

  1161.          if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
    1162. 

  1162.            nsacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2)    &
    1163. 

  1163.                        *ncr(i,k,2)*denfac(i,k),ncr(i,k,2)/dtcld)
    1164. 

  1164.          endif
    1165. 

  1165. !-------------------------------------------------------------
    1166. 

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

  1167. !        (T<T0: QC->QG, and T>=T0: QC->QR)
    1168. 

  1168. !-------------------------------------------------------------
    1169. 

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

  1170.             pgacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3)*qci(i,k,1)    &
    1171. 

  1171.                         *denfac(i,k),qci(i,k,1)/dtcld)
    1172. 

  1172.           endif
    1173. 

  1173. !-------------------------------------------------------------
    1174. 

  1174. ! ngacw: Accretion of cloud water by graupel [LH A13]
    1175. 

  1175. !        (NC-> 
    1176. 

  1176. !-------------------------------------------------------------
    1177. 

  1177.           if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
    1178. 

  1178.             ngacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3)*ncr(i,k,2)    &
    1179. 

  1179.                         *denfac(i,k),ncr(i,k,2)/dtcld)
    1180. 

  1180.           endif
    1181. 

  1181. !-------------------------------------------------------------
    1182. 

  1182. ! paacw: Accretion of cloud water by averaged snow/graupel 
    1183. 

  1183. !        (T<T0: QC->QG or QS, and T>=T0: QC->QR) 
    1184. 

  1184. !-------------------------------------------------------------
    1185. 

  1185.           if(qrs(i,k,2).gt.qcrmin .and. qrs(i,k,3).gt.qcrmin) then
    1186. 

  1186.             paacw(i,k) = (qrs(i,k,2)*psacw(i,k)+qrs(i,k,3)*pgacw(i,k))/(qsum(i,k))
    1187. 

  1187. !-------------------------------------------------------------
    1188. 

  1188. ! naacw: Accretion of cloud water by averaged snow/graupel
    1189. 

  1189. !        (Nc->)
    1190. 

  1190. !-------------------------------------------------------------
    1191. 

  1191.             naacw(i,k) = (qrs(i,k,2)*nsacw(i,k)+qrs(i,k,3)*ngacw(i,k))/(qsum(i,k))
    1192. 

  1192.           endif      
    1193. 

  1193. !-------------------------------------------------------------
    1194. 

  1194. ! pracs: Accretion of snow by rain [HL A11] [LFO 27]
    1195. 

  1195. !         (T<T0: QS->QG)
    1196. 

  1196. !-------------------------------------------------------------
    1197. 

  1197.           if(qrs(i,k,2).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin) then
    1198. 

  1198.             if(supcol.gt.0) then
    1199. 

  1199.               acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2)                        &
    1200. 

  1200.                       + 4.*rslope3(i,k,2)*rslope2(i,k,2)*rslope(i,k,1)         &
    1201. 

  1201.                       + 1.5*rslope2(i,k,2)*rslope2(i,k,2)*rslope2(i,k,1)
    1202. 

  1202.               pracs(i,k) = pi*pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2r-vt2ave)   &
    1203. 

  1203.                           *(dens/den(i,k))*acrfac
    1204. 

  1204.               pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
    1205. 

  1205.             endif
    1206. 

  1206. !-------------------------------------------------------------
    1207. 

  1207. ! psacr: Accretion of rain by snow [HL A10] [LFO 28]
    1208. 

  1208. !         (T<T0:QR->QS or QR->QG) (T>=T0: enhance melting of snow)
    1209. 

  1209. !-------------------------------------------------------------
    1210. 

  1210.             acrfac = 30.*rslope3(i,k,1)*rslope2(i,k,1)*rslope(i,k,2)           &
    1211. 

  1211.                      + 5.*rslope2(i,k,1)*rslope2(i,k,1)*rslope2(i,k,2)         &
    1212. 

  1212.                      + 2.*rslope3(i,k,1)*rslope3(i,k,2)
    1213. 

  1213.             psacr(i,k) = pi*pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2ave-vt2r)     &
    1214. 

  1214.                         *(denr/den(i,k))*acrfac
    1215. 

  1215.             psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
    1216. 

  1216.           endif
    1217. 

  1217.           if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,3).gt.nrmin) then
    1218. 

  1218. !-------------------------------------------------------------
    1219. 

  1219. ! nsacr: Accretion of rain by snow  [LH A23] 
    1220. 

  1220. !       (T<T0:NR->)
    1221. 

  1221. !-------------------------------------------------------------
    1222. 

  1222.             acrfac = 1.5*rslope2(i,k,1)*rslope(i,k,2)                          &
    1223. 

  1223.                     + 1.0*rslope(i,k,1)*rslope2(i,k,2)+.5*rslope3(i,k,2)        
    1224. 

  1224.             nsacr(i,k) = pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2ave-vt2r)        &
    1225. 

  1225.                         *acrfac
    1226. 

  1226.             nsacr(i,k) = min(nsacr(i,k),ncr(i,k,3)/dtcld)
    1227. 

  1227.           endif
    1228. 

  1228. !-------------------------------------------------------------
    1229. 

  1229. ! pgacr: Accretion of rain by graupel [HL A12] [LFO 42]
    1230. 

  1230. !         (T<T0: QR->QG) (T>=T0: enhance melting of graupel)
    1231. 

  1231. !-------------------------------------------------------------
    1232. 

  1232.           if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin) then
    1233. 

  1233.             acrfac = 30.*rslope3(i,k,1)*rslope2(i,k,1)*rslope(i,k,3)           &
    1234. 

  1234.                     + 5.*rslope2(i,k,1)*rslope2(i,k,1)*rslope2(i,k,3)          &
    1235. 

  1235.                     + 2.*rslope3(i,k,1)*rslope3(i,k,3)
    1236. 

  1236.             pgacr(i,k) = pi*pi*ncr(i,k,3)*n0g*abs(vt2ave-vt2r)*(denr/den(i,k)) &
    1237. 

  1237.                         *acrfac
    1238. 

  1238.             pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
    1239. 

  1239.           endif
    1240. 

  1240. !-------------------------------------------------------------
    1241. 

  1241. ! ngacr: Accretion of rain by graupel  [LH A24] 
    1242. 

  1242. !        (T<T0: NR->) 
    1243. 

  1243. !-------------------------------------------------------------
    1244. 

  1244.           if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,3).gt.nrmin) then
    1245. 

  1245.             acrfac = 1.5*rslope2(i,k,1)*rslope(i,k,3)                          &
    1246. 

  1246.                     + 1.0*rslope(i,k,1)*rslope2(i,k,3) + .5*rslope3(i,k,3)   
    1247. 

  1247.             ngacr(i,k) = pi*ncr(i,k,3)*n0g*abs(vt2ave-vt2r)*acrfac
    1248. 

  1248.             ngacr(i,k) = min(ngacr(i,k),ncr(i,k,3)/dtcld)
    1249. 

  1249.           endif
    1250. 

  1250. !
    1251. 

  1251. !-------------------------------------------------------------
    1252. 

  1252. ! pgacs: Accretion of snow by graupel [HL A13] [LFO 29]
    1253. 

  1253. !        (QS->QG) : This process is eliminated in V3.0 with the
    1254. 

  1254. !        new combined snow/graupel fall speeds
    1255. 

  1255. !-------------------------------------------------------------
    1256. 

  1256.           if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,2).gt.qcrmin) then
    1257. 

  1257.             pgacs(i,k) = 0. 
    1258. 

  1258.           endif
    1259. 

  1259.           if(supcol.le.0) then
    1260. 

  1260.             xlf = xlf0
    1261. 

  1261. !-------------------------------------------------------------
    1262. 

  1262. ! pseml: Enhanced melting of snow by accretion of water [HL A34]
    1263. 

  1263. !        (T>=T0: QS->QR)
    1264. 

  1264. !-------------------------------------------------------------
    1265. 

  1265.             if(qrs(i,k,2).gt.0.)                                               & 
    1266. 

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

  1267.                           /xlf,-qrs(i,k,2)/dtcld),0.)
    1268. 

  1268. !--------------------------------------------------------------
    1269. 

  1269. ! nseml: Enhanced melting of snow by accretion of water    [LH A29]
    1270. 

  1270. !        (T>=T0: ->NR)
    1271. 

  1271. !--------------------------------------------------------------
    1272. 

  1272.               if  (qrs(i,k,2).gt.qcrmin) then
    1273. 

  1273.                 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
    1274. 

  1274.                 nseml(i,k) = -sfac*pseml(i,k)
    1275. 

  1275.               endif
    1276. 

  1276. !-------------------------------------------------------------
    1277. 

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

  1278. !        (T>=T0: QG->QR)
    1279. 

  1279. !-------------------------------------------------------------
    1280. 

  1280.             if(qrs(i,k,3).gt.0.)                                               &
    1281. 

  1281.               pgeml(i,k) = min(max(cliq*supcol*(paacw(i,k)+pgacr(i,k))/xlf     &
    1282. 

  1282.                           ,-qrs(i,k,3)/dtcld),0.)
    1283. 

  1283. !--------------------------------------------------------------
    1284. 

  1284. ! ngeml: Enhanced melting of graupel by accretion of water [LH A30] 
    1285. 

  1285. !         (T>=T0: -> NR)
    1286. 

  1286. !--------------------------------------------------------------
    1287. 

  1287.               if (qrs(i,k,3).gt.qcrmin) then
    1288. 

  1288.                 gfac = rslope(i,k,3)*n0g/qrs(i,k,3)
    1289. 

  1289.                 ngeml(i,k) = -gfac*pgeml(i,k)
    1290. 

  1290.               endif
    1291. 

  1291.           endif
    1292. 

  1292.           if(supcol.gt.0) then
    1293. 

  1293. !-------------------------------------------------------------
    1294. 

  1294. ! pidep: Deposition/Sublimation rate of ice [HDC 9]
    1295. 

  1295. !       (T<T0: QV->QI or QI->QV)
    1296. 

  1296. !-------------------------------------------------------------
    1297. 

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

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

  1299.               supice = satdt-prevp(i,k)
    1300. 

  1300.               if(pidep(i,k).lt.0.) then
    1301. 

  1301.                 pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
    1302. 

  1302.                 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
    1303. 

  1303.               else
    1304. 

  1304.                 pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
    1305. 

  1305.               endif
    1306. 

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

  1307.             endif
    1308. 

  1308. !-------------------------------------------------------------
    1309. 

  1309. ! psdep: deposition/sublimation rate of snow [HDC 14]
    1310. 

  1310. !        (T<T0: QV->QS or QS->QV)
    1311. 

  1311. !-------------------------------------------------------------
    1312. 

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

  1313.               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    1314. 

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

  1315.                            + precs2*work2(i,k)*coeres)/work1(i,k,2)
    1316. 

  1316.               supice = satdt-prevp(i,k)-pidep(i,k)
    1317. 

  1317.               if(psdep(i,k).lt.0.) then
    1318. 

  1318.                 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
    1319. 

  1319.                 psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
    1320. 

  1320.               else
    1321. 

  1321.                 psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
    1322. 

  1322.               endif
    1323. 

  1323.               if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) ifsat = 1
    1324. 

  1324.             endif
    1325. 

  1325. !-------------------------------------------------------------
    1326. 

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

  1327. !        (T<T0: QV->QG or QG->QV)
    1328. 

  1328. !-------------------------------------------------------------
    1329. 

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

  1330.               coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
    1331. 

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

  1332.                           + precg2*work2(i,k)*coeres)/work1(i,k,2)
    1333. 

  1333.               supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
    1334. 

  1334.               if(pgdep(i,k).lt.0.) then
    1335. 

  1335.                 pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
    1336. 

  1336.                 pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
    1337. 

  1337.               else
    1338. 

  1338.                 pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
    1339. 

  1339.               endif
    1340. 

  1340.               if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge.          &
    1341. 

  1341.                 abs(satdt)) ifsat = 1
    1342. 

  1342.             endif
    1343. 

  1343. !-------------------------------------------------------------
    1344. 

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

  1345. !       (T<T0: QV->QI)
    1346. 

  1346. !-------------------------------------------------------------
    1347. 

  1347.             if(supsat.gt.0. .and. ifsat.ne.1) then
    1348. 

  1348.               supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
    1349. 

  1349.               xni0 = 1.e3*exp(0.1*supcol)
    1350. 

  1350.               roqi0 = 4.92e-11*xni0**1.33
    1351. 

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

  1352.               pigen(i,k) = min(min(pigen(i,k),satdt),supice)
    1353. 

  1353.             endif
    1354. 

  1354. !
    1355. 

  1355. !-------------------------------------------------------------
    1356. 

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

  1357. !        (T<T0: QI->QS)
    1358. 

  1358. !-------------------------------------------------------------
    1359. 

  1359.             if(qci(i,k,2).gt.0.) then
    1360. 

  1360.               qimax = roqimax/den(i,k)
    1361. 

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

  1362.             endif
    1363. 

  1363. !
    1364. 

  1364. !-------------------------------------------------------------
    1365. 

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

  1366. !        (T<T0: QS->QG)
    1367. 

  1367. !-------------------------------------------------------------
    1368. 

  1368.             if(qrs(i,k,2).gt.0.) then
    1369. 

  1369.               alpha2 = 1.e-3*exp(0.09*(-supcol))
    1370. 

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

  1371.             endif
    1372. 

  1372.           endif
    1373. 

  1373. !
    1374. 

  1374. !-------------------------------------------------------------
    1375. 

  1375. ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
    1376. 

  1376. !       (T>=T0: QS->QV)
    1377. 

  1377. !-------------------------------------------------------------
    1378. 

  1378.           if(supcol.lt.0.) then
    1379. 

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

  1380.               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    1381. 

  1381.               psevp(i,k) = (rh(i,k,1)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2)   &
    1382. 

  1382.                            +precs2*work2(i,k)*coeres)/work1(i,k,1)
    1383. 

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

  1384.             endif
    1385. 

  1385. !-------------------------------------------------------------
    1386. 

  1386. ! pgevp: Evaporation of melting graupel [HL A25] [RH84 A19]
    1387. 

  1387. !       (T>=T0: QG->QV)
    1388. 

  1388. !-------------------------------------------------------------
    1389. 

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

  1390.               coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
    1391. 

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

  1392.                          + precg2*work2(i,k)*coeres)/work1(i,k,1)
    1393. 

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

  1394.             endif
    1395. 

  1395.           endif
    1396. 

  1396.         enddo
    1397. 

  1397.       enddo
    1398. 

  1398. !
    1399. 

  1399. !
    1400. 

  1400. !----------------------------------------------------------------
    1401. 

  1401. !     check mass conservation of generation terms and feedback to the
    1402. 

  1402. !     large scale
    1403. 

  1403. !
    1404. 

  1404.       do k = kts, kte
    1405. 

  1405.         do i = its, ite
    1406. 

  1406. !
    1407. 

  1407.           delta2=0.
    1408. 

  1408.           delta3=0.
    1409. 

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

  1410.           if(qrs(i,k,1).lt.1.e-4) delta3=1.
    1411. 

  1411.           if(t(i,k).le.t0c) then
    1412. 

  1412. !
    1413. 

  1413. !     cloud water
    1414. 

  1414. !
    1415. 

  1415.             value = max(qmin,qci(i,k,1))
    1416. 

  1416.             source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))             &
    1417. 

  1417.                     *dtcld
    1418. 

  1418.             if (source.gt.value) then
    1419. 

  1419.               factor = value/source
    1420. 

  1420.               praut(i,k) = praut(i,k)*factor
    1421. 

  1421.               pracw(i,k) = pracw(i,k)*factor
    1422. 

  1422.               paacw(i,k) = paacw(i,k)*factor
    1423. 

  1423.             endif
    1424. 

  1424. !
    1425. 

  1425. !     cloud ice
    1426. 

  1426. !
    1427. 

  1427.             value = max(qmin,qci(i,k,2))
    1428. 

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

  1429.                     +pgaci(i,k))*dtcld
    1430. 

  1430.             if (source.gt.value) then
    1431. 

  1431.               factor = value/source
    1432. 

  1432.               psaut(i,k) = psaut(i,k)*factor
    1433. 

  1433.               pigen(i,k) = pigen(i,k)*factor
    1434. 

  1434.               pidep(i,k) = pidep(i,k)*factor
    1435. 

  1435.               praci(i,k) = praci(i,k)*factor
    1436. 

  1436.               psaci(i,k) = psaci(i,k)*factor
    1437. 

  1437.               pgaci(i,k) = pgaci(i,k)*factor
    1438. 

  1438.             endif
    1439. 

  1439. !
    1440. 

  1440. !     rain
    1441. 

  1441. !
    1442. 

  1442.             value = max(qmin,qrs(i,k,1))
    1443. 

  1443.             source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k)             &
    1444. 

  1444.                     +psacr(i,k)+pgacr(i,k))*dtcld
    1445. 

  1445.             if (source.gt.value) then
    1446. 

  1446.               factor = value/source
    1447. 

  1447.               praut(i,k) = praut(i,k)*factor
    1448. 

  1448.               prevp(i,k) = prevp(i,k)*factor
    1449. 

  1449.               pracw(i,k) = pracw(i,k)*factor
    1450. 

  1450.               piacr(i,k) = piacr(i,k)*factor
    1451. 

  1451.               psacr(i,k) = psacr(i,k)*factor
    1452. 

  1452.               pgacr(i,k) = pgacr(i,k)*factor
    1453. 

  1453.             endif
    1454. 

  1454. !
    1455. 

  1455. !     snow
    1456. 

  1456. !
    1457. 

  1457.             value = max(qmin,qrs(i,k,2))
    1458. 

  1458.             source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k)             &
    1459. 

  1459.                      +piacr(i,k)*delta3+praci(i,k)*delta3                      &
    1460. 

  1460.                      -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2                 &
    1461. 

  1461.                      +psaci(i,k)-pgacs(i,k) )*dtcld
    1462. 

  1462.             if (source.gt.value) then
    1463. 

  1463.               factor = value/source
    1464. 

  1464.               psdep(i,k) = psdep(i,k)*factor
    1465. 

  1465.               psaut(i,k) = psaut(i,k)*factor
    1466. 

  1466.               pgaut(i,k) = pgaut(i,k)*factor
    1467. 

  1467.               paacw(i,k) = paacw(i,k)*factor
    1468. 

  1468.               piacr(i,k) = piacr(i,k)*factor
    1469. 

  1469.               praci(i,k) = praci(i,k)*factor
    1470. 

  1470.               psaci(i,k) = psaci(i,k)*factor
    1471. 

  1471.               pracs(i,k) = pracs(i,k)*factor
    1472. 

  1472.               psacr(i,k) = psacr(i,k)*factor
    1473. 

  1473.               pgacs(i,k) = pgacs(i,k)*factor
    1474. 

  1474.             endif
    1475. 

  1475. !
    1476. 

  1476. !     graupel
    1477. 

  1477. !
    1478. 

  1478.             value = max(qmin,qrs(i,k,3))
    1479. 

  1479.             source = -(pgdep(i,k)+pgaut(i,k)                                   &
    1480. 

  1480.                      +piacr(i,k)*(1.-delta3)+praci(i,k)*(1.-delta3)            &
    1481. 

  1481.                      +psacr(i,k)*(1.-delta2)+pracs(i,k)*(1.-delta2)            &
    1482. 

  1482.                      +pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
    1483. 

  1483.             if (source.gt.value) then
    1484. 

  1484.               factor = value/source
    1485. 

  1485.               pgdep(i,k) = pgdep(i,k)*factor
    1486. 

  1486.               pgaut(i,k) = pgaut(i,k)*factor
    1487. 

  1487.               piacr(i,k) = piacr(i,k)*factor
    1488. 

  1488.               praci(i,k) = praci(i,k)*factor
    1489. 

  1489.               psacr(i,k) = psacr(i,k)*factor
    1490. 

  1490.               pracs(i,k) = pracs(i,k)*factor
    1491. 

  1491.               paacw(i,k) = paacw(i,k)*factor
    1492. 

  1492.               pgaci(i,k) = pgaci(i,k)*factor
    1493. 

  1493.               pgacr(i,k) = pgacr(i,k)*factor
    1494. 

  1494.               pgacs(i,k) = pgacs(i,k)*factor
    1495. 

  1495.             endif
    1496. 

  1496. !
    1497. 

  1497. !     cloud
    1498. 

  1498. !
    1499. 

  1499.             value = max(ncmin,ncr(i,k,2))
    1500. 

  1500.             source = (nraut(i,k)+nccol(i,k)+nracw(i,k)                         &
    1501. 

  1501.                     +naacw(i,k)+naacw(i,k))*dtcld
    1502. 

  1502.             if (source.gt.value) then
    1503. 

  1503.               factor = value/source
    1504. 

  1504.               nraut(i,k) = nraut(i,k)*factor
    1505. 

  1505.               nccol(i,k) = nccol(i,k)*factor
    1506. 

  1506.               nracw(i,k) = nracw(i,k)*factor
    1507. 

  1507.               naacw(i,k) = naacw(i,k)*factor
    1508. 

  1508.             endif
    1509. 

  1509. !
    1510. 

  1510. !     rain
    1511. 

  1511. !
    1512. 

  1512.             value = max(nrmin,ncr(i,k,3))
    1513. 

  1513.             source = (-nraut(i,k)+nrcol(i,k)+niacr(i,k)+nsacr(i,k)+ngacr(i,k)  &
    1514. 

  1514.                      )*dtcld
    1515. 

  1515.             if (source.gt.value) then
    1516. 

  1516.               factor = value/source
    1517. 

  1517.               nraut(i,k) = nraut(i,k)*factor
    1518. 

  1518.               nrcol(i,k) = nrcol(i,k)*factor
    1519. 

  1519.               niacr(i,k) = niacr(i,k)*factor
    1520. 

  1520.               nsacr(i,k) = nsacr(i,k)*factor
    1521. 

  1521.               ngacr(i,k) = ngacr(i,k)*factor
    1522. 

  1522.             endif
    1523. 

  1523. !
    1524. 

  1524.             work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k))
    1525. 

  1525. !     update
    1526. 

  1526.             q(i,k) = q(i,k)+work2(i,k)*dtcld
    1527. 

  1527.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
    1528. 

  1528.                            +paacw(i,k)+paacw(i,k))*dtcld,0.)
    1529. 

  1529.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
    1530. 

  1530.                            +prevp(i,k)-piacr(i,k)-pgacr(i,k)                   &
    1531. 

  1531.                            -psacr(i,k))*dtcld,0.)
    1532. 

  1532.             qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k)                 &
    1533. 

  1533.                            +psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k))       &
    1534. 

  1534.                            *dtcld,0.)
    1535. 

  1535.             qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k)      &
    1536. 

  1536.                            -pgaut(i,k)+piacr(i,k)*delta3                       &
    1537. 

  1537.                            +praci(i,k)*delta3+psaci(i,k)-pgacs(i,k)            &
    1538. 

  1538.                            -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2)          &
    1539. 

  1539.                            *dtcld,0.)
    1540. 

  1540.             qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k)                 &
    1541. 

  1541.                            +piacr(i,k)*(1.-delta3)                             &
    1542. 

  1542.                            +praci(i,k)*(1.-delta3)+psacr(i,k)*(1.-delta2)      &
    1543. 

  1543.                            +pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k)       &
    1544. 

  1544.                            +pgacr(i,k)+pgacs(i,k))*dtcld,0.)
    1545. 

  1545.             ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k)     &
    1546. 

  1546.                            -naacw(i,k)-naacw(i,k))*dtcld,0.)
    1547. 

  1547.             ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)-niacr(i,k)      &
    1548. 

  1548.                            -nsacr(i,k)-ngacr(i,k))*dtcld,0.)
    1549. 

  1549.             xlf = xls-xl(i,k)
    1550. 

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

  1551.                       -xl(i,k)*prevp(i,k)-xlf*(piacr(i,k)+paacw(i,k)           &
    1552. 

  1552.                       +paacw(i,k)+pgacr(i,k)+psacr(i,k))
    1553. 

  1553.             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
    1554. 

  1554.           else
    1555. 

  1555. !
    1556. 

  1556. !     cloud water
    1557. 

  1557. !
    1558. 

  1558.             value = max(qmin,qci(i,k,1))
    1559. 

  1559.             source= (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))              &
    1560. 

  1560.                    *dtcld
    1561. 

  1561.             if (source.gt.value) then
    1562. 

  1562.               factor = value/source
    1563. 

  1563.               praut(i,k) = praut(i,k)*factor
    1564. 

  1564.               pracw(i,k) = pracw(i,k)*factor
    1565. 

  1565.               paacw(i,k) = paacw(i,k)*factor
    1566. 

  1566.             endif
    1567. 

  1567. !
    1568. 

  1568. !     rain
    1569. 

  1569. !
    1570. 

  1570.             value = max(qmin,qrs(i,k,1))
    1571. 

  1571.             source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k)             &
    1572. 

  1572.                      -pracw(i,k)-paacw(i,k)-prevp(i,k))*dtcld
    1573. 

  1573.             if (source.gt.value) then
    1574. 

  1574.               factor = value/source
    1575. 

  1575.               praut(i,k) = praut(i,k)*factor
    1576. 

  1576.               prevp(i,k) = prevp(i,k)*factor
    1577. 

  1577.               pracw(i,k) = pracw(i,k)*factor
    1578. 

  1578.               paacw(i,k) = paacw(i,k)*factor
    1579. 

  1579.               pseml(i,k) = pseml(i,k)*factor
    1580. 

  1580.               pgeml(i,k) = pgeml(i,k)*factor
    1581. 

  1581.             endif
    1582. 

  1582. !
    1583. 

  1583. !     snow
    1584. 

  1584. !
    1585. 

  1585.             value = max(qcrmin,qrs(i,k,2))
    1586. 

  1586.             source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
    1587. 

  1587.             if (source.gt.value) then
    1588. 

  1588.               factor = value/source
    1589. 

  1589.               pgacs(i,k) = pgacs(i,k)*factor
    1590. 

  1590.               psevp(i,k) = psevp(i,k)*factor
    1591. 

  1591.               pseml(i,k) = pseml(i,k)*factor
    1592. 

  1592.             endif
    1593. 

  1593. !
    1594. 

  1594. !     graupel
    1595. 

  1595. !
    1596. 

  1596.             value = max(qcrmin,qrs(i,k,3))
    1597. 

  1597.             source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
    1598. 

  1598.             if (source.gt.value) then
    1599. 

  1599.               factor = value/source
    1600. 

  1600.               pgacs(i,k) = pgacs(i,k)*factor
    1601. 

  1601.               pgevp(i,k) = pgevp(i,k)*factor
    1602. 

  1602.               pgeml(i,k) = pgeml(i,k)*factor
    1603. 

  1603.             endif
    1604. 

  1604. !
    1605. 

  1605. !     cloud
    1606. 

  1606. !
    1607. 

  1607.             value = max(ncmin,ncr(i,k,2))
    1608. 

  1608.             source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+naacw(i,k)             &
    1609. 

  1609.                      +naacw(i,k))*dtcld
    1610. 

  1610.             if (source.gt.value) then
    1611. 

  1611.               factor = value/source
    1612. 

  1612.               nraut(i,k) = nraut(i,k)*factor
    1613. 

  1613.               nccol(i,k) = nccol(i,k)*factor
    1614. 

  1614.               nracw(i,k) = nracw(i,k)*factor
    1615. 

  1615.               naacw(i,k) = naacw(i,k)*factor
    1616. 

  1616.             endif
    1617. 

  1617. !
    1618. 

  1618. !     rain
    1619. 

  1619. !
    1620. 

  1620.             value = max(nrmin,ncr(i,k,3))
    1621. 

  1621.             source = (-nraut(i,k)+nrcol(i,k)-nseml(i,k)-ngeml(i,k)             &
    1622. 

  1622.                       )*dtcld
    1623. 

  1623.             if (source.gt.value) then
    1624. 

  1624.               factor = value/source
    1625. 

  1625.               nraut(i,k) = nraut(i,k)*factor
    1626. 

  1626.               nrcol(i,k) = nrcol(i,k)*factor
    1627. 

  1627.               nseml(i,k) = nseml(i,k)*factor
    1628. 

  1628.               ngeml(i,k) = ngeml(i,k)*factor
    1629. 

  1629.             endif
    1630. 

  1630. !
    1631. 

  1631.             work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
    1632. 

  1632. !     update
    1633. 

  1633.             q(i,k) = q(i,k)+work2(i,k)*dtcld
    1634. 

  1634.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
    1635. 

  1635.                     +paacw(i,k)+paacw(i,k))*dtcld,0.)
    1636. 

  1636.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
    1637. 

  1637.                     +prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k)               &
    1638. 

  1638.                     -pgeml(i,k))*dtcld,0.)
    1639. 

  1639.             qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k)                 &
    1640. 

  1640.                     +pseml(i,k))*dtcld,0.)
    1641. 

  1641.             qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k)                 &
    1642. 

  1642.                     +pgeml(i,k))*dtcld,0.)
    1643. 

  1643.             ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k)     &
    1644. 

  1644.                    -naacw(i,k)-naacw(i,k))*dtcld,0.)
    1645. 

  1645.             ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)+nseml(i,k)      &
    1646. 

  1646.                            +ngeml(i,k))*dtcld,0.)
    1647. 

  1647.             xlf = xls-xl(i,k)
    1648. 

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

  1649.                       -xlf*(pseml(i,k)+pgeml(i,k))
    1650. 

  1650.             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
    1651. 

  1651.           endif
    1652. 

  1652.         enddo
    1653. 

  1653.       enddo
    1654. 

  1654. !
    1655. 

  1655. ! Inline expansion for fpvs
    1656. 

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

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

  1658.       hsub = xls
    1659. 

  1659.       hvap = xlv0
    1660. 

  1660.       cvap = cpv
    1661. 

  1661.       ttp=t0c+0.01
    1662. 

  1662.       dldt=cvap-cliq
    1663. 

  1663.       xa=-dldt/rv
    1664. 

  1664.       xb=xa+hvap/(rv*ttp)
    1665. 

  1665.       dldti=cvap-cice
    1666. 

  1666.       xai=-dldti/rv
    1667. 

  1667.       xbi=xai+hsub/(rv*ttp)
    1668. 

  1668.       do k = kts, kte
    1669. 

  1669.         do i = its, ite
    1670. 

  1670.           tr=ttp/t(i,k)
    1671. 

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

  1672.           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
    1673. 

  1673.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    1674. 

  1674.           qs(i,k,1) = max(qs(i,k,1),qmin)
    1675. 

  1675.           tr=ttp/t(i,k)
    1676. 

  1676.           if(t(i,k).lt.ttp) then
    1677. 

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

  1678.           else
    1679. 

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

  1680.           endif
    1681. 

  1681.           qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
    1682. 

  1682.           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
    1683. 

  1683.           qs(i,k,2) = max(qs(i,k,2),qmin)
    1684. 

  1684.           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
    1685. 

  1685.         enddo
    1686. 

  1686.       enddo
    1687. 

  1687. !
    1688. 

  1688.       call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
    1689. 

  1689.                      rslope3,work1,workn,its,ite,kts,kte)
    1690. 

  1690.       do k = kts, kte
    1691. 

  1691.         do i = its, ite
    1692. 

  1692. !-----------------------------------------------------------------
    1693. 

  1693. ! re-compute the mean-volume drop diameter       [LH A10]
    1694. 

  1694. ! for raindrop distribution
    1695. 

  1695. !-----------------------------------------------------------------
    1696. 

  1696.           avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
    1697. 

  1697. !----------------------------------------------------------------
    1698. 

  1698. ! Nrevp_s: evaporation/condensation rate of rain   [LH A14]
    1699. 

  1699. !        (NR->NC)
    1700. 

  1700. !----------------------------------------------------------------
    1701. 

  1701.           if(avedia(i,k,2).le.di82) then
    1702. 

  1702.             ncr(i,k,2) = ncr(i,k,2)+ncr(i,k,3)
    1703. 

  1703.             ncr(i,k,3) = 0.
    1704. 

  1704. !----------------------------------------------------------------
    1705. 

  1705. ! Prevp_s: evaporation/condensation rate of rain [LH A15] [KK 23]
    1706. 

  1706. !        (QR->QC)
    1707. 

  1707. !----------------------------------------------------------------
    1708. 

  1708.             qci(i,k,1) = qci(i,k,1)+qrs(i,k,1)
    1709. 

  1709.             qrs(i,k,1) = 0.
    1710. 

  1710.           endif
    1711. 

  1711.         enddo
    1712. 

  1712.       enddo
    1713. 

  1713. !
    1714. 

  1714.       do k = kts, kte
    1715. 

  1715.         do i = its, ite
    1716. 

  1716. !---------------------------------------------------------------
    1717. 

  1717. ! rate of change of cloud drop concentration due to CCN activation
    1718. 

  1718. ! pcact: QV -> QC                                 [LH 8]  [KK 14]
    1719. 

  1719. ! ncact: NCCN -> NC                               [LH A2] [KK 12]
    1720. 

  1720. !---------------------------------------------------------------
    1721. 

  1721.           if(rh(i,k,1).gt.1.) then
    1722. 

  1722.             ncact(i,k) = max(0.,((ncr(i,k,1)+ncr(i,k,2))                       &
    1723. 

  1723.                        *min(1.,(rh(i,k,1)/satmax)**actk) - ncr(i,k,2)))/dtcld
    1724. 

  1724.             ncact(i,k) =min(ncact(i,k),max(ncr(i,k,1),0.)/dtcld)
    1725. 

  1725.             pcact(i,k) = min(4.*pi*denr*(actr*1.E-6)**3*ncact(i,k)/            &
    1726. 

  1726.                          (3.*den(i,k)),max(q(i,k),0.)/dtcld)
    1727. 

  1727.             q(i,k) = max(q(i,k)-pcact(i,k)*dtcld,0.)
    1728. 

  1728.             qci(i,k,1) = max(qci(i,k,1)+pcact(i,k)*dtcld,0.)
    1729. 

  1729.             ncr(i,k,1) = max(ncr(i,k,1)-ncact(i,k)*dtcld,0.)
    1730. 

  1730.             ncr(i,k,2) = max(ncr(i,k,2)+ncact(i,k)*dtcld,0.)
    1731. 

  1731.             t(i,k) = t(i,k)+pcact(i,k)*xl(i,k)/cpm(i,k)*dtcld
    1732. 

  1732.           endif  
    1733. 

  1733. !---------------------------------------------------------------
    1734. 

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

  1735. !     if there exists additional water vapor condensated/if
    1736. 

  1736. !     evaporation of cloud water is not enough to remove subsaturation
    1737. 

  1737. !  (QV->QC or QC->QV)     
    1738. 

  1738. !---------------------------------------------------------------
    1739. 

  1739.           tr=ttp/t(i,k)
    1740. 

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

  1741.           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
    1742. 

  1742.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    1743. 

  1743.           qs(i,k,1) = max(qs(i,k,1),qmin)
    1744. 

  1744.           work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
    1745. 

  1745.           work2(i,k) = qci(i,k,1)+work1(i,k,1)
    1746. 

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

  1747.           if(qci(i,k,1).gt.0. .and. work1(i,k,1).lt.0.)                        & 
    1748. 

  1748.             pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
    1749. 

  1749. !----------------------------------------------------------------
    1750. 

  1750. ! ncevp: evpration of Cloud number concentration  [LH A3] 
    1751. 

  1751. ! (NC->NCCN) 
    1752. 

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

  1753.           if(pcond(i,k).eq.-qci(i,k,1)/dtcld) then
    1754. 

  1754.             ncr(i,k,2) = 0.
    1755. 

  1755.             ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,2)
    1756. 

  1756.           endif
    1757. 

  1757. !
    1758. 

  1758.           q(i,k) = q(i,k)-pcond(i,k)*dtcld
    1759. 

  1759.           qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
    1760. 

  1760.           t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
    1761. 

  1761.         enddo
    1762. 

  1762.       enddo
    1763. 

  1763. !
    1764. 

  1764. !----------------------------------------------------------------
    1765. 

  1765. !     padding for small values
    1766. 

  1766. !
    1767. 

  1767.       do k = kts, kte
    1768. 

  1768.         do i = its, ite
    1769. 

  1769.           if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
    1770. 

  1770.           if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
    1771. 

  1771.         enddo
    1772. 

  1772.       enddo
    1773. 

  1773.       enddo                  ! big loops
    1774. 

  1774.   END SUBROUTINE wdm62d
    1775. 

  1775. ! ...................................................................
    1776. 

  1776.       REAL FUNCTION rgmma(x)
    1777. 

  1777. !-------------------------------------------------------------------
    1778. 

  1778.   IMPLICIT NONE
    1779. 

  1779. !-------------------------------------------------------------------
    1780. 

  1780. !     rgmma function:  use infinite product form
    1781. 

  1781.       REAL :: euler
    1782. 

  1782.       PARAMETER (euler=0.577215664901532)
    1783. 

  1783.       REAL :: x, y
    1784. 

  1784.       INTEGER :: i
    1785. 

  1785.       if(x.eq.1.)then
    1786. 

  1786.         rgmma=0.
    1787. 

  1787.           else
    1788. 

  1788.         rgmma=x*exp(euler*x)
    1789. 

  1789.         do i=1,10000
    1790. 

  1790.           y=float(i)
    1791. 

  1791.           rgmma=rgmma*(1.000+x/y)*exp(-x/y)
    1792. 

  1792.         enddo
    1793. 

  1793.         rgmma=1./rgmma
    1794. 

  1794.       endif
    1795. 

  1795.       END FUNCTION rgmma
    1796. 

  1796. !
    1797. 

  1797. !--------------------------------------------------------------------------
    1798. 

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

  1799. !--------------------------------------------------------------------------
    1800. 

  1800.       IMPLICIT NONE
    1801. 

  1801. !--------------------------------------------------------------------------
    1802. 

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

  1803.            xai,xbi,ttp,tr
    1804. 

  1804.       INTEGER ice
    1805. 

  1805. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1806. 

  1806.       ttp=t0c+0.01
    1807. 

  1807.       dldt=cvap-cliq
    1808. 

  1808.       xa=-dldt/rv
    1809. 

  1809.       xb=xa+hvap/(rv*ttp)
    1810. 

  1810.       dldti=cvap-cice
    1811. 

  1811.       xai=-dldti/rv
    1812. 

  1812.       xbi=xai+hsub/(rv*ttp)
    1813. 

  1813.       tr=ttp/t
    1814. 

  1814.       if(t.lt.ttp .and. ice.eq.1) then
    1815. 

  1815.         fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
    1816. 

  1816.       else
    1817. 

  1817.         fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
    1818. 

  1818.       endif
    1819. 

  1819. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1820. 

  1820.       END FUNCTION fpvs
    1821. 

  1821. !-------------------------------------------------------------------
    1822. 

  1822.   SUBROUTINE wdm6init(den0,denr,dens,cl,cpv, ccn0, allowed_to_read)
    1823. 

  1823. !-------------------------------------------------------------------
    1824. 

  1824.   IMPLICIT NONE
    1825. 

  1825. !-------------------------------------------------------------------
    1826. 

  1826. !.... constants which may not be tunable
    1827. 

  1827.    REAL, INTENT(IN) :: den0,denr,dens,cl,cpv,ccn0
    1828. 

  1828.    LOGICAL, INTENT(IN) :: allowed_to_read
    1829. 

  1829. !
    1830. 

  1830.    pi = 4.*atan(1.)
    1831. 

  1831.    xlv1 = cl-cpv
    1832. 

  1832. !
    1833. 

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

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

  1835.    pidnc = pi*denr/6.
    1836. 

  1836. !
    1837. 

  1837.    bvtr1 = 1.+bvtr
    1838. 

  1838.    bvtr2 = 2.+bvtr
    1839. 

  1839.    bvtr3 = 3.+bvtr
    1840. 

  1840.    bvtr4 = 4.+bvtr
    1841. 

  1841.    bvtr5 = 5.+bvtr
    1842. 

  1842.    bvtr6 = 6.+bvtr
    1843. 

  1843.    bvtr7 = 7.+bvtr
    1844. 

  1844.    bvtr2o5 = 2.5+.5*bvtr
    1845. 

  1845.    bvtr3o5 = 3.5+.5*bvtr
    1846. 

  1846.    g1pbr = rgmma(bvtr1)
    1847. 

  1847.    g2pbr = rgmma(bvtr2)
    1848. 

  1848.    g3pbr = rgmma(bvtr3)
    1849. 

  1849.    g4pbr = rgmma(bvtr4)            ! 17.837825
    1850. 

  1850.    g5pbr = rgmma(bvtr5)
    1851. 

  1851.    g6pbr = rgmma(bvtr6)
    1852. 

  1852.    g7pbr = rgmma(bvtr7)
    1853. 

  1853.    g5pbro2 = rgmma(bvtr2o5) 
    1854. 

  1854.    g7pbro2 = rgmma(bvtr3o5)
    1855. 

  1855.    pvtr = avtr*g5pbr/24.
    1856. 

  1856.    pvtrn = avtr*g2pbr
    1857. 

  1857.    eacrr = 1.0
    1858. 

  1858.    pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
    1859. 

  1859.    precr1 = 2.*pi*1.56
    1860. 

  1860.    precr2 = 2.*pi*.31*avtr**.5*g7pbro2
    1861. 

  1861.    pidn0r =  pi*denr*n0r
    1862. 

  1862.    pidnr = 4.*pi*denr
    1863. 

  1863. !
    1864. 

  1864.    xmmax = (dimax/dicon)**2
    1865. 

  1865.    roqimax = 2.08e22*dimax**8
    1866. 

  1866. !
    1867. 

  1867.    bvts1 = 1.+bvts
    1868. 

  1868.    bvts2 = 2.5+.5*bvts
    1869. 

  1869.    bvts3 = 3.+bvts
    1870. 

  1870.    bvts4 = 4.+bvts
    1871. 

  1871.    g1pbs = rgmma(bvts1)    !.8875
    1872. 

  1872.    g3pbs = rgmma(bvts3)
    1873. 

  1873.    g4pbs = rgmma(bvts4)    ! 12.0786
    1874. 

  1874.    g5pbso2 = rgmma(bvts2)
    1875. 

  1875.    pvts = avts*g4pbs/6.
    1876. 

  1876.    pacrs = pi*n0s*avts*g3pbs*.25
    1877. 

  1877.    precs1 = 4.*n0s*.65
    1878. 

  1878.    precs2 = 4.*n0s*.44*avts**.5*g5pbso2
    1879. 

  1879.    pidn0s =  pi*dens*n0s
    1880. 

  1880. !
    1881. 

  1881.    pacrc = pi*n0s*avts*g3pbs*.25*eacrc
    1882. 

  1882. !
    1883. 

  1883.    bvtg1 = 1.+bvtg
    1884. 

  1884.    bvtg2 = 2.5+.5*bvtg
    1885. 

  1885.    bvtg3 = 3.+bvtg
    1886. 

  1886.    bvtg4 = 4.+bvtg
    1887. 

  1887.    g1pbg = rgmma(bvtg1)
    1888. 

  1888.    g3pbg = rgmma(bvtg3)
    1889. 

  1889.    g4pbg = rgmma(bvtg4)
    1890. 

  1890.    g5pbgo2 = rgmma(bvtg2)
    1891. 

  1891.    pacrg = pi*n0g*avtg*g3pbg*.25
    1892. 

  1892.    pvtg = avtg*g4pbg/6.
    1893. 

  1893.    precg1 = 2.*pi*n0g*.78
    1894. 

  1894.    precg2 = 2.*pi*n0g*.31*avtg**.5*g5pbgo2
    1895. 

  1895.    pidn0g =  pi*deng*n0g
    1896. 

  1896. !
    1897. 

  1897.    rslopecmax = 1./lamdacmax
    1898. 

  1898.    rslopermax = 1./lamdarmax
    1899. 

  1899.    rslopesmax = 1./lamdasmax
    1900. 

  1900.    rslopegmax = 1./lamdagmax
    1901. 

  1901.    rsloperbmax = rslopermax ** bvtr
    1902. 

  1902.    rslopesbmax = rslopesmax ** bvts
    1903. 

  1903.    rslopegbmax = rslopegmax ** bvtg
    1904. 

  1904.    rslopec2max = rslopecmax * rslopecmax
    1905. 

  1905.    rsloper2max = rslopermax * rslopermax
    1906. 

  1906.    rslopes2max = rslopesmax * rslopesmax
    1907. 

  1907.    rslopeg2max = rslopegmax * rslopegmax
    1908. 

  1908.    rslopec3max = rslopec2max * rslopecmax
    1909. 

  1909.    rsloper3max = rsloper2max * rslopermax
    1910. 

  1910.    rslopes3max = rslopes2max * rslopesmax
    1911. 

  1911.    rslopeg3max = rslopeg2max * rslopegmax
    1912. 

  1912. !
    1913. 

  1913.   END SUBROUTINE wdm6init
    1914. 

  1914. !------------------------------------------------------------------------------
    1915. 

  1915.       subroutine slope_wdm6(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
    1916. 

  1916.                             vt,vtn,its,ite,kts,kte)
    1917. 

  1917.   IMPLICIT NONE
    1918. 

  1918.   INTEGER       ::               its,ite, jts,jte, kts,kte
    1919. 

  1919.   REAL, DIMENSION( its:ite , kts:kte,3) ::                                     &
    1920. 

  1920.                                                                           qrs, &
    1921. 

  1921.                                                                        rslope, &
    1922. 

  1922.                                                                       rslopeb, &
    1923. 

  1923.                                                                       rslope2, &
    1924. 

  1924.                                                                       rslope3, & 
    1925. 

  1925.                                                                            vt
    1926. 

  1926.   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
    1927. 

  1927.                                                                           ncr, & 
    1928. 

  1928.                                                                           vtn, & 
    1929. 

  1929.                                                                           den, &
    1930. 

  1930.                                                                        denfac, &
    1931. 

  1931.                                                                             t
    1932. 

  1932.   REAL, PARAMETER  :: t0c = 273.15
    1933. 

  1933.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    1934. 

  1934.                                                                        n0sfac
    1935. 

  1935.   REAL       ::  lamdar, lamdas, lamdag, x, y, z, supcol
    1936. 

  1936.   integer :: i, j, k
    1937. 

  1937. !----------------------------------------------------------------
    1938. 

  1938. !     size distributions: (x=mixing ratio, y=air density):
    1939. 

  1939. !     valid for mixing ratio > 1.e-9 kg/kg.
    1940. 

  1940. !
    1941. 

  1941. !  Optimizatin : A**B => exp(log(A)*(B))
    1942. 

  1942.       lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
    1943. 

  1943.       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
    1944. 

  1944.       lamdag(x,y)=   sqrt(sqrt(pidn0g/(x*y)))      ! (pidn0g/(x*y))**.25
    1945. 

  1945. !
    1946. 

  1946.       do k = kts, kte
    1947. 

  1947.         do i = its, ite
    1948. 

  1948.           supcol = t0c-t(i,k)
    1949. 

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

  1950. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
    1951. 

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

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

  1953.           if(qrs(i,k,1).le.qcrmin .or. ncr(i,k).le.nrmin ) then
    1954. 

  1954.             rslope(i,k,1) = rslopermax
    1955. 

  1955.             rslopeb(i,k,1) = rsloperbmax
    1956. 

  1956.             rslope2(i,k,1) = rsloper2max
    1957. 

  1957.             rslope3(i,k,1) = rsloper3max
    1958. 

  1958.           else
    1959. 

  1959.             rslope(i,k,1) = min(1./lamdar(qrs(i,k,1),den(i,k),ncr(i,k)),1.e-3)
    1960. 

  1960.             rslopeb(i,k,1) = rslope(i,k,1)**bvtr
    1961. 

  1961.             rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
    1962. 

  1962.             rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
    1963. 

  1963.           endif
    1964. 

  1964.           if(qrs(i,k,2).le.qcrmin) then
    1965. 

  1965.             rslope(i,k,2) = rslopesmax
    1966. 

  1966.             rslopeb(i,k,2) = rslopesbmax
    1967. 

  1967.             rslope2(i,k,2) = rslopes2max
    1968. 

  1968.             rslope3(i,k,2) = rslopes3max
    1969. 

  1969.           else
    1970. 

  1970.             rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
    1971. 

  1971.             rslopeb(i,k,2) = rslope(i,k,2)**bvts
    1972. 

  1972.             rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
    1973. 

  1973.             rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
    1974. 

  1974.           endif
    1975. 

  1975.           if(qrs(i,k,3).le.qcrmin) then
    1976. 

  1976.             rslope(i,k,3) = rslopegmax
    1977. 

  1977.             rslopeb(i,k,3) = rslopegbmax
    1978. 

  1978.             rslope2(i,k,3) = rslopeg2max
    1979. 

  1979.             rslope3(i,k,3) = rslopeg3max
    1980. 

  1980.           else
    1981. 

  1981.             rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
    1982. 

  1982.             rslopeb(i,k,3) = rslope(i,k,3)**bvtg
    1983. 

  1983.             rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
    1984. 

  1984.             rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
    1985. 

  1985.           endif
    1986. 

  1986.           vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
    1987. 

  1987.           vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
    1988. 

  1988.           vt(i,k,3) = pvtg*rslopeb(i,k,3)*denfac(i,k)
    1989. 

  1989.           vtn(i,k) = pvtrn*rslopeb(i,k,1)*denfac(i,k)
    1990. 

  1990.           if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0 
    1991. 

  1991.           if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
    1992. 

  1992.           if(qrs(i,k,3).le.0.0) vt(i,k,3) = 0.0
    1993. 

  1993.           if(ncr(i,k).le.0.0) vtn(i,k) = 0.0 
    1994. 

  1994.         enddo
    1995. 

  1995.       enddo
    1996. 

  1996.   END subroutine slope_wdm6
    1997. 

  1997. !-----------------------------------------------------------------------------
    1998. 

  1998.       subroutine slope_rain(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
    1999. 

  1999.                             vt,vtn,its,ite,kts,kte)
    2000. 

  2000.   IMPLICIT NONE
    2001. 

  2001.   INTEGER       ::               its,ite, jts,jte, kts,kte
    2002. 

  2002.   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
    2003. 

  2003.                                                                           qrs, &
    2004. 

  2004.                                                                           ncr, & 
    2005. 

  2005.                                                                        rslope, &
    2006. 

  2006.                                                                       rslopeb, &
    2007. 

  2007.                                                                       rslope2, &
    2008. 

  2008.                                                                       rslope3, &
    2009. 

  2009.                                                                            vt, &
    2010. 

  2010.                                                                           vtn, &
    2011. 

  2011.                                                                           den, &
    2012. 

  2012.                                                                        denfac, &
    2013. 

  2013.                                                                             t
    2014. 

  2014.   REAL, PARAMETER  :: t0c = 273.15
    2015. 

  2015.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    2016. 

  2016.                                                                        n0sfac
    2017. 

  2017.   REAL       ::  lamdar, x, y, z, supcol
    2018. 

  2018.   integer :: i, j, k
    2019. 

  2019. !----------------------------------------------------------------
    2020. 

  2020. !     size distributions: (x=mixing ratio, y=air density):
    2021. 

  2021. !     valid for mixing ratio > 1.e-9 kg/kg.
    2022. 

  2022.       lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
    2023. 

  2023. !
    2024. 

  2024.       do k = kts, kte
    2025. 

  2025.         do i = its, ite
    2026. 

  2026.           if(qrs(i,k).le.qcrmin .or. ncr(i,k).le.nrmin) then
    2027. 

  2027.             rslope(i,k) = rslopermax
    2028. 

  2028.             rslopeb(i,k) = rsloperbmax
    2029. 

  2029.             rslope2(i,k) = rsloper2max
    2030. 

  2030.             rslope3(i,k) = rsloper3max
    2031. 

  2031.           else
    2032. 

  2032.             rslope(i,k) = min(1./lamdar(qrs(i,k),den(i,k),ncr(i,k)),1.e-3)
    2033. 

  2033.             rslopeb(i,k) = rslope(i,k)**bvtr
    2034. 

  2034.             rslope2(i,k) = rslope(i,k)*rslope(i,k)
    2035. 

  2035.             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
    2036. 

  2036.           endif
    2037. 

  2037.           vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
    2038. 

  2038.           vtn(i,k) = pvtrn*rslopeb(i,k)*denfac(i,k)
    2039. 

  2039.           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
    2040. 

  2040.           if(ncr(i,k).le.0.0) vtn(i,k) = 0.0 
    2041. 

  2041.         enddo
    2042. 

  2042.       enddo
    2043. 

  2043.   END subroutine slope_rain
    2044. 

  2044. !------------------------------------------------------------------------------
    2045. 

  2045.       subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
    2046. 

  2046.                             vt,its,ite,kts,kte)
    2047. 

  2047.   IMPLICIT NONE
    2048. 

  2048.   INTEGER       ::               its,ite, jts,jte, kts,kte
    2049. 

  2049.   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
    2050. 

  2050.                                                                           qrs, &
    2051. 

  2051.                                                                        rslope, &
    2052. 

  2052.                                                                       rslopeb, &
    2053. 

  2053.                                                                       rslope2, &
    2054. 

  2054.                                                                       rslope3, &
    2055. 

  2055.                                                                            vt, &
    2056. 

  2056.                                                                           den, &
    2057. 

  2057.                                                                        denfac, &
    2058. 

  2058.                                                                             t
    2059. 

  2059.   REAL, PARAMETER  :: t0c = 273.15
    2060. 

  2060.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    2061. 

  2061.                                                                        n0sfac
    2062. 

  2062.   REAL       ::  lamdas, x, y, z, supcol
    2063. 

  2063.   integer :: i, j, k
    2064. 

  2064. !----------------------------------------------------------------
    2065. 

  2065. !     size distributions: (x=mixing ratio, y=air density):
    2066. 

  2066. !     valid for mixing ratio > 1.e-9 kg/kg.
    2067. 

  2067.       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
    2068. 

  2068. !
    2069. 

  2069.       do k = kts, kte
    2070. 

  2070.         do i = its, ite
    2071. 

  2071.           supcol = t0c-t(i,k)
    2072. 

  2072. !---------------------------------------------------------------
    2073. 

  2073. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
    2074. 

  2074. !---------------------------------------------------------------
    2075. 

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

  2076.           if(qrs(i,k).le.qcrmin)then
    2077. 

  2077.             rslope(i,k) = rslopesmax
    2078. 

  2078.             rslopeb(i,k) = rslopesbmax
    2079. 

  2079.             rslope2(i,k) = rslopes2max
    2080. 

  2080.             rslope3(i,k) = rslopes3max
    2081. 

  2081.           else
    2082. 

  2082.             rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
    2083. 

  2083.             rslopeb(i,k) = rslope(i,k)**bvts
    2084. 

  2084.             rslope2(i,k) = rslope(i,k)*rslope(i,k)
    2085. 

  2085.             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
    2086. 

  2086.           endif
    2087. 

  2087.           vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
    2088. 

  2088.           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
    2089. 

  2089.         enddo
    2090. 

  2090.       enddo
    2091. 

  2091.   END subroutine slope_snow
    2092. 

  2092. !----------------------------------------------------------------------------------
    2093. 

  2093.       subroutine slope_graup(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
    2094. 

  2094.                             vt,its,ite,kts,kte)
    2095. 

  2095.   IMPLICIT NONE
    2096. 

  2096.   INTEGER       ::               its,ite, jts,jte, kts,kte
    2097. 

  2097.   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
    2098. 

  2098.                                                                           qrs, &
    2099. 

  2099.                                                                        rslope, &
    2100. 

  2100.                                                                       rslopeb, &
    2101. 

  2101.                                                                       rslope2, &
    2102. 

  2102.                                                                       rslope3, &
    2103. 

  2103.                                                                            vt, &
    2104. 

  2104.                                                                           den, &
    2105. 

  2105.                                                                        denfac, &
    2106. 

  2106.                                                                             t
    2107. 

  2107.   REAL, PARAMETER  :: t0c = 273.15
    2108. 

  2108.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    2109. 

  2109.                                                                        n0sfac
    2110. 

  2110.   REAL       ::  lamdag, x, y, z, supcol
    2111. 

  2111.   integer :: i, j, k
    2112. 

  2112. !----------------------------------------------------------------
    2113. 

  2113. !     size distributions: (x=mixing ratio, y=air density):
    2114. 

  2114. !     valid for mixing ratio > 1.e-9 kg/kg.
    2115. 

  2115.       lamdag(x,y)=   sqrt(sqrt(pidn0g/(x*y)))      ! (pidn0g/(x*y))**.25
    2116. 

  2116. !
    2117. 

  2117.       do k = kts, kte
    2118. 

  2118.         do i = its, ite
    2119. 

  2119. !---------------------------------------------------------------
    2120. 

  2120. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
    2121. 

  2121. !---------------------------------------------------------------
    2122. 

  2122.           if(qrs(i,k).le.qcrmin)then
    2123. 

  2123.             rslope(i,k) = rslopegmax
    2124. 

  2124.             rslopeb(i,k) = rslopegbmax
    2125. 

  2125.             rslope2(i,k) = rslopeg2max
    2126. 

  2126.             rslope3(i,k) = rslopeg3max
    2127. 

  2127.           else
    2128. 

  2128.             rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k))
    2129. 

  2129.             rslopeb(i,k) = rslope(i,k)**bvtg
    2130. 

  2130.             rslope2(i,k) = rslope(i,k)*rslope(i,k)
    2131. 

  2131.             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
    2132. 

  2132.           endif
    2133. 

  2133.           vt(i,k) = pvtg*rslopeb(i,k)*denfac(i,k)
    2134. 

  2134.           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
    2135. 

  2135.         enddo
    2136. 

  2136.       enddo
    2137. 

  2137.   END subroutine slope_graup
    2138. 

  2138. !---------------------------------------------------------------------------------
    2139. 

  2139. !-------------------------------------------------------------------
    2140. 

  2140.       SUBROUTINE nislfv_rain_plmr(im,km,denl,denfacl,tkl,dzl,wwl,rql,rnl,precip,dt,id,iter,rid)
    2141. 

  2141. !-------------------------------------------------------------------
    2142. 

  2142. !
    2143. 

  2143. ! for non-iteration semi-Lagrangain forward advection for cloud
    2144. 

  2144. ! with mass conservation and positive definite advection
    2145. 

  2145. ! 2nd order interpolation with monotonic piecewise linear method
    2146. 

  2146. ! this routine is under assumption of decfl < 1 for semi_Lagrangian
    2147. 

  2147. !
    2148. 

  2148. ! dzl    depth of model layer in meter
    2149. 

  2149. ! wwl    terminal velocity at model layer m/s
    2150. 

  2150. ! rql    cloud density*mixing ration
    2151. 

  2151. ! precip precipitation
    2152. 

  2152. ! dt     time step
    2153. 

  2153. ! id     kind of precip: 0 test case; 1 raindrop
    2154. 

  2154. ! iter   how many time to guess mean terminal velocity: 0 pure forward.
    2155. 

  2155. !        0 : use departure wind for advection
    2156. 

  2156. !        1 : use mean wind for advection
    2157. 

  2157. !        > 1 : use mean wind after iter-1 iterations
    2158. 

  2158. !        rid : 1 for number 0 for mixing ratio
    2159. 

  2159. !
    2160. 

  2160. ! author: hann-ming henry juang <henry.juang@noaa.gov>
    2161. 

  2161. !         implemented by song-you hong
    2162. 

  2162. !
    2163. 

  2163.       implicit none
    2164. 

  2164.       integer  im,km,id
    2165. 

  2165.       real  dt
    2166. 

  2166.       real  dzl(im,km),wwl(im,km),rql(im,km),rnl(im,km),precip(im)
    2167. 

  2167.       real  denl(im,km),denfacl(im,km),tkl(im,km)
    2168. 

  2168. !
    2169. 

  2169.       integer  i,k,n,m,kk,kb,kt,iter,rid
    2170. 

  2170.       real  tl,tl2,qql,dql,qqd
    2171. 

  2171.       real  th,th2,qqh,dqh
    2172. 

  2172.       real  zsum,qsum,dim,dip,c1,con1,fa1,fa2
    2173. 

  2173.       real  allold, allnew, zz, dzamin, cflmax, decfl
    2174. 

  2174.       real  dz(km), ww(km), qq(km), nr(km), wd(km), wa(km), wa2(km), was(km)
    2175. 

  2175.       real  den(km), denfac(km), tk(km)
    2176. 

  2176.       real  wi(km+1), zi(km+1), za(km+1)
    2177. 

  2177.       real  qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
    2178. 

  2178.       real  dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
    2179. 

  2179. !
    2180. 

  2180.       precip(:) = 0.0
    2181. 

  2181. !
    2182. 

  2182.       i_loop : do i=1,im
    2183. 

  2183. ! -----------------------------------
    2184. 

  2184.       dz(:) = dzl(i,:)
    2185. 

  2185.       qq(:) = rql(i,:)
    2186. 

  2186.       nr(:) = rnl(i,:)
    2187. 

  2187.       if(rid .eq. 1) nr(:) = rnl(i,:)/denl(i,:)
    2188. 

  2188.       ww(:) = wwl(i,:)
    2189. 

  2189.       den(:) = denl(i,:)
    2190. 

  2190.       denfac(:) = denfacl(i,:)
    2191. 

  2191.       tk(:) = tkl(i,:)
    2192. 

  2192. ! skip for no precipitation for all layers
    2193. 

  2193.       allold = 0.0
    2194. 

  2194.       do k=1,km
    2195. 

  2195.         allold = allold + qq(k)
    2196. 

  2196.       enddo
    2197. 

  2197.       if(allold.le.0.0) then
    2198. 

  2198.         cycle i_loop
    2199. 

  2199.       endif
    2200. 

  2200. !
    2201. 

  2201. ! compute interface values
    2202. 

  2202.       zi(1)=0.0
    2203. 

  2203.       do k=1,km
    2204. 

  2204.         zi(k+1) = zi(k)+dz(k)
    2205. 

  2205.       enddo
    2206. 

  2206. !
    2207. 

  2207. ! save departure wind
    2208. 

  2208.       wd(:) = ww(:)
    2209. 

  2209.       n=1
    2210. 

  2210.  100  continue
    2211. 

  2211. ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
    2212. 

  2212. ! 2nd order interpolation to get wi
    2213. 

  2213.       wi(1) = ww(1)
    2214. 

  2214.       wi(km+1) = ww(km)
    2215. 

  2215.       do k=2,km
    2216. 

  2216.         wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
    2217. 

  2217.       enddo
    2218. 

  2218. ! 3rd order interpolation to get wi
    2219. 

  2219.       fa1 = 9./16.
    2220. 

  2220.       fa2 = 1./16.
    2221. 

  2221.       wi(1) = ww(1)
    2222. 

  2222.       wi(2) = 0.5*(ww(2)+ww(1))
    2223. 

  2223.       do k=3,km-1
    2224. 

  2224.         wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
    2225. 

  2225.       enddo
    2226. 

  2226.       wi(km) = 0.5*(ww(km)+ww(km-1))
    2227. 

  2227.       wi(km+1) = ww(km)
    2228. 

  2228. !
    2229. 

  2229. ! terminate of top of raingroup
    2230. 

  2230.       do k=2,km
    2231. 

  2231.         if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
    2232. 

  2232.       enddo
    2233. 

  2233. !
    2234. 

  2234. ! diffusivity of wi
    2235. 

  2235.       con1 = 0.05
    2236. 

  2236.       do k=km,1,-1
    2237. 

  2237.         decfl = (wi(k+1)-wi(k))*dt/dz(k)
    2238. 

  2238.         if( decfl .gt. con1 ) then
    2239. 

  2239.           wi(k) = wi(k+1) - con1*dz(k)/dt
    2240. 

  2240.         endif
    2241. 

  2241.       enddo
    2242. 

  2242. ! compute arrival point
    2243. 

  2243.       do k=1,km+1
    2244. 

  2244.         za(k) = zi(k) - wi(k)*dt
    2245. 

  2245.       enddo
    2246. 

  2246. !
    2247. 

  2247.       do k=1,km
    2248. 

  2248.         dza(k) = za(k+1)-za(k)
    2249. 

  2249.       enddo
    2250. 

  2250.       dza(km+1) = zi(km+1) - za(km+1)
    2251. 

  2251. !     
    2252. 

  2252. ! computer deformation at arrival point
    2253. 

  2253.       do k=1,km
    2254. 

  2254.         qa(k) = qq(k)*dz(k)/dza(k)
    2255. 

  2255.         qr(k) = qa(k)/den(k)
    2256. 

  2256.         if(rid .eq. 1) qr(k) = qa(K)
    2257. 

  2257.       enddo
    2258. 

  2258.       qa(km+1) = 0.0
    2259. 

  2259. !     call maxmin(km,1,qa,' arrival points ')
    2260. 

  2260. !     
    2261. 

  2261. ! compute arrival terminal velocity, and estimate mean terminal velocity
    2262. 

  2262. ! then back to use mean terminal velocity
    2263. 

  2263.       if( n.le.iter ) then
    2264. 

  2264.         if(rid.eq.1) then
    2265. 

  2265.         call slope_rain(nr,qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,wa2,1,1,1,km)
    2266. 

  2266.         else
    2267. 

  2267.         call slope_rain(qr,nr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,wa2,1,1,1,km)
    2268. 

  2268.         endif
    2269. 

  2269.         if(rid.eq.1) wa(:) = wa2(:)
    2270. 

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

  2271.         do k=1,km
    2272. 

  2272. !#ifdef DEBUG 
    2273. 

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

  2274. !#endif
    2275. 

  2275. ! mean wind is average of departure and new arrival winds
    2276. 

  2276.           ww(k) = 0.5* ( wd(k)+wa(k) )
    2277. 

  2277.         enddo
    2278. 

  2278.         was(:) = wa(:)
    2279. 

  2279.         n=n+1
    2280. 

  2280.         go to 100
    2281. 

  2281.       endif
    2282. 

  2282. !     
    2283. 

  2283. ! estimate values at arrival cell interface with monotone
    2284. 

  2284.       do k=2,km
    2285. 

  2285.         dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
    2286. 

  2286.         dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
    2287. 

  2287.         if( dip*dim.le.0.0 ) then
    2288. 

  2288.           qmi(k)=qa(k)
    2289. 

  2289.           qpi(k)=qa(k)
    2290. 

  2290.         else
    2291. 

  2291.           qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
    2292. 

  2292.           qmi(k)=2.0*qa(k)-qpi(k)
    2293. 

  2293.           if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
    2294. 

  2294.             qpi(k) = qa(k)
    2295. 

  2295.             qmi(k) = qa(k)
    2296. 

  2296.           endif
    2297. 

  2297.         endif
    2298. 

  2298.       enddo
    2299. 

  2299.       qpi(1)=qa(1)
    2300. 

  2300.       qmi(1)=qa(1)
    2301. 

  2301.       qmi(km+1)=qa(km+1)
    2302. 

  2302.       qpi(km+1)=qa(km+1)
    2303. 

  2303. !
    2304. 

  2304. ! interpolation to regular point
    2305. 

  2305.       qn = 0.0
    2306. 

  2306.       kb=1
    2307. 

  2307.       kt=1
    2308. 

  2308.       intp : do k=1,km
    2309. 

  2309.              kb=max(kb-1,1)
    2310. 

  2310.              kt=max(kt-1,1)
    2311. 

  2311. ! find kb and kt
    2312. 

  2312.              if( zi(k).ge.za(km+1) ) then
    2313. 

  2313.                exit intp
    2314. 

  2314.              else
    2315. 

  2315.                find_kb : do kk=kb,km
    2316. 

  2316.                          if( zi(k).le.za(kk+1) ) then
    2317. 

  2317.                            kb = kk
    2318. 

  2318.                            exit find_kb
    2319. 

  2319.                          else
    2320. 

  2320.                            cycle find_kb
    2321. 

  2321.                          endif
    2322. 

  2322.                enddo find_kb
    2323. 

  2323.                find_kt : do kk=kt,km
    2324. 

  2324.                          if( zi(k+1).le.za(kk) ) then
    2325. 

  2325.                            kt = kk
    2326. 

  2326.                            exit find_kt
    2327. 

  2327.                          else
    2328. 

  2328.                            cycle find_kt
    2329. 

  2329.                          endif
    2330. 

  2330.                enddo find_kt
    2331. 

  2331.                kt = kt - 1
    2332. 

  2332. ! compute q with piecewise constant method
    2333. 

  2333.                if( kt.eq.kb ) then
    2334. 

  2334.                  tl=(zi(k)-za(kb))/dza(kb)
    2335. 

  2335.                  th=(zi(k+1)-za(kb))/dza(kb)
    2336. 

  2336.                  tl2=tl*tl
    2337. 

  2337.                  th2=th*th
    2338. 

  2338.                  qqd=0.5*(qpi(kb)-qmi(kb))
    2339. 

  2339.                  qqh=qqd*th2+qmi(kb)*th
    2340. 

  2340.                  qql=qqd*tl2+qmi(kb)*tl
    2341. 

  2341.                  qn(k) = (qqh-qql)/(th-tl)
    2342. 

  2342.                else if( kt.gt.kb ) then
    2343. 

  2343.                  tl=(zi(k)-za(kb))/dza(kb)
    2344. 

  2344.                  tl2=tl*tl
    2345. 

  2345.                  qqd=0.5*(qpi(kb)-qmi(kb))
    2346. 

  2346.                  qql=qqd*tl2+qmi(kb)*tl
    2347. 

  2347.                  dql = qa(kb)-qql
    2348. 

  2348.                  zsum  = (1.-tl)*dza(kb)
    2349. 

  2349.                  qsum  = dql*dza(kb)
    2350. 

  2350.                  if( kt-kb.gt.1 ) then
    2351. 

  2351.                  do m=kb+1,kt-1
    2352. 

  2352.                    zsum = zsum + dza(m)
    2353. 

  2353.                    qsum = qsum + qa(m) * dza(m)
    2354. 

  2354.                  enddo
    2355. 

  2355.                  endif
    2356. 

  2356.                  th=(zi(k+1)-za(kt))/dza(kt)
    2357. 

  2357.                  th2=th*th
    2358. 

  2358.                  qqd=0.5*(qpi(kt)-qmi(kt))
    2359. 

  2359.                  dqh=qqd*th2+qmi(kt)*th
    2360. 

  2360.                  zsum  = zsum + th*dza(kt)
    2361. 

  2361.                  qsum  = qsum + dqh*dza(kt)
    2362. 

  2362.                  qn(k) = qsum/zsum
    2363. 

  2363.                endif
    2364. 

  2364.                cycle intp
    2365. 

  2365.              endif
    2366. 

  2366. !     
    2367. 

  2367.        enddo intp
    2368. 

  2368. !            
    2369. 

  2369. ! rain out
    2370. 

  2370.       sum_precip: do k=1,km
    2371. 

  2371.                     if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
    2372. 

  2372.                       precip(i) = precip(i) + qa(k)*dza(k)
    2373. 

  2373.                       cycle sum_precip
    2374. 

  2374.                     else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
    2375. 

  2375.                       precip(i) = precip(i) + qa(k)*(0.0-za(k))
    2376. 

  2376.                       exit sum_precip
    2377. 

  2377.                     endif
    2378. 

  2378.                     exit sum_precip
    2379. 

  2379.       enddo sum_precip
    2380. 

  2380. !              
    2381. 

  2381. ! replace the new values 
    2382. 

  2382.       rql(i,:) = qn(:)     
    2383. 

  2383. !                          
    2384. 

  2384. ! ----------------------------------
    2385. 

  2385.       enddo i_loop         
    2386. 

  2386. !                        
    2387. 

  2387.   END SUBROUTINE nislfv_rain_plmr
    2388. 

  2388. !-------------------------------------------------------------------
    2389. 

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

  2390. !-------------------------------------------------------------------
    2391. 

  2391. !
    2392. 

  2392. ! for non-iteration semi-Lagrangain forward advection for cloud
    2393. 

  2393. ! with mass conservation and positive definite advection
    2394. 

  2394. ! 2nd order interpolation with monotonic piecewise linear method
    2395. 

  2395. ! this routine is under assumption of decfl < 1 for semi_Lagrangian
    2396. 

  2396. !
    2397. 

  2397. ! dzl    depth of model layer in meter
    2398. 

  2398. ! wwl    terminal velocity at model layer m/s
    2399. 

  2399. ! rql    cloud density*mixing ration
    2400. 

  2400. ! precip precipitation
    2401. 

  2401. ! dt     time step
    2402. 

  2402. ! id     kind of precip: 0 test case; 1 raindrop
    2403. 

  2403. ! iter   how many time to guess mean terminal velocity: 0 pure forward.
    2404. 

  2404. !        0 : use departure wind for advection
    2405. 

  2405. !        1 : use mean wind for advection
    2406. 

  2406. !        > 1 : use mean wind after iter-1 iterations
    2407. 

  2407. !
    2408. 

  2408. ! author: hann-ming henry juang <henry.juang@noaa.gov>
    2409. 

  2409. !         implemented by song-you hong
    2410. 

  2410. !
    2411. 

  2411.       implicit none
    2412. 

  2412.       integer  im,km,id
    2413. 

  2413.       real  dt
    2414. 

  2414.       real  dzl(im,km),wwl(im,km),rql(im,km),rql2(im,km),precip(im),precip1(im),precip2(im)
    2415. 

  2415.       real  denl(im,km),denfacl(im,km),tkl(im,km)
    2416. 

  2416. !
    2417. 

  2417.       integer  i,k,n,m,kk,kb,kt,iter,ist
    2418. 

  2418.       real  tl,tl2,qql,dql,qqd
    2419. 

  2419.       real  th,th2,qqh,dqh
    2420. 

  2420.       real  zsum,qsum,dim,dip,c1,con1,fa1,fa2
    2421. 

  2421.       real  allold, allnew, zz, dzamin, cflmax, decfl
    2422. 

  2422.       real  dz(km), ww(km), qq(km), qq2(km), wd(km), wa(km), wa2(km), was(km)
    2423. 

  2423.       real  den(km), denfac(km), tk(km)
    2424. 

  2424.       real  wi(km+1), zi(km+1), za(km+1)
    2425. 

  2425.       real  qn(km), qr(km),qr2(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
    2426. 

  2426.       real  dza(km+1), qa(km+1), qa2(km+1),qmi(km+1), qpi(km+1)
    2427. 

  2427. !
    2428. 

  2428.       precip(:) = 0.0
    2429. 

  2429.       precip1(:) = 0.0
    2430. 

  2430.       precip2(:) = 0.0
    2431. 

  2431. !
    2432. 

  2432.       i_loop : do i=1,im
    2433. 

  2433. ! -----------------------------------
    2434. 

  2434.       dz(:) = dzl(i,:)
    2435. 

  2435.       qq(:) = rql(i,:)
    2436. 

  2436.       qq2(:) = rql2(i,:)
    2437. 

  2437.       ww(:) = wwl(i,:)
    2438. 

  2438.       den(:) = denl(i,:)
    2439. 

  2439.       denfac(:) = denfacl(i,:)
    2440. 

  2440.       tk(:) = tkl(i,:)
    2441. 

  2441. ! skip for no precipitation for all layers
    2442. 

  2442.       allold = 0.0
    2443. 

  2443.       do k=1,km
    2444. 

  2444.         allold = allold + qq(k)
    2445. 

  2445.       enddo
    2446. 

  2446.       if(allold.le.0.0) then
    2447. 

  2447.         cycle i_loop
    2448. 

  2448.       endif
    2449. 

  2449. !
    2450. 

  2450. ! compute interface values
    2451. 

  2451.       zi(1)=0.0
    2452. 

  2452.       do k=1,km
    2453. 

  2453.         zi(k+1) = zi(k)+dz(k)
    2454. 

  2454.       enddo
    2455. 

  2455. !
    2456. 

  2456. ! save departure wind
    2457. 

  2457.       wd(:) = ww(:)
    2458. 

  2458.       n=1
    2459. 

  2459.  100  continue
    2460. 

  2460. ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
    2461. 

  2461. ! 2nd order interpolation to get wi
    2462. 

  2462.       wi(1) = ww(1)
    2463. 

  2463.       wi(km+1) = ww(km)
    2464. 

  2464.       do k=2,km
    2465. 

  2465.         wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
    2466. 

  2466.       enddo
    2467. 

  2467. ! 3rd order interpolation to get wi
    2468. 

  2468.       fa1 = 9./16.
    2469. 

  2469.       fa2 = 1./16.
    2470. 

  2470.       wi(1) = ww(1)
    2471. 

  2471.       wi(2) = 0.5*(ww(2)+ww(1))
    2472. 

  2472.       do k=3,km-1
    2473. 

  2473.         wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
    2474. 

  2474.       enddo
    2475. 

  2475.       wi(km) = 0.5*(ww(km)+ww(km-1))
    2476. 

  2476.       wi(km+1) = ww(km)
    2477. 

  2477. !
    2478. 

  2478. ! terminate of top of raingroup
    2479. 

  2479.       do k=2,km
    2480. 

  2480.         if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
    2481. 

  2481.       enddo
    2482. 

  2482. !
    2483. 

  2483. ! diffusivity of wi
    2484. 

  2484.       con1 = 0.05
    2485. 

  2485.       do k=km,1,-1
    2486. 

  2486.         decfl = (wi(k+1)-wi(k))*dt/dz(k)
    2487. 

  2487.         if( decfl .gt. con1 ) then
    2488. 

  2488.           wi(k) = wi(k+1) - con1*dz(k)/dt
    2489. 

  2489.         endif
    2490. 

  2490.       enddo
    2491. 

  2491. ! compute arrival point
    2492. 

  2492.       do k=1,km+1
    2493. 

  2493.         za(k) = zi(k) - wi(k)*dt
    2494. 

  2494.       enddo
    2495. 

  2495. !
    2496. 

  2496.       do k=1,km
    2497. 

  2497.         dza(k) = za(k+1)-za(k)
    2498. 

  2498.       enddo
    2499. 

  2499.       dza(km+1) = zi(km+1) - za(km+1)
    2500. 

  2500. !
    2501. 

  2501. ! computer deformation at arrival point
    2502. 

  2502.       do k=1,km
    2503. 

  2503.         qa(k) = qq(k)*dz(k)/dza(k)
    2504. 

  2504.         qa2(k) = qq2(k)*dz(k)/dza(k)
    2505. 

  2505.         qr(k) = qa(k)/den(k)
    2506. 

  2506.         qr2(k) = qa2(k)/den(k)
    2507. 

  2507.       enddo
    2508. 

  2508.       qa(km+1) = 0.0
    2509. 

  2509.       qa2(km+1) = 0.0
    2510. 

  2510. !     call maxmin(km,1,qa,' arrival points ')
    2511. 

  2511. !
    2512. 

  2512. ! compute arrival terminal velocity, and estimate mean terminal velocity
    2513. 

  2513. ! then back to use mean terminal velocity
    2514. 

  2514.       if( n.le.iter ) then
    2515. 

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

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

  2517.         do k = 1, km
    2518. 

  2518.           tmp(k) = max((qr(k)+qr2(k)), 1.E-15)
    2519. 

  2519.           IF ( tmp(k) .gt. 1.e-15 ) THEN
    2520. 

  2520.             wa(k) = (wa(k)*qr(k) + wa2(k)*qr2(k))/tmp(k)
    2521. 

  2521.           ELSE
    2522. 

  2522.             wa(k) = 0.
    2523. 

  2523.           ENDIF
    2524. 

  2524.         enddo
    2525. 

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

  2526.         do k=1,km
    2527. 

  2527. !#ifdef DEBUG
    2528. 

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

  2529. !           ww(k),wa(k)
    2530. 

  2530. !#endif
    2531. 

  2531. ! mean wind is average of departure and new arrival winds
    2532. 

  2532.           ww(k) = 0.5* ( wd(k)+wa(k) )
    2533. 

  2533.         enddo
    2534. 

  2534.         was(:) = wa(:)
    2535. 

  2535.         n=n+1
    2536. 

  2536.         go to 100
    2537. 

  2537.       endif
    2538. 

  2538.       ist_loop : do ist = 1, 2
    2539. 

  2539.       if (ist.eq.2) then
    2540. 

  2540.        qa(:) = qa2(:)
    2541. 

  2541.       endif
    2542. 

  2542. !
    2543. 

  2543.       precip(i) = 0.
    2544. 

  2544. !
    2545. 

  2545. ! estimate values at arrival cell interface with monotone
    2546. 

  2546.       do k=2,km
    2547. 

  2547.         dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
    2548. 

  2548.         dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
    2549. 

  2549.         if( dip*dim.le.0.0 ) then
    2550. 

  2550.           qmi(k)=qa(k)
    2551. 

  2551.           qpi(k)=qa(k)
    2552. 

  2552.         else
    2553. 

  2553.           qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
    2554. 

  2554.           qmi(k)=2.0*qa(k)-qpi(k)
    2555. 

  2555.           if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
    2556. 

  2556.             qpi(k) = qa(k)
    2557. 

  2557.             qmi(k) = qa(k)
    2558. 

  2558.           endif
    2559. 

  2559.         endif
    2560. 

  2560.       enddo
    2561. 

  2561.       qpi(1)=qa(1)
    2562. 

  2562.       qmi(1)=qa(1)
    2563. 

  2563.       qmi(km+1)=qa(km+1)
    2564. 

  2564.       qpi(km+1)=qa(km+1)
    2565. 

  2565. !
    2566. 

  2566. ! interpolation to regular point
    2567. 

  2567.       qn = 0.0
    2568. 

  2568.       kb=1
    2569. 

  2569.       kt=1
    2570. 

  2570.       intp : do k=1,km
    2571. 

  2571.              kb=max(kb-1,1)
    2572. 

  2572.              kt=max(kt-1,1)
    2573. 

  2573. ! find kb and kt
    2574. 

  2574.              if( zi(k).ge.za(km+1) ) then
    2575. 

  2575.                exit intp
    2576. 

  2576.              else
    2577. 

  2577.                find_kb : do kk=kb,km
    2578. 

  2578.                          if( zi(k).le.za(kk+1) ) then
    2579. 

  2579.                            kb = kk
    2580. 

  2580.                            exit find_kb
    2581. 

  2581.                          else
    2582. 

  2582.                            cycle find_kb
    2583. 

  2583.                          endif
    2584. 

  2584.                enddo find_kb
    2585. 

  2585.                find_kt : do kk=kt,km
    2586. 

  2586.                          if( zi(k+1).le.za(kk) ) then
    2587. 

  2587.                            kt = kk
    2588. 

  2588.                            exit find_kt
    2589. 

  2589.                          else
    2590. 

  2590.                            cycle find_kt
    2591. 

  2591.                          endif
    2592. 

  2592.                enddo find_kt
    2593. 

  2593.                kt = kt - 1
    2594. 

  2594. ! compute q with piecewise constant method
    2595. 

  2595.                if( kt.eq.kb ) then
    2596. 

  2596.                  tl=(zi(k)-za(kb))/dza(kb)
    2597. 

  2597.                  th=(zi(k+1)-za(kb))/dza(kb)
    2598. 

  2598.                  tl2=tl*tl
    2599. 

  2599.                  th2=th*th
    2600. 

  2600.                  qqd=0.5*(qpi(kb)-qmi(kb))
    2601. 

  2601.                  qqh=qqd*th2+qmi(kb)*th
    2602. 

  2602.                  qql=qqd*tl2+qmi(kb)*tl
    2603. 

  2603.                  qn(k) = (qqh-qql)/(th-tl)
    2604. 

  2604.                else if( kt.gt.kb ) then
    2605. 

  2605.                  tl=(zi(k)-za(kb))/dza(kb)
    2606. 

  2606.                  tl2=tl*tl
    2607. 

  2607.                  qqd=0.5*(qpi(kb)-qmi(kb))
    2608. 

  2608.                  qql=qqd*tl2+qmi(kb)*tl
    2609. 

  2609.                  dql = qa(kb)-qql
    2610. 

  2610.                  zsum  = (1.-tl)*dza(kb)
    2611. 

  2611.                  qsum  = dql*dza(kb)
    2612. 

  2612.                  if( kt-kb.gt.1 ) then
    2613. 

  2613.                  do m=kb+1,kt-1
    2614. 

  2614.                    zsum = zsum + dza(m)
    2615. 

  2615.                    qsum = qsum + qa(m) * dza(m)
    2616. 

  2616.                  enddo
    2617. 

  2617.                  endif
    2618. 

  2618.                  th=(zi(k+1)-za(kt))/dza(kt)
    2619. 

  2619.                  th2=th*th
    2620. 

  2620.                  qqd=0.5*(qpi(kt)-qmi(kt))
    2621. 

  2621.                  dqh=qqd*th2+qmi(kt)*th
    2622. 

  2622.                  zsum  = zsum + th*dza(kt)
    2623. 

  2623.                  qsum  = qsum + dqh*dza(kt)
    2624. 

  2624.                  qn(k) = qsum/zsum
    2625. 

  2625.                endif
    2626. 

  2626.                cycle intp
    2627. 

  2627.              endif
    2628. 

  2628. !
    2629. 

  2629.        enddo intp
    2630. 

  2630. !
    2631. 

  2631. ! rain out
    2632. 

  2632.       sum_precip: do k=1,km
    2633. 

  2633.                     if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
    2634. 

  2634.                       precip(i) = precip(i) + qa(k)*dza(k)
    2635. 

  2635.                       cycle sum_precip
    2636. 

  2636.                     else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
    2637. 

  2637.                       precip(i) = precip(i) + qa(k)*(0.0-za(k))
    2638. 

  2638.                       exit sum_precip
    2639. 

  2639.                     endif
    2640. 

  2640.                     exit sum_precip
    2641. 

  2641.       enddo sum_precip
    2642. 

  2642. !
    2643. 

  2643. ! replace the new values
    2644. 

  2644.       if(ist.eq.1) then
    2645. 

  2645.         rql(i,:) = qn(:)
    2646. 

  2646.         precip1(i) = precip(i)
    2647. 

  2647.       else
    2648. 

  2648.         rql2(i,:) = qn(:)
    2649. 

  2649.         precip2(i) = precip(i)
    2650. 

  2650.       endif
    2651. 

  2651.       enddo ist_loop
    2652. 

  2652. !
    2653. 

  2653. ! ----------------------------------
    2654. 

  2654.       enddo i_loop
    2655. 

  2655. !
    2656. 

  2656.   END SUBROUTINE nislfv_rain_plm6
    2657. 

  2657. END MODULE module_mp_wdm6
    2658. 

  2658.