PageRenderTime 45ms CodeModel.GetById 34ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_mp_wsm5_accel.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 3093 lines | 2052 code | 59 blank | 982 comment | 153 complexity | 4a961fc46c3dd28793c37a174419a9c5 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. !Including inline expansion statistical function 
    10. 

  10. MODULE module_mp_wsm5
    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 :: lamdarmax = 8.e4   ! limited maximum value for slope parameter of rain
    25. 

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

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

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

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

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

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

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

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

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

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

  35.    REAL, SAVE ::                                      &
    36. 

  36.              qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
    37. 

  37.              g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr,    &
    38. 

  38.              precr1,precr2,xmmax,roqimax,bvts1,       &
    39. 

  39.              bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs,     &
    40. 

  40.              g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
    41. 

  41.              pidn0s,xlv1,pacrc,                       &
    42. 

  42.              rslopermax,rslopesmax,rslopegmax,        &
    43. 

  43.              rsloperbmax,rslopesbmax,rslopegbmax,     &
    44. 

  44.              rsloper2max,rslopes2max,rslopeg2max,     &
    45. 

  45.              rsloper3max,rslopes3max,rslopeg3max
    46. 

  46. !
    47. 

  47. ! Specifies code-inlining of fpvs function in WSM52D below. JM 20040507
    48. 

  48. !
    49. 

  49. CONTAINS
    50. 

  50. !===================================================================
    51. 

  51. !
    52. 

  52.   SUBROUTINE wsm5(th, q, qc, qr, qi, qs                            &
    53. 

  53.                  ,den, pii, p, delz                                &
    54. 

  54.                  ,delt,g, cpd, cpv, rd, rv, t0c                    &
    55. 

  55.                  ,ep1, ep2, qmin                                   &
    56. 

  56.                  ,XLS, XLV0, XLF0, den0, denr                      &
    57. 

  57.                  ,cliq,cice,psat                                   &
    58. 

  58.                  ,rain, rainncv                                    &
    59. 

  59.                  ,snow, snowncv                                    &
    60. 

  60.                  ,sr                                               &
    61. 

  61.                  ,ids,ide, jds,jde, kds,kde                        &
    62. 

  62.                  ,ims,ime, jms,jme, kms,kme                        &
    63. 

  63.                  ,its,ite, jts,jte, kts,kte                        &
    64. 

  64.                                                                    )
    65. 

  65. #ifdef _OPENMP
    66. 

  66.   use omp_lib
    67. 

  67. #endif
    68. 

  68. !-------------------------------------------------------------------
    69. 

  69.   IMPLICIT NONE
    70. 

  70. !-------------------------------------------------------------------
    71. 

  71. !
    72. 

  72. !  This code is a 5-class mixed ice microphyiscs scheme (WSM5) of the WRF
    73. 

  73. !  Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
    74. 

  74. !  number concentration is a function of temperature, and seperate assumption
    75. 

  75. !  is developed, in which ice crystal number concentration is a function
    76. 

  76. !  of ice amount. A theoretical background of the ice-microphysics and related
    77. 

  77. !  processes in the WSMMPs are described in Hong et al. (2004).
    78. 

  78. !  Production terms in the WSM6 scheme are described in Hong and Lim (2006).
    79. 

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

  80. !
    81. 

  81. !  WSM5 cloud scheme
    82. 

  82. !
    83. 

  83. !  Coded by Song-You Hong (Yonsei Univ.)
    84. 

  84. !             Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
    85. 

  85. !             Summer 2002
    86. 

  86. !
    87. 

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

  88. !             Summer 2003
    89. 

  89. !
    90. 

  90. !  Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
    91. 

  91. !             Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
    92. 

  92. !             Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
    93. 

  93. !
    94. 

  94.   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
    95. 

  95.                                       ims,ime, jms,jme, kms,kme , &
    96. 

  96.                                       its,ite, jts,jte, kts,kte
    97. 

  97.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    98. 

  98.         INTENT(INOUT) ::                                          &
    99. 

  99.                                                              th,  &
    100. 

  100.                                                               q,  &
    101. 

  101.                                                               qc, &
    102. 

  102.                                                               qi, &
    103. 

  103.                                                               qr, &
    104. 

  104.                                                               qs
    105. 

  105.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    106. 

  106.         INTENT(IN   ) ::                                          &
    107. 

  107.                                                              den, &
    108. 

  108.                                                              pii, &
    109. 

  109.                                                                p, &
    110. 

  110.                                                             delz
    111. 

  111.   REAL, INTENT(IN   ) ::                                    delt, &
    112. 

  112.                                                                g, &
    113. 

  113.                                                               rd, &
    114. 

  114.                                                               rv, &
    115. 

  115.                                                              t0c, &
    116. 

  116.                                                             den0, &
    117. 

  117.                                                              cpd, &
    118. 

  118.                                                              cpv, &
    119. 

  119.                                                              ep1, &
    120. 

  120.                                                              ep2, &
    121. 

  121.                                                             qmin, &
    122. 

  122.                                                              XLS, &
    123. 

  123.                                                             XLV0, &
    124. 

  124.                                                             XLF0, &
    125. 

  125.                                                             cliq, &
    126. 

  126.                                                             cice, &
    127. 

  127.                                                             psat, &
    128. 

  128.                                                             denr
    129. 

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

  130.         INTENT(INOUT) ::                                    rain, &
    131. 

  131.                                                          rainncv, &
    132. 

  132.                                                               sr
    133. 

  133. 

  134.   REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                 &
    135. 

  135.         INTENT(INOUT) ::                                    snow, &
    136. 

  136.                                                          snowncv
    137. 

  137. 

  138. ! LOCAL VAR
    139. 

  139.   REAL, DIMENSION( its:ite , kts:kte ) ::   t
    140. 

  140.   REAL, DIMENSION( its:ite , kts:kte, 2 ) ::   qci, qrs
    141. 

  141.   CHARACTER*256 :: emess
    142. 

  142.   INTEGER :: mkx_test
    143. 

  143.   INTEGER ::               i,j,k
    144. 

  144. 

  145. #ifdef _ACCEL_PROF
    146. 

  146.   INTEGER :: l
    147. 

  147.   real*8 wsm3_t(8,256), wsm5_t(8,256), t1, t2
    148. 

  148.   common /wsm_times/ wsm3_t(8,256), wsm5_t(8,256)
    149. 

  149. #endif
    150. 

  150. 

  151. 

  152. !-------------------------------------------------------------------
    153. 

  153. 

  154. #ifdef _ACCEL_PROF
    155. 

  155. call cpu_time(t1)
    156. 

  156. #endif
    157. 

  157. 

  158. #ifndef RUN_ON_GPU
    159. 

  159. 

  160. #ifdef _ACCEL
    161. 

  161. 

  162.       !  Need to send th, pii, qc, qi, qr, qs
    163. 

  163.       !  Don't send t
    164. 

    165. 

  165.       CALL wsm52D(th, pii, q, qc, qr, qi, qs                       &
    166. 

  166.                   ,den                                             &
    167. 

  167.                   ,p, delz                                         &
    168. 

  168.                   ,delt,g, cpd, cpv, rd, rv, t0c                   &
    169. 

  169.                   ,ep1, ep2, qmin                                  &
    170. 

  170.                   ,XLS, XLV0, XLF0, den0, denr                     &
    171. 

  171.                   ,cliq,cice,psat                                  &
    172. 

  172.                   ,rain,rainncv                                    &
    173. 

  173.                   ,sr                                              &
    174. 

  174.                   ,ids,ide, jds,jde, kds,kde                       &
    175. 

  175.                   ,ims,ime, jms,jme, kms,kme                       &
    176. 

  176.                   ,its,ite, jts,jte, kts,kte                       &
    177. 

  177.                   ,snow,snowncv                                    &
    178. 

  178.                                                                    )
    179. 

  179. 

  180. #else
    181. 

  181. 

  182.       DO j=jts,jte
    183. 

  183.          DO k=kts,kte
    184. 

  184.          DO i=its,ite
    185. 

  185.             t(i,k)=th(i,k,j)*pii(i,k,j)
    186. 

  186.             qci(i,k,1) = qc(i,k,j)
    187. 

  187.             qci(i,k,2) = qi(i,k,j)
    188. 

  188.             qrs(i,k,1) = qr(i,k,j)
    189. 

  189.             qrs(i,k,2) = qs(i,k,j)
    190. 

  190.          ENDDO
    191. 

  191.          ENDDO
    192. 

  192. 

  193.          !  Sending array starting locations of optional variables may cause
    194. 

  194.          !  troubles, so we explicitly change the call.
    195. 

  195. 

  196.          CALL wsm52D(t, q(ims,kms,j), qci, qrs                     &
    197. 

  197.                     ,den(ims,kms,j)                                &
    198. 

  198.                     ,p(ims,kms,j), delz(ims,kms,j)                 &
    199. 

  199.                     ,delt,g, cpd, cpv, rd, rv, t0c                 &
    200. 

  200.                     ,ep1, ep2, qmin                                &
    201. 

  201.                     ,XLS, XLV0, XLF0, den0, denr                   &
    202. 

  202.                     ,cliq,cice,psat                                &
    203. 

  203.                     ,j                                             &
    204. 

  204.                     ,rain(ims,j),rainncv(ims,j)                    &
    205. 

  205.                     ,sr(ims,j)                                     &
    206. 

  206.                     ,ids,ide, jds,jde, kds,kde                     &
    207. 

  207.                     ,ims,ime, jms,jme, kms,kme                     &
    208. 

  208.                     ,its,ite, jts,jte, kts,kte                     &
    209. 

  209.                     ,snow,snowncv                                  &
    210. 

  210.                                                                    )
    211. 

  211. 

  212.          DO K=kts,kte
    213. 

  213.          DO I=its,ite
    214. 

  214.             th(i,k,j)=t(i,k)/pii(i,k,j)
    215. 

  215.             qc(i,k,j) = qci(i,k,1)
    216. 

  216.             qi(i,k,j) = qci(i,k,2)
    217. 

  217.             qr(i,k,j) = qrs(i,k,1)
    218. 

  218.             qs(i,k,j) = qrs(i,k,2)
    219. 

  219.          ENDDO
    220. 

  220.          ENDDO
    221. 

  221.       ENDDO
    222. 

  222. #endif
    223. 

  223. #else
    224. 

  224.       CALL get_wsm5_gpu_levels ( mkx_test )
    225. 

  225.       IF ( mkx_test .LT. kte ) THEN
    226. 

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

  227.                       mkx_test,' < ',kte
    228. 

  228.         CALL wrf_error_fatal(emess)
    229. 

  229.       ENDIF
    230. 

  230.       CALL wsm5_host (                                                           &
    231. 

  231.                     th(its:ite,kts:kte,jts:jte), pii(its:ite,kts:kte,jts:jte)    &
    232. 

  232.                    ,q(its:ite,kts:kte,jts:jte), qc(its:ite,kts:kte,jts:jte)      &
    233. 

  233.                    ,qi(its:ite,kts:kte,jts:jte), qr(its:ite,kts:kte,jts:jte)     &
    234. 

  234.                    ,qs(its:ite,kts:kte,jts:jte), den(its:ite,kts:kte,jts:jte)    &
    235. 

  235.                    ,p(its:ite,kts:kte,jts:jte), delz(its:ite,kts:kte,jts:jte)    &
    236. 

  236.                    ,delt                                                         &
    237. 

  237.                    ,rain(its:ite,jts:jte),rainncv(its:ite,jts:jte)               &
    238. 

  238.                    ,snow(its:ite,jts:jte),snowncv(its:ite,jts:jte)               &
    239. 

  239.                    ,sr(its:ite,jts:jte)                                          &
    240. 

  240.                    ,its, ite,  jts, jte,  kts, kte                               &
    241. 

  241.                    ,its, ite,  jts, jte,  kts, kte                               &
    242. 

  242.                    ,its, ite,  jts, jte,  kts, kte                               &
    243. 

  243.           )
    244. 

  244. #endif
    245. 

  245. 

  246. 

  247. #ifdef _ACCEL_PROF
    248. 

  248.   call cpu_time(t2)
    249. 

  249. #ifdef _OPENMP
    250. 

  250.   l = omp_get_thread_num() + 1
    251. 

  251. #else
    252. 

  252.   l = 1
    253. 

  253. #endif
    254. 

  254.   wsm5_t(1,l) = wsm5_t(1,l) + (t2 - t1)
    255. 

  255. #endif
    256. 

  256. 

  257.   END SUBROUTINE wsm5
    258. 

  258. 

  259. 

  260. #ifdef _ACCEL
    261. 

  261. 

  262. !===================================================================
    263. 

  263. !
    264. 

  264.   SUBROUTINE wsm52D(th, pii, q, qc, qr, qi, qqs, den, p, delz     &
    265. 

  265.                    ,delt,g, cpd, cpv, rd, rv, t0c                 &
    266. 

  266.                    ,ep1, ep2, qmin                                &
    267. 

  267.                    ,XLS, XLV0, XLF0, den0, denr                   &
    268. 

  268.                    ,cliq,cice,psat                                &
    269. 

  269.                    ,rain,rainncv                                  &
    270. 

  270.                    ,sr                                            &
    271. 

  271.                    ,ids,ide, jds,jde, kds,kde                     &
    272. 

  272.                    ,ims,ime, jms,jme, kms,kme                     &
    273. 

  273.                    ,its,ite, jts,jte, kts,kte                     &
    274. 

  274.                    ,snow,snowncv                                  &
    275. 

  275.                                                                   )
    276. 

  276. !-------------------------------------------------------------------
    277. 

  277.   IMPLICIT NONE
    278. 

  278. !-------------------------------------------------------------------
    279. 

  279.   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
    280. 

  280.                                       ims,ime, jms,jme, kms,kme , &
    281. 

  281.                                       its,ite, jts,jte, kts,kte
    282. 

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

  283.         INTENT(INOUT) ::                                          &
    284. 

  284.                                                               th
    285. 

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

  286.         INTENT(IN) ::                                             &
    287. 

  287.                                                              pii
    288. 

  288.   REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),                  &
    289. 

  289.         INTENT(INOUT) ::                                          &
    290. 

  290.                                                               qc, &
    291. 

  291.                                                               qr, &
    292. 

  292.                                                               qi, &
    293. 

  293.                                                              qqs
    294. 

  294.   REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),                  &
    295. 

  295.         INTENT(INOUT) ::                                          &
    296. 

  296.                                                                q
    297. 

  297.   REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),                  &
    298. 

  298.         INTENT(IN   ) ::                                          &
    299. 

  299.                                                              den, &
    300. 

  300.                                                                p, &
    301. 

  301.                                                             delz
    302. 

  302.   REAL, INTENT(IN   ) ::                                    delt, &
    303. 

  303.                                                                g, &
    304. 

  304.                                                              cpd, &
    305. 

  305.                                                              cpv, &
    306. 

  306.                                                              t0c, &
    307. 

  307.                                                             den0, &
    308. 

  308.                                                               rd, &
    309. 

  309.                                                               rv, &
    310. 

  310.                                                              ep1, &
    311. 

  311.                                                              ep2, &
    312. 

  312.                                                             qmin, &
    313. 

  313.                                                              XLS, &
    314. 

  314.                                                             XLV0, &
    315. 

  315.                                                             XLF0, &
    316. 

  316.                                                             cliq, &
    317. 

  317.                                                             cice, &
    318. 

  318.                                                             psat, &
    319. 

  319.                                                             denr
    320. 

  320.   REAL, DIMENSION( ims:ime, jms:jme ),                            &
    321. 

  321.         INTENT(INOUT) ::                                    rain, &
    322. 

  322.                                                          rainncv, &
    323. 

  323.                                                               sr
    324. 

  324. 

  325.   REAL, DIMENSION( ims:ime, jms:jme ),     OPTIONAL,              &
    326. 

  326.         INTENT(INOUT) ::                                    snow, &
    327. 

  327.                                                          snowncv
    328. 

  328. 

  329. ! LOCAL VAR
    330. 

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

  331.                                                               rh, &
    332. 

  332.                                                               qs, &
    333. 

  333.                                                           rslope, &
    334. 

  334.                                                          rslope2, &
    335. 

  335.                                                          rslope3, &
    336. 

  336.                                                          rslopeb, &
    337. 

  337.                                                             falk, &
    338. 

  338.                                                             fall, &
    339. 

  339.                                                            work1
    340. 

  340.   REAL, DIMENSION( its:ite , kts:kte, jts:jte ) ::                &
    341. 

  341.                                                                t
    342. 

  342.   REAL, DIMENSION( its:ite , kts:kte , 2 ) ::                     &
    343. 

  343.                                                              qci, &
    344. 

  344.                                                              qrs
    345. 

  345.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    346. 

  346.                                                            falkc, &
    347. 

  347.                                                            fallc, &
    348. 

  348.                                                               xl, &
    349. 

  349.                                                              cpm, &
    350. 

  350.                                                           denfac, &
    351. 

  351.                                                              xni, &
    352. 

  352.                                                           n0sfac, &
    353. 

  353.                                                            work2, &
    354. 

  354.                                                           work1c, &
    355. 

  355.                                                           work2c
    356. 

  356. 

  357.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    358. 

  358.                                                            pigen, &
    359. 

  359.                                                            pidep, &
    360. 

  360.                                                            psdep, &
    361. 

  361.                                                            praut, &
    362. 

  362.                                                            psaut, &
    363. 

  363.                                                            prevp, &
    364. 

  364.                                                            psevp, &
    365. 

  365.                                                            pracw, &
    366. 

  366.                                                            psacw, &
    367. 

  367.                                                            psaci, &
    368. 

  368.                                                            pcond, &
    369. 

  369.                                                            psmlt
    370. 

  370.   INTEGER ::                                                      &
    371. 

  371.                                                            mstep, &
    372. 

  372.                                                            numdt
    373. 

  373.   REAL ::                                                 rmstep
    374. 

  374.   REAL dtcldden, rdelz, rdtcld
    375. 

  375.   LOGICAL ::                                              flgcld
    376. 

  376. 

  377. #define WSM_NO_CONDITIONAL_IN_VECTOR
    378. 

  378. #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
    379. 

  379.   REAL ::                     xal, xbl
    380. 

  380. #endif
    381. 

  381. 

  382.   REAL  ::  pi,                                                   &
    383. 

  383.             cpmcal, xlcal, lamdar, lamdas, diffus,                &
    384. 

  384.             viscos, xka, venfac, conden, diffac,                  &
    385. 

  385.             x, y, z, a, b, c, d, e,                               &
    386. 

  386.             qdt, holdrr, holdrs, supcol, supcolt, pvt,            &
    387. 

  387.             coeres, supsat, dtcld, xmi, eacrs, satdt,             &
    388. 

  388.             vt2i,vt2s,acrfac,                                     &
    389. 

  389.             qimax, diameter, xni0, roqi0,                         &
    390. 

  390.             fallsum, fallsum_qsi, xlwork2, factor, source,        &
    391. 

  391.             value, xlf, pfrzdtc, pfrzdtr, supice,  holdc, holdci
    392. 

  392. ! variables for optimization
    393. 

  393.   REAL, DIMENSION( its:ite )           ::                    tvec1
    394. 

  394.   REAL ::                                                    temp 
    395. 

  395.   INTEGER :: i, j, k,                                             &
    396. 

  396.             iprt, latd, lond, loop, loops, ifsat, n
    397. 

  397. ! Temporaries used for inlining fpvs function
    398. 

  398.   REAL  :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
    399. 

  399.   REAL  :: logtr
    400. 

  400. !
    401. 

  401. !=================================================================
    402. 

  402. !   compute internal functions
    403. 

  403. !
    404. 

  404.       cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
    405. 

  405.       xlcal(x) = xlv0-xlv1*(x-t0c)
    406. 

  406. !----------------------------------------------------------------
    407. 

  407. !     size distributions: (x=mixing ratio, y=air density):
    408. 

  408. !     valid for mixing ratio > 1.e-9 kg/kg.
    409. 

  409. !
    410. 

  410. ! Optimizatin : A**B => exp(log(A)*(B))
    411. 

  411.       lamdar(x,y)=   sqrt(sqrt(pidn0r/(x*y)))      ! (pidn0r/(x*y))**.25
    412. 

  412.       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
    413. 

  413. !
    414. 

  414. !----------------------------------------------------------------
    415. 

  415. !     diffus: diffusion coefficient of the water vapor
    416. 

  416. !     viscos: kinematic viscosity(m2s-1)
    417. 

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

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

  419. !     xka(x,y) = 1.414e3*viscos(x,y)*y
    420. 

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

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

  422. !                    /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
    423. 

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

  424. !
    425. 

  425. !
    426. 

  426.       pi = 4. * atan(1.)
    427. 

  427. !
    428. 

  428. !----------------------------------------------------------------
    429. 

  429. !     paddint 0 for negative values generated by dynamics
    430. 

  430. !
    431. 

  431. 

  432. !
    433. 

  433. ! Moved outside of accelerator region
    434. 

  434. !
    435. 

  435.       loops = max(nint(delt/dtcldcr),1)
    436. 

  436.       dtcld = delt/loops
    437. 

  437.       if(delt.le.dtcldcr) dtcld = delt
    438. 

  438. !
    439. 

  439. 

  440. !....!$acc        local(t) &
    441. 

  441. 

  442.       IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
    443. 

  443. 

  444. !$acc region &
    445. 

  445. !$acc        local(t) &
    446. 

  446. !$acc        copyin(delz(:,:,:),p(:,:,:),den(:,:,:),pii(:,:,:)) &
    447. 

  447. !$acc        copyout(snowncv(:,:),rainncv(:,:),sr(:,:)) &
    448. 

  448. !$acc        copy(qqs(:,:,:),qr(:,:,:),qi(:,:,:),qc(:,:,:)) &
    449. 

  449. !$acc        copy(th(:,:,:),q(:,:,:),snow(:,:),rain(:,:))
    450. 

  450. !$acc do &
    451. 

  451. !$acc        private(rh,qs,rslope,rslope2,rslope3,rslopeb,falk,fall) &
    452. 

  452. !$acc        private(work1,qci,qrs,falkc,fallc,xl,cpm,denfac,xni) &
    453. 

  453. !$acc        private(n0sfac,work2,work1c,work2c,pigen,pidep,psdep) &
    454. 

  454. !$acc        private(praut,psaut,prevp,psevp) &
    455. 

  455. !$acc        private(pracw,psacw,psaci,pcond,psmlt) &
    456. 

  456. !$acc        parallel
    457. 

  457.       do j = jts, jte
    458. 

  458. !$acc do &
    459. 

  459. !$acc        private(numdt,mstep) &
    460. 

  460. !$acc        kernel vector
    461. 

  461.         do i = its, ite
    462. 

  462.         do k = kts, kte
    463. 

  463.           t(i,k,j)=th(i,k,j)*pii(i,k,j)
    464. 

  464.           qci(i,k,1) = max(qc(i,k,j),0.0)
    465. 

  465.           qci(i,k,2) = max(qi(i,k,j),0.0)
    466. 

  466.           qrs(i,k,1) = max(qr(i,k,j),0.0)
    467. 

  467.           qrs(i,k,2) = max(qqs(i,k,j),0.0)
    468. 

  468.         enddo
    469. 

  469. !
    470. 

  470. !----------------------------------------------------------------
    471. 

  471. !     latent heat for phase changes and heat capacity. neglect the
    472. 

  472. !     changes during microphysical process calculation
    473. 

  473. !     emanuel(1994)
    474. 

  474. !
    475. 

  475.       do k = kts, kte
    476. 

  476.           cpm(i,k) = cpmcal(q(i,k,j))
    477. 

  477.           xl(i,k) = xlcal(t(i,k,j))
    478. 

  478.       enddo
    479. 

  479. !
    480. 

  480. !----------------------------------------------------------------
    481. 

  481. !     compute the minor time steps.
    482. 

  482. !
    483. 

  483. !     loops = max(nint(delt/dtcldcr),1)
    484. 

  484. !     dtcld = delt/loops
    485. 

  485. !     if(delt.le.dtcldcr) dtcld = delt
    486. 

  486. !
    487. 

  487.       do loop = 1,loops
    488. 

  488. !
    489. 

  489. !----------------------------------------------------------------
    490. 

  490. !     initialize the large scale variables
    491. 

  491. !
    492. 

  492.         mstep = 1
    493. 

  493.         flgcld = .true.
    494. 

  494. !
    495. 

  495.       do k = kts, kte
    496. 

  496.           denfac(i,k) = sqrt(den0/den(i,k,j))
    497. 

  497.       enddo
    498. 

  498. !     do k = kts, kte
    499. 

  499. !       CALL VREC( tvec1(its), den(its,k,j), ite-its+1)
    500. 

  500. !       do i = its, ite
    501. 

  501. !         tvec1(i) = tvec1(i)*den0
    502. 

  502. !       enddo
    503. 

  503. !       CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
    504. 

  504. !     enddo
    505. 

  505. !
    506. 

  506. ! Inline expansion for fpvs
    507. 

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

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

  509.       hsub = xls
    510. 

  510.       hvap = xlv0
    511. 

  511.       cvap = cpv
    512. 

  512.       ttp=t0c+0.01
    513. 

  513.       dldt=cvap-cliq
    514. 

  514.       xa=-dldt/rv
    515. 

  515.       xb=xa+hvap/(rv*ttp)
    516. 

  516.       dldti=cvap-cice
    517. 

  517.       xai=-dldti/rv
    518. 

  518.       xbi=xai+hsub/(rv*ttp)
    519. 

  519. 

  520. ! this is for compilers where the conditional inhibits vectorization
    521. 

  521. #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
    522. 

  522.       do k = kts, kte
    523. 

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

  524.             xal = xai
    525. 

  525.             xbl = xbi
    526. 

  526.           else
    527. 

  527.             xal = xa
    528. 

  528.             xbl = xb
    529. 

  529.           endif
    530. 

  530.           tr=ttp/t(i,k,j)
    531. 

  531.           logtr=log(tr)
    532. 

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

  533.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
    534. 

  534.           qs(i,k,1) = max(qs(i,k,1),qmin)
    535. 

  535.           rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
    536. 

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

  537.           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
    538. 

  538.           qs(i,k,2) = max(qs(i,k,2),qmin)
    539. 

  539.           rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
    540. 

  540.       enddo
    541. 

  541. #else
    542. 

  542.       do k = kts, kte
    543. 

  543.           tr=ttp/t(i,k,j)
    544. 

  544.           logtr=log(tr)
    545. 

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

  546.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
    547. 

  547.           qs(i,k,1) = max(qs(i,k,1),qmin)
    548. 

  548.           rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
    549. 

  549.           if(t(i,k,j).lt.ttp) then
    550. 

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

  551.           else
    552. 

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

  553.           endif
    554. 

  554.           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
    555. 

  555.           qs(i,k,2) = max(qs(i,k,2),qmin)
    556. 

  556.           rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
    557. 

  557.       enddo
    558. 

  558. #endif
    559. 

  559. 

  560. !
    561. 

  561. !----------------------------------------------------------------
    562. 

  562. !     initialize the variables for microphysical physics
    563. 

  563. !
    564. 

  564. !
    565. 

  565.       do k = kts, kte
    566. 

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

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

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

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

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

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

  572.           psacw(i,k) = 0.
    573. 

  573.           pigen(i,k) = 0.
    574. 

  574.           pidep(i,k) = 0.
    575. 

  575.           pcond(i,k) = 0.
    576. 

  576.           psmlt(i,k) = 0.
    577. 

  577.           psevp(i,k) = 0.
    578. 

  578.           falk(i,k,1) = 0.
    579. 

  579.           falk(i,k,2) = 0.
    580. 

  580.           fall(i,k,1) = 0.
    581. 

  581.           fall(i,k,2) = 0.
    582. 

  582.           fallc(i,k) = 0.
    583. 

  583.           falkc(i,k) = 0.
    584. 

  584.           xni(i,k) = 1.e3
    585. 

  585.       enddo
    586. 

  586. !
    587. 

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

  588. !     compute the fallout term:
    589. 

  589. !     first, vertical terminal velosity for minor loops
    590. 

  590. !
    591. 

  591.       do k = kts, kte
    592. 

  592.           supcol = t0c-t(i,k,j)
    593. 

  593. !---------------------------------------------------------------
    594. 

  594. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
    595. 

  595. !---------------------------------------------------------------
    596. 

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

  597.           if(qrs(i,k,1).le.qcrmin)then
    598. 

  598.             rslope(i,k,1) = rslopermax
    599. 

  599.             rslopeb(i,k,1) = rsloperbmax
    600. 

  600.             rslope2(i,k,1) = rsloper2max
    601. 

  601.             rslope3(i,k,1) = rsloper3max
    602. 

  602.           else
    603. 

  603.             rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
    604. 

  604.             rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
    605. 

  605.             rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
    606. 

  606.             rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
    607. 

  607.           endif
    608. 

  608.           if(qrs(i,k,2).le.qcrmin)then
    609. 

  609.             rslope(i,k,2) = rslopesmax
    610. 

  610.             rslopeb(i,k,2) = rslopesbmax
    611. 

  611.             rslope2(i,k,2) = rslopes2max
    612. 

  612.             rslope3(i,k,2) = rslopes3max
    613. 

  613.           else
    614. 

  614.             rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
    615. 

  615.             rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
    616. 

  616.             rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
    617. 

  617.             rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
    618. 

  618.           endif
    619. 

  619. !-------------------------------------------------------------
    620. 

  620. ! Ni: ice crystal number concentraiton   [HDC 5c]
    621. 

  621. !-------------------------------------------------------------
    622. 

  622. !         xni(i,k) = min(max(5.38e7*(den(i,k,j)                                  &
    623. 

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

  624.           temp = (den(i,k,j)*max(qci(i,k,2),qmin))
    625. 

  625.           temp = sqrt(sqrt(temp*temp*temp))
    626. 

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

  627.       enddo
    628. 

  628. !
    629. 

  629.       numdt = 1
    630. 

  630.       do k = kte, kts, -1
    631. 

  631.           work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k,j)
    632. 

  632.           work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k,j)
    633. 

  633.           numdt = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
    634. 

  634.           if(numdt.ge.mstep) mstep = numdt
    635. 

  635.       enddo
    636. 

  636.         rmstep = 1./mstep
    637. 

  637. !
    638. 

  638.       do n = 1, mstep
    639. 

  639.         k = kte
    640. 

  640. !             falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
    641. 

  641. !             falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
    642. 

  642.               falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
    643. 

  643.               falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
    644. 

  644.               fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
    645. 

  645.               fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
    646. 

  646. !             qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k,j),0.)
    647. 

  647. !             qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k,j),0.)
    648. 

  648.               dtcldden = dtcld/den(i,k,j)
    649. 

  649.               qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
    650. 

  650.               qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
    651. 

  651. !           endif
    652. 

  652.         do k = kte-1, kts, -1
    653. 

  653.               falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
    654. 

  654.               falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
    655. 

  655.               fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
    656. 

  656.               fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
    657. 

  657.               dtcldden = dtcld/den(i,k,j)
    658. 

  658.               rdelz = 1./delz(i,k,j)
    659. 

  659.               qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1)           &
    660. 

  660.                           *delz(i,k+1,j)*rdelz)*dtcldden,0.)
    661. 

  661.               qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2)           &
    662. 

  662.                           *delz(i,k+1,j)*rdelz)*dtcldden,0.)
    663. 

  663.         enddo
    664. 

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

  665.               if(t(i,k,j).gt.t0c.and.qrs(i,k,2).gt.0.) then
    666. 

  666. !----------------------------------------------------------------
    667. 

  667. ! psmlt: melting of snow [HL A33] [RH83 A25]
    668. 

  668. !       (T>T0: S->R)
    669. 

  669. !----------------------------------------------------------------
    670. 

  670.                 xlf = xlf0
    671. 

  671. !               work2(i,k)= venfac(p(i,k),t(i,k,j),den(i,k,j))
    672. 

  672.                 work2(i,k)= (exp(log(((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j)))        &
    673. 

  673.                             /((t(i,k,j))+120.)/(den(i,k,j)))/(8.794e-5             &
    674. 

  674.                             *exp(log(t(i,k,j))*(1.81))/p(i,k,j))))                 &
    675. 

  675.                             *((.3333333)))/sqrt((1.496e-6*((t(i,k,j))            &
    676. 

  676.                             *sqrt(t(i,k,j)))/((t(i,k,j))+120.)/(den(i,k,j))))        &
    677. 

  677.                             *sqrt(sqrt(den0/(den(i,k,j)))))
    678. 

  678.                 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    679. 

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

  680. !                           *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
    681. 

  681. !                           *work2(i,k)*coeres)
    682. 

  682.                 psmlt(i,k) = (1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j)))        &         
    683. 

  683.                             /((t(i,k,j))+120.)/(den(i,k,j)) )*(den(i,k,j)))          &
    684. 

  684.                             /xlf*(t0c-t(i,k,j))*pi/2.                            &
    685. 

  685.                             *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
    686. 

  686.                             *work2(i,k)*coeres)
    687. 

  687.                 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep,                &
    688. 

  688.                             -qrs(i,k,2)/mstep),0.)
    689. 

  689.                 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
    690. 

  690.                 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
    691. 

  691.                 t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*psmlt(i,k)
    692. 

  692.               endif
    693. 

  693.         enddo
    694. 

  694.       enddo
    695. 

  695. 

  696. 

  697. !---------------------------------------------------------------
    698. 

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

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

  700.       mstep = 1
    701. 

  701.       numdt = 1
    702. 

  702.       do k = kte, kts, -1
    703. 

  703.           if(qci(i,k,2).le.0.) then
    704. 

  704.             work2c(i,k) = 0.
    705. 

  705.           else
    706. 

  706.             xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
    707. 

  707. !           diameter  = min(dicon * sqrt(xmi),dimax)
    708. 

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

  709.             work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
    710. 

  710.             work2c(i,k) = work1c(i,k)/delz(i,k,j)
    711. 

  711.           endif
    712. 

  712.           numdt = max(nint(work2c(i,k)*dtcld+.5),1)
    713. 

  713.           if(numdt.ge.mstep) mstep = numdt
    714. 

  714.       enddo
    715. 

  715. !
    716. 

  716.       do n = 1, mstep
    717. 

  717.         k = kte
    718. 

  718.             falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
    719. 

  719.             holdc = falkc(i,k)
    720. 

  720.             fallc(i,k) = fallc(i,k)+falkc(i,k)
    721. 

  721.             holdci = qci(i,k,2)
    722. 

  722.             qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k,j),0.)
    723. 

  723. !         endif
    724. 

  724.         do k = kte-1, kts, -1
    725. 

  725.               falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
    726. 

  726.               holdc = falkc(i,k)
    727. 

  727.               fallc(i,k) = fallc(i,k)+falkc(i,k)
    728. 

  728.               holdci = qci(i,k,2)
    729. 

  729.               qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1)             &
    730. 

  730.                           *delz(i,k+1,j)/delz(i,k,j))*dtcld/den(i,k,j),0.)
    731. 

  731. !           endif
    732. 

  732.         enddo
    733. 

  733.       enddo
    734. 

  734. !
    735. 

  735. !
    736. 

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

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

  738. !
    739. 

  739.         fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
    740. 

  740.         fallsum_qsi = fall(i,1,2)+fallc(i,1)
    741. 

  741.         rainncv(i,j) = 0.
    742. 

  742.         if(fallsum.gt.0.) then
    743. 

  743.           rainncv(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000.
    744. 

  744.           rain(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000. + rain(i,j)
    745. 

  745.         endif
    746. 

  746.         snowncv(i,j) = 0.
    747. 

  747.         if(fallsum_qsi.gt.0.) then
    748. 

  748.           snowncv(i,j) = fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000.
    749. 

  749.           snow(i,j) = fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000. + snow(i,j)
    750. 

  750.         endif
    751. 

  751.         sr(i,j) = 0.
    752. 

  752.         if(fallsum.gt.0.)sr(i,j)=fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000.        &    
    753. 

  753.         /(rainncv(i,j)+1.e-12)
    754. 

  754. !
    755. 

  755. !---------------------------------------------------------------
    756. 

  756. ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
    757. 

  757. !       (T>T0: I->C)
    758. 

  758. !---------------------------------------------------------------
    759. 

  759.       do k = kts, kte
    760. 

  760.           supcol = t0c-t(i,k,j)
    761. 

  761.           xlf = xls-xl(i,k)
    762. 

  762.           if(supcol.lt.0.) xlf = xlf0
    763. 

  763.           if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
    764. 

  764.             qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
    765. 

  765.             t(i,k,j) = t(i,k,j) - xlf/cpm(i,k)*qci(i,k,2)
    766. 

  766.             qci(i,k,2) = 0.
    767. 

  767.           endif
    768. 

  768. !---------------------------------------------------------------
    769. 

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

  770. !        (T<-40C: C->I)
    771. 

  771. !---------------------------------------------------------------
    772. 

  772.           if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
    773. 

  773.             qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
    774. 

  774.             t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*qci(i,k,1)
    775. 

  775.             qci(i,k,1) = 0.
    776. 

  776.           endif
    777. 

  777. !---------------------------------------------------------------
    778. 

  778. ! pihtf: heterogeneous freezing of cloud water [HL A44]
    779. 

  779. !        (T0>T>-40C: C->I)
    780. 

  780. !---------------------------------------------------------------
    781. 

  781.           if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
    782. 

  782.             supcolt=min(supcol,50.)
    783. 

  783. !           pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.)                         &
    784. 

  784. !              *den(i,k,j)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
    785. 

  785.             pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.)                        &
    786. 

  786.             *den(i,k,j)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
    787. 

  787.             qci(i,k,2) = qci(i,k,2) + pfrzdtc
    788. 

  788.             t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtc
    789. 

  789.             qci(i,k,1) = qci(i,k,1)-pfrzdtc
    790. 

  790.           endif
    791. 

  791. !---------------------------------------------------------------
    792. 

  792. ! psfrz: freezing of rain water [HL A20] [LFO 45]
    793. 

  793. !        (T<T0, R->S)
    794. 

  794. !---------------------------------------------------------------
    795. 

  795.           if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
    796. 

  796.             supcolt=min(supcol,50.)
    797. 

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

  798. !                 *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld,              &
    799. 

  799. !                 qrs(i,k,1))
    800. 

  800.             temp = rslope(i,k,1)
    801. 

  801.             temp = temp*temp*temp*temp*temp*temp*temp
    802. 

  802.             pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k,j)                  &
    803. 

  803.                   *(exp(pfrz2*supcolt)-1.)*temp*dtcld,                         &
    804. 

  804.                   qrs(i,k,1))
    805. 

  805.             qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
    806. 

  806.             t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtr
    807. 

  807.             qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
    808. 

  808.           endif
    809. 

  809.       enddo
    810. 

  810. !
    811. 

  811. !----------------------------------------------------------------
    812. 

  812. !     rsloper: reverse of the slope parameter of the rain(m)
    813. 

  813. !     xka:    thermal conductivity of air(jm-1s-1k-1)
    814. 

  814. !     work1:  the thermodynamic term in the denominator associated with
    815. 

  815. !             heat conduction and vapor diffusion
    816. 

  816. !             (ry88, y93, h85)
    817. 

  817. !     work2: parameter associated with the ventilation effects(y93)
    818. 

  818. !
    819. 

  819.       do k = kts, kte
    820. 

  820.           if(qrs(i,k,1).le.qcrmin)then
    821. 

  821.             rslope(i,k,1) = rslopermax
    822. 

  822.             rslopeb(i,k,1) = rsloperbmax
    823. 

  823.             rslope2(i,k,1) = rsloper2max
    824. 

  824.             rslope3(i,k,1) = rsloper3max
    825. 

  825.           else
    826. 

  826. !           rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
    827. 

  827.             rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k,j))))))
    828. 

  828.             rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
    829. 

  829.             rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
    830. 

  830.             rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
    831. 

  831.           endif
    832. 

  832.           if(qrs(i,k,2).le.qcrmin)then
    833. 

  833.             rslope(i,k,2) = rslopesmax
    834. 

  834.             rslopeb(i,k,2) = rslopesbmax
    835. 

  835.             rslope2(i,k,2) = rslopes2max
    836. 

  836.             rslope3(i,k,2) = rslopes3max
    837. 

  837.           else
    838. 

  838. !            rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
    839. 

  839.             rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2))   &  
    840. 

  840.                           *(den(i,k,j))))))
    841. 

  841.             rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
    842. 

  842.             rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
    843. 

  843.             rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
    844. 

  844.           endif
    845. 

  845.       enddo
    846. 

  846. !
    847. 

  847.       do k = kts, kte
    848. 

  848. !         work1(i,k,1) = diffac(xl(i,k),p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,1))
    849. 

  849.           work1(i,k,1) = ((((den(i,k,j))*(xl(i,k))*(xl(i,k)))*((t(i,k,j))+120.)    &
    850. 

  850.                        *(den(i,k,j)))/(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))&
    851. 

  851.                        *(den(i,k,j))*(rv*(t(i,k,j))*(t(i,k,j)))))                    &
    852. 

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

  853. !         work1(i,k,2) = diffac(xls,p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,2))
    854. 

  854.           work1(i,k,2) = ((((den(i,k,j))*(xls)*(xls))*((t(i,k,j))+120.)*(den(i,k,j)))&
    855. 

  855.                        /(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))*(den(i,k,j)) &
    856. 

  856.                        *(rv*(t(i,k,j))*(t(i,k,j))))                                &
    857. 

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

  858. !         work2(i,k) = venfac(p(i,k,j),t(i,k,j),den(i,k,j))
    859. 

  859.           work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k,j))*sqrt(t(i,k,j)))) &
    860. 

  860.                      *p(i,k,j))/(((t(i,k,j))+120.)*den(i,k,j)*(8.794e-5              &
    861. 

  861.                      *exp(log(t(i,k,j))*(1.81))))))*sqrt(sqrt(den0/(den(i,k,j))))) &
    862. 

  862.                       /sqrt((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))                 &
    863. 

  863.                       /(((t(i,k,j))+120.)*den(i,k,j)))
    864. 

  864.       enddo 
    865. 

  865. !
    866. 

  866. !===============================================================
    867. 

  867. !
    868. 

  868. ! warm rain processes
    869. 

  869. !
    870. 

  870. ! - follows the processes in RH83 and LFO except for autoconcersion
    871. 

  871. !
    872. 

  872. !===============================================================
    873. 

  873. !
    874. 

  874.       do k = kts, kte
    875. 

  875.           supsat = max(q(i,k,j),qmin)-qs(i,k,1)
    876. 

  876.           satdt = supsat/dtcld
    877. 

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

  878. ! praut: auto conversion rate from cloud to rain [HDC 16]
    879. 

  879. !        (C->R)
    880. 

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

  881.           if(qci(i,k,1).gt.qc0) then
    882. 

  882.             praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
    883. 

  883.             praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
    884. 

  884.           endif
    885. 

  885. !---------------------------------------------------------------
    886. 

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

  887. !        (C->R)
    888. 

  888. !---------------------------------------------------------------
    889. 

  889.           if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
    890. 

  890.             pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1)               & 
    891. 

  891.                          *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
    892. 

  892.           endif
    893. 

  893. !---------------------------------------------------------------
    894. 

  894. ! prevp: evaporation/condensation rate of rain [HDC 14]
    895. 

  895. !        (V->R or R->V)
    896. 

  896. !---------------------------------------------------------------
    897. 

  897.           if(qrs(i,k,1).gt.0.) then
    898. 

  898.             coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
    899. 

  899.             prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1)                 &
    900. 

  900.                          +precr2*work2(i,k)*coeres)/work1(i,k,1)
    901. 

  901.             if(prevp(i,k).lt.0.) then
    902. 

  902.               prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
    903. 

  903.               prevp(i,k) = max(prevp(i,k),satdt/2)
    904. 

  904.             else
    905. 

  905.               prevp(i,k) = min(prevp(i,k),satdt/2)
    906. 

  906.             endif
    907. 

  907.           endif
    908. 

  908.       enddo
    909. 

  909. !
    910. 

  910. !===============================================================
    911. 

  911. !
    912. 

  912. ! cold rain processes
    913. 

  913. !
    914. 

  914. ! - follows the revised ice microphysics processes in HDC
    915. 

  915. ! - the processes same as in RH83 and RH84  and LFO behave
    916. 

  916. !   following ice crystal hapits defined in HDC, inclduing
    917. 

  917. !   intercept parameter for snow (n0s), ice crystal number
    918. 

  918. !   concentration (ni), ice nuclei number concentration
    919. 

  919. !   (n0i), ice diameter (d)
    920. 

  920. !
    921. 

  921. !===============================================================
    922. 

  922. !
    923. 

  923.       rdtcld = 1./dtcld
    924. 

  924.       do k = kts, kte
    925. 

  925.           supcol = t0c-t(i,k,j)
    926. 

  926.           supsat = max(q(i,k,j),qmin)-qs(i,k,2)
    927. 

  927.           satdt = supsat/dtcld
    928. 

  928.           ifsat = 0
    929. 

  929. !-------------------------------------------------------------
    930. 

  930. ! Ni: ice crystal number concentraiton   [HDC 5c]
    931. 

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

  932. !         xni(i,k) = min(max(5.38e7*(den(i,k,j)                                  &
    933. 

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

  934.           temp = (den(i,k,j)*max(qci(i,k,2),qmin))
    935. 

  935.           temp = sqrt(sqrt(temp*temp*temp))
    936. 

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

  937.           eacrs = exp(0.07*(-supcol))
    938. 

  938. !
    939. 

  939.           if(supcol.gt.0) then
    940. 

  940.             if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
    941. 

  941.               xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
    942. 

  942.               diameter  = min(dicon * sqrt(xmi),dimax)
    943. 

  943.               vt2i = 1.49e4*diameter**1.31
    944. 

  944.               vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
    945. 

  945. !-------------------------------------------------------------
    946. 

  946. ! psaci: Accretion of cloud ice by rain [HDC 10]
    947. 

  947. !        (T<T0: I->S)
    948. 

  948. !-------------------------------------------------------------
    949. 

  949.               acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2)            &
    950. 

  950.                       +diameter**2*rslope(i,k,2)
    951. 

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

  952.                            *abs(vt2s-vt2i)*acrfac/4.
    953. 

  953.             endif
    954. 

  954.           endif
    955. 

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

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

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

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

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

  960.             psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)                  &
    961. 

  961.                            *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k)              &
    962. 

  962. !                          ,qci(i,k,1)/dtcld)
    963. 

  963.                            ,qci(i,k,1)*rdtcld)
    964. 

  964.           endif
    965. 

  965.           if(supcol .gt. 0) then
    966. 

  966. !-------------------------------------------------------------
    967. 

  967. ! pidep: Deposition/Sublimation rate of ice [HDC 9]
    968. 

  968. !       (T<T0: V->I or I->V)
    969. 

  969. !-------------------------------------------------------------
    970. 

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

  971.               xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
    972. 

  972.               diameter = dicon * sqrt(xmi)
    973. 

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

  974.               supice = satdt-prevp(i,k)
    975. 

  975.               if(pidep(i,k).lt.0.) then
    976. 

  976. !               pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
    977. 

  977. !               pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
    978. 

  978.                 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
    979. 

  979.                 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
    980. 

  980.               else
    981. 

  981. !               pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
    982. 

  982.                 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
    983. 

  983.               endif
    984. 

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

  985.             endif
    986. 

  986. !-------------------------------------------------------------
    987. 

  987. ! psdep: deposition/sublimation rate of snow [HDC 14]
    988. 

  988. !        (V->S or S->V)
    989. 

  989. !-------------------------------------------------------------
    990. 

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

  991.               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    992. 

  992.               psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)                          &
    993. 

  993.                            *(precs1*rslope2(i,k,2)+precs2                      &
    994. 

  994.                            *work2(i,k)*coeres)/work1(i,k,2)
    995. 

  995.               supice = satdt-prevp(i,k)-pidep(i,k)
    996. 

  996.               if(psdep(i,k).lt.0.) then
    997. 

  997. !               psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
    998. 

  998. !               psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
    999. 

  999.                 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
    1000. 

  1000.                 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
    1001. 

  1001.               else
    1002. 

  1002. !             psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
    1003. 

  1003.                 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
    1004. 

  1004.               endif
    1005. 

  1005.               if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt))          &
    1006. 

  1006.                 ifsat = 1
    1007. 

  1007.             endif
    1008. 

  1008. !-------------------------------------------------------------
    1009. 

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

  1010. !       (T<T0: V->I)
    1011. 

  1011. !-------------------------------------------------------------
    1012. 

  1012.             if(supsat.gt.0.and.ifsat.ne.1) then
    1013. 

  1013.               supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
    1014. 

  1014.               xni0 = 1.e3*exp(0.1*supcol)
    1015. 

  1015.               roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
    1016. 

  1016.               pigen(i,k) = max(0.,(roqi0/den(i,k,j)-max(qci(i,k,2),0.))          &
    1017. 

  1017. !                        /dtcld)
    1018. 

  1018.                          *rdtcld)
    1019. 

  1019.               pigen(i,k) = min(min(pigen(i,k),satdt),supice)
    1020. 

  1020.             endif
    1021. 

  1021. !
    1022. 

  1022. !-------------------------------------------------------------
    1023. 

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

  1024. !       (T<T0: I->S)
    1025. 

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

  1026.             if(qci(i,k,2).gt.0.) then
    1027. 

  1027.               qimax = roqimax/den(i,k,j)
    1028. 

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

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

  1030.             endif
    1031. 

  1031.           endif
    1032. 

  1032. !-------------------------------------------------------------
    1033. 

  1033. ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
    1034. 

  1034. !       (T>T0: S->V)
    1035. 

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

  1036.           if(supcol.lt.0.) then
    1037. 

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

  1038.               psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
    1039. 

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

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

  1041.           endif
    1042. 

  1042.       enddo
    1043. 

  1043. !
    1044. 

  1044. !
    1045. 

  1045. !----------------------------------------------------------------
    1046. 

  1046. !     check mass conservation of generation terms and feedback to the
    1047. 

  1047. !     large scale
    1048. 

  1048. !
    1049. 

  1049.       do k = kts, kte
    1050. 

  1050.           if(t(i,k,j).le.t0c) then
    1051. 

  1051. !
    1052. 

  1052. !     cloud water
    1053. 

  1053. !
    1054. 

  1054.             value = max(qmin,qci(i,k,1))
    1055. 

  1055.             source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
    1056. 

  1056.             if (source.gt.value) then
    1057. 

  1057.               factor = value/source
    1058. 

  1058.               praut(i,k) = praut(i,k)*factor
    1059. 

  1059.               pracw(i,k) = pracw(i,k)*factor
    1060. 

  1060.               psacw(i,k) = psacw(i,k)*factor
    1061. 

  1061.             endif
    1062. 

  1062. !
    1063. 

  1063. !     cloud ice
    1064. 

  1064. !
    1065. 

  1065.             value = max(qmin,qci(i,k,2))
    1066. 

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

  1067.             if (source.gt.value) then
    1068. 

  1068.               factor = value/source
    1069. 

  1069.               psaut(i,k) = psaut(i,k)*factor
    1070. 

  1070.               psaci(i,k) = psaci(i,k)*factor
    1071. 

  1071.               pigen(i,k) = pigen(i,k)*factor
    1072. 

  1072.               pidep(i,k) = pidep(i,k)*factor
    1073. 

  1073.             endif
    1074. 

  1074. !
    1075. 

  1075. !     rain
    1076. 

  1076. !
    1077. 

  1077. !
    1078. 

  1078.             value = max(qmin,qrs(i,k,1))
    1079. 

  1079.             source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
    1080. 

  1080.             if (source.gt.value) then
    1081. 

  1081.               factor = value/source
    1082. 

  1082.               praut(i,k) = praut(i,k)*factor
    1083. 

  1083.               pracw(i,k) = pracw(i,k)*factor
    1084. 

  1084.               prevp(i,k) = prevp(i,k)*factor
    1085. 

  1085.             endif
    1086. 

  1086. !
    1087. 

  1087. !    snow
    1088. 

  1088. !
    1089. 

  1089.             value = max(qmin,qrs(i,k,2))
    1090. 

  1090.             source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld  
    1091. 

  1091.             if (source.gt.value) then
    1092. 

  1092.               factor = value/source
    1093. 

  1093.               psdep(i,k) = psdep(i,k)*factor
    1094. 

  1094.               psaut(i,k) = psaut(i,k)*factor
    1095. 

  1095.               psaci(i,k) = psaci(i,k)*factor
    1096. 

  1096.               psacw(i,k) = psacw(i,k)*factor
    1097. 

  1097.             endif
    1098. 

  1098. !
    1099. 

  1099.             work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
    1100. 

  1100. !     update
    1101. 

  1101.             q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
    1102. 

  1102.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
    1103. 

  1103.                         +psacw(i,k))*dtcld,0.)
    1104. 

  1104.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
    1105. 

  1105.                         +prevp(i,k))*dtcld,0.)
    1106. 

  1106.             qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k)                 &
    1107. 

  1107.                         -pigen(i,k)-pidep(i,k))*dtcld,0.)
    1108. 

  1108.             qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)                 &
    1109. 

  1109.                         +psaci(i,k)+psacw(i,k))*dtcld,0.)
    1110. 

  1110.             xlf = xls-xl(i,k)
    1111. 

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

  1112.                       -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
    1113. 

  1113.             t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
    1114. 

  1114.           else
    1115. 

  1115. !
    1116. 

  1116. !     cloud water
    1117. 

  1117. !
    1118. 

  1118.             value = max(qmin,qci(i,k,1))
    1119. 

  1119.             source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
    1120. 

  1120.             if (source.gt.value) then
    1121. 

  1121.               factor = value/source
    1122. 

  1122.               praut(i,k) = praut(i,k)*factor
    1123. 

  1123.               pracw(i,k) = pracw(i,k)*factor
    1124. 

  1124.               psacw(i,k) = psacw(i,k)*factor
    1125. 

  1125.             endif
    1126. 

  1126. !
    1127. 

  1127. !     rain
    1128. 

  1128. !
    1129. 

  1129.             value = max(qmin,qrs(i,k,1))
    1130. 

  1130.             source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
    1131. 

  1131.             if (source.gt.value) then
    1132. 

  1132.               factor = value/source
    1133. 

  1133.               praut(i,k) = praut(i,k)*factor
    1134. 

  1134.               pracw(i,k) = pracw(i,k)*factor
    1135. 

  1135.               prevp(i,k) = prevp(i,k)*factor
    1136. 

  1136.               psacw(i,k) = psacw(i,k)*factor
    1137. 

  1137.             endif  
    1138. 

  1138. !
    1139. 

  1139. !     snow
    1140. 

  1140. !
    1141. 

  1141.             value = max(qcrmin,qrs(i,k,2))
    1142. 

  1142.             source=(-psevp(i,k))*dtcld
    1143. 

  1143.             if (source.gt.value) then
    1144. 

  1144.               factor = value/source
    1145. 

  1145.               psevp(i,k) = psevp(i,k)*factor
    1146. 

  1146.             endif
    1147. 

  1147.             work2(i,k)=-(prevp(i,k)+psevp(i,k))
    1148. 

  1148. !     update
    1149. 

  1149.             q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
    1150. 

  1150.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
    1151. 

  1151.                         +psacw(i,k))*dtcld,0.)
    1152. 

  1152.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
    1153. 

  1153.                         +prevp(i,k) +psacw(i,k))*dtcld,0.)
    1154. 

  1154.             qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
    1155. 

  1155.             xlf = xls-xl(i,k)
    1156. 

  1156.             xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
    1157. 

  1157.             t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
    1158. 

  1158.           endif
    1159. 

  1159.       enddo
    1160. 

  1160. !
    1161. 

  1161. ! Inline expansion for fpvs
    1162. 

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

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

  1164.       hsub = xls
    1165. 

  1165.       hvap = xlv0
    1166. 

  1166.       cvap = cpv
    1167. 

  1167.       ttp=t0c+0.01
    1168. 

  1168.       dldt=cvap-cliq
    1169. 

  1169.       xa=-dldt/rv
    1170. 

  1170.       xb=xa+hvap/(rv*ttp)
    1171. 

  1171.       dldti=cvap-cice
    1172. 

  1172.       xai=-dldti/rv
    1173. 

  1173.       xbi=xai+hsub/(rv*ttp)
    1174. 

  1174.       do k = kts, kte
    1175. 

  1175.           tr=ttp/t(i,k,j)
    1176. 

  1176.           logtr = log(tr)
    1177. 

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

  1178.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
    1179. 

  1179.           qs(i,k,1) = max(qs(i,k,1),qmin)
    1180. 

  1180.       enddo
    1181. 

  1181. !
    1182. 

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

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

  1184. !     if there exists additional water vapor condensated/if
    1185. 

  1185. !     evaporation of cloud water is not enough to remove subsaturation
    1186. 

  1186. !
    1187. 

  1187.       do k = kts, kte
    1188. 

  1188. !         work1(i,k,1) = conden(t(i,k,j),q(i,k,j),qs(i,k,1),xl(i,k),cpm(i,k))
    1189. 

  1189.           work1(i,k,1) = ((max(q(i,k,j),qmin)-(qs(i,k,1)))/(1.+(xl(i,k))         &   
    1190. 

  1190.                         *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))                 &
    1191. 

  1191.                         /((t(i,k,j))*(t(i,k,j)))))
    1192. 

  1192.           work2(i,k) = qci(i,k,1)+work1(i,k,1)
    1193. 

  1193.           pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k,j),0.)/dtcld)
    1194. 

  1194.           if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.)                          &
    1195. 

  1195.             pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
    1196. 

  1196.           q(i,k,j) = q(i,k,j)-pcond(i,k)*dtcld
    1197. 

  1197.           qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
    1198. 

  1198.           t(i,k,j) = t(i,k,j)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
    1199. 

  1199.       enddo
    1200. 

  1200. !
    1201. 

  1201. !
    1202. 

  1202. !----------------------------------------------------------------
    1203. 

  1203. !     padding for small values
    1204. 

  1204. !
    1205. 

  1205.       do k = kts, kte
    1206. 

  1206.           if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
    1207. 

  1207.           if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
    1208. 

  1208.       enddo
    1209. 

  1209.       enddo                  ! big loops
    1210. 

  1210. 

  1211.          DO K=kts,kte
    1212. 

  1212.             th(i,k,j)=t(i,k,j)/pii(i,k,j)
    1213. 

  1213.             qc(i,k,j) = qci(i,k,1)
    1214. 

  1214.             qi(i,k,j) = qci(i,k,2)
    1215. 

  1215.             qr(i,k,j) = qrs(i,k,1)
    1216. 

  1216.             qqs(i,k,j) = qrs(i,k,2)
    1217. 

  1217.          ENDDO
    1218. 

  1218. 

  1219. 

  1220. 

  1221.       ENDDO     !  i loop
    1222. 

  1222.       enddo     !  j loop
    1223. 

  1223. !$acc end region
    1224. 

  1224. 

  1225.     ELSE
    1226. 

  1226. 

  1227. !
    1228. 

  1228. ! Moved outside of accelerator region
    1229. 

  1229. !
    1230. 

  1230.       loops = max(nint(delt/dtcldcr),1)
    1231. 

  1231.       dtcld = delt/loops
    1232. 

  1232.       if(delt.le.dtcldcr) dtcld = delt
    1233. 

  1233. 

  1234. !$acc region &
    1235. 

  1235. !$acc        local(t) &
    1236. 

  1236. !$acc        copyin(delz(:,:,:),p(:,:,:),den(:,:,:),pii(:,:,:)) &
    1237. 

  1237. !$acc        copyout(rainncv(:,:),sr(:,:)) &
    1238. 

  1238. !$acc        copy(qqs(:,:,:),qr(:,:,:),qi(:,:,:),qc(:,:,:)) &
    1239. 

  1239. !$acc        copy(th(:,:,:),q(:,:,:),rain(:,:))
    1240. 

  1240. !$acc do &
    1241. 

  1241. !$acc        private(rh,qs,rslope,rslope2,rslope3,rslopeb,falk,fall) &
    1242. 

  1242. !$acc        private(work1,qci,qrs,falkc,fallc,xl,cpm,denfac,xni) &
    1243. 

  1243. !$acc        private(n0sfac,work2,work1c,work2c,pigen,pidep,psdep) &
    1244. 

  1244. !$acc        private(praut,psaut,prevp,psevp) &
    1245. 

  1245. !$acc        private(pracw,psacw,psaci,pcond,psmlt) &
    1246. 

  1246. !$acc        parallel
    1247. 

  1247.       do j = jts, jte
    1248. 

  1248. !$acc do &
    1249. 

  1249. !$acc        private(numdt,mstep) &
    1250. 

  1250. !$acc        kernel vector
    1251. 

  1251.         do i = its, ite
    1252. 

  1252.         do k = kts, kte
    1253. 

  1253.           t(i,k,j)=th(i,k,j)*pii(i,k,j)
    1254. 

  1254.           qci(i,k,1) = max(qc(i,k,j),0.0)
    1255. 

  1255.           qci(i,k,2) = max(qi(i,k,j),0.0)
    1256. 

  1256.           qrs(i,k,1) = max(qr(i,k,j),0.0)
    1257. 

  1257.           qrs(i,k,2) = max(qqs(i,k,j),0.0)
    1258. 

  1258.         enddo
    1259. 

  1259. !
    1260. 

  1260. !----------------------------------------------------------------
    1261. 

  1261. !     latent heat for phase changes and heat capacity. neglect the
    1262. 

  1262. !     changes during microphysical process calculation
    1263. 

  1263. !     emanuel(1994)
    1264. 

  1264. !
    1265. 

  1265.       do k = kts, kte
    1266. 

  1266.           cpm(i,k) = cpmcal(q(i,k,j))
    1267. 

  1267.           xl(i,k) = xlcal(t(i,k,j))
    1268. 

  1268.       enddo
    1269. 

  1269. !
    1270. 

  1270. !----------------------------------------------------------------
    1271. 

  1271. !     compute the minor time steps.
    1272. 

  1272. !
    1273. 

  1273. !     loops = max(nint(delt/dtcldcr),1)
    1274. 

  1274. !     dtcld = delt/loops
    1275. 

  1275. !     if(delt.le.dtcldcr) dtcld = delt
    1276. 

  1276. !
    1277. 

  1277.       do loop = 1,loops
    1278. 

  1278. !
    1279. 

  1279. !----------------------------------------------------------------
    1280. 

  1280. !     initialize the large scale variables
    1281. 

  1281. !
    1282. 

  1282.         mstep = 1
    1283. 

  1283.         flgcld = .true.
    1284. 

  1284. !
    1285. 

  1285.       do k = kts, kte
    1286. 

  1286.           denfac(i,k) = sqrt(den0/den(i,k,j))
    1287. 

  1287.       enddo
    1288. 

  1288. !     do k = kts, kte
    1289. 

  1289. !       CALL VREC( tvec1(its), den(its,k,j), ite-its+1)
    1290. 

  1290. !       do i = its, ite
    1291. 

  1291. !         tvec1(i) = tvec1(i)*den0
    1292. 

  1292. !       enddo
    1293. 

  1293. !       CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
    1294. 

  1294. !     enddo
    1295. 

  1295. !
    1296. 

  1296. ! Inline expansion for fpvs
    1297. 

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

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

  1299.       hsub = xls
    1300. 

  1300.       hvap = xlv0
    1301. 

  1301.       cvap = cpv
    1302. 

  1302.       ttp=t0c+0.01
    1303. 

  1303.       dldt=cvap-cliq
    1304. 

  1304.       xa=-dldt/rv
    1305. 

  1305.       xb=xa+hvap/(rv*ttp)
    1306. 

  1306.       dldti=cvap-cice
    1307. 

  1307.       xai=-dldti/rv
    1308. 

  1308.       xbi=xai+hsub/(rv*ttp)
    1309. 

  1309. 

  1310. ! this is for compilers where the conditional inhibits vectorization
    1311. 

  1311. #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
    1312. 

  1312.       do k = kts, kte
    1313. 

  1313.           if(t(i,k,j).lt.ttp) then
    1314. 

  1314.             xal = xai
    1315. 

  1315.             xbl = xbi
    1316. 

  1316.           else
    1317. 

  1317.             xal = xa
    1318. 

  1318.             xbl = xb
    1319. 

  1319.           endif
    1320. 

  1320.           tr=ttp/t(i,k,j)
    1321. 

  1321.           logtr=log(tr)
    1322. 

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

  1323.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
    1324. 

  1324.           qs(i,k,1) = max(qs(i,k,1),qmin)
    1325. 

  1325.           rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
    1326. 

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

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

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

  1329.           rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
    1330. 

  1330.       enddo
    1331. 

  1331. #else
    1332. 

  1332.       do k = kts, kte
    1333. 

  1333.           tr=ttp/t(i,k,j)
    1334. 

  1334.           logtr=log(tr)
    1335. 

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

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

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

  1338.           rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
    1339. 

  1339.           if(t(i,k,j).lt.ttp) then
    1340. 

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

  1341.           else
    1342. 

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

  1343.           endif
    1344. 

  1344.           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
    1345. 

  1345.           qs(i,k,2) = max(qs(i,k,2),qmin)
    1346. 

  1346.           rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
    1347. 

  1347.       enddo
    1348. 

  1348. #endif
    1349. 

  1349. !
    1350. 

  1350. !----------------------------------------------------------------
    1351. 

  1351. !     initialize the variables for microphysical physics
    1352. 

  1352. !
    1353. 

  1353. !
    1354. 

  1354.       do k = kts, kte
    1355. 

  1355.           prevp(i,k) = 0.
    1356. 

  1356.           psdep(i,k) = 0.
    1357. 

  1357.           praut(i,k) = 0.
    1358. 

  1358.           psaut(i,k) = 0.
    1359. 

  1359.           pracw(i,k) = 0.
    1360. 

  1360.           psaci(i,k) = 0.
    1361. 

  1361.           psacw(i,k) = 0.
    1362. 

  1362.           pigen(i,k) = 0.
    1363. 

  1363.           pidep(i,k) = 0.
    1364. 

  1364.           pcond(i,k) = 0.
    1365. 

  1365.           psmlt(i,k) = 0.
    1366. 

  1366.           psevp(i,k) = 0.
    1367. 

  1367.           falk(i,k,1) = 0.
    1368. 

  1368.           falk(i,k,2) = 0.
    1369. 

  1369.           fall(i,k,1) = 0.
    1370. 

  1370.           fall(i,k,2) = 0.
    1371. 

  1371.           fallc(i,k) = 0.
    1372. 

  1372.           falkc(i,k) = 0.
    1373. 

  1373.           xni(i,k) = 1.e3
    1374. 

  1374.       enddo
    1375. 

  1375. !
    1376. 

  1376. !----------------------------------------------------------------
    1377. 

  1377. !     compute the fallout term:
    1378. 

  1378. !     first, vertical terminal velosity for minor loops
    1379. 

  1379. !
    1380. 

  1380.       do k = kts, kte
    1381. 

  1381.           supcol = t0c-t(i,k,j)
    1382. 

  1382. !---------------------------------------------------------------
    1383. 

  1383. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
    1384. 

  1384. !---------------------------------------------------------------
    1385. 

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

  1386.           if(qrs(i,k,1).le.qcrmin)then
    1387. 

  1387.             rslope(i,k,1) = rslopermax
    1388. 

  1388.             rslopeb(i,k,1) = rsloperbmax
    1389. 

  1389.             rslope2(i,k,1) = rsloper2max
    1390. 

  1390.             rslope3(i,k,1) = rsloper3max
    1391. 

  1391.           else
    1392. 

  1392.             rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
    1393. 

  1393.             rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
    1394. 

  1394.             rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
    1395. 

  1395.             rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
    1396. 

  1396.           endif
    1397. 

  1397.           if(qrs(i,k,2).le.qcrmin)then
    1398. 

  1398.             rslope(i,k,2) = rslopesmax
    1399. 

  1399.             rslopeb(i,k,2) = rslopesbmax
    1400. 

  1400.             rslope2(i,k,2) = rslopes2max
    1401. 

  1401.             rslope3(i,k,2) = rslopes3max
    1402. 

  1402.           else
    1403. 

  1403.             rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
    1404. 

  1404.             rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
    1405. 

  1405.             rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
    1406. 

  1406.             rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
    1407. 

  1407.           endif
    1408. 

  1408. !-------------------------------------------------------------
    1409. 

  1409. ! Ni: ice crystal number concentraiton   [HDC 5c]
    1410. 

  1410. !-------------------------------------------------------------
    1411. 

  1411. !         xni(i,k) = min(max(5.38e7*(den(i,k,j)                                  &
    1412. 

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

  1413.           temp = (den(i,k,j)*max(qci(i,k,2),qmin))
    1414. 

  1414.           temp = sqrt(sqrt(temp*temp*temp))
    1415. 

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

  1416.       enddo
    1417. 

  1417. !
    1418. 

  1418.       numdt = 1
    1419. 

  1419.       do k = kte, kts, -1
    1420. 

  1420.           work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k,j)
    1421. 

  1421.           work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k,j)
    1422. 

  1422.           numdt = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
    1423. 

  1423.           if(numdt.ge.mstep) mstep = numdt
    1424. 

  1424.       enddo
    1425. 

  1425.         rmstep = 1./mstep
    1426. 

  1426. !
    1427. 

  1427.       do n = 1, mstep
    1428. 

  1428.         k = kte
    1429. 

  1429. !             falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
    1430. 

  1430. !             falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
    1431. 

  1431.               falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
    1432. 

  1432.               falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
    1433. 

  1433.               fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
    1434. 

  1434.               fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
    1435. 

  1435. !             qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k,j),0.)
    1436. 

  1436. !             qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k,j),0.)
    1437. 

  1437.               dtcldden = dtcld/den(i,k,j)
    1438. 

  1438.               qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
    1439. 

  1439.               qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
    1440. 

  1440. !           endif
    1441. 

  1441.         do k = kte-1, kts, -1
    1442. 

  1442.               falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
    1443. 

  1443.               falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
    1444. 

  1444.               fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
    1445. 

  1445.               fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
    1446. 

  1446.               dtcldden = dtcld/den(i,k,j)
    1447. 

  1447.               rdelz = 1./delz(i,k,j)
    1448. 

  1448.               qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1)           &
    1449. 

  1449.                           *delz(i,k+1,j)*rdelz)*dtcldden,0.)
    1450. 

  1450.               qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2)           &
    1451. 

  1451.                           *delz(i,k+1,j)*rdelz)*dtcldden,0.)
    1452. 

  1452.         enddo
    1453. 

  1453.         do k = kte, kts, -1
    1454. 

  1454.               if(t(i,k,j).gt.t0c.and.qrs(i,k,2).gt.0.) then
    1455. 

  1455. !----------------------------------------------------------------
    1456. 

  1456. ! psmlt: melting of snow [HL A33] [RH83 A25]
    1457. 

  1457. !       (T>T0: S->R)
    1458. 

  1458. !----------------------------------------------------------------
    1459. 

  1459.                 xlf = xlf0
    1460. 

  1460. !               work2(i,k)= venfac(p(i,k),t(i,k,j),den(i,k,j))
    1461. 

  1461.                 work2(i,k)= (exp(log(((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j)))        &
    1462. 

  1462.                             /((t(i,k,j))+120.)/(den(i,k,j)))/(8.794e-5             &
    1463. 

  1463.                             *exp(log(t(i,k,j))*(1.81))/p(i,k,j))))                 &
    1464. 

  1464.                             *((.3333333)))/sqrt((1.496e-6*((t(i,k,j))            &
    1465. 

  1465.                             *sqrt(t(i,k,j)))/((t(i,k,j))+120.)/(den(i,k,j))))        &
    1466. 

  1466.                             *sqrt(sqrt(den0/(den(i,k,j)))))
    1467. 

  1467.                 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    1468. 

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

  1469. !                           *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
    1470. 

  1470. !                           *work2(i,k)*coeres)
    1471. 

  1471.                 psmlt(i,k) = (1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j)))        &         
    1472. 

  1472.                             /((t(i,k,j))+120.)/(den(i,k,j)) )*(den(i,k,j)))          &
    1473. 

  1473.                             /xlf*(t0c-t(i,k,j))*pi/2.                            &
    1474. 

  1474.                             *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
    1475. 

  1475.                             *work2(i,k)*coeres)
    1476. 

  1476.                 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep,                &
    1477. 

  1477.                             -qrs(i,k,2)/mstep),0.)
    1478. 

  1478.                 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
    1479. 

  1479.                 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
    1480. 

  1480.                 t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*psmlt(i,k)
    1481. 

  1481.               endif
    1482. 

  1482.         enddo
    1483. 

  1483.       enddo
    1484. 

  1484. !---------------------------------------------------------------
    1485. 

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

  1486. !---------------------------------------------------------------
    1487. 

  1487.       mstep = 1
    1488. 

  1488.       numdt = 1
    1489. 

  1489.       do k = kte, kts, -1
    1490. 

  1490.           if(qci(i,k,2).le.0.) then
    1491. 

  1491.             work2c(i,k) = 0.
    1492. 

  1492.           else
    1493. 

  1493.             xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
    1494. 

  1494. !           diameter  = min(dicon * sqrt(xmi),dimax)
    1495. 

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

  1496.             work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
    1497. 

  1497.             work2c(i,k) = work1c(i,k)/delz(i,k,j)
    1498. 

  1498.           endif
    1499. 

  1499.           numdt = max(nint(work2c(i,k)*dtcld+.5),1)
    1500. 

  1500.           if(numdt.ge.mstep) mstep = numdt
    1501. 

  1501.       enddo
    1502. 

  1502. !
    1503. 

  1503.       do n = 1, mstep
    1504. 

  1504.         k = kte
    1505. 

  1505.             falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
    1506. 

  1506.             holdc = falkc(i,k)
    1507. 

  1507.             fallc(i,k) = fallc(i,k)+falkc(i,k)
    1508. 

  1508.             holdci = qci(i,k,2)
    1509. 

  1509.             qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k,j),0.)
    1510. 

  1510.         do k = kte-1, kts, -1
    1511. 

  1511.               falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
    1512. 

  1512.               holdc = falkc(i,k)
    1513. 

  1513.               fallc(i,k) = fallc(i,k)+falkc(i,k)
    1514. 

  1514.               holdci = qci(i,k,2)
    1515. 

  1515.               qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1)             &
    1516. 

  1516.                           *delz(i,k+1,j)/delz(i,k,j))*dtcld/den(i,k,j),0.)
    1517. 

  1517.         enddo
    1518. 

  1518.       enddo
    1519. 

  1519. !
    1520. 

  1520. !
    1521. 

  1521. !----------------------------------------------------------------
    1522. 

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

  1523. !
    1524. 

  1524.         fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
    1525. 

  1525.         fallsum_qsi = fall(i,1,2)+fallc(i,1)
    1526. 

  1526.         rainncv(i,j) = 0.
    1527. 

  1527.         if(fallsum.gt.0.) then
    1528. 

  1528.           rainncv(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000.
    1529. 

  1529.           rain(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000. + rain(i,j)
    1530. 

  1530.         endif
    1531. 

  1531.         sr(i,j) = 0.
    1532. 

  1532.         if(fallsum.gt.0.)sr(i,j)=fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000.        &    
    1533. 

  1533.         /(rainncv(i,j)+1.e-12)
    1534. 

  1534. !
    1535. 

  1535. !---------------------------------------------------------------
    1536. 

  1536. ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
    1537. 

  1537. !       (T>T0: I->C)
    1538. 

  1538. !---------------------------------------------------------------
    1539. 

  1539.       do k = kts, kte
    1540. 

  1540.           supcol = t0c-t(i,k,j)
    1541. 

  1541.           xlf = xls-xl(i,k)
    1542. 

  1542.           if(supcol.lt.0.) xlf = xlf0
    1543. 

  1543.           if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
    1544. 

  1544.             qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
    1545. 

  1545.             t(i,k,j) = t(i,k,j) - xlf/cpm(i,k)*qci(i,k,2)
    1546. 

  1546.             qci(i,k,2) = 0.
    1547. 

  1547.           endif
    1548. 

  1548. !---------------------------------------------------------------
    1549. 

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

  1550. !        (T<-40C: C->I)
    1551. 

  1551. !---------------------------------------------------------------
    1552. 

  1552.           if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
    1553. 

  1553.             qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
    1554. 

  1554.             t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*qci(i,k,1)
    1555. 

  1555.             qci(i,k,1) = 0.
    1556. 

  1556.           endif
    1557. 

  1557. !---------------------------------------------------------------
    1558. 

  1558. ! pihtf: heterogeneous freezing of cloud water [HL A44]
    1559. 

  1559. !        (T0>T>-40C: C->I)
    1560. 

  1560. !---------------------------------------------------------------
    1561. 

  1561.           if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
    1562. 

  1562.             supcolt=min(supcol,50.)
    1563. 

  1563. !           pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.)                         &
    1564. 

  1564. !              *den(i,k,j)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
    1565. 

  1565.             pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.)                        &
    1566. 

  1566.             *den(i,k,j)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
    1567. 

  1567.             qci(i,k,2) = qci(i,k,2) + pfrzdtc
    1568. 

  1568.             t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtc
    1569. 

  1569.             qci(i,k,1) = qci(i,k,1)-pfrzdtc
    1570. 

  1570.           endif
    1571. 

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

  1572. ! psfrz: freezing of rain water [HL A20] [LFO 45]
    1573. 

  1573. !        (T<T0, R->S)
    1574. 

  1574. !---------------------------------------------------------------
    1575. 

  1575.           if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
    1576. 

  1576.             supcolt=min(supcol,50.)
    1577. 

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

  1578. !                 *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld,              &
    1579. 

  1579. !                 qrs(i,k,1))
    1580. 

  1580.             temp = rslope(i,k,1)
    1581. 

  1581.             temp = temp*temp*temp*temp*temp*temp*temp
    1582. 

  1582.             pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k,j)                  &
    1583. 

  1583.                   *(exp(pfrz2*supcolt)-1.)*temp*dtcld,                         &
    1584. 

  1584.                   qrs(i,k,1))
    1585. 

  1585.             qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
    1586. 

  1586.             t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtr
    1587. 

  1587.             qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
    1588. 

  1588.           endif
    1589. 

  1589.       enddo
    1590. 

  1590. !
    1591. 

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

  1592. !     rsloper: reverse of the slope parameter of the rain(m)
    1593. 

  1593. !     xka:    thermal conductivity of air(jm-1s-1k-1)
    1594. 

  1594. !     work1:  the thermodynamic term in the denominator associated with
    1595. 

  1595. !             heat conduction and vapor diffusion
    1596. 

  1596. !             (ry88, y93, h85)
    1597. 

  1597. !     work2: parameter associated with the ventilation effects(y93)
    1598. 

  1598. !
    1599. 

  1599.       do k = kts, kte
    1600. 

  1600.           if(qrs(i,k,1).le.qcrmin)then
    1601. 

  1601.             rslope(i,k,1) = rslopermax
    1602. 

  1602.             rslopeb(i,k,1) = rsloperbmax
    1603. 

  1603.             rslope2(i,k,1) = rsloper2max
    1604. 

  1604.             rslope3(i,k,1) = rsloper3max
    1605. 

  1605.           else
    1606. 

  1606. !           rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
    1607. 

  1607.             rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k,j))))))
    1608. 

  1608.             rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
    1609. 

  1609.             rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
    1610. 

  1610.             rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
    1611. 

  1611.           endif
    1612. 

  1612.           if(qrs(i,k,2).le.qcrmin)then
    1613. 

  1613.             rslope(i,k,2) = rslopesmax
    1614. 

  1614.             rslopeb(i,k,2) = rslopesbmax
    1615. 

  1615.             rslope2(i,k,2) = rslopes2max
    1616. 

  1616.             rslope3(i,k,2) = rslopes3max
    1617. 

  1617.           else
    1618. 

  1618. !            rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
    1619. 

  1619.             rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2))   &  
    1620. 

  1620.                           *(den(i,k,j))))))
    1621. 

  1621.             rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
    1622. 

  1622.             rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
    1623. 

  1623.             rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
    1624. 

  1624.           endif
    1625. 

  1625.       enddo
    1626. 

  1626. !
    1627. 

  1627.       do k = kts, kte
    1628. 

  1628. !         work1(i,k,1) = diffac(xl(i,k),p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,1))
    1629. 

  1629.           work1(i,k,1) = ((((den(i,k,j))*(xl(i,k))*(xl(i,k)))*((t(i,k,j))+120.)    &
    1630. 

  1630.                        *(den(i,k,j)))/(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))&
    1631. 

  1631.                        *(den(i,k,j))*(rv*(t(i,k,j))*(t(i,k,j)))))                    &
    1632. 

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

  1633. !         work1(i,k,2) = diffac(xls,p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,2))
    1634. 

  1634.           work1(i,k,2) = ((((den(i,k,j))*(xls)*(xls))*((t(i,k,j))+120.)*(den(i,k,j)))&
    1635. 

  1635.                        /(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))*(den(i,k,j)) &
    1636. 

  1636.                        *(rv*(t(i,k,j))*(t(i,k,j))))                                &
    1637. 

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

  1638. !         work2(i,k) = venfac(p(i,k,j),t(i,k,j),den(i,k,j))
    1639. 

  1639.           work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k,j))*sqrt(t(i,k,j)))) &
    1640. 

  1640.                      *p(i,k,j))/(((t(i,k,j))+120.)*den(i,k,j)*(8.794e-5              &
    1641. 

  1641.                      *exp(log(t(i,k,j))*(1.81))))))*sqrt(sqrt(den0/(den(i,k,j))))) &
    1642. 

  1642.                       /sqrt((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))                 &
    1643. 

  1643.                       /(((t(i,k,j))+120.)*den(i,k,j)))
    1644. 

  1644.       enddo 
    1645. 

  1645. !
    1646. 

  1646. !===============================================================
    1647. 

  1647. !
    1648. 

  1648. ! warm rain processes
    1649. 

  1649. !
    1650. 

  1650. ! - follows the processes in RH83 and LFO except for autoconcersion
    1651. 

  1651. !
    1652. 

  1652. !===============================================================
    1653. 

  1653. !
    1654. 

  1654.       do k = kts, kte
    1655. 

  1655.           supsat = max(q(i,k,j),qmin)-qs(i,k,1)
    1656. 

  1656.           satdt = supsat/dtcld
    1657. 

  1657. !---------------------------------------------------------------
    1658. 

  1658. ! praut: auto conversion rate from cloud to rain [HDC 16]
    1659. 

  1659. !        (C->R)
    1660. 

  1660. !---------------------------------------------------------------
    1661. 

  1661.           if(qci(i,k,1).gt.qc0) then
    1662. 

  1662.             praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
    1663. 

  1663.             praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
    1664. 

  1664.           endif
    1665. 

  1665. !---------------------------------------------------------------
    1666. 

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

  1667. !        (C->R)
    1668. 

  1668. !---------------------------------------------------------------
    1669. 

  1669.           if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
    1670. 

  1670.             pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1)               & 
    1671. 

  1671.                          *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
    1672. 

  1672.           endif
    1673. 

  1673. !---------------------------------------------------------------
    1674. 

  1674. ! prevp: evaporation/condensation rate of rain [HDC 14]
    1675. 

  1675. !        (V->R or R->V)
    1676. 

  1676. !---------------------------------------------------------------
    1677. 

  1677.           if(qrs(i,k,1).gt.0.) then
    1678. 

  1678.             coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
    1679. 

  1679.             prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1)                 &
    1680. 

  1680.                          +precr2*work2(i,k)*coeres)/work1(i,k,1)
    1681. 

  1681.             if(prevp(i,k).lt.0.) then
    1682. 

  1682.               prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
    1683. 

  1683.               prevp(i,k) = max(prevp(i,k),satdt/2)
    1684. 

  1684.             else
    1685. 

  1685.               prevp(i,k) = min(prevp(i,k),satdt/2)
    1686. 

  1686.             endif
    1687. 

  1687.           endif
    1688. 

  1688.       enddo
    1689. 

  1689. !
    1690. 

  1690. !===============================================================
    1691. 

  1691. !
    1692. 

  1692. ! cold rain processes
    1693. 

  1693. !
    1694. 

  1694. ! - follows the revised ice microphysics processes in HDC
    1695. 

  1695. ! - the processes same as in RH83 and RH84  and LFO behave
    1696. 

  1696. !   following ice crystal hapits defined in HDC, inclduing
    1697. 

  1697. !   intercept parameter for snow (n0s), ice crystal number
    1698. 

  1698. !   concentration (ni), ice nuclei number concentration
    1699. 

  1699. !   (n0i), ice diameter (d)
    1700. 

  1700. !
    1701. 

  1701. !===============================================================
    1702. 

  1702. !
    1703. 

  1703.       rdtcld = 1./dtcld
    1704. 

  1704.       do k = kts, kte
    1705. 

  1705.           supcol = t0c-t(i,k,j)
    1706. 

  1706.           supsat = max(q(i,k,j),qmin)-qs(i,k,2)
    1707. 

  1707.           satdt = supsat/dtcld
    1708. 

  1708.           ifsat = 0
    1709. 

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

  1710. ! Ni: ice crystal number concentraiton   [HDC 5c]
    1711. 

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

  1712. !         xni(i,k) = min(max(5.38e7*(den(i,k,j)                                  &
    1713. 

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

  1714.           temp = (den(i,k,j)*max(qci(i,k,2),qmin))
    1715. 

  1715.           temp = sqrt(sqrt(temp*temp*temp))
    1716. 

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

  1717.           eacrs = exp(0.07*(-supcol))
    1718. 

  1718. !
    1719. 

  1719.           if(supcol.gt.0) then
    1720. 

  1720.             if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
    1721. 

  1721.               xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
    1722. 

  1722.               diameter  = min(dicon * sqrt(xmi),dimax)
    1723. 

  1723.               vt2i = 1.49e4*diameter**1.31
    1724. 

  1724.               vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
    1725. 

  1725. !-------------------------------------------------------------
    1726. 

  1726. ! psaci: Accretion of cloud ice by rain [HDC 10]
    1727. 

  1727. !        (T<T0: I->S)
    1728. 

  1728. !-------------------------------------------------------------
    1729. 

  1729.               acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2)            &
    1730. 

  1730.                       +diameter**2*rslope(i,k,2)
    1731. 

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

  1732.                            *abs(vt2s-vt2i)*acrfac/4.
    1733. 

  1733.             endif
    1734. 

  1734.           endif
    1735. 

  1735. !-------------------------------------------------------------
    1736. 

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

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

  1738. !-------------------------------------------------------------
    1739. 

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

  1740.             psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)                  &
    1741. 

  1741.                            *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k)              &
    1742. 

  1742. !                          ,qci(i,k,1)/dtcld)
    1743. 

  1743.                            ,qci(i,k,1)*rdtcld)
    1744. 

  1744.           endif
    1745. 

  1745.           if(supcol .gt. 0) then
    1746. 

  1746. !-------------------------------------------------------------
    1747. 

  1747. ! pidep: Deposition/Sublimation rate of ice [HDC 9]
    1748. 

  1748. !       (T<T0: V->I or I->V)
    1749. 

  1749. !-------------------------------------------------------------
    1750. 

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

  1751.               xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
    1752. 

  1752.               diameter = dicon * sqrt(xmi)
    1753. 

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

  1754.               supice = satdt-prevp(i,k)
    1755. 

  1755.               if(pidep(i,k).lt.0.) then
    1756. 

  1756. !               pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
    1757. 

  1757. !               pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
    1758. 

  1758.                 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
    1759. 

  1759.                 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
    1760. 

  1760.               else
    1761. 

  1761. !               pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
    1762. 

  1762.                 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
    1763. 

  1763.               endif
    1764. 

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

  1765.             endif
    1766. 

  1766. !-------------------------------------------------------------
    1767. 

  1767. ! psdep: deposition/sublimation rate of snow [HDC 14]
    1768. 

  1768. !        (V->S or S->V)
    1769. 

  1769. !-------------------------------------------------------------
    1770. 

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

  1771.               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    1772. 

  1772.               psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)                          &
    1773. 

  1773.                            *(precs1*rslope2(i,k,2)+precs2                      &
    1774. 

  1774.                            *work2(i,k)*coeres)/work1(i,k,2)
    1775. 

  1775.               supice = satdt-prevp(i,k)-pidep(i,k)
    1776. 

  1776.               if(psdep(i,k).lt.0.) then
    1777. 

  1777. !               psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
    1778. 

  1778. !               psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
    1779. 

  1779.                 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
    1780. 

  1780.                 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
    1781. 

  1781.               else
    1782. 

  1782. !             psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
    1783. 

  1783.                 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
    1784. 

  1784.               endif
    1785. 

  1785.               if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt))          &
    1786. 

  1786.                 ifsat = 1
    1787. 

  1787.             endif
    1788. 

  1788. !-------------------------------------------------------------
    1789. 

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

  1790. !       (T<T0: V->I)
    1791. 

  1791. !-------------------------------------------------------------
    1792. 

  1792.             if(supsat.gt.0.and.ifsat.ne.1) then
    1793. 

  1793.               supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
    1794. 

  1794.               xni0 = 1.e3*exp(0.1*supcol)
    1795. 

  1795.               roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
    1796. 

  1796.               pigen(i,k) = max(0.,(roqi0/den(i,k,j)-max(qci(i,k,2),0.))          &
    1797. 

  1797. !                        /dtcld)
    1798. 

  1798.                          *rdtcld)
    1799. 

  1799.               pigen(i,k) = min(min(pigen(i,k),satdt),supice)
    1800. 

  1800.             endif
    1801. 

  1801. !
    1802. 

  1802. !-------------------------------------------------------------
    1803. 

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

  1804. !       (T<T0: I->S)
    1805. 

  1805. !-------------------------------------------------------------
    1806. 

  1806.             if(qci(i,k,2).gt.0.) then
    1807. 

  1807.               qimax = roqimax/den(i,k,j)
    1808. 

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

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

  1810.             endif
    1811. 

  1811.           endif
    1812. 

  1812. !-------------------------------------------------------------
    1813. 

  1813. ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
    1814. 

  1814. !       (T>T0: S->V)
    1815. 

  1815. !-------------------------------------------------------------
    1816. 

  1816.           if(supcol.lt.0.) then
    1817. 

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

  1818.               psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
    1819. 

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

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

  1821.           endif
    1822. 

  1822.       enddo
    1823. 

  1823. !
    1824. 

  1824. !
    1825. 

  1825. !----------------------------------------------------------------
    1826. 

  1826. !     check mass conservation of generation terms and feedback to the
    1827. 

  1827. !     large scale
    1828. 

  1828. !
    1829. 

  1829.       do k = kts, kte
    1830. 

  1830.           if(t(i,k,j).le.t0c) then
    1831. 

  1831. !
    1832. 

  1832. !     cloud water
    1833. 

  1833. !
    1834. 

  1834.             value = max(qmin,qci(i,k,1))
    1835. 

  1835.             source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
    1836. 

  1836.             if (source.gt.value) then
    1837. 

  1837.               factor = value/source
    1838. 

  1838.               praut(i,k) = praut(i,k)*factor
    1839. 

  1839.               pracw(i,k) = pracw(i,k)*factor
    1840. 

  1840.               psacw(i,k) = psacw(i,k)*factor
    1841. 

  1841.             endif
    1842. 

  1842. !
    1843. 

  1843. !     cloud ice
    1844. 

  1844. !
    1845. 

  1845.             value = max(qmin,qci(i,k,2))
    1846. 

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

  1847.             if (source.gt.value) then
    1848. 

  1848.               factor = value/source
    1849. 

  1849.               psaut(i,k) = psaut(i,k)*factor
    1850. 

  1850.               psaci(i,k) = psaci(i,k)*factor
    1851. 

  1851.               pigen(i,k) = pigen(i,k)*factor
    1852. 

  1852.               pidep(i,k) = pidep(i,k)*factor
    1853. 

  1853.             endif
    1854. 

  1854. !
    1855. 

  1855. !     rain
    1856. 

  1856. !
    1857. 

  1857. !
    1858. 

  1858.             value = max(qmin,qrs(i,k,1))
    1859. 

  1859.             source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
    1860. 

  1860.             if (source.gt.value) then
    1861. 

  1861.               factor = value/source
    1862. 

  1862.               praut(i,k) = praut(i,k)*factor
    1863. 

  1863.               pracw(i,k) = pracw(i,k)*factor
    1864. 

  1864.               prevp(i,k) = prevp(i,k)*factor
    1865. 

  1865.             endif
    1866. 

  1866. !
    1867. 

  1867. !    snow
    1868. 

  1868. !
    1869. 

  1869.             value = max(qmin,qrs(i,k,2))
    1870. 

  1870.             source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld  
    1871. 

  1871.             if (source.gt.value) then
    1872. 

  1872.               factor = value/source
    1873. 

  1873.               psdep(i,k) = psdep(i,k)*factor
    1874. 

  1874.               psaut(i,k) = psaut(i,k)*factor
    1875. 

  1875.               psaci(i,k) = psaci(i,k)*factor
    1876. 

  1876.               psacw(i,k) = psacw(i,k)*factor
    1877. 

  1877.             endif
    1878. 

  1878. !
    1879. 

  1879.             work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
    1880. 

  1880. !     update
    1881. 

  1881.             q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
    1882. 

  1882.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
    1883. 

  1883.                         +psacw(i,k))*dtcld,0.)
    1884. 

  1884.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
    1885. 

  1885.                         +prevp(i,k))*dtcld,0.)
    1886. 

  1886.             qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k)                 &
    1887. 

  1887.                         -pigen(i,k)-pidep(i,k))*dtcld,0.)
    1888. 

  1888.             qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)                 &
    1889. 

  1889.                         +psaci(i,k)+psacw(i,k))*dtcld,0.)
    1890. 

  1890.             xlf = xls-xl(i,k)
    1891. 

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

  1892.                       -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
    1893. 

  1893.             t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
    1894. 

  1894.           else
    1895. 

  1895. !
    1896. 

  1896. !     cloud water
    1897. 

  1897. !
    1898. 

  1898.             value = max(qmin,qci(i,k,1))
    1899. 

  1899.             source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
    1900. 

  1900.             if (source.gt.value) then
    1901. 

  1901.               factor = value/source
    1902. 

  1902.               praut(i,k) = praut(i,k)*factor
    1903. 

  1903.               pracw(i,k) = pracw(i,k)*factor
    1904. 

  1904.               psacw(i,k) = psacw(i,k)*factor
    1905. 

  1905.             endif
    1906. 

  1906. !
    1907. 

  1907. !     rain
    1908. 

  1908. !
    1909. 

  1909.             value = max(qmin,qrs(i,k,1))
    1910. 

  1910.             source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
    1911. 

  1911.             if (source.gt.value) then
    1912. 

  1912.               factor = value/source
    1913. 

  1913.               praut(i,k) = praut(i,k)*factor
    1914. 

  1914.               pracw(i,k) = pracw(i,k)*factor
    1915. 

  1915.               prevp(i,k) = prevp(i,k)*factor
    1916. 

  1916.               psacw(i,k) = psacw(i,k)*factor
    1917. 

  1917.             endif  
    1918. 

  1918. !
    1919. 

  1919. !     snow
    1920. 

  1920. !
    1921. 

  1921.             value = max(qcrmin,qrs(i,k,2))
    1922. 

  1922.             source=(-psevp(i,k))*dtcld
    1923. 

  1923.             if (source.gt.value) then
    1924. 

  1924.               factor = value/source
    1925. 

  1925.               psevp(i,k) = psevp(i,k)*factor
    1926. 

  1926.             endif
    1927. 

  1927.             work2(i,k)=-(prevp(i,k)+psevp(i,k))
    1928. 

  1928. !     update
    1929. 

  1929.             q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
    1930. 

  1930.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
    1931. 

  1931.                         +psacw(i,k))*dtcld,0.)
    1932. 

  1932.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
    1933. 

  1933.                         +prevp(i,k) +psacw(i,k))*dtcld,0.)
    1934. 

  1934.             qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
    1935. 

  1935.             xlf = xls-xl(i,k)
    1936. 

  1936.             xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
    1937. 

  1937.             t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
    1938. 

  1938.           endif
    1939. 

  1939.       enddo
    1940. 

  1940. !
    1941. 

  1941. ! Inline expansion for fpvs
    1942. 

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

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

  1944.       hsub = xls
    1945. 

  1945.       hvap = xlv0
    1946. 

  1946.       cvap = cpv
    1947. 

  1947.       ttp=t0c+0.01
    1948. 

  1948.       dldt=cvap-cliq
    1949. 

  1949.       xa=-dldt/rv
    1950. 

  1950.       xb=xa+hvap/(rv*ttp)
    1951. 

  1951.       dldti=cvap-cice
    1952. 

  1952.       xai=-dldti/rv
    1953. 

  1953.       xbi=xai+hsub/(rv*ttp)
    1954. 

  1954.       do k = kts, kte
    1955. 

  1955.           tr=ttp/t(i,k,j)
    1956. 

  1956.           logtr = log(tr)
    1957. 

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

  1958.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
    1959. 

  1959.           qs(i,k,1) = max(qs(i,k,1),qmin)
    1960. 

  1960.       enddo
    1961. 

  1961. !
    1962. 

  1962. !----------------------------------------------------------------
    1963. 

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

  1964. !     if there exists additional water vapor condensated/if
    1965. 

  1965. !     evaporation of cloud water is not enough to remove subsaturation
    1966. 

  1966. !
    1967. 

  1967.       do k = kts, kte
    1968. 

  1968. !         work1(i,k,1) = conden(t(i,k,j),q(i,k,j),qs(i,k,1),xl(i,k),cpm(i,k))
    1969. 

  1969.           work1(i,k,1) = ((max(q(i,k,j),qmin)-(qs(i,k,1)))/(1.+(xl(i,k))         &   
    1970. 

  1970.                         *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))                 &
    1971. 

  1971.                         /((t(i,k,j))*(t(i,k,j)))))
    1972. 

  1972.           work2(i,k) = qci(i,k,1)+work1(i,k,1)
    1973. 

  1973.           pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k,j),0.)/dtcld)
    1974. 

  1974.           if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.)                          &
    1975. 

  1975.             pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
    1976. 

  1976.           q(i,k,j) = q(i,k,j)-pcond(i,k)*dtcld
    1977. 

  1977.           qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
    1978. 

  1978.           t(i,k,j) = t(i,k,j)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
    1979. 

  1979.       enddo
    1980. 

  1980. !
    1981. 

  1981. !
    1982. 

  1982. !----------------------------------------------------------------
    1983. 

  1983. !     padding for small values
    1984. 

  1984. !
    1985. 

  1985.       do k = kts, kte
    1986. 

  1986.           if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
    1987. 

  1987.           if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
    1988. 

  1988.       enddo
    1989. 

  1989.       enddo                  ! big loops
    1990. 

  1990. 

  1991.          DO K=kts,kte
    1992. 

  1992.             th(i,k,j)=t(i,k,j)/pii(i,k,j)
    1993. 

  1993.             qc(i,k,j) = qci(i,k,1)
    1994. 

  1994.             qi(i,k,j) = qci(i,k,2)
    1995. 

  1995.             qr(i,k,j) = qrs(i,k,1)
    1996. 

  1996.             qqs(i,k,j) = qrs(i,k,2)
    1997. 

  1997.          ENDDO
    1998. 

  1998. 

  1999.       ENDDO     !  i loop
    2000. 

  2000.       enddo     !  j loop
    2001. 

  2001. !$acc end region
    2002. 

  2002. 

  2003.     ENDIF
    2004. 

  2004. 

  2005.   END SUBROUTINE wsm52d
    2006. 

  2006. 

  2007. 

  2008. #else
    2009. 

  2009. 

  2010. !===================================================================
    2011. 

  2011. !
    2012. 

  2012.   SUBROUTINE wsm52D(t, q, qci, qrs, den, p, delz                  &
    2013. 

  2013.                    ,delt,g, cpd, cpv, rd, rv, t0c                 &
    2014. 

  2014.                    ,ep1, ep2, qmin                                &
    2015. 

  2015.                    ,XLS, XLV0, XLF0, den0, denr                   &
    2016. 

  2016.                    ,cliq,cice,psat                                &
    2017. 

  2017.                    ,lat                                           &
    2018. 

  2018.                    ,rain,rainncv                                  &
    2019. 

  2019.                    ,sr                                            &
    2020. 

  2020.                    ,ids,ide, jds,jde, kds,kde                     &
    2021. 

  2021.                    ,ims,ime, jms,jme, kms,kme                     &
    2022. 

  2022.                    ,its,ite, jts,jte, kts,kte                     &
    2023. 

  2023.                    ,snow,snowncv                                  &
    2024. 

  2024.                                                                   )
    2025. 

  2025. !-------------------------------------------------------------------
    2026. 

  2026.   IMPLICIT NONE
    2027. 

  2027. !-------------------------------------------------------------------
    2028. 

  2028.   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
    2029. 

  2029.                                       ims,ime, jms,jme, kms,kme , &
    2030. 

  2030.                                       its,ite, jts,jte, kts,kte,  &
    2031. 

  2031.                                       lat
    2032. 

  2032.   REAL, DIMENSION( its:ite , kts:kte ),                           &
    2033. 

  2033.         INTENT(INOUT) ::                                          &
    2034. 

  2034.                                                                t
    2035. 

  2035.   REAL, DIMENSION( its:ite , kts:kte, 2 ),                        &
    2036. 

  2036.         INTENT(INOUT) ::                                          &
    2037. 

  2037.                                                              qci, &
    2038. 

  2038.                                                              qrs
    2039. 

  2039.   REAL, DIMENSION( ims:ime , kms:kme ),                           &
    2040. 

  2040.         INTENT(INOUT) ::                                          &
    2041. 

  2041.                                                                q
    2042. 

  2042.   REAL, DIMENSION( ims:ime , kms:kme ),                           &
    2043. 

  2043.         INTENT(IN   ) ::                                          &
    2044. 

  2044.                                                              den, &
    2045. 

  2045.                                                                p, &
    2046. 

  2046.                                                             delz
    2047. 

  2047.   REAL, INTENT(IN   ) ::                                    delt, &
    2048. 

  2048.                                                                g, &
    2049. 

  2049.                                                              cpd, &
    2050. 

  2050.                                                              cpv, &
    2051. 

  2051.                                                              t0c, &
    2052. 

  2052.                                                             den0, &
    2053. 

  2053.                                                               rd, &
    2054. 

  2054.                                                               rv, &
    2055. 

  2055.                                                              ep1, &
    2056. 

  2056.                                                              ep2, &
    2057. 

  2057.                                                             qmin, &
    2058. 

  2058.                                                              XLS, &
    2059. 

  2059.                                                             XLV0, &
    2060. 

  2060.                                                             XLF0, &
    2061. 

  2061.                                                             cliq, &
    2062. 

  2062.                                                             cice, &
    2063. 

  2063.                                                             psat, &
    2064. 

  2064.                                                             denr
    2065. 

  2065.   REAL, DIMENSION( ims:ime ),                                     &
    2066. 

  2066.         INTENT(INOUT) ::                                    rain, &
    2067. 

  2067.                                                          rainncv, &
    2068. 

  2068.                                                               sr
    2069. 

  2069. 

  2070.   REAL, DIMENSION( ims:ime, jms:jme ),     OPTIONAL,              &
    2071. 

  2071.         INTENT(INOUT) ::                                    snow, &
    2072. 

  2072.                                                          snowncv
    2073. 

  2073. 

  2074. ! LOCAL VAR
    2075. 

  2075.   REAL, DIMENSION( its:ite , kts:kte , 2) ::                      &
    2076. 

  2076.                                                               rh, &
    2077. 

  2077.                                                               qs, &
    2078. 

  2078.                                                           rslope, &
    2079. 

  2079.                                                          rslope2, &
    2080. 

  2080.                                                          rslope3, &
    2081. 

  2081.                                                          rslopeb, &
    2082. 

  2082.                                                             falk, &
    2083. 

  2083.                                                             fall, &
    2084. 

  2084.                                                            work1
    2085. 

  2085. 

  2086.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    2087. 

  2087.                                                            falkc, &
    2088. 

  2088.                                                            fallc, &
    2089. 

  2089.                                                               xl, &
    2090. 

  2090.                                                              cpm, &
    2091. 

  2091.                                                           denfac, &
    2092. 

  2092.                                                              xni, &
    2093. 

  2093.                                                           n0sfac, &
    2094. 

  2094.                                                            work2, &
    2095. 

  2095.                                                           work1c, &
    2096. 

  2096.                                                           work2c
    2097. 

  2097. 

  2098.   REAL, DIMENSION( its:ite , kts:kte ) ::                         &
    2099. 

  2099.                                                            pigen, &
    2100. 

  2100.                                                            pidep, &
    2101. 

  2101.                                                            psdep, &
    2102. 

  2102.                                                            praut, &
    2103. 

  2103.                                                            psaut, &
    2104. 

  2104.                                                            prevp, &
    2105. 

  2105.                                                            psevp, &
    2106. 

  2106.                                                            pracw, &
    2107. 

  2107.                                                            psacw, &
    2108. 

  2108.                                                            psaci, &
    2109. 

  2109.                                                            pcond, &
    2110. 

  2110.                                                            psmlt
    2111. 

  2111.   INTEGER, DIMENSION( its:ite ) ::                                &
    2112. 

  2112.                                                            mstep, &
    2113. 

  2113.                                                            numdt
    2114. 

  2114.   REAL, DIMENSION(its:ite) ::                             rmstep
    2115. 

  2115.   REAL dtcldden, rdelz, rdtcld
    2116. 

  2116.   LOGICAL, DIMENSION( its:ite ) ::                        flgcld
    2117. 

  2117. 

  2118. #define WSM_NO_CONDITIONAL_IN_VECTOR
    2119. 

  2119. #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
    2120. 

  2120.   REAL, DIMENSION(its:ite) :: xal, xbl
    2121. 

  2121. #endif
    2122. 

  2122. 

  2123.   REAL  ::  pi,                                                   &
    2124. 

  2124.             cpmcal, xlcal, lamdar, lamdas, diffus,                &
    2125. 

  2125.             viscos, xka, venfac, conden, diffac,                  &
    2126. 

  2126.             x, y, z, a, b, c, d, e,                               &
    2127. 

  2127.             qdt, holdrr, holdrs, supcol, supcolt, pvt,            &
    2128. 

  2128.             coeres, supsat, dtcld, xmi, eacrs, satdt,             &
    2129. 

  2129.             vt2i,vt2s,acrfac,                                     &
    2130. 

  2130.             qimax, diameter, xni0, roqi0,                         &
    2131. 

  2131.             fallsum, fallsum_qsi, xlwork2, factor, source,        &
    2132. 

  2132.             value, xlf, pfrzdtc, pfrzdtr, supice,  holdc, holdci
    2133. 

  2133. ! variables for optimization
    2134. 

  2134.   REAL, DIMENSION( its:ite )           ::                    tvec1
    2135. 

  2135.   REAL ::                                                    temp 
    2136. 

  2136.   INTEGER :: i, j, k, mstepmax,                                   &
    2137. 

  2137.             iprt, latd, lond, loop, loops, ifsat, n
    2138. 

  2138. ! Temporaries used for inlining fpvs function
    2139. 

  2139.   REAL  :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
    2140. 

  2140.   REAL  :: logtr
    2141. 

  2141. !
    2142. 

  2142. !=================================================================
    2143. 

  2143. !   compute internal functions
    2144. 

  2144. !
    2145. 

  2145.       cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
    2146. 

  2146.       xlcal(x) = xlv0-xlv1*(x-t0c)
    2147. 

  2147. !----------------------------------------------------------------
    2148. 

  2148. !     size distributions: (x=mixing ratio, y=air density):
    2149. 

  2149. !     valid for mixing ratio > 1.e-9 kg/kg.
    2150. 

  2150. !
    2151. 

  2151. ! Optimizatin : A**B => exp(log(A)*(B))
    2152. 

  2152.       lamdar(x,y)=   sqrt(sqrt(pidn0r/(x*y)))      ! (pidn0r/(x*y))**.25
    2153. 

  2153.       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
    2154. 

  2154. !
    2155. 

  2155. !----------------------------------------------------------------
    2156. 

  2156. !     diffus: diffusion coefficient of the water vapor
    2157. 

  2157. !     viscos: kinematic viscosity(m2s-1)
    2158. 

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

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

  2160. !     xka(x,y) = 1.414e3*viscos(x,y)*y
    2161. 

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

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

  2163. !                    /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
    2164. 

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

  2165. !
    2166. 

  2166. !
    2167. 

  2167.       pi = 4. * atan(1.)
    2168. 

  2168. !
    2169. 

  2169. !----------------------------------------------------------------
    2170. 

  2170. !     paddint 0 for negative values generated by dynamics
    2171. 

  2171. !
    2172. 

  2172.       do k = kts, kte
    2173. 

  2173.         do i = its, ite
    2174. 

  2174.           qci(i,k,1) = max(qci(i,k,1),0.0)
    2175. 

  2175.           qrs(i,k,1) = max(qrs(i,k,1),0.0)
    2176. 

  2176.           qci(i,k,2) = max(qci(i,k,2),0.0)
    2177. 

  2177.           qrs(i,k,2) = max(qrs(i,k,2),0.0)
    2178. 

  2178.         enddo
    2179. 

  2179.       enddo
    2180. 

  2180. !
    2181. 

  2181. !----------------------------------------------------------------
    2182. 

  2182. !     latent heat for phase changes and heat capacity. neglect the
    2183. 

  2183. !     changes during microphysical process calculation
    2184. 

  2184. !     emanuel(1994)
    2185. 

  2185. !
    2186. 

  2186.       do k = kts, kte
    2187. 

  2187.         do i = its, ite
    2188. 

  2188.           cpm(i,k) = cpmcal(q(i,k))
    2189. 

  2189.           xl(i,k) = xlcal(t(i,k))
    2190. 

  2190.         enddo
    2191. 

  2191.       enddo
    2192. 

  2192. !
    2193. 

  2193. !----------------------------------------------------------------
    2194. 

  2194. !     compute the minor time steps.
    2195. 

  2195. !
    2196. 

  2196.       loops = max(nint(delt/dtcldcr),1)
    2197. 

  2197.       dtcld = delt/loops
    2198. 

  2198.       if(delt.le.dtcldcr) dtcld = delt
    2199. 

  2199. !
    2200. 

  2200.       do loop = 1,loops
    2201. 

  2201. !
    2202. 

  2202. !----------------------------------------------------------------
    2203. 

  2203. !     initialize the large scale variables
    2204. 

  2204. !
    2205. 

  2205.       do i = its, ite
    2206. 

  2206.         mstep(i) = 1
    2207. 

  2207.         flgcld(i) = .true.
    2208. 

  2208.       enddo
    2209. 

  2209. !
    2210. 

  2210. !     do k = kts, kte
    2211. 

  2211. !       do i = its, ite
    2212. 

  2212. !         denfac(i,k) = sqrt(den0/den(i,k))
    2213. 

  2213. !       enddo
    2214. 

  2214. !     enddo
    2215. 

  2215.       do k = kts, kte
    2216. 

  2216.         CALL VREC( tvec1(its), den(its,k), ite-its+1)
    2217. 

  2217.         do i = its, ite
    2218. 

  2218.           tvec1(i) = tvec1(i)*den0
    2219. 

  2219.         enddo
    2220. 

  2220.         CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
    2221. 

  2221.       enddo
    2222. 

  2222. !
    2223. 

  2223. ! Inline expansion for fpvs
    2224. 

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

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

  2226.       hsub = xls
    2227. 

  2227.       hvap = xlv0
    2228. 

  2228.       cvap = cpv
    2229. 

  2229.       ttp=t0c+0.01
    2230. 

  2230.       dldt=cvap-cliq
    2231. 

  2231.       xa=-dldt/rv
    2232. 

  2232.       xb=xa+hvap/(rv*ttp)
    2233. 

  2233.       dldti=cvap-cice
    2234. 

  2234.       xai=-dldti/rv
    2235. 

  2235.       xbi=xai+hsub/(rv*ttp)
    2236. 

  2236. 

  2237. ! this is for compilers where the conditional inhibits vectorization
    2238. 

  2238. #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
    2239. 

  2239.       do k = kts, kte
    2240. 

  2240.         do i = its, ite
    2241. 

  2241.           if(t(i,k).lt.ttp) then
    2242. 

  2242.             xal(i) = xai
    2243. 

  2243.             xbl(i) = xbi
    2244. 

  2244.           else
    2245. 

  2245.             xal(i) = xa
    2246. 

  2246.             xbl(i) = xb
    2247. 

  2247.           endif
    2248. 

  2248.         enddo
    2249. 

  2249.         do i = its, ite
    2250. 

  2250.           tr=ttp/t(i,k)
    2251. 

  2251.           logtr=log(tr)
    2252. 

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

  2253.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    2254. 

  2254.           qs(i,k,1) = max(qs(i,k,1),qmin)
    2255. 

  2255.           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
    2256. 

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

  2257.           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
    2258. 

  2258.           qs(i,k,2) = max(qs(i,k,2),qmin)
    2259. 

  2259.           rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
    2260. 

  2260.         enddo
    2261. 

  2261.       enddo
    2262. 

  2262. #else
    2263. 

  2263.       do k = kts, kte
    2264. 

  2264.         do i = its, ite
    2265. 

  2265.           tr=ttp/t(i,k)
    2266. 

  2266.           logtr=log(tr)
    2267. 

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

  2268.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    2269. 

  2269.           qs(i,k,1) = max(qs(i,k,1),qmin)
    2270. 

  2270.           rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
    2271. 

  2271.           if(t(i,k).lt.ttp) then
    2272. 

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

  2273.           else
    2274. 

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

  2275.           endif
    2276. 

  2276.           qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
    2277. 

  2277.           qs(i,k,2) = max(qs(i,k,2),qmin)
    2278. 

  2278.           rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
    2279. 

  2279.         enddo
    2280. 

  2280.       enddo
    2281. 

  2281. #endif
    2282. 

  2282. !
    2283. 

  2283. !----------------------------------------------------------------
    2284. 

  2284. !     initialize the variables for microphysical physics
    2285. 

  2285. !
    2286. 

  2286. !
    2287. 

  2287.       do k = kts, kte
    2288. 

  2288.         do i = its, ite
    2289. 

  2289.           prevp(i,k) = 0.
    2290. 

  2290.           psdep(i,k) = 0.
    2291. 

  2291.           praut(i,k) = 0.
    2292. 

  2292.           psaut(i,k) = 0.
    2293. 

  2293.           pracw(i,k) = 0.
    2294. 

  2294.           psaci(i,k) = 0.
    2295. 

  2295.           psacw(i,k) = 0.
    2296. 

  2296.           pigen(i,k) = 0.
    2297. 

  2297.           pidep(i,k) = 0.
    2298. 

  2298.           pcond(i,k) = 0.
    2299. 

  2299.           psmlt(i,k) = 0.
    2300. 

  2300.           psevp(i,k) = 0.
    2301. 

  2301.           falk(i,k,1) = 0.
    2302. 

  2302.           falk(i,k,2) = 0.
    2303. 

  2303.           fall(i,k,1) = 0.
    2304. 

  2304.           fall(i,k,2) = 0.
    2305. 

  2305.           fallc(i,k) = 0.
    2306. 

  2306.           falkc(i,k) = 0.
    2307. 

  2307.           xni(i,k) = 1.e3
    2308. 

  2308.         enddo
    2309. 

  2309.       enddo
    2310. 

  2310. !
    2311. 

  2311. !----------------------------------------------------------------
    2312. 

  2312. !     compute the fallout term:
    2313. 

  2313. !     first, vertical terminal velosity for minor loops
    2314. 

  2314. !
    2315. 

  2315.       do k = kts, kte
    2316. 

  2316.         do i = its, ite
    2317. 

  2317.           supcol = t0c-t(i,k)
    2318. 

  2318. !---------------------------------------------------------------
    2319. 

  2319. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
    2320. 

  2320. !---------------------------------------------------------------
    2321. 

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

  2322.           if(qrs(i,k,1).le.qcrmin)then
    2323. 

  2323.             rslope(i,k,1) = rslopermax
    2324. 

  2324.             rslopeb(i,k,1) = rsloperbmax
    2325. 

  2325.             rslope2(i,k,1) = rsloper2max
    2326. 

  2326.             rslope3(i,k,1) = rsloper3max
    2327. 

  2327.           else
    2328. 

  2328.             rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
    2329. 

  2329.             rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
    2330. 

  2330.             rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
    2331. 

  2331.             rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
    2332. 

  2332.           endif
    2333. 

  2333.           if(qrs(i,k,2).le.qcrmin)then
    2334. 

  2334.             rslope(i,k,2) = rslopesmax
    2335. 

  2335.             rslopeb(i,k,2) = rslopesbmax
    2336. 

  2336.             rslope2(i,k,2) = rslopes2max
    2337. 

  2337.             rslope3(i,k,2) = rslopes3max
    2338. 

  2338.           else
    2339. 

  2339.             rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
    2340. 

  2340.             rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
    2341. 

  2341.             rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
    2342. 

  2342.             rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
    2343. 

  2343.           endif
    2344. 

  2344. !-------------------------------------------------------------
    2345. 

  2345. ! Ni: ice crystal number concentraiton   [HDC 5c]
    2346. 

  2346. !-------------------------------------------------------------
    2347. 

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

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

  2349.           temp = (den(i,k)*max(qci(i,k,2),qmin))
    2350. 

  2350.           temp = sqrt(sqrt(temp*temp*temp))
    2351. 

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

  2352.         enddo
    2353. 

  2353.       enddo
    2354. 

  2354. !
    2355. 

  2355.       mstepmax = 1
    2356. 

  2356.       numdt = 1
    2357. 

  2357.       do k = kte, kts, -1
    2358. 

  2358.         do i = its, ite
    2359. 

  2359.           work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k)
    2360. 

  2360.           work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k)
    2361. 

  2361.           numdt(i) = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
    2362. 

  2362.           if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
    2363. 

  2363.         enddo
    2364. 

  2364.       enddo
    2365. 

  2365.       do i = its, ite
    2366. 

  2366.         if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
    2367. 

  2367.         rmstep(i) = 1./mstep(i)
    2368. 

  2368.       enddo
    2369. 

  2369. !
    2370. 

  2370.       do n = 1, mstepmax
    2371. 

  2371.         k = kte
    2372. 

  2372.         do i = its, ite
    2373. 

  2373.           if(n.le.mstep(i)) then
    2374. 

  2374. !             falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
    2375. 

  2375. !             falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
    2376. 

  2376.               falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)*rmstep(i)
    2377. 

  2377.               falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)*rmstep(i)
    2378. 

  2378.               fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
    2379. 

  2379.               fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
    2380. 

  2380. !             qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
    2381. 

  2381. !             qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k),0.)
    2382. 

  2382.               dtcldden = dtcld/den(i,k)
    2383. 

  2383.               qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
    2384. 

  2384.               qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
    2385. 

  2385.             endif
    2386. 

  2386.           enddo
    2387. 

  2387.         do k = kte-1, kts, -1
    2388. 

  2388.           do i = its, ite
    2389. 

  2389.             if(n.le.mstep(i)) then
    2390. 

  2390.               falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)*rmstep(i)
    2391. 

  2391.               falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)*rmstep(i)
    2392. 

  2392.               fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
    2393. 

  2393.               fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
    2394. 

  2394.               dtcldden = dtcld/den(i,k)
    2395. 

  2395.               rdelz = 1./delz(i,k)
    2396. 

  2396.               qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1)           &
    2397. 

  2397.                           *delz(i,k+1)*rdelz)*dtcldden,0.)
    2398. 

  2398.               qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2)           &
    2399. 

  2399.                           *delz(i,k+1)*rdelz)*dtcldden,0.)
    2400. 

  2400.             endif
    2401. 

  2401.           enddo
    2402. 

  2402.         enddo
    2403. 

  2403.         do k = kte, kts, -1
    2404. 

  2404.           do i = its, ite
    2405. 

  2405.             if(n.le.mstep(i)) then
    2406. 

  2406.               if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
    2407. 

  2407. !----------------------------------------------------------------
    2408. 

  2408. ! psmlt: melting of snow [HL A33] [RH83 A25]
    2409. 

  2409. !       (T>T0: S->R)
    2410. 

  2410. !----------------------------------------------------------------
    2411. 

  2411.                 xlf = xlf0
    2412. 

  2412. !               work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
    2413. 

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

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

  2415.                             *exp(log(t(i,k))*(1.81))/p(i,k))))                 &
    2416. 

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

  2417.                             *sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k))))        &
    2418. 

  2418.                             *sqrt(sqrt(den0/(den(i,k)))))
    2419. 

  2419.                 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    2420. 

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

  2421. !                           *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
    2422. 

  2422. !                           *work2(i,k)*coeres)
    2423. 

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

  2424.                             /((t(i,k))+120.)/(den(i,k)) )*(den(i,k)))          &
    2425. 

  2425.                             /xlf*(t0c-t(i,k))*pi/2.                            &
    2426. 

  2426.                             *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
    2427. 

  2427.                             *work2(i,k)*coeres)
    2428. 

  2428.                 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),                &
    2429. 

  2429.                             -qrs(i,k,2)/mstep(i)),0.)
    2430. 

  2430.                 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
    2431. 

  2431.                 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
    2432. 

  2432.                 t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
    2433. 

  2433.               endif
    2434. 

  2434.             endif
    2435. 

  2435.           enddo
    2436. 

  2436.         enddo
    2437. 

  2437.       enddo
    2438. 

  2438. !---------------------------------------------------------------
    2439. 

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

  2440. !---------------------------------------------------------------
    2441. 

  2441.       mstepmax = 1
    2442. 

  2442.       mstep = 1
    2443. 

  2443.       numdt = 1
    2444. 

  2444.       do k = kte, kts, -1
    2445. 

  2445.         do i = its, ite
    2446. 

  2446.           if(qci(i,k,2).le.0.) then
    2447. 

  2447.             work2c(i,k) = 0.
    2448. 

  2448.           else
    2449. 

  2449.             xmi = den(i,k)*qci(i,k,2)/xni(i,k)
    2450. 

  2450. !           diameter  = min(dicon * sqrt(xmi),dimax)
    2451. 

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

  2452.             work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
    2453. 

  2453.             work2c(i,k) = work1c(i,k)/delz(i,k)
    2454. 

  2454.           endif
    2455. 

  2455.           numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
    2456. 

  2456.           if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
    2457. 

  2457.         enddo
    2458. 

  2458.       enddo
    2459. 

  2459.       do i = its, ite
    2460. 

  2460.         if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
    2461. 

  2461.       enddo
    2462. 

  2462. !
    2463. 

  2463.       do n = 1, mstepmax
    2464. 

  2464.         k = kte
    2465. 

  2465.         do i = its, ite
    2466. 

  2466.           if(n.le.mstep(i)) then
    2467. 

  2467.             falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
    2468. 

  2468.             holdc = falkc(i,k)
    2469. 

  2469.             fallc(i,k) = fallc(i,k)+falkc(i,k)
    2470. 

  2470.             holdci = qci(i,k,2)
    2471. 

  2471.             qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
    2472. 

  2472.           endif
    2473. 

  2473.         enddo
    2474. 

  2474.         do k = kte-1, kts, -1
    2475. 

  2475.           do i = its, ite
    2476. 

  2476.             if(n.le.mstep(i)) then
    2477. 

  2477.               falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
    2478. 

  2478.               holdc = falkc(i,k)
    2479. 

  2479.               fallc(i,k) = fallc(i,k)+falkc(i,k)
    2480. 

  2480.               holdci = qci(i,k,2)
    2481. 

  2481.               qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1)             &
    2482. 

  2482.                           *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
    2483. 

  2483.             endif
    2484. 

  2484.           enddo
    2485. 

  2485.         enddo
    2486. 

  2486.       enddo
    2487. 

  2487. !
    2488. 

  2488. !
    2489. 

  2489. !----------------------------------------------------------------
    2490. 

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

  2491. !
    2492. 

  2492.       do i = its, ite
    2493. 

  2493.         fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
    2494. 

  2494.         fallsum_qsi = fall(i,1,2)+fallc(i,1)
    2495. 

  2495.         rainncv(i) = 0.
    2496. 

  2496.         if(fallsum.gt.0.) then
    2497. 

  2497.           rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000.
    2498. 

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

  2499.         endif
    2500. 

  2500.         IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
    2501. 

  2501.         snowncv(i,lat) = 0.
    2502. 

  2502.         if(fallsum_qsi.gt.0.) then
    2503. 

  2503.           snowncv(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
    2504. 

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

  2505.         endif
    2506. 

  2506.         ENDIF
    2507. 

  2507.         sr(i) = 0.
    2508. 

  2508.         if(fallsum.gt.0.)sr(i)=fallsum_qsi*delz(i,kts)/denr*dtcld*1000.        &    
    2509. 

  2509.         /(rainncv(i)+1.e-12)
    2510. 

  2510.       enddo
    2511. 

  2511. !
    2512. 

  2512. !---------------------------------------------------------------
    2513. 

  2513. ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
    2514. 

  2514. !       (T>T0: I->C)
    2515. 

  2515. !---------------------------------------------------------------
    2516. 

  2516.       do k = kts, kte
    2517. 

  2517.         do i = its, ite
    2518. 

  2518.           supcol = t0c-t(i,k)
    2519. 

  2519.           xlf = xls-xl(i,k)
    2520. 

  2520.           if(supcol.lt.0.) xlf = xlf0
    2521. 

  2521.           if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
    2522. 

  2522.             qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
    2523. 

  2523.             t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
    2524. 

  2524.             qci(i,k,2) = 0.
    2525. 

  2525.           endif
    2526. 

  2526. !---------------------------------------------------------------
    2527. 

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

  2528. !        (T<-40C: C->I)
    2529. 

  2529. !---------------------------------------------------------------
    2530. 

  2530.           if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
    2531. 

  2531.             qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
    2532. 

  2532.             t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
    2533. 

  2533.             qci(i,k,1) = 0.
    2534. 

  2534.           endif
    2535. 

  2535. !---------------------------------------------------------------
    2536. 

  2536. ! pihtf: heterogeneous freezing of cloud water [HL A44]
    2537. 

  2537. !        (T0>T>-40C: C->I)
    2538. 

  2538. !---------------------------------------------------------------
    2539. 

  2539.           if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
    2540. 

  2540.             supcolt=min(supcol,50.)
    2541. 

  2541. !           pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.)                         &
    2542. 

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

  2543.             pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.)                        &
    2544. 

  2544.             *den(i,k)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
    2545. 

  2545.             qci(i,k,2) = qci(i,k,2) + pfrzdtc
    2546. 

  2546.             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
    2547. 

  2547.             qci(i,k,1) = qci(i,k,1)-pfrzdtc
    2548. 

  2548.           endif
    2549. 

  2549. !---------------------------------------------------------------
    2550. 

  2550. ! psfrz: freezing of rain water [HL A20] [LFO 45]
    2551. 

  2551. !        (T<T0, R->S)
    2552. 

  2552. !---------------------------------------------------------------
    2553. 

  2553.           if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
    2554. 

  2554.             supcolt=min(supcol,50.)
    2555. 

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

  2556. !                 *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld,              &
    2557. 

  2557. !                 qrs(i,k,1))
    2558. 

  2558.             temp = rslope(i,k,1)
    2559. 

  2559.             temp = temp*temp*temp*temp*temp*temp*temp
    2560. 

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

  2561.                   *(exp(pfrz2*supcolt)-1.)*temp*dtcld,                         &
    2562. 

  2562.                   qrs(i,k,1))
    2563. 

  2563.             qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
    2564. 

  2564.             t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
    2565. 

  2565.             qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
    2566. 

  2566.           endif
    2567. 

  2567.         enddo
    2568. 

  2568.       enddo
    2569. 

  2569. !
    2570. 

  2570. !----------------------------------------------------------------
    2571. 

  2571. !     rsloper: reverse of the slope parameter of the rain(m)
    2572. 

  2572. !     xka:    thermal conductivity of air(jm-1s-1k-1)
    2573. 

  2573. !     work1:  the thermodynamic term in the denominator associated with
    2574. 

  2574. !             heat conduction and vapor diffusion
    2575. 

  2575. !             (ry88, y93, h85)
    2576. 

  2576. !     work2: parameter associated with the ventilation effects(y93)
    2577. 

  2577. !
    2578. 

  2578.       do k = kts, kte
    2579. 

  2579.         do i = its, ite
    2580. 

  2580.           if(qrs(i,k,1).le.qcrmin)then
    2581. 

  2581.             rslope(i,k,1) = rslopermax
    2582. 

  2582.             rslopeb(i,k,1) = rsloperbmax
    2583. 

  2583.             rslope2(i,k,1) = rsloper2max
    2584. 

  2584.             rslope3(i,k,1) = rsloper3max
    2585. 

  2585.           else
    2586. 

  2586. !           rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
    2587. 

  2587.             rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k))))))
    2588. 

  2588.             rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
    2589. 

  2589.             rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
    2590. 

  2590.             rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
    2591. 

  2591.           endif
    2592. 

  2592.           if(qrs(i,k,2).le.qcrmin)then
    2593. 

  2593.             rslope(i,k,2) = rslopesmax
    2594. 

  2594.             rslopeb(i,k,2) = rslopesbmax
    2595. 

  2595.             rslope2(i,k,2) = rslopes2max
    2596. 

  2596.             rslope3(i,k,2) = rslopes3max
    2597. 

  2597.           else
    2598. 

  2598. !            rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
    2599. 

  2599.             rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2))   &  
    2600. 

  2600.                           *(den(i,k))))))
    2601. 

  2601.             rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
    2602. 

  2602.             rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
    2603. 

  2603.             rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
    2604. 

  2604.           endif
    2605. 

  2605.         enddo
    2606. 

  2606.       enddo
    2607. 

  2607. !
    2608. 

  2608.       do k = kts, kte
    2609. 

  2609.         do i = its, ite
    2610. 

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

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

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

  2613.                        *(den(i,k))*(rv*(t(i,k))*(t(i,k)))))                    &
    2614. 

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

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

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

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

  2618.                        *(rv*(t(i,k))*(t(i,k))))                                &
    2619. 

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

  2620. !         work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
    2621. 

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

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

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

  2624.                       /sqrt((1.496e-6*((t(i,k))*sqrt(t(i,k))))                 &
    2625. 

  2625.                       /(((t(i,k))+120.)*den(i,k)))
    2626. 

  2626.         enddo 
    2627. 

  2627.       enddo 
    2628. 

  2628. !
    2629. 

  2629. !===============================================================
    2630. 

  2630. !
    2631. 

  2631. ! warm rain processes
    2632. 

  2632. !
    2633. 

  2633. ! - follows the processes in RH83 and LFO except for autoconcersion
    2634. 

  2634. !
    2635. 

  2635. !===============================================================
    2636. 

  2636. !
    2637. 

  2637.       do k = kts, kte
    2638. 

  2638.         do i = its, ite
    2639. 

  2639.           supsat = max(q(i,k),qmin)-qs(i,k,1)
    2640. 

  2640.           satdt = supsat/dtcld
    2641. 

  2641. !---------------------------------------------------------------
    2642. 

  2642. ! praut: auto conversion rate from cloud to rain [HDC 16]
    2643. 

  2643. !        (C->R)
    2644. 

  2644. !---------------------------------------------------------------
    2645. 

  2645.           if(qci(i,k,1).gt.qc0) then
    2646. 

  2646.             praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
    2647. 

  2647.             praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
    2648. 

  2648.           endif
    2649. 

  2649. !---------------------------------------------------------------
    2650. 

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

  2651. !        (C->R)
    2652. 

  2652. !---------------------------------------------------------------
    2653. 

  2653.           if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
    2654. 

  2654.             pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1)               & 
    2655. 

  2655.                          *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
    2656. 

  2656.           endif
    2657. 

  2657. !---------------------------------------------------------------
    2658. 

  2658. ! prevp: evaporation/condensation rate of rain [HDC 14]
    2659. 

  2659. !        (V->R or R->V)
    2660. 

  2660. !---------------------------------------------------------------
    2661. 

  2661.           if(qrs(i,k,1).gt.0.) then
    2662. 

  2662.             coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
    2663. 

  2663.             prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1)                 &
    2664. 

  2664.                          +precr2*work2(i,k)*coeres)/work1(i,k,1)
    2665. 

  2665.             if(prevp(i,k).lt.0.) then
    2666. 

  2666.               prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
    2667. 

  2667.               prevp(i,k) = max(prevp(i,k),satdt/2)
    2668. 

  2668.             else
    2669. 

  2669.               prevp(i,k) = min(prevp(i,k),satdt/2)
    2670. 

  2670.             endif
    2671. 

  2671.           endif
    2672. 

  2672.         enddo
    2673. 

  2673.       enddo
    2674. 

  2674. !
    2675. 

  2675. !===============================================================
    2676. 

  2676. !
    2677. 

  2677. ! cold rain processes
    2678. 

  2678. !
    2679. 

  2679. ! - follows the revised ice microphysics processes in HDC
    2680. 

  2680. ! - the processes same as in RH83 and RH84  and LFO behave
    2681. 

  2681. !   following ice crystal hapits defined in HDC, inclduing
    2682. 

  2682. !   intercept parameter for snow (n0s), ice crystal number
    2683. 

  2683. !   concentration (ni), ice nuclei number concentration
    2684. 

  2684. !   (n0i), ice diameter (d)
    2685. 

  2685. !
    2686. 

  2686. !===============================================================
    2687. 

  2687. !
    2688. 

  2688.       rdtcld = 1./dtcld
    2689. 

  2689.       do k = kts, kte
    2690. 

  2690.         do i = its, ite
    2691. 

  2691.           supcol = t0c-t(i,k)
    2692. 

  2692.           supsat = max(q(i,k),qmin)-qs(i,k,2)
    2693. 

  2693.           satdt = supsat/dtcld
    2694. 

  2694.           ifsat = 0
    2695. 

  2695. !-------------------------------------------------------------
    2696. 

  2696. ! Ni: ice crystal number concentraiton   [HDC 5c]
    2697. 

  2697. !-------------------------------------------------------------
    2698. 

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

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

  2700.           temp = (den(i,k)*max(qci(i,k,2),qmin))
    2701. 

  2701.           temp = sqrt(sqrt(temp*temp*temp))
    2702. 

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

  2703.           eacrs = exp(0.07*(-supcol))
    2704. 

  2704. !
    2705. 

  2705.           if(supcol.gt.0) then
    2706. 

  2706.             if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
    2707. 

  2707.               xmi = den(i,k)*qci(i,k,2)/xni(i,k)
    2708. 

  2708.               diameter  = min(dicon * sqrt(xmi),dimax)
    2709. 

  2709.               vt2i = 1.49e4*diameter**1.31
    2710. 

  2710.               vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
    2711. 

  2711. !-------------------------------------------------------------
    2712. 

  2712. ! psaci: Accretion of cloud ice by rain [HDC 10]
    2713. 

  2713. !        (T<T0: I->S)
    2714. 

  2714. !-------------------------------------------------------------
    2715. 

  2715.               acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2)            &
    2716. 

  2716.                       +diameter**2*rslope(i,k,2)
    2717. 

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

  2718.                            *abs(vt2s-vt2i)*acrfac/4.
    2719. 

  2719.             endif
    2720. 

  2720.           endif
    2721. 

  2721. !-------------------------------------------------------------
    2722. 

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

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

  2724. !-------------------------------------------------------------
    2725. 

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

  2726.             psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)                  &
    2727. 

  2727.                            *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k)              &
    2728. 

  2728. !                          ,qci(i,k,1)/dtcld)
    2729. 

  2729.                            ,qci(i,k,1)*rdtcld)
    2730. 

  2730.           endif
    2731. 

  2731.           if(supcol .gt. 0) then
    2732. 

  2732. !-------------------------------------------------------------
    2733. 

  2733. ! pidep: Deposition/Sublimation rate of ice [HDC 9]
    2734. 

  2734. !       (T<T0: V->I or I->V)
    2735. 

  2735. !-------------------------------------------------------------
    2736. 

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

  2737.               xmi = den(i,k)*qci(i,k,2)/xni(i,k)
    2738. 

  2738.               diameter = dicon * sqrt(xmi)
    2739. 

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

  2740.               supice = satdt-prevp(i,k)
    2741. 

  2741.               if(pidep(i,k).lt.0.) then
    2742. 

  2742. !               pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
    2743. 

  2743. !               pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
    2744. 

  2744.                 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
    2745. 

  2745.                 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
    2746. 

  2746.               else
    2747. 

  2747. !               pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
    2748. 

  2748.                 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
    2749. 

  2749.               endif
    2750. 

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

  2751.             endif
    2752. 

  2752. !-------------------------------------------------------------
    2753. 

  2753. ! psdep: deposition/sublimation rate of snow [HDC 14]
    2754. 

  2754. !        (V->S or S->V)
    2755. 

  2755. !-------------------------------------------------------------
    2756. 

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

  2757.               coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
    2758. 

  2758.               psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)                          &
    2759. 

  2759.                            *(precs1*rslope2(i,k,2)+precs2                      &
    2760. 

  2760.                            *work2(i,k)*coeres)/work1(i,k,2)
    2761. 

  2761.               supice = satdt-prevp(i,k)-pidep(i,k)
    2762. 

  2762.               if(psdep(i,k).lt.0.) then
    2763. 

  2763. !               psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
    2764. 

  2764. !               psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
    2765. 

  2765.                 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
    2766. 

  2766.                 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
    2767. 

  2767.               else
    2768. 

  2768. !             psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
    2769. 

  2769.                 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
    2770. 

  2770.               endif
    2771. 

  2771.               if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt))          &
    2772. 

  2772.                 ifsat = 1
    2773. 

  2773.             endif
    2774. 

  2774. !-------------------------------------------------------------
    2775. 

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

  2776. !       (T<T0: V->I)
    2777. 

  2777. !-------------------------------------------------------------
    2778. 

  2778.             if(supsat.gt.0.and.ifsat.ne.1) then
    2779. 

  2779.               supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
    2780. 

  2780.               xni0 = 1.e3*exp(0.1*supcol)
    2781. 

  2781.               roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
    2782. 

  2782.               pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))          &
    2783. 

  2783. !                        /dtcld)
    2784. 

  2784.                          *rdtcld)
    2785. 

  2785.               pigen(i,k) = min(min(pigen(i,k),satdt),supice)
    2786. 

  2786.             endif
    2787. 

  2787. !
    2788. 

  2788. !-------------------------------------------------------------
    2789. 

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

  2790. !       (T<T0: I->S)
    2791. 

  2791. !-------------------------------------------------------------
    2792. 

  2792.             if(qci(i,k,2).gt.0.) then
    2793. 

  2793.               qimax = roqimax/den(i,k)
    2794. 

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

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

  2796.             endif
    2797. 

  2797.           endif
    2798. 

  2798. !-------------------------------------------------------------
    2799. 

  2799. ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
    2800. 

  2800. !       (T>T0: S->V)
    2801. 

  2801. !-------------------------------------------------------------
    2802. 

  2802.           if(supcol.lt.0.) then
    2803. 

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

  2804.               psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
    2805. 

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

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

  2807.           endif
    2808. 

  2808.         enddo
    2809. 

  2809.       enddo
    2810. 

  2810. !
    2811. 

  2811. !
    2812. 

  2812. !----------------------------------------------------------------
    2813. 

  2813. !     check mass conservation of generation terms and feedback to the
    2814. 

  2814. !     large scale
    2815. 

  2815. !
    2816. 

  2816.       do k = kts, kte
    2817. 

  2817.         do i = its, ite
    2818. 

  2818.           if(t(i,k).le.t0c) then
    2819. 

  2819. !
    2820. 

  2820. !     cloud water
    2821. 

  2821. !
    2822. 

  2822.             value = max(qmin,qci(i,k,1))
    2823. 

  2823.             source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
    2824. 

  2824.             if (source.gt.value) then
    2825. 

  2825.               factor = value/source
    2826. 

  2826.               praut(i,k) = praut(i,k)*factor
    2827. 

  2827.               pracw(i,k) = pracw(i,k)*factor
    2828. 

  2828.               psacw(i,k) = psacw(i,k)*factor
    2829. 

  2829.             endif
    2830. 

  2830. !
    2831. 

  2831. !     cloud ice
    2832. 

  2832. !
    2833. 

  2833.             value = max(qmin,qci(i,k,2))
    2834. 

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

  2835.             if (source.gt.value) then
    2836. 

  2836.               factor = value/source
    2837. 

  2837.               psaut(i,k) = psaut(i,k)*factor
    2838. 

  2838.               psaci(i,k) = psaci(i,k)*factor
    2839. 

  2839.               pigen(i,k) = pigen(i,k)*factor
    2840. 

  2840.               pidep(i,k) = pidep(i,k)*factor
    2841. 

  2841.             endif
    2842. 

  2842. !
    2843. 

  2843. !     rain
    2844. 

  2844. !
    2845. 

  2845. !
    2846. 

  2846.             value = max(qmin,qrs(i,k,1))
    2847. 

  2847.             source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
    2848. 

  2848.             if (source.gt.value) then
    2849. 

  2849.               factor = value/source
    2850. 

  2850.               praut(i,k) = praut(i,k)*factor
    2851. 

  2851.               pracw(i,k) = pracw(i,k)*factor
    2852. 

  2852.               prevp(i,k) = prevp(i,k)*factor
    2853. 

  2853.             endif
    2854. 

  2854. !
    2855. 

  2855. !    snow
    2856. 

  2856. !
    2857. 

  2857.             value = max(qmin,qrs(i,k,2))
    2858. 

  2858.             source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld  
    2859. 

  2859.             if (source.gt.value) then
    2860. 

  2860.               factor = value/source
    2861. 

  2861.               psdep(i,k) = psdep(i,k)*factor
    2862. 

  2862.               psaut(i,k) = psaut(i,k)*factor
    2863. 

  2863.               psaci(i,k) = psaci(i,k)*factor
    2864. 

  2864.               psacw(i,k) = psacw(i,k)*factor
    2865. 

  2865.             endif
    2866. 

  2866. !
    2867. 

  2867.             work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
    2868. 

  2868. !     update
    2869. 

  2869.             q(i,k) = q(i,k)+work2(i,k)*dtcld
    2870. 

  2870.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
    2871. 

  2871.                         +psacw(i,k))*dtcld,0.)
    2872. 

  2872.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
    2873. 

  2873.                         +prevp(i,k))*dtcld,0.)
    2874. 

  2874.             qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k)                 &
    2875. 

  2875.                         -pigen(i,k)-pidep(i,k))*dtcld,0.)
    2876. 

  2876.             qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)                 &
    2877. 

  2877.                         +psaci(i,k)+psacw(i,k))*dtcld,0.)
    2878. 

  2878.             xlf = xls-xl(i,k)
    2879. 

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

  2880.                       -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
    2881. 

  2881.             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
    2882. 

  2882.           else
    2883. 

  2883. !
    2884. 

  2884. !     cloud water
    2885. 

  2885. !
    2886. 

  2886.             value = max(qmin,qci(i,k,1))
    2887. 

  2887.             source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
    2888. 

  2888.             if (source.gt.value) then
    2889. 

  2889.               factor = value/source
    2890. 

  2890.               praut(i,k) = praut(i,k)*factor
    2891. 

  2891.               pracw(i,k) = pracw(i,k)*factor
    2892. 

  2892.               psacw(i,k) = psacw(i,k)*factor
    2893. 

  2893.             endif
    2894. 

  2894. !
    2895. 

  2895. !     rain
    2896. 

  2896. !
    2897. 

  2897.             value = max(qmin,qrs(i,k,1))
    2898. 

  2898.             source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
    2899. 

  2899.             if (source.gt.value) then
    2900. 

  2900.               factor = value/source
    2901. 

  2901.               praut(i,k) = praut(i,k)*factor
    2902. 

  2902.               pracw(i,k) = pracw(i,k)*factor
    2903. 

  2903.               prevp(i,k) = prevp(i,k)*factor
    2904. 

  2904.               psacw(i,k) = psacw(i,k)*factor
    2905. 

  2905.             endif  
    2906. 

  2906. !
    2907. 

  2907. !     snow
    2908. 

  2908. !
    2909. 

  2909.             value = max(qcrmin,qrs(i,k,2))
    2910. 

  2910.             source=(-psevp(i,k))*dtcld
    2911. 

  2911.             if (source.gt.value) then
    2912. 

  2912.               factor = value/source
    2913. 

  2913.               psevp(i,k) = psevp(i,k)*factor
    2914. 

  2914.             endif
    2915. 

  2915.             work2(i,k)=-(prevp(i,k)+psevp(i,k))
    2916. 

  2916. !     update
    2917. 

  2917.             q(i,k) = q(i,k)+work2(i,k)*dtcld
    2918. 

  2918.             qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
    2919. 

  2919.                         +psacw(i,k))*dtcld,0.)
    2920. 

  2920.             qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
    2921. 

  2921.                         +prevp(i,k) +psacw(i,k))*dtcld,0.)
    2922. 

  2922.             qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
    2923. 

  2923.             xlf = xls-xl(i,k)
    2924. 

  2924.             xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
    2925. 

  2925.             t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
    2926. 

  2926.           endif
    2927. 

  2927.         enddo
    2928. 

  2928.       enddo
    2929. 

  2929. !
    2930. 

  2930. ! Inline expansion for fpvs
    2931. 

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

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

  2933.       hsub = xls
    2934. 

  2934.       hvap = xlv0
    2935. 

  2935.       cvap = cpv
    2936. 

  2936.       ttp=t0c+0.01
    2937. 

  2937.       dldt=cvap-cliq
    2938. 

  2938.       xa=-dldt/rv
    2939. 

  2939.       xb=xa+hvap/(rv*ttp)
    2940. 

  2940.       dldti=cvap-cice
    2941. 

  2941.       xai=-dldti/rv
    2942. 

  2942.       xbi=xai+hsub/(rv*ttp)
    2943. 

  2943.       do k = kts, kte
    2944. 

  2944.       do i = its, ite
    2945. 

  2945.           tr=ttp/t(i,k)
    2946. 

  2946.           logtr = log(tr)
    2947. 

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

  2948.           qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
    2949. 

  2949.           qs(i,k,1) = max(qs(i,k,1),qmin)
    2950. 

  2950.         enddo
    2951. 

  2951.       enddo
    2952. 

  2952. !
    2953. 

  2953. !----------------------------------------------------------------
    2954. 

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

  2955. !     if there exists additional water vapor condensated/if
    2956. 

  2956. !     evaporation of cloud water is not enough to remove subsaturation
    2957. 

  2957. !
    2958. 

  2958.       do k = kts, kte
    2959. 

  2959.         do i = its, ite
    2960. 

  2960. !         work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
    2961. 

  2961.           work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/(1.+(xl(i,k))         &   
    2962. 

  2962.                         *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))                 &
    2963. 

  2963.                         /((t(i,k))*(t(i,k)))))
    2964. 

  2964.           work2(i,k) = qci(i,k,1)+work1(i,k,1)
    2965. 

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

  2966.           if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.)                          &
    2967. 

  2967.             pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
    2968. 

  2968.           q(i,k) = q(i,k)-pcond(i,k)*dtcld
    2969. 

  2969.           qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
    2970. 

  2970.           t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
    2971. 

  2971.         enddo
    2972. 

  2972.       enddo
    2973. 

  2973. !
    2974. 

  2974. !
    2975. 

  2975. !----------------------------------------------------------------
    2976. 

  2976. !     padding for small values
    2977. 

  2977. !
    2978. 

  2978.       do k = kts, kte
    2979. 

  2979.         do i = its, ite
    2980. 

  2980.           if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
    2981. 

  2981.           if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
    2982. 

  2982.         enddo
    2983. 

  2983.       enddo
    2984. 

  2984.       enddo                  ! big loops
    2985. 

  2985.   END SUBROUTINE wsm52d
    2986. 

  2986. 

  2987. #endif
    2988. 

  2988. 

  2989. ! ...................................................................
    2990. 

  2990.       REAL FUNCTION rgmma(x)
    2991. 

  2991. !-------------------------------------------------------------------
    2992. 

  2992.   IMPLICIT NONE
    2993. 

  2993. !-------------------------------------------------------------------
    2994. 

  2994. !     rgmma function:  use infinite product form
    2995. 

  2995.       REAL :: euler
    2996. 

  2996.       PARAMETER (euler=0.577215664901532)
    2997. 

  2997.       REAL :: x, y
    2998. 

  2998.       INTEGER :: i
    2999. 

  2999.       if(x.eq.1.)then
    3000. 

  3000.         rgmma=0.
    3001. 

  3001.           else
    3002. 

  3002.         rgmma=x*exp(euler*x)
    3003. 

  3003.         do i=1,10000
    3004. 

  3004.           y=float(i)
    3005. 

  3005.           rgmma=rgmma*(1.000+x/y)*exp(-x/y)
    3006. 

  3006.         enddo
    3007. 

  3007.         rgmma=1./rgmma
    3008. 

  3008.       endif
    3009. 

  3009.       END FUNCTION rgmma
    3010. 

  3010. !
    3011. 

  3011. !--------------------------------------------------------------------------
    3012. 

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

  3013. !--------------------------------------------------------------------------
    3014. 

  3014.       IMPLICIT NONE
    3015. 

  3015. !--------------------------------------------------------------------------
    3016. 

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

  3017.            xai,xbi,ttp,tr
    3018. 

  3018.       INTEGER ice
    3019. 

  3019. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    3020. 

  3020.       ttp=t0c+0.01
    3021. 

  3021.       dldt=cvap-cliq
    3022. 

  3022.       xa=-dldt/rv
    3023. 

  3023.       xb=xa+hvap/(rv*ttp)
    3024. 

  3024.       dldti=cvap-cice
    3025. 

  3025.       xai=-dldti/rv
    3026. 

  3026.       xbi=xai+hsub/(rv*ttp)
    3027. 

  3027.       tr=ttp/t
    3028. 

  3028.       if(t.lt.ttp.and.ice.eq.1) then
    3029. 

  3029.         fpvs=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
    3030. 

  3030.       else
    3031. 

  3031.         fpvs=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
    3032. 

  3032.       endif
    3033. 

  3033. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    3034. 

  3034.       END FUNCTION fpvs
    3035. 

  3035. !-------------------------------------------------------------------
    3036. 

  3036.   SUBROUTINE wsm5init(den0,denr,dens,cl,cpv,allowed_to_read)
    3037. 

  3037. !-------------------------------------------------------------------
    3038. 

  3038.   IMPLICIT NONE
    3039. 

  3039. !-------------------------------------------------------------------
    3040. 

  3040. !.... constants which may not be tunable
    3041. 

  3041.    REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
    3042. 

  3042.    LOGICAL, INTENT(IN) :: allowed_to_read
    3043. 

  3043.    REAL :: pi
    3044. 

  3044. !
    3045. 

  3045.    pi = 4.*atan(1.)
    3046. 

  3046.    xlv1 = cl-cpv
    3047. 

  3047. !
    3048. 

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

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

  3050. !
    3051. 

  3051.    bvtr1 = 1.+bvtr
    3052. 

  3052.    bvtr2 = 2.5+.5*bvtr
    3053. 

  3053.    bvtr3 = 3.+bvtr
    3054. 

  3054.    bvtr4 = 4.+bvtr
    3055. 

  3055.    g1pbr = rgmma(bvtr1)
    3056. 

  3056.    g3pbr = rgmma(bvtr3)
    3057. 

  3057.    g4pbr = rgmma(bvtr4)            ! 17.837825
    3058. 

  3058.    g5pbro2 = rgmma(bvtr2)          ! 1.8273
    3059. 

  3059.    pvtr = avtr*g4pbr/6.
    3060. 

  3060.    eacrr = 1.0
    3061. 

  3061.    pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
    3062. 

  3062.    precr1 = 2.*pi*n0r*.78
    3063. 

  3063.    precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
    3064. 

  3064.    xmmax = (dimax/dicon)**2
    3065. 

  3065.    roqimax = 2.08e22*dimax**8
    3066. 

  3066. !
    3067. 

  3067.    bvts1 = 1.+bvts
    3068. 

  3068.    bvts2 = 2.5+.5*bvts
    3069. 

  3069.    bvts3 = 3.+bvts
    3070. 

  3070.    bvts4 = 4.+bvts
    3071. 

  3071.    g1pbs = rgmma(bvts1)    !.8875
    3072. 

  3072.    g3pbs = rgmma(bvts3)
    3073. 

  3073.    g4pbs = rgmma(bvts4)    ! 12.0786
    3074. 

  3074.    g5pbso2 = rgmma(bvts2)
    3075. 

  3075.    pvts = avts*g4pbs/6.
    3076. 

  3076.    pacrs = pi*n0s*avts*g3pbs*.25
    3077. 

  3077.    precs1 = 4.*n0s*.65
    3078. 

  3078.    precs2 = 4.*n0s*.44*avts**.5*g5pbso2
    3079. 

  3079.    pidn0r =  pi*denr*n0r
    3080. 

  3080.    pidn0s =  pi*dens*n0s
    3081. 

  3081.    pacrc = pi*n0s*avts*g3pbs*.25*eacrc
    3082. 

  3082. !
    3083. 

  3083.    rslopermax = 1./lamdarmax
    3084. 

  3084.    rslopesmax = 1./lamdasmax
    3085. 

  3085.    rsloperbmax = rslopermax ** bvtr
    3086. 

  3086.    rslopesbmax = rslopesmax ** bvts
    3087. 

  3087.    rsloper2max = rslopermax * rslopermax
    3088. 

  3088.    rslopes2max = rslopesmax * rslopesmax
    3089. 

  3089.    rsloper3max = rsloper2max * rslopermax
    3090. 

  3090.    rslopes3max = rslopes2max * rslopesmax
    3091. 

  3091. !
    3092. 

  3092.   END SUBROUTINE wsm5init
    3093. 

  3093. END MODULE module_mp_wsm5
    3094. 

  3094.