PageRenderTime 69ms CodeModel.GetById 18ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_mp_wdm5.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2015 lines | 1436 code | 0 blank | 579 comment | 101 complexity | a595cff214ff08031fdb6a3c0c1a78e6 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. 

  9. !Including inline expansion statistical function 
    10. 

  10. MODULE module_mp_wdm5
    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 :: avtr = 841.9       ! a constant for terminal velocity of rain  
    16. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  31.    REAL, PARAMETER, PRIVATE :: alpha = .12        ! .122 exponen factor for n0s
    32. 

  32.    REAL, PARAMETER, PRIVATE :: pfrz1 = 100.       ! constant in Biggs freezing  
    33. 

  33.    REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66       ! constant in Biggs freezing
    34. 

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

  35.    REAL, PARAMETER, PRIVATE :: ncmin = 1.e1       ! minimum value for Nc 
    36. 

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

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

  38. !
    39. 

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

  40.                                                   ! 1.008 for maritime air mass /1.0048 for conti
    41. 

  41.    REAL, PARAMETER, PRIVATE :: actk = 0.6         ! parameter for the CCN activation  
    42. 

  42.    REAL, PARAMETER, PRIVATE :: actr = 1.5         ! radius of activated CCN drops
    43. 

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

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

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

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

  47.    REAL, PARAMETER, PRIVATE :: di2000 = 20.e-4    ! parameter related with accretion and collection of cloud drops
    48. 

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

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

  50.    REAL, SAVE ::                                      &
    51. 

  51.              qc0, qck1,pidnc,bvtr1,bvtr2,bvtr3,bvtr4, &
    52. 

  52.              bvtr5,bvtr7,bvtr2o5,bvtr3o5,g1pbr,g2pbr, &
    53. 

  53.              g3pbr,g4pbr,g5pbr,g7pbr,g5pbro2,g7pbro2, &
    54. 

  54.              pvtr,pvtrn,eacrr,pacrr, pi,              &
    55. 

  55.              precr1,precr2,xmmax,roqimax,bvts1,       &
    56. 

  56.              bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs,     &
    57. 

  57.              g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
    58. 

  58.              pidn0s,pidnr,xlv1,pacrc,                 &
    59. 

  59.              rslopecmax,rslopec2max,rslopec3max,      &
    60. 

  60.              rslopermax,rslopesmax,rslopegmax,        &
    61. 

  61.              rsloperbmax,rslopesbmax,rslopegbmax,     &
    62. 

  62.              rsloper2max,rslopes2max,rslopeg2max,     &
    63. 

  63.              rsloper3max,rslopes3max,rslopeg3max
    64. 

  64. !
    65. 

  65. ! Specifies code-inlining of fpvs function in WDM52D below. JM 20040507
    66. 

  66. !
    67. 

  67. CONTAINS
    68. 

  68. !===================================================================
    69. 

  69. !
    70. 

  70.   SUBROUTINE wdm5(th, q, qc, qr, qi, qs                            &
    71. 

  71.                  ,nn, nc, nr                                       &
    72. 

  72.                  ,den, pii, p, delz                                &
    73. 

  73.                  ,delt,g, cpd, cpv, ccn0, rd, rv, t0c              &
    74. 

  74.                  ,ep1, ep2, qmin                                   &
    75. 

  75.                  ,XLS, XLV0, XLF0, den0, denr                      &
    76. 

  76.                  ,cliq,cice,psat                                   &
    77. 

  77.                  ,rain, rainncv                                    &
    78. 

  78.                  ,snow, snowncv                                    &
    79. 

  79.                  ,sr                                               &
    80. 

  80.                  ,itimestep                                        &
    81. 

  81.                  ,ids,ide, jds,jde, kds,kde                        &
    82. 

  82.                  ,ims,ime, jms,jme, kms,kme                        &
    83. 

  83.                  ,its,ite, jts,jte, kts,kte                        &
    84. 

  84.                                                                    )
    85. 

  85. !-------------------------------------------------------------------
    86. 

  86.   IMPLICIT NONE
    87. 

  87. !-------------------------------------------------------------------
    88. 

  88. !
    89. 

  89. !  This code is a WRF double-moment 5-class  mixed ice
    90. 

  90. !  microphyiscs scheme (WDM5). The WDM microphysics scheme predicts
    91. 

  91. !  number concentrations for warm rain species including clouds and
    92. 

  92. !  rain. cloud condensation nuclei (CCN) is also predicted.
    93. 

  93. !  The cold rain species including ice, snow, graupel follow the
    94. 

  94. !  WRF single-moment 5-class microphysics (WSM5)
    95. 

  95. !  in which theoretical background for WSM ice phase microphysics is
    96. 

  96. !  based on Hong et al. (2004). 
    97. 

  97. !  The WDM scheme is described in Lim and Hong (2010).
    98. 

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

  99. !
    100. 

  100. !  WDM5 cloud scheme
    101. 

  101. !
    102. 

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

  103. !
    104. 

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

  105. !
    106. 

  106. !  Reference) Lim and Hong (LH, 2010) Mon. Wea. Rev. 
    107. 

  107. !             Juang and Hong (JH, 2010) Mon. Wea. Rev.
    108. 

  108. !             Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
    109. 

  109. !             Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
    110. 

  110. !             Cohard and Pinty (CP, 2000) Quart. J. Roy. Meteor. Soc.
    111. 

  111. !             Khairoutdinov and Kogan (KK, 2000) Mon. Wea. Rev.
    112. 

  112. !             Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
    113. 

  113. !
    114. 

  114. !             Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
    115. 

  115. !             Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
    116. 

  116. !             Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
    117. 

  117. !
    118. 

  118.   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
    119. 

  119.                                       ims,ime, jms,jme, kms,kme , &
    120. 

  120.                                       its,ite, jts,jte, kts,kte
    121. 

  121.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    122. 

  122.         INTENT(INOUT) ::                                          &
    123. 

  123.                                                              th,  &
    124. 

  124.                                                               q,  &
    125. 

  125.                                                               qc, &
    126. 

  126.                                                               qi, &
    127. 

  127.                                                               qr, &
    128. 

  128.                                                               qs, &
    129. 

  129.                                                               nn, & 
    130. 

  130.                                                               nc, &
    131. 

  131.                                                               nr   
    132. 

  132.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    133. 

  133.         INTENT(IN   ) ::                                          &
    134. 

  134.                                                              den, &
    135. 

  135.                                                              pii, &
    136. 

  136.                                                                p, &
    137. 

  137.                                                             delz
    138. 

  138.   REAL, INTENT(IN   ) ::                                    delt, &
    139. 

  139.                                                                g, &
    140. 

  140.                                                               rd, &
    141. 

  141.                                                               rv, &
    142. 

  142.                                                              t0c, &
    143. 

  143.                                                             den0, &
    144. 

  144.                                                              cpd, &
    145. 

  145.                                                              cpv, &
    146. 

  146.                                                             ccn0, &
    147. 

  147.                                                              ep1, &
    148. 

  148.                                                              ep2, &
    149. 

  149.                                                             qmin, &
    150. 

  150.                                                              XLS, &
    151. 

  151.                                                             XLV0, &
    152. 

  152.                                                             XLF0, &
    153. 

  153.                                                             cliq, &
    154. 

  154.                                                             cice, &
    155. 

  155.                                                             psat, &
    156. 

  156.                                                             denr
    157. 

  157.   INTEGER, INTENT(IN   ) ::                            itimestep
    158. 

  158.   REAL, DIMENSION( ims:ime , jms:jme ),                           &
    159. 

  159.         INTENT(INOUT) ::                                    rain, &
    160. 

  160.                                                          rainncv, &
    161. 

  161.                                                               sr
    162. 

  162.   REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                 &
    163. 

  163.         INTENT(INOUT) ::                                    snow, &
    164. 

  164.                                                          snowncv
    165. 

  165. ! LOCAL VAR
    166. 

  166.   REAL, DIMENSION( its:ite , kts:kte ) ::   t
    167. 

  167.   REAL, DIMENSION( its:ite , kts:kte, 2 ) ::   qci, qrs
    168. 

  168.   REAL, DIMENSION( its:ite , kts:kte, 3 ) ::   ncr
    169. 

  169.   CHARACTER*256 :: emess
    170. 

  170.   INTEGER :: mkx_test
    171. 

  171.   INTEGER ::               i,j,k
    172. 

  172. !-------------------------------------------------------------------
    173. 

  173. #ifndef RUN_ON_GPU
    174. 

  174.       IF (itimestep .eq. 1) THEN
    175. 

  175.         DO j=jms,jme
    176. 

  176.            DO k=kms,kme
    177. 

  177.            DO i=ims,ime
    178. 

  178.               nn(i,k,j) = ccn0
    179. 

  179.            ENDDO
    180. 

  180.            ENDDO
    181. 

  181.         ENDDO
    182. 

  182.       ENDIF
    183. 

  183. !
    184. 

  184.       DO j=jts,jte
    185. 

  185.          DO k=kts,kte
    186. 

  186.          DO i=its,ite
    187. 

  187.             t(i,k)=th(i,k,j)*pii(i,k,j)
    188. 

  188.             qci(i,k,1) = qc(i,k,j)
    189. 

  189.             qci(i,k,2) = qi(i,k,j)
    190. 

  190.             qrs(i,k,1) = qr(i,k,j)
    191. 

  191.             qrs(i,k,2) = qs(i,k,j)
    192. 

  192.             ncr(i,k,1) = nn(i,k,j)                          
    193. 

  193.             ncr(i,k,2) = nc(i,k,j)                         
    194. 

  194.             ncr(i,k,3) = nr(i,k,j)                          
    195. 

  195.          ENDDO
    196. 

  196.          ENDDO
    197. 

  197.          !  Sending array starting locations of optional variables may cause
    198. 

  198.          !  troubles, so we explicitly change the call.
    199. 

  199.          CALL wdm52D(t, q(ims,kms,j), qci, qrs, ncr               &
    200. 

  200.                     ,den(ims,kms,j)                               &
    201. 

  201.                     ,p(ims,kms,j), delz(ims,kms,j)                &
    202. 

  202.                     ,delt,g, cpd, cpv, ccn0, rd, rv, t0c          &
    203. 

  203.                     ,ep1, ep2, qmin                               &
    204. 

  204.                     ,XLS, XLV0, XLF0, den0, denr                  &
    205. 

  205.                     ,cliq,cice,psat                               &
    206. 

  206.                     ,j                                            &
    207. 

  207.                     ,rain(ims,j),rainncv(ims,j)                   &
    208. 

  208.                     ,sr(ims,j)                                    &
    209. 

  209.                     ,ids,ide, jds,jde, kds,kde                    &
    210. 

  210.                     ,ims,ime, jms,jme, kms,kme                    &
    211. 

  211.                     ,its,ite, jts,jte, kts,kte                    &
    212. 

  212.                     ,snow(ims,j),snowncv(ims,j)                   &
    213. 

  213.                                                                   )
    214. 

  214.          DO K=kts,kte
    215. 

  215.          DO I=its,ite
    216. 

  216.             th(i,k,j)=t(i,k)/pii(i,k,j)
    217. 

  217.             qc(i,k,j) = qci(i,k,1)
    218. 

  218.             qi(i,k,j) = qci(i,k,2)
    219. 

  219.             qr(i,k,j) = qrs(i,k,1)
    220. 

  220.             qs(i,k,j) = qrs(i,k,2)
    221. 

  221.             nn(i,k,j) = ncr(i,k,1)
    222. 

  222.             nc(i,k,j) = ncr(i,k,2)
    223. 

  223.             nr(i,k,j) = ncr(i,k,3)
    224. 

  224.          ENDDO
    225. 

  225.          ENDDO
    226. 

  226.       ENDDO
    227. 

  227. #else
    228. 

  228.       CALL get_wsm5_gpu_levels ( mkx_test )
    229. 

  229.       IF ( mkx_test .LT. kte ) THEN
    230. 

  230.         WRITE(emess,*)'Number of levels compiled for GPU WSM5 too small. ',    &
    231. 

  231.                       mkx_test,' < ',kte
    232. 

  232.         CALL wrf_error_fatal(emess)
    233. 

  233.       ENDIF
    234. 

  234.       CALL wsm5_host (                                                         &
    235. 

  235.                     th(its:ite,kts:kte,jts:jte), pii(its:ite,kts:kte,jts:jte)  &
    236. 

  236.                    ,q(its:ite,kts:kte,jts:jte), qc(its:ite,kts:kte,jts:jte)    &
    237. 

  237.                    ,qi(its:ite,kts:kte,jts:jte), qr(its:ite,kts:kte,jts:jte)   &
    238. 

  238.                    ,qs(its:ite,kts:kte,jts:jte), den(its:ite,kts:kte,jts:jte)  &
    239. 

  239.                    ,p(its:ite,kts:kte,jts:jte), delz(its:ite,kts:kte,jts:jte)  &
    240. 

  240.                    ,delt                                                       &
    241. 

  241.                    ,rain(its:ite,jts:jte),rainncv(its:ite,jts:jte)             &
    242. 

  242.                    ,snow(its:ite,jts:jte),snowncv(its:ite,jts:jte)             &
    243. 

  243.                    ,sr(its:ite,jts:jte)                                        &
    244. 

  244.                    ,its, ite,  jts, jte,  kts, kte                             &
    245. 

  245.                    ,its, ite,  jts, jte,  kts, kte                             &
    246. 

  246.                    ,its, ite,  jts, jte,  kts, kte                             &
    247. 

  247.           )
    248. 

  248. #endif
    249. 

  249.   END SUBROUTINE wdm5
    250. 

  250. !===================================================================
    251. 

  251. !
    252. 

  252.   SUBROUTINE wdm52D(t, q, qci, qrs, ncr, den, p, delz             &
    253. 

  253.                    ,delt,g, cpd, cpv, ccn0, rd, rv, t0c           &
    254. 

  254.                    ,ep1, ep2, qmin                                &
    255. 

  255.                    ,XLS, XLV0, XLF0, den0, denr                   &
    256. 

  256.                    ,cliq,cice,psat                                &
    257. 

  257.                    ,lat                                           &
    258. 

  258.                    ,rain,rainncv                                  &
    259. 

  259.                    ,sr                                            &
    260. 

  260.                    ,ids,ide, jds,jde, kds,kde                     &
    261. 

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

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

  263.                    ,snow,snowncv                                  &
    264. 

  264.                                                                   )
    265. 

  265. !-------------------------------------------------------------------
    266. 

  266.   IMPLICIT NONE
    267. 

  267. !-------------------------------------------------------------------
    268. 

  268.   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
    269. 

  269.                                       ims,ime, jms,jme, kms,kme , &
    270. 

  270.                                       its,ite, jts,jte, kts,kte,  &
    271. 

  271.                                       lat
    272. 

  272.   REAL, DIMENSION( its:ite , kts:kte ),                           &
    273. 

  273.         INTENT(INOUT) ::                                          &
    274. 

  274.                                                                t
    275. 

  275.   REAL, DIMENSION( its:ite , kts:kte, 2 ),                        &
    276. 

  276.         INTENT(INOUT) ::                                          &
    277. 

  277.                                                              qci, &
    278. 

  278.                                                              qrs 
    279. 

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

  280.         INTENT(INOUT) ::                                          &
    281. 

  281.                                                              ncr
    282. 

  282.   REAL, DIMENSION( ims:ime , kms:kme ),                           &
    283. 

  283.         INTENT(INOUT) ::                                          &
    284. 

  284.                                                                q
    285. 

  285.   REAL, DIMENSION( ims:ime , kms:kme ),                           &
    286. 

  286.         INTENT(IN   ) ::                                          &
    287. 

  287.                                                              den, &
    288. 

  288.                                                                p, &
    289. 

  289.                                                             delz
    290. 

  290.   REAL, INTENT(IN   ) ::                                    delt, &
    291. 

  291.                                                                g, &
    292. 

  292.                                                              cpd, &
    293. 

  293.                                                              cpv, &
    294. 

  294.                                                             ccn0, &
    295. 

  295.                                                              t0c, &
    296. 

  296.                                                             den0, &
    297. 

  297.                                                               rd, &
    298. 

  298.                                                               rv, &
    299. 

  299.                                                              ep1, &
    300. 

  300.                                                              ep2, &
    301. 

  301.                                                             qmin, &
    302. 

  302.                                                              XLS, &
    303. 

  303.                                                             XLV0, &
    304. 

  304.                                                             XLF0, &
    305. 

  305.                                                             cliq, &
    306. 

  306.                                                             cice, &
    307. 

  307.                                                             psat, &
    308. 

  308.                                                             denr
    309. 

  309.   REAL, DIMENSION( ims:ime ),                                     &
    310. 

  310.         INTENT(INOUT) ::                                    rain, &
    311. 

  311.                                                          rainncv, &
    312. 

  312.                                                               sr
    313. 

  313.   REAL, DIMENSION( ims:ime ),     OPTIONAL,                       &
    314. 

  314.         INTENT(INOUT) ::                                    snow, &
    315. 

  315.                                                          snowncv
    316. 

  316. ! LOCAL VAR
    317. 

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

  318.         rh, qs, rslope, rslope2, rslope3, rslopeb,                &
    319. 

  319.         falk, fall, work1, qrs_tmp
    320. 

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

  321.         rslopec, rslopec2,rslopec3                           
    322. 

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

  323.         avedia 
    324. 

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

  325.         workn,falln,falkn
    326. 

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

  327.         works
    328. 

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

  329.         den_tmp, delz_tmp, ncr_tmp
    330. 

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

  331.               falkc, work1c, work2c, fallc
    332. 

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

  333.         pcact, praut, psaut, prevp, psdep, pracw, psaci, psacw,   &  
    334. 

  334.         pigen, pidep, pcond,                                      &
    335. 

  335.         xl, cpm, work2, psmlt, psevp, denfac, xni,                &
    336. 

  336.         n0sfac, denqrs2, denqci
    337. 

  337.   REAL, DIMENSION( its:ite ) ::                                   &
    338. 

  338.         delqrs2, delqi
    339. 

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

  340.         nraut, nracw, ncevp, nccol, nrcol,                        &
    341. 

  341.         nsacw, nseml, ncact 
    342. 

  342.   REAL :: ifac, sfac
    343. 

  343.   REAL, DIMENSION(its:ite) :: tstepsnow
    344. 

  344. !
    345. 

  345. #define WSM_NO_CONDITIONAL_IN_VECTOR
    346. 

  346. #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
    347. 

  347.   REAL, DIMENSION(its:ite) :: xal, xbl
    348. 

  348. #endif
    349. 

  349. ! variables for optimization
    350. 

  350.   REAL, DIMENSION( its:ite )           :: tvec1
    351. 

  351.   INTEGER, DIMENSION( its:ite ) :: mnstep, numndt
    352. 

  352.   INTEGER, DIMENSION( its:ite ) :: mstep, numdt
    353. 

  353.   REAL, DIMENSION(its:ite) :: rmstep
    354. 

  354.   REAL dtcldden, rdelz, rdtcld
    355. 

  355.   LOGICAL, DIMENSION( its:ite ) :: flgcld
    356. 

  356.   REAL  ::                                                        &
    357. 

  357.             cpmcal, xlcal, lamdac, diffus,                        &
    358. 

  358.             viscos, xka, venfac, conden, diffac,                  &
    359. 

  359.             x, y, z, a, b, c, d, e,                               &
    360. 

  360.             ndt, qdt, holdrr, holdrs, supcol, supcolt, pvt,       &
    361. 

  361.             coeres, supsat, dtcld, xmi, eacrs, satdt,             &
    362. 

  362.             vt2i,vt2s,acrfac, coecol,                             &
    363. 

  363.             nfrzdtr, nfrzdtc,                                     &
    364. 

  364.             taucon, lencon, lenconcr,                             &
    365. 

  365.             qimax, diameter, xni0, roqi0,                         &
    366. 

  366.             fallsum, fallsum_qsi, xlwork2, factor, source,        &
    367. 

  367.             value, xlf, pfrzdtc, pfrzdtr, supice
    368. 

  368.   REAL :: temp 
    369. 

  369.   REAL  :: holdc, holdci
    370. 

  370.   INTEGER :: i, j, k, mstepmax,                                                &
    371. 

  371.             iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
    372. 

  372. ! Temporaries used for inlining fpvs function
    373. 

  373.   REAL  :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
    374. 

  374.   REAL  :: logtr
    375. 

  375. !
    376. 

  376. !=================================================================
    377. 

  377. !   compute internal functions
    378. 

  378. !
    379. 

  379.       cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
    380. 

  380.       xlcal(x) = xlv0-xlv1*(x-t0c)
    381. 

  381. !----------------------------------------------------------------
    382. 

  382. !     size distributions: (x=mixing ratio, y=air density):
    383. 

  383. !     valid for mixing ratio > 1.e-9 kg/kg.
    384. 

  384. !
    385. 

  385. ! Optimizatin : A**B => exp(log(A)*(B))
    386. 

  386.       lamdac(x,y,z)= exp(log(((pidnc*z)/(x*y)))*((.33333333)))
    387. 

  387. !
    388. 

  388. !----------------------------------------------------------------
    389. 

  389. !     diffus: diffusion coefficient of the water vapor
    390. 

  390. !     viscos: kinematic viscosity(m2s-1)
    391. 

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

  392. !     viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y  
    393. 

  393. !     xka(x,y) = 1.414e3*viscos(x,y)*y
    394. 

  394. !     diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
    395. 

  395. !     venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333)))         &
    396. 

  396. !                    /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
    397. 

  397. !     conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
    398. 

  398. !
    399. 

  399. !
    400. 

  400.       idim = ite-its+1
    401. 

  401.       kdim = kte-kts+1
    402. 

  402. !
    403. 

  403. !----------------------------------------------------------------
    404. 

  404. !     paddint 0 for negative values generated by dynamics
    405. 

  405. !
    406. 

  406.       do k = kts, kte
    407. 

  407.         do i = its, ite
    408. 

  408.           qci(i,k,1) = max(qci(i,k,1),0.0)
    409. 

  409.           qrs(i,k,1) = max(qrs(i,k,1),0.0)
    410. 

  410.           qci(i,k,2) = max(qci(i,k,2),0.0)
    411. 

  411.           qrs(i,k,2) = max(qrs(i,k,2),0.0)
    412. 

  412.           ncr(i,k,1) = max(ncr(i,k,1),0.)
    413. 

  413.           ncr(i,k,2) = max(ncr(i,k,2),0.)
    414. 

  414.           ncr(i,k,3) = max(ncr(i,k,3),0.) 
    415. 

  415.         enddo
    416. 

  416.       enddo
    417. 

  417. !
    418. 

  418. !     latent heat for phase changes and heat capacity. neglect the
    419. 

  419. !     changes during microphysical process calculation
    420. 

  420. !     emanuel(1994)
    421. 

  421. !
    422. 

  422.       do k = kts, kte
    423. 

  423.         do i = its, ite
    424. 

  424.           cpm(i,k) = cpmcal(q(i,k))
    425. 

  425.           xl(i,k) = xlcal(t(i,k))
    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
    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.         sr(i) = 0.
    438. 

  438. ! new local array to catch step snow
    439. 

  439.         tstepsnow(i) = 0.
    440. 

  440.       enddo
    441. 

  441. !
    442. 

  442. !----------------------------------------------------------------
    443. 

  443. !     compute the minor time steps.
    444. 

  444. !
    445. 

  445.       loops = max(nint(delt/dtcldcr),1)
    446. 

  446.       dtcld = delt/loops
    447. 

  447.       if(delt.le.dtcldcr) dtcld = delt
    448. 

  448. !
    449. 

  449.       do loop = 1,loops
    450. 

  450. !
    451. 

  451. !----------------------------------------------------------------
    452. 

  452. !     initialize the large scale variables
    453. 

  453. !
    454. 

  454.       do i = its, ite
    455. 

  455.         mstep(i) = 1
    456. 

  456.         mnstep(i) = 1
    457. 

  457.         flgcld(i) = .true.
    458. 

  458.       enddo
    459. 

  459. !
    460. 

  460. !     do k = kts, kte
    461. 

  461. !       do i = its, ite
    462. 

  462. !         denfac(i,k) = sqrt(den0/den(i,k))
    463. 

  463. !       enddo
    464. 

  464. !     enddo
    465. 

  465.       do k = kts, kte
    466. 

  466.         CALL VREC( tvec1(its), den(its,k), ite-its+1)
    467. 

  467.         do i = its, ite
    468. 

  468.           tvec1(i) = tvec1(i)*den0
    469. 

  469.         enddo
    470. 

  470.         CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
    471. 

  471.       enddo
    472. 

  472. !
    473. 

  473. ! Inline expansion for fpvs
    474. 

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

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

  476.       hsub = xls
    477. 

  477.       hvap = xlv0
    478. 

  478.       cvap = cpv
    479. 

  479.       ttp=t0c+0.01
    480. 

  480.       dldt=cvap-cliq
    481. 

  481.       xa=-dldt/rv
    482. 

  482.       xb=xa+hvap/(rv*ttp)
    483. 

  483.       dldti=cvap-cice
    484. 

  484.       xai=-dldti/rv
    485. 

  485.       xbi=xai+hsub/(rv*ttp)
    486. 

  486. ! this is for compilers where the conditional inhibits vectorization
    487. 

  487. #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
    488. 

  488.       do k = kts, kte
    489. 

  489.         do i = its, ite
    490. 

  490.           if(t(i,k).lt.ttp) then
    491. 

  491.             xal(i) = xai
    492. 

  492.             xbl(i) = xbi
    493. 

  493.           else
    494. 

  494.             xal(i) = xa
    495. 

  495.             xbl(i) = xb
    496. 

  496.           endif
    497. 

  497.         enddo
    498. 

  498.         do i = its, ite
    499. 

  499.           tr=ttp/t(i,k)
    500. 

  500.           logtr=log(tr)
    501. 

  501.           qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
    502. 

  502.           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
    503. 

  503.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    504. 

  504.           qs(i,k,1) = max(qs(i,k,1),qmin)
    505. 

  505.           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
    506. 

  506.           qs(i,k,2)=psat*exp(logtr*(xal(i))+xbl(i)*(1.-tr))
    507. 

  507.           qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
    508. 

  508.           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
    509. 

  509.           qs(i,k,2) = max(qs(i,k,2),qmin)
    510. 

  510.           rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
    511. 

  511.         enddo
    512. 

  512.       enddo
    513. 

  513. #else
    514. 

  514.       do k = kts, kte
    515. 

  515.         do i = its, ite
    516. 

  516.           tr=ttp/t(i,k)
    517. 

  517.           logtr=log(tr)
    518. 

  518.           qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
    519. 

  519.           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
    520. 

  520.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    521. 

  521.           qs(i,k,1) = max(qs(i,k,1),qmin)
    522. 

  522.           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
    523. 

  523.           if(t(i,k).lt.ttp) then
    524. 

  524.             qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
    525. 

  525.           else
    526. 

  526.             qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
    527. 

  527.           endif
    528. 

  528.           qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
    529. 

  529.           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
    530. 

  530.           qs(i,k,2) = max(qs(i,k,2),qmin)
    531. 

  531.           rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
    532. 

  532.         enddo
    533. 

  533.       enddo
    534. 

  534. #endif
    535. 

  535. !
    536. 

  536. !----------------------------------------------------------------
    537. 

  537. !     initialize the variables for microphysical physics
    538. 

  538. !
    539. 

  539. !
    540. 

  540.       do k = kts, kte
    541. 

  541.         do i = its, ite
    542. 

  542.           prevp(i,k) = 0.
    543. 

  543.           psdep(i,k) = 0.
    544. 

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

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

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

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

  548.           psacw(i,k) = 0.
    549. 

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

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

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

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

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

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

  555.           falk(i,k,1) = 0.
    556. 

  556.           falk(i,k,2) = 0.
    557. 

  557.           fall(i,k,1) = 0.
    558. 

  558.           fall(i,k,2) = 0.
    559. 

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

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

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

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

  563.           xni(i,k) = 1.e3
    564. 

  564.           nsacw(i,k) = 0.
    565. 

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

  566.           nracw(i,k) = 0.
    567. 

  567.           nccol(i,k) = 0.
    568. 

  568.           nrcol(i,k) = 0.
    569. 

  569.           ncact(i,k) = 0.
    570. 

  570.           nraut(i,k) = 0.
    571. 

  571.           ncevp(i,k) = 0.
    572. 

  572.         enddo
    573. 

  573.       enddo
    574. 

  574. !
    575. 

  575. !----------------------------------------------------------------
    576. 

  576. !     compute the fallout term:
    577. 

  577. !     first, vertical terminal velosity for minor loops
    578. 

  578. !
    579. 

  579.       do k = kts, kte
    580. 

  580.         do i = its, ite
    581. 

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

  582.             rslopec(i,k) = rslopecmax
    583. 

  583.             rslopec2(i,k) = rslopec2max
    584. 

  584.             rslopec3(i,k) = rslopec3max
    585. 

  585.           else
    586. 

  586.             rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
    587. 

  587.             rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
    588. 

  588.             rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
    589. 

  589.           endif
    590. 

  590. !-------------------------------------------------------------
    591. 

  591. ! Ni: ice crystal number concentraiton   [HDC 5c]
    592. 

  592. !-------------------------------------------------------------
    593. 

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

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

  595.           temp = (den(i,k)*max(qci(i,k,2),qmin))
    596. 

  596.           temp = sqrt(sqrt(temp*temp*temp))
    597. 

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

  598.         enddo
    599. 

  599.       enddo
    600. 

  600.       do k = kts, kte
    601. 

  601.         do i = its, ite
    602. 

  602.           qrs_tmp(i,k,1) = qrs(i,k,1)
    603. 

  603.           qrs_tmp(i,k,2) = qrs(i,k,2)
    604. 

  604.           ncr_tmp(i,k) = ncr(i,k,3)
    605. 

  605.         enddo
    606. 

  606.       enddo
    607. 

  607.       call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
    608. 

  608.                      rslope3,work1,workn,its,ite,kts,kte)
    609. 

  609. !----------------------------------------------------------------
    610. 

  610. !     compute the fallout term:
    611. 

  611. !     first, vertical terminal velosity for minor loops
    612. 

  612. !----------------------------------------------------------------
    613. 

  613. !
    614. 

  614. ! vt update for qr and nr
    615. 

  615.       mstepmax = 1
    616. 

  616.       numdt = 1
    617. 

  617.       do k = kte, kts, -1
    618. 

  618.         do i = its, ite
    619. 

  619.           work1(i,k,1) = work1(i,k,1)/delz(i,k)
    620. 

  620.           workn(i,k) = workn(i,k)/delz(i,k)
    621. 

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

  622.           if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
    623. 

  623.         enddo
    624. 

  624.       enddo
    625. 

  625.       do i = its, ite
    626. 

  626.         if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
    627. 

  627.       enddo
    628. 

  628. !
    629. 

  629.       do n = 1, mstepmax
    630. 

  630.         k = kte
    631. 

  631.         do i = its, ite
    632. 

  632.           if(n.le.mstep(i)) then
    633. 

  633.             falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
    634. 

  634.             falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
    635. 

  635.             fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
    636. 

  636.             falln(i,k) = falln(i,k)+falkn(i,k)
    637. 

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

  638.             ncr(i,k,3) = max(ncr(i,k,3)-falkn(i,k)*dtcld,0.)
    639. 

  639.           endif
    640. 

  640.         enddo
    641. 

  641.         do k = kte-1, kts, -1
    642. 

  642.           do i = its, ite
    643. 

  643.             if(n.le.mstep(i)) then
    644. 

  644.               falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
    645. 

  645.               falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
    646. 

  646.               fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
    647. 

  647.               falln(i,k) = falln(i,k)+falkn(i,k)
    648. 

  648.               qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1)           &
    649. 

  649.                           *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
    650. 

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

  651.                           /delz(i,k))*dtcld,0.)
    652. 

  652.             endif
    653. 

  653.           enddo
    654. 

  654.         enddo
    655. 

  655.         do k = kts, kte
    656. 

  656.           do i = its, ite
    657. 

  657.             qrs_tmp(i,k,1) = qrs(i,k,1)
    658. 

  658.             ncr_tmp(i,k) = ncr(i,k,3)
    659. 

  659.           enddo
    660. 

  660.         enddo
    661. 

  661.         call slope_rain(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
    662. 

  662.                      rslope3,work1,workn,its,ite,kts,kte)
    663. 

  663.         do k = kte, kts, -1
    664. 

  664.           do i = its, ite
    665. 

  665.             work1(i,k,1) = work1(i,k,1)/delz(i,k)
    666. 

  666.             workn(i,k) = workn(i,k)/delz(i,k)
    667. 

  667.           enddo
    668. 

  668.         enddo
    669. 

  669.       enddo
    670. 

  670. ! for semi
    671. 

  671.       do k = kte, kts, -1
    672. 

  672.         do i = its, ite
    673. 

  673.           works(i,k) = work1(i,k,2)
    674. 

  674.           denqrs2(i,k) = den(i,k)*qrs(i,k,2)
    675. 

  675.           if(qrs(i,k,2).le.0.0) works(i,k) = 0.0
    676. 

  676.         enddo
    677. 

  677.       enddo
    678. 

  678.       call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,works,denqrs2,  &
    679. 

  679.                            delqrs2,dtcld,2,1)
    680. 

  680.       do k = kts, kte
    681. 

  681.         do i = its, ite
    682. 

  682.           qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
    683. 

  683.           fall(i,k,2) = denqrs2(i,k)*works(i,k)/delz(i,k)
    684. 

  684.         enddo
    685. 

  685.       enddo
    686. 

  686.       do i = its, ite
    687. 

  687.         fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
    688. 

  688.       enddo
    689. 

  689.       do k = kts, kte
    690. 

  690.         do i = its, ite
    691. 

  691.           qrs_tmp(i,k,1) = qrs(i,k,1)
    692. 

  692.           qrs_tmp(i,k,2) = qrs(i,k,2)
    693. 

  693.           ncr_tmp(i,k) = ncr(i,k,3)
    694. 

  694.         enddo
    695. 

  695.       enddo
    696. 

  696.       call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
    697. 

  697.                      rslope3,work1,workn,its,ite,kts,kte)
    698. 

  698. !
    699. 

  699.       do k = kte, kts, -1
    700. 

  700.         do i = its, ite
    701. 

  701.           supcol = t0c-t(i,k)
    702. 

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

  703.           if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
    704. 

  704. !----------------------------------------------------------------
    705. 

  705. ! psmlt: melting of snow [HL A33] [RH83 A25]
    706. 

  706. !       (T>T0: QS->QR)
    707. 

  707. !----------------------------------------------------------------
    708. 

  708.             xlf = xlf0
    709. 

  709. !           work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
    710. 

  710.             work2(i,k)= (exp(log(((1.496e-6*((t(i,k))*sqrt(t(i,k)))        &
    711. 

  711.                         /((t(i,k))+120.)/(den(i,k)))/(8.794e-5             &
    712. 

  712.                         *exp(log(t(i,k))*(1.81))/p(i,k))))                 &
    713. 

  713.                         *((.3333333)))/sqrt((1.496e-6*((t(i,k))            &
    714. 

  714.                         *sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k))))        &
    715. 

  715.                         *sqrt(sqrt(den0/(den(i,k)))))
    716. 

  716.             coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    717. 

  717. !           psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2.       &
    718. 

  718. !                       *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
    719. 

  719. !                       *work2(i,k)*coeres)
    720. 

  720.             psmlt(i,k) = (1.414e3*(1.496e-6 * ((t(i,k))*sqrt(t(i,k)))      &
    721. 

  721.                          /((t(i,k))+120.)/(den(i,k)))*(den(i,k)))/xlf      &
    722. 

  722.                         *(t0c-t(i,k))*pi/2.*n0sfac(i,k)                    &
    723. 

  723.                         *(precs1*rslope2(i,k,2)+precs2*work2(i,k)*coeres)
    724. 

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

  725.                           /mstep(i)),0.)
    726. 

  726. !-------------------------------------------------------------------
    727. 

  727. ! nsmlt: melgin of snow  [LH A27]
    728. 

  728. !       (T>T0: ->NR)
    729. 

  729. !-------------------------------------------------------------------
    730. 

  730.             if(qrs(i,k,2).gt.qcrmin) then
    731. 

  731.               sfac = rslope(i,k,2)*n0s*n0sfac(i,k)*mstep(i)/qrs(i,k,2)
    732. 

  732.               ncr(i,k,3) = ncr(i,k,3) - sfac*psmlt(i,k)
    733. 

  733.             endif
    734. 

  734.             qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
    735. 

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

  736.             t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
    737. 

  737.           endif
    738. 

  738.         enddo
    739. 

  739.       enddo
    740. 

  740. !---------------------------------------------------------------
    741. 

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

  742. !---------------------------------------------------------------
    743. 

  743.       do k = kte, kts, -1
    744. 

  744.         do i = its, ite
    745. 

  745.           if(qci(i,k,2).le.0.) then
    746. 

  746.             work1c(i,k) = 0.
    747. 

  747.           else
    748. 

  748.             xmi = den(i,k)*qci(i,k,2)/xni(i,k)
    749. 

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

  750.             work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
    751. 

  751.           endif
    752. 

  752.         enddo
    753. 

  753.       enddo
    754. 

  754. !
    755. 

  755. !  forward semi-laglangian scheme (JH), PCM (piecewise constant),  (linear)
    756. 

  756. !
    757. 

  757.       do k = kte, kts, -1
    758. 

  758.         do i = its, ite
    759. 

  759.           denqci(i,k) = den(i,k)*qci(i,k,2)
    760. 

  760.         enddo
    761. 

  761.       enddo
    762. 

  762.       call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci,  &
    763. 

  763.                            delqi,dtcld,1,0)
    764. 

  764.       do k = kts, kte
    765. 

  765.         do i = its, ite
    766. 

  766.           qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
    767. 

  767.         enddo
    768. 

  768.       enddo
    769. 

  769.       do i = its, ite
    770. 

  770.         fallc(i,1) = delqi(i)/delz(i,1)/dtcld
    771. 

  771.       enddo
    772. 

  772. !
    773. 

  773. !----------------------------------------------------------------
    774. 

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

  775. !
    776. 

  776.       do i = its, ite
    777. 

  777.         fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
    778. 

  778.         fallsum_qsi = fall(i,1,2)+fallc(i,1)
    779. 

  779.         if(fallsum.gt.0.) then
    780. 

  780.           rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rainncv(i)
    781. 

  781.           rain(i) = fallsum*delz(i,1)/denr*dtcld*1000.+rain(i)
    782. 

  782.         endif
    783. 

  783.        if(fallsum_qsi.gt.0.) then
    784. 

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

  785.           if (PRESENT (snowncv) .and. PRESENT (snow)) then 
    786. 

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

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

  788.           endif
    789. 

  789.         endif 
    790. 

  790.         if(fallsum.gt.0.)sr(i)= tstepsnow(i)/(rainncv(i)+1.e-12)
    791. 

  791.       enddo
    792. 

  792. !
    793. 

  793. !---------------------------------------------------------------
    794. 

  794. ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
    795. 

  795. !       (T>T0: QI->QC)
    796. 

  796. !---------------------------------------------------------------
    797. 

  797.       do k = kts, kte
    798. 

  798.         do i = its, ite
    799. 

  799.           supcol = t0c-t(i,k)
    800. 

  800.           xlf = xls-xl(i,k)
    801. 

  801.           if(supcol.lt.0.) xlf = xlf0
    802. 

  802.           if(supcol.lt.0 .and. qci(i,k,2).gt.0.) then
    803. 

  803.             qci(i,k,1) = qci(i,k,1)+qci(i,k,2)
    804. 

  804. !---------------------------------------------------------------
    805. 

  805. ! nimlt: instantaneous melting of cloud ice  [LH A18]
    806. 

  806. !        (T>T0: ->NC)
    807. 

  807. !--------------------------------------------------------------          
    808. 

  808.             ncr(i,k,2) = ncr(i,k,2) + xni(i,k)
    809. 

  809.             t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
    810. 

  810.             qci(i,k,2) = 0.
    811. 

  811.           endif
    812. 

  812. !---------------------------------------------------------------
    813. 

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

  814. !        (T<-40C: QC->QI)
    815. 

  815. !---------------------------------------------------------------
    816. 

  816.           if(supcol.gt.40. .and. qci(i,k,1).gt.0.) then
    817. 

  817.             qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
    818. 

  818. !---------------------------------------------------------------
    819. 

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

  820. !        (T<-40C: NC->)
    821. 

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

  822.             if(ncr(i,k,2).gt.0.) ncr(i,k,2) = 0.
    823. 

  823.             t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
    824. 

  824.             qci(i,k,1) = 0.
    825. 

  825.           endif
    826. 

  826. !---------------------------------------------------------------
    827. 

  827. ! pihtf: heterogeneous freezing of cloud water [HL A44]
    828. 

  828. !        (T0>T>-40C: QC->QI)
    829. 

  829. !---------------------------------------------------------------
    830. 

  830.           if(supcol.gt.0. .and. qci(i,k,1).gt.0.) then
    831. 

  831.             supcolt=min(supcol,70.)
    832. 

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

  833.                    *ncr(i,k,2)*rslopec3(i,k)*rslopec3(i,k)/18.*dtcld,qci(i,k,1))       
    834. 

  834. !---------------------------------------------------------------
    835. 

  835. ! nihtf: heterogeneous  of cloud water  [LH A16]
    836. 

  836. !         (T0>T>-40C: NC->)
    837. 

  837. !---------------------------------------------------------------
    838. 

  838.             if(ncr(i,k,2).gt.ncmin) then
    839. 

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

  840.                       *rslopec3(i,k)/6.*dtcld,ncr(i,k,2))
    841. 

  841.               ncr(i,k,2) = ncr(i,k,2) - nfrzdtc
    842. 

  842.             endif
    843. 

  843.             qci(i,k,2) = qci(i,k,2) + pfrzdtc
    844. 

  844.             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
    845. 

  845.             qci(i,k,1) = qci(i,k,1)-pfrzdtc
    846. 

  846.            endif
    847. 

  847. !---------------------------------------------------------------
    848. 

  848. ! psfrz: freezing of rain water [HL A20] [LFO 45]
    849. 

  849. !        (T<T0, QR->QS)
    850. 

  850. !---------------------------------------------------------------
    851. 

  851.           if(supcol.gt.0. .and. qrs(i,k,1).gt.0.) then
    852. 

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

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

  854.                   *(exp(pfrz2*supcolt)-1.)*rslope3(i,k,1)*rslope3(i,k,1)       &
    855. 

  855.                   *dtcld,qrs(i,k,1)) 
    856. 

  856. !---------------------------------------------------------------
    857. 

  857. ! nsfrz: freezing of rain water [LH A26]
    858. 

  858. !        (T<T0, NR-> )
    859. 

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

  860.             if(ncr(i,k,3).gt.nrmin) then
    861. 

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

  862.                        *rslope3(i,k,1)*dtcld,ncr(i,k,3))
    863. 

  863.               ncr(i,k,3) = ncr(i,k,3)-nfrzdtr
    864. 

  864.             endif
    865. 

  865.             qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
    866. 

  866.             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
    867. 

  867.             qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
    868. 

  868.           endif
    869. 

  869.         enddo
    870. 

  870.       enddo
    871. 

  871. !
    872. 

  872.       do k = kts, kte
    873. 

  873.         do i = its, ite
    874. 

  874.           ncr(i,k,2) = max(ncr(i,k,2),0.0)
    875. 

  875.           ncr(i,k,3) = max(ncr(i,k,3),0.0)
    876. 

  876.         enddo
    877. 

  877.       enddo
    878. 

  878. !----------------------------------------------------------------
    879. 

  879. !     update the slope parameters for microphysics computation
    880. 

  880. !
    881. 

  881.       do k = kts, kte
    882. 

  882.         do i = its, ite
    883. 

  883.           qrs_tmp(i,k,1) = qrs(i,k,1)
    884. 

  884.           qrs_tmp(i,k,2) = qrs(i,k,2)
    885. 

  885.           ncr_tmp(i,k) = ncr(i,k,3)
    886. 

  886.         enddo
    887. 

  887.       enddo
    888. 

  888.       call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
    889. 

  889.                      rslope3,work1,workn,its,ite,kts,kte)
    890. 

  890.       do k = kts, kte
    891. 

  891.         do i = its, ite
    892. 

  892. !-----------------------------------------------------------------
    893. 

  893. ! compute the mean-volume drop diameter                  [LH A10]
    894. 

  894. ! for raindrop distribution
    895. 

  895. !-----------------------------------------------------------------
    896. 

  896.           avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
    897. 

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

  898.             rslopec(i,k) = rslopecmax
    899. 

  899.             rslopec2(i,k) = rslopec2max
    900. 

  900.             rslopec3(i,k) = rslopec3max
    901. 

  901.           else
    902. 

  902.             rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
    903. 

  903.             rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
    904. 

  904.             rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
    905. 

  905.           endif
    906. 

  906. !--------------------------------------------------------------------
    907. 

  907. ! compute the mean-volume drop diameter                   [LH A7]
    908. 

  908. ! for cloud-droplet distribution
    909. 

  909. !--------------------------------------------------------------------
    910. 

  910.           avedia(i,k,1) = rslopec(i,k)
    911. 

  911.         enddo
    912. 

  912.       enddo
    913. 

  913. !----------------------------------------------------------------
    914. 

  914. !     work1:  the thermodynamic term in the denominator associated with
    915. 

  915. !             heat conduction and vapor diffusion
    916. 

  916. !             (ry88, y93, h85)
    917. 

  917. !     work2: parameter associated with the ventilation effects(y93)
    918. 

  918. !
    919. 

  919.       do k = kts, kte
    920. 

  920.         do i = its, ite
    921. 

  921. !         work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
    922. 

  922.           work1(i,k,1) = ((((den(i,k))*(xl(i,k))*(xl(i,k)))*((t(i,k))+120.)    &                          
    923. 

  923.                        *(den(i,k)))/(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))&
    924. 

  924.                        *(den(i,k))*(rv*(t(i,k))*(t(i,k)))))                    &
    925. 

  925.                       +  p(i,k)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k))*(1.81))))
    926. 

  926. !         work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
    927. 

  927.           work1(i,k,2) = ((((den(i,k))*(xls)*(xls))*((t(i,k))+120.)*(den(i,k)))&
    928. 

  928.                        /(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))*(den(i,k)) & 
    929. 

  929.                        *(rv*(t(i,k))*(t(i,k))))                                & 
    930. 

  930.                       + p(i,k)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k))*(1.81)))))
    931. 

  931. !         work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
    932. 

  932.           work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) &
    933. 

  933.                      *p(i,k))/(((t(i,k))+120.)*den(i,k)*(8.794e-5              &
    934. 

  934.                      *exp(log(t(i,k))*(1.81))))))*sqrt(sqrt(den0/(den(i,k))))) & 
    935. 

  935.                       /sqrt((1.496e-6*((t(i,k))*sqrt(t(i,k))))                 &
    936. 

  936.                       /(((t(i,k))+120.)*den(i,k)))
    937. 

  937.         enddo 
    938. 

  938.       enddo 
    939. 

  939. !
    940. 

  940. !===============================================================
    941. 

  941. !
    942. 

  942. ! warm rain processes
    943. 

  943. !
    944. 

  944. ! - follows the processes in RH83 and LFO except for autoconcersion
    945. 

  945. !
    946. 

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

  947. !
    948. 

  948.       do k = kts, kte
    949. 

  949.         do i = its, ite
    950. 

  950.           supsat = max(q(i,k),qmin)-qs(i,k,1)
    951. 

  951.           satdt = supsat/dtcld
    952. 

  952. !---------------------------------------------------------------
    953. 

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

  954. !        (QC->QR)
    955. 

  955. !---------------------------------------------------------------
    956. 

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

  957.                    *rslopec2(i,k)-0.4)
    958. 

  958.           lenconcr = max(1.2*lencon,qcrmin)
    959. 

  959.           if(avedia(i,k,1).gt.di15) then
    960. 

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

  961.             taucon = max(taucon, 1.e-12)
    962. 

  962.             praut(i,k) = lencon/(taucon*den(i,k))
    963. 

  963.             praut(i,k) = min(max(praut(i,k),0.),qci(i,k,1)/dtcld)
    964. 

  964. !---------------------------------------------------------------
    965. 

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

  966. !        (NC->NR)
    967. 

  967. !---------------------------------------------------------------
    968. 

  968.             nraut(i,k) = 3.5e9*den(i,k)*praut(i,k)
    969. 

  969.             if(qrs(i,k,1).gt.lenconcr)                                         &
    970. 

  970.             nraut(i,k) = ncr(i,k,3)/qrs(i,k,1)*praut(i,k)
    971. 

  971.             nraut(i,k) = min(nraut(i,k),ncr(i,k,2)/dtcld)
    972. 

  972.           endif
    973. 

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

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

  975. !        (QC->QR)
    976. 

  976. ! nracw: accretion of cloud water by rain   [LH A9]
    977. 

  977. !        (NC->)
    978. 

  978. !---------------------------------------------------------------
    979. 

  979.           if(qrs(i,k,1).ge.lenconcr) then
    980. 

  980.             if(avedia(i,k,2).ge.di100) then
    981. 

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

  982.                          + 24.*rslope3(i,k,1)),ncr(i,k,2)/dtcld)
    983. 

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

  984.                          *ncr(i,k,3)*rslopec3(i,k)*(2.*rslopec3(i,k)           &
    985. 

  985.                          + 24.*rslope3(i,k,1)),qci(i,k,1)/dtcld)
    986. 

  986.             else
    987. 

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

  988.                          *rslopec3(i,k)+5040.*rslope3(i,k,1)                   &
    989. 

  989.                          *rslope3(i,k,1)),ncr(i,k,2)/dtcld)
    990. 

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

  991.                          *ncr(i,k,3)*rslopec3(i,k)*(6.*rslopec3(i,k)           &
    992. 

  992.                          *rslopec3(i,k)+5040.*rslope3(i,k,1)                   &
    993. 

  993.                          *rslope3(i,k,1)),qci(i,k,1)/dtcld)
    994. 

  994.             endif
    995. 

  995.           endif
    996. 

  996. !----------------------------------------------------------------
    997. 

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

  998. !        (NC->)
    999. 

  999. !----------------------------------------------------------------
    1000. 

  1000.           if(avedia(i,k,1).ge.di100) then
    1001. 

  1001.             nccol(i,k) = ncrk1*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)
    1002. 

  1002.           else
    1003. 

  1003.             nccol(i,k) = 2.*ncrk2*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)        & 
    1004. 

  1004.                         *rslopec3(i,k)
    1005. 

  1005.           endif
    1006. 

  1006. !----------------------------------------------------------------
    1007. 

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

  1008. !        (NR->)
    1009. 

  1009. !----------------------------------------------------------------
    1010. 

  1010.           if(qrs(i,k,1).ge.lenconcr) then
    1011. 

  1011.             if(avedia(i,k,2).lt.di100) then
    1012. 

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

  1013.                           *rslope3(i,k,1)
    1014. 

  1014.             elseif(avedia(i,k,2).ge.di100 .and. avedia(i,k,2).lt.di600) then
    1015. 

  1015.               nrcol(i,k) = 24.*ncrk1*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1)
    1016. 

  1016.             elseif(avedia(i,k,2).ge.di600 .and. avedia(i,k,2).lt.di2000) then
    1017. 

  1017.               coecol = -2.5e3*(avedia(i,k,2)-di600)
    1018. 

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

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

  1020.             else
    1021. 

  1021.               nrcol(i,k) = 0.
    1022. 

  1022.             endif
    1023. 

  1023.           endif
    1024. 

  1024. !---------------------------------------------------------------
    1025. 

  1025. ! prevp: evaporation/condensation rate of rain  [HL A41]
    1026. 

  1026. !        (QV->QR or QR->QV)
    1027. 

  1027. !---------------------------------------------------------------
    1028. 

  1028.           if(qrs(i,k,1).gt.0.) then
    1029. 

  1029.             coeres = rslope(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
    1030. 

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

  1031.                          +precr2*work2(i,k)*coeres)/work1(i,k,1)
    1032. 

  1032.             if(prevp(i,k).lt.0.) then
    1033. 

  1033.               prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
    1034. 

  1034.               prevp(i,k) = max(prevp(i,k),satdt/2)
    1035. 

  1035. !----------------------------------------------------------------
    1036. 

  1036. ! Nrevp: evaporation/condensation rate of rain  [LH A14] 
    1037. 

  1037. !        (NR->NCCN)
    1038. 

  1038. !----------------------------------------------------------------
    1039. 

  1039.               if(prevp(i,k).eq.-qrs(i,k,1)/dtcld) then
    1040. 

  1040.                  ncr(i,k,1) = ncr(i,k,1) + ncr(i,k,3)
    1041. 

  1041.                  ncr(i,k,3) = 0.
    1042. 

  1042.               endif
    1043. 

  1043.             else
    1044. 

  1044. !
    1045. 

  1045.               prevp(i,k) = min(prevp(i,k),satdt/2)
    1046. 

  1046.             endif
    1047. 

  1047.           endif
    1048. 

  1048.         enddo
    1049. 

  1049.       enddo
    1050. 

  1050. !
    1051. 

  1051. !===============================================================
    1052. 

  1052. !
    1053. 

  1053. ! cold rain processes
    1054. 

  1054. !
    1055. 

  1055. ! - follows the revised ice microphysics processes in HDC
    1056. 

  1056. ! - the processes same as in RH83 and RH84  and LFO behave
    1057. 

  1057. !   following ice crystal hapits defined in HDC, inclduing
    1058. 

  1058. !   intercept parameter for snow (n0s), ice crystal number
    1059. 

  1059. !   concentration (ni), ice nuclei number concentration
    1060. 

  1060. !   (n0i), ice diameter (d)
    1061. 

  1061. !
    1062. 

  1062. !===============================================================
    1063. 

  1063. !
    1064. 

  1064.       rdtcld = 1./dtcld
    1065. 

  1065.       do k = kts, kte
    1066. 

  1066.         do i = its, ite
    1067. 

  1067.           supcol = t0c-t(i,k)
    1068. 

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

  1069.           supsat = max(q(i,k),qmin)-qs(i,k,2)
    1070. 

  1070.           satdt = supsat/dtcld
    1071. 

  1071.           ifsat = 0
    1072. 

  1072. !-------------------------------------------------------------
    1073. 

  1073. ! Ni: ice crystal number concentraiton   [HDC 5c]
    1074. 

  1074. !-------------------------------------------------------------
    1075. 

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

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

  1077.           temp = (den(i,k)*max(qci(i,k,2),qmin))
    1078. 

  1078.           temp = sqrt(sqrt(temp*temp*temp))
    1079. 

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

  1080.           eacrs = exp(0.07*(-supcol))
    1081. 

  1081. !
    1082. 

  1082.           if(supcol.gt.0) then
    1083. 

  1083.             if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,2).gt.qmin) then
    1084. 

  1084.               xmi = den(i,k)*qci(i,k,2)/xni(i,k)
    1085. 

  1085.               diameter  = min(dicon * sqrt(xmi),dimax)
    1086. 

  1086.               vt2i = 1.49e4*diameter**1.31
    1087. 

  1087.               vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
    1088. 

  1088. !-------------------------------------------------------------
    1089. 

  1089. ! psaci: Accretion of cloud ice by rain [HDC 10]
    1090. 

  1090. !        (T<T0: QI->QS)
    1091. 

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

  1092.               acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2)            &
    1093. 

  1093.                       + diameter**2*rslope(i,k,2)
    1094. 

  1094.               psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k)*abs(vt2s-vt2i)  &
    1095. 

  1095.                          *acrfac/4.
    1096. 

  1096.             endif
    1097. 

  1097.           endif
    1098. 

  1098. !-------------------------------------------------------------
    1099. 

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

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

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

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

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

  1104.                         *qci(i,k,1)*denfac(i,k),qci(i,k,1)*rdtcld)
    1105. 

  1105.           endif
    1106. 

  1106. !-------------------------------------------------------------
    1107. 

  1107. ! nsacw: Accretion of cloud water by snow  [LH A12]
    1108. 

  1108. !        (NC ->)
    1109. 

  1109. !-------------------------------------------------------------
    1110. 

  1110.          if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
    1111. 

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

  1112.                        *ncr(i,k,2)*denfac(i,k),ncr(i,k,2)/dtcld)
    1113. 

  1113.          endif
    1114. 

  1114.          if(supcol.le.0) then
    1115. 

  1115.            xlf = xlf0
    1116. 

  1116. !--------------------------------------------------------------
    1117. 

  1117. ! nseml: Enhanced melting of snow by accretion of water  [LH A29]
    1118. 

  1118. !        (T>=T0: ->NR)
    1119. 

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

  1120.               if  (qrs(i,k,2).gt.qcrmin) then
    1121. 

  1121.                 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
    1122. 

  1122.                 nseml(i,k) = -sfac*min(max(cliq*supcol*(psacw(i,k))/xlf        &
    1123. 

  1123.                             ,-qrs(i,k,2)/dtcld),0.)
    1124. 

  1124.               endif   
    1125. 

  1125.          endif
    1126. 

  1126.          if(supcol.gt.0) then
    1127. 

  1127. !-------------------------------------------------------------
    1128. 

  1128. ! pidep: Deposition/Sublimation rate of ice [HDC 9]
    1129. 

  1129. !       (T<T0: QV->QI or QI->QV)
    1130. 

  1130. !-------------------------------------------------------------
    1131. 

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

  1132.               xmi = den(i,k)*qci(i,k,2)/xni(i,k)
    1133. 

  1133.               diameter = dicon * sqrt(xmi)
    1134. 

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

  1135.               supice = satdt-prevp(i,k)
    1136. 

  1136.               if(pidep(i,k).lt.0.) then
    1137. 

  1137. !               pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
    1138. 

  1138. !               pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
    1139. 

  1139.                 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
    1140. 

  1140.                 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
    1141. 

  1141.               else
    1142. 

  1142. !               pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
    1143. 

  1143.                 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
    1144. 

  1144.               endif
    1145. 

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

  1146.             endif
    1147. 

  1147. !-------------------------------------------------------------
    1148. 

  1148. ! psdep: deposition/sublimation rate of snow [HDC 14]
    1149. 

  1149. !        (QV->QS or QS->QV)
    1150. 

  1150. !-------------------------------------------------------------
    1151. 

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

  1152.               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    1153. 

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

  1154.                            + precs2*work2(i,k)*coeres)/work1(i,k,2)
    1155. 

  1155.               supice = satdt-prevp(i,k)-pidep(i,k)
    1156. 

  1156.               if(psdep(i,k).lt.0.) then
    1157. 

  1157. !               psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
    1158. 

  1158. !               psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
    1159. 

  1159.                 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
    1160. 

  1160.                 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
    1161. 

  1161.               else
    1162. 

  1162. !             psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
    1163. 

  1163.                 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
    1164. 

  1164.               endif
    1165. 

  1165.               if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) ifsat = 1
    1166. 

  1166.             endif
    1167. 

  1167. !-------------------------------------------------------------
    1168. 

  1168. ! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
    1169. 

  1169. !       (T<T0: QV->QI)
    1170. 

  1170. !-------------------------------------------------------------
    1171. 

  1171.             if(supsat.gt.0 .and. ifsat.ne.1) then
    1172. 

  1172.               supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
    1173. 

  1173.               xni0 = 1.e3*exp(0.1*supcol)
    1174. 

  1174.               roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
    1175. 

  1175.               pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))*rdtcld)
    1176. 

  1176.               pigen(i,k) = min(min(pigen(i,k),satdt),supice)
    1177. 

  1177.             endif
    1178. 

  1178. !
    1179. 

  1179. !-------------------------------------------------------------
    1180. 

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

  1181. !       (T<T0: QI->QS)
    1182. 

  1182. !-------------------------------------------------------------
    1183. 

  1183.             if(qci(i,k,2).gt.0.) then
    1184. 

  1184.               qimax = roqimax/den(i,k)
    1185. 

  1185. !             psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
    1186. 

  1186.               psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
    1187. 

  1187.             endif
    1188. 

  1188.           endif
    1189. 

  1189. !-------------------------------------------------------------
    1190. 

  1190. ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
    1191. 

  1191. !       (T>T0: QS->QV)
    1192. 

  1192. !-------------------------------------------------------------
    1193. 

  1193.           if(supcol.lt.0.) then
    1194. 

  1194.             if(qrs(i,k,2).gt.0. .and. rh(i,k,1).lt.1.)                         &
    1195. 

  1195.               psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
    1196. 

  1196. !              psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
    1197. 

  1197.               psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
    1198. 

  1198.           endif
    1199. 

  1199.         enddo
    1200. 

  1200.       enddo
    1201. 

  1201. !
    1202. 

  1202. !
    1203. 

  1203. !----------------------------------------------------------------
    1204. 

  1204. !     check mass conservation of generation terms and feedback to the
    1205. 

  1205. !     large scale
    1206. 

  1206. !
    1207. 

  1207.       do k = kts, kte
    1208. 

  1208.         do i = its, ite
    1209. 

  1209.           if(t(i,k).le.t0c) then
    1210. 

  1210. !
    1211. 

  1211. !     Q_cloud water
    1212. 

  1212. !
    1213. 

  1213.             value = max(qmin,qci(i,k,1))
    1214. 

  1214.             source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
    1215. 

  1215.             if (source.gt.value) then
    1216. 

  1216.               factor = value/source
    1217. 

  1217.               praut(i,k) = praut(i,k)*factor
    1218. 

  1218.               pracw(i,k) = pracw(i,k)*factor
    1219. 

  1219.               psacw(i,k) = psacw(i,k)*factor
    1220. 

  1220.             endif
    1221. 

  1221. !
    1222. 

  1222. !     Q_cloud ice
    1223. 

  1223. !
    1224. 

  1224.             value = max(qmin,qci(i,k,2))
    1225. 

  1225.             source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
    1226. 

  1226.             if (source.gt.value) then
    1227. 

  1227.               factor = value/source
    1228. 

  1228.               psaut(i,k) = psaut(i,k)*factor
    1229. 

  1229.               psaci(i,k) = psaci(i,k)*factor
    1230. 

  1230.               pigen(i,k) = pigen(i,k)*factor
    1231. 

  1231.               pidep(i,k) = pidep(i,k)*factor
    1232. 

  1232.             endif
    1233. 

  1233. !
    1234. 

  1234. !     Q_rain
    1235. 

  1235. !
    1236. 

  1236. !
    1237. 

  1237.             value = max(qmin,qrs(i,k,1))
    1238. 

  1238.             source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
    1239. 

  1239.             if (source.gt.value) then
    1240. 

  1240.               factor = value/source
    1241. 

  1241.               praut(i,k) = praut(i,k)*factor
    1242. 

  1242.               pracw(i,k) = pracw(i,k)*factor
    1243. 

  1243.               prevp(i,k) = prevp(i,k)*factor
    1244. 

  1244.             endif
    1245. 

  1245. !
    1246. 

  1246. !    Q_snow
    1247. 

  1247. !
    1248. 

  1248.             value = max(qmin,qrs(i,k,2))
    1249. 

  1249.             source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld  
    1250. 

  1250.             if (source.gt.value) then
    1251. 

  1251.               factor = value/source
    1252. 

  1252.               psdep(i,k) = psdep(i,k)*factor
    1253. 

  1253.               psaut(i,k) = psaut(i,k)*factor
    1254. 

  1254.               psaci(i,k) = psaci(i,k)*factor
    1255. 

  1255.               psacw(i,k) = psacw(i,k)*factor
    1256. 

  1256.             endif
    1257. 

  1257. !
    1258. 

  1258. !     N_cloud
    1259. 

  1259. !
    1260. 

  1260.             value = max(ncmin,ncr(i,k,2))
    1261. 

  1261.             source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+nsacw(i,k))*dtcld
    1262. 

  1262.             if (source.gt.value) then
    1263. 

  1263.               factor = value/source
    1264. 

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

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

  1266.               nracw(i,k) = nracw(i,k)*factor
    1267. 

  1267.               nsacw(i,k) = nsacw(i,k)*factor
    1268. 

  1268.             endif
    1269. 

  1269. !
    1270. 

  1270. !     N_rain
    1271. 

  1271. !
    1272. 

  1272.             value = max(nrmin,ncr(i,k,3))
    1273. 

  1273.             source = (-nraut(i,k)+nrcol(i,k))*dtcld
    1274. 

  1274.             if (source.gt.value) then
    1275. 

  1275.               factor = value/source
    1276. 

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

  1277.               nrcol(i,k) = nrcol(i,k)*factor
    1278. 

  1278.             endif
    1279. 

  1279. !
    1280. 

  1280.             work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
    1281. 

  1281. !     update
    1282. 

  1282.             q(i,k) = q(i,k)+work2(i,k)*dtcld
    1283. 

  1283.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)+psacw(i,k)      &  
    1284. 

  1284.                         )*dtcld,0.)
    1285. 

  1285.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)+prevp(i,k)      &    
    1286. 

  1286.                         )*dtcld,0.)
    1287. 

  1287.             qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k)-pigen(i,k)      & 
    1288. 

  1288.                         -pidep(i,k))*dtcld,0.)
    1289. 

  1289.             qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+psaci(i,k)      &
    1290. 

  1290.                         +psacw(i,k))*dtcld,0.)
    1291. 

  1291.             ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k)     &
    1292. 

  1292.                          -nsacw(i,k))*dtcld,0.)
    1293. 

  1293.             ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k))                &
    1294. 

  1294.                         *dtcld,0.)
    1295. 

  1295.             xlf = xls-xl(i,k)
    1296. 

  1296.             xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k))                  &
    1297. 

  1297.                       -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
    1298. 

  1298.             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
    1299. 

  1299.           else
    1300. 

  1300. !
    1301. 

  1301. !     Q_cloud water
    1302. 

  1302. !
    1303. 

  1303.             value = max(qmin,qci(i,k,1))
    1304. 

  1304.             source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
    1305. 

  1305.             if (source.gt.value) then
    1306. 

  1306.               factor = value/source
    1307. 

  1307.               praut(i,k) = praut(i,k)*factor
    1308. 

  1308.               pracw(i,k) = pracw(i,k)*factor
    1309. 

  1309.               psacw(i,k) = psacw(i,k)*factor
    1310. 

  1310.             endif
    1311. 

  1311. !
    1312. 

  1312. !     Q_rain
    1313. 

  1313. !
    1314. 

  1314.             value = max(qmin,qrs(i,k,1))
    1315. 

  1315.             source = (-praut(i,k)-pracw(i,k)-prevp(i,k)                        &         
    1316. 

  1316.                     -psacw(i,k))*dtcld
    1317. 

  1317.             if (source.gt.value) then
    1318. 

  1318.               factor = value/source
    1319. 

  1319.               praut(i,k) = praut(i,k)*factor
    1320. 

  1320.               pracw(i,k) = pracw(i,k)*factor
    1321. 

  1321.               prevp(i,k) = prevp(i,k)*factor
    1322. 

  1322.               psacw(i,k) = psacw(i,k)*factor
    1323. 

  1323.             endif  
    1324. 

  1324. !
    1325. 

  1325. !     Q_snow
    1326. 

  1326. !
    1327. 

  1327.             value = max(qcrmin,qrs(i,k,2))
    1328. 

  1328.             source=(-psevp(i,k))*dtcld
    1329. 

  1329.             if (source.gt.value) then
    1330. 

  1330.               factor = value/source
    1331. 

  1331.               psevp(i,k) = psevp(i,k)*factor
    1332. 

  1332.             endif
    1333. 

  1333. !
    1334. 

  1334. !     N_cloud
    1335. 

  1335. !
    1336. 

  1336.             value = max(ncmin,ncr(i,k,2))
    1337. 

  1337.             source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+nsacw(i,k))*dtcld
    1338. 

  1338.             if (source.gt.value) then
    1339. 

  1339.               factor = value/source
    1340. 

  1340.               nraut(i,k) = nraut(i,k)*factor
    1341. 

  1341.               nccol(i,k) = nccol(i,k)*factor
    1342. 

  1342.               nracw(i,k) = nracw(i,k)*factor
    1343. 

  1343.               nsacw(i,k) = nsacw(i,k)*factor
    1344. 

  1344.             endif
    1345. 

  1345. !
    1346. 

  1346. !     N_rain
    1347. 

  1347. !
    1348. 

  1348.             value = max(nrmin,ncr(i,k,3))
    1349. 

  1349.             source = (-nraut(i,k)-nseml(i,k)+nrcol(i,k))*dtcld
    1350. 

  1350.             if (source.gt.value) then
    1351. 

  1351.               factor = value/source
    1352. 

  1352.               nraut(i,k) = nraut(i,k)*factor
    1353. 

  1353.               nseml(i,k) = nseml(i,k)*factor
    1354. 

  1354.               nrcol(i,k) = nrcol(i,k)*factor
    1355. 

  1355.             endif
    1356. 

  1356.             work2(i,k)=-(prevp(i,k)+psevp(i,k))
    1357. 

  1357. !     update
    1358. 

  1358.             q(i,k) = q(i,k)+work2(i,k)*dtcld
    1359. 

  1359.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)+psacw(i,k)      &         
    1360. 

  1360.                         )*dtcld,0.)
    1361. 

  1361.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)+prevp(i,k)      & 
    1362. 

  1362.                         +psacw(i,k))*dtcld,0.)
    1363. 

  1363.             qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
    1364. 

  1364.             ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k)     &
    1365. 

  1365.                         -nsacw(i,k))*dtcld,0.)
    1366. 

  1366.             ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)+nseml(i,k)-nrcol(i,k)      &
    1367. 

  1367.                         )*dtcld,0.)
    1368. 

  1368.             xlf = xls-xl(i,k)
    1369. 

  1369.             xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
    1370. 

  1370.             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
    1371. 

  1371.           endif
    1372. 

  1372.         enddo
    1373. 

  1373.       enddo
    1374. 

  1374. !
    1375. 

  1375. ! Inline expansion for fpvs
    1376. 

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

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

  1378.       hsub = xls
    1379. 

  1379.       hvap = xlv0
    1380. 

  1380.       cvap = cpv
    1381. 

  1381.       ttp=t0c+0.01
    1382. 

  1382.       dldt=cvap-cliq
    1383. 

  1383.       xa=-dldt/rv
    1384. 

  1384.       xb=xa+hvap/(rv*ttp)
    1385. 

  1385.       dldti=cvap-cice
    1386. 

  1386.       xai=-dldti/rv
    1387. 

  1387.       xbi=xai+hsub/(rv*ttp)
    1388. 

  1388.       do k = kts, kte
    1389. 

  1389.       do i = its, ite
    1390. 

  1390.           tr=ttp/t(i,k)
    1391. 

  1391.           logtr = log(tr)
    1392. 

  1392.           qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
    1393. 

  1393.           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
    1394. 

  1394.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    1395. 

  1395.           qs(i,k,1) = max(qs(i,k,1),qmin)
    1396. 

  1396.           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
    1397. 

  1397.         enddo
    1398. 

  1398.       enddo
    1399. 

  1399. !
    1400. 

  1400.       call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
    1401. 

  1401.                      rslope3,work1,workn,its,ite,kts,kte)
    1402. 

  1402.       do k = kts, kte
    1403. 

  1403.         do i = its, ite
    1404. 

  1404. !-----------------------------------------------------------------
    1405. 

  1405. ! re-compute the mean-volume drop diameter            [LH A10]
    1406. 

  1406. ! for raindrop distribution
    1407. 

  1407. !-----------------------------------------------------------------
    1408. 

  1408.           avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
    1409. 

  1409. !----------------------------------------------------------------
    1410. 

  1410. ! Nrevp_s: evaporation/condensation rate of rain   [LH A14]
    1411. 

  1411. !        (NR->NC)
    1412. 

  1412. !----------------------------------------------------------------
    1413. 

  1413.           if(avedia(i,k,2).le.di82) then
    1414. 

  1414.             ncr(i,k,2) = ncr(i,k,2)+ncr(i,k,3)
    1415. 

  1415.             ncr(i,k,3) = 0.
    1416. 

  1416. !----------------------------------------------------------------
    1417. 

  1417. ! Prevp_s: evaporation/condensation rate of rain [LH A15] [KK 23]
    1418. 

  1418. !        (QR->QC)
    1419. 

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

  1420.             qci(i,k,1) = qci(i,k,1)+qrs(i,k,1)
    1421. 

  1421.             qrs(i,k,1) = 0.
    1422. 

  1422.           endif
    1423. 

  1423.         enddo
    1424. 

  1424.       enddo
    1425. 

  1425. !
    1426. 

  1426.       do k = kts, kte
    1427. 

  1427.         do i = its, ite
    1428. 

  1428. !-------------------------------------------------------------------
    1429. 

  1429. ! rate of change of cloud drop concentration due to CCN activation
    1430. 

  1430. ! pcact: QV -> QC                                 [LH 8]  [KK 14]
    1431. 

  1431. ! ncact: NCCN -> NC                               [LH A2] [KK 12]
    1432. 

  1432. !-------------------------------------------------------------------
    1433. 

  1433.           if(rh(i,k,1).gt.1.) then
    1434. 

  1434.             ncact(i,k) = max(0.,((ncr(i,k,1)+ncr(i,k,2))                       &
    1435. 

  1435.                        *min(1.,(rh(i,k,1)/satmax)**actk) - ncr(i,k,2)))/dtcld
    1436. 

  1436.             ncact(i,k) =min(ncact(i,k),max(ncr(i,k,1),0.)/dtcld)
    1437. 

  1437.             pcact(i,k) = min(4.*pi*denr*(actr*1.E-6)**3*ncact(i,k)/            &
    1438. 

  1438.                          (3.*den(i,k)),max(q(i,k),0.)/dtcld)
    1439. 

  1439.             q(i,k) = max(q(i,k)-pcact(i,k)*dtcld,0.)
    1440. 

  1440.             qci(i,k,1) = max(qci(i,k,1)+pcact(i,k)*dtcld,0.)
    1441. 

  1441.             ncr(i,k,1) = max(ncr(i,k,1)-ncact(i,k)*dtcld,0.)
    1442. 

  1442.             ncr(i,k,2) = max(ncr(i,k,2)+ncact(i,k)*dtcld,0.)
    1443. 

  1443.             t(i,k) = t(i,k)+pcact(i,k)*xl(i,k)/cpm(i,k)*dtcld
    1444. 

  1444.           endif
    1445. 

  1445. !---------------------------------------------------------------------
    1446. 

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

  1447. !     if there exists additional water vapor condensated/if
    1448. 

  1448. !     evaporation of cloud water is not enough to remove subsaturation
    1449. 

  1449. !  (QV->QC or QC->QV)
    1450. 

  1450. !---------------------------------------------------------------------
    1451. 

  1451.           tr=ttp/t(i,k)
    1452. 

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

  1453.           qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
    1454. 

  1454.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    1455. 

  1455.           qs(i,k,1) = max(qs(i,k,1),qmin)
    1456. 

  1456.           work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/(1.+(xl(i,k))         &
    1457. 

  1457.                        *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))/((t(i,k))        &
    1458. 

  1458.                        *(t(i,k)))))
    1459. 

  1459.           work2(i,k) = qci(i,k,1)+work1(i,k,1)
    1460. 

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

  1461.           if(qci(i,k,1).gt.0. .and. work1(i,k,1).lt.0.)                        &
    1462. 

  1462.             pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
    1463. 

  1463. !----------------------------------------------------------------------
    1464. 

  1464. ! ncevp: evpration of Cloud number concentration  [LH A3]
    1465. 

  1465. ! (NC->NCCN)
    1466. 

  1466. !----------------------------------------------------------------------
    1467. 

  1467.           if(pcond(i,k).eq.-qci(i,k,1)/dtcld) then
    1468. 

  1468.             ncr(i,k,2) = 0.
    1469. 

  1469.             ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,2)
    1470. 

  1470.           endif
    1471. 

  1471. !
    1472. 

  1472.           q(i,k) = q(i,k)-pcond(i,k)*dtcld
    1473. 

  1473.           qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
    1474. 

  1474.           t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
    1475. 

  1475.         enddo
    1476. 

  1476.       enddo
    1477. 

  1477. !
    1478. 

  1478. !----------------------------------------------------------------
    1479. 

  1479. !     padding for small values
    1480. 

  1480. !
    1481. 

  1481.       do k = kts, kte
    1482. 

  1482.         do i = its, ite
    1483. 

  1483.           if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
    1484. 

  1484.           if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
    1485. 

  1485.         enddo
    1486. 

  1486.       enddo
    1487. 

  1487.       enddo                  ! big loops
    1488. 

  1488.   END SUBROUTINE wdm52d
    1489. 

  1489. ! ...................................................................
    1490. 

  1490.       REAL FUNCTION rgmma(x)
    1491. 

  1491. !-------------------------------------------------------------------
    1492. 

  1492.   IMPLICIT NONE
    1493. 

  1493. !-------------------------------------------------------------------
    1494. 

  1494. !     rgmma function:  use infinite product form
    1495. 

  1495.       REAL :: euler
    1496. 

  1496.       PARAMETER (euler=0.577215664901532)
    1497. 

  1497.       REAL :: x, y
    1498. 

  1498.       INTEGER :: i
    1499. 

  1499.       if(x.eq.1.)then
    1500. 

  1500.         rgmma=0.
    1501. 

  1501.           else
    1502. 

  1502.         rgmma=x*exp(euler*x)
    1503. 

  1503.         do i=1,10000
    1504. 

  1504.           y=float(i)
    1505. 

  1505.           rgmma=rgmma*(1.000+x/y)*exp(-x/y)
    1506. 

  1506.         enddo
    1507. 

  1507.         rgmma=1./rgmma
    1508. 

  1508.       endif
    1509. 

  1509.       END FUNCTION rgmma
    1510. 

  1510. !
    1511. 

  1511. !--------------------------------------------------------------------------
    1512. 

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

  1513. !--------------------------------------------------------------------------
    1514. 

  1514.       IMPLICIT NONE
    1515. 

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

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

  1517.            xai,xbi,ttp,tr
    1518. 

  1518.       INTEGER ice
    1519. 

  1519. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1520. 

  1520.       ttp=t0c+0.01
    1521. 

  1521.       dldt=cvap-cliq
    1522. 

  1522.       xa=-dldt/rv
    1523. 

  1523.       xb=xa+hvap/(rv*ttp)
    1524. 

  1524.       dldti=cvap-cice
    1525. 

  1525.       xai=-dldti/rv
    1526. 

  1526.       xbi=xai+hsub/(rv*ttp)
    1527. 

  1527.       tr=ttp/t
    1528. 

  1528.       if(t.lt.ttp .and. ice.eq.1) then
    1529. 

  1529.         fpvs=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
    1530. 

  1530.       else
    1531. 

  1531.         fpvs=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
    1532. 

  1532.       endif
    1533. 

  1533. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1534. 

  1534.       END FUNCTION fpvs
    1535. 

  1535. !-------------------------------------------------------------------
    1536. 

  1536.   SUBROUTINE wdm5init(den0,denr,dens,cl,cpv,ccn0,allowed_to_read)
    1537. 

  1537. !-------------------------------------------------------------------
    1538. 

  1538.   IMPLICIT NONE
    1539. 

  1539. !-------------------------------------------------------------------
    1540. 

  1540. !.... constants which may not be tunable
    1541. 

  1541.    REAL, INTENT(IN) :: den0,denr,dens,cl,cpv,ccn0
    1542. 

  1542.    LOGICAL, INTENT(IN) :: allowed_to_read
    1543. 

  1543. !
    1544. 

  1544.    pi = 4.*atan(1.)
    1545. 

  1545.    xlv1 = cl-cpv
    1546. 

  1546. !
    1547. 

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

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

  1549.    pidnc = pi*denr/6.
    1550. 

  1550. !
    1551. 

  1551.    bvtr1 = 1.+bvtr
    1552. 

  1552.    bvtr2 = 2.+bvtr
    1553. 

  1553.    bvtr3 = 3.+bvtr
    1554. 

  1554.    bvtr4 = 4.+bvtr
    1555. 

  1555.    bvtr5 = 5.+bvtr
    1556. 

  1556.    bvtr7 = 7.+bvtr
    1557. 

  1557.    bvtr2o5 = 2.5+.5*bvtr
    1558. 

  1558.    bvtr3o5 = 3.5+.5*bvtr
    1559. 

  1559.    g1pbr = rgmma(bvtr1)
    1560. 

  1560.    g2pbr = rgmma(bvtr2)
    1561. 

  1561.    g3pbr = rgmma(bvtr3)
    1562. 

  1562.    g4pbr = rgmma(bvtr4)            ! 17.837825
    1563. 

  1563.    g5pbr = rgmma(bvtr5)
    1564. 

  1564.    g7pbr = rgmma(bvtr7)
    1565. 

  1565.    g5pbro2 = rgmma(bvtr2o5)
    1566. 

  1566.    g7pbro2 = rgmma(bvtr3o5)
    1567. 

  1567.    pvtr = avtr*g5pbr/24.
    1568. 

  1568.    pvtrn = avtr*g2pbr
    1569. 

  1569.    eacrr = 1.0
    1570. 

  1570.    pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
    1571. 

  1571.    precr1 = 2.*pi*1.56
    1572. 

  1572.    precr2 = 2.*pi*.31*avtr**.5*g7pbro2
    1573. 

  1573.    pidn0r =  pi*denr*n0r
    1574. 

  1574.    pidnr = 4.*pi*denr
    1575. 

  1575.    xmmax = (dimax/dicon)**2
    1576. 

  1576.    roqimax = 2.08e22*dimax**8
    1577. 

  1577. !
    1578. 

  1578.    bvts1 = 1.+bvts
    1579. 

  1579.    bvts2 = 2.5+.5*bvts
    1580. 

  1580.    bvts3 = 3.+bvts
    1581. 

  1581.    bvts4 = 4.+bvts
    1582. 

  1582.    g1pbs = rgmma(bvts1)    !.8875
    1583. 

  1583.    g3pbs = rgmma(bvts3)
    1584. 

  1584.    g4pbs = rgmma(bvts4)    ! 12.0786
    1585. 

  1585.    g5pbso2 = rgmma(bvts2)
    1586. 

  1586.    pvts = avts*g4pbs/6.
    1587. 

  1587.    pacrs = pi*n0s*avts*g3pbs*.25
    1588. 

  1588.    precs1 = 4.*n0s*.65
    1589. 

  1589.    precs2 = 4.*n0s*.44*avts**.5*g5pbso2
    1590. 

  1590.    pidn0s =  pi*dens*n0s
    1591. 

  1591.    pacrc = pi*n0s*avts*g3pbs*.25*eacrc
    1592. 

  1592. !
    1593. 

  1593.    rslopecmax = 1./lamdacmax
    1594. 

  1594.    rslopermax = 1./lamdarmax
    1595. 

  1595.    rslopesmax = 1./lamdasmax
    1596. 

  1596.    rsloperbmax = rslopermax ** bvtr
    1597. 

  1597.    rslopesbmax = rslopesmax ** bvts
    1598. 

  1598.    rslopec2max = rslopecmax * rslopecmax
    1599. 

  1599.    rsloper2max = rslopermax * rslopermax
    1600. 

  1600.    rslopes2max = rslopesmax * rslopesmax
    1601. 

  1601.    rslopec3max = rslopec2max * rslopecmax
    1602. 

  1602.    rsloper3max = rsloper2max * rslopermax
    1603. 

  1603.    rslopes3max = rslopes2max * rslopesmax
    1604. 

  1604. !
    1605. 

  1605.   END SUBROUTINE wdm5init
    1606. 

  1606. !------------------------------------------------------------------------------
    1607. 

  1607.       subroutine slope_wdm5(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
    1608. 

  1608.                             vt,vtn,its,ite,kts,kte)
    1609. 

  1609.   IMPLICIT NONE
    1610. 

  1610.   INTEGER       ::               its,ite, jts,jte, kts,kte
    1611. 

  1611.   REAL, DIMENSION( its:ite , kts:kte,2) ::                                     &
    1612. 

  1612.                                                                           qrs, &
    1613. 

  1613.                                                                        rslope, &
    1614. 

  1614.                                                                       rslopeb, &
    1615. 

  1615.                                                                       rslope2, &
    1616. 

  1616.                                                                       rslope3, &
    1617. 

  1617.                                                                            vt
    1618. 

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

  1619.                                                                           ncr, &
    1620. 

  1620.                                                                           vtn, &
    1621. 

  1621.                                                                           den, &
    1622. 

  1622.                                                                        denfac, &
    1623. 

  1623.                                                                             t
    1624. 

  1624.   REAL, PARAMETER  :: t0c = 273.15
    1625. 

  1625.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    1626. 

  1626.                                                                        n0sfac
    1627. 

  1627.   REAL       ::  lamdar, lamdas, x, y, z, supcol
    1628. 

  1628.   integer :: i, j, k
    1629. 

  1629. !----------------------------------------------------------------
    1630. 

  1630. !     size distributions: (x=mixing ratio, y=air density):
    1631. 

  1631. !     valid for mixing ratio > 1.e-9 kg/kg.
    1632. 

  1632. !
    1633. 

  1633. !  Optimizatin : A**B => exp(log(A)*(B))
    1634. 

  1634.       lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
    1635. 

  1635.       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
    1636. 

  1636. !
    1637. 

  1637.       do k = kts, kte
    1638. 

  1638.         do i = its, ite
    1639. 

  1639.           supcol = t0c-t(i,k)
    1640. 

  1640. !---------------------------------------------------------------
    1641. 

  1641. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
    1642. 

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

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

  1644.           if(qrs(i,k,1).le.qcrmin .or. ncr(i,k).le.nrmin ) then
    1645. 

  1645.             rslope(i,k,1) = rslopermax
    1646. 

  1646.             rslopeb(i,k,1) = rsloperbmax
    1647. 

  1647.             rslope2(i,k,1) = rsloper2max
    1648. 

  1648.             rslope3(i,k,1) = rsloper3max
    1649. 

  1649.           else
    1650. 

  1650.             rslope(i,k,1) = min(1./lamdar(qrs(i,k,1),den(i,k),ncr(i,k)),1.e-3)
    1651. 

  1651.             rslopeb(i,k,1) = rslope(i,k,1)**bvtr
    1652. 

  1652.             rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
    1653. 

  1653.             rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
    1654. 

  1654.           endif
    1655. 

  1655.           if(qrs(i,k,2).le.qcrmin) then
    1656. 

  1656.             rslope(i,k,2) = rslopesmax
    1657. 

  1657.             rslopeb(i,k,2) = rslopesbmax
    1658. 

  1658.             rslope2(i,k,2) = rslopes2max
    1659. 

  1659.             rslope3(i,k,2) = rslopes3max
    1660. 

  1660.           else
    1661. 

  1661.             rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
    1662. 

  1662.             rslopeb(i,k,2) = rslope(i,k,2)**bvts
    1663. 

  1663.             rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
    1664. 

  1664.             rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
    1665. 

  1665.           endif
    1666. 

  1666.           vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
    1667. 

  1667.           vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
    1668. 

  1668.           vtn(i,k) = pvtrn*rslopeb(i,k,1)*denfac(i,k)
    1669. 

  1669.           if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
    1670. 

  1670.           if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
    1671. 

  1671.           if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
    1672. 

  1672.         enddo
    1673. 

  1673.       enddo
    1674. 

  1674.   END subroutine slope_wdm5
    1675. 

  1675. !-----------------------------------------------------------------------------
    1676. 

  1676.       subroutine slope_rain(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
    1677. 

  1677.                             vt,vtn,its,ite,kts,kte)
    1678. 

  1678.   IMPLICIT NONE
    1679. 

  1679.   INTEGER       ::               its,ite, jts,jte, kts,kte
    1680. 

  1680.   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
    1681. 

  1681.                                                                           qrs, &
    1682. 

  1682.                                                                           ncr, &
    1683. 

  1683.                                                                        rslope, &
    1684. 

  1684.                                                                       rslopeb, &
    1685. 

  1685.                                                                       rslope2, &
    1686. 

  1686.                                                                       rslope3, &
    1687. 

  1687.                                                                            vt, &
    1688. 

  1688.                                                                           vtn, &
    1689. 

  1689.                                                                           den, &
    1690. 

  1690.                                                                        denfac, &
    1691. 

  1691.                                                                             t
    1692. 

  1692.   REAL, PARAMETER  :: t0c = 273.15
    1693. 

  1693.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    1694. 

  1694.                                                                        n0sfac
    1695. 

  1695.   REAL       ::  lamdar, x, y, z, supcol
    1696. 

  1696.   integer :: i, j, k
    1697. 

  1697. !----------------------------------------------------------------
    1698. 

  1698. !     size distributions: (x=mixing ratio, y=air density):
    1699. 

  1699. !     valid for mixing ratio > 1.e-9 kg/kg.
    1700. 

  1700.       lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
    1701. 

  1701. !
    1702. 

  1702.       do k = kts, kte
    1703. 

  1703.         do i = its, ite
    1704. 

  1704.           if(qrs(i,k).le.qcrmin .or. ncr(i,k).le.nrmin) then
    1705. 

  1705.             rslope(i,k) = rslopermax
    1706. 

  1706.             rslopeb(i,k) = rsloperbmax
    1707. 

  1707.             rslope2(i,k) = rsloper2max
    1708. 

  1708.             rslope3(i,k) = rsloper3max
    1709. 

  1709.           else
    1710. 

  1710.             rslope(i,k) = min(1./lamdar(qrs(i,k),den(i,k),ncr(i,k)),1.e-3)
    1711. 

  1711.             rslopeb(i,k) = rslope(i,k)**bvtr
    1712. 

  1712.             rslope2(i,k) = rslope(i,k)*rslope(i,k)
    1713. 

  1713.             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
    1714. 

  1714.           endif
    1715. 

  1715.           vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
    1716. 

  1716.           vtn(i,k) = pvtrn*rslopeb(i,k)*denfac(i,k)
    1717. 

  1717.           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
    1718. 

  1718.           if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
    1719. 

  1719.         enddo
    1720. 

  1720.       enddo
    1721. 

  1721.   END subroutine slope_rain
    1722. 

  1722. !------------------------------------------------------------------------------
    1723. 

  1723.       subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,   &
    1724. 

  1724.                             vt,its,ite,kts,kte)
    1725. 

  1725.   IMPLICIT NONE
    1726. 

  1726.   INTEGER       ::               its,ite, jts,jte, kts,kte
    1727. 

  1727.   REAL, DIMENSION( its:ite , kts:kte) ::                                       &
    1728. 

  1728.                                                                           qrs, &
    1729. 

  1729.                                                                        rslope, &
    1730. 

  1730.                                                                       rslopeb, &
    1731. 

  1731.                                                                       rslope2, &
    1732. 

  1732.                                                                       rslope3, &
    1733. 

  1733.                                                                            vt, &
    1734. 

  1734.                                                                           den, &
    1735. 

  1735.                                                                        denfac, &
    1736. 

  1736.                                                                             t
    1737. 

  1737.   REAL, PARAMETER  :: t0c = 273.15
    1738. 

  1738.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    1739. 

  1739.                                                                        n0sfac
    1740. 

  1740.   REAL       ::  lamdas, x, y, z, supcol
    1741. 

  1741.   integer :: i, j, k
    1742. 

  1742. !----------------------------------------------------------------
    1743. 

  1743. !     size distributions: (x=mixing ratio, y=air density):
    1744. 

  1744. !     valid for mixing ratio > 1.e-9 kg/kg.
    1745. 

  1745.       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
    1746. 

  1746. !
    1747. 

  1747.       do k = kts, kte
    1748. 

  1748.         do i = its, ite
    1749. 

  1749.           supcol = t0c-t(i,k)
    1750. 

  1750. !---------------------------------------------------------------
    1751. 

  1751. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
    1752. 

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

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

  1754.           if(qrs(i,k).le.qcrmin)then
    1755. 

  1755.             rslope(i,k) = rslopesmax
    1756. 

  1756.             rslopeb(i,k) = rslopesbmax
    1757. 

  1757.             rslope2(i,k) = rslopes2max
    1758. 

  1758.             rslope3(i,k) = rslopes3max
    1759. 

  1759.           else
    1760. 

  1760.             rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
    1761. 

  1761.             rslopeb(i,k) = rslope(i,k)**bvts
    1762. 

  1762.             rslope2(i,k) = rslope(i,k)*rslope(i,k)
    1763. 

  1763.             rslope3(i,k) = rslope2(i,k)*rslope(i,k)
    1764. 

  1764.           endif
    1765. 

  1765.           vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
    1766. 

  1766.           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
    1767. 

  1767.         enddo
    1768. 

  1768.       enddo
    1769. 

  1769.   END subroutine slope_snow
    1770. 

  1770. !-------------------------------------------------------------------
    1771. 

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

  1772. !-------------------------------------------------------------------
    1773. 

  1773. !
    1774. 

  1774. ! for non-iteration semi-Lagrangain forward advection for cloud
    1775. 

  1775. ! with mass conservation and positive definite advection
    1776. 

  1776. ! 2nd order interpolation with monotonic piecewise linear method
    1777. 

  1777. ! this routine is under assumption of decfl < 1 for semi_Lagrangian
    1778. 

  1778. !
    1779. 

  1779. ! dzl    depth of model layer in meter
    1780. 

  1780. ! wwl    terminal velocity at model layer m/s
    1781. 

  1781. ! rql    cloud density*mixing ration
    1782. 

  1782. ! precip precipitation
    1783. 

  1783. ! dt     time step
    1784. 

  1784. ! id     kind of precip: 0 test case; 1 raindrop  2: snow
    1785. 

  1785. ! iter   how many time to guess mean terminal velocity: 0 pure forward.
    1786. 

  1786. !        0 : use departure wind for advection
    1787. 

  1787. !        1 : use mean wind for advection
    1788. 

  1788. !        > 1 : use mean wind after iter-1 iterations
    1789. 

  1789. !
    1790. 

  1790. ! author: hann-ming henry juang <henry.juang@noaa.gov>
    1791. 

  1791. !         implemented by song-you hong
    1792. 

  1792. !
    1793. 

  1793.       implicit none
    1794. 

  1794.       integer  im,km,id
    1795. 

  1795.       real  dt
    1796. 

  1796.       real  dzl(im,km),wwl(im,km),rql(im,km),precip(im)
    1797. 

  1797.       real  denl(im,km),denfacl(im,km),tkl(im,km)
    1798. 

  1798. !
    1799. 

  1799.       integer  i,k,n,m,kk,kb,kt,iter
    1800. 

  1800.       real  tl,tl2,qql,dql,qqd
    1801. 

  1801.       real  th,th2,qqh,dqh
    1802. 

  1802.       real  zsum,qsum,dim,dip,c1,con1,fa1,fa2
    1803. 

  1803.       real  allold, allnew, zz, dzamin, cflmax, decfl
    1804. 

  1804.       real  dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
    1805. 

  1805.       real  den(km), denfac(km), tk(km)
    1806. 

  1806.       real  wi(km+1), zi(km+1), za(km+1)
    1807. 

  1807.       real  qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
    1808. 

  1808.       real  dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
    1809. 

  1809. !
    1810. 

  1810.       precip(:) = 0.0
    1811. 

  1811. !
    1812. 

  1812.       i_loop : do i=1,im
    1813. 

  1813. ! -----------------------------------
    1814. 

  1814.       dz(:) = dzl(i,:)
    1815. 

  1815.       qq(:) = rql(i,:)
    1816. 

  1816.       ww(:) = wwl(i,:)
    1817. 

  1817.       den(:) = denl(i,:)
    1818. 

  1818.       denfac(:) = denfacl(i,:)
    1819. 

  1819.       tk(:) = tkl(i,:)
    1820. 

  1820. ! skip for no precipitation for all layers
    1821. 

  1821.       allold = 0.0
    1822. 

  1822.       do k=1,km
    1823. 

  1823.         allold = allold + qq(k)
    1824. 

  1824.       enddo
    1825. 

  1825.       if(allold.le.0.0) then
    1826. 

  1826.         cycle i_loop
    1827. 

  1827.       endif
    1828. 

  1828. !
    1829. 

  1829. ! compute interface values
    1830. 

  1830.       zi(1)=0.0
    1831. 

  1831.       do k=1,km
    1832. 

  1832.         zi(k+1) = zi(k)+dz(k)
    1833. 

  1833.       enddo
    1834. 

  1834. !
    1835. 

  1835. ! save departure wind
    1836. 

  1836.       wd(:) = ww(:)
    1837. 

  1837.       n=1
    1838. 

  1838.  100  continue
    1839. 

  1839. ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
    1840. 

  1840. ! 2nd order interpolation to get wi
    1841. 

  1841.       wi(1) = ww(1)
    1842. 

  1842.       wi(km+1) = ww(km)
    1843. 

  1843.       do k=2,km
    1844. 

  1844.         wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
    1845. 

  1845.       enddo
    1846. 

  1846. ! 3rd order interpolation to get wi
    1847. 

  1847.       fa1 = 9./16.
    1848. 

  1848.       fa2 = 1./16.
    1849. 

  1849.       wi(1) = ww(1)
    1850. 

  1850.       wi(2) = 0.5*(ww(2)+ww(1))
    1851. 

  1851.       do k=3,km-1
    1852. 

  1852.         wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
    1853. 

  1853.       enddo
    1854. 

  1854.       wi(km) = 0.5*(ww(km)+ww(km-1))
    1855. 

  1855.       wi(km+1) = ww(km)
    1856. 

  1856. !
    1857. 

  1857. ! terminate of top of raingroup
    1858. 

  1858.       do k=2,km
    1859. 

  1859.         if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
    1860. 

  1860.       enddo
    1861. 

  1861. !
    1862. 

  1862. ! diffusivity of wi
    1863. 

  1863.       con1 = 0.05
    1864. 

  1864.       do k=km,1,-1
    1865. 

  1865.         decfl = (wi(k+1)-wi(k))*dt/dz(k)
    1866. 

  1866.         if( decfl .gt. con1 ) then
    1867. 

  1867.           wi(k) = wi(k+1) - con1*dz(k)/dt
    1868. 

  1868.         endif
    1869. 

  1869.       enddo
    1870. 

  1870. ! compute arrival point
    1871. 

  1871.       do k=1,km+1
    1872. 

  1872.         za(k) = zi(k) - wi(k)*dt
    1873. 

  1873.       enddo
    1874. 

  1874. !
    1875. 

  1875.       do k=1,km
    1876. 

  1876.         dza(k) = za(k+1)-za(k)
    1877. 

  1877.       enddo
    1878. 

  1878.       dza(km+1) = zi(km+1) - za(km+1)
    1879. 

  1879. !
    1880. 

  1880. ! computer deformation at arrival point
    1881. 

  1881.       do k=1,km
    1882. 

  1882.         qa(k) = qq(k)*dz(k)/dza(k)
    1883. 

  1883.         qr(k) = qa(k)/den(k)
    1884. 

  1884.       enddo
    1885. 

  1885.       qa(km+1) = 0.0
    1886. 

  1886. !     call maxmin(km,1,qa,' arrival points ')
    1887. 

  1887. !
    1888. 

  1888. ! compute arrival terminal velocity, and estimate mean terminal velocity
    1889. 

  1889. ! then back to use mean terminal velocity
    1890. 

  1890.       if( n.le.iter ) then
    1891. 

  1891. !        if (id.eq.1) then
    1892. 

  1892. !
    1893. 

  1893. !          call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
    1894. 

  1894. !        else
    1895. 

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

  1896. !        endif
    1897. 

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

  1898.         do k=1,km
    1899. 

  1899. !#ifdef DEBUG
    1900. 

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

  1901. !#endif
    1902. 

  1902. ! mean wind is average of departure and new arrival winds
    1903. 

  1903.           ww(k) = 0.5* ( wd(k)+wa(k) )
    1904. 

  1904.         enddo
    1905. 

  1905.         was(:) = wa(:)
    1906. 

  1906.         n=n+1
    1907. 

  1907.         go to 100
    1908. 

  1908.       endif
    1909. 

  1909. !
    1910. 

  1910. ! estimate values at arrival cell interface with monotone
    1911. 

  1911.       do k=2,km
    1912. 

  1912.         dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
    1913. 

  1913.         dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
    1914. 

  1914.         if( dip*dim.le.0.0 ) then
    1915. 

  1915.           qmi(k)=qa(k)
    1916. 

  1916.           qpi(k)=qa(k)
    1917. 

  1917.         else
    1918. 

  1918.           qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
    1919. 

  1919.           qmi(k)=2.0*qa(k)-qpi(k)
    1920. 

  1920.           if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
    1921. 

  1921.             qpi(k) = qa(k)
    1922. 

  1922.             qmi(k) = qa(k)
    1923. 

  1923.           endif
    1924. 

  1924.         endif
    1925. 

  1925.       enddo
    1926. 

  1926.       qpi(1)=qa(1)
    1927. 

  1927.       qmi(1)=qa(1)
    1928. 

  1928.       qmi(km+1)=qa(km+1)
    1929. 

  1929.       qpi(km+1)=qa(km+1)
    1930. 

  1930. !
    1931. 

  1931. ! interpolation to regular point
    1932. 

  1932.       qn = 0.0
    1933. 

  1933.       kb=1
    1934. 

  1934.       kt=1
    1935. 

  1935.       intp : do k=1,km
    1936. 

  1936.              kb=max(kb-1,1)
    1937. 

  1937.              kt=max(kt-1,1)
    1938. 

  1938. ! find kb and kt
    1939. 

  1939.              if( zi(k).ge.za(km+1) ) then
    1940. 

  1940.                exit intp
    1941. 

  1941.              else
    1942. 

  1942.                find_kb : do kk=kb,km
    1943. 

  1943.                          if( zi(k).le.za(kk+1) ) then
    1944. 

  1944.                            kb = kk
    1945. 

  1945.                            exit find_kb
    1946. 

  1946.                          else
    1947. 

  1947.                            cycle find_kb
    1948. 

  1948.                          endif
    1949. 

  1949.                enddo find_kb
    1950. 

  1950.                find_kt : do kk=kt,km
    1951. 

  1951.                          if( zi(k+1).le.za(kk) ) then
    1952. 

  1952.                            kt = kk
    1953. 

  1953.                            exit find_kt
    1954. 

  1954.                          else
    1955. 

  1955.                            cycle find_kt
    1956. 

  1956.                          endif
    1957. 

  1957.                enddo find_kt
    1958. 

  1958.                kt = kt - 1
    1959. 

  1959. ! compute q with piecewise constant method
    1960. 

  1960.                if( kt.eq.kb ) then
    1961. 

  1961.                  tl=(zi(k)-za(kb))/dza(kb)
    1962. 

  1962.                  th=(zi(k+1)-za(kb))/dza(kb)
    1963. 

  1963.                  tl2=tl*tl
    1964. 

  1964.                  th2=th*th
    1965. 

  1965.                  qqd=0.5*(qpi(kb)-qmi(kb))
    1966. 

  1966.                  qqh=qqd*th2+qmi(kb)*th
    1967. 

  1967.                  qql=qqd*tl2+qmi(kb)*tl
    1968. 

  1968.                  qn(k) = (qqh-qql)/(th-tl)
    1969. 

  1969.                else if( kt.gt.kb ) then
    1970. 

  1970.                  tl=(zi(k)-za(kb))/dza(kb)
    1971. 

  1971.                  tl2=tl*tl
    1972. 

  1972.                  qqd=0.5*(qpi(kb)-qmi(kb))
    1973. 

  1973.                  qql=qqd*tl2+qmi(kb)*tl
    1974. 

  1974.                  dql = qa(kb)-qql
    1975. 

  1975.                  zsum  = (1.-tl)*dza(kb)
    1976. 

  1976.                  qsum  = dql*dza(kb)
    1977. 

  1977.                  if( kt-kb.gt.1 ) then
    1978. 

  1978.                  do m=kb+1,kt-1
    1979. 

  1979.                    zsum = zsum + dza(m)
    1980. 

  1980.                    qsum = qsum + qa(m) * dza(m)
    1981. 

  1981.                  enddo
    1982. 

  1982.                  endif
    1983. 

  1983.                  th=(zi(k+1)-za(kt))/dza(kt)
    1984. 

  1984.                  th2=th*th
    1985. 

  1985.                  qqd=0.5*(qpi(kt)-qmi(kt))
    1986. 

  1986.                  dqh=qqd*th2+qmi(kt)*th
    1987. 

  1987.                  zsum  = zsum + th*dza(kt)
    1988. 

  1988.                  qsum  = qsum + dqh*dza(kt)
    1989. 

  1989.                  qn(k) = qsum/zsum
    1990. 

  1990.                endif
    1991. 

  1991.                cycle intp
    1992. 

  1992.              endif
    1993. 

  1993. !
    1994. 

  1994.        enddo intp
    1995. 

  1995. !
    1996. 

  1996. ! rain out
    1997. 

  1997.       sum_precip: do k=1,km
    1998. 

  1998.                     if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
    1999. 

  1999.                       precip(i) = precip(i) + qa(k)*dza(k)
    2000. 

  2000.                       cycle sum_precip
    2001. 

  2001.                     else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
    2002. 

  2002.                       precip(i) = precip(i) + qa(k)*(0.0-za(k))
    2003. 

  2003.                       exit sum_precip
    2004. 

  2004.                     endif
    2005. 

  2005.                     exit sum_precip
    2006. 

  2006.       enddo sum_precip
    2007. 

  2007. !
    2008. 

  2008. ! replace the new values
    2009. 

  2009.       rql(i,:) = qn(:)
    2010. 

  2010. !
    2011. 

  2011. ! ----------------------------------
    2012. 

  2012.       enddo i_loop
    2013. 

  2013. !
    2014. 

  2014.   END SUBROUTINE nislfv_rain_plm
    2015. 

  2015. END MODULE module_mp_wdm5
    2016. 

  2016.