PageRenderTime 61ms CodeModel.GetById 20ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_mp_wdm6.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2657 lines | 1935 code | 1 blank | 721 comment | 156 complexity | 487eded20e93a0325d996a56a299f8ab MD5 | raw file
Possible License(s): AGPL-1.0

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. MODULE module_mp_wdm6
    10. 

  10. !
    11. 

  11. !
    12. 

  12. !
    13. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  44. !
    45. 

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

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

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

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

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

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

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

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

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

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

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

  56. !
    57. 

  57.    REAL, SAVE ::                                           &
    58. 

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

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

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

  61.              g5pbro2,g7pbro2,pi,                           &
    62. 

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

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

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

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

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

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

  68.              rslopecmax,rslopec2max,rslopec3max,           &
    69. 

  69.              rslopermax,rslopesmax,rslopegmax,             &
    70. 

  70.              rsloperbmax,rslopesbmax,rslopegbmax,          &
    71. 

  71.              rsloper2max,rslopes2max,rslopeg2max,          &
    72. 

  72.              rsloper3max,rslopes3max,rslopeg3max
    73. 

  73. CONTAINS
    74. 

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

  75. !
    76. 

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

  77.                     nn, nc, nr,                            &
    78. 

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

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

  80.                     ep1, ep2, qmin,                        &
    81. 

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

  82.                     cliq,cice,psat,                        &
    83. 

  83.                     rain, rainncv,                         &
    84. 

  84.                     snow, snowncv,                         &
    85. 

  85.                     sr,                                    &
    86. 

  86.                     graupel, graupelncv,                   &
    87. 

  87.                     itimestep,                             &
    88. 

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

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

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

  91.                                                            )
    92. 

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

  93.   IMPLICIT NONE
    94. 

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

  95. !
    96. 

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

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

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

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

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

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

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

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

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

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

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

  107. !
    108. 

  108. !  WDM6 cloud scheme
    109. 

  109. !
    110. 

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

  111. !
    112. 

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

  113. !
    114. 

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

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

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

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

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

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

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

  121. !             
    122. 

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

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

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

  125. !
    126. 

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

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

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

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

  130.         INTENT(INOUT) ::                                          &
    131. 

  131.                                                               th, &
    132. 

  132.                                                                q, &
    133. 

  133.                                                               qc, &
    134. 

  134.                                                               qi, &
    135. 

  135.                                                               qr, &
    136. 

  136.                                                               qs, &
    137. 

  137.                                                               qg, &
    138. 

  138.                                                               nn, & 
    139. 

  139.                                                               nc, &
    140. 

  140.                                                               nr
    141. 

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

  142.         INTENT(IN   ) ::                                          &
    143. 

  143.                                                              den, &
    144. 

  144.                                                              pii, &
    145. 

  145.                                                                p, &
    146. 

  146.                                                             delz
    147. 

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

  148.                                                                g, &
    149. 

  149.                                                               rd, &
    150. 

  150.                                                               rv, &
    151. 

  151.                                                              t0c, &
    152. 

  152.                                                             den0, &
    153. 

  153.                                                              cpd, &
    154. 

  154.                                                              cpv, &
    155. 

  155.                                                             ccn0, &
    156. 

  156.                                                              ep1, &
    157. 

  157.                                                              ep2, &
    158. 

  158.                                                             qmin, &
    159. 

  159.                                                              XLS, &
    160. 

  160.                                                             XLV0, &
    161. 

  161.                                                             XLF0, &
    162. 

  162.                                                             cliq, &
    163. 

  163.                                                             cice, &
    164. 

  164.                                                             psat, &
    165. 

  165.                                                             denr
    166. 

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

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

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

  169.                                                          rainncv, &
    170. 

  170.                                                               sr
    171. 

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

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

  173.                                                          snowncv
    174. 

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

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

  176.                                                         graupelncv
    177. 

  177. ! LOCAL VAR
    178. 

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

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

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

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

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

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

  184.         DO j=jms,jme
    185. 

  185.            DO k=kms,kme    
    186. 

  186.            DO i=ims,ime
    187. 

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

  188.            ENDDO
    189. 

  189.            ENDDO
    190. 

  190.         ENDDO
    191. 

  191.       ENDIF
    192. 

  192. !
    193. 

  193.       DO j=jts,jte
    194. 

  194.          DO k=kts,kte
    195. 

  195.          DO i=its,ite
    196. 

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

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

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

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

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

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

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

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

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

  205.          ENDDO
    206. 

  206.          ENDDO
    207. 

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

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

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

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

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

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

  213.                     ,ep1, ep2, qmin                               &
    214. 

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

  215.                     ,cliq,cice,psat                               &
    216. 

  216.                     ,j                                            &
    217. 

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

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

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

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

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

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

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

  224.                                                                    )
    225. 

  225.          DO K=kts,kte
    226. 

  226.          DO I=its,ite
    227. 

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

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

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

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

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

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

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

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

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

  236.          ENDDO
    237. 

  237.          ENDDO
    238. 

  238.       ENDDO
    239. 

  239.   END SUBROUTINE wdm6
    240. 

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

  241. !
    242. 

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

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

  244.                    ,ep1, ep2, qmin                                &
    245. 

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

  246.                    ,cliq,cice,psat                                &
    247. 

  247.                    ,lat                                           &
    248. 

  248.                    ,rain,rainncv                                  &
    249. 

  249.                    ,sr                                            &
    250. 

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

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

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

  253.                    ,snow,snowncv                                  &
    254. 

  254.                    ,graupel,graupelncv                            &
    255. 

  255.                                                                   )
    256. 

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

  257.   IMPLICIT NONE
    258. 

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

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

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

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

  262.                                       lat
    263. 

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

  264.         INTENT(INOUT) ::                                          &
    265. 

  265.                                                                t
    266. 

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

  267.         INTENT(INOUT) ::                                          &
    268. 

  268.                                                              qci
    269. 

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

  270.         INTENT(INOUT) ::                                          &
    271. 

  271.                                                              qrs, &
    272. 

  272.                                                              ncr 
    273. 

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

  274.         INTENT(INOUT) ::                                          &
    275. 

  275.                                                                q
    276. 

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

  277.         INTENT(IN   ) ::                                          &
    278. 

  278.                                                              den, &
    279. 

  279.                                                                p, &
    280. 

  280.                                                             delz
    281. 

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

  282.                                                                g, &
    283. 

  283.                                                              cpd, &
    284. 

  284.                                                              cpv, &
    285. 

  285.                                                             ccn0, &
    286. 

  286.                                                              t0c, &
    287. 

  287.                                                             den0, &
    288. 

  288.                                                               rd, &
    289. 

  289.                                                               rv, &
    290. 

  290.                                                              ep1, &
    291. 

  291.                                                              ep2, &
    292. 

  292.                                                             qmin, &
    293. 

  293.                                                              XLS, &
    294. 

  294.                                                             XLV0, &
    295. 

  295.                                                             XLF0, &
    296. 

  296.                                                             cliq, &
    297. 

  297.                                                             cice, &
    298. 

  298.                                                             psat, &
    299. 

  299.                                                             denr
    300. 

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

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

  302.                                                          rainncv, &
    303. 

  303.                                                               sr
    304. 

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

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

  306.                                                          snowncv
    307. 

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

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

  309.                                                       graupelncv
    310. 

  310. ! LOCAL VAR
    311. 

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

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

  313.         falk, fall, work1, qrs_tmp
    314. 

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

  315.         rslopec, rslopec2,rslopec3 
    316. 

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

  317.         avedia 
    318. 

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

  319.         workn,falln,falkn
    320. 

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

  321.         worka,workr
    322. 

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

  323.         den_tmp, delz_tmp, ncr_tmp 
    324. 

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

  325.         falkc, work1c, work2c, fallc
    326. 

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

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

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

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

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

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

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

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

  334.         nseml, ngeml, ncact 
    335. 

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

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

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

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

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

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

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

  342.   REAL :: gfac, sfac
    343. 

  343. ! variables for optimization
    344. 

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

  345.   REAL :: temp
    346. 

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

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

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

  349.   REAL  ::                                                        &
    350. 

  350.             cpmcal, xlcal, lamdac,                                &
    351. 

  351.             diffus,                                               &
    352. 

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

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

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

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

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

  357.             fallsum, fallsum_qsi, fallsum_qg,                     &
    358. 

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

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

  360.             nfrzdtr, nfrzdtc,                                     &
    361. 

  361.             taucon, lencon, lenconcr,                       &
    362. 

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

  363.   REAL  :: vt2ave
    364. 

  364.   REAL  :: holdc, holdci
    365. 

  365. !
    366. 

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

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

  368. ! Temporaries used for inlining fpvs function
    369. 

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

  370. !
    371. 

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

  372. !   compute internal functions
    373. 

  373. !
    374. 

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

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

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

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

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

  379. !
    380. 

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

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

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

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

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

  385. !
    386. 

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

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

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

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

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

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

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

  393. !
    394. 

  394.       idim = ite-its+1
    395. 

  395.       kdim = kte-kts+1
    396. 

  396. !
    397. 

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

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

  399. !
    400. 

  400.       do k = kts, kte
    401. 

  401.         do i = its, ite
    402. 

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

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

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

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

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

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

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

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

  410.         enddo
    411. 

  411.       enddo
    412. 

  412. !
    413. 

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

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

  415. !     changes during microphysical process calculation
    416. 

  416. !     emanuel(1994)
    417. 

  417. !
    418. 

  418.       do k = kts, kte
    419. 

  419.         do i = its, ite
    420. 

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

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

  422.         enddo
    423. 

  423.       enddo
    424. 

  424.       do k = kts, kte
    425. 

  425.         do i = its, ite
    426. 

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

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

  428.         enddo
    429. 

  429.       enddo
    430. 

  430. !
    431. 

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

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

  433. !
    434. 

  434.       do i = its, ite
    435. 

  435.         rainncv(i) = 0.
    436. 

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

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

  438.         sr(i) = 0.
    439. 

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

  440.         tstepsnow(i) = 0.
    441. 

  441.         tstepgraup(i) = 0.
    442. 

  442.       enddo
    443. 

  443. !
    444. 

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

  445. !     compute the minor time steps.
    446. 

  446. !
    447. 

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

  448.       dtcld = delt/loops
    449. 

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

  450. !
    451. 

  451.       do loop = 1,loops
    452. 

  452. !
    453. 

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

  454. !     initialize the large scale variables
    455. 

  455. !
    456. 

  456.       do i = its, ite
    457. 

  457.         mstep(i) = 1
    458. 

  458.         mnstep(i) = 1
    459. 

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

  460.       enddo
    461. 

  461. !
    462. 

  462.       do k = kts, kte
    463. 

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

  464.         do i = its, ite
    465. 

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

  466.         enddo
    467. 

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

  468.       enddo
    469. 

  469. !
    470. 

  470. ! Inline expansion for fpvs
    471. 

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

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

  473.       hsub = xls
    474. 

  474.       hvap = xlv0
    475. 

  475.       cvap = cpv
    476. 

  476.       ttp=t0c+0.01
    477. 

  477.       dldt=cvap-cliq
    478. 

  478.       xa=-dldt/rv
    479. 

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

  480.       dldti=cvap-cice
    481. 

  481.       xai=-dldti/rv
    482. 

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

  483.       do k = kts, kte
    484. 

  484.         do i = its, ite
    485. 

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

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

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

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

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

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

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

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

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

  494.           else
    495. 

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

  496.           endif
    497. 

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

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

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

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

  501.         enddo
    502. 

  502.       enddo
    503. 

  503. !
    504. 

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

  505. !     initialize the variables for microphysical physics
    506. 

  506. !
    507. 

  507. !
    508. 

  508.       do k = kts, kte
    509. 

  509.         do i = its, ite
    510. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  563.         enddo
    564. 

  564.       enddo
    565. 

  565.       do k = kts, kte
    566. 

  566.         do i = its, ite
    567. 

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

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

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

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

  571.           else
    572. 

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

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

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

  575.           endif
    576. 

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

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

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

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

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

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

  582.         enddo
    583. 

  583.       enddo
    584. 

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

  585. !     compute the fallout term:
    586. 

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

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

  588.       do k = kts, kte
    589. 

  589.         do i = its, ite
    590. 

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

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

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

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

  594.         enddo
    595. 

  595.       enddo
    596. 

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

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

  598. !
    599. 

  599. ! vt update for qr and nr
    600. 

  600.       mstepmax = 1
    601. 

  601.       numdt = 1
    602. 

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

  603.         do i = its, ite
    604. 

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

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

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

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

  608.         enddo
    609. 

  609.       enddo
    610. 

  610.       do i = its, ite
    611. 

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

  612.       enddo
    613. 

  613. !
    614. 

  614.       do n = 1, mstepmax
    615. 

  615.         k = kte
    616. 

  616.         do i = its, ite
    617. 

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

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

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

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

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

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

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

  624.           endif
    625. 

  625.         enddo
    626. 

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

  627.           do i = its, ite
    628. 

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

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

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

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

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

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

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

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

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

  637.             endif
    638. 

  638.           enddo
    639. 

  639.         enddo
    640. 

  640.         do k = kts, kte
    641. 

  641.           do i = its, ite
    642. 

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

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

  644.           enddo
    645. 

  645.         enddo
    646. 

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

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

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

  649.           do i = its, ite
    650. 

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

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

  652.           enddo
    653. 

  653.         enddo
    654. 

  654.       enddo
    655. 

  655. ! for semi
    656. 

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

  657.         do i = its, ite
    658. 

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

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

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

  661.                       /qsum(i,k)
    662. 

  662.           else
    663. 

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

  664.           endif
    665. 

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

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

  667.         enddo
    668. 

  668.       enddo
    669. 

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

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

  671.       do k = kts, kte
    672. 

  672.         do i = its, ite
    673. 

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

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

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

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

  677.         enddo
    678. 

  678.       enddo
    679. 

  679.       do i = its, ite
    680. 

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

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

  682.       enddo
    683. 

  683.       do k = kts, kte
    684. 

  684.         do i = its, ite
    685. 

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

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

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

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

  689.         enddo
    690. 

  690.       enddo
    691. 

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

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

  693. !
    694. 

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

  695.         do i = its, ite
    696. 

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

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

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

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

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

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

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

  703.             xlf = xlf0
    704. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  721.               endif
    722. 

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

  723.             endif
    724. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  743.               endif
    744. 

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

  745.             endif
    746. 

  746.           endif
    747. 

  747.         enddo
    748. 

  748.       enddo
    749. 

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

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

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

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

  753.         do i = its, ite
    754. 

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

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

  756.           else
    757. 

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

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

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

  760.           endif
    761. 

  761.         enddo
    762. 

  762.       enddo
    763. 

  763. !
    764. 

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

  765. !
    766. 

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

  767.         do i = its, ite
    768. 

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

  769.         enddo
    770. 

  770.       enddo
    771. 

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

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

  773.       do k = kts, kte
    774. 

  774.         do i = its, ite
    775. 

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

  776.         enddo
    777. 

  777.       enddo
    778. 

  778.       do i = its, ite
    779. 

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

  780.       enddo
    781. 

  781. !
    782. 

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

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

  784. !
    785. 

  785.       do i = its, ite
    786. 

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

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

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

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

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

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

  792.         endif
    793. 

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

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

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

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

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

  798.           ENDIF
    799. 

  799.         ENDIF
    800. 

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

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

  802.                             + tstepgraup(i)
    803. 

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

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

  805.                             + graupelncv(i)
    806. 

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

  807.           ENDIF
    808. 

  808.         ENDIF
    809. 

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

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

  811.       enddo
    812. 

  812. !
    813. 

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

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

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

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

  817.       do k = kts, kte
    818. 

  818.         do i = its, ite
    819. 

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

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

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

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

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

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

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

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

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

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

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

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

  831.           endif
    832. 

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

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

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

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

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

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

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

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

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

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

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

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

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

  845.           endif
    846. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  863.             endif
    864. 

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

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

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

  867.           endif 
    868. 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  885.             endif
    886. 

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

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

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

  889.           endif
    890. 

  890.         enddo
    891. 

  891.       enddo
    892. 

  892. !
    893. 

  893.       do k = kts, kte
    894. 

  894.         do i = its, ite
    895. 

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

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

  897.         enddo
    898. 

  898.       enddo
    899. 

  899. !
    900. 

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

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

  902. !
    903. 

  903.       do k = kts, kte
    904. 

  904.         do i = its, ite
    905. 

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

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

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

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

  909.         enddo
    910. 

  910.       enddo
    911. 

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

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

  913.       do k = kts, kte
    914. 

  914.         do i = its, ite
    915. 

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

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

  917. ! for raindrop distribution                          
    918. 

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

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

  920. !
    921. 

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

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

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

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

  925.           else
    926. 

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

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

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

  929.           endif
    930. 

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

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

  932. ! for cloud-droplet distribution
    933. 

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

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

  935.         enddo
    936. 

  936.       enddo
    937. 

  937. !
    938. 

  938.       do k = kts, kte
    939. 

  939.         do i = its, ite
    940. 

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

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

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

  943.         enddo
    944. 

  944.       enddo
    945. 

  945. !
    946. 

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

  947. !
    948. 

  948. ! warm rain processes
    949. 

  949. !
    950. 

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

  951. !
    952. 

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

  953. !
    954. 

  954.       do k = kts, kte
    955. 

  955.         do i = its, ite
    956. 

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

  957.           satdt = supsat/dtcld
    958. 

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

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

  960. !        (QC->QR)
    961. 

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

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

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

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

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

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

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

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

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

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

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

  972. !        (NC->NR)
    973. 

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

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

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

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

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

  978.           endif
    979. 

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

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

  981. !        (QC->QR)
    982. 

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

  983. !        (NC->)
    984. 

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

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

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

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

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

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

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

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

  992.             else
    993. 

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

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

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

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

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

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

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

  1000.             endif
    1001. 

  1001.           endif 
    1002. 

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

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

  1004. !        (NC->)
    1005. 

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

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

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

  1008.           else
    1009. 

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

  1010.                          *rslopec3(i,k)
    1011. 

  1011.           endif
    1012. 

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

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

  1014. !        (NR->)
    1015. 

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

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

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

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

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

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

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

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

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

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

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

  1026.             else
    1027. 

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

  1028.             endif
    1029. 

  1029.           endif
    1030. 

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

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

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

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

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

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

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

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

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

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

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

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

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

  1043. !        (NR->NCCN) 
    1044. 

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

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

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

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

  1048.               endif
    1049. 

  1049.             else
    1050. 

  1050. !
    1051. 

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

  1052.             endif
    1053. 

  1053.           endif
    1054. 

  1054.         enddo
    1055. 

  1055.       enddo
    1056. 

  1056. !
    1057. 

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

  1058. !
    1059. 

  1059. ! cold rain processes
    1060. 

  1060. !
    1061. 

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

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

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

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

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

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

  1067. !
    1068. 

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

  1069. !
    1070. 

  1070.       do k = kts, kte
    1071. 

  1071.         do i = its, ite
    1072. 

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

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

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

  1075.           satdt = supsat/dtcld
    1076. 

  1076.           ifsat = 0
    1077. 

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

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

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

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

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

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

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

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

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

  1086. !
    1087. 

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

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

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

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

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

  1092.           vt2g…
    1093.

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