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

/wrfv2_fire/phys/module_mp_wsm3.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1420 lines | 1023 code | 1 blank | 396 comment | 74 complexity | 0129a2020876822b6ce3d11a8e777a7a 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. #ifdef _ACCEL
    2. 

  2. #  include "module_mp_wsm3_accel.F"
    3. 

  3. #else
    4. 

  4. #if ( RWORDSIZE == 4 )
    5. 

  5. #  define VREC vsrec
    6. 

  6. #  define VSQRT vssqrt
    7. 

  7. #else
    8. 

  8. #  define VREC vrec
    9. 

  9. #  define VSQRT vsqrt
    10. 

  10. #endif
    11. 

  11. 

  12. MODULE module_mp_wsm3
    13. 

  13. !
    14. 

  14. !
    15. 

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

  16.    REAL, PARAMETER, PRIVATE :: n0r = 8.e6         ! intercept parameter rain
    17. 

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

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

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

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

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

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

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

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

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

  26.    REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4   ! limited maximum value for slope parameter of rain
    27. 

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

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

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

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

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

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

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

  34.    REAL, SAVE ::                                      &
    35. 

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

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

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

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

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

  40.              pidn0s,xlv1,pi,                          &
    41. 

  41.              rslopermax,rslopesmax,rslopegmax,        &
    42. 

  42.              rsloperbmax,rslopesbmax,rslopegbmax,     &
    43. 

  43.              rsloper2max,rslopes2max,rslopeg2max,     &
    44. 

  44.              rsloper3max,rslopes3max,rslopeg3max
    45. 

  45. !
    46. 

  46. ! Specifies code-inlining of fpvs function in WSM32D below. JM 20040507
    47. 

  47. !
    48. 

  48. CONTAINS
    49. 

  49. !===================================================================
    50. 

  50. !
    51. 

  51.   SUBROUTINE wsm3(th, q, qci, qrs                     &
    52. 

  52.                    , w, den, pii, p, delz             &
    53. 

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

  54.                    , ep1, ep2, qmin                   &
    55. 

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

  56.                    , cliq,cice,psat                   &
    57. 

  57.                    , rain, rainncv                    &
    58. 

  58.                    , snow, snowncv                    &
    59. 

  59.                    , sr                               &
    60. 

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

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

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

  63.                                                       )
    64. 

  64. !-------------------------------------------------------------------
    65. 

  65.   IMPLICIT NONE
    66. 

  66. !-------------------------------------------------------------------
    67. 

  67. !
    68. 

  68.   INTEGER,      INTENT(IN   )    ::                ids,ide, jds,jde, kds,kde , &
    69. 

  69.                                                    ims,ime, jms,jme, kms,kme , &
    70. 

  70.                                                    its,ite, jts,jte, kts,kte
    71. 

  71.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                              &
    72. 

  72.         INTENT(INOUT) ::                                                       &
    73. 

  73.                                                                           th,  &
    74. 

  74.                                                                            q,  &
    75. 

  75.                                                                           qci, &
    76. 

  76.                                                                           qrs
    77. 

  77.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                              &
    78. 

  78.         INTENT(IN   ) ::                                                    w, &
    79. 

  79.                                                                           den, &
    80. 

  80.                                                                           pii, &
    81. 

  81.                                                                             p, &
    82. 

  82.                                                                          delz
    83. 

  83.   REAL, INTENT(IN   ) ::                                                 delt, &
    84. 

  84.                                                                             g, &
    85. 

  85.                                                                            rd, &
    86. 

  86.                                                                            rv, &
    87. 

  87.                                                                           t0c, &
    88. 

  88.                                                                          den0, &
    89. 

  89.                                                                           cpd, &
    90. 

  90.                                                                           cpv, &
    91. 

  91.                                                                           ep1, &
    92. 

  92.                                                                           ep2, &
    93. 

  93.                                                                          qmin, &
    94. 

  94.                                                                           XLS, &
    95. 

  95.                                                                          XLV0, &
    96. 

  96.                                                                          XLF0, &
    97. 

  97.                                                                          cliq, &
    98. 

  98.                                                                          cice, &
    99. 

  99.                                                                          psat, &
    100. 

  100.                                                                          denr
    101. 

  101.   REAL, DIMENSION( ims:ime , jms:jme ),                                        &
    102. 

  102.         INTENT(INOUT) ::                                                 rain, &
    103. 

  103.                                                                       rainncv
    104. 

  104.   REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                              &
    105. 

  105.         INTENT(INOUT) ::                                                 snow, &
    106. 

  106.                                                                       snowncv, &
    107. 

  107.                                                                            sr
    108. 

  108. ! LOCAL VAR
    109. 

  109.   REAL, DIMENSION( its:ite , kts:kte ) ::                                   t
    110. 

  110.   INTEGER ::                                                            i,j,k
    111. 

  111. !-------------------------------------------------------------------
    112. 

  112.       DO j=jts,jte
    113. 

  113.          DO k=kts,kte
    114. 

  114.          DO i=its,ite
    115. 

  115.             t(i,k)=th(i,k,j)*pii(i,k,j)
    116. 

  116.          ENDDO
    117. 

  117.          ENDDO
    118. 

  118.          CALL wsm32D(t, q(ims,kms,j), qci(ims,kms,j)                           &
    119. 

  119.                     ,qrs(ims,kms,j),w(ims,kms,j), den(ims,kms,j)               &
    120. 

  120.                     ,p(ims,kms,j), delz(ims,kms,j)                             &
    121. 

  121.                     ,delt,g, cpd, cpv, rd, rv, t0c                             &
    122. 

  122.                     ,ep1, ep2, qmin                                            &
    123. 

  123.                     ,XLS, XLV0, XLF0, den0, denr                               &
    124. 

  124.                     ,cliq,cice,psat                                            &
    125. 

  125.                     ,j                                                         &
    126. 

  126.                     ,rain(ims,j), rainncv(ims,j)                               &
    127. 

  127.                     ,snow(ims,j),snowncv(ims,j)                                &
    128. 

  128.                     ,sr(ims,j)                                                 &
    129. 

  129.                     ,ids,ide, jds,jde, kds,kde                                 &
    130. 

  130.                     ,ims,ime, jms,jme, kms,kme                                 &
    131. 

  131.                     ,its,ite, jts,jte, kts,kte                                 &
    132. 

  132.                                                                                )
    133. 

  133.          DO K=kts,kte
    134. 

  134.          DO I=its,ite
    135. 

  135.             th(i,k,j)=t(i,k)/pii(i,k,j)
    136. 

  136.          ENDDO
    137. 

  137.          ENDDO
    138. 

  138.       ENDDO
    139. 

  139.   END SUBROUTINE wsm3
    140. 

  140. !===================================================================
    141. 

  141. !
    142. 

  142.   SUBROUTINE wsm32D(t, q                                                       &
    143. 

  143.                    ,qci, qrs,w, den, p, delz                                   &
    144. 

  144.                    ,delt,g, cpd, cpv, rd, rv, t0c                              &
    145. 

  145.                    ,ep1, ep2, qmin                                             &
    146. 

  146.                    ,XLS, XLV0, XLF0, den0, denr                                &
    147. 

  147.                    ,cliq,cice,psat                                             &
    148. 

  148.                    ,lat                                                        &
    149. 

  149.                    ,rain, rainncv                                              &
    150. 

  150.                    ,snow,snowncv                                               &
    151. 

  151.                    ,sr                                                         &
    152. 

  152.                    ,ids,ide, jds,jde, kds,kde                                  &
    153. 

  153.                    ,ims,ime, jms,jme, kms,kme                                  &
    154. 

  154.                    ,its,ite, jts,jte, kts,kte                                  &
    155. 

  155.                                                                                )
    156. 

  156. !-------------------------------------------------------------------
    157. 

  157.   IMPLICIT NONE
    158. 

  158. !-------------------------------------------------------------------
    159. 

  159. !
    160. 

  160. !  This code is a 3-class simple ice microphyiscs scheme (WSM3) of the 
    161. 

  161. !  Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
    162. 

  162. !  number concentration is a function of temperature, and seperate assumption
    163. 

  163. !  is developed, in which ice crystal number concentration is a function
    164. 

  164. !  of ice amount. A theoretical background of the ice-microphysics and related
    165. 

  165. !  processes in the WSMMPs are described in Hong et al. (2004).
    166. 

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

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

  168. !
    169. 

  169. !  WSM3 cloud scheme
    170. 

  170. !
    171. 

  171. !  Developed by Song-You Hong (Yonsei Univ.), Jimy Dudhia (NCAR) 
    172. 

  172. !             and Shu-Hua Chen (UC Davis) 
    173. 

  173. !             Summer 2002
    174. 

  174. !
    175. 

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

  176. !             Summer 2003
    177. 

  177. !
    178. 

  178. !  History :  semi-lagrangian scheme sedimentation(JH), and clean up
    179. 

  179. !             Hong, August 2009
    180. 

  180. !
    181. 

  181. !  Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
    182. 

  182. !             Dudhia (D89, 1989) J. Atmos. Sci.
    183. 

  183. !             Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
    184. 

  184. !             Juang and Hong (JH, 2010) Mon. Wea. Rev.
    185. 

  185. !
    186. 

  186.   INTEGER,      INTENT(IN   )    ::                 ids,ide, jds,jde, kds,kde, &
    187. 

  187.                                                     ims,ime, jms,jme, kms,kme, &
    188. 

  188.                                                     its,ite, jts,jte, kts,kte, &
    189. 

  189.                                                     lat
    190. 

  190.   REAL, DIMENSION( its:ite , kts:kte ),                                        &
    191. 

  191.         INTENT(INOUT) ::                                                       &
    192. 

  192.                                                                             t
    193. 

  193.   REAL, DIMENSION( ims:ime , kms:kme ),                                        &
    194. 

  194.         INTENT(INOUT) ::                                                       &
    195. 

  195.                                                                             q, &
    196. 

  196.                                                                           qci, &
    197. 

  197.                                                                           qrs
    198. 

  198.   REAL, DIMENSION( ims:ime , kms:kme ),                                        &
    199. 

  199.         INTENT(IN   ) ::                                                    w, &
    200. 

  200.                                                                           den, &
    201. 

  201.                                                                             p, &
    202. 

  202.                                                                          delz
    203. 

  203.   REAL, INTENT(IN   ) ::                                                 delt, &
    204. 

  204.                                                                             g, &
    205. 

  205.                                                                           cpd, &
    206. 

  206.                                                                           cpv, &
    207. 

  207.                                                                           t0c, &
    208. 

  208.                                                                          den0, &
    209. 

  209.                                                                            rd, &
    210. 

  210.                                                                            rv, &
    211. 

  211.                                                                           ep1, &
    212. 

  212.                                                                           ep2, &
    213. 

  213.                                                                          qmin, &
    214. 

  214.                                                                           XLS, &
    215. 

  215.                                                                          XLV0, &
    216. 

  216.                                                                          XLF0, &
    217. 

  217.                                                                          cliq, &
    218. 

  218.                                                                          cice, &
    219. 

  219.                                                                          psat, &
    220. 

  220.                                                                          denr
    221. 

  221.   REAL, DIMENSION( ims:ime ),                                                  &
    222. 

  222.         INTENT(INOUT) ::                                                 rain, &
    223. 

  223.                                                                       rainncv
    224. 

  224.   REAL, DIMENSION( ims:ime ),     OPTIONAL,                                    &
    225. 

  225.         INTENT(INOUT) ::                                                 snow, &
    226. 

  226.                                                                       snowncv, &
    227. 

  227.                                                                            sr
    228. 

  228. ! LOCAL VAR
    229. 

  229.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    230. 

  230.                                                                            rh, &
    231. 

  231.                                                                            qs, &
    232. 

  232.                                                                        denfac, &
    233. 

  233.                                                                        rslope, &
    234. 

  234.                                                                       rslope2, &
    235. 

  235.                                                                       rslope3, &
    236. 

  236.                                                                       qrs_tmp, &
    237. 

  237.                                                                       den_tmp, &
    238. 

  238.                                                                      delz_tmp, &
    239. 

  239.                                                                       rslopeb
    240. 

  240.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    241. 

  241.                                                                          pgen, &
    242. 

  242.                                                                          pisd, &
    243. 

  243.                                                                          paut, &
    244. 

  244.                                                                          pacr, &
    245. 

  245.                                                                          pres, &
    246. 

  246.                                                                          pcon
    247. 

  247.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    248. 

  248.                                                                          fall, &
    249. 

  249.                                                                          falk, &
    250. 

  250.                                                                            xl, &
    251. 

  251.                                                                           cpm, &
    252. 

  252.                                                                         work1, &
    253. 

  253.                                                                         work2, &
    254. 

  254.                                                                           xni, &
    255. 

  255.                                                                           qs0, &
    256. 

  256.                                                                        denqci, &
    257. 

  257.                                                                        denqrs, &
    258. 

  258.                                                                        n0sfac, &
    259. 

  259.                                                                         falkc, &
    260. 

  260.                                                                        work1c, &
    261. 

  261.                                                                        work2c, &
    262. 

  262.                                                                         fallc
    263. 

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

  264.                                                                          delqi
    265. 

  265.   REAL, DIMENSION(its:ite) ::                                        tstepsnow
    266. 

  266.   INTEGER, DIMENSION( its:ite ) ::                                      kwork1,&
    267. 

  267.                                                                         kwork2
    268. 

  268.   INTEGER, DIMENSION( its:ite ) ::                                      mstep, &
    269. 

  269.                                                                         numdt
    270. 

  270.   LOGICAL, DIMENSION( its:ite ) :: flgcld
    271. 

  271.   REAL  ::                                                                     &
    272. 

  272.             cpmcal, xlcal, diffus,                                             &
    273. 

  273.             viscos, xka, venfac, conden, diffac,                               &
    274. 

  274.             x, y, z, a, b, c, d, e,                                            &
    275. 

  275.             fallsum, fallsum_qsi, vt2i,vt2s,acrfac,                            &      
    276. 

  276.             qdt, pvt, qik, delq, facq, qrsci, frzmlt,                          &
    277. 

  277.             snomlt, hold, holdrs, facqci, supcol, coeres,                      &
    278. 

  278.             supsat, dtcld, xmi, qciik, delqci, eacrs, satdt,                   &
    279. 

  279.             qimax, diameter, xni0, roqi0, supice,holdc, holdci
    280. 

  280.   INTEGER :: i, j, k, mstepmax,                                                &
    281. 

  281.             iprt, latd, lond, loop, loops, ifsat, kk, n, idim, kdim
    282. 

  282. ! Temporaries used for inlining fpvs function
    283. 

  283.   REAL  :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
    284. 

  284. ! variables for optimization
    285. 

  285.   REAL, DIMENSION( its:ite )    :: tvec1
    286. 

  286. !
    287. 

  287. !=================================================================
    288. 

  288. !   compute internal functions
    289. 

  289. !
    290. 

  290.       cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
    291. 

  291.       xlcal(x) = xlv0-xlv1*(x-t0c)
    292. 

  292. !----------------------------------------------------------------
    293. 

  293. !     diffus: diffusion coefficient of the water vapor
    294. 

  294. !     viscos: kinematic viscosity(m2s-1)
    295. 

  295. !     Optimizatin : A**B => exp(log(A)*(B))
    296. 

  296. !
    297. 

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

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

  299.       xka(x,y) = 1.414e3*viscos(x,y)*y
    300. 

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

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

  302.                      /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
    303. 

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

  304. !
    305. 

  305.       idim = ite-its+1
    306. 

  306.       kdim = kte-kts+1
    307. 

  307. !
    308. 

  308. !----------------------------------------------------------------
    309. 

  309. !     paddint 0 for negative values generated by dynamics
    310. 

  310. !
    311. 

  311.       do k = kts, kte
    312. 

  312.         do i = its, ite
    313. 

  313.           qci(i,k) = max(qci(i,k),0.0)
    314. 

  314.           qrs(i,k) = max(qrs(i,k),0.0)
    315. 

  315.         enddo
    316. 

  316.       enddo
    317. 

  317. !
    318. 

  318. !----------------------------------------------------------------
    319. 

  319. !     latent heat for phase changes and heat capacity. neglect the
    320. 

  320. !     changes during microphysical process calculation
    321. 

  321. !     emanuel(1994)
    322. 

  322. !
    323. 

  323.       do k = kts, kte
    324. 

  324.         do i = its, ite
    325. 

  325.           cpm(i,k) = cpmcal(q(i,k))
    326. 

  326.           xl(i,k) = xlcal(t(i,k))
    327. 

  327.         enddo
    328. 

  328.       enddo
    329. 

  329.       do k = kts, kte
    330. 

  330.         do i = its, ite
    331. 

  331.           delz_tmp(i,k) = delz(i,k)
    332. 

  332.           den_tmp(i,k) = den(i,k)
    333. 

  333.         enddo
    334. 

  334.       enddo
    335. 

  335. !
    336. 

  336. !----------------------------------------------------------------
    337. 

  337. !    initialize the surface rain, snow
    338. 

  338. !
    339. 

  339.       do i = its, ite
    340. 

  340.         rainncv(i) = 0.
    341. 

  341.         if(PRESENT (snowncv) .AND. PRESENT (snow)) snowncv(i) = 0.
    342. 

  342.         sr(i) = 0.
    343. 

  343. ! new local array to catch step snow
    344. 

  344.         tstepsnow(i) = 0.
    345. 

  345.       enddo
    346. 

  346. !
    347. 

  347. !----------------------------------------------------------------
    348. 

  348. !     compute the minor time steps.
    349. 

  349. !
    350. 

  350.       loops = max(nint(delt/dtcldcr),1)
    351. 

  351.       dtcld = delt/loops
    352. 

  352.       if(delt.le.dtcldcr) dtcld = delt
    353. 

  353. !
    354. 

  354.       do loop = 1,loops
    355. 

  355. !
    356. 

  356. !----------------------------------------------------------------
    357. 

  357. !     initialize the large scale variables
    358. 

  358. !
    359. 

  359.       do i = its, ite
    360. 

  360.         flgcld(i) = .true.
    361. 

  361.       enddo
    362. 

  362. !
    363. 

  363.       do k = kts, kte
    364. 

  364.         CALL VREC( tvec1(its), den(its,k), ite-its+1)
    365. 

  365.         do i = its, ite
    366. 

  366.           tvec1(i) = tvec1(i)*den0
    367. 

  367.         enddo
    368. 

  368.         CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
    369. 

  369.       enddo
    370. 

  370. !
    371. 

  371. ! Inline expansion for fpvs
    372. 

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

  373. !         qs0(i,k) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
    374. 

  374.       cvap = cpv
    375. 

  375.       hvap=xlv0
    376. 

  376.       hsub=xls
    377. 

  377.       ttp=t0c+0.01
    378. 

  378.       dldt=cvap-cliq
    379. 

  379.       xa=-dldt/rv
    380. 

  380.       xb=xa+hvap/(rv*ttp)
    381. 

  381.       dldti=cvap-cice
    382. 

  382.       xai=-dldti/rv
    383. 

  383.       xbi=xai+hsub/(rv*ttp)
    384. 

  384.       do k = kts, kte
    385. 

  385.         do i = its, ite
    386. 

  386.           tr=ttp/t(i,k)
    387. 

  387.           if(t(i,k).lt.ttp) then
    388. 

  388.             qs(i,k) =psat*(exp(log(tr)*(xai)))*exp(xbi*(1.-tr))
    389. 

  389.           else
    390. 

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

  391.           endif
    392. 

  392.           qs0(i,k)  =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
    393. 

  393.           qs0(i,k) = (qs0(i,k)-qs(i,k))/qs(i,k)
    394. 

  394.           qs(i,k) = min(qs(i,k),0.99*p(i,k))
    395. 

  395.           qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
    396. 

  396.           qs(i,k) = max(qs(i,k),qmin)
    397. 

  397.           rh(i,k) = max(q(i,k) / qs(i,k),qmin)
    398. 

  398.         enddo
    399. 

  399.       enddo
    400. 

  400. !
    401. 

  401. !----------------------------------------------------------------
    402. 

  402. !     initialize the variables for microphysical physics
    403. 

  403. !
    404. 

  404. !
    405. 

  405.       do k = kts, kte
    406. 

  406.         do i = its, ite
    407. 

  407.           pres(i,k) = 0.
    408. 

  408.           paut(i,k) = 0.
    409. 

  409.           pacr(i,k) = 0.
    410. 

  410.           pgen(i,k) = 0.
    411. 

  411.           pisd(i,k) = 0.
    412. 

  412.           pcon(i,k) = 0.
    413. 

  413.           fall(i,k) = 0.
    414. 

  414.           falk(i,k) = 0.
    415. 

  415.           fallc(i,k) = 0.
    416. 

  416.           falkc(i,k) = 0.
    417. 

  417.           xni(i,k) = 1.e3
    418. 

  418.         enddo
    419. 

  419.       enddo
    420. 

  420. !-------------------------------------------------------------
    421. 

  421. ! Ni: ice crystal number concentraiton   [HDC 5c]
    422. 

  422. !-------------------------------------------------------------
    423. 

  423.       do k = kts, kte
    424. 

  424.         do i = its, ite
    425. 

  425.           xni(i,k) = min(max(5.38e7                                            &
    426. 

  426.                     *exp(log((den(i,k)*max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6)
    427. 

  427.         enddo
    428. 

  428.       enddo
    429. 

  429. !
    430. 

  430. !----------------------------------------------------------------
    431. 

  431. !     compute the fallout term:
    432. 

  432. !     first, vertical terminal velosity for minor loops
    433. 

  433. !---------------------------------------------------------------
    434. 

  434.       do k = kts, kte
    435. 

  435.         do i = its, ite
    436. 

  436.           qrs_tmp(i,k) = qrs(i,k)
    437. 

  437.         enddo
    438. 

  438.       enddo
    439. 

  439.       call slope_wsm3(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
    440. 

  440.                       work1,its,ite,kts,kte)
    441. 

  441. !
    442. 

  442. !
    443. 

  443. !  forward semi-laglangian scheme (JH), PCM (piecewise constant),  (linear)
    444. 

  444. !
    445. 

  445.       do k = kte, kts, -1
    446. 

  446.         do i = its, ite
    447. 

  447.           denqrs(i,k) = den(i,k)*qrs(i,k)
    448. 

  448.         enddo
    449. 

  449.       enddo
    450. 

  450.       call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1,denqrs,   &
    451. 

  451.                            delqrs,dtcld,1,1)
    452. 

  452.       do k = kts, kte
    453. 

  453.         do i = its, ite
    454. 

  454.           qrs(i,k) = max(denqrs(i,k)/den(i,k),0.)
    455. 

  455.           fall(i,k) = denqrs(i,k)*work1(i,k)/delz(i,k)
    456. 

  456.         enddo
    457. 

  457.       enddo
    458. 

  458.       do i = its, ite
    459. 

  459.         fall(i,1) = delqrs(i)/delz(i,1)/dtcld
    460. 

  460.       enddo
    461. 

  461. !---------------------------------------------------------------
    462. 

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

  463. !---------------------------------------------------------------
    464. 

  464.       do k = kte, kts, -1
    465. 

  465.         do i = its, ite
    466. 

  466.           if(t(i,k).lt.t0c.and.qci(i,k).gt.0.) then
    467. 

  467.             xmi = den(i,k)*qci(i,k)/xni(i,k)
    468. 

  468.             diameter  = max(dicon * sqrt(xmi), 1.e-25)
    469. 

  469.             work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
    470. 

  470.           else
    471. 

  471.             work1c(i,k) = 0.
    472. 

  472.           endif
    473. 

  473.         enddo
    474. 

  474.       enddo
    475. 

  475. !
    476. 

  476. !  forward semi-laglangian scheme (JH), PCM (piecewise constant),  (linear)
    477. 

  477. !
    478. 

  478.       do k = kte, kts, -1
    479. 

  479.         do i = its, ite
    480. 

  480.           denqci(i,k) = den(i,k)*qci(i,k)
    481. 

  481.         enddo
    482. 

  482.       enddo
    483. 

  483.       call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci,  &
    484. 

  484.                            delqi,dtcld,1,0)
    485. 

  485.       do k = kts, kte
    486. 

  486.         do i = its, ite
    487. 

  487.           qci(i,k) = max(denqci(i,k)/den(i,k),0.)
    488. 

  488.         enddo
    489. 

  489.       enddo
    490. 

  490.       do i = its, ite
    491. 

  491.         fallc(i,1) = delqi(i)/delz(i,1)/dtcld
    492. 

  492.       enddo
    493. 

  493. !
    494. 

  494. !----------------------------------------------------------------
    495. 

  495. !     compute the freezing/melting term. [D89 B16-B17]
    496. 

  496. !     freezing occurs one layer above the melting level
    497. 

  497. !
    498. 

  498.       do i = its, ite
    499. 

  499.         mstep(i) = 0
    500. 

  500.       enddo
    501. 

  501.       do k = kts, kte
    502. 

  502. !
    503. 

  503.         do i = its, ite
    504. 

  504.           if(t(i,k).ge.t0c) then
    505. 

  505.             mstep(i) = k
    506. 

  506.           endif
    507. 

  507.         enddo
    508. 

  508.       enddo
    509. 

  509. !
    510. 

  510.       do i = its, ite
    511. 

  511.         kwork2(i) = mstep(i)
    512. 

  512.         kwork1(i) = mstep(i)
    513. 

  513.         if(mstep(i).ne.0) then
    514. 

  514.           if (w(i,mstep(i)).gt.0.) then
    515. 

  515.             kwork1(i) = mstep(i) + 1
    516. 

  516.           endif
    517. 

  517.         endif
    518. 

  518.       enddo
    519. 

  519. !
    520. 

  520.       do i = its, ite
    521. 

  521.         k  = kwork1(i)
    522. 

  522.         kk = kwork2(i)
    523. 

  523.         if(k*kk.ge.1) then
    524. 

  524.           qrsci = qrs(i,k) + qci(i,k)
    525. 

  525.           if(qrsci.gt.0..or.fall(i,kk).gt.0.) then
    526. 

  526.             frzmlt = min(max(-w(i,k)*qrsci/delz(i,k),-qrsci/dtcld),            &
    527. 

  527.                     qrsci/dtcld)
    528. 

  528.             snomlt = min(max(fall(i,kk)/den(i,kk),-qrs(i,k)/dtcld),            &
    529. 

  529.                     qrs(i,k)/dtcld)
    530. 

  530.             if(k.eq.kk) then
    531. 

  531.               t(i,k) = t(i,k) - xlf0/cpm(i,k)*(frzmlt+snomlt)*dtcld
    532. 

  532.             else
    533. 

  533.               t(i,k) = t(i,k) - xlf0/cpm(i,k)*frzmlt*dtcld
    534. 

  534.               t(i,kk) = t(i,kk) - xlf0/cpm(i,kk)*snomlt*dtcld
    535. 

  535.             endif
    536. 

  536.           endif
    537. 

  537.         endif
    538. 

  538.       enddo
    539. 

  539. !
    540. 

  540. !----------------------------------------------------------------
    541. 

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

  542. !
    543. 

  543.       do i = its, ite
    544. 

  544.         fallsum = fall(i,1)
    545. 

  545.         fallsum_qsi = 0.
    546. 

  546.         if((t0c-t(i,1)).gt.0) then
    547. 

  547.         fallsum = fallsum+fallc(i,1)
    548. 

  548.         fallsum_qsi = fall(i,1)+fallc(i,1)
    549. 

  549.         endif
    550. 

  550.         if(fallsum.gt.0.) then
    551. 

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

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

  553.         endif
    554. 

  554.         if(fallsum_qsi.gt.0.) then
    555. 

  555.           tstepsnow(i)   = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.            &
    556. 

  556.                            +tstepsnow(i)
    557. 

  557.         IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
    558. 

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

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

  560.         ENDIF
    561. 

  561.         endif
    562. 

  562. !       if(fallsum.gt.0.) sr(i) = snowncv(i)/(rainncv(i)+1.e-12)
    563. 

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

  564.       enddo
    565. 

  565. !
    566. 

  566. !----------------------------------------------------------------
    567. 

  567. !     update the slope parameters for microphysics computation 
    568. 

  568. !
    569. 

  569.       do k = kts, kte
    570. 

  570.         do i = its, ite
    571. 

  571.           qrs_tmp(i,k) = qrs(i,k)
    572. 

  572.         enddo
    573. 

  573.       enddo
    574. 

  574.       call slope_wsm3(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
    575. 

  575.                       work1,its,ite,kts,kte)
    576. 

  576. !
    577. 

  577. !     work1: the thermodynamic term in the denominator associated with
    578. 

  578. !             heat conduction and vapor diffusion
    579. 

  579. !     work2: parameter associated with the ventilation effects(y93)
    580. 

  580. !
    581. 

  581.       do k = kts, kte
    582. 

  582.         do i = its, ite
    583. 

  583.           if(t(i,k).ge.t0c) then
    584. 

  584.             work1(i,k) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k))
    585. 

  585.           else
    586. 

  586.             work1(i,k) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k))
    587. 

  587.           endif
    588. 

  588.           work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
    589. 

  589.         enddo
    590. 

  590.       enddo
    591. 

  591. !
    592. 

  592.       do k = kts, kte
    593. 

  593.         do i = its, ite
    594. 

  594.           supsat = max(q(i,k),qmin)-qs(i,k)
    595. 

  595.           satdt = supsat/dtcld
    596. 

  596.           if(t(i,k).ge.t0c) then
    597. 

  597. !
    598. 

  598. !===============================================================
    599. 

  599. !
    600. 

  600. ! warm rain processes
    601. 

  601. !
    602. 

  602. ! - follows the processes in RH83 and LFO except for autoconcersion
    603. 

  603. !
    604. 

  604. !===============================================================
    605. 

  605. !---------------------------------------------------------------
    606. 

  606. ! praut: auto conversion rate from cloud to rain [HDC 16]
    607. 

  607. !        (C->R)
    608. 

  608. !---------------------------------------------------------------
    609. 

  609.             if(qci(i,k).gt.qc0) then
    610. 

  610. !             paut(i,k) = qck1*qci(i,k)**(7./3.)
    611. 

  611.               paut(i,k) = qck1*exp(log(qci(i,k))*((7./3.)))
    612. 

  612.               paut(i,k) = min(paut(i,k),qci(i,k)/dtcld)
    613. 

  613.             endif
    614. 

  614. !---------------------------------------------------------------
    615. 

  615. ! pracw: accretion of cloud water by rain [HL A40] [D89 B15]
    616. 

  616. !        (C->R)
    617. 

  617. !---------------------------------------------------------------
    618. 

  618.             if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then
    619. 

  619.                 pacr(i,k) = min(pacrr*rslope3(i,k)*rslopeb(i,k)                &
    620. 

  620.                      *qci(i,k)*denfac(i,k),qci(i,k)/dtcld)
    621. 

  621.             endif
    622. 

  622. !---------------------------------------------------------------
    623. 

  623. ! prevp: evaporation/condensation rate of rain [HDC 14]
    624. 

  624. !        (V->R or R->V)
    625. 

  625. !---------------------------------------------------------------
    626. 

  626.             if(qrs(i,k).gt.0.) then
    627. 

  627.                 coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k))
    628. 

  628.                 pres(i,k) = (rh(i,k)-1.)*(precr1*rslope2(i,k)                  &
    629. 

  629.                          +precr2*work2(i,k)*coeres)/work1(i,k)
    630. 

  630.               if(pres(i,k).lt.0.) then
    631. 

  631.                 pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
    632. 

  632.                 pres(i,k) = max(pres(i,k),satdt/2)
    633. 

  633.               else
    634. 

  634.                 pres(i,k) = min(pres(i,k),satdt/2)
    635. 

  635.               endif
    636. 

  636.             endif
    637. 

  637.           else
    638. 

  638. !
    639. 

  639. !===============================================================
    640. 

  640. !
    641. 

  641. ! cold rain processes
    642. 

  642. !
    643. 

  643. ! - follows the revised ice microphysics processes in HDC
    644. 

  644. ! - the processes same as in RH83 and LFO behave
    645. 

  645. !   following ice crystal hapits defined in HDC, inclduing
    646. 

  646. !   intercept parameter for snow (n0s), ice crystal number
    647. 

  647. !   concentration (ni), ice nuclei number concentration
    648. 

  648. !   (n0i), ice diameter (d)
    649. 

  649. !
    650. 

  650. !===============================================================
    651. 

  651. !
    652. 

  652.             supcol = t0c-t(i,k)
    653. 

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

  654.             ifsat = 0
    655. 

  655. !-------------------------------------------------------------
    656. 

  656. ! Ni: ice crystal number concentraiton   [HDC 5c]
    657. 

  657. !-------------------------------------------------------------
    658. 

  658.             xni(i,k) = min(max(5.38e7                                          &
    659. 

  659.                     *exp(log((den(i,k)*max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6)
    660. 

  660.             eacrs = exp(0.07*(-supcol))
    661. 

  661.             if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then
    662. 

  662.               xmi = den(i,k)*qci(i,k)/xni(i,k)
    663. 

  663.               diameter  = min(dicon * sqrt(xmi),dimax)
    664. 

  664.               vt2i = 1.49e4*diameter**1.31
    665. 

  665.               vt2s = pvts*rslopeb(i,k)*denfac(i,k)
    666. 

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

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

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

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

  670.               acrfac = 2.*rslope3(i,k)+2.*diameter*rslope2(i,k)                &
    671. 

  671.                       +diameter**2*rslope(i,k)
    672. 

  672.               pacr(i,k) = min(pi*qci(i,k)*eacrs*n0s*n0sfac(i,k)                &
    673. 

  673.                        *abs(vt2s-vt2i)*acrfac/4.,qci(i,k)/dtcld)
    674. 

  674.             endif
    675. 

  675. !-------------------------------------------------------------
    676. 

  676. ! pidep: Deposition/Sublimation rate of ice [HDC 9]
    677. 

  677. !       (T<T0: V->I or I->V)
    678. 

  678. !-------------------------------------------------------------
    679. 

  679.             if(qci(i,k).gt.0.) then
    680. 

  680.               xmi = den(i,k)*qci(i,k)/xni(i,k)
    681. 

  681.               diameter = dicon * sqrt(xmi)
    682. 

  682.               pisd(i,k) = 4.*diameter*xni(i,k)*(rh(i,k)-1.)/work1(i,k)
    683. 

  683.               if(pisd(i,k).lt.0.) then
    684. 

  684.                 pisd(i,k) = max(pisd(i,k),satdt/2)
    685. 

  685.                 pisd(i,k) = max(pisd(i,k),-qci(i,k)/dtcld)
    686. 

  686.               else
    687. 

  687.                 pisd(i,k) = min(pisd(i,k),satdt/2)
    688. 

  688.               endif
    689. 

  689.               if(abs(pisd(i,k)).ge.abs(satdt)) ifsat = 1
    690. 

  690.             endif
    691. 

  691. !-------------------------------------------------------------
    692. 

  692. ! psdep: deposition/sublimation rate of snow [HDC 14]
    693. 

  693. !        (V->S or S->V)
    694. 

  694. !-------------------------------------------------------------
    695. 

  695.             if(qrs(i,k).gt.0..and.ifsat.ne.1) then
    696. 

  696.               coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k))
    697. 

  697.               pres(i,k) = (rh(i,k)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k)        &
    698. 

  698.                         +precs2*work2(i,k)*coeres)/work1(i,k)
    699. 

  699.               supice = satdt-pisd(i,k)
    700. 

  700.               if(pres(i,k).lt.0.) then
    701. 

  701.                 pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
    702. 

  702.                 pres(i,k) = max(max(pres(i,k),satdt/2),supice)
    703. 

  703.               else
    704. 

  704.                 pres(i,k) = min(min(pres(i,k),satdt/2),supice)
    705. 

  705.               endif
    706. 

  706.               if(abs(pisd(i,k)+pres(i,k)).ge.abs(satdt)) ifsat = 1
    707. 

  707.             endif
    708. 

  708. !-------------------------------------------------------------
    709. 

  709. ! pigen: generation(nucleation) of ice from vapor [HDC 7-8]
    710. 

  710. !       (T<T0: V->I)
    711. 

  711. !-------------------------------------------------------------
    712. 

  712.             if(supsat.gt.0.and.ifsat.ne.1) then
    713. 

  713.               supice = satdt-pisd(i,k)-pres(i,k)
    714. 

  714.               xni0 = 1.e3*exp(0.1*supcol)
    715. 

  715.               roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
    716. 

  716.               pgen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k),0.))/dtcld)
    717. 

  717.               pgen(i,k) = min(min(pgen(i,k),satdt),supice)
    718. 

  718.             endif
    719. 

  719. !-------------------------------------------------------------
    720. 

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

  721. !       (T<T0: I->S)
    722. 

  722. !-------------------------------------------------------------
    723. 

  723.             if(qci(i,k).gt.0.) then
    724. 

  724.               qimax = roqimax/den(i,k)
    725. 

  725.               paut(i,k) = max(0.,(qci(i,k)-qimax)/dtcld)
    726. 

  726.             endif
    727. 

  727.           endif
    728. 

  728.         enddo
    729. 

  729.       enddo
    730. 

  730. !
    731. 

  731. !----------------------------------------------------------------
    732. 

  732. !     check mass conservation of generation terms and feedback to the
    733. 

  733. !     large scale
    734. 

  734. !
    735. 

  735.       do k = kts, kte
    736. 

  736.         do i = its, ite
    737. 

  737.           qciik = max(qmin,qci(i,k))
    738. 

  738.           delqci = (paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k))*dtcld
    739. 

  739.           if(delqci.ge.qciik) then
    740. 

  740.             facqci = qciik/delqci
    741. 

  741.             paut(i,k) = paut(i,k)*facqci
    742. 

  742.             pacr(i,k) = pacr(i,k)*facqci
    743. 

  743.             pgen(i,k) = pgen(i,k)*facqci
    744. 

  744.             pisd(i,k) = pisd(i,k)*facqci
    745. 

  745.           endif
    746. 

  746.           qik = max(qmin,q(i,k))
    747. 

  747.           delq = (pres(i,k)+pgen(i,k)+pisd(i,k))*dtcld
    748. 

  748.           if(delq.ge.qik) then
    749. 

  749.             facq = qik/delq
    750. 

  750.             pres(i,k) = pres(i,k)*facq
    751. 

  751.             pgen(i,k) = pgen(i,k)*facq
    752. 

  752.             pisd(i,k) = pisd(i,k)*facq
    753. 

  753.           endif
    754. 

  754.           work2(i,k) = -pres(i,k)-pgen(i,k)-pisd(i,k)
    755. 

  755.           q(i,k) = q(i,k)+work2(i,k)*dtcld
    756. 

  756.           qci(i,k) = max(qci(i,k)-(paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k))    &
    757. 

  757.                     *dtcld,0.)
    758. 

  758.           qrs(i,k) = max(qrs(i,k)+(paut(i,k)+pacr(i,k)+pres(i,k))*dtcld,0.)
    759. 

  759.           if(t(i,k).lt.t0c) then
    760. 

  760.             t(i,k) = t(i,k)-xls*work2(i,k)/cpm(i,k)*dtcld
    761. 

  761.           else
    762. 

  762.             t(i,k) = t(i,k)-xl(i,k)*work2(i,k)/cpm(i,k)*dtcld
    763. 

  763.           endif
    764. 

  764.         enddo
    765. 

  765.       enddo
    766. 

  766. !
    767. 

  767.       cvap = cpv
    768. 

  768.       hvap = xlv0
    769. 

  769.       hsub = xls
    770. 

  770.       ttp=t0c+0.01
    771. 

  771.       dldt=cvap-cliq
    772. 

  772.       xa=-dldt/rv
    773. 

  773.       xb=xa+hvap/(rv*ttp)
    774. 

  774.       dldti=cvap-cice
    775. 

  775.       xai=-dldti/rv
    776. 

  776.       xbi=xai+hsub/(rv*ttp)
    777. 

  777.       do k = kts, kte
    778. 

  778.         do i = its, ite
    779. 

  779.           tr=ttp/t(i,k)
    780. 

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

  781.           qs(i,k) = min(qs(i,k),0.99*p(i,k))
    782. 

  782.           qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
    783. 

  783.           qs(i,k) = max(qs(i,k),qmin)
    784. 

  784.           denfac(i,k) = sqrt(den0/den(i,k))
    785. 

  785.         enddo
    786. 

  786.       enddo
    787. 

  787. !
    788. 

  788. !----------------------------------------------------------------
    789. 

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

  790. !     if there exists additional water vapor condensated/if
    791. 

  791. !     evaporation of cloud water is not enough to remove subsaturation
    792. 

  792. !
    793. 

  793.       do k = kts, kte
    794. 

  794.         do i = its, ite
    795. 

  795.           work1(i,k) = conden(t(i,k),q(i,k),qs(i,k),xl(i,k),cpm(i,k))
    796. 

  796.           work2(i,k) = qci(i,k)+work1(i,k)
    797. 

  797.           pcon(i,k) = min(max(work1(i,k),0.),max(q(i,k),0.))/dtcld
    798. 

  798.           if(qci(i,k).gt.0..and.work1(i,k).lt.0.and.t(i,k).gt.t0c)             &
    799. 

  799.             pcon(i,k) = max(work1(i,k),-qci(i,k))/dtcld
    800. 

  800.           q(i,k) = q(i,k)-pcon(i,k)*dtcld
    801. 

  801.           qci(i,k) = max(qci(i,k)+pcon(i,k)*dtcld,0.)
    802. 

  802.           t(i,k) = t(i,k)+pcon(i,k)*xl(i,k)/cpm(i,k)*dtcld
    803. 

  803.         enddo
    804. 

  804.       enddo
    805. 

  805. !
    806. 

  806. !----------------------------------------------------------------
    807. 

  807. !     padding for small values
    808. 

  808. !
    809. 

  809.       do k = kts, kte
    810. 

  810.         do i = its, ite
    811. 

  811.           if(qci(i,k).le.qmin) qci(i,k) = 0.0
    812. 

  812.           if(qrs(i,k).le.qcrmin) qrs(i,k) = 0.0
    813. 

  813.         enddo
    814. 

  814.       enddo
    815. 

  815. !
    816. 

  816.       enddo                  ! big loops
    817. 

  817.   END SUBROUTINE wsm32D
    818. 

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

  819.       REAL FUNCTION rgmma(x)
    820. 

  820. !-------------------------------------------------------------------
    821. 

  821.   IMPLICIT NONE
    822. 

  822. !-------------------------------------------------------------------
    823. 

  823. !     rgmma function:  use infinite product form
    824. 

  824.       REAL :: euler
    825. 

  825.       PARAMETER (euler=0.577215664901532)
    826. 

  826.       REAL :: x, y
    827. 

  827.       INTEGER :: i
    828. 

  828.       if(x.eq.1.)then
    829. 

  829.         rgmma=0.
    830. 

  830.           else
    831. 

  831.         rgmma=x*exp(euler*x)
    832. 

  832.         do i=1,10000
    833. 

  833.           y=float(i)
    834. 

  834.           rgmma=rgmma*(1.000+x/y)*exp(-x/y)
    835. 

  835.         enddo
    836. 

  836.         rgmma=1./rgmma
    837. 

  837.       endif
    838. 

  838.       END FUNCTION rgmma
    839. 

  839. !
    840. 

  840. !--------------------------------------------------------------------------
    841. 

  841.       REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
    842. 

  842. !--------------------------------------------------------------------------
    843. 

  843.       IMPLICIT NONE
    844. 

  844. !--------------------------------------------------------------------------
    845. 

  845.       REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti,         &
    846. 

  846.            xai,xbi,ttp,tr
    847. 

  847.       INTEGER ice
    848. 

  848. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    849. 

  849.       ttp=t0c+0.01
    850. 

  850.       dldt=cvap-cliq
    851. 

  851.       xa=-dldt/rv
    852. 

  852.       xb=xa+hvap/(rv*ttp)
    853. 

  853.       dldti=cvap-cice
    854. 

  854.       xai=-dldti/rv
    855. 

  855.       xbi=xai+hsub/(rv*ttp)
    856. 

  856.       tr=ttp/t
    857. 

  857.       if(t.lt.ttp.and.ice.eq.1) then
    858. 

  858.         fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
    859. 

  859.       else
    860. 

  860.         fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
    861. 

  861.       endif
    862. 

  862. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    863. 

  863.       END FUNCTION fpvs
    864. 

  864. !-------------------------------------------------------------------
    865. 

  865.   SUBROUTINE wsm3init(den0,denr,dens,cl,cpv,allowed_to_read)
    866. 

  866. !-------------------------------------------------------------------
    867. 

  867.   IMPLICIT NONE
    868. 

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

  869. !.... constants which may not be tunable
    870. 

  870.    REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
    871. 

  871.    LOGICAL, INTENT(IN) :: allowed_to_read
    872. 

  872. !  
    873. 

  873.    pi = 4.*atan(1.)
    874. 

  874.    xlv1 = cl-cpv
    875. 

  875. !  
    876. 

  876.    qc0  = 4./3.*pi*denr*r0**3*xncr/den0  ! 0.419e-3 -- .61e-3
    877. 

  877.    qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
    878. 

  878. !  
    879. 

  879.    bvtr1 = 1.+bvtr
    880. 

  880.    bvtr2 = 2.5+.5*bvtr
    881. 

  881.    bvtr3 = 3.+bvtr
    882. 

  882.    bvtr4 = 4.+bvtr
    883. 

  883.    g1pbr = rgmma(bvtr1)
    884. 

  884.    g3pbr = rgmma(bvtr3)
    885. 

  885.    g4pbr = rgmma(bvtr4)            ! 17.837825
    886. 

  886.    g5pbro2 = rgmma(bvtr2)          ! 1.8273
    887. 

  887.    pvtr = avtr*g4pbr/6.
    888. 

  888.    eacrr = 1.0
    889. 

  889.    pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
    890. 

  890.    precr1 = 2.*pi*n0r*.78
    891. 

  891.    precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
    892. 

  892.    xmmax = (dimax/dicon)**2
    893. 

  893.    roqimax = 2.08e22*dimax**8
    894. 

  894. !
    895. 

  895.    bvts1 = 1.+bvts
    896. 

  896.    bvts2 = 2.5+.5*bvts
    897. 

  897.    bvts3 = 3.+bvts
    898. 

  898.    bvts4 = 4.+bvts
    899. 

  899.    g1pbs = rgmma(bvts1)    !.8875
    900. 

  900.    g3pbs = rgmma(bvts3)
    901. 

  901.    g4pbs = rgmma(bvts4)    ! 12.0786
    902. 

  902.    g5pbso2 = rgmma(bvts2)
    903. 

  903.    pvts = avts*g4pbs/6.
    904. 

  904.    pacrs = pi*n0s*avts*g3pbs*.25
    905. 

  905.    precs1 = 4.*n0s*.65
    906. 

  906.    precs2 = 4.*n0s*.44*avts**.5*g5pbso2
    907. 

  907.    pidn0r =  pi*denr*n0r
    908. 

  908.    pidn0s =  pi*dens*n0s
    909. 

  909. !
    910. 

  910.    rslopermax = 1./lamdarmax
    911. 

  911.    rslopesmax = 1./lamdasmax
    912. 

  912.    rsloperbmax = rslopermax ** bvtr
    913. 

  913.    rslopesbmax = rslopesmax ** bvts
    914. 

  914.    rsloper2max = rslopermax * rslopermax
    915. 

  915.    rslopes2max = rslopesmax * rslopesmax
    916. 

  916.    rsloper3max = rsloper2max * rslopermax
    917. 

  917.    rslopes3max = rslopes2max * rslopesmax
    918. 

  918. !
    919. 

  919.   END SUBROUTINE wsm3init
    920. 

  920. !
    921. 

  921.       subroutine slope_wsm3(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,vt,its,ite,kts,kte)
    922. 

  922.   IMPLICIT NONE
    923. 

  923.   INTEGER       ::               its,ite, jts,jte, kts,kte
    924. 

  924.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    925. 

  925.                                                                           qrs, &
    926. 

  926.                                                                           den, &
    927. 

  927.                                                                        denfac, &
    928. 

  928.                                                                             t, &
    929. 

  929.                                                                        rslope, &
    930. 

  930.                                                                       rslopeb, &
    931. 

  931.                                                                       rslope2, &
    932. 

  932.                                                                       rslope3, &
    933. 

  933.                                                                            vt
    934. 

  934.   REAL, PARAMETER  :: t0c = 273.15
    935. 

  935.   REAL, DIMENSION( its:ite , kts:kte ) ::                                      &
    936. 

  936.                                                                        n0sfac
    937. 

  937.   REAL       ::  lamdar,lamdas,x, y, z, supcol, pvt
    938. 

  938.   integer :: i, j, k
    939. 

  939. !----------------------------------------------------------------
    940. 

  940. !     size distributions: (x=mixing ratio, y=air density):
    941. 

  941. !     valid for mixing ratio > 1.e-9 kg/kg.
    942. 

  942. !
    943. 

  943.       lamdar(x,y)=   sqrt(sqrt(pidn0r/(x*y)))      ! (pidn0r/(x*y))**.25
    944. 

  944.       lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
    945. 

  945. !
    946. 

  946.       do k = kts, kte
    947. 

  947.         do i = its, ite
    948. 

  948.           if(t(i,k).ge.t0c) then
    949. 

  949.             pvt = pvtr
    950. 

  950.             if(qrs(i,k).le.qcrmin)then
    951. 

  951.               rslope(i,k) = rslopermax
    952. 

  952.               rslopeb(i,k) = rsloperbmax
    953. 

  953.               rslope2(i,k) = rsloper2max
    954. 

  954.               rslope3(i,k) = rsloper3max
    955. 

  955.             else
    956. 

  956.               rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
    957. 

  957.               rslopeb(i,k) = exp(log(rslope(i,k))*(bvtr))
    958. 

  958.               rslope2(i,k) = rslope(i,k)*rslope(i,k)
    959. 

  959.               rslope3(i,k) = rslope2(i,k)*rslope(i,k)
    960. 

  960.             endif
    961. 

  961.           else
    962. 

  962.             supcol = t0c-t(i,k)
    963. 

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

  964.             pvt = pvts
    965. 

  965.             if(qrs(i,k).le.qcrmin)then
    966. 

  966.               rslope(i,k) = rslopesmax
    967. 

  967.               rslopeb(i,k) = rslopesbmax
    968. 

  968.               rslope2(i,k) = rslopes2max
    969. 

  969.               rslope3(i,k) = rslopes3max
    970. 

  970.             else
    971. 

  971.               rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
    972. 

  972.               rslopeb(i,k) = exp(log(rslope(i,k))*(bvts))
    973. 

  973.               rslope2(i,k) = rslope(i,k)*rslope(i,k)
    974. 

  974.               rslope3(i,k) = rslope2(i,k)*rslope(i,k)
    975. 

  975.             endif
    976. 

  976.           endif
    977. 

  977.           vt(i,k) = pvt*rslopeb(i,k)*denfac(i,k)
    978. 

  978.           if(qrs(i,k).le.0.0) vt(i,k) = 0.0
    979. 

  979.         enddo
    980. 

  980.       enddo
    981. 

  981.   END subroutine slope_wsm3
    982. 

  982. !-------------------------------------------------------------------
    983. 

  983.       SUBROUTINE nislfv_rain_pcm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
    984. 

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

  985. !
    986. 

  986. ! for non-iteration semi-Lagrangain forward advection for cloud
    987. 

  987. ! with mass conservation and positive definite advection
    988. 

  988. ! 2nd order interpolation with monotonic piecewise linear method
    989. 

  989. ! this routine is under assumption of decfl < 1 for semi_Lagrangian
    990. 

  990. !
    991. 

  991. ! dzl    depth of model layer in meter
    992. 

  992. ! wwl    terminal velocity at model layer m/s
    993. 

  993. ! rql    cloud density*mixing ration
    994. 

  994. ! precip precipitation
    995. 

  995. ! dt     time step
    996. 

  996. ! id     kind of precip: 0 test case; 1 raindrop
    997. 

  997. ! iter   how many time to guess mean terminal velocity: 0 pure forward.
    998. 

  998. !        0 : use departure wind for advection 
    999. 

  999. !        1 : use mean wind for advection 
    1000. 

  1000. !        > 1 : use mean wind after iter-1 iterations
    1001. 

  1001. !
    1002. 

  1002. ! author: hann-ming henry juang <henry.juang@noaa.gov>
    1003. 

  1003. !         implemented by song-you hong
    1004. 

  1004. !
    1005. 

  1005.       implicit none
    1006. 

  1006.       integer  im,km,id
    1007. 

  1007.       real  dt
    1008. 

  1008.       real  dzl(im,km),wwl(im,km),rql(im,km),precip(im)
    1009. 

  1009.       real  denl(im,km),denfacl(im,km),tkl(im,km)
    1010. 

  1010. !
    1011. 

  1011.       integer  i,k,n,m,kk,kb,kt,iter
    1012. 

  1012.       real  tl,tl2,qql,dql,qqd
    1013. 

  1013.       real  th,th2,qqh,dqh
    1014. 

  1014.       real  zsum,qsum,dim,dip,c1,con1,fa1,fa2
    1015. 

  1015.       real  zsumt,qsumt,zsumb,qsumb
    1016. 

  1016.       real  allold, allnew, zz, dzamin, cflmax, decfl
    1017. 

  1017.       real  dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
    1018. 

  1018.       real  den(km), denfac(km), tk(km)
    1019. 

  1019.       real  wi(km+1), zi(km+1), za(km+1)
    1020. 

  1020.       real  qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
    1021. 

  1021.       real  dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
    1022. 

  1022. !
    1023. 

  1023.       precip(:) = 0.0
    1024. 

  1024. !
    1025. 

  1025.       i_loop : do i=1,im
    1026. 

  1026. ! -----------------------------------
    1027. 

  1027.       dz(:) = dzl(i,:)
    1028. 

  1028.       qq(:) = rql(i,:)
    1029. 

  1029.       ww(:) = wwl(i,:)
    1030. 

  1030.       den(:) = denl(i,:)
    1031. 

  1031.       denfac(:) = denfacl(i,:)
    1032. 

  1032.       tk(:) = tkl(i,:)
    1033. 

  1033. ! skip for no precipitation for all layers
    1034. 

  1034.       allold = 0.0
    1035. 

  1035.       do k=1,km
    1036. 

  1036.         allold = allold + qq(k)
    1037. 

  1037.       enddo
    1038. 

  1038.       if(allold.le.0.0) then
    1039. 

  1039.         cycle i_loop
    1040. 

  1040.       endif
    1041. 

  1041. !
    1042. 

  1042. ! compute interface values
    1043. 

  1043.       zi(1)=0.0
    1044. 

  1044.       do k=1,km
    1045. 

  1045.         zi(k+1) = zi(k)+dz(k)
    1046. 

  1046.       enddo
    1047. 

  1047. !
    1048. 

  1048. ! save departure wind
    1049. 

  1049.       wd(:) = ww(:)
    1050. 

  1050.       n=1
    1051. 

  1051.  100  continue
    1052. 

  1052. ! pcm is 1st order, we should use 2nd order wi
    1053. 

  1053. ! 2nd order interpolation to get wi
    1054. 

  1054.       wi(1) = ww(1)
    1055. 

  1055.       do k=2,km
    1056. 

  1056.         wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
    1057. 

  1057.       enddo
    1058. 

  1058.       wi(km+1) = ww(km)
    1059. 

  1059. !
    1060. 

  1060. ! terminate of top of raingroup
    1061. 

  1061.       do k=2,km
    1062. 

  1062.         if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
    1063. 

  1063.       enddo
    1064. 

  1064. !
    1065. 

  1065. ! diffusivity of wi
    1066. 

  1066.       con1 = 0.05
    1067. 

  1067.       do k=km,1,-1
    1068. 

  1068.         decfl = (wi(k+1)-wi(k))*dt/dz(k)
    1069. 

  1069.         if( decfl .gt. con1 ) then
    1070. 

  1070.           wi(k) = wi(k+1) - con1*dz(k)/dt
    1071. 

  1071.         endif
    1072. 

  1072.       enddo
    1073. 

  1073. ! compute arrival point
    1074. 

  1074.       do k=1,km+1
    1075. 

  1075.         za(k) = zi(k) - wi(k)*dt
    1076. 

  1076.       enddo
    1077. 

  1077. !
    1078. 

  1078.       do k=1,km
    1079. 

  1079.         dza(k) = za(k+1)-za(k)
    1080. 

  1080.       enddo
    1081. 

  1081.       dza(km+1) = zi(km+1) - za(km+1)
    1082. 

  1082. !
    1083. 

  1083. ! computer deformation at arrival point
    1084. 

  1084.       do k=1,km
    1085. 

  1085.         qa(k) = qq(k)*dz(k)/dza(k)
    1086. 

  1086.         qr(k) = qa(k)/den(k)
    1087. 

  1087.       enddo
    1088. 

  1088.       qa(km+1) = 0.0
    1089. 

  1089. !     call maxmin(km,1,qa,' arrival points ')
    1090. 

  1090. !
    1091. 

  1091. ! compute arrival terminal velocity, and estimate mean terminal velocity
    1092. 

  1092. ! then back to use mean terminal velocity
    1093. 

  1093.       if( n.le.iter ) then
    1094. 

  1094.         call slope_wsm3(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
    1095. 

  1095.         if( n.eq.2 ) wa(1:km) = 0.5*(wa(1:km)+was(1:km))
    1096. 

  1096.         do k=1,km
    1097. 

  1097. !#ifdef DEBUG
    1098. 

  1098. !      print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
    1099. 

  1099. !#endif
    1100. 

  1100. ! mean wind is average of departure and new arrival winds
    1101. 

  1101.           ww(k) = 0.5* ( wd(k)+wa(k) )
    1102. 

  1102.         enddo
    1103. 

  1103.         was(:) = wa(:)
    1104. 

  1104.         n=n+1
    1105. 

  1105.         go to 100
    1106. 

  1106.       endif
    1107. 

  1107. !
    1108. 

  1108. !
    1109. 

  1109. ! interpolation to regular point
    1110. 

  1110.       qn = 0.0
    1111. 

  1111.       kb=1
    1112. 

  1112.       kt=1
    1113. 

  1113.       intp : do k=1,km
    1114. 

  1114.              kb=max(kb-1,1)
    1115. 

  1115.              kt=max(kt-1,1)
    1116. 

  1116. ! find kb and kt
    1117. 

  1117.              if( zi(k).ge.za(km+1) ) then
    1118. 

  1118.                exit intp
    1119. 

  1119.              else
    1120. 

  1120.                find_kb : do kk=kb,km
    1121. 

  1121.                          if( zi(k).le.za(kk+1) ) then
    1122. 

  1122.                            kb = kk
    1123. 

  1123.                            exit find_kb
    1124. 

  1124.                          else
    1125. 

  1125.                            cycle find_kb
    1126. 

  1126.                          endif
    1127. 

  1127.                enddo find_kb
    1128. 

  1128.                find_kt : do kk=kt,km
    1129. 

  1129.                          if( zi(k+1).le.za(kk) ) then
    1130. 

  1130.                            kt = kk
    1131. 

  1131.                            exit find_kt
    1132. 

  1132.                          else
    1133. 

  1133.                            cycle find_kt
    1134. 

  1134.                          endif
    1135. 

  1135.                enddo find_kt
    1136. 

  1136. ! compute q with piecewise constant method
    1137. 

  1137.                if( kt-kb.eq.1 ) then
    1138. 

  1138.                  qn(k) = qa(kb)
    1139. 

  1139.                else if( kt-kb.ge.2 ) then
    1140. 

  1140.                  zsumb = za(kb+1)-zi(k)
    1141. 

  1141.                  qsumb = qa(kb) * zsumb
    1142. 

  1142.                  zsumt = zi(k+1)-za(kt-1)
    1143. 

  1143.                  qsumt = qa(kt-1) * zsumt
    1144. 

  1144.                  qsum = 0.0
    1145. 

  1145.                  zsum = 0.0
    1146. 

  1146.                  if( kt-kb.ge.3 ) then
    1147. 

  1147.                  do m=kb+1,kt-2
    1148. 

  1148.                    qsum = qsum + qa(m) * dza(m)
    1149. 

  1149.                    zsum = zsum + dza(m)
    1150. 

  1150.                  enddo
    1151. 

  1151.                  endif
    1152. 

  1152.                  qn(k) = (qsumb+qsum+qsumt)/(zsumb+zsum+zsumt)
    1153. 

  1153.                endif
    1154. 

  1154.                cycle intp
    1155. 

  1155.              endif
    1156. 

  1156. !
    1157. 

  1157.        enddo intp
    1158. 

  1158. !
    1159. 

  1159. ! rain out
    1160. 

  1160.       sum_precip: do k=1,km
    1161. 

  1161.                     if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
    1162. 

  1162.                       precip(i) = precip(i) + qa(k)*dza(k)
    1163. 

  1163.                       cycle sum_precip
    1164. 

  1164.                     else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
    1165. 

  1165.                       precip(i) = precip(i) + qa(k)*(0.0-za(k))
    1166. 

  1166.                       exit sum_precip
    1167. 

  1167.                     endif
    1168. 

  1168.                     exit sum_precip
    1169. 

  1169.       enddo sum_precip
    1170. 

  1170. !
    1171. 

  1171. ! replace the new values
    1172. 

  1172.       rql(i,:) = qn(:)
    1173. 

  1173. !
    1174. 

  1174. ! -------------------------------…
    1175.

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