PageRenderTime 74ms CodeModel.GetById 35ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_mp_wdm5.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2015 lines | 1436 code | 0 blank | 579 comment | 101 complexity | a595cff214ff08031fdb6a3c0c1a78e6 MD5 | raw file
Possible License(s): AGPL-1.0

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

    

  1. #if ( RWORDSIZE == 4 )
    2. 

  2. #  define VREC vsrec
    3. 

  3. #  define VSQRT vssqrt
    4. 

  4. #else
    5. 

  5. #  define VREC vrec
    6. 

  6. #  define VSQRT vsqrt
    7. 

  7. #endif
    8. 

  8. !
    9. 

  9. !Including inline expansion statistical function 
    10. 

  10. MODULE module_mp_wdm5
    11. 

  11. !
    12. 

  12. !
    13. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  38. !
    39. 

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

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

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

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

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

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

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

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

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

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

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

  50.    REAL, SAVE ::                                      &
    51. 

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

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

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

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

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

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

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

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

  59.              rslopecmax,rslopec2max,rslopec3max,      &
    60. 

  60.              rslopermax,rslopesmax,rslopegmax,        &
    61. 

  61.              rsloperbmax,rslopesbmax,rslopegbmax,     &
    62. 

  62.              rsloper2max,rslopes2max,rslopeg2max,     &
    63. 

  63.              rsloper3max,rslopes3max,rslopeg3max
    64. 

  64. !
    65. 

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

  66. !
    67. 

  67. CONTAINS
    68. 

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

  69. !
    70. 

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

  71.                  ,nn, nc, nr                                       &
    72. 

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

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

  74.                  ,ep1, ep2, qmin                                   &
    75. 

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

  76.                  ,cliq,cice,psat                                   &
    77. 

  77.                  ,rain, rainncv                                    &
    78. 

  78.                  ,snow, snowncv                                    &
    79. 

  79.                  ,sr                                               &
    80. 

  80.                  ,itimestep                                        &
    81. 

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

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

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

  84.                                                                    )
    85. 

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

  86.   IMPLICIT NONE
    87. 

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

  88. !
    89. 

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

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

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

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

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

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

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

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

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

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

  99. !
    100. 

  100. !  WDM5 cloud scheme
    101. 

  101. !
    102. 

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

  103. !
    104. 

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

  105. !
    106. 

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

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

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

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

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

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

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

  113. !
    114. 

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

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

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

  117. !
    118. 

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

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

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

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

  122.         INTENT(INOUT) ::                                          &
    123. 

  123.                                                              th,  &
    124. 

  124.                                                               q,  &
    125. 

  125.                                                               qc, &
    126. 

  126.                                                               qi, &
    127. 

  127.                                                               qr, &
    128. 

  128.                                                               qs, &
    129. 

  129.                                                               nn, & 
    130. 

  130.                                                               nc, &
    131. 

  131.                                                               nr   
    132. 

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

  133.         INTENT(IN   ) ::                                          &
    134. 

  134.                                                              den, &
    135. 

  135.                                                              pii, &
    136. 

  136.                                                                p, &
    137. 

  137.                                                             delz
    138. 

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

  139.                                                                g, &
    140. 

  140.                                                               rd, &
    141. 

  141.                                                               rv, &
    142. 

  142.                                                              t0c, &
    143. 

  143.                                                             den0, &
    144. 

  144.                                                              cpd, &
    145. 

  145.                                                              cpv, &
    146. 

  146.                                                             ccn0, &
    147. 

  147.                                                              ep1, &
    148. 

  148.                                                              ep2, &
    149. 

  149.                                                             qmin, &
    150. 

  150.                                                              XLS, &
    151. 

  151.                                                             XLV0, &
    152. 

  152.                                                             XLF0, &
    153. 

  153.                                                             cliq, &
    154. 

  154.                                                             cice, &
    155. 

  155.                                                             psat, &
    156. 

  156.                                                             denr
    157. 

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

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

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

  160.                                                          rainncv, &
    161. 

  161.                                                               sr
    162. 

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

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

  164.                                                          snowncv
    165. 

  165. ! LOCAL VAR
    166. 

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

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

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

  169.   CHARACTER*256 :: emess
    170. 

  170.   INTEGER :: mkx_test
    171. 

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

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

  173. #ifndef RUN_ON_GPU
    174. 

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

  175.         DO j=jms,jme
    176. 

  176.            DO k=kms,kme
    177. 

  177.            DO i=ims,ime
    178. 

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

  179.            ENDDO
    180. 

  180.            ENDDO
    181. 

  181.         ENDDO
    182. 

  182.       ENDIF
    183. 

  183. !
    184. 

  184.       DO j=jts,jte
    185. 

  185.          DO k=kts,kte
    186. 

  186.          DO i=its,ite
    187. 

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

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

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

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

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

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

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

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

  195.          ENDDO
    196. 

  196.          ENDDO
    197. 

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

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

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

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

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

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

  203.                     ,ep1, ep2, qmin                               &
    204. 

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

  205.                     ,cliq,cice,psat                               &
    206. 

  206.                     ,j                                            &
    207. 

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

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

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

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

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

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

  213.                                                                   )
    214. 

  214.          DO K=kts,kte
    215. 

  215.          DO I=its,ite
    216. 

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

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

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

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

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

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

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

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

  224.          ENDDO
    225. 

  225.          ENDDO
    226. 

  226.       ENDDO
    227. 

  227. #else
    228. 

  228.       CALL get_wsm5_gpu_levels ( mkx_test )
    229. 

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

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

  231.                       mkx_test,' < ',kte
    232. 

  232.         CALL wrf_error_fatal(emess)
    233. 

  233.       ENDIF
    234. 

  234.       CALL wsm5_host (                                                         &
    235. 

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

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

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

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

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

  240.                    ,delt                                                       &
    241. 

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

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

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

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

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

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

  247.           )
    248. 

  248. #endif
    249. 

  249.   END SUBROUTINE wdm5
    250. 

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

  251. !
    252. 

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

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

  254.                    ,ep1, ep2, qmin                                &
    255. 

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

  256.                    ,cliq,cice,psat                                &
    257. 

  257.                    ,lat                                           &
    258. 

  258.                    ,rain,rainncv                                  &
    259. 

  259.                    ,sr                                            &
    260. 

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

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

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

  263.                    ,snow,snowncv                                  &
    264. 

  264.                                                                   )
    265. 

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

  266.   IMPLICIT NONE
    267. 

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

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

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

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

  271.                                       lat
    272. 

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

  273.         INTENT(INOUT) ::                                          &
    274. 

  274.                                                                t
    275. 

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

  276.         INTENT(INOUT) ::                                          &
    277. 

  277.                                                              qci, &
    278. 

  278.                                                              qrs 
    279. 

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

  280.         INTENT(INOUT) ::                                          &
    281. 

  281.                                                              ncr
    282. 

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

  283.         INTENT(INOUT) ::                                          &
    284. 

  284.                                                                q
    285. 

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

  286.         INTENT(IN   ) ::                                          &
    287. 

  287.                                                              den, &
    288. 

  288.                                                                p, &
    289. 

  289.                                                             delz
    290. 

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

  291.                                                                g, &
    292. 

  292.                                                              cpd, &
    293. 

  293.                                                              cpv, &
    294. 

  294.                                                             ccn0, &
    295. 

  295.                                                              t0c, &
    296. 

  296.                                                             den0, &
    297. 

  297.                                                               rd, &
    298. 

  298.                                                               rv, &
    299. 

  299.                                                              ep1, &
    300. 

  300.                                                              ep2, &
    301. 

  301.                                                             qmin, &
    302. 

  302.                                                              XLS, &
    303. 

  303.                                                             XLV0, &
    304. 

  304.                                                             XLF0, &
    305. 

  305.                                                             cliq, &
    306. 

  306.                                                             cice, &
    307. 

  307.                                                             psat, &
    308. 

  308.                                                             denr
    309. 

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

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

  311.                                                          rainncv, &
    312. 

  312.                                                               sr
    313. 

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

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

  315.                                                          snowncv
    316. 

  316. ! LOCAL VAR
    317. 

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

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

  319.         falk, fall, work1, qrs_tmp
    320. 

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

  321.         rslopec, rslopec2,rslopec3                           
    322. 

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

  323.         avedia 
    324. 

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

  325.         workn,falln,falkn
    326. 

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

  327.         works
    328. 

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

  329.         den_tmp, delz_tmp, ncr_tmp
    330. 

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

  331.               falkc, work1c, work2c, fallc
    332. 

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

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

  334.         pigen, pidep, pcond,                                      &
    335. 

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

  336.         n0sfac, denqrs2, denqci
    337. 

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

  338.         delqrs2, delqi
    339. 

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

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

  341.         nsacw, nseml, ncact 
    342. 

  342.   REAL :: ifac, sfac
    343. 

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

  344. !
    345. 

  345. #define WSM_NO_CONDITIONAL_IN_VECTOR
    346. 

  346. #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
    347. 

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

  348. #endif
    349. 

  349. ! variables for optimization
    350. 

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

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

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

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

  354.   REAL dtcldden, rdelz, rdtcld
    355. 

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

  356.   REAL  ::                                                        &
    357. 

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

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

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

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

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

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

  363.             nfrzdtr, nfrzdtc,                                     &
    364. 

  364.             taucon, lencon, lenconcr,                             &
    365. 

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

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

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

  368.   REAL :: temp 
    369. 

  369.   REAL  :: holdc, holdci
    370. 

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

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

  372. ! Temporaries used for inlining fpvs function
    373. 

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

  374.   REAL  :: logtr
    375. 

  375. !
    376. 

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

  377. !   compute internal functions
    378. 

  378. !
    379. 

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

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

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

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

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

  384. !
    385. 

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

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

  387. !
    388. 

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

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

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

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

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

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

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

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

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

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

  398. !
    399. 

  399. !
    400. 

  400.       idim = ite-its+1
    401. 

  401.       kdim = kte-kts+1
    402. 

  402. !
    403. 

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

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

  405. !
    406. 

  406.       do k = kts, kte
    407. 

  407.         do i = its, ite
    408. 

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

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

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

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

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

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

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

  415.         enddo
    416. 

  416.       enddo
    417. 

  417. !
    418. 

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

  419. !     changes during microphysical process calculation
    420. 

  420. !     emanuel(1994)
    421. 

  421. !
    422. 

  422.       do k = kts, kte
    423. 

  423.         do i = its, ite
    424. 

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

  425.           xl(i,k) = xlcal(t(i,k))
    426. 

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

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

  428.         enddo
    429. 

  429.       enddo
    430. 

  430. !
    431. 

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

  432. !    initialize the surface rain, snow
    433. 

  433. !
    434. 

  434.       do i = its, ite
    435. 

  435.         rainncv(i) = 0.
    436. 

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

  437.         sr(i) = 0.
    438. 

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

  439.         tstepsnow(i) = 0.
    440. 

  440.       enddo
    441. 

  441. !
    442. 

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

  443. !     compute the minor time steps.
    444. 

  444. !
    445. 

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

  446.       dtcld = delt/loops
    447. 

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

  448. !
    449. 

  449.       do loop = 1,loops
    450. 

  450. !
    451. 

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

  452. !     initialize the large scale variables
    453. 

  453. !
    454. 

  454.       do i = its, ite
    455. 

  455.         mstep(i) = 1
    456. 

  456.         mnstep(i) = 1
    457. 

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

  458.       enddo
    459. 

  459. !
    460. 

  460. !     do k = kts, kte
    461. 

  461. !       do i = its, ite
    462. 

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

  463. !       enddo
    464. 

  464. !     enddo
    465. 

  465.       do k = kts, kte
    466. 

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

  467.         do i = its, ite
    468. 

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

  469.         enddo
    470. 

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

  471.       enddo
    472. 

  472. !
    473. 

  473. ! Inline expansion for fpvs
    474. 

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

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

  476.       hsub = xls
    477. 

  477.       hvap = xlv0
    478. 

  478.       cvap = cpv
    479. 

  479.       ttp=t0c+0.01
    480. 

  480.       dldt=cvap-cliq
    481. 

  481.       xa=-dldt/rv
    482. 

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

  483.       dldti=cvap-cice
    484. 

  484.       xai=-dldti/rv
    485. 

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

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

  487. #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
    488. 

  488.       do k = kts, kte
    489. 

  489.         do i = its, ite
    490. 

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

  491.             xal(i) = xai
    492. 

  492.             xbl(i) = xbi
    493. 

  493.           else
    494. 

  494.             xal(i) = xa
    495. 

  495.             xbl(i) = xb
    496. 

  496.           endif
    497. 

  497.         enddo
    498. 

  498.         do i = its, ite
    499. 

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

  500.           logtr=log(tr)
    501. 

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

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

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

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

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

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

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

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

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

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

  511.         enddo
    512. 

  512.       enddo
    513. 

  513. #else
    514. 

  514.       do k = kts, kte
    515. 

  515.         do i = its, ite
    516. 

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

  517.           logtr=log(tr)
    518. 

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

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

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

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

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

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

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

  525.           else
    526. 

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

  527.           endif
    528. 

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

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

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

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

  532.         enddo
    533. 

  533.       enddo
    534. 

  534. #endif
    535. 

  535. !
    536. 

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

  537. !     initialize the variables for microphysical physics
    538. 

  538. !
    539. 

  539. !
    540. 

  540.       do k = kts, kte
    541. 

  541.         do i = its, ite
    542. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  572.         enddo
    573. 

  573.       enddo
    574. 

  574. !
    575. 

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

  576. !     compute the fallout term:
    577. 

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

  578. !
    579. 

  579.       do k = kts, kte
    580. 

  580.         do i = its, ite
    581. 

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

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

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

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

  585.           else
    586. 

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

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

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

  589.           endif
    590. 

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

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

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

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

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

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

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

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

  598.         enddo
    599. 

  599.       enddo
    600. 

  600.       do k = kts, kte
    601. 

  601.         do i = its, ite
    602. 

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

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

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

  605.         enddo
    606. 

  606.       enddo
    607. 

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

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

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

  610. !     compute the fallout term:
    611. 

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

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

  613. !
    614. 

  614. ! vt update for qr and nr
    615. 

  615.       mstepmax = 1
    616. 

  616.       numdt = 1
    617. 

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

  618.         do i = its, ite
    619. 

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

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

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

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

  623.         enddo
    624. 

  624.       enddo
    625. 

  625.       do i = its, ite
    626. 

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

  627.       enddo
    628. 

  628. !
    629. 

  629.       do n = 1, mstepmax
    630. 

  630.         k = kte
    631. 

  631.         do i = its, ite
    632. 

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

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

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

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

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

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

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

  639.           endif
    640. 

  640.         enddo
    641. 

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

  642.           do i = its, ite
    643. 

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

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

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

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

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

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

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

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

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

  652.             endif
    653. 

  653.           enddo
    654. 

  654.         enddo
    655. 

  655.         do k = kts, kte
    656. 

  656.           do i = its, ite
    657. 

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

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

  659.           enddo
    660. 

  660.         enddo
    661. 

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

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

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

  664.           do i = its, ite
    665. 

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

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

  667.           enddo
    668. 

  668.         enddo
    669. 

  669.       enddo
    670. 

  670. ! for semi
    671. 

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

  672.         do i = its, ite
    673. 

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

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

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

  676.         enddo
    677. 

  677.       enddo
    678. 

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

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

  680.       do k = kts, kte
    681. 

  681.         do i = its, ite
    682. 

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

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

  684.         enddo
    685. 

  685.       enddo
    686. 

  686.       do i = its, ite
    687. 

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

  688.       enddo
    689. 

  689.       do k = kts, kte
    690. 

  690.         do i = its, ite
    691. 

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

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

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

  694.         enddo
    695. 

  695.       enddo
    696. 

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

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

  698. !
    699. 

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

  700.         do i = its, ite
    701. 

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

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

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

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

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

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

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

  708.             xlf = xlf0
    709. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  733.             endif
    734. 

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

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

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

  737.           endif
    738. 

  738.         enddo
    739. 

  739.       enddo
    740. 

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

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

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

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

  744.         do i = its, ite
    745. 

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

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

  747.           else
    748. 

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

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

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

  751.           endif
    752. 

  752.         enddo
    753. 

  753.       enddo
    754. 

  754. !
    755. 

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

  756. !
    757. 

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

  758.         do i = its, ite
    759. 

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

  760.         enddo
    761. 

  761.       enddo
    762. 

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

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

  764.       do k = kts, kte
    765. 

  765.         do i = its, ite
    766. 

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

  767.         enddo
    768. 

  768.       enddo
    769. 

  769.       do i = its, ite
    770. 

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

  771.       enddo
    772. 

  772. !
    773. 

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

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

  775. !
    776. 

  776.       do i = its, ite
    777. 

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

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

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

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

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

  782.         endif
    783. 

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

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

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

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

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

  788.           endif
    789. 

  789.         endif 
    790. 

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

  791.       enddo
    792. 

  792. !
    793. 

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

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

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

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

  797.       do k = kts, kte
    798. 

  798.         do i = its, ite
    799. 

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

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

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

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

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

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

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

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

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

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

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

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

  811.           endif
    812. 

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

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

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

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

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

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

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

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

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

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

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

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

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

  825.           endif
    826. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  842.             endif
    843. 

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

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

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

  846.            endif
    847. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  864.             endif
    865. 

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

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

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

  868.           endif
    869. 

  869.         enddo
    870. 

  870.       enddo
    871. 

  871. !
    872. 

  872.       do k = kts, kte
    873. 

  873.         do i = its, ite
    874. 

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

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

  876.         enddo
    877. 

  877.       enddo
    878. 

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

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

  880. !
    881. 

  881.       do k = kts, kte
    882. 

  882.         do i = its, ite
    883. 

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

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

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

  886.         enddo
    887. 

  887.       enddo
    888. 

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

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

  890.       do k = kts, kte
    891. 

  891.         do i = its, ite
    892. 

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

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

  894. ! for raindrop distribution
    895. 

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

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

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

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

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

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

  901.           else
    902. 

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

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

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

  905.           endif
    906. 

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

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

  908. ! for cloud-droplet distribution
    909. 

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

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

  911.         enddo
    912. 

  912.       enddo
    913. 

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

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

  915. !             heat conduction and vapor diffusion
    916. 

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

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

  918. !
    919. 

  919.       do k = kts, kte
    920. 

  920.         do i = its, ite
    921. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  937.         enddo 
    938. 

  938.       enddo 
    939. 

  939. !
    940. 

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

  941. !
    942. 

  942. ! warm rain processes
    943. 

  943. !
    944. 

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

  945. !
    946. 

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

  947. !
    948. 

  948.       do k = kts, kte
    949. 

  949.         do i = its, ite
    950. 

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

  951.           satdt = supsat/dtcld
    952. 

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

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

  954. !        (QC->QR)
    955. 

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

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

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

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

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

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

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

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

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

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

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

  966. !        (NC->NR)
    967. 

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

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

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

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

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

  972.           endif
    973. 

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

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

  975. !        (QC->QR)
    976. 

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

  977. !        (NC->)
    978. 

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

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

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

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

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

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

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

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

  986.             else
    987. 

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

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

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

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

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

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

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

  994.             endif
    995. 

  995.           endif
    996. 

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

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

  998. !        (NC->)
    999. 

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

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

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

  1002.           else
    1003. 

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

  1004.                         *rslopec3(i,k)
    1005. 

  1005.           endif
    1006. 

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

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

  1008. !        (NR->)
    1009. 

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

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

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

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

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

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

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

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

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

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

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

  1020.             else
    1021. 

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

  1022.             endif
    1023. 

  1023.           endif
    1024. 

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

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

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

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

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

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

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

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

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

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

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

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

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

  1037. !        (NR->NCCN)
    1038. 

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

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

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

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

  1042.               endif
    1043. 

  1043.             else
    1044. 

  1044. !
    1045. 

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

  1046.             endif
    1047. 

  1047.           endif
    1048. 

  1048.         enddo
    1049. 

  1049.       enddo
    1050. 

  1050. !
    1051. 

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

  1052. !
    1053. 

  1053. ! cold rain processes
    1054. 

  1054. !
    1055. 

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

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

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

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

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

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

  1061. !
    1062. 

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

  1063. !
    1064. 

  1064.       rdtcld = 1./dtcld
    1065. 

  1065.       do k = kts, kte
    1066. 

  1066.         do i = its, ite
    1067. 

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

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

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

  1070.           satdt = supsat/dtcld
    1071. 

  1071.           ifsat = 0
    1072. 

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

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

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

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

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

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

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

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

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

  1081. !
    1082. 

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

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

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

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

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

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

  1088. !---…
    1089.

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