PageRenderTime 33ms CodeModel.GetById 31ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_mp_thompson.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 3629 lines | 2525 code | 368 blank | 736 comment | 188 complexity | 3d748c85f72e2afd42fc1cc9f03b6927 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. !+---+-----------------------------------------------------------------+
  2. !.. This subroutine computes the moisture tendencies of water vapor,
  3. !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
  4. !.. Prior to WRFv2.2 this code was based on Reisner et al (1998), but
  5. !.. few of those pieces remain. A complete description is now found in
  6. !.. Thompson, G., P. R. Field, R. M. Rasmussen, and W. D. Hall, 2008:
  7. !.. Explicit Forecasts of winter precipitation using an improved bulk
  8. !.. microphysics scheme. Part II: Implementation of a new snow
  9. !.. parameterization. Mon. Wea. Rev., 136, 5095-5115.
  10. !.. Prior to WRFv3.1, this code was single-moment rain prediction as
  11. !.. described in the reference above, but in v3.1 and higher, the
  12. !.. scheme is two-moment rain (predicted rain number concentration).
  13. !..
  14. !.. Most importantly, users may wish to modify the prescribed number of
  15. !.. cloud droplets (Nt_c; see guidelines mentioned below). Otherwise,
  16. !.. users may alter the rain and graupel size distribution parameters
  17. !.. to use exponential (Marshal-Palmer) or generalized gamma shape.
  18. !.. The snow field assumes a combination of two gamma functions (from
  19. !.. Field et al. 2005) and would require significant modifications
  20. !.. throughout the entire code to alter its shape as well as accretion
  21. !.. rates. Users may also alter the constants used for density of rain,
  22. !.. graupel, ice, and snow, but the latter is not constant when using
  23. !.. Paul Field's snow distribution and moments methods. Other values
  24. !.. users can modify include the constants for mass and/or velocity
  25. !.. power law relations and assumed capacitances used in deposition/
  26. !.. sublimation/evaporation/melting.
  27. !.. Remaining values should probably be left alone.
  28. !..
  29. !..Author: Greg Thompson, NCAR-RAL, gthompsn@ucar.edu, 303-497-2805
  30. !..Last modified: 20 Mar 2012
  31. !+---+-----------------------------------------------------------------+
  32. !wrft:model_layer:physics
  33. !+---+-----------------------------------------------------------------+
  34. !
  35. MODULE module_mp_thompson
  36. USE module_wrf_error
  37. ! USE module_mp_radar
  38. ! USE module_utility, ONLY: WRFU_Clock, WRFU_Alarm
  39. ! USE module_domain, ONLY : HISTORY_ALARM, Is_alarm_tstep
  40. IMPLICIT NONE
  41. LOGICAL, PARAMETER, PRIVATE:: iiwarm = .false.
  42. INTEGER, PARAMETER, PRIVATE:: IFDRY = 0
  43. REAL, PARAMETER, PRIVATE:: T_0 = 273.15
  44. REAL, PARAMETER, PRIVATE:: PI = 3.1415926536
  45. !..Densities of rain, snow, graupel, and cloud ice.
  46. REAL, PARAMETER, PRIVATE:: rho_w = 1000.0
  47. REAL, PARAMETER, PRIVATE:: rho_s = 100.0
  48. REAL, PARAMETER, PRIVATE:: rho_g = 500.0
  49. REAL, PARAMETER, PRIVATE:: rho_i = 890.0
  50. !..Prescribed number of cloud droplets. Set according to known data or
  51. !.. roughly 100 per cc (100.E6 m^-3) for Maritime cases and
  52. !.. 300 per cc (300.E6 m^-3) for Continental. Gamma shape parameter,
  53. !.. mu_c, calculated based on Nt_c is important in autoconversion
  54. !.. scheme.
  55. REAL, PARAMETER, PRIVATE:: Nt_c = 100.E6
  56. !..Generalized gamma distributions for rain, graupel and cloud ice.
  57. !.. N(D) = N_0 * D**mu * exp(-lamda*D); mu=0 is exponential.
  58. REAL, PARAMETER, PRIVATE:: mu_r = 0.0
  59. REAL, PARAMETER, PRIVATE:: mu_g = 0.0
  60. REAL, PARAMETER, PRIVATE:: mu_i = 0.0
  61. REAL, PRIVATE:: mu_c
  62. !..Sum of two gamma distrib for snow (Field et al. 2005).
  63. !.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3)
  64. !.. + Kap1*(M2/M3)**mu_s * D**mu_s * exp(-M2*Lam1*D/M3)]
  65. !.. M2 and M3 are the (bm_s)th and (bm_s+1)th moments respectively
  66. !.. calculated as function of ice water content and temperature.
  67. REAL, PARAMETER, PRIVATE:: mu_s = 0.6357
  68. REAL, PARAMETER, PRIVATE:: Kap0 = 490.6
  69. REAL, PARAMETER, PRIVATE:: Kap1 = 17.46
  70. REAL, PARAMETER, PRIVATE:: Lam0 = 20.78
  71. REAL, PARAMETER, PRIVATE:: Lam1 = 3.29
  72. !..Y-intercept parameter for graupel is not constant and depends on
  73. !.. mixing ratio. Also, when mu_g is non-zero, these become equiv
  74. !.. y-intercept for an exponential distrib and proper values are
  75. !.. computed based on same mixing ratio and total number concentration.
  76. REAL, PARAMETER, PRIVATE:: gonv_min = 2.E4
  77. REAL, PARAMETER, PRIVATE:: gonv_max = 2.E6
  78. !..Mass power law relations: mass = am*D**bm
  79. !.. Snow from Field et al. (2005), others assume spherical form.
  80. REAL, PARAMETER, PRIVATE:: am_r = PI*rho_w/6.0
  81. REAL, PARAMETER, PRIVATE:: bm_r = 3.0
  82. REAL, PARAMETER, PRIVATE:: am_s = 0.069
  83. REAL, PARAMETER, PRIVATE:: bm_s = 2.0
  84. REAL, PARAMETER, PRIVATE:: am_g = PI*rho_g/6.0
  85. REAL, PARAMETER, PRIVATE:: bm_g = 3.0
  86. REAL, PARAMETER, PRIVATE:: am_i = PI*rho_i/6.0
  87. REAL, PARAMETER, PRIVATE:: bm_i = 3.0
  88. !..Fallspeed power laws relations: v = (av*D**bv)*exp(-fv*D)
  89. !.. Rain from Ferrier (1994), ice, snow, and graupel from
  90. !.. Thompson et al (2008). Coefficient fv is zero for graupel/ice.
  91. REAL, PARAMETER, PRIVATE:: av_r = 4854.0
  92. REAL, PARAMETER, PRIVATE:: bv_r = 1.0
  93. REAL, PARAMETER, PRIVATE:: fv_r = 195.0
  94. REAL, PARAMETER, PRIVATE:: av_s = 40.0
  95. REAL, PARAMETER, PRIVATE:: bv_s = 0.55
  96. REAL, PARAMETER, PRIVATE:: fv_s = 100.0
  97. REAL, PARAMETER, PRIVATE:: av_g = 442.0
  98. REAL, PARAMETER, PRIVATE:: bv_g = 0.89
  99. REAL, PARAMETER, PRIVATE:: av_i = 1847.5
  100. REAL, PARAMETER, PRIVATE:: bv_i = 1.0
  101. !..Capacitance of sphere and plates/aggregates: D**3, D**2
  102. REAL, PARAMETER, PRIVATE:: C_cube = 0.5
  103. REAL, PARAMETER, PRIVATE:: C_sqrd = 0.3
  104. !..Collection efficiencies. Rain/snow/graupel collection of cloud
  105. !.. droplets use variables (Ef_rw, Ef_sw, Ef_gw respectively) and
  106. !.. get computed elsewhere because they are dependent on stokes
  107. !.. number.
  108. REAL, PARAMETER, PRIVATE:: Ef_si = 0.05
  109. REAL, PARAMETER, PRIVATE:: Ef_rs = 0.95
  110. REAL, PARAMETER, PRIVATE:: Ef_rg = 0.75
  111. REAL, PARAMETER, PRIVATE:: Ef_ri = 0.95
  112. !..Minimum microphys values
  113. !.. R1 value, 1.E-12, cannot be set lower because of numerical
  114. !.. problems with Paul Field's moments and should not be set larger
  115. !.. because of truncation problems in snow/ice growth.
  116. REAL, PARAMETER, PRIVATE:: R1 = 1.E-12
  117. REAL, PARAMETER, PRIVATE:: R2 = 1.E-6
  118. REAL, PARAMETER, PRIVATE:: eps = 1.E-15
  119. !..Constants in Cooper curve relation for cloud ice number.
  120. REAL, PARAMETER, PRIVATE:: TNO = 5.0
  121. REAL, PARAMETER, PRIVATE:: ATO = 0.304
  122. !..Rho_not used in fallspeed relations (rho_not/rho)**.5 adjustment.
  123. REAL, PARAMETER, PRIVATE:: rho_not = 101325.0/(287.05*298.0)
  124. !..Schmidt number
  125. REAL, PARAMETER, PRIVATE:: Sc = 0.632
  126. REAL, PRIVATE:: Sc3
  127. !..Homogeneous freezing temperature
  128. REAL, PARAMETER, PRIVATE:: HGFR = 235.16
  129. !..Water vapor and air gas constants at constant pressure
  130. REAL, PARAMETER, PRIVATE:: Rv = 461.5
  131. REAL, PARAMETER, PRIVATE:: oRv = 1./Rv
  132. REAL, PARAMETER, PRIVATE:: R = 287.04
  133. REAL, PARAMETER, PRIVATE:: Cp = 1004.0
  134. !..Enthalpy of sublimation, vaporization, and fusion at 0C.
  135. REAL, PARAMETER, PRIVATE:: lsub = 2.834E6
  136. REAL, PARAMETER, PRIVATE:: lvap0 = 2.5E6
  137. REAL, PARAMETER, PRIVATE:: lfus = lsub - lvap0
  138. REAL, PARAMETER, PRIVATE:: olfus = 1./lfus
  139. !..Ice initiates with this mass (kg), corresponding diameter calc.
  140. !..Min diameters and mass of cloud, rain, snow, and graupel (m, kg).
  141. REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12
  142. REAL, PARAMETER, PRIVATE:: D0c = 1.E-6
  143. REAL, PARAMETER, PRIVATE:: D0r = 50.E-6
  144. REAL, PARAMETER, PRIVATE:: D0s = 200.E-6
  145. REAL, PARAMETER, PRIVATE:: D0g = 250.E-6
  146. REAL, PRIVATE:: D0i, xm0s, xm0g
  147. !..Lookup table dimensions
  148. INTEGER, PARAMETER, PRIVATE:: nbins = 100
  149. INTEGER, PARAMETER, PRIVATE:: nbc = nbins
  150. INTEGER, PARAMETER, PRIVATE:: nbi = nbins
  151. INTEGER, PARAMETER, PRIVATE:: nbr = nbins
  152. INTEGER, PARAMETER, PRIVATE:: nbs = nbins
  153. INTEGER, PARAMETER, PRIVATE:: nbg = nbins
  154. INTEGER, PARAMETER, PRIVATE:: ntb_c = 37
  155. INTEGER, PARAMETER, PRIVATE:: ntb_i = 64
  156. INTEGER, PARAMETER, PRIVATE:: ntb_r = 37
  157. INTEGER, PARAMETER, PRIVATE:: ntb_s = 28
  158. INTEGER, PARAMETER, PRIVATE:: ntb_g = 28
  159. INTEGER, PARAMETER, PRIVATE:: ntb_g1 = 28
  160. INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37
  161. INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55
  162. INTEGER, PARAMETER, PRIVATE:: ntb_t = 9
  163. INTEGER, PRIVATE:: nic2, nii2, nii3, nir2, nir3, nis2, nig2, nig3
  164. DOUBLE PRECISION, DIMENSION(nbins+1):: xDx
  165. DOUBLE PRECISION, DIMENSION(nbc):: Dc, dtc
  166. DOUBLE PRECISION, DIMENSION(nbi):: Di, dti
  167. DOUBLE PRECISION, DIMENSION(nbr):: Dr, dtr
  168. DOUBLE PRECISION, DIMENSION(nbs):: Ds, dts
  169. DOUBLE PRECISION, DIMENSION(nbg):: Dg, dtg
  170. !..Lookup tables for cloud water content (kg/m**3).
  171. REAL, DIMENSION(ntb_c), PARAMETER, PRIVATE:: &
  172. r_c = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, &
  173. 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
  174. 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
  175. 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, &
  176. 1.e-2/)
  177. !..Lookup tables for cloud ice content (kg/m**3).
  178. REAL, DIMENSION(ntb_i), PARAMETER, PRIVATE:: &
  179. r_i = (/1.e-10,2.e-10,3.e-10,4.e-10, &
  180. 5.e-10,6.e-10,7.e-10,8.e-10,9.e-10, &
  181. 1.e-9,2.e-9,3.e-9,4.e-9,5.e-9,6.e-9,7.e-9,8.e-9,9.e-9, &
  182. 1.e-8,2.e-8,3.e-8,4.e-8,5.e-8,6.e-8,7.e-8,8.e-8,9.e-8, &
  183. 1.e-7,2.e-7,3.e-7,4.e-7,5.e-7,6.e-7,7.e-7,8.e-7,9.e-7, &
  184. 1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, &
  185. 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
  186. 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
  187. 1.e-3/)
  188. !..Lookup tables for rain content (kg/m**3).
  189. REAL, DIMENSION(ntb_r), PARAMETER, PRIVATE:: &
  190. r_r = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, &
  191. 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
  192. 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
  193. 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, &
  194. 1.e-2/)
  195. !..Lookup tables for graupel content (kg/m**3).
  196. REAL, DIMENSION(ntb_g), PARAMETER, PRIVATE:: &
  197. r_g = (/1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
  198. 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
  199. 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, &
  200. 1.e-2/)
  201. !..Lookup tables for snow content (kg/m**3).
  202. REAL, DIMENSION(ntb_s), PARAMETER, PRIVATE:: &
  203. r_s = (/1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
  204. 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
  205. 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, &
  206. 1.e-2/)
  207. !..Lookup tables for rain y-intercept parameter (/m**4).
  208. REAL, DIMENSION(ntb_r1), PARAMETER, PRIVATE:: &
  209. N0r_exp = (/1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, &
  210. 1.e7,2.e7,3.e7,4.e7,5.e7,6.e7,7.e7,8.e7,9.e7, &
  211. 1.e8,2.e8,3.e8,4.e8,5.e8,6.e8,7.e8,8.e8,9.e8, &
  212. 1.e9,2.e9,3.e9,4.e9,5.e9,6.e9,7.e9,8.e9,9.e9, &
  213. 1.e10/)
  214. !..Lookup tables for graupel y-intercept parameter (/m**4).
  215. REAL, DIMENSION(ntb_g1), PARAMETER, PRIVATE:: &
  216. N0g_exp = (/1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
  217. 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
  218. 1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, &
  219. 1.e7/)
  220. !..Lookup tables for ice number concentration (/m**3).
  221. REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: &
  222. Nt_i = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, &
  223. 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, &
  224. 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, &
  225. 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
  226. 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
  227. 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
  228. 1.e6/)
  229. !..For snow moments conversions (from Field et al. 2005)
  230. REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
  231. sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, &
  232. 0.31255, 0.000204, 0.003199, 0.0, -0.015952/)
  233. REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
  234. sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, &
  235. 0.060366, 0.000079, 0.000594, 0.0, -0.003577/)
  236. !..Temperatures (5 C interval 0 to -40) used in lookup tables.
  237. REAL, DIMENSION(ntb_t), PARAMETER, PRIVATE:: &
  238. Tc = (/-0.01, -5., -10., -15., -20., -25., -30., -35., -40./)
  239. !..Lookup tables for various accretion/collection terms.
  240. !.. ntb_x refers to the number of elements for rain, snow, graupel,
  241. !.. and temperature array indices. Variables beginning with t-p/c/m/n
  242. !.. represent lookup tables. Save compile-time memory by making
  243. !.. allocatable (2009Jun12, J. Michalakes).
  244. INTEGER, PARAMETER, PRIVATE:: R8SIZE = 8
  245. REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
  246. tcg_racg, tmr_racg, tcr_gacr, tmg_gacr, &
  247. tnr_racg, tnr_gacr
  248. REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
  249. tcs_racs1, tmr_racs1, tcs_racs2, tmr_racs2, &
  250. tcr_sacr1, tms_sacr1, tcr_sacr2, tms_sacr2, &
  251. tnr_racs1, tnr_racs2, tnr_sacr1, tnr_sacr2
  252. REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: &
  253. tpi_qcfz, tni_qcfz
  254. REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: &
  255. tpi_qrfz, tpg_qrfz, tni_qrfz, tnr_qrfz
  256. REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: &
  257. tps_iaus, tni_iaus, tpi_ide
  258. REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efrw
  259. REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efsw
  260. REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: tnr_rev
  261. !..Variables holding a bunch of exponents and gamma values (cloud water,
  262. !.. cloud ice, rain, snow, then graupel).
  263. REAL, DIMENSION(3), PRIVATE:: cce, ccg
  264. REAL, PRIVATE:: ocg1, ocg2
  265. REAL, DIMENSION(7), PRIVATE:: cie, cig
  266. REAL, PRIVATE:: oig1, oig2, obmi
  267. REAL, DIMENSION(13), PRIVATE:: cre, crg
  268. REAL, PRIVATE:: ore1, org1, org2, org3, obmr
  269. REAL, DIMENSION(18), PRIVATE:: cse, csg
  270. REAL, PRIVATE:: oams, obms, ocms
  271. REAL, DIMENSION(12), PRIVATE:: cge, cgg
  272. REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, oamg, obmg, ocmg
  273. !..Declaration of precomputed constants in various rate eqns.
  274. REAL:: t1_qr_qc, t1_qr_qi, t2_qr_qi, t1_qg_qc, t1_qs_qc, t1_qs_qi
  275. REAL:: t1_qr_ev, t2_qr_ev
  276. REAL:: t1_qs_sd, t2_qs_sd, t1_qg_sd, t2_qg_sd
  277. REAL:: t1_qs_me, t2_qs_me, t1_qg_me, t2_qg_me
  278. CHARACTER*256:: mp_debug
  279. !+---+
  280. !+---+-----------------------------------------------------------------+
  281. !..END DECLARATIONS
  282. !+---+-----------------------------------------------------------------+
  283. !+---+
  284. !ctrlL
  285. CONTAINS
  286. SUBROUTINE thompson_init
  287. IMPLICIT NONE
  288. INTEGER:: i, j, k, m, n
  289. LOGICAL:: micro_init
  290. !..Allocate space for lookup tables (J. Michalakes 2009Jun08).
  291. micro_init = .FALSE.
  292. if (.NOT. ALLOCATED(tcg_racg) ) then
  293. ALLOCATE(tcg_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
  294. micro_init = .TRUE.
  295. endif
  296. if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
  297. if (.NOT. ALLOCATED(tcr_gacr)) ALLOCATE(tcr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
  298. if (.NOT. ALLOCATED(tmg_gacr)) ALLOCATE(tmg_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
  299. if (.NOT. ALLOCATED(tnr_racg)) ALLOCATE(tnr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
  300. if (.NOT. ALLOCATED(tnr_gacr)) ALLOCATE(tnr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
  301. if (.NOT. ALLOCATED(tcs_racs1)) ALLOCATE(tcs_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
  302. if (.NOT. ALLOCATED(tmr_racs1)) ALLOCATE(tmr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
  303. if (.NOT. ALLOCATED(tcs_racs2)) ALLOCATE(tcs_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
  304. if (.NOT. ALLOCATED(tmr_racs2)) ALLOCATE(tmr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
  305. if (.NOT. ALLOCATED(tcr_sacr1)) ALLOCATE(tcr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
  306. if (.NOT. ALLOCATED(tms_sacr1)) ALLOCATE(tms_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
  307. if (.NOT. ALLOCATED(tcr_sacr2)) ALLOCATE(tcr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
  308. if (.NOT. ALLOCATED(tms_sacr2)) ALLOCATE(tms_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
  309. if (.NOT. ALLOCATED(tnr_racs1)) ALLOCATE(tnr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
  310. if (.NOT. ALLOCATED(tnr_racs2)) ALLOCATE(tnr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
  311. if (.NOT. ALLOCATED(tnr_sacr1)) ALLOCATE(tnr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
  312. if (.NOT. ALLOCATED(tnr_sacr2)) ALLOCATE(tnr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
  313. if (.NOT. ALLOCATED(tpi_qcfz)) ALLOCATE(tpi_qcfz(ntb_c,45))
  314. if (.NOT. ALLOCATED(tni_qcfz)) ALLOCATE(tni_qcfz(ntb_c,45))
  315. if (.NOT. ALLOCATED(tpi_qrfz)) ALLOCATE(tpi_qrfz(ntb_r,ntb_r1,45))
  316. if (.NOT. ALLOCATED(tpg_qrfz)) ALLOCATE(tpg_qrfz(ntb_r,ntb_r1,45))
  317. if (.NOT. ALLOCATED(tni_qrfz)) ALLOCATE(tni_qrfz(ntb_r,ntb_r1,45))
  318. if (.NOT. ALLOCATED(tnr_qrfz)) ALLOCATE(tnr_qrfz(ntb_r,ntb_r1,45))
  319. if (.NOT. ALLOCATED(tps_iaus)) ALLOCATE(tps_iaus(ntb_i,ntb_i1))
  320. if (.NOT. ALLOCATED(tni_iaus)) ALLOCATE(tni_iaus(ntb_i,ntb_i1))
  321. if (.NOT. ALLOCATED(tpi_ide)) ALLOCATE(tpi_ide(ntb_i,ntb_i1))
  322. if (.NOT. ALLOCATED(t_Efrw)) ALLOCATE(t_Efrw(nbr,nbc))
  323. if (.NOT. ALLOCATED(t_Efsw)) ALLOCATE(t_Efsw(nbs,nbc))
  324. if (.NOT. ALLOCATED(tnr_rev)) ALLOCATE(tnr_rev(nbr, ntb_r1, ntb_r))
  325. if (micro_init) then
  326. !..From Martin et al. (1994), assign gamma shape parameter mu for cloud
  327. !.. drops according to general dispersion characteristics (disp=~0.25
  328. !.. for Maritime and 0.45 for Continental).
  329. !.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime
  330. !.. to 2 for really dirty air.
  331. mu_c = MIN(15., (1000.E6/Nt_c + 2.))
  332. !..Schmidt number to one-third used numerous times.
  333. Sc3 = Sc**(1./3.)
  334. !..Compute min ice diam from mass, min snow/graupel mass from diam.
  335. D0i = (xm0i/am_i)**(1./bm_i)
  336. xm0s = am_s * D0s**bm_s
  337. xm0g = am_g * D0g**bm_g
  338. !..These constants various exponents and gamma() assoc with cloud,
  339. !.. rain, snow, and graupel.
  340. cce(1) = mu_c + 1.
  341. cce(2) = bm_r + mu_c + 1.
  342. cce(3) = bm_r + mu_c + 4.
  343. ccg(1) = WGAMMA(cce(1))
  344. ccg(2) = WGAMMA(cce(2))
  345. ccg(3) = WGAMMA(cce(3))
  346. ocg1 = 1./ccg(1)
  347. ocg2 = 1./ccg(2)
  348. cie(1) = mu_i + 1.
  349. cie(2) = bm_i + mu_i + 1.
  350. cie(3) = bm_i + mu_i + bv_i + 1.
  351. cie(4) = mu_i + bv_i + 1.
  352. cie(5) = mu_i + 2.
  353. cie(6) = bm_i*0.5 + mu_i + bv_i + 1.
  354. cie(7) = bm_i*0.5 + mu_i + 1.
  355. cig(1) = WGAMMA(cie(1))
  356. cig(2) = WGAMMA(cie(2))
  357. cig(3) = WGAMMA(cie(3))
  358. cig(4) = WGAMMA(cie(4))
  359. cig(5) = WGAMMA(cie(5))
  360. cig(6) = WGAMMA(cie(6))
  361. cig(7) = WGAMMA(cie(7))
  362. oig1 = 1./cig(1)
  363. oig2 = 1./cig(2)
  364. obmi = 1./bm_i
  365. cre(1) = bm_r + 1.
  366. cre(2) = mu_r + 1.
  367. cre(3) = bm_r + mu_r + 1.
  368. cre(4) = bm_r*2. + mu_r + 1.
  369. cre(5) = mu_r + bv_r + 1.
  370. cre(6) = bm_r + mu_r + bv_r + 1.
  371. cre(7) = bm_r*0.5 + mu_r + bv_r + 1.
  372. cre(8) = bm_r + mu_r + bv_r + 3.
  373. cre(9) = mu_r + bv_r + 3.
  374. cre(10) = mu_r + 2.
  375. cre(11) = 0.5*(bv_r + 5. + 2.*mu_r)
  376. cre(12) = bm_r*0.5 + mu_r + 1.
  377. cre(13) = bm_r*2. + mu_r + bv_r + 1.
  378. do n = 1, 13
  379. crg(n) = WGAMMA(cre(n))
  380. enddo
  381. obmr = 1./bm_r
  382. ore1 = 1./cre(1)
  383. org1 = 1./crg(1)
  384. org2 = 1./crg(2)
  385. org3 = 1./crg(3)
  386. cse(1) = bm_s + 1.
  387. cse(2) = bm_s + 2.
  388. cse(3) = bm_s*2.
  389. cse(4) = bm_s + bv_s + 1.
  390. cse(5) = bm_s*2. + bv_s + 1.
  391. cse(6) = bm_s*2. + 1.
  392. cse(7) = bm_s + mu_s + 1.
  393. cse(8) = bm_s + mu_s + 2.
  394. cse(9) = bm_s + mu_s + 3.
  395. cse(10) = bm_s + mu_s + bv_s + 1.
  396. cse(11) = bm_s*2. + mu_s + bv_s + 1.
  397. cse(12) = bm_s*2. + mu_s + 1.
  398. cse(13) = bv_s + 2.
  399. cse(14) = bm_s + bv_s
  400. cse(15) = mu_s + 1.
  401. cse(16) = 1.0 + (1.0 + bv_s)/2.
  402. cse(17) = cse(16) + mu_s + 1.
  403. cse(18) = bv_s + mu_s + 3.
  404. do n = 1, 18
  405. csg(n) = WGAMMA(cse(n))
  406. enddo
  407. oams = 1./am_s
  408. obms = 1./bm_s
  409. ocms = oams**obms
  410. cge(1) = bm_g + 1.
  411. cge(2) = mu_g + 1.
  412. cge(3) = bm_g + mu_g + 1.
  413. cge(4) = bm_g*2. + mu_g + 1.
  414. cge(5) = bm_g*2. + mu_g + bv_g + 1.
  415. cge(6) = bm_g + mu_g + bv_g + 1.
  416. cge(7) = bm_g + mu_g + bv_g + 2.
  417. cge(8) = bm_g + mu_g + bv_g + 3.
  418. cge(9) = mu_g + bv_g + 3.
  419. cge(10) = mu_g + 2.
  420. cge(11) = 0.5*(bv_g + 5. + 2.*mu_g)
  421. cge(12) = 0.5*(bv_g + 5.) + mu_g
  422. do n = 1, 12
  423. cgg(n) = WGAMMA(cge(n))
  424. enddo
  425. oamg = 1./am_g
  426. obmg = 1./bm_g
  427. ocmg = oamg**obmg
  428. oge1 = 1./cge(1)
  429. ogg1 = 1./cgg(1)
  430. ogg2 = 1./cgg(2)
  431. ogg3 = 1./cgg(3)
  432. !+---+-----------------------------------------------------------------+
  433. !..Simplify various rate eqns the best we can now.
  434. !+---+-----------------------------------------------------------------+
  435. !..Rain collecting cloud water and cloud ice
  436. t1_qr_qc = PI*.25*av_r * crg(9)
  437. t1_qr_qi = PI*.25*av_r * crg(9)
  438. t2_qr_qi = PI*.25*am_r*av_r * crg(8)
  439. !..Graupel collecting cloud water
  440. t1_qg_qc = PI*.25*av_g * cgg(9)
  441. !..Snow collecting cloud water
  442. t1_qs_qc = PI*.25*av_s
  443. !..Snow collecting cloud ice
  444. t1_qs_qi = PI*.25*av_s
  445. !..Evaporation of rain; ignore depositional growth of rain.
  446. t1_qr_ev = 0.78 * crg(10)
  447. t2_qr_ev = 0.308*Sc3*SQRT(av_r) * crg(11)
  448. !..Sublimation/depositional growth of snow
  449. t1_qs_sd = 0.86
  450. t2_qs_sd = 0.28*Sc3*SQRT(av_s)
  451. !..Melting of snow
  452. t1_qs_me = PI*4.*C_sqrd*olfus * 0.86
  453. t2_qs_me = PI*4.*C_sqrd*olfus * 0.28*Sc3*SQRT(av_s)
  454. !..Sublimation/depositional growth of graupel
  455. t1_qg_sd = 0.86 * cgg(10)
  456. t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11)
  457. !..Melting of graupel
  458. t1_qg_me = PI*4.*C_cube*olfus * 0.86 * cgg(10)
  459. t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11)
  460. !..Constants for helping find lookup table indexes.
  461. nic2 = NINT(ALOG10(r_c(1)))
  462. nii2 = NINT(ALOG10(r_i(1)))
  463. nii3 = NINT(ALOG10(Nt_i(1)))
  464. nir2 = NINT(ALOG10(r_r(1)))
  465. nir3 = NINT(ALOG10(N0r_exp(1)))
  466. nis2 = NINT(ALOG10(r_s(1)))
  467. nig2 = NINT(ALOG10(r_g(1)))
  468. nig3 = NINT(ALOG10(N0g_exp(1)))
  469. !..Create bins of cloud water (from min diameter up to 100 microns).
  470. Dc(1) = D0c*1.0d0
  471. dtc(1) = D0c*1.0d0
  472. do n = 2, nbc
  473. Dc(n) = Dc(n-1) + 1.0D-6
  474. dtc(n) = (Dc(n) - Dc(n-1))
  475. enddo
  476. !..Create bins of cloud ice (from min diameter up to 5x min snow size).
  477. xDx(1) = D0i*1.0d0
  478. xDx(nbi+1) = 5.0d0*D0s
  479. do n = 2, nbi
  480. xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbi) &
  481. *DLOG(xDx(nbi+1)/xDx(1)) +DLOG(xDx(1)))
  482. enddo
  483. do n = 1, nbi
  484. Di(n) = DSQRT(xDx(n)*xDx(n+1))
  485. dti(n) = xDx(n+1) - xDx(n)
  486. enddo
  487. !..Create bins of rain (from min diameter up to 5 mm).
  488. xDx(1) = D0r*1.0d0
  489. xDx(nbr+1) = 0.005d0
  490. do n = 2, nbr
  491. xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbr) &
  492. *DLOG(xDx(nbr+1)/xDx(1)) +DLOG(xDx(1)))
  493. enddo
  494. do n = 1, nbr
  495. Dr(n) = DSQRT(xDx(n)*xDx(n+1))
  496. dtr(n) = xDx(n+1) - xDx(n)
  497. enddo
  498. !..Create bins of snow (from min diameter up to 2 cm).
  499. xDx(1) = D0s*1.0d0
  500. xDx(nbs+1) = 0.02d0
  501. do n = 2, nbs
  502. xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbs) &
  503. *DLOG(xDx(nbs+1)/xDx(1)) +DLOG(xDx(1)))
  504. enddo
  505. do n = 1, nbs
  506. Ds(n) = DSQRT(xDx(n)*xDx(n+1))
  507. dts(n) = xDx(n+1) - xDx(n)
  508. enddo
  509. !..Create bins of graupel (from min diameter up to 5 cm).
  510. xDx(1) = D0g*1.0d0
  511. xDx(nbg+1) = 0.05d0
  512. do n = 2, nbg
  513. xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbg) &
  514. *DLOG(xDx(nbg+1)/xDx(1)) +DLOG(xDx(1)))
  515. enddo
  516. do n = 1, nbg
  517. Dg(n) = DSQRT(xDx(n)*xDx(n+1))
  518. dtg(n) = xDx(n+1) - xDx(n)
  519. enddo
  520. !+---+-----------------------------------------------------------------+
  521. !..Create lookup tables for most costly calculations.
  522. !+---+-----------------------------------------------------------------+
  523. do m = 1, ntb_r
  524. do k = 1, ntb_r1
  525. do j = 1, ntb_g
  526. do i = 1, ntb_g1
  527. tcg_racg(i,j,k,m) = 0.0d0
  528. tmr_racg(i,j,k,m) = 0.0d0
  529. tcr_gacr(i,j,k,m) = 0.0d0
  530. tmg_gacr(i,j,k,m) = 0.0d0
  531. tnr_racg(i,j,k,m) = 0.0d0
  532. tnr_gacr(i,j,k,m) = 0.0d0
  533. enddo
  534. enddo
  535. enddo
  536. enddo
  537. do m = 1, ntb_r
  538. do k = 1, ntb_r1
  539. do j = 1, ntb_t
  540. do i = 1, ntb_s
  541. tcs_racs1(i,j,k,m) = 0.0d0
  542. tmr_racs1(i,j,k,m) = 0.0d0
  543. tcs_racs2(i,j,k,m) = 0.0d0
  544. tmr_racs2(i,j,k,m) = 0.0d0
  545. tcr_sacr1(i,j,k,m) = 0.0d0
  546. tms_sacr1(i,j,k,m) = 0.0d0
  547. tcr_sacr2(i,j,k,m) = 0.0d0
  548. tms_sacr2(i,j,k,m) = 0.0d0
  549. tnr_racs1(i,j,k,m) = 0.0d0
  550. tnr_racs2(i,j,k,m) = 0.0d0
  551. tnr_sacr1(i,j,k,m) = 0.0d0
  552. tnr_sacr2(i,j,k,m) = 0.0d0
  553. enddo
  554. enddo
  555. enddo
  556. enddo
  557. do k = 1, 45
  558. do j = 1, ntb_r1
  559. do i = 1, ntb_r
  560. tpi_qrfz(i,j,k) = 0.0d0
  561. tni_qrfz(i,j,k) = 0.0d0
  562. tpg_qrfz(i,j,k) = 0.0d0
  563. tnr_qrfz(i,j,k) = 0.0d0
  564. enddo
  565. enddo
  566. do i = 1, ntb_c
  567. tpi_qcfz(i,k) = 0.0d0
  568. tni_qcfz(i,k) = 0.0d0
  569. enddo
  570. enddo
  571. do j = 1, ntb_i1
  572. do i = 1, ntb_i
  573. tps_iaus(i,j) = 0.0d0
  574. tni_iaus(i,j) = 0.0d0
  575. tpi_ide(i,j) = 0.0d0
  576. enddo
  577. enddo
  578. do j = 1, nbc
  579. do i = 1, nbr
  580. t_Efrw(i,j) = 0.0
  581. enddo
  582. do i = 1, nbs
  583. t_Efsw(i,j) = 0.0
  584. enddo
  585. enddo
  586. do k = 1, ntb_r
  587. do j = 1, ntb_r1
  588. do i = 1, nbr
  589. tnr_rev(i,j,k) = 0.0d0
  590. enddo
  591. enddo
  592. enddo
  593. CALL wrf_debug(150, 'CREATING MICROPHYSICS LOOKUP TABLES ... ')
  594. WRITE (wrf_err_message, '(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') &
  595. ' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g
  596. CALL wrf_debug(150, wrf_err_message)
  597. !..Collision efficiency between rain/snow and cloud water.
  598. CALL wrf_debug(200, ' creating qc collision eff tables')
  599. call table_Efrw
  600. call table_Efsw
  601. !..Drop evaporation.
  602. ! CALL wrf_debug(200, ' creating rain evap table')
  603. ! call table_dropEvap
  604. !..Initialize various constants for computing radar reflectivity.
  605. ! xam_r = am_r
  606. ! xbm_r = bm_r
  607. ! xmu_r = mu_r
  608. ! xam_s = am_s
  609. ! xbm_s = bm_s
  610. ! xmu_s = mu_s
  611. ! xam_g = am_g
  612. ! xbm_g = bm_g
  613. ! xmu_g = mu_g
  614. ! call radar_init
  615. if (.not. iiwarm) then
  616. !..Rain collecting graupel & graupel collecting rain.
  617. CALL wrf_debug(200, ' creating rain collecting graupel table')
  618. call qr_acr_qg
  619. !..Rain collecting snow & snow collecting rain.
  620. CALL wrf_debug(200, ' creating rain collecting snow table')
  621. call qr_acr_qs
  622. !..Cloud water and rain freezing (Bigg, 1953).
  623. CALL wrf_debug(200, ' creating freezing of water drops table')
  624. call freezeH2O
  625. !..Conversion of some ice mass into snow category.
  626. CALL wrf_debug(200, ' creating ice converting to snow table')
  627. call qi_aut_qs
  628. endif
  629. CALL wrf_debug(150, ' ... DONE microphysical lookup tables')
  630. endif
  631. END SUBROUTINE thompson_init
  632. !+---+-----------------------------------------------------------------+
  633. !ctrlL
  634. !+---+-----------------------------------------------------------------+
  635. !..This is a wrapper routine designed to transfer values from 3D to 1D.
  636. !+---+-----------------------------------------------------------------+
  637. SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, &
  638. th, pii, p, dz, dt_in, itimestep, &
  639. RAINNC, RAINNCV, &
  640. SNOWNC, SNOWNCV, &
  641. GRAUPELNC, GRAUPELNCV, SR, &
  642. #ifdef WRF_CHEM
  643. rainprod, evapprod, &
  644. #endif
  645. ! refl_10cm, grid_clock, grid_alarms, &
  646. ids,ide, jds,jde, kds,kde, & ! domain dims
  647. ims,ime, jms,jme, kms,kme, & ! memory dims
  648. its,ite, jts,jte, kts,kte) ! tile dims
  649. implicit none
  650. !..Subroutine arguments
  651. INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, &
  652. ims,ime, jms,jme, kms,kme, &
  653. its,ite, jts,jte, kts,kte
  654. REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
  655. qv, qc, qr, qi, qs, qg, ni, nr, th
  656. #ifdef WRF_CHEM
  657. REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
  658. rainprod, evapprod
  659. #endif
  660. REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
  661. pii, p, dz
  662. REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
  663. RAINNC, RAINNCV, SR
  664. REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(INOUT):: &
  665. SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV
  666. ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
  667. ! refl_10cm
  668. REAL, INTENT(IN):: dt_in
  669. INTEGER, INTENT(IN):: itimestep
  670. ! TYPE (WRFU_Clock):: grid_clock
  671. ! TYPE (WRFU_Alarm), POINTER:: grid_alarms(:)
  672. !..Local variables
  673. REAL, DIMENSION(kts:kte):: &
  674. qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
  675. nr1d, t1d, p1d, dz1d, dBZ
  676. #ifdef WRF_CHEM
  677. REAL, DIMENSION(kts:kte):: &
  678. rainprod1d, evapprod1d
  679. #endif
  680. REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
  681. REAL:: dt, pptrain, pptsnow, pptgraul, pptice
  682. REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
  683. INTEGER:: i, j, k
  684. INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr
  685. INTEGER:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr
  686. INTEGER:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr
  687. INTEGER:: i_start, j_start, i_end, j_end
  688. LOGICAL:: dBZ_tstep
  689. !+---+
  690. dBZ_tstep = .false.
  691. ! if ( Is_alarm_tstep(grid_clock, grid_alarms(HISTORY_ALARM)) ) then
  692. ! dBZ_tstep = .true.
  693. ! endif
  694. i_start = its
  695. j_start = jts
  696. i_end = ite
  697. j_end = jte
  698. !..For idealized testing by developer.
  699. ! if ( (ide-ids+1).gt.4 .and. (jde-jds+1).lt.4 .and. &
  700. ! ids.eq.its.and.ide.eq.ite.and.jds.eq.jts.and.jde.eq.jte) then
  701. ! i_start = its + 2
  702. ! i_end = ite
  703. ! j_start = jts
  704. ! j_end = jte
  705. ! endif
  706. dt = dt_in
  707. qc_max = 0.
  708. qr_max = 0.
  709. qs_max = 0.
  710. qi_max = 0.
  711. qg_max = 0
  712. ni_max = 0.
  713. nr_max = 0.
  714. imax_qc = 0
  715. imax_qr = 0
  716. imax_qi = 0
  717. imax_qs = 0
  718. imax_qg = 0
  719. imax_ni = 0
  720. imax_nr = 0
  721. jmax_qc = 0
  722. jmax_qr = 0
  723. jmax_qi = 0
  724. jmax_qs = 0
  725. jmax_qg = 0
  726. jmax_ni = 0
  727. jmax_nr = 0
  728. kmax_qc = 0
  729. kmax_qr = 0
  730. kmax_qi = 0
  731. kmax_qs = 0
  732. kmax_qg = 0
  733. kmax_ni = 0
  734. kmax_nr = 0
  735. do i = 1, 256
  736. mp_debug(i:i) = char(0)
  737. enddo
  738. j_loop: do j = j_start, j_end
  739. i_loop: do i = i_start, i_end
  740. pptrain = 0.
  741. pptsnow = 0.
  742. pptgraul = 0.
  743. pptice = 0.
  744. RAINNCV(i,j) = 0.
  745. IF ( PRESENT (snowncv) ) THEN
  746. SNOWNCV(i,j) = 0.
  747. ENDIF
  748. IF ( PRESENT (graupelncv) ) THEN
  749. GRAUPELNCV(i,j) = 0.
  750. ENDIF
  751. SR(i,j) = 0.
  752. do k = kts, kte
  753. t1d(k) = th(i,k,j)*pii(i,k,j)
  754. p1d(k) = p(i,k,j)
  755. dz1d(k) = dz(i,k,j)
  756. qv1d(k) = qv(i,k,j)
  757. qc1d(k) = qc(i,k,j)
  758. qi1d(k) = qi(i,k,j)
  759. qr1d(k) = qr(i,k,j)
  760. qs1d(k) = qs(i,k,j)
  761. qg1d(k) = qg(i,k,j)
  762. ni1d(k) = ni(i,k,j)
  763. nr1d(k) = nr(i,k,j)
  764. enddo
  765. call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
  766. nr1d, t1d, p1d, dz1d, &
  767. pptrain, pptsnow, pptgraul, pptice, &
  768. #ifdef WRF_CHEM
  769. rainprod1d, evapprod1d, &
  770. #endif
  771. kts, kte, dt, i, j)
  772. pcp_ra(i,j) = pptrain
  773. pcp_sn(i,j) = pptsnow
  774. pcp_gr(i,j) = pptgraul
  775. pcp_ic(i,j) = pptice
  776. RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice
  777. RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
  778. IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN
  779. SNOWNCV(i,j) = pptsnow + pptice
  780. SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice
  781. ENDIF
  782. IF ( PRESENT(graupelncv) .AND. PRESENT(graupelnc) ) THEN
  783. GRAUPELNCV(i,j) = pptgraul
  784. GRAUPELNC(i,j) = GRAUPELNC(i,j) + pptgraul
  785. ENDIF
  786. SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12)
  787. do k = kts, kte
  788. qv(i,k,j) = qv1d(k)
  789. qc(i,k,j) = qc1d(k)
  790. qi(i,k,j) = qi1d(k)
  791. qr(i,k,j) = qr1d(k)
  792. qs(i,k,j) = qs1d(k)
  793. qg(i,k,j) = qg1d(k)
  794. ni(i,k,j) = ni1d(k)
  795. nr(i,k,j) = nr1d(k)
  796. #ifdef WRF_CHEM
  797. rainprod(i,k,j) = rainprod1d(k)
  798. evapprod(i,k,j) = evapprod1d(k)
  799. #endif
  800. th(i,k,j) = t1d(k)/pii(i,k,j)
  801. if (qc1d(k) .gt. qc_max) then
  802. imax_qc = i
  803. jmax_qc = j
  804. kmax_qc = k
  805. qc_max = qc1d(k)
  806. elseif (qc1d(k) .lt. 0.0) then
  807. write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), &
  808. ' at i,j,k=', i,j,k
  809. CALL wrf_debug(150, mp_debug)
  810. endif
  811. if (qr1d(k) .gt. qr_max) then
  812. imax_qr = i
  813. jmax_qr = j
  814. kmax_qr = k
  815. qr_max = qr1d(k)
  816. elseif (qr1d(k) .lt. 0.0) then
  817. write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), &
  818. ' at i,j,k=', i,j,k
  819. CALL wrf_debug(150, mp_debug)
  820. endif
  821. if (nr1d(k) .gt. nr_max) then
  822. imax_nr = i
  823. jmax_nr = j
  824. kmax_nr = k
  825. nr_max = nr1d(k)
  826. elseif (nr1d(k) .lt. 0.0) then
  827. write(mp_debug,*) 'WARNING, negative nr ', nr1d(k), &
  828. ' at i,j,k=', i,j,k
  829. CALL wrf_debug(150, mp_debug)
  830. endif
  831. if (qs1d(k) .gt. qs_max) then
  832. imax_qs = i
  833. jmax_qs = j
  834. kmax_qs = k
  835. qs_max = qs1d(k)
  836. elseif (qs1d(k) .lt. 0.0) then
  837. write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), &
  838. ' at i,j,k=', i,j,k
  839. CALL wrf_debug(150, mp_debug)
  840. endif
  841. if (qi1d(k) .gt. qi_max) then
  842. imax_qi = i
  843. jmax_qi = j
  844. kmax_qi = k
  845. qi_max = qi1d(k)
  846. elseif (qi1d(k) .lt. 0.0) then
  847. write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), &
  848. ' at i,j,k=', i,j,k
  849. CALL wrf_debug(150, mp_debug)
  850. endif
  851. if (qg1d(k) .gt. qg_max) then
  852. imax_qg = i
  853. jmax_qg = j
  854. kmax_qg = k
  855. qg_max = qg1d(k)
  856. elseif (qg1d(k) .lt. 0.0) then
  857. write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), &
  858. ' at i,j,k=', i,j,k
  859. CALL wrf_debug(150, mp_debug)
  860. endif
  861. if (ni1d(k) .gt. ni_max) then
  862. imax_ni = i
  863. jmax_ni = j
  864. kmax_ni = k
  865. ni_max = ni1d(k)
  866. elseif (ni1d(k) .lt. 0.0) then
  867. write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), &
  868. ' at i,j,k=', i,j,k
  869. CALL wrf_debug(150, mp_debug)
  870. endif
  871. if (qv1d(k) .lt. 0.0) then
  872. if (k.lt.kte-2 .and. k.gt.kts+1) then
  873. qv(i,k,j) = 0.5*(qv(i,k-1,j) + qv(i,k+1,j))
  874. else
  875. qv(i,k,j) = 1.E-7
  876. endif
  877. write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), &
  878. ' at i,j,k=', i,j,k
  879. CALL wrf_debug(150, mp_debug)
  880. endif
  881. enddo
  882. ! if (dBZ_tstep) then
  883. ! call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
  884. ! t1d, p1d, dBZ, kts, kte, i, j)
  885. ! do k = kts, kte
  886. ! refl_10cm(i,k,j) = MAX(-35., dBZ(k))
  887. ! enddo
  888. ! endif
  889. enddo i_loop
  890. enddo j_loop
  891. ! DEBUG - GT
  892. write(mp_debug,'(a,7(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', &
  893. 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', &
  894. 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', &
  895. 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', &
  896. 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', &
  897. 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', &
  898. 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', &
  899. 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')'
  900. CALL wrf_debug(150, mp_debug)
  901. ! END DEBUG - GT
  902. do i = 1, 256
  903. mp_debug(i:i) = char(0)
  904. enddo
  905. END SUBROUTINE mp_gt_driver
  906. !+---+-----------------------------------------------------------------+
  907. !ctrlL
  908. !+---+-----------------------------------------------------------------+
  909. !+---+-----------------------------------------------------------------+
  910. !.. This subroutine computes the moisture tendencies of water vapor,
  911. !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
  912. !.. Previously this code was based on Reisner et al (1998), but few of
  913. !.. those pieces remain. A complete description is now found in
  914. !.. Thompson et al. (2004, 2008).
  915. !+---+-----------------------------------------------------------------+
  916. !
  917. subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
  918. nr1d, t1d, p1d, dzq, &
  919. pptrain, pptsnow, pptgraul, pptice, &
  920. #ifdef WRF_CHEM
  921. rainprod, evapprod, &
  922. #endif
  923. kts, kte, dt, ii, jj)
  924. implicit none
  925. !..Sub arguments
  926. INTEGER, INTENT(IN):: kts, kte, ii, jj
  927. REAL, DIMENSION(kts:kte), INTENT(INOUT):: &
  928. qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
  929. nr1d, t1d, p1d
  930. #ifdef WRF_CHEM
  931. REAL, DIMENSION(kts:kte), INTENT(INOUT):: &
  932. rainprod, evapprod
  933. #endif
  934. REAL, DIMENSION(kts:kte), INTENT(IN):: dzq
  935. REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice
  936. REAL, INTENT(IN):: dt
  937. !..Local variables
  938. REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, &
  939. qrten, qsten, qgten, niten, nrten
  940. DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd
  941. DOUBLE PRECISION, DIMENSION(kts:kte):: prr_wau, prr_rcw, prr_rcs, &
  942. prr_rcg, prr_sml, prr_gml, &
  943. prr_rci, prv_rev, &
  944. pnr_wau, pnr_rcs, pnr_rcg, &
  945. pnr_rci, pnr_sml, pnr_gml, &
  946. pnr_rev, pnr_rcr, pnr_rfz
  947. DOUBLE PRECISION, DIMENSION(kts:kte):: pri_inu, pni_inu, pri_ihm, &
  948. pni_ihm, pri_wfz, pni_wfz, &
  949. pri_rfz, pni_rfz, pri_ide, &
  950. pni_ide, pri_rci, pni_rci, &
  951. pni_sci, pni_iau
  952. DOUBLE PRECISION, DIMENSION(kts:kte):: prs_iau, prs_sci, prs_rcs, &
  953. prs_scw, prs_sde, prs_ihm, &
  954. prs_ide
  955. DOUBLE PRECISION, DIMENSION(kts:kte):: prg_scw, prg_rfz, prg_gde, &
  956. prg_gcw, prg_rci, prg_rcs, &
  957. prg_rcg, prg_ihm
  958. DOUBLE PRECISION, PARAMETER:: zeroD0 = 0.0d0
  959. REAL, DIMENSION(kts:kte):: temp, pres, qv
  960. REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr
  961. REAL, DIMENSION(kts:kte):: rho, rhof, rhof2
  962. REAL, DIMENSION(kts:kte):: qvs, qvsi
  963. REAL, DIMENSION(kts:kte):: satw, sati, ssatw, ssati
  964. REAL, DIMENSION(kts:kte):: diffu, visco, vsc2, &
  965. tcond, lvap, ocp, lvt2
  966. DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g
  967. REAL, DIMENSION(kts:kte):: mvd_r, mvd_c
  968. REAL, DIMENSION(kts:kte):: smob, smo2, smo1, smo0, &
  969. smoc, smod, smoe, smof
  970. REAL, DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n
  971. REAL:: rgvm, delta_tp, orho, lfus2
  972. REAL, DIMENSION(4):: onstep
  973. DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg
  974. DOUBLE PRECISION:: lami, ilami
  975. REAL:: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m
  976. DOUBLE PRECISION:: Dr_star
  977. REAL:: zeta1, zeta, taud, tau
  978. REAL:: stoke_r, stoke_s, stoke_g, stoke_i
  979. REAL:: vti, vtr, vts, vtg
  980. REAL, DIMENSION(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk
  981. REAL, DIMENSION(kts:kte):: vts_boost
  982. REAL:: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow
  983. REAL:: a_, b_, loga_, A1, A2, tf
  984. REAL:: tempc, tc0, r_mvd1, r_mvd2, xkrat
  985. REAL:: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr
  986. REAL:: xsat, rate_max, sump, ratio
  987. REAL:: clap, fcd, dfcd
  988. REAL:: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl
  989. REAL:: r_frac, g_frac
  990. REAL:: Ef_rw, Ef_sw, Ef_gw, Ef_rr
  991. REAL:: dtsave, odts, odt, odzq
  992. REAL:: xslw1, ygra1, zans1
  993. INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq
  994. INTEGER, DIMENSION(4):: ksed1
  995. INTEGER:: nir, nis, nig, nii, nic
  996. INTEGER:: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, &
  997. idx_i1, idx_i, idx_c, idx, idx_d
  998. LOGICAL:: melti, no_micro
  999. LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg
  1000. LOGICAL:: debug_flag
  1001. !+---+
  1002. debug_flag = .false.
  1003. ! if (ii.eq.315 .and. jj.eq.2) debug_flag = .true.
  1004. no_micro = .true.
  1005. dtsave = dt
  1006. odt = 1./dt
  1007. odts = 1./dtsave
  1008. iexfrq = 1
  1009. !+---+-----------------------------------------------------------------+
  1010. !.. Source/sink terms. First 2 chars: "pr" represents source/sink of
  1011. !.. mass while "pn" represents source/sink of number. Next char is one
  1012. !.. of "v" for water vapor, "r" for rain, "i" for cloud ice, "w" for
  1013. !.. cloud water, "s" for snow, and "g" for graupel. Next chars
  1014. !.. represent processes: "de" for sublimation/deposition, "ev" for
  1015. !.. evaporation, "fz" for freezing, "ml" for melting, "au" for
  1016. !.. autoconversion, "nu" for ice nucleation, "hm" for Hallet/Mossop
  1017. !.. secondary ice production, and "c" for collection followed by the
  1018. !.. character for the species being collected. ALL of these terms are
  1019. !.. positive (except for deposition/sublimation terms which can switch
  1020. !.. signs based on super/subsaturation) and are treated as negatives
  1021. !.. where necessary in the tendency equations.
  1022. !+---+-----------------------------------------------------------------+
  1023. do k = kts, kte
  1024. tten(k) = 0.
  1025. qvten(k) = 0.
  1026. qcten(k) = 0.
  1027. qiten(k) = 0.
  1028. qrten(k) = 0.
  1029. qsten(k) = 0.
  1030. qgten(k) = 0.
  1031. niten(k) = 0.
  1032. nrten(k) = 0.
  1033. prw_vcd(k) = 0.
  1034. prv_rev(k) = 0.
  1035. prr_wau(k) = 0.
  1036. prr_rcw(k) = 0.
  1037. prr_rcs(k) = 0.
  1038. prr_rcg(k) = 0.
  1039. prr_sml(k) = 0.
  1040. prr_gml(k) = 0.
  1041. prr_rci(k) = 0.
  1042. pnr_wau(k) = 0.
  1043. pnr_rcs(k) = 0.
  1044. pnr_rcg(k) = 0.
  1045. pnr_rci(k) = 0.
  1046. pnr_sml(k) = 0.
  1047. pnr_gml(k) = 0.
  1048. pnr_rev(k) = 0.
  1049. pnr_rcr(k) = 0.
  1050. pnr_rfz(k) = 0.
  1051. pri_inu(k) = 0.
  1052. pni_inu(k) = 0.
  1053. pri_ihm(k) = 0.
  1054. pni_ihm(k) = 0.
  1055. pri_wfz(k) = 0.
  1056. pni_wfz(k) = 0.
  1057. pri_rfz(k) = 0.
  1058. pni_rfz(k) = 0.
  1059. pri_ide(k) = 0.
  1060. pni_ide(k) = 0.
  1061. pri_rci(k) = 0.
  1062. pni_rci(k) = 0.
  1063. pni_sci(k) = 0.
  1064. pni_iau(k) = 0.
  1065. prs_iau(k) = 0.
  1066. prs_sci(k) = 0.
  1067. prs_rcs(k) = 0.
  1068. prs_scw(k) = 0.
  1069. prs_sde(k) = 0.
  1070. prs_ihm(k) = 0.
  1071. prs_ide(k) = 0.
  1072. prg_scw(k) = 0.
  1073. prg_rfz(k) = 0.
  1074. prg_gde(k) = 0.
  1075. prg_gcw(k) = 0.
  1076. prg_rci(k) = 0.
  1077. prg_rcs(k) = 0.
  1078. prg_rcg(k) = 0.
  1079. prg_ihm(k) = 0.
  1080. enddo
  1081. #ifdef WRF_CHEM
  1082. do k = kts, kte
  1083. rainprod(k) = 0.
  1084. evapprod(k) = 0.
  1085. enddo
  1086. #endif
  1087. !+---+-----------------------------------------------------------------+
  1088. !..Put column of data into local arrays.
  1089. !+---+-----------------------------------------------------------------+
  1090. do k = kts, kte
  1091. temp(k) = t1d(k)
  1092. qv(k) = MAX(1.E-10, qv1d(k))
  1093. pres(k) = p1d(k)
  1094. rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
  1095. if (qc1d(k) .gt. R1) then
  1096. no_micro = .false.
  1097. rc(k) = qc1d(k)*rho(k)
  1098. L_qc(k) = .true.
  1099. else
  1100. qc1d(k) = 0.0
  1101. rc(k) = R1
  1102. L_qc(k) = .false.
  1103. endif
  1104. if (qi1d(k) .gt. R1) then
  1105. no_micro = .false.
  1106. ri(k) = qi1d(k)*rho(k)
  1107. ni(k) = MAX(R2, ni1d(k)*rho(k))
  1108. L_qi(k) = .true.
  1109. lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
  1110. ilami = 1./lami
  1111. xDi = (bm_i + mu_i + 1.) * ilami
  1112. if (xDi.lt. 20.E-6) then
  1113. lami = cie(2)/20.E-6
  1114. ni(k) = MIN(250.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
  1115. elseif (xDi.gt. 300.E-6) then
  1116. lami = cie(2)/300.E-6
  1117. ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i
  1118. endif
  1119. else
  1120. qi1d(k) = 0.0
  1121. ni1d(k) = 0.0
  1122. ri(k) = R1
  1123. ni(k) = R2
  1124. L_qi(k) = .false.
  1125. endif
  1126. if (qr1d(k) .gt. R1) then
  1127. no_micro = .false.
  1128. rr(k) = qr1d(k)*rho(k)
  1129. nr(k) = MAX(R2, nr1d(k)*rho(k))
  1130. L_qr(k) = .true.
  1131. lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
  1132. mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
  1133. if (mvd_r(k) .gt. 2.5E-3) then
  1134. mvd_r(k) = 2.5E-3
  1135. lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
  1136. nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
  1137. elseif (mvd_r(k) .lt. D0r*0.75) then
  1138. mvd_r(k) = D0r*0.75
  1139. lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
  1140. nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
  1141. endif
  1142. else
  1143. qr1d(k) = 0.0
  1144. nr1d(k) = 0.0
  1145. rr(k) = R1
  1146. nr(k) = R2
  1147. L_qr(k) = .false.
  1148. endif
  1149. if (qs1d(k) .gt. R1) then
  1150. no_micro = .false.
  1151. rs(k) = qs1d(k)*rho(k)
  1152. L_qs(k) = .true.
  1153. else
  1154. qs1d(k) = 0.0
  1155. rs(k) = R1
  1156. L_qs(k) = .false.
  1157. endif
  1158. if (qg1d(k) .gt. R1) then
  1159. no_micro = .false.
  1160. rg(k) = qg1d(k)*rho(k)
  1161. L_qg(k) = .true.
  1162. else
  1163. qg1d(k) = 0.0
  1164. rg(k) = R1
  1165. L_qg(k) = .false.
  1166. endif
  1167. enddo
  1168. !+---+-----------------------------------------------------------------+
  1169. !..Derive various thermodynamic variables frequently used.
  1170. !.. Saturation vapor pressure (mixing ratio) over liquid/ice comes from
  1171. !.. Flatau et al. 1992; enthalpy (latent heat) of vaporization from
  1172. !.. Bohren & Albrecht 1998; others from Pruppacher & Klett 1978.
  1173. !+---+-----------------------------------------------------------------+
  1174. do k = kts, kte
  1175. tempc = temp(k) - 273.15
  1176. rhof(k) = SQRT(RHO_NOT/rho(k))
  1177. rhof2(k) = SQRT(rhof(k))
  1178. qvs(k) = rslf(pres(k), temp(k))
  1179. if (tempc .le. 0.0) then
  1180. qvsi(k) = rsif(pres(k), temp(k))
  1181. else
  1182. qvsi(k) = qvs(k)
  1183. endif
  1184. satw(k) = qv(k)/qvs(k)
  1185. sati(k) = qv(k)/qvsi(k)
  1186. ssatw(k) = satw(k) - 1.
  1187. ssati(k) = sati(k) - 1.
  1188. if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
  1189. if (abs(ssati(k)).lt. eps) ssati(k) = 0.0
  1190. if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false.
  1191. diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
  1192. if (tempc .ge. 0.0) then
  1193. visco(k) = (1.718+0.0049*tempc)*1.0E-5
  1194. else
  1195. visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
  1196. endif
  1197. ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
  1198. vsc2(k) = SQRT(rho(k)/visco(k))
  1199. lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
  1200. tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
  1201. enddo
  1202. !+---+-----------------------------------------------------------------+
  1203. !..If no existing hydrometeor species and no chance to initiate ice or
  1204. !.. condense cloud water, just exit quickly!
  1205. !+---+-----------------------------------------------------------------+
  1206. if (no_micro) return
  1207. !+---+-----------------------------------------------------------------+
  1208. !..Calculate y-intercept, slope, and useful moments for snow.
  1209. !+---+-------

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