Line data Source code
1 : module micro_mg1_0
2 :
3 : !---------------------------------------------------------------------------------
4 : ! Purpose:
5 : ! MG microphysics
6 : !
7 : ! Author: Andrew Gettelman, Hugh Morrison.
8 : ! Contributions from: Xiaohong Liu and Steve Ghan
9 : ! December 2005-May 2010
10 : ! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
11 : ! Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)
12 : ! for questions contact Hugh Morrison, Andrew Gettelman
13 : ! e-mail: morrison@ucar.edu, andrew@ucar.edu
14 : !
15 : ! NOTE: Modified to allow other microphysics packages (e.g. CARMA) to do ice
16 : ! microphysics in cooperation with the MG liquid microphysics. This is
17 : ! controlled by the do_cldice variable.
18 : !
19 : ! NOTE: If do_cldice is false, then MG microphysics should not update CLDICE
20 : ! or NUMICE; however, it is assumed that the other microphysics scheme will have
21 : ! updated CLDICE and NUMICE. The other microphysics should handle the following
22 : ! processes that would have been done by MG:
23 : ! - Detrainment (liquid and ice)
24 : ! - Homogeneous ice nucleation
25 : ! - Heterogeneous ice nucleation
26 : ! - Bergeron process
27 : ! - Melting of ice
28 : ! - Freezing of cloud drops
29 : ! - Autoconversion (ice -> snow)
30 : ! - Growth/Sublimation of ice
31 : ! - Sedimentation of ice
32 : !---------------------------------------------------------------------------------
33 : ! modification for sub-columns, HM, (orig 8/11/10)
34 : ! This is done using the logical 'microp_uniform' set to .true. = uniform for subcolumns
35 : !---------------------------------------------------------------------------------
36 :
37 : ! Procedures required:
38 : ! 1) An implementation of the gamma function (if not intrinsic).
39 : ! 2) saturation vapor pressure to specific humidity formula
40 : ! 3) svp over water
41 : ! 4) svp over ice
42 :
43 : #ifndef HAVE_GAMMA_INTRINSICS
44 : use shr_spfn_mod, only: gamma => shr_spfn_gamma
45 : #endif
46 :
47 : use wv_sat_methods, only: &
48 : svp_water => wv_sat_svp_water, &
49 : svp_ice => wv_sat_svp_ice, &
50 : svp_to_qsat => wv_sat_svp_to_qsat
51 :
52 : use phys_control, only: phys_getopts
53 :
54 : implicit none
55 : private
56 : save
57 :
58 : ! Note: The liu_in option has been removed, as there was a serious bug with this
59 : ! option being set to false. The code now behaves as if the default liu_in=.true.
60 : ! is always on. Addition/reinstatement of ice nucleation options will likely be
61 : ! done outside of this module.
62 :
63 : public :: &
64 : micro_mg_init, &
65 : micro_mg_get_cols, &
66 : micro_mg_tend
67 :
68 : integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
69 :
70 : real(r8) :: g !gravity
71 : real(r8) :: r !Dry air Gas constant
72 : real(r8) :: rv !water vapor gas contstant
73 : real(r8) :: cpp !specific heat of dry air
74 : real(r8) :: rhow !density of liquid water
75 : real(r8) :: tmelt ! Freezing point of water (K)
76 : real(r8) :: xxlv ! latent heat of vaporization
77 : real(r8) :: xlf !latent heat of freezing
78 : real(r8) :: xxls !latent heat of sublimation
79 :
80 : real(r8) :: rhosn ! bulk density snow
81 : real(r8) :: rhoi ! bulk density ice
82 :
83 : real(r8) :: ac,bc,as,bs,ai,bi,ar,br !fall speed parameters
84 : real(r8) :: ci,di !ice mass-diameter relation parameters
85 : real(r8) :: cs,ds !snow mass-diameter relation parameters
86 : real(r8) :: cr,dr !drop mass-diameter relation parameters
87 : real(r8) :: f1s,f2s !ventilation param for snow
88 : real(r8) :: Eii !collection efficiency aggregation of ice
89 : real(r8) :: Ecr !collection efficiency cloud droplets/rain
90 : real(r8) :: f1r,f2r !ventilation param for rain
91 : real(r8) :: DCS !autoconversion size threshold
92 : real(r8) :: qsmall !min mixing ratio
93 : real(r8) :: bimm,aimm !immersion freezing
94 : real(r8) :: rhosu !typical 850mn air density
95 : real(r8) :: mi0 ! new crystal mass
96 : real(r8) :: rin ! radius of contact nuclei
97 : real(r8) :: pi ! pi
98 :
99 : ! Additional constants to help speed up code
100 :
101 : real(r8) :: cons1
102 : real(r8) :: cons4
103 : real(r8) :: cons5
104 : real(r8) :: cons6
105 : real(r8) :: cons7
106 : real(r8) :: cons8
107 : real(r8) :: cons11
108 : real(r8) :: cons13
109 : real(r8) :: cons14
110 : real(r8) :: cons16
111 : real(r8) :: cons17
112 : real(r8) :: cons22
113 : real(r8) :: cons23
114 : real(r8) :: cons24
115 : real(r8) :: cons25
116 : real(r8) :: cons27
117 : real(r8) :: cons28
118 :
119 : real(r8) :: lammini
120 : real(r8) :: lammaxi
121 : real(r8) :: lamminr
122 : real(r8) :: lammaxr
123 : real(r8) :: lammins
124 : real(r8) :: lammaxs
125 :
126 : ! parameters for snow/rain fraction for convective clouds
127 : real(r8) :: tmax_fsnow ! max temperature for transition to convective snow
128 : real(r8) :: tmin_fsnow ! min temperature for transition to convective snow
129 :
130 : !needed for findsp
131 : real(r8) :: tt0 ! Freezing temperature
132 :
133 : real(r8) :: csmin,csmax,minrefl,mindbz
134 :
135 : real(r8) :: rhmini ! Minimum rh for ice cloud fraction > 0.
136 :
137 : logical :: use_hetfrz_classnuc ! option to use heterogeneous freezing
138 :
139 : character(len=16) :: micro_mg_precip_frac_method ! type of precipitation fraction method
140 : real(r8) :: micro_mg_berg_eff_factor ! berg efficiency factor
141 :
142 : ! Switches for specification rather than prediction of droplet and crystal number
143 : ! note: number will be adjusted as needed to keep mean size within bounds,
144 : ! even when specified droplet or ice number is used
145 : !
146 : ! If constant cloud ice number is set (nicons = .true.),
147 : ! then all microphysical processes except mass transfer due to ice nucleation
148 : ! (mnuccd) are based on the fixed cloud ice number. Calculation of
149 : ! mnuccd follows from the prognosed ice crystal number ni.
150 : logical :: nccons ! nccons=.true. to specify constant cloud droplet number
151 : logical :: nicons ! nicons=.true. to specify constant cloud ice number
152 :
153 : ! parameters for specified ice and droplet number concentration
154 : ! note: these are local in-cloud values, not grid-mean
155 : real(r8) :: ncnst ! droplet num concentration when nccons=.true. (m-3)
156 : real(r8) :: ninst ! ice num concentration when nicons=.true. (m-3)
157 :
158 : !===============================================================================
159 : contains
160 : !===============================================================================
161 :
162 0 : subroutine micro_mg_init( &
163 : kind, gravit, rair, rh2o, cpair, &
164 : rhoh2o, tmelt_in, latvap, latice, &
165 : rhmini_in, micro_mg_dcs, use_hetfrz_classnuc_in, &
166 0 : micro_mg_precip_frac_method_in, micro_mg_berg_eff_factor_in, &
167 0 : nccons_in, nicons_in, ncnst_in, ninst_in, errstring)
168 :
169 : !-----------------------------------------------------------------------
170 : !
171 : ! Purpose:
172 : ! initialize constants for the morrison microphysics
173 : !
174 : ! Author: Andrew Gettelman Dec 2005
175 : !
176 : !-----------------------------------------------------------------------
177 :
178 : integer, intent(in) :: kind ! Kind used for reals
179 : real(r8), intent(in) :: gravit
180 : real(r8), intent(in) :: rair
181 : real(r8), intent(in) :: rh2o
182 : real(r8), intent(in) :: cpair
183 : real(r8), intent(in) :: rhoh2o
184 : real(r8), intent(in) :: tmelt_in ! Freezing point of water (K)
185 : real(r8), intent(in) :: latvap
186 : real(r8), intent(in) :: latice
187 : real(r8), intent(in) :: rhmini_in ! Minimum rh for ice cloud fraction > 0.
188 : real(r8), intent(in) :: micro_mg_dcs
189 : logical, intent(in) :: use_hetfrz_classnuc_in
190 : character(len=16),intent(in) :: micro_mg_precip_frac_method_in ! type of precipitation fraction method
191 : real(r8), intent(in) :: micro_mg_berg_eff_factor_in ! berg efficiency factor
192 : logical, intent(in) :: nccons_in
193 : logical, intent(in) :: nicons_in
194 : real(r8), intent(in) :: ncnst_in
195 : real(r8), intent(in) :: ninst_in
196 :
197 : character(128), intent(out) :: errstring ! Output status (non-blank for error return)
198 :
199 : integer k
200 :
201 : integer l,m, iaer
202 : real(r8) surften ! surface tension of water w/respect to air (N/m)
203 : real(r8) arg
204 : !-----------------------------------------------------------------------
205 :
206 0 : errstring = ' '
207 :
208 0 : if( kind .ne. r8 ) then
209 0 : errstring = 'micro_mg_init: KIND of reals does not match'
210 0 : return
211 : end if
212 :
213 : !declarations for morrison codes (transforms variable names)
214 :
215 0 : g= gravit !gravity
216 0 : r= rair !Dry air Gas constant: note units(phys_constants are in J/K/kmol)
217 0 : rv= rh2o !water vapor gas contstant
218 0 : cpp = cpair !specific heat of dry air
219 0 : rhow = rhoh2o !density of liquid water
220 0 : tmelt = tmelt_in
221 0 : rhmini = rhmini_in
222 0 : micro_mg_precip_frac_method = micro_mg_precip_frac_method_in
223 0 : micro_mg_berg_eff_factor = micro_mg_berg_eff_factor_in
224 :
225 0 : nccons = nccons_in
226 0 : nicons = nicons_in
227 0 : ncnst = ncnst_in
228 0 : ninst = ninst_in
229 :
230 : ! latent heats
231 :
232 0 : xxlv = latvap ! latent heat vaporization
233 0 : xlf = latice ! latent heat freezing
234 0 : xxls = xxlv + xlf ! latent heat of sublimation
235 :
236 : ! flags
237 0 : use_hetfrz_classnuc = use_hetfrz_classnuc_in
238 :
239 : ! parameters for snow/rain fraction for convective clouds
240 :
241 0 : tmax_fsnow = tmelt
242 0 : tmin_fsnow = tmelt-5._r8
243 :
244 : ! parameters below from Reisner et al. (1998)
245 : ! density parameters (kg/m3)
246 :
247 0 : rhosn = 250._r8 ! bulk density snow (++ ceh)
248 0 : rhoi = 500._r8 ! bulk density ice
249 0 : rhow = 1000._r8 ! bulk density liquid
250 :
251 :
252 : ! fall speed parameters, V = aD^b
253 : ! V is in m/s
254 :
255 : ! droplets
256 0 : ac = 3.e7_r8
257 0 : bc = 2._r8
258 :
259 : ! snow
260 0 : as = 11.72_r8
261 0 : bs = 0.41_r8
262 :
263 : ! cloud ice
264 0 : ai = 700._r8
265 0 : bi = 1._r8
266 :
267 : ! rain
268 0 : ar = 841.99667_r8
269 0 : br = 0.8_r8
270 :
271 : ! particle mass-diameter relationship
272 : ! currently we assume spherical particles for cloud ice/snow
273 : ! m = cD^d
274 :
275 0 : pi= 3.1415927_r8
276 :
277 : ! cloud ice mass-diameter relationship
278 :
279 0 : ci = rhoi*pi/6._r8
280 0 : di = 3._r8
281 :
282 : ! snow mass-diameter relationship
283 :
284 0 : cs = rhosn*pi/6._r8
285 0 : ds = 3._r8
286 :
287 : ! drop mass-diameter relationship
288 :
289 0 : cr = rhow*pi/6._r8
290 0 : dr = 3._r8
291 :
292 : ! ventilation parameters for snow
293 : ! hall and prupacher
294 :
295 0 : f1s = 0.86_r8
296 0 : f2s = 0.28_r8
297 :
298 : ! collection efficiency, aggregation of cloud ice and snow
299 :
300 0 : Eii = 0.1_r8
301 :
302 : ! collection efficiency, accretion of cloud water by rain
303 :
304 0 : Ecr = 1.0_r8
305 :
306 : ! ventilation constants for rain
307 :
308 0 : f1r = 0.78_r8
309 0 : f2r = 0.32_r8
310 :
311 : ! autoconversion size threshold for cloud ice to snow (m)
312 :
313 0 : Dcs = micro_mg_dcs
314 :
315 : ! smallest mixing ratio considered in microphysics
316 :
317 0 : qsmall = 1.e-18_r8
318 :
319 : ! immersion freezing parameters, bigg 1953
320 :
321 0 : bimm = 100._r8
322 0 : aimm = 0.66_r8
323 :
324 : ! typical air density at 850 mb
325 :
326 0 : rhosu = 85000._r8/(rair * tmelt)
327 :
328 : ! mass of new crystal due to aerosol freezing and growth (kg)
329 :
330 0 : mi0 = 4._r8/3._r8*pi*rhoi*(10.e-6_r8)*(10.e-6_r8)*(10.e-6_r8)
331 :
332 : ! radius of contact nuclei aerosol (m)
333 :
334 0 : rin = 0.1e-6_r8
335 :
336 : ! freezing temperature
337 0 : tt0=273.15_r8
338 :
339 0 : pi=4._r8*atan(1.0_r8)
340 :
341 : !Range of cloudsat reflectivities (dBz) for analytic simulator
342 0 : csmin= -30._r8
343 0 : csmax= 26._r8
344 0 : mindbz = -99._r8
345 : ! minrefl = 10._r8**(mindbz/10._r8)
346 0 : minrefl = 1.26e-10_r8
347 :
348 : ! Define constants to help speed up code (limit calls to gamma function)
349 :
350 0 : cons1=gamma(1._r8+di)
351 0 : cons4=gamma(1._r8+br)
352 0 : cons5=gamma(4._r8+br)
353 0 : cons6=gamma(1._r8+ds)
354 0 : cons7=gamma(1._r8+bs)
355 0 : cons8=gamma(4._r8+bs)
356 0 : cons11=gamma(3._r8+bs)
357 0 : cons13=gamma(5._r8/2._r8+br/2._r8)
358 0 : cons14=gamma(5._r8/2._r8+bs/2._r8)
359 0 : cons16=gamma(1._r8+bi)
360 0 : cons17=gamma(4._r8+bi)
361 0 : cons22=(4._r8/3._r8*pi*rhow*(25.e-6_r8)**3)
362 0 : cons23=dcs**3
363 0 : cons24=dcs**2
364 0 : cons25=dcs**bs
365 0 : cons27=xxlv**2
366 0 : cons28=xxls**2
367 :
368 0 : lammaxi = 1._r8/10.e-6_r8
369 0 : lammini = 1._r8/(2._r8*dcs)
370 0 : lammaxr = 1._r8/20.e-6_r8
371 0 : lamminr = 1._r8/500.e-6_r8
372 0 : lammaxs = 1._r8/10.e-6_r8
373 0 : lammins = 1._r8/2000.e-6_r8
374 :
375 : end subroutine micro_mg_init
376 :
377 : !===============================================================================
378 : !microphysics routine for each timestep goes here...
379 :
380 0 : subroutine micro_mg_tend ( &
381 : microp_uniform, pcols, pver, ncol, top_lev, deltatin,&
382 0 : tn, qn, qc, qi, nc, &
383 0 : ni, p, pdel, cldn, liqcldf, &
384 0 : relvar, accre_enhan, &
385 0 : icecldf, rate1ord_cw2pr_st, naai, npccnin, &
386 0 : rndst, nacon, tlat, qvlat, qctend, &
387 0 : qitend, nctend, nitend, effc, effc_fn, &
388 0 : effi, prect, preci, nevapr, evapsnow, am_evp_st, &
389 0 : prain, prodsnow, cmeout, deffi, pgamrad, &
390 0 : lamcrad, qsout, dsout, rflx, sflx, &
391 0 : qrout, reff_rain, reff_snow, qcsevap, qisevap, &
392 0 : qvres, cmeiout, vtrmc, vtrmi, qcsedten, &
393 0 : qisedten, prao, prco, mnuccco, mnuccto, &
394 0 : msacwio, psacwso, bergso, bergo, melto, &
395 0 : homoo, qcreso, prcio, praio, qireso, &
396 0 : mnuccro, pracso, meltsdt, frzrdt, mnuccdo, &
397 0 : nrout, nsout, refl, arefl, areflz, &
398 0 : frefl, csrfl, acsrfl, fcsrfl, rercld, &
399 0 : ncai, ncal, qrout2, qsout2, nrout2, &
400 0 : nsout2, drout2, dsout2, freqs, freqr, &
401 0 : nfice, prer_evap, do_cldice, errstring, &
402 0 : tnd_qsnow, tnd_nsnow, re_ice, &
403 0 : frzimm, frzcnt, frzdep)
404 :
405 : ! input arguments
406 : logical, intent(in) :: microp_uniform ! True = configure uniform for sub-columns False = use w/o sub-columns (standard)
407 : integer, intent(in) :: pcols ! size of column (first) index
408 : integer, intent(in) :: pver ! number of layers in columns
409 : integer, intent(in) :: ncol ! number of columns
410 : integer, intent(in) :: top_lev ! top level microphys is applied
411 : real(r8), intent(in) :: deltatin ! time step (s)
412 : real(r8), intent(in) :: tn(pcols,pver) ! input temperature (K)
413 : real(r8), intent(in) :: qn(pcols,pver) ! input h20 vapor mixing ratio (kg/kg)
414 : real(r8), intent(in) :: relvar(pcols,pver) ! relative variance of cloud water (-)
415 : real(r8), intent(in) :: accre_enhan(pcols,pver) ! optional accretion enhancement factor (-)
416 :
417 : ! note: all input cloud variables are grid-averaged
418 : real(r8), intent(inout) :: qc(pcols,pver) ! cloud water mixing ratio (kg/kg)
419 : real(r8), intent(inout) :: qi(pcols,pver) ! cloud ice mixing ratio (kg/kg)
420 : real(r8), intent(inout) :: nc(pcols,pver) ! cloud water number conc (1/kg)
421 : real(r8), intent(inout) :: ni(pcols,pver) ! cloud ice number conc (1/kg)
422 : real(r8), intent(in) :: p(pcols,pver) ! air pressure (pa)
423 : real(r8), intent(in) :: pdel(pcols,pver) ! pressure difference across level (pa)
424 : real(r8), intent(in) :: cldn(pcols,pver) ! cloud fraction
425 : real(r8), intent(in) :: icecldf(pcols,pver) ! ice cloud fraction
426 : real(r8), intent(in) :: liqcldf(pcols,pver) ! liquid cloud fraction
427 :
428 : real(r8), intent(out) :: rate1ord_cw2pr_st(pcols,pver) ! 1st order rate for direct cw to precip conversion
429 : ! used for scavenging
430 : ! Inputs for aerosol activation
431 : real(r8), intent(in) :: naai(pcols,pver) ! ice nulceation number (from microp_aero_ts)
432 : real(r8), intent(in) :: npccnin(pcols,pver) ! ccn activated number tendency (from microp_aero_ts)
433 : real(r8), intent(in) :: rndst(pcols,pver,4) ! radius of 4 dust bins for contact freezing (from microp_aero_ts)
434 : real(r8), intent(in) :: nacon(pcols,pver,4) ! number in 4 dust bins for contact freezing (from microp_aero_ts)
435 :
436 : ! Used with CARMA cirrus microphysics
437 : ! (or similar external microphysics model)
438 : logical, intent(in) :: do_cldice ! Prognosing cldice
439 :
440 : ! output arguments
441 :
442 : real(r8), intent(out) :: tlat(pcols,pver) ! latent heating rate (W/kg)
443 : real(r8), intent(out) :: qvlat(pcols,pver) ! microphysical tendency qv (1/s)
444 : real(r8), intent(out) :: qctend(pcols,pver) ! microphysical tendency qc (1/s)
445 : real(r8), intent(out) :: qitend(pcols,pver) ! microphysical tendency qi (1/s)
446 : real(r8), intent(out) :: nctend(pcols,pver) ! microphysical tendency nc (1/(kg*s))
447 : real(r8), intent(out) :: nitend(pcols,pver) ! microphysical tendency ni (1/(kg*s))
448 : real(r8), intent(out) :: effc(pcols,pver) ! droplet effective radius (micron)
449 : real(r8), intent(out) :: effc_fn(pcols,pver) ! droplet effective radius, assuming nc = 1.e8 kg-1
450 : real(r8), intent(out) :: effi(pcols,pver) ! cloud ice effective radius (micron)
451 : real(r8), intent(out) :: prect(pcols) ! surface precip rate (m/s)
452 : real(r8), intent(out) :: preci(pcols) ! cloud ice/snow precip rate (m/s)
453 : real(r8), intent(out) :: nevapr(pcols,pver) ! evaporation rate of rain + snow
454 : real(r8), intent(out) :: evapsnow(pcols,pver)! sublimation rate of snow
455 : real(r8), intent(out) :: am_evp_st(pcols,pver)! stratiform evaporation area
456 : real(r8), intent(out) :: prain(pcols,pver) ! production of rain + snow
457 : real(r8), intent(out) :: prodsnow(pcols,pver)! production of snow
458 : real(r8), intent(out) :: cmeout(pcols,pver) ! evap/sub of cloud
459 : real(r8), intent(out) :: deffi(pcols,pver) ! ice effective diameter for optics (radiation)
460 : real(r8), intent(out) :: pgamrad(pcols,pver) ! ice gamma parameter for optics (radiation)
461 : real(r8), intent(out) :: lamcrad(pcols,pver) ! slope of droplet distribution for optics (radiation)
462 : real(r8), intent(out) :: qsout(pcols,pver) ! snow mixing ratio (kg/kg)
463 : real(r8), intent(out) :: dsout(pcols,pver) ! snow diameter (m)
464 : real(r8), intent(out) :: rflx(pcols,pver+1) ! grid-box average rain flux (kg m^-2 s^-1)
465 : real(r8), intent(out) :: sflx(pcols,pver+1) ! grid-box average snow flux (kg m^-2 s^-1)
466 : real(r8), intent(out) :: qrout(pcols,pver) ! grid-box average rain mixing ratio (kg/kg)
467 : real(r8), intent(inout) :: reff_rain(pcols,pver) ! rain effective radius (micron)
468 : real(r8), intent(inout) :: reff_snow(pcols,pver) ! snow effective radius (micron)
469 : real(r8), intent(out) :: qcsevap(pcols,pver) ! cloud water evaporation due to sedimentation
470 : real(r8), intent(out) :: qisevap(pcols,pver) ! cloud ice sublimation due to sublimation
471 : real(r8), intent(out) :: qvres(pcols,pver) ! residual condensation term to ensure RH < 100%
472 : real(r8), intent(out) :: cmeiout(pcols,pver) ! grid-mean cloud ice sub/dep
473 : real(r8), intent(out) :: vtrmc(pcols,pver) ! mass-weighted cloud water fallspeed
474 : real(r8), intent(out) :: vtrmi(pcols,pver) ! mass-weighted cloud ice fallspeed
475 : real(r8), intent(out) :: qcsedten(pcols,pver) ! qc sedimentation tendency
476 : real(r8), intent(out) :: qisedten(pcols,pver) ! qi sedimentation tendency
477 : ! microphysical process rates for output (mixing ratio tendencies)
478 : real(r8), intent(out) :: prao(pcols,pver) ! accretion of cloud by rain
479 : real(r8), intent(out) :: prco(pcols,pver) ! autoconversion of cloud to rain
480 : real(r8), intent(out) :: mnuccco(pcols,pver) ! mixing rat tend due to immersion freezing
481 : real(r8), intent(out) :: mnuccto(pcols,pver) ! mixing ratio tend due to contact freezing
482 : real(r8), intent(out) :: msacwio(pcols,pver) ! mixing ratio tend due to H-M splintering
483 : real(r8), intent(out) :: psacwso(pcols,pver) ! collection of cloud water by snow
484 : real(r8), intent(out) :: bergso(pcols,pver) ! bergeron process on snow
485 : real(r8), intent(out) :: bergo(pcols,pver) ! bergeron process on cloud ice
486 : real(r8), intent(out) :: melto(pcols,pver) ! melting of cloud ice
487 : real(r8), intent(out) :: homoo(pcols,pver) ! homogeneos freezign cloud water
488 : real(r8), intent(out) :: qcreso(pcols,pver) ! residual cloud condensation due to removal of excess supersat
489 : real(r8), intent(out) :: prcio(pcols,pver) ! autoconversion of cloud ice to snow
490 : real(r8), intent(out) :: praio(pcols,pver) ! accretion of cloud ice by snow
491 : real(r8), intent(out) :: qireso(pcols,pver) ! residual ice deposition due to removal of excess supersat
492 : real(r8), intent(out) :: mnuccro(pcols,pver) ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s)
493 : real(r8), intent(out) :: pracso (pcols,pver) ! mixing ratio tendency due to accretion of rain by snow (1/s)
494 : real(r8), intent(out) :: meltsdt(pcols,pver) ! latent heating rate due to melting of snow (W/kg)
495 : real(r8), intent(out) :: frzrdt (pcols,pver) ! latent heating rate due to homogeneous freezing of rain (W/kg)
496 : real(r8), intent(out) :: mnuccdo(pcols,pver) ! mass tendency from ice nucleation
497 : real(r8), intent(out) :: nrout(pcols,pver) ! rain number concentration (1/m3)
498 : real(r8), intent(out) :: nsout(pcols,pver) ! snow number concentration (1/m3)
499 : real(r8), intent(out) :: refl(pcols,pver) ! analytic radar reflectivity
500 : real(r8), intent(out) :: arefl(pcols,pver) !average reflectivity will zero points outside valid range
501 : real(r8), intent(out) :: areflz(pcols,pver) !average reflectivity in z.
502 : real(r8), intent(out) :: frefl(pcols,pver)
503 : real(r8), intent(out) :: csrfl(pcols,pver) !cloudsat reflectivity
504 : real(r8), intent(out) :: acsrfl(pcols,pver) !cloudsat average
505 : real(r8), intent(out) :: fcsrfl(pcols,pver)
506 : real(r8), intent(out) :: rercld(pcols,pver) ! effective radius calculation for rain + cloud
507 : real(r8), intent(out) :: ncai(pcols,pver) ! output number conc of ice nuclei available (1/m3)
508 : real(r8), intent(out) :: ncal(pcols,pver) ! output number conc of CCN (1/m3)
509 : real(r8), intent(out) :: qrout2(pcols,pver)
510 : real(r8), intent(out) :: qsout2(pcols,pver)
511 : real(r8), intent(out) :: nrout2(pcols,pver)
512 : real(r8), intent(out) :: nsout2(pcols,pver)
513 : real(r8), intent(out) :: drout2(pcols,pver) ! mean rain particle diameter (m)
514 : real(r8), intent(out) :: dsout2(pcols,pver) ! mean snow particle diameter (m)
515 : real(r8), intent(out) :: freqs(pcols,pver)
516 : real(r8), intent(out) :: freqr(pcols,pver)
517 : real(r8), intent(out) :: nfice(pcols,pver)
518 : real(r8), intent(out) :: prer_evap(pcols,pver)
519 :
520 0 : real(r8) :: nevapr2(pcols,pver)
521 :
522 : character(128), intent(out) :: errstring ! Output status (non-blank for error return)
523 :
524 : ! Tendencies calculated by external schemes that can replace MG's native
525 : ! process tendencies.
526 :
527 : ! Used with CARMA cirrus microphysics
528 : ! (or similar external microphysics model)
529 : real(r8), intent(in) :: tnd_qsnow(:,:) ! snow mass tendency (kg/kg/s)
530 : real(r8), intent(in) :: tnd_nsnow(:,:) ! snow number tendency (#/kg/s)
531 : real(r8), intent(in) :: re_ice(:,:) ! ice effective radius (m)
532 :
533 : ! From external ice nucleation.
534 : real(r8), intent(in) :: frzimm(:,:) ! Number tendency due to immersion freezing (1/cm3)
535 : real(r8), intent(in) :: frzcnt(:,:) ! Number tendency due to contact freezing (1/cm3)
536 : real(r8), intent(in) :: frzdep(:,:) ! Number tendency due to deposition nucleation (1/cm3)
537 :
538 : ! local workspace
539 : ! all units mks unless otherwise stated
540 :
541 : ! Additional constants to help speed up code
542 : real(r8) :: cons2
543 : real(r8) :: cons3
544 : real(r8) :: cons9
545 : real(r8) :: cons10
546 : real(r8) :: cons12
547 : real(r8) :: cons15
548 : real(r8) :: cons18
549 : real(r8) :: cons19
550 : real(r8) :: cons20
551 :
552 : ! temporary variables for sub-stepping
553 0 : real(r8) :: t1(pcols,pver)
554 0 : real(r8) :: q1(pcols,pver)
555 0 : real(r8) :: qc1(pcols,pver)
556 0 : real(r8) :: qi1(pcols,pver)
557 0 : real(r8) :: nc1(pcols,pver)
558 0 : real(r8) :: ni1(pcols,pver)
559 0 : real(r8) :: tlat1(pcols,pver)
560 0 : real(r8) :: qvlat1(pcols,pver)
561 0 : real(r8) :: qctend1(pcols,pver)
562 0 : real(r8) :: qitend1(pcols,pver)
563 0 : real(r8) :: nctend1(pcols,pver)
564 0 : real(r8) :: nitend1(pcols,pver)
565 0 : real(r8) :: prect1(pcols)
566 0 : real(r8) :: preci1(pcols)
567 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
568 :
569 : real(r8) :: deltat ! sub-time step (s)
570 : real(r8) :: omsm ! number near unity for round-off issues
571 : real(r8) :: dto2 ! dt/2 (s)
572 : real(r8) :: mincld ! minimum allowed cloud fraction
573 0 : real(r8) :: q(pcols,pver) ! water vapor mixing ratio (kg/kg)
574 0 : real(r8) :: t(pcols,pver) ! temperature (K)
575 0 : real(r8) :: rho(pcols,pver) ! air density (kg m-3)
576 0 : real(r8) :: dv(pcols,pver) ! diffusivity of water vapor in air
577 0 : real(r8) :: mu(pcols,pver) ! viscocity of air
578 0 : real(r8) :: sc(pcols,pver) ! schmidt number
579 0 : real(r8) :: kap(pcols,pver) ! thermal conductivity of air
580 0 : real(r8) :: rhof(pcols,pver) ! air density correction factor for fallspeed
581 0 : real(r8) :: cldmax(pcols,pver) ! precip fraction assuming maximum overlap
582 0 : real(r8) :: cldm(pcols,pver) ! cloud fraction
583 0 : real(r8) :: icldm(pcols,pver) ! ice cloud fraction
584 0 : real(r8) :: lcldm(pcols,pver) ! liq cloud fraction
585 : real(r8) :: icwc(pcols) ! in cloud water content (liquid+ice)
586 : real(r8) :: calpha(pcols) ! parameter for cond/evap (Zhang et al. 2003)
587 : real(r8) :: cbeta(pcols) ! parameter for cond/evap (Zhang et al. 2003)
588 : real(r8) :: cbetah(pcols) ! parameter for cond/evap (Zhang et al. 2003)
589 : real(r8) :: cgamma(pcols) ! parameter for cond/evap (Zhang et al. 2003)
590 : real(r8) :: cgamah(pcols) ! parameter for cond/evap (Zhang et al. 2003)
591 : real(r8) :: rcgama(pcols) ! parameter for cond/evap (Zhang et al. 2003)
592 : real(r8) :: cmec1(pcols) ! parameter for cond/evap (Zhang et al. 2003)
593 : real(r8) :: cmec2(pcols) ! parameter for cond/evap (Zhang et al. 2003)
594 : real(r8) :: cmec3(pcols) ! parameter for cond/evap (Zhang et al. 2003)
595 : real(r8) :: cmec4(pcols) ! parameter for cond/evap (Zhang et al. 2003)
596 : real(r8) :: qtmp ! dummy qv
597 : real(r8) :: dum ! temporary dummy variable
598 :
599 0 : real(r8) :: cme(pcols,pver) ! total (liquid+ice) cond/evap rate of cloud
600 :
601 0 : real(r8) :: cmei(pcols,pver) ! dep/sublimation rate of cloud ice
602 0 : real(r8) :: cwml(pcols,pver) ! cloud water mixing ratio
603 0 : real(r8) :: cwmi(pcols,pver) ! cloud ice mixing ratio
604 0 : real(r8) :: nnuccd(pver) ! ice nucleation rate from deposition/cond.-freezing
605 0 : real(r8) :: mnuccd(pver) ! mass tendency from ice nucleation
606 : real(r8) :: qcld ! total cloud water
607 : real(r8) :: lcldn(pcols,pver) ! fractional coverage of new liquid cloud
608 : real(r8) :: lcldo(pcols,pver) ! fractional coverage of old liquid cloud
609 : real(r8) :: nctend_mixnuc(pcols,pver)
610 : real(r8) :: arg ! argument of erfc
611 :
612 : ! for calculation of rate1ord_cw2pr_st
613 0 : real(r8) :: qcsinksum_rate1ord(pver) ! sum over iterations of cw to precip sink
614 0 : real(r8) :: qcsum_rate1ord(pver) ! sum over iterations of cloud water
615 :
616 : real(r8) :: alpha
617 :
618 : real(r8) :: dum1,dum2 !general dummy variables
619 :
620 0 : real(r8) :: npccn(pver) ! droplet activation rate
621 0 : real(r8) :: qcic(pcols,pver) ! in-cloud cloud liquid mixing ratio
622 0 : real(r8) :: qiic(pcols,pver) ! in-cloud cloud ice mixing ratio
623 0 : real(r8) :: qniic(pcols,pver) ! in-precip snow mixing ratio
624 0 : real(r8) :: qric(pcols,pver) ! in-precip rain mixing ratio
625 0 : real(r8) :: ncic(pcols,pver) ! in-cloud droplet number conc
626 0 : real(r8) :: niic(pcols,pver) ! in-cloud cloud ice number conc
627 0 : real(r8) :: nsic(pcols,pver) ! in-precip snow number conc
628 0 : real(r8) :: nric(pcols,pver) ! in-precip rain number conc
629 0 : real(r8) :: lami(pver) ! slope of cloud ice size distr
630 0 : real(r8) :: n0i(pver) ! intercept of cloud ice size distr
631 0 : real(r8) :: lamc(pver) ! slope of cloud liquid size distr
632 : real(r8) :: n0c(pver) ! intercept of cloud liquid size distr
633 0 : real(r8) :: lams(pver) ! slope of snow size distr
634 0 : real(r8) :: n0s(pver) ! intercept of snow size distr
635 0 : real(r8) :: lamr(pver) ! slope of rain size distr
636 0 : real(r8) :: n0r(pver) ! intercept of rain size distr
637 0 : real(r8) :: cdist1(pver) ! size distr parameter to calculate droplet freezing
638 : ! combined size of precip & cloud drops
639 0 : real(r8) :: arcld(pcols,pver) ! averaging control flag
640 : real(r8) :: Actmp !area cross section of drops
641 : real(r8) :: Artmp !area cross section of rain
642 :
643 0 : real(r8) :: pgam(pver) ! spectral width parameter of droplet size distr
644 : real(r8) :: lammax ! maximum allowed slope of size distr
645 : real(r8) :: lammin ! minimum allowed slope of size distr
646 : real(r8) :: nacnt ! number conc of contact ice nuclei
647 0 : real(r8) :: mnuccc(pver) ! mixing ratio tendency due to freezing of cloud water
648 0 : real(r8) :: nnuccc(pver) ! number conc tendency due to freezing of cloud water
649 :
650 0 : real(r8) :: mnucct(pver) ! mixing ratio tendency due to contact freezing of cloud water
651 0 : real(r8) :: nnucct(pver) ! number conc tendency due to contact freezing of cloud water
652 0 : real(r8) :: msacwi(pver) ! mixing ratio tendency due to HM ice multiplication
653 0 : real(r8) :: nsacwi(pver) ! number conc tendency due to HM ice multiplication
654 :
655 0 : real(r8) :: prc(pver) ! qc tendency due to autoconversion of cloud droplets
656 0 : real(r8) :: nprc(pver) ! number conc tendency due to autoconversion of cloud droplets
657 0 : real(r8) :: nprc1(pver) ! qr tendency due to autoconversion of cloud droplets
658 0 : real(r8) :: nsagg(pver) ! ns tendency due to self-aggregation of snow
659 : real(r8) :: dc0 ! mean size droplet size distr
660 : real(r8) :: ds0 ! mean size snow size distr (area weighted)
661 : real(r8) :: eci ! collection efficiency for riming of snow by droplets
662 0 : real(r8) :: psacws(pver) ! mixing rat tendency due to collection of droplets by snow
663 0 : real(r8) :: npsacws(pver) ! number conc tendency due to collection of droplets by snow
664 : real(r8) :: uni ! number-weighted cloud ice fallspeed
665 : real(r8) :: umi ! mass-weighted cloud ice fallspeed
666 0 : real(r8) :: uns(pver) ! number-weighted snow fallspeed
667 0 : real(r8) :: ums(pver) ! mass-weighted snow fallspeed
668 0 : real(r8) :: unr(pver) ! number-weighted rain fallspeed
669 0 : real(r8) :: umr(pver) ! mass-weighted rain fallspeed
670 : real(r8) :: unc ! number-weighted cloud droplet fallspeed
671 : real(r8) :: umc ! mass-weighted cloud droplet fallspeed
672 0 : real(r8) :: pracs(pver) ! mixing rat tendency due to collection of rain by snow
673 0 : real(r8) :: npracs(pver) ! number conc tendency due to collection of rain by snow
674 0 : real(r8) :: mnuccr(pver) ! mixing rat tendency due to freezing of rain
675 0 : real(r8) :: nnuccr(pver) ! number conc tendency due to freezing of rain
676 0 : real(r8) :: pra(pver) ! mixing rat tendnency due to accretion of droplets by rain
677 0 : real(r8) :: npra(pver) ! nc tendnency due to accretion of droplets by rain
678 0 : real(r8) :: nragg(pver) ! nr tendency due to self-collection of rain
679 0 : real(r8) :: prci(pver) ! mixing rat tendency due to autoconversion of cloud ice to snow
680 0 : real(r8) :: nprci(pver) ! number conc tendency due to autoconversion of cloud ice to snow
681 0 : real(r8) :: prai(pver) ! mixing rat tendency due to accretion of cloud ice by snow
682 0 : real(r8) :: nprai(pver) ! number conc tendency due to accretion of cloud ice by snow
683 : real(r8) :: qvs ! liquid saturation vapor mixing ratio
684 : real(r8) :: qvi ! ice saturation vapor mixing ratio
685 : real(r8) :: dqsdt ! change of sat vapor mixing ratio with temperature
686 : real(r8) :: dqsidt ! change of ice sat vapor mixing ratio with temperature
687 : real(r8) :: ab ! correction factor for rain evap to account for latent heat
688 : real(r8) :: qclr ! water vapor mixing ratio in clear air
689 : real(r8) :: abi ! correction factor for snow sublimation to account for latent heat
690 : real(r8) :: epss ! 1/ sat relaxation timescale for snow
691 : real(r8) :: epsr ! 1/ sat relaxation timescale for rain
692 0 : real(r8) :: pre(pver) ! rain mixing rat tendency due to evaporation
693 0 : real(r8) :: prds(pver) ! snow mixing rat tendency due to sublimation
694 : real(r8) :: qce ! dummy qc for conservation check
695 : real(r8) :: qie ! dummy qi for conservation check
696 : real(r8) :: nce ! dummy nc for conservation check
697 : real(r8) :: nie ! dummy ni for conservation check
698 : real(r8) :: ratio ! parameter for conservation check
699 0 : real(r8) :: dumc(pcols,pver) ! dummy in-cloud qc
700 0 : real(r8) :: dumnc(pcols,pver) ! dummy in-cloud nc
701 0 : real(r8) :: dumi(pcols,pver) ! dummy in-cloud qi
702 0 : real(r8) :: dumni(pcols,pver) ! dummy in-cloud ni
703 : real(r8) :: dums(pcols,pver) ! dummy in-cloud snow mixing rat
704 : real(r8) :: dumns(pcols,pver) ! dummy in-cloud snow number conc
705 : real(r8) :: dumr(pcols,pver) ! dummy in-cloud rain mixing rat
706 : real(r8) :: dumnr(pcols,pver) ! dummy in-cloud rain number conc
707 : ! below are parameters for cloud water and cloud ice sedimentation calculations
708 : real(r8) :: fr(pver)
709 : real(r8) :: fnr(pver)
710 0 : real(r8) :: fc(pver)
711 0 : real(r8) :: fnc(pver)
712 0 : real(r8) :: fi(pver)
713 0 : real(r8) :: fni(pver)
714 : real(r8) :: fs(pver)
715 : real(r8) :: fns(pver)
716 : real(r8) :: faloutr(pver)
717 : real(r8) :: faloutnr(pver)
718 0 : real(r8) :: faloutc(pver)
719 0 : real(r8) :: faloutnc(pver)
720 0 : real(r8) :: falouti(pver)
721 0 : real(r8) :: faloutni(pver)
722 : real(r8) :: falouts(pver)
723 : real(r8) :: faloutns(pver)
724 : real(r8) :: faltndr
725 : real(r8) :: faltndnr
726 : real(r8) :: faltndc
727 : real(r8) :: faltndnc
728 : real(r8) :: faltndi
729 : real(r8) :: faltndni
730 : real(r8) :: faltnds
731 : real(r8) :: faltndns
732 : real(r8) :: faltndqie
733 : real(r8) :: faltndqce
734 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
735 0 : real(r8) :: relhum(pcols,pver) ! relative humidity
736 : real(r8) :: csigma(pcols) ! parameter for cond/evap of cloud water/ice
737 : real(r8) :: rgvm ! max fallspeed for all species
738 0 : real(r8) :: arn(pcols,pver) ! air density corrected rain fallspeed parameter
739 0 : real(r8) :: asn(pcols,pver) ! air density corrected snow fallspeed parameter
740 0 : real(r8) :: acn(pcols,pver) ! air density corrected cloud droplet fallspeed parameter
741 0 : real(r8) :: ain(pcols,pver) ! air density corrected cloud ice fallspeed parameter
742 0 : real(r8) :: nsubi(pver) ! evaporation of cloud ice number
743 0 : real(r8) :: nsubc(pver) ! evaporation of droplet number
744 0 : real(r8) :: nsubs(pver) ! evaporation of snow number
745 0 : real(r8) :: nsubr(pver) ! evaporation of rain number
746 : real(r8) :: mtime ! factor to account for droplet activation timescale
747 0 : real(r8) :: dz(pcols,pver) ! height difference across model vertical level
748 :
749 :
750 : !! add precip flux variables for sub-stepping
751 0 : real(r8) :: rflx1(pcols,pver+1)
752 0 : real(r8) :: sflx1(pcols,pver+1)
753 :
754 : ! returns from function/subroutine calls
755 : real(r8) :: tsp(pcols,pver) ! saturation temp (K)
756 : real(r8) :: qsp(pcols,pver) ! saturation mixing ratio (kg/kg)
757 : real(r8) :: qsphy(pcols,pver) ! saturation mixing ratio (kg/kg): hybrid rh
758 0 : real(r8) :: qs(pcols) ! liquid-ice weighted sat mixing rat (kg/kg)
759 0 : real(r8) :: es(pcols) ! liquid-ice weighted sat vapor press (pa)
760 0 : real(r8) :: esl(pcols,pver) ! liquid sat vapor pressure (pa)
761 0 : real(r8) :: esi(pcols,pver) ! ice sat vapor pressure (pa)
762 :
763 : ! sum of source/sink terms for diagnostic precip
764 :
765 0 : real(r8) :: qnitend(pcols,pver) ! snow mixing ratio source/sink term
766 0 : real(r8) :: nstend(pcols,pver) ! snow number concentration source/sink term
767 0 : real(r8) :: qrtend(pcols,pver) ! rain mixing ratio source/sink term
768 0 : real(r8) :: nrtend(pcols,pver) ! rain number concentration source/sink term
769 : real(r8) :: qrtot ! vertically-integrated rain mixing rat source/sink term
770 : real(r8) :: nrtot ! vertically-integrated rain number conc source/sink term
771 : real(r8) :: qstot ! vertically-integrated snow mixing rat source/sink term
772 : real(r8) :: nstot ! vertically-integrated snow number conc source/sink term
773 :
774 : ! new terms for Bergeron process
775 :
776 : real(r8) :: dumnnuc ! provisional ice nucleation rate (for calculating bergeron)
777 : real(r8) :: ninew ! provisional cloud ice number conc (for calculating bergeron)
778 : real(r8) :: qinew ! provisional cloud ice mixing ratio (for calculating bergeron)
779 : real(r8) :: qvl ! liquid sat mixing ratio
780 : real(r8) :: epsi ! 1/ sat relaxation timecale for cloud ice
781 : real(r8) :: prd ! provisional deposition rate of cloud ice at water sat
782 0 : real(r8) :: berg(pcols,pver) ! mixing rat tendency due to bergeron process for cloud ice
783 0 : real(r8) :: bergs(pver) ! mixing rat tendency due to bergeron process for snow
784 :
785 : !bergeron terms
786 : real(r8) :: bergtsf !bergeron timescale to remove all liquid
787 : real(r8) :: rhin !modified RH for vapor deposition
788 :
789 : ! diagnostic rain/snow for output to history
790 : ! values are in-precip (local) !!!!
791 :
792 0 : real(r8) :: drout(pcols,pver) ! rain diameter (m)
793 :
794 : !averageed rain/snow for history
795 : real(r8) :: dumfice
796 :
797 : !ice nucleation, droplet activation
798 0 : real(r8) :: dum2i(pcols,pver) ! number conc of ice nuclei available (1/kg)
799 0 : real(r8) :: dum2l(pcols,pver) ! number conc of CCN (1/kg)
800 : real(r8) :: ncmax
801 : real(r8) :: nimax
802 :
803 : real(r8) :: qcvar ! 1/relative variance of sub-grid qc
804 :
805 : ! loop array variables
806 : integer i,k,nstep,n, l
807 : integer ii,kk, m
808 :
809 : ! loop variables for sub-step solution
810 0 : integer iter,it,ltrue(pcols)
811 :
812 : ! used in contact freezing via dust particles
813 : real(r8) tcnt, viscosity, mfp
814 : real(r8) slip1, slip2, slip3, slip4
815 : ! real(r8) dfaer1, dfaer2, dfaer3, dfaer4
816 : ! real(r8) nacon1,nacon2,nacon3,nacon4
817 : real(r8) ndfaer1, ndfaer2, ndfaer3, ndfaer4
818 : real(r8) nslip1, nslip2, nslip3, nslip4
819 :
820 : ! used in ice effective radius
821 : real(r8) bbi, cci, ak, iciwc, rvi
822 :
823 : ! used in Bergeron processe and water vapor deposition
824 : real(r8) Tk, deles, Aprpr, Bprpr, Cice, qi0, Crate, qidep
825 :
826 : ! mean cloud fraction over the time step
827 0 : real(r8) cldmw(pcols,pver)
828 :
829 : ! used in secondary ice production
830 : real(r8) ni_secp
831 :
832 : ! variabels to check for RH after rain evap
833 :
834 : real(r8) :: esn
835 : real(r8) :: qsn
836 : real(r8) :: ttmp
837 :
838 :
839 :
840 0 : real(r8) :: rainrt(pcols,pver) ! rain rate for reflectivity calculation
841 0 : real(r8) :: rainrt1(pcols,pver)
842 : real(r8) :: tmp
843 :
844 : real(r8) dmc,ssmc,dstrn ! variables for modal scheme.
845 :
846 : real(r8), parameter :: cdnl = 0.e6_r8 ! cloud droplet number limiter
847 :
848 : ! heterogeneous freezing
849 0 : real(r8) :: mnudep(pver) ! mixing ratio tendency due to deposition of water vapor
850 0 : real(r8) :: nnudep(pver) ! number conc tendency due to deposition of water vapor
851 : real(r8) :: con1 ! work cnstant
852 : real(r8) :: r3lx ! Mean volume radius (m)
853 : real(r8) :: mi0l
854 : real(r8) :: frztmp
855 :
856 : logical :: do_clubb_sgs
857 :
858 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
859 :
860 : ! Return error message
861 0 : errstring = ' '
862 :
863 0 : call phys_getopts(do_clubb_sgs_out = do_clubb_sgs)
864 :
865 : ! initialize output fields for number conc qand ice nucleation
866 0 : ncai(1:ncol,1:pver)=0._r8
867 0 : ncal(1:ncol,1:pver)=0._r8
868 :
869 : !Initialize rain size
870 0 : rercld(1:ncol,1:pver)=0._r8
871 0 : arcld(1:ncol,1:pver)=0._r8
872 :
873 : !initialize radiation output variables
874 0 : pgamrad(1:ncol,1:pver)=0._r8 ! liquid gamma parameter for optics (radiation)
875 0 : lamcrad(1:ncol,1:pver)=0._r8 ! slope of droplet distribution for optics (radiation)
876 0 : deffi (1:ncol,1:pver)=0._r8 ! slope of droplet distribution for optics (radiation)
877 : !initialize radiation output variables
878 : !initialize water vapor tendency term output
879 0 : qcsevap(1:ncol,1:pver)=0._r8
880 0 : qisevap(1:ncol,1:pver)=0._r8
881 0 : qvres (1:ncol,1:pver)=0._r8
882 0 : cmeiout (1:ncol,1:pver)=0._r8
883 0 : vtrmc (1:ncol,1:pver)=0._r8
884 0 : vtrmi (1:ncol,1:pver)=0._r8
885 0 : qcsedten (1:ncol,1:pver)=0._r8
886 0 : qisedten (1:ncol,1:pver)=0._r8
887 :
888 0 : prao(1:ncol,1:pver)=0._r8
889 0 : prco(1:ncol,1:pver)=0._r8
890 0 : mnuccco(1:ncol,1:pver)=0._r8
891 0 : mnuccto(1:ncol,1:pver)=0._r8
892 0 : msacwio(1:ncol,1:pver)=0._r8
893 0 : psacwso(1:ncol,1:pver)=0._r8
894 0 : bergso(1:ncol,1:pver)=0._r8
895 0 : bergo(1:ncol,1:pver)=0._r8
896 0 : melto(1:ncol,1:pver)=0._r8
897 0 : homoo(1:ncol,1:pver)=0._r8
898 0 : qcreso(1:ncol,1:pver)=0._r8
899 0 : prcio(1:ncol,1:pver)=0._r8
900 0 : praio(1:ncol,1:pver)=0._r8
901 0 : qireso(1:ncol,1:pver)=0._r8
902 0 : mnuccro(1:ncol,1:pver)=0._r8
903 0 : pracso (1:ncol,1:pver)=0._r8
904 0 : meltsdt(1:ncol,1:pver)=0._r8
905 0 : frzrdt (1:ncol,1:pver)=0._r8
906 0 : mnuccdo(1:ncol,1:pver)=0._r8
907 :
908 0 : rflx(:,:)=0._r8
909 0 : sflx(:,:)=0._r8
910 0 : effc(:,:)=0._r8
911 0 : effc_fn(:,:)=0._r8
912 0 : effi(:,:)=0._r8
913 :
914 : ! assign variable deltat for sub-stepping...
915 0 : deltat=deltatin
916 :
917 : ! parameters for scheme
918 :
919 0 : omsm=0.99999_r8
920 0 : dto2=0.5_r8*deltat
921 0 : mincld=0.0001_r8
922 :
923 : ! initialize multi-level fields
924 0 : q(1:ncol,1:pver)=qn(1:ncol,1:pver)
925 0 : t(1:ncol,1:pver)=tn(1:ncol,1:pver)
926 :
927 : ! initialize time-varying parameters
928 :
929 0 : do k=1,pver
930 0 : do i=1,ncol
931 0 : rho(i,k)=p(i,k)/(r*t(i,k))
932 0 : dv(i,k) = 8.794E-5_r8*t(i,k)**1.81_r8/p(i,k)
933 0 : mu(i,k) = 1.496E-6_r8*t(i,k)**1.5_r8/(t(i,k)+120._r8)
934 0 : sc(i,k) = mu(i,k)/(rho(i,k)*dv(i,k))
935 0 : kap(i,k) = 1.414e3_r8*1.496e-6_r8*t(i,k)**1.5_r8/(t(i,k)+120._r8)
936 :
937 : ! air density adjustment for fallspeed parameters
938 : ! includes air density correction factor to the
939 : ! power of 0.54 following Heymsfield and Bansemer 2007
940 :
941 0 : rhof(i,k)=(rhosu/rho(i,k))**0.54_r8
942 :
943 0 : arn(i,k)=ar*rhof(i,k)
944 0 : asn(i,k)=as*rhof(i,k)
945 0 : acn(i,k)=ac*rhof(i,k)
946 0 : ain(i,k)=ai*rhof(i,k)
947 :
948 : ! get dz from dp and hydrostatic approx
949 : ! keep dz positive (define as layer k-1 - layer k)
950 :
951 0 : dz(i,k)= pdel(i,k)/(rho(i,k)*g)
952 :
953 : end do
954 : end do
955 :
956 : ! initialization
957 0 : qc(1:ncol,1:top_lev-1) = 0._r8
958 0 : qi(1:ncol,1:top_lev-1) = 0._r8
959 0 : nc(1:ncol,1:top_lev-1) = 0._r8
960 0 : ni(1:ncol,1:top_lev-1) = 0._r8
961 0 : t1(1:ncol,1:pver) = t(1:ncol,1:pver)
962 0 : q1(1:ncol,1:pver) = q(1:ncol,1:pver)
963 0 : qc1(1:ncol,1:pver) = qc(1:ncol,1:pver)
964 0 : qi1(1:ncol,1:pver) = qi(1:ncol,1:pver)
965 0 : nc1(1:ncol,1:pver) = nc(1:ncol,1:pver)
966 0 : ni1(1:ncol,1:pver) = ni(1:ncol,1:pver)
967 :
968 : ! initialize tendencies to zero
969 0 : tlat1(1:ncol,1:pver)=0._r8
970 0 : qvlat1(1:ncol,1:pver)=0._r8
971 0 : qctend1(1:ncol,1:pver)=0._r8
972 0 : qitend1(1:ncol,1:pver)=0._r8
973 0 : nctend1(1:ncol,1:pver)=0._r8
974 0 : nitend1(1:ncol,1:pver)=0._r8
975 :
976 : ! initialize precip output
977 0 : qrout(1:ncol,1:pver)=0._r8
978 0 : qsout(1:ncol,1:pver)=0._r8
979 0 : nrout(1:ncol,1:pver)=0._r8
980 0 : nsout(1:ncol,1:pver)=0._r8
981 0 : dsout(1:ncol,1:pver)=0._r8
982 :
983 0 : drout(1:ncol,1:pver)=0._r8
984 :
985 0 : reff_rain(1:ncol,1:pver)=0._r8
986 0 : reff_snow(1:ncol,1:pver)=0._r8
987 :
988 : ! initialize variables for trop_mozart
989 0 : nevapr(1:ncol,1:pver) = 0._r8
990 0 : nevapr2(1:ncol,1:pver) = 0._r8
991 0 : evapsnow(1:ncol,1:pver) = 0._r8
992 0 : prain(1:ncol,1:pver) = 0._r8
993 0 : prodsnow(1:ncol,1:pver) = 0._r8
994 0 : cmeout(1:ncol,1:pver) = 0._r8
995 :
996 0 : am_evp_st(1:ncol,1:pver) = 0._r8
997 :
998 : ! for refl calc
999 0 : rainrt1(1:ncol,1:pver) = 0._r8
1000 :
1001 : ! initialize precip fraction and output tendencies
1002 0 : cldmax(1:ncol,1:pver)=mincld
1003 :
1004 : !initialize aerosol number
1005 : ! naer2(1:ncol,1:pver,:)=0._r8
1006 0 : dum2l(1:ncol,1:pver)=0._r8
1007 0 : dum2i(1:ncol,1:pver)=0._r8
1008 :
1009 : ! initialize avg precip rate
1010 0 : prect1(1:ncol)=0._r8
1011 0 : preci1(1:ncol)=0._r8
1012 :
1013 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1014 : !Get humidity and saturation vapor pressures
1015 :
1016 0 : do k=top_lev,pver
1017 :
1018 0 : do i=1,ncol
1019 :
1020 : ! find wet bulk temperature and saturation value for provisional t and q without
1021 : ! condensation
1022 :
1023 0 : es(i) = svp_water(t(i,k))
1024 0 : qs(i) = svp_to_qsat(es(i), p(i,k))
1025 :
1026 : ! Prevents negative values.
1027 0 : if (qs(i) < 0.0_r8) then
1028 0 : qs(i) = 1.0_r8
1029 0 : es(i) = p(i,k)
1030 : end if
1031 :
1032 0 : esl(i,k)=svp_water(t(i,k))
1033 0 : esi(i,k)=svp_ice(t(i,k))
1034 :
1035 : ! hm fix, make sure when above freezing that esi=esl, not active yet
1036 0 : if (t(i,k).gt.tmelt)esi(i,k)=esl(i,k)
1037 :
1038 0 : relhum(i,k)=q(i,k)/qs(i)
1039 :
1040 : ! get cloud fraction, check for minimum
1041 :
1042 0 : cldm(i,k)=max(cldn(i,k),mincld)
1043 0 : cldmw(i,k)=max(cldn(i,k),mincld)
1044 :
1045 0 : icldm(i,k)=max(icecldf(i,k),mincld)
1046 0 : lcldm(i,k)=max(liqcldf(i,k),mincld)
1047 :
1048 : ! subcolumns, set cloud fraction variables to one
1049 : ! if cloud water or ice is present, if not present
1050 : ! set to mincld (mincld used instead of zero, to prevent
1051 : ! possible division by zero errors
1052 :
1053 0 : if (microp_uniform) then
1054 :
1055 0 : cldm(i,k)=mincld
1056 0 : cldmw(i,k)=mincld
1057 0 : icldm(i,k)=mincld
1058 0 : lcldm(i,k)=mincld
1059 :
1060 0 : if (qc(i,k).ge.qsmall) then
1061 0 : lcldm(i,k)=1._r8
1062 0 : cldm(i,k)=1._r8
1063 0 : cldmw(i,k)=1._r8
1064 : end if
1065 :
1066 0 : if (qi(i,k).ge.qsmall) then
1067 0 : cldm(i,k)=1._r8
1068 0 : icldm(i,k)=1._r8
1069 : end if
1070 :
1071 : end if ! sub-columns
1072 :
1073 : ! calculate nfice based on liquid and ice mmr (no rain and snow mmr available yet)
1074 :
1075 0 : nfice(i,k)=0._r8
1076 0 : dumfice=qc(i,k)+qi(i,k)
1077 0 : if (dumfice.gt.qsmall .and. qi(i,k).gt.qsmall) then
1078 0 : nfice(i,k)=qi(i,k)/dumfice
1079 : endif
1080 :
1081 0 : if (do_cldice .and. (t(i,k).lt.tmelt - 5._r8)) then
1082 :
1083 : ! if aerosols interact with ice set number of activated ice nuclei
1084 0 : dum2=naai(i,k)
1085 :
1086 0 : dumnnuc=(dum2-ni(i,k)/icldm(i,k))/deltat*icldm(i,k)
1087 0 : dumnnuc=max(dumnnuc,0._r8)
1088 : ! get provisional ni and qi after nucleation in order to calculate
1089 : ! Bergeron process below
1090 0 : ninew=ni(i,k)+dumnnuc*deltat
1091 0 : qinew=qi(i,k)+dumnnuc*deltat*mi0
1092 :
1093 : !T>268
1094 : else
1095 0 : ninew=ni(i,k)
1096 0 : qinew=qi(i,k)
1097 : end if
1098 :
1099 : ! Initialize CME components
1100 :
1101 0 : cme(i,k) = 0._r8
1102 0 : cmei(i,k)=0._r8
1103 :
1104 :
1105 : !-------------------------------------------------------------------
1106 : !Bergeron process
1107 :
1108 : ! make sure to initialize bergeron process to zero
1109 0 : berg(i,k)=0._r8
1110 0 : prd = 0._r8
1111 :
1112 : !condensation loop.
1113 :
1114 : ! get in-cloud qi and ni after nucleation
1115 0 : if (icldm(i,k) .gt. 0._r8) then
1116 0 : qiic(i,k)=qinew/icldm(i,k)
1117 0 : niic(i,k)=ninew/icldm(i,k)
1118 : else
1119 0 : qiic(i,k)=0._r8
1120 0 : niic(i,k)=0._r8
1121 : endif
1122 :
1123 0 : if (nicons) then
1124 0 : niic(i,k) = ninst/rho(i,k)
1125 : end if
1126 :
1127 : !if T < 0 C then bergeron.
1128 0 : if (do_cldice .and. (t(i,k).lt.273.15_r8)) then
1129 :
1130 : !if ice exists
1131 0 : if (qi(i,k).gt.qsmall) then
1132 :
1133 0 : bergtsf = 0._r8 ! bergeron time scale (fraction of timestep)
1134 :
1135 0 : qvi = svp_to_qsat(esi(i,k), p(i,k))
1136 0 : qvl = svp_to_qsat(esl(i,k), p(i,k))
1137 :
1138 0 : dqsidt = xxls*qvi/(rv*t(i,k)**2)
1139 0 : abi = 1._r8+dqsidt*xxls/cpp
1140 :
1141 : ! get ice size distribution parameters
1142 :
1143 0 : if (qiic(i,k).ge.qsmall) then
1144 : lami(k) = (cons1*ci* &
1145 0 : niic(i,k)/qiic(i,k))**(1._r8/di)
1146 0 : n0i(k) = niic(i,k)*lami(k)
1147 :
1148 : ! check for slope
1149 : ! adjust vars
1150 0 : if (lami(k).lt.lammini) then
1151 :
1152 0 : lami(k) = lammini
1153 0 : n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
1154 0 : else if (lami(k).gt.lammaxi) then
1155 0 : lami(k) = lammaxi
1156 0 : n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
1157 : end if
1158 :
1159 0 : epsi = 2._r8*pi*n0i(k)*rho(i,k)*Dv(i,k)/(lami(k)*lami(k))
1160 :
1161 : !if liquid exists
1162 0 : if (qc(i,k).gt. qsmall) then
1163 :
1164 : !begin bergeron process
1165 : ! do bergeron (vapor deposition with RHw=1)
1166 : ! code to find berg (a rate) goes here
1167 :
1168 : ! calculate Bergeron process
1169 :
1170 0 : prd = epsi*(qvl-qvi)/abi
1171 :
1172 : else
1173 : prd = 0._r8
1174 : end if
1175 :
1176 : ! multiply by cloud fraction
1177 :
1178 0 : prd = prd*min(icldm(i,k),lcldm(i,k))
1179 :
1180 : ! transfer of existing cloud liquid to ice
1181 :
1182 0 : berg(i,k)=max(0._r8,prd)
1183 :
1184 : end if !end liquid exists bergeron
1185 :
1186 0 : if (berg(i,k).gt.0._r8) then
1187 0 : bergtsf=max(0._r8,(qc(i,k)/berg(i,k))/deltat)
1188 :
1189 0 : if(bergtsf.lt.1._r8) berg(i,k) = max(0._r8,qc(i,k)/deltat)
1190 :
1191 : endif
1192 :
1193 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1194 :
1195 0 : if (bergtsf.lt.1._r8.or.icldm(i,k).gt.lcldm(i,k)) then
1196 :
1197 0 : if (qiic(i,k).ge.qsmall) then
1198 :
1199 : ! first case is for case when liquid water is present, but is completely depleted
1200 : ! in time step, i.e., bergrsf > 0 but < 1
1201 :
1202 0 : if (qc(i,k).ge.qsmall) then
1203 0 : rhin = (1.0_r8 + relhum(i,k)) / 2._r8
1204 0 : if ((rhin*esl(i,k)/esi(i,k)) > 1._r8) then
1205 0 : prd = epsi*(rhin*qvl-qvi)/abi
1206 :
1207 : ! multiply by cloud fraction assuming liquid/ice maximum overlap
1208 0 : prd = prd*min(icldm(i,k),lcldm(i,k))
1209 :
1210 : ! add to cmei
1211 0 : cmei(i,k) = cmei(i,k) + (prd * (1._r8- bergtsf))
1212 :
1213 : end if ! rhin
1214 : end if ! qc > qsmall
1215 :
1216 : ! second case is for pure ice cloud, either no liquid, or icldm > lcldm
1217 :
1218 0 : if (qc(i,k).lt.qsmall.or.icldm(i,k).gt.lcldm(i,k)) then
1219 :
1220 : ! note: for case of no liquid, need to set liquid cloud fraction to zero
1221 : ! store liquid cloud fraction in 'dum'
1222 :
1223 0 : if (qc(i,k).lt.qsmall) then
1224 : dum=0._r8
1225 : else
1226 0 : dum=lcldm(i,k)
1227 : end if
1228 :
1229 : ! set RH to grid-mean value for pure ice cloud
1230 0 : rhin = relhum(i,k)
1231 :
1232 0 : if ((rhin*esl(i,k)/esi(i,k)) > 1._r8) then
1233 :
1234 0 : prd = epsi*(rhin*qvl-qvi)/abi
1235 :
1236 : ! multiply by relevant cloud fraction for pure ice cloud
1237 : ! assuming maximum overlap of liquid/ice
1238 0 : prd = prd*max((icldm(i,k)-dum),0._r8)
1239 0 : cmei(i,k) = cmei(i,k) + prd
1240 :
1241 : end if ! rhin
1242 : end if ! qc or icldm > lcldm
1243 : end if ! qiic
1244 : end if ! bergtsf or icldm > lcldm
1245 :
1246 : ! if deposition, it should not reduce grid mean rhi below 1.0
1247 0 : if(cmei(i,k) > 0.0_r8 .and. (relhum(i,k)*esl(i,k)/esi(i,k)) > 1._r8 ) &
1248 0 : cmei(i,k)=min(cmei(i,k),(q(i,k)-qs(i)*esi(i,k)/esl(i,k))/abi/deltat)
1249 :
1250 : end if !end ice exists loop
1251 : !this ends temperature < 0. loop
1252 :
1253 : !-------------------------------------------------------------------
1254 : end if !
1255 : !..............................................................
1256 :
1257 : ! evaporation should not exceed available water
1258 :
1259 0 : if ((-berg(i,k)).lt.-qc(i,k)/deltat) berg(i,k) = max(qc(i,k)/deltat,0._r8)
1260 :
1261 : !sublimation process...
1262 0 : if (do_cldice .and. ((relhum(i,k)*esl(i,k)/esi(i,k)).lt.1._r8 .and. qiic(i,k).ge.qsmall )) then
1263 :
1264 0 : qvi = svp_to_qsat(esi(i,k), p(i,k))
1265 0 : qvl = svp_to_qsat(esl(i,k), p(i,k))
1266 0 : dqsidt = xxls*qvi/(rv*t(i,k)**2)
1267 0 : abi = 1._r8+dqsidt*xxls/cpp
1268 :
1269 : ! get ice size distribution parameters
1270 :
1271 : lami(k) = (cons1*ci* &
1272 0 : niic(i,k)/qiic(i,k))**(1._r8/di)
1273 0 : n0i(k) = niic(i,k)*lami(k)
1274 :
1275 : ! check for slope
1276 : ! adjust vars
1277 0 : if (lami(k).lt.lammini) then
1278 :
1279 0 : lami(k) = lammini
1280 0 : n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
1281 0 : else if (lami(k).gt.lammaxi) then
1282 0 : lami(k) = lammaxi
1283 0 : n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
1284 : end if
1285 :
1286 0 : epsi = 2._r8*pi*n0i(k)*rho(i,k)*Dv(i,k)/(lami(k)*lami(k))
1287 :
1288 : ! modify for ice fraction below
1289 0 : prd = epsi*(relhum(i,k)*qvl-qvi)/abi * icldm(i,k)
1290 0 : cmei(i,k)=min(prd,0._r8)
1291 :
1292 : endif
1293 :
1294 : ! sublimation should not exceed available ice
1295 0 : if (cmei(i,k).lt.-qi(i,k)/deltat) cmei(i,k)=-qi(i,k)/deltat
1296 :
1297 : ! sublimation should not increase grid mean rhi above 1.0
1298 0 : if(cmei(i,k) < 0.0_r8 .and. (relhum(i,k)*esl(i,k)/esi(i,k)) < 1._r8 ) &
1299 0 : cmei(i,k)=min(0._r8,max(cmei(i,k),(q(i,k)-qs(i)*esi(i,k)/esl(i,k))/abi/deltat))
1300 :
1301 : ! limit cmei due for roundoff error
1302 :
1303 0 : cmei(i,k)=cmei(i,k)*omsm
1304 :
1305 : ! conditional for ice nucleation
1306 0 : if (do_cldice .and. (t(i,k).lt.(tmelt - 5._r8))) then
1307 :
1308 : ! using Liu et al. (2007) ice nucleation with hooks into simulated aerosol
1309 : ! ice nucleation rate (dum2) has already been calculated and read in (naai)
1310 :
1311 0 : dum2i(i,k)=naai(i,k)
1312 : else
1313 0 : dum2i(i,k)=0._r8
1314 : end if
1315 :
1316 : end do ! i loop
1317 : end do ! k loop
1318 :
1319 :
1320 : !! initialize sub-step precip flux variables
1321 0 : do i=1,ncol
1322 : !! flux is zero at top interface, so these should stay as 0.
1323 0 : rflx1(i,1)=0._r8
1324 0 : sflx1(i,1)=0._r8
1325 0 : do k=top_lev,pver
1326 :
1327 : ! initialize normal and sub-step precip flux variables
1328 0 : rflx1(i,k+1)=0._r8
1329 0 : sflx1(i,k+1)=0._r8
1330 : end do ! i loop
1331 : end do ! k loop
1332 : !! initialize final precip flux variables.
1333 0 : do i=1,ncol
1334 : !! flux is zero at top interface, so these should stay as 0.
1335 0 : rflx(i,1)=0._r8
1336 0 : sflx(i,1)=0._r8
1337 0 : do k=top_lev,pver
1338 : ! initialize normal and sub-step precip flux variables
1339 0 : rflx(i,k+1)=0._r8
1340 0 : sflx(i,k+1)=0._r8
1341 : end do ! i loop
1342 : end do ! k loop
1343 :
1344 0 : do i=1,ncol
1345 0 : ltrue(i)=0
1346 0 : do k=top_lev,pver
1347 : ! skip microphysical calculations if no cloud water
1348 :
1349 0 : if (qc(i,k).ge.qsmall.or.qi(i,k).ge.qsmall.or.cmei(i,k).ge.qsmall) ltrue(i)=1
1350 : end do
1351 : end do
1352 :
1353 : ! assign number of sub-steps to iter
1354 : ! use 2 sub-steps, following tests described in MG2008
1355 : iter = 2
1356 :
1357 : ! get sub-step time step
1358 0 : deltat=deltat/real(iter)
1359 :
1360 : ! since activation/nucleation processes are fast, need to take into account
1361 : ! factor mtime = mixing timescale in cloud / model time step
1362 : ! mixing time can be interpreted as cloud depth divided by sub-grid vertical velocity
1363 : ! for now mixing timescale is assumed to be 1 timestep for modal aerosols, 20 min bulk
1364 :
1365 : ! note: mtime for bulk aerosols was set to: mtime=deltat/1200._r8
1366 :
1367 : mtime=1._r8
1368 0 : rate1ord_cw2pr_st(:,:)=0._r8 ! rce 2010/05/01
1369 :
1370 : !!!! skip calculations if no cloud water
1371 0 : do i=1,ncol
1372 0 : if (ltrue(i).eq.0) then
1373 0 : tlat(i,1:pver)=0._r8
1374 0 : qvlat(i,1:pver)=0._r8
1375 0 : qctend(i,1:pver)=0._r8
1376 0 : qitend(i,1:pver)=0._r8
1377 0 : qnitend(i,1:pver)=0._r8
1378 0 : qrtend(i,1:pver)=0._r8
1379 0 : nctend(i,1:pver)=0._r8
1380 0 : nitend(i,1:pver)=0._r8
1381 0 : nrtend(i,1:pver)=0._r8
1382 0 : nstend(i,1:pver)=0._r8
1383 0 : prect(i)=0._r8
1384 0 : preci(i)=0._r8
1385 0 : rflx(i,1:pver+1)=0._r8
1386 0 : sflx(i,1:pver+1)=0._r8
1387 0 : qniic(i,1:pver)=0._r8
1388 0 : qric(i,1:pver)=0._r8
1389 0 : nsic(i,1:pver)=0._r8
1390 0 : nric(i,1:pver)=0._r8
1391 0 : rainrt(i,1:pver)=0._r8
1392 : goto 300
1393 : end if
1394 :
1395 0 : qcsinksum_rate1ord(1:pver)=0._r8
1396 0 : qcsum_rate1ord(1:pver)=0._r8
1397 :
1398 :
1399 : !!!!!!!!! begin sub-step!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1400 : !.....................................................................................................
1401 0 : do it=1,iter
1402 :
1403 : ! initialize sub-step microphysical tendencies
1404 :
1405 0 : tlat(i,1:pver)=0._r8
1406 0 : qvlat(i,1:pver)=0._r8
1407 0 : qctend(i,1:pver)=0._r8
1408 0 : qitend(i,1:pver)=0._r8
1409 0 : qnitend(i,1:pver)=0._r8
1410 0 : qrtend(i,1:pver)=0._r8
1411 0 : nctend(i,1:pver)=0._r8
1412 0 : nitend(i,1:pver)=0._r8
1413 0 : nrtend(i,1:pver)=0._r8
1414 0 : nstend(i,1:pver)=0._r8
1415 :
1416 : ! initialize diagnostic precipitation to zero
1417 :
1418 0 : qniic(i,1:pver)=0._r8
1419 0 : qric(i,1:pver)=0._r8
1420 0 : nsic(i,1:pver)=0._r8
1421 0 : nric(i,1:pver)=0._r8
1422 :
1423 0 : rainrt(i,1:pver)=0._r8
1424 :
1425 :
1426 : ! begin new i,k loop, calculate new cldmax after adjustment to cldm above
1427 :
1428 : ! initialize vertically-integrated rain and snow tendencies
1429 :
1430 0 : qrtot = 0._r8
1431 0 : nrtot = 0._r8
1432 0 : qstot = 0._r8
1433 0 : nstot = 0._r8
1434 :
1435 : ! initialize precip at surface
1436 :
1437 0 : prect(i)=0._r8
1438 0 : preci(i)=0._r8
1439 :
1440 : ! initialize fluxes
1441 0 : rflx(i,1:pver+1)=0._r8
1442 0 : sflx(i,1:pver+1)=0._r8
1443 :
1444 0 : do k=top_lev,pver
1445 :
1446 0 : qcvar=relvar(i,k)
1447 0 : cons2=gamma(qcvar+2.47_r8)
1448 0 : cons3=gamma(qcvar)
1449 0 : cons9=gamma(qcvar+2._r8)
1450 0 : cons10=gamma(qcvar+1._r8)
1451 0 : cons12=gamma(qcvar+1.15_r8)
1452 0 : cons15=gamma(qcvar+bc/3._r8)
1453 0 : cons18=qcvar**2.47_r8
1454 0 : cons19=qcvar**2
1455 0 : cons20=qcvar**1.15_r8
1456 :
1457 : ! set cwml and cwmi to current qc and qi
1458 :
1459 0 : cwml(i,k)=qc(i,k)
1460 0 : cwmi(i,k)=qi(i,k)
1461 :
1462 : ! initialize precip fallspeeds to zero
1463 :
1464 0 : ums(k)=0._r8
1465 0 : uns(k)=0._r8
1466 0 : umr(k)=0._r8
1467 0 : unr(k)=0._r8
1468 :
1469 : ! calculate precip fraction based on maximum overlap assumption
1470 :
1471 : ! for sub-columns cldm has already been set to 1 if cloud
1472 : ! water or ice is present, so cldmax will be correctly set below
1473 : ! and nothing extra needs to be done here
1474 :
1475 0 : if (k.eq.top_lev) then
1476 0 : cldmax(i,k)=cldm(i,k)
1477 : else
1478 : ! if rain or snow mix ratio is smaller than
1479 : ! threshold, then set cldmax to cloud fraction at current level
1480 :
1481 0 : if (do_clubb_sgs) then
1482 0 : if (qc(i,k).ge.qsmall.or.qi(i,k).ge.qsmall) then
1483 0 : cldmax(i,k)=cldm(i,k)
1484 : else
1485 0 : cldmax(i,k)=cldmax(i,k-1)
1486 : end if
1487 : else
1488 :
1489 0 : if (qric(i,k-1).ge.qsmall.or.qniic(i,k-1).ge.qsmall) then
1490 0 : cldmax(i,k)=max(cldmax(i,k-1),cldm(i,k))
1491 : else
1492 0 : cldmax(i,k)=cldm(i,k)
1493 : end if
1494 : endif
1495 : end if
1496 :
1497 : ! decrease in number concentration due to sublimation/evap
1498 : ! divide by cloud fraction to get in-cloud decrease
1499 : ! don't reduce Nc due to bergeron process
1500 :
1501 0 : if (cmei(i,k) < 0._r8 .and. qi(i,k) > qsmall .and. cldm(i,k) > mincld) then
1502 0 : nsubi(k)=cmei(i,k)/qi(i,k)*ni(i,k)/cldm(i,k)
1503 : else
1504 0 : nsubi(k)=0._r8
1505 : end if
1506 0 : nsubc(k)=0._r8
1507 :
1508 :
1509 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1510 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1511 :
1512 : ! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5%
1513 :
1514 0 : if (do_cldice .and. dum2i(i,k).gt.0._r8.and.t(i,k).lt.(tmelt - 5._r8).and. &
1515 : relhum(i,k)*esl(i,k)/esi(i,k).gt. rhmini+0.05_r8) then
1516 :
1517 : !if NCAI > 0. then set numice = ncai (as before)
1518 : !note: this is gridbox averaged
1519 :
1520 0 : nnuccd(k)=(dum2i(i,k)-ni(i,k)/icldm(i,k))/deltat*icldm(i,k)
1521 0 : nnuccd(k)=max(nnuccd(k),0._r8)
1522 0 : nimax = dum2i(i,k)*icldm(i,k)
1523 :
1524 : !Calc mass of new particles using new crystal mass...
1525 : !also this will be multiplied by mtime as nnuccd is...
1526 :
1527 0 : mnuccd(k) = nnuccd(k) * mi0
1528 :
1529 : ! add mnuccd to cmei....
1530 0 : cmei(i,k)= cmei(i,k) + mnuccd(k) * mtime
1531 :
1532 : ! limit cmei
1533 :
1534 0 : qvi = svp_to_qsat(esi(i,k), p(i,k))
1535 0 : dqsidt = xxls*qvi/(rv*t(i,k)**2)
1536 0 : abi = 1._r8+dqsidt*xxls/cpp
1537 0 : cmei(i,k)=min(cmei(i,k),(q(i,k)-qvi)/abi/deltat)
1538 :
1539 : ! limit for roundoff error
1540 0 : cmei(i,k)=cmei(i,k)*omsm
1541 :
1542 : else
1543 0 : nnuccd(k)=0._r8
1544 0 : nimax = 0._r8
1545 0 : mnuccd(k) = 0._r8
1546 : end if
1547 :
1548 : !c............................................................................
1549 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1550 : ! obtain in-cloud values of cloud water/ice mixing ratios and number concentrations
1551 : ! for microphysical process calculations
1552 : ! units are kg/kg for mixing ratio, 1/kg for number conc
1553 :
1554 : ! limit in-cloud values to 0.005 kg/kg
1555 :
1556 0 : qcic(i,k)=min(cwml(i,k)/lcldm(i,k),5.e-3_r8)
1557 0 : qiic(i,k)=min(cwmi(i,k)/icldm(i,k),5.e-3_r8)
1558 0 : ncic(i,k)=max(nc(i,k)/lcldm(i,k),0._r8)
1559 0 : niic(i,k)=max(ni(i,k)/icldm(i,k),0._r8)
1560 :
1561 0 : if (nccons) then
1562 0 : ncic(i,k) = ncnst/rho(i,k)
1563 : end if
1564 0 : if (nicons) then
1565 0 : niic(i,k) = ninst/rho(i,k)
1566 : end if
1567 :
1568 0 : if (qc(i,k) - berg(i,k)*deltat.lt.qsmall) then
1569 0 : qcic(i,k)=0._r8
1570 0 : ncic(i,k)=0._r8
1571 0 : if (qc(i,k)-berg(i,k)*deltat.lt.0._r8) then
1572 0 : berg(i,k)=qc(i,k)/deltat*omsm
1573 : end if
1574 : end if
1575 :
1576 0 : if (do_cldice .and. qi(i,k)+(cmei(i,k)+berg(i,k))*deltat.lt.qsmall) then
1577 0 : qiic(i,k)=0._r8
1578 0 : niic(i,k)=0._r8
1579 0 : if (qi(i,k)+(cmei(i,k)+berg(i,k))*deltat.lt.0._r8) then
1580 0 : cmei(i,k)=(-qi(i,k)/deltat-berg(i,k))*omsm
1581 : end if
1582 : end if
1583 :
1584 : ! add to cme output
1585 :
1586 0 : cmeout(i,k) = cmeout(i,k)+cmei(i,k)
1587 :
1588 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1589 : ! droplet activation
1590 : ! calculate potential for droplet activation if cloud water is present
1591 : ! formulation from Abdul-Razzak and Ghan (2000) and Abdul-Razzak et al. (1998), AR98
1592 : ! number tendency (npccnin) is read in from companion routine
1593 :
1594 : ! assume aerosols already activated are equal to number of existing droplets for simplicity
1595 : ! multiply by cloud fraction to obtain grid-average tendency
1596 :
1597 0 : if (qcic(i,k).ge.qsmall) then
1598 0 : npccn(k) = max(0._r8,npccnin(i,k))
1599 0 : dum2l(i,k)=(nc(i,k)+npccn(k)*deltat)/lcldm(i,k)
1600 0 : dum2l(i,k)=max(dum2l(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3
1601 0 : ncmax = dum2l(i,k)*lcldm(i,k)
1602 : else
1603 0 : npccn(k)=0._r8
1604 0 : dum2l(i,k)=0._r8
1605 0 : ncmax = 0._r8
1606 : end if
1607 :
1608 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1609 : ! get size distribution parameters based on in-cloud cloud water/ice
1610 : ! these calculations also ensure consistency between number and mixing ratio
1611 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1612 :
1613 : !......................................................................
1614 : ! cloud ice
1615 :
1616 0 : if (qiic(i,k).ge.qsmall) then
1617 :
1618 : ! add upper limit to in-cloud number concentration to prevent numerical error
1619 0 : niic(i,k)=min(niic(i,k),qiic(i,k)*1.e20_r8)
1620 :
1621 0 : lami(k) = (cons1*ci*niic(i,k)/qiic(i,k))**(1._r8/di)
1622 0 : n0i(k) = niic(i,k)*lami(k)
1623 :
1624 : ! check for slope
1625 : ! adjust vars
1626 :
1627 0 : if (lami(k).lt.lammini) then
1628 :
1629 0 : lami(k) = lammini
1630 0 : n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
1631 0 : niic(i,k) = n0i(k)/lami(k)
1632 0 : else if (lami(k).gt.lammaxi) then
1633 0 : lami(k) = lammaxi
1634 0 : n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
1635 0 : niic(i,k) = n0i(k)/lami(k)
1636 : end if
1637 :
1638 : else
1639 0 : lami(k) = 0._r8
1640 0 : n0i(k) = 0._r8
1641 : end if
1642 :
1643 0 : if (qcic(i,k).ge.qsmall) then
1644 :
1645 : ! add upper limit to in-cloud number concentration to prevent numerical error
1646 0 : ncic(i,k)=min(ncic(i,k),qcic(i,k)*1.e20_r8)
1647 :
1648 0 : ncic(i,k)=max(ncic(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm
1649 :
1650 : ! get pgam from fit to observations of martin et al. 1994
1651 :
1652 0 : pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
1653 0 : pgam(k)=1._r8/(pgam(k)**2)-1._r8
1654 0 : pgam(k)=max(pgam(k),2._r8)
1655 0 : pgam(k)=min(pgam(k),15._r8)
1656 :
1657 : ! calculate lamc
1658 :
1659 : lamc(k) = (pi/6._r8*rhow*ncic(i,k)*gamma(pgam(k)+4._r8)/ &
1660 0 : (qcic(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
1661 :
1662 : ! lammin, 50 micron diameter max mean size
1663 :
1664 0 : lammin = (pgam(k)+1._r8)/50.e-6_r8
1665 0 : lammax = (pgam(k)+1._r8)/2.e-6_r8
1666 :
1667 0 : if (lamc(k).lt.lammin) then
1668 0 : lamc(k) = lammin
1669 : ncic(i,k) = 6._r8*lamc(k)**3*qcic(i,k)* &
1670 0 : gamma(pgam(k)+1._r8)/(pi*rhow*gamma(pgam(k)+4._r8))
1671 0 : else if (lamc(k).gt.lammax) then
1672 0 : lamc(k) = lammax
1673 : ncic(i,k) = 6._r8*lamc(k)**3*qcic(i,k)* &
1674 0 : gamma(pgam(k)+1._r8)/(pi*rhow*gamma(pgam(k)+4._r8))
1675 : end if
1676 :
1677 : ! parameter to calculate droplet freezing
1678 :
1679 0 : cdist1(k) = ncic(i,k)/gamma(pgam(k)+1._r8)
1680 :
1681 : else
1682 0 : lamc(k) = 0._r8
1683 0 : cdist1(k) = 0._r8
1684 : end if
1685 :
1686 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1687 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1688 : ! begin micropysical process calculations
1689 : !.................................................................
1690 : ! autoconversion of cloud liquid water to rain
1691 : ! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
1692 : ! minimum qc of 1 x 10^-8 prevents floating point error
1693 :
1694 0 : if (qcic(i,k).ge.1.e-8_r8) then
1695 :
1696 : ! nprc is increase in rain number conc due to autoconversion
1697 : ! nprc1 is decrease in cloud droplet conc due to autoconversion
1698 :
1699 : ! assume exponential sub-grid distribution of qc, resulting in additional
1700 : ! factor related to qcvar below
1701 :
1702 : ! hm switch for sub-columns, don't include sub-grid qc
1703 0 : if (microp_uniform) then
1704 :
1705 : prc(k) = 1350._r8*qcic(i,k)**2.47_r8* &
1706 0 : (ncic(i,k)/1.e6_r8*rho(i,k))**(-1.79_r8)
1707 0 : nprc(k) = prc(k)/(4._r8/3._r8*pi*rhow*(25.e-6_r8)**3)
1708 0 : nprc1(k) = prc(k)/(qcic(i,k)/ncic(i,k))
1709 :
1710 : else
1711 :
1712 : prc(k) = cons2/(cons3*cons18)*1350._r8*qcic(i,k)**2.47_r8* &
1713 0 : (ncic(i,k)/1.e6_r8*rho(i,k))**(-1.79_r8)
1714 0 : nprc(k) = prc(k)/cons22
1715 0 : nprc1(k) = prc(k)/(qcic(i,k)/ncic(i,k))
1716 :
1717 : end if ! sub-column switch
1718 :
1719 : else
1720 0 : prc(k)=0._r8
1721 0 : nprc(k)=0._r8
1722 0 : nprc1(k)=0._r8
1723 : end if
1724 :
1725 : ! add autoconversion to precip from above to get provisional rain mixing ratio
1726 : ! and number concentration (qric and nric)
1727 :
1728 : ! 0.45 m/s is fallspeed of new rain drop (80 micron diameter)
1729 :
1730 0 : dum=0.45_r8
1731 0 : dum1=0.45_r8
1732 :
1733 0 : if (k.eq.top_lev) then
1734 0 : qric(i,k)=prc(k)*lcldm(i,k)*dz(i,k)/cldmax(i,k)/dum
1735 0 : nric(i,k)=nprc(k)*lcldm(i,k)*dz(i,k)/cldmax(i,k)/dum
1736 : else
1737 0 : if (qric(i,k-1).ge.qsmall) then
1738 0 : dum=umr(k-1)
1739 0 : dum1=unr(k-1)
1740 : end if
1741 :
1742 : ! no autoconversion of rain number if rain/snow falling from above
1743 : ! this assumes that new drizzle drops formed by autoconversion are rapidly collected
1744 : ! by the existing rain/snow particles from above
1745 :
1746 0 : if (qric(i,k-1).ge.1.e-9_r8.or.qniic(i,k-1).ge.1.e-9_r8) then
1747 0 : nprc(k)=0._r8
1748 : end if
1749 :
1750 : qric(i,k) = (rho(i,k-1)*umr(k-1)*qric(i,k-1)*cldmax(i,k-1)+ &
1751 : (rho(i,k)*dz(i,k)*((pra(k-1)+prc(k))*lcldm(i,k)+(pre(k-1)-pracs(k-1)-mnuccr(k-1))*cldmax(i,k))))&
1752 0 : /(dum*rho(i,k)*cldmax(i,k))
1753 : nric(i,k) = (rho(i,k-1)*unr(k-1)*nric(i,k-1)*cldmax(i,k-1)+ &
1754 : (rho(i,k)*dz(i,k)*(nprc(k)*lcldm(i,k)+(nsubr(k-1)-npracs(k-1)-nnuccr(k-1)+nragg(k-1))*cldmax(i,k))))&
1755 0 : /(dum1*rho(i,k)*cldmax(i,k))
1756 :
1757 : end if
1758 :
1759 : !.......................................................................
1760 : ! Autoconversion of cloud ice to snow
1761 : ! similar to Ferrier (1994)
1762 :
1763 0 : if (do_cldice) then
1764 0 : if (t(i,k).le.273.15_r8.and.qiic(i,k).ge.qsmall) then
1765 :
1766 : ! note: assumes autoconversion timescale of 180 sec
1767 :
1768 0 : nprci(k) = n0i(k)/(lami(k)*180._r8)*exp(-lami(k)*dcs)
1769 :
1770 : prci(k) = pi*rhoi*n0i(k)/(6._r8*180._r8)* &
1771 : (cons23/lami(k)+3._r8*cons24/lami(k)**2+ &
1772 0 : 6._r8*dcs/lami(k)**3+6._r8/lami(k)**4)*exp(-lami(k)*dcs)
1773 : else
1774 0 : prci(k)=0._r8
1775 0 : nprci(k)=0._r8
1776 : end if
1777 : else
1778 : ! Add in the particles that we have already converted to snow, and
1779 : ! don't do any further autoconversion of ice.
1780 0 : prci(k) = tnd_qsnow(i, k) / cldm(i,k)
1781 0 : nprci(k) = tnd_nsnow(i, k) / cldm(i,k)
1782 : end if
1783 :
1784 : ! add autoconversion to flux from level above to get provisional snow mixing ratio
1785 : ! and number concentration (qniic and nsic)
1786 :
1787 0 : dum=(asn(i,k)*cons25)
1788 0 : dum1=(asn(i,k)*cons25)
1789 :
1790 0 : if (k.eq.top_lev) then
1791 0 : qniic(i,k)=prci(k)*icldm(i,k)*dz(i,k)/cldmax(i,k)/dum
1792 0 : nsic(i,k)=nprci(k)*icldm(i,k)*dz(i,k)/cldmax(i,k)/dum
1793 : else
1794 0 : if (qniic(i,k-1).ge.qsmall) then
1795 0 : dum=ums(k-1)
1796 0 : dum1=uns(k-1)
1797 : end if
1798 :
1799 : qniic(i,k) = (rho(i,k-1)*ums(k-1)*qniic(i,k-1)*cldmax(i,k-1)+ &
1800 : (rho(i,k)*dz(i,k)*((prci(k)+prai(k-1)+psacws(k-1)+bergs(k-1))*icldm(i,k)+(prds(k-1)+ &
1801 : pracs(k-1)+mnuccr(k-1))*cldmax(i,k))))&
1802 0 : /(dum*rho(i,k)*cldmax(i,k))
1803 :
1804 : nsic(i,k) = (rho(i,k-1)*uns(k-1)*nsic(i,k-1)*cldmax(i,k-1)+ &
1805 : (rho(i,k)*dz(i,k)*(nprci(k)*icldm(i,k)+(nsubs(k-1)+nsagg(k-1)+nnuccr(k-1))*cldmax(i,k))))&
1806 0 : /(dum1*rho(i,k)*cldmax(i,k))
1807 :
1808 : end if
1809 :
1810 : ! if precip mix ratio is zero so should number concentration
1811 :
1812 0 : if (qniic(i,k).lt.qsmall) then
1813 0 : qniic(i,k)=0._r8
1814 0 : nsic(i,k)=0._r8
1815 : end if
1816 :
1817 0 : if (qric(i,k).lt.qsmall) then
1818 0 : qric(i,k)=0._r8
1819 0 : nric(i,k)=0._r8
1820 : end if
1821 :
1822 : ! make sure number concentration is a positive number to avoid
1823 : ! taking root of negative later
1824 :
1825 0 : nric(i,k)=max(nric(i,k),0._r8)
1826 0 : nsic(i,k)=max(nsic(i,k),0._r8)
1827 :
1828 : !.......................................................................
1829 : ! get size distribution parameters for precip
1830 : !......................................................................
1831 : ! rain
1832 :
1833 0 : if (qric(i,k).ge.qsmall) then
1834 0 : lamr(k) = (pi*rhow*nric(i,k)/qric(i,k))**(1._r8/3._r8)
1835 0 : n0r(k) = nric(i,k)*lamr(k)
1836 :
1837 : ! check for slope
1838 : ! adjust vars
1839 :
1840 0 : if (lamr(k).lt.lamminr) then
1841 :
1842 0 : lamr(k) = lamminr
1843 :
1844 0 : n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
1845 0 : nric(i,k) = n0r(k)/lamr(k)
1846 0 : else if (lamr(k).gt.lammaxr) then
1847 0 : lamr(k) = lammaxr
1848 0 : n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
1849 0 : nric(i,k) = n0r(k)/lamr(k)
1850 : end if
1851 :
1852 : ! provisional rain number and mass weighted mean fallspeed (m/s)
1853 :
1854 0 : unr(k) = min(arn(i,k)*cons4/lamr(k)**br,9.1_r8*rhof(i,k))
1855 0 : umr(k) = min(arn(i,k)*cons5/(6._r8*lamr(k)**br),9.1_r8*rhof(i,k))
1856 :
1857 : else
1858 0 : lamr(k) = 0._r8
1859 0 : n0r(k) = 0._r8
1860 0 : umr(k) = 0._r8
1861 0 : unr(k) = 0._r8
1862 : end if
1863 :
1864 : !......................................................................
1865 : ! snow
1866 :
1867 0 : if (qniic(i,k).ge.qsmall) then
1868 0 : lams(k) = (cons6*cs*nsic(i,k)/qniic(i,k))**(1._r8/ds)
1869 0 : n0s(k) = nsic(i,k)*lams(k)
1870 :
1871 : ! check for slope
1872 : ! adjust vars
1873 :
1874 0 : if (lams(k).lt.lammins) then
1875 0 : lams(k) = lammins
1876 0 : n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
1877 0 : nsic(i,k) = n0s(k)/lams(k)
1878 :
1879 0 : else if (lams(k).gt.lammaxs) then
1880 0 : lams(k) = lammaxs
1881 0 : n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
1882 0 : nsic(i,k) = n0s(k)/lams(k)
1883 : end if
1884 :
1885 : ! provisional snow number and mass weighted mean fallspeed (m/s)
1886 :
1887 0 : ums(k) = min(asn(i,k)*cons8/(6._r8*lams(k)**bs),1.2_r8*rhof(i,k))
1888 0 : uns(k) = min(asn(i,k)*cons7/lams(k)**bs,1.2_r8*rhof(i,k))
1889 :
1890 : else
1891 0 : lams(k) = 0._r8
1892 0 : n0s(k) = 0._r8
1893 0 : ums(k) = 0._r8
1894 0 : uns(k) = 0._r8
1895 : end if
1896 :
1897 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1898 :
1899 : ! heterogeneous freezing of cloud water
1900 :
1901 0 : if (.not. use_hetfrz_classnuc) then
1902 :
1903 0 : if (do_cldice .and. qcic(i,k).ge.qsmall .and. t(i,k).lt.269.15_r8) then
1904 :
1905 : ! immersion freezing (Bigg, 1953)
1906 :
1907 :
1908 : ! subcolumns
1909 :
1910 0 : if (microp_uniform) then
1911 :
1912 : mnuccc(k) = &
1913 : pi*pi/36._r8*rhow* &
1914 : cdist1(k)*gamma(7._r8+pgam(k))* &
1915 : bimm*(exp(aimm*(273.15_r8-t(i,k)))-1._r8)/ &
1916 0 : lamc(k)**3/lamc(k)**3
1917 :
1918 : nnuccc(k) = &
1919 : pi/6._r8*cdist1(k)*gamma(pgam(k)+4._r8) &
1920 : *bimm* &
1921 0 : (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamc(k)**3
1922 :
1923 : else
1924 :
1925 : mnuccc(k) = cons9/(cons3*cons19)* &
1926 : pi*pi/36._r8*rhow* &
1927 : cdist1(k)*gamma(7._r8+pgam(k))* &
1928 : bimm*(exp(aimm*(273.15_r8-t(i,k)))-1._r8)/ &
1929 0 : lamc(k)**3/lamc(k)**3
1930 :
1931 : nnuccc(k) = cons10/(cons3*qcvar)* &
1932 : pi/6._r8*cdist1(k)*gamma(pgam(k)+4._r8) &
1933 : *bimm* &
1934 0 : (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamc(k)**3
1935 : end if ! sub-columns
1936 :
1937 :
1938 : ! contact freezing (-40<T<-3 C) (Young, 1974) with hooks into simulated dust
1939 : ! dust size and number in 4 bins are read in from companion routine
1940 :
1941 0 : tcnt=(270.16_r8-t(i,k))**1.3_r8
1942 0 : viscosity=1.8e-5_r8*(t(i,k)/298.0_r8)**0.85_r8 ! Viscosity (kg/m/s)
1943 : mfp=2.0_r8*viscosity/(p(i,k) & ! Mean free path (m)
1944 0 : *sqrt(8.0_r8*28.96e-3_r8/(pi*8.314409_r8*t(i,k))))
1945 :
1946 0 : nslip1=1.0_r8+(mfp/rndst(i,k,1))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,1)/mfp))))! Slip correction factor
1947 0 : nslip2=1.0_r8+(mfp/rndst(i,k,2))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,2)/mfp))))
1948 0 : nslip3=1.0_r8+(mfp/rndst(i,k,3))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,3)/mfp))))
1949 0 : nslip4=1.0_r8+(mfp/rndst(i,k,4))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,4)/mfp))))
1950 :
1951 0 : ndfaer1=1.381e-23_r8*t(i,k)*nslip1/(6._r8*pi*viscosity*rndst(i,k,1)) ! aerosol diffusivity (m2/s)
1952 0 : ndfaer2=1.381e-23_r8*t(i,k)*nslip2/(6._r8*pi*viscosity*rndst(i,k,2))
1953 0 : ndfaer3=1.381e-23_r8*t(i,k)*nslip3/(6._r8*pi*viscosity*rndst(i,k,3))
1954 0 : ndfaer4=1.381e-23_r8*t(i,k)*nslip4/(6._r8*pi*viscosity*rndst(i,k,4))
1955 :
1956 :
1957 0 : if (microp_uniform) then
1958 :
1959 : mnucct(k) = &
1960 : (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ &
1961 : ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*pi*pi/3._r8*rhow* &
1962 0 : cdist1(k)*gamma(pgam(k)+5._r8)/lamc(k)**4
1963 :
1964 : nnucct(k) = (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ &
1965 : ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*2._r8*pi* &
1966 0 : cdist1(k)*gamma(pgam(k)+2._r8)/lamc(k)
1967 :
1968 : else
1969 :
1970 : mnucct(k) = gamma(qcvar+4._r8/3._r8)/(cons3*qcvar**(4._r8/3._r8))* &
1971 : (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ &
1972 : ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*pi*pi/3._r8*rhow* &
1973 0 : cdist1(k)*gamma(pgam(k)+5._r8)/lamc(k)**4
1974 :
1975 : nnucct(k) = gamma(qcvar+1._r8/3._r8)/(cons3*qcvar**(1._r8/3._r8))* &
1976 : (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ &
1977 : ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*2._r8*pi* &
1978 0 : cdist1(k)*gamma(pgam(k)+2._r8)/lamc(k)
1979 :
1980 : end if ! sub-column switch
1981 :
1982 : ! make sure number of droplets frozen does not exceed available ice nuclei concentration
1983 : ! this prevents 'runaway' droplet freezing
1984 :
1985 0 : if (nnuccc(k)*lcldm(i,k).gt.nnuccd(k)) then
1986 0 : dum=(nnuccd(k)/(nnuccc(k)*lcldm(i,k)))
1987 : ! scale mixing ratio of droplet freezing with limit
1988 0 : mnuccc(k)=mnuccc(k)*dum
1989 0 : nnuccc(k)=nnuccd(k)/lcldm(i,k)
1990 : end if
1991 :
1992 : else
1993 0 : mnuccc(k)=0._r8
1994 0 : nnuccc(k)=0._r8
1995 0 : mnucct(k)=0._r8
1996 0 : nnucct(k)=0._r8
1997 : end if
1998 :
1999 : else
2000 0 : if (do_cldice .and. qcic(i,k) >= qsmall) then
2001 0 : con1 = 1._r8/(1.333_r8*pi)**0.333_r8
2002 0 : r3lx = con1*(rho(i,k)*qcic(i,k)/(rhow*max(ncic(i,k)*rho(i,k), 1.0e6_r8)))**0.333_r8 ! in m
2003 0 : r3lx = max(4.e-6_r8, r3lx)
2004 0 : mi0l = 4._r8/3._r8*pi*rhow*r3lx**3_r8
2005 :
2006 0 : nnuccc(k) = frzimm(i,k)*1.0e6_r8/rho(i,k)
2007 0 : mnuccc(k) = nnuccc(k)*mi0l
2008 :
2009 0 : nnucct(k) = frzcnt(i,k)*1.0e6_r8/rho(i,k)
2010 0 : mnucct(k) = nnucct(k)*mi0l
2011 :
2012 0 : nnudep(k) = frzdep(i,k)*1.0e6_r8/rho(i,k)
2013 0 : mnudep(k) = nnudep(k)*mi0
2014 : else
2015 0 : nnuccc(k) = 0._r8
2016 0 : mnuccc(k) = 0._r8
2017 :
2018 0 : nnucct(k) = 0._r8
2019 0 : mnucct(k) = 0._r8
2020 :
2021 0 : nnudep(k) = 0._r8
2022 0 : mnudep(k) = 0._r8
2023 : end if
2024 : endif
2025 :
2026 :
2027 : !.......................................................................
2028 : ! snow self-aggregation from passarelli, 1978, used by reisner, 1998
2029 : ! this is hard-wired for bs = 0.4 for now
2030 : ! ignore self-collection of cloud ice
2031 :
2032 0 : if (qniic(i,k).ge.qsmall .and. t(i,k).le.273.15_r8) then
2033 : nsagg(k) = -1108._r8*asn(i,k)*Eii* &
2034 : pi**((1._r8-bs)/3._r8)*rhosn**((-2._r8-bs)/3._r8)*rho(i,k)** &
2035 : ((2._r8+bs)/3._r8)*qniic(i,k)**((2._r8+bs)/3._r8)* &
2036 : (nsic(i,k)*rho(i,k))**((4._r8-bs)/3._r8)/ &
2037 0 : (4._r8*720._r8*rho(i,k))
2038 : else
2039 0 : nsagg(k)=0._r8
2040 : end if
2041 :
2042 : !.......................................................................
2043 : ! accretion of cloud droplets onto snow/graupel
2044 : ! here use continuous collection equation with
2045 : ! simple gravitational collection kernel
2046 : ! ignore collisions between droplets/cloud ice
2047 : ! since minimum size ice particle for accretion is 50 - 150 micron
2048 :
2049 : ! ignore collision of snow with droplets above freezing
2050 :
2051 0 : if (qniic(i,k).ge.qsmall .and. t(i,k).le.tmelt .and. &
2052 : qcic(i,k).ge.qsmall) then
2053 :
2054 : ! put in size dependent collection efficiency
2055 : ! mean diameter of snow is area-weighted, since
2056 : ! accretion is function of crystal geometric area
2057 : ! collection efficiency is approximation based on stoke's law (Thompson et al. 2004)
2058 :
2059 0 : dc0 = (pgam(k)+1._r8)/lamc(k)
2060 0 : ds0 = 1._r8/lams(k)
2061 0 : dum = dc0*dc0*uns(k)*rhow/(9._r8*mu(i,k)*ds0)
2062 0 : eci = dum*dum/((dum+0.4_r8)*(dum+0.4_r8))
2063 :
2064 0 : eci = max(eci,0._r8)
2065 0 : eci = min(eci,1._r8)
2066 :
2067 :
2068 : ! no impact of sub-grid distribution of qc since psacws
2069 : ! is linear in qc
2070 :
2071 : psacws(k) = pi/4._r8*asn(i,k)*qcic(i,k)*rho(i,k)* &
2072 : n0s(k)*Eci*cons11/ &
2073 0 : lams(k)**(bs+3._r8)
2074 : npsacws(k) = pi/4._r8*asn(i,k)*ncic(i,k)*rho(i,k)* &
2075 : n0s(k)*Eci*cons11/ &
2076 0 : lams(k)**(bs+3._r8)
2077 : else
2078 0 : psacws(k)=0._r8
2079 0 : npsacws(k)=0._r8
2080 : end if
2081 :
2082 : ! add secondary ice production due to accretion of droplets by snow
2083 : ! (Hallet-Mossop process) (from Cotton et al., 1986)
2084 :
2085 0 : if (.not. do_cldice) then
2086 0 : ni_secp = 0.0_r8
2087 0 : nsacwi(k) = 0.0_r8
2088 0 : msacwi(k) = 0.0_r8
2089 0 : else if((t(i,k).lt.270.16_r8) .and. (t(i,k).ge.268.16_r8)) then
2090 0 : ni_secp = 3.5e8_r8*(270.16_r8-t(i,k))/2.0_r8*psacws(k)
2091 0 : nsacwi(k) = ni_secp
2092 0 : msacwi(k) = min(ni_secp*mi0,psacws(k))
2093 0 : else if((t(i,k).lt.268.16_r8) .and. (t(i,k).ge.265.16_r8)) then
2094 0 : ni_secp = 3.5e8_r8*(t(i,k)-265.16_r8)/3.0_r8*psacws(k)
2095 0 : nsacwi(k) = ni_secp
2096 0 : msacwi(k) = min(ni_secp*mi0,psacws(k))
2097 : else
2098 0 : ni_secp = 0.0_r8
2099 0 : nsacwi(k) = 0.0_r8
2100 0 : msacwi(k) = 0.0_r8
2101 : endif
2102 0 : psacws(k) = max(0.0_r8,psacws(k)-ni_secp*mi0)
2103 :
2104 : !.......................................................................
2105 : ! accretion of rain water by snow
2106 : ! formula from ikawa and saito, 1991, used by reisner et al., 1998
2107 :
2108 0 : if (qric(i,k).ge.1.e-8_r8 .and. qniic(i,k).ge.1.e-8_r8 .and. &
2109 : t(i,k).le.273.15_r8) then
2110 :
2111 : pracs(k) = pi*pi*ecr*(((1.2_r8*umr(k)-0.95_r8*ums(k))**2+ &
2112 : 0.08_r8*ums(k)*umr(k))**0.5_r8*rhow*rho(i,k)* &
2113 : n0r(k)*n0s(k)* &
2114 : (5._r8/(lamr(k)**6*lams(k))+ &
2115 : 2._r8/(lamr(k)**5*lams(k)**2)+ &
2116 0 : 0.5_r8/(lamr(k)**4*lams(k)**3)))
2117 :
2118 : npracs(k) = pi/2._r8*rho(i,k)*ecr*(1.7_r8*(unr(k)-uns(k))**2+ &
2119 : 0.3_r8*unr(k)*uns(k))**0.5_r8*n0r(k)*n0s(k)* &
2120 : (1._r8/(lamr(k)**3*lams(k))+ &
2121 : 1._r8/(lamr(k)**2*lams(k)**2)+ &
2122 0 : 1._r8/(lamr(k)*lams(k)**3))
2123 :
2124 : else
2125 0 : pracs(k)=0._r8
2126 0 : npracs(k)=0._r8
2127 : end if
2128 :
2129 : !.......................................................................
2130 : ! heterogeneous freezing of rain drops
2131 : ! follows from Bigg (1953)
2132 :
2133 0 : if (t(i,k).lt.269.15_r8 .and. qric(i,k).ge.qsmall) then
2134 :
2135 : mnuccr(k) = 20._r8*pi*pi*rhow*nric(i,k)*bimm* &
2136 : (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamr(k)**3 &
2137 0 : /lamr(k)**3
2138 :
2139 : nnuccr(k) = pi*nric(i,k)*bimm* &
2140 0 : (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamr(k)**3
2141 : else
2142 0 : mnuccr(k)=0._r8
2143 0 : nnuccr(k)=0._r8
2144 : end if
2145 :
2146 : !.......................................................................
2147 : ! accretion of cloud liquid water by rain
2148 : ! formula from Khrouditnov and Kogan (2000)
2149 : ! gravitational collection kernel, droplet fall speed neglected
2150 :
2151 0 : if (qric(i,k).ge.qsmall .and. qcic(i,k).ge.qsmall) then
2152 :
2153 : ! include sub-grid distribution of cloud water
2154 :
2155 : ! add sub-column switch
2156 :
2157 0 : if (microp_uniform) then
2158 :
2159 0 : pra(k) = 67._r8*(qcic(i,k)*qric(i,k))**1.15_r8
2160 0 : npra(k) = pra(k)/(qcic(i,k)/ncic(i,k))
2161 :
2162 : else
2163 :
2164 0 : pra(k) = accre_enhan(i,k)*(cons12/(cons3*cons20)*67._r8*(qcic(i,k)*qric(i,k))**1.15_r8)
2165 0 : npra(k) = pra(k)/(qcic(i,k)/ncic(i,k))
2166 :
2167 : end if ! sub-column switch
2168 :
2169 : else
2170 0 : pra(k)=0._r8
2171 0 : npra(k)=0._r8
2172 : end if
2173 :
2174 : !.......................................................................
2175 : ! Self-collection of rain drops
2176 : ! from Beheng(1994)
2177 :
2178 0 : if (qric(i,k).ge.qsmall) then
2179 0 : nragg(k) = -8._r8*nric(i,k)*qric(i,k)*rho(i,k)
2180 : else
2181 0 : nragg(k)=0._r8
2182 : end if
2183 :
2184 : !.......................................................................
2185 : ! Accretion of cloud ice by snow
2186 : ! For this calculation, it is assumed that the Vs >> Vi
2187 : ! and Ds >> Di for continuous collection
2188 :
2189 : if (do_cldice .and. qniic(i,k).ge.qsmall.and.qiic(i,k).ge.qsmall &
2190 0 : .and.t(i,k).le.273.15_r8) then
2191 :
2192 : prai(k) = pi/4._r8*asn(i,k)*qiic(i,k)*rho(i,k)* &
2193 : n0s(k)*Eii*cons11/ &
2194 0 : lams(k)**(bs+3._r8)
2195 : nprai(k) = pi/4._r8*asn(i,k)*niic(i,k)* &
2196 : rho(i,k)*n0s(k)*Eii*cons11/ &
2197 0 : lams(k)**(bs+3._r8)
2198 : else
2199 0 : prai(k)=0._r8
2200 0 : nprai(k)=0._r8
2201 : end if
2202 :
2203 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2204 : ! calculate evaporation/sublimation of rain and snow
2205 : ! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
2206 : ! in-cloud condensation/deposition of rain and snow is neglected
2207 : ! except for transfer of cloud water to snow through bergeron process
2208 :
2209 : ! initialize evap/sub tendncies
2210 0 : pre(k)=0._r8
2211 0 : prds(k)=0._r8
2212 :
2213 : ! evaporation of rain
2214 : ! only calculate if there is some precip fraction > cloud fraction
2215 :
2216 0 : if (qcic(i,k)+qiic(i,k).lt.1.e-6_r8.or.cldmax(i,k).gt.lcldm(i,k)) then
2217 :
2218 : ! set temporary cloud fraction to zero if cloud water + ice is very small
2219 : ! this will ensure that evaporation/sublimation of precip occurs over
2220 : ! entire grid cell, since min cloud fraction is specified otherwise
2221 0 : if (qcic(i,k)+qiic(i,k).lt.1.e-6_r8) then
2222 : dum=0._r8
2223 : else
2224 0 : dum=lcldm(i,k)
2225 : end if
2226 :
2227 : ! saturation vapor pressure
2228 0 : esn=svp_water(t(i,k))
2229 0 : qsn=svp_to_qsat(esn, p(i,k))
2230 :
2231 : ! recalculate saturation vapor pressure for liquid and ice
2232 0 : esl(i,k)=esn
2233 0 : esi(i,k)=svp_ice(t(i,k))
2234 : ! hm fix, make sure when above freezing that esi=esl, not active yet
2235 0 : if (t(i,k).gt.tmelt)esi(i,k)=esl(i,k)
2236 :
2237 : ! calculate q for out-of-cloud region
2238 0 : qclr=(q(i,k)-dum*qsn)/(1._r8-dum)
2239 :
2240 0 : if (qric(i,k).ge.qsmall) then
2241 :
2242 0 : qvs=svp_to_qsat(esl(i,k), p(i,k))
2243 0 : dqsdt = xxlv*qvs/(rv*t(i,k)**2)
2244 0 : ab = 1._r8+dqsdt*xxlv/cpp
2245 : epsr = 2._r8*pi*n0r(k)*rho(i,k)*Dv(i,k)* &
2246 : (f1r/(lamr(k)*lamr(k))+ &
2247 : f2r*(arn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
2248 : sc(i,k)**(1._r8/3._r8)*cons13/ &
2249 0 : (lamr(k)**(5._r8/2._r8+br/2._r8)))
2250 :
2251 0 : pre(k) = epsr*(qclr-qvs)/ab
2252 :
2253 : ! only evaporate in out-of-cloud region
2254 : ! and distribute across cldmax
2255 0 : pre(k)=min(pre(k)*(cldmax(i,k)-dum),0._r8)
2256 0 : pre(k)=pre(k)/cldmax(i,k)
2257 0 : am_evp_st(i,k) = max(cldmax(i,k)-dum, 0._r8)
2258 : end if
2259 :
2260 : ! sublimation of snow
2261 0 : if (qniic(i,k).ge.qsmall) then
2262 0 : qvi=svp_to_qsat(esi(i,k), p(i,k))
2263 0 : dqsidt = xxls*qvi/(rv*t(i,k)**2)
2264 0 : abi = 1._r8+dqsidt*xxls/cpp
2265 : epss = 2._r8*pi*n0s(k)*rho(i,k)*Dv(i,k)* &
2266 : (f1s/(lams(k)*lams(k))+ &
2267 : f2s*(asn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
2268 : sc(i,k)**(1._r8/3._r8)*cons14/ &
2269 0 : (lams(k)**(5._r8/2._r8+bs/2._r8)))
2270 0 : prds(k) = epss*(qclr-qvi)/abi
2271 :
2272 : ! only sublimate in out-of-cloud region and distribute over cldmax
2273 0 : prds(k)=min(prds(k)*(cldmax(i,k)-dum),0._r8)
2274 0 : prds(k)=prds(k)/cldmax(i,k)
2275 0 : am_evp_st(i,k) = max(cldmax(i,k)-dum, 0._r8)
2276 : end if
2277 :
2278 : ! make sure RH not pushed above 100% due to rain evaporation/snow sublimation
2279 : ! get updated RH at end of time step based on cloud water/ice condensation/evap
2280 :
2281 0 : qtmp=q(i,k)-(cmei(i,k)+(pre(k)+prds(k))*cldmax(i,k))*deltat
2282 : ttmp=t(i,k)+((pre(k)*cldmax(i,k))*xxlv+ &
2283 0 : (cmei(i,k)+prds(k)*cldmax(i,k))*xxls)*deltat/cpp
2284 :
2285 : !limit range of temperatures!
2286 0 : ttmp=max(180._r8,min(ttmp,323._r8))
2287 :
2288 0 : esn=svp_water(ttmp) ! use rhw to allow ice supersaturation
2289 0 : qsn=svp_to_qsat(esn, p(i,k))
2290 :
2291 : ! modify precip evaporation rate if q > qsat
2292 0 : if (qtmp.gt.qsn) then
2293 0 : if (pre(k)+prds(k).lt.-1.e-20_r8) then
2294 0 : dum1=pre(k)/(pre(k)+prds(k))
2295 : ! recalculate q and t after cloud water cond but without precip evap
2296 0 : qtmp=q(i,k)-(cmei(i,k))*deltat
2297 0 : ttmp=t(i,k)+(cmei(i,k)*xxls)*deltat/cpp
2298 0 : esn=svp_water(ttmp) ! use rhw to allow ice supersaturation
2299 0 : qsn=svp_to_qsat(esn, p(i,k))
2300 0 : dum=(qtmp-qsn)/(1._r8 + cons27*qsn/(cpp*rv*ttmp**2))
2301 0 : dum=min(dum,0._r8)
2302 :
2303 : ! modify rates if needed, divide by cldmax to get local (in-precip) value
2304 0 : pre(k)=dum*dum1/deltat/cldmax(i,k)
2305 :
2306 : ! do separately using RHI for prds....
2307 0 : esn=svp_ice(ttmp) ! use rhi to allow ice supersaturation
2308 0 : qsn=svp_to_qsat(esn, p(i,k))
2309 0 : dum=(qtmp-qsn)/(1._r8 + cons28*qsn/(cpp*rv*ttmp**2))
2310 0 : dum=min(dum,0._r8)
2311 :
2312 : ! modify rates if needed, divide by cldmax to get local (in-precip) value
2313 0 : prds(k)=dum*(1._r8-dum1)/deltat/cldmax(i,k)
2314 : end if
2315 : end if
2316 : end if
2317 :
2318 : ! bergeron process - evaporation of droplets and deposition onto snow
2319 :
2320 0 : if (qniic(i,k).ge.qsmall.and.qcic(i,k).ge.qsmall.and.t(i,k).lt.tmelt) then
2321 0 : qvi=svp_to_qsat(esi(i,k), p(i,k))
2322 0 : qvs=svp_to_qsat(esl(i,k), p(i,k))
2323 0 : dqsidt = xxls*qvi/(rv*t(i,k)**2)
2324 0 : abi = 1._r8+dqsidt*xxls/cpp
2325 : epss = 2._r8*pi*n0s(k)*rho(i,k)*Dv(i,k)* &
2326 : (f1s/(lams(k)*lams(k))+ &
2327 : f2s*(asn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
2328 : sc(i,k)**(1._r8/3._r8)*cons14/ &
2329 0 : (lams(k)**(5._r8/2._r8+bs/2._r8)))
2330 0 : bergs(k)=epss*(qvs-qvi)/abi
2331 : else
2332 0 : bergs(k)=0._r8
2333 : end if
2334 :
2335 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2336 : ! conservation to ensure no negative values of cloud water/precipitation
2337 : ! in case microphysical process rates are large
2338 :
2339 : ! make sure and use end-of-time step values for cloud water, ice, due
2340 : ! condensation/deposition
2341 :
2342 : ! note: for check on conservation, processes are multiplied by omsm
2343 : ! to prevent problems due to round off error
2344 :
2345 : ! include mixing timescale (mtime)
2346 :
2347 0 : qce=(qc(i,k) - berg(i,k)*deltat)
2348 0 : nce=(nc(i,k)+npccn(k)*deltat*mtime)
2349 0 : qie=(qi(i,k)+(cmei(i,k)+berg(i,k))*deltat)
2350 0 : nie=(ni(i,k)+nnuccd(k)*deltat*mtime)
2351 :
2352 : ! conservation of qc
2353 :
2354 : dum = (prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+ &
2355 0 : psacws(k)+bergs(k))*lcldm(i,k)*deltat
2356 :
2357 0 : if (dum.gt.qce) then
2358 0 : ratio = qce/deltat/lcldm(i,k)/(prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+psacws(k)+bergs(k))*omsm
2359 :
2360 0 : prc(k) = prc(k)*ratio
2361 0 : pra(k) = pra(k)*ratio
2362 0 : mnuccc(k) = mnuccc(k)*ratio
2363 0 : mnucct(k) = mnucct(k)*ratio
2364 0 : msacwi(k) = msacwi(k)*ratio
2365 0 : psacws(k) = psacws(k)*ratio
2366 0 : bergs(k) = bergs(k)*ratio
2367 : end if
2368 :
2369 : ! conservation of nc
2370 :
2371 : dum = (nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+ &
2372 0 : npsacws(k)-nsubc(k))*lcldm(i,k)*deltat
2373 :
2374 0 : if (dum.gt.nce) then
2375 : ratio = nce/deltat/((nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+&
2376 0 : npsacws(k)-nsubc(k))*lcldm(i,k))*omsm
2377 :
2378 0 : nprc1(k) = nprc1(k)*ratio
2379 0 : npra(k) = npra(k)*ratio
2380 0 : nnuccc(k) = nnuccc(k)*ratio
2381 0 : nnucct(k) = nnucct(k)*ratio
2382 0 : npsacws(k) = npsacws(k)*ratio
2383 0 : nsubc(k)=nsubc(k)*ratio
2384 : end if
2385 :
2386 : ! conservation of qi
2387 :
2388 0 : if (do_cldice) then
2389 :
2390 0 : frztmp = -mnuccc(k) - mnucct(k) - msacwi(k)
2391 0 : if (use_hetfrz_classnuc) frztmp = -mnuccc(k)-mnucct(k)-mnudep(k)-msacwi(k)
2392 0 : dum = ( frztmp*lcldm(i,k) + (prci(k)+prai(k))*icldm(i,k) )*deltat
2393 :
2394 0 : if (dum.gt.qie) then
2395 :
2396 0 : frztmp = mnuccc(k) + mnucct(k) + msacwi(k)
2397 0 : if (use_hetfrz_classnuc) frztmp = mnuccc(k) + mnucct(k) + mnudep(k) + msacwi(k)
2398 0 : ratio = (qie/deltat + frztmp*lcldm(i,k))/((prci(k)+prai(k))*icldm(i,k))*omsm
2399 0 : prci(k) = prci(k)*ratio
2400 0 : prai(k) = prai(k)*ratio
2401 : end if
2402 :
2403 : ! conservation of ni
2404 0 : frztmp = -nnucct(k) - nsacwi(k)
2405 0 : if (use_hetfrz_classnuc) frztmp = -nnucct(k) - nnuccc(k) - nnudep(k) - nsacwi(k)
2406 0 : dum = ( frztmp*lcldm(i,k) + (nprci(k)+nprai(k)-nsubi(k))*icldm(i,k) )*deltat
2407 :
2408 0 : if (dum.gt.nie) then
2409 :
2410 0 : frztmp = nnucct(k) + nsacwi(k)
2411 0 : if (use_hetfrz_classnuc) frztmp = nnucct(k) + nnuccc(k) + nnudep(k) + nsacwi(k)
2412 : ratio = (nie/deltat + frztmp*lcldm(i,k))/ &
2413 0 : ((nprci(k)+nprai(k)-nsubi(k))*icldm(i,k))*omsm
2414 0 : nprci(k) = nprci(k)*ratio
2415 0 : nprai(k) = nprai(k)*ratio
2416 0 : nsubi(k) = nsubi(k)*ratio
2417 : end if
2418 : end if
2419 :
2420 : ! for precipitation conservation, use logic that vertical integral
2421 : ! of tendency from current level to top of model (i.e., qrtot) cannot be negative
2422 :
2423 : ! conservation of rain mixing rat
2424 :
2425 0 : if (((prc(k)+pra(k))*lcldm(i,k)+(-mnuccr(k)+pre(k)-pracs(k))*&
2426 : cldmax(i,k))*dz(i,k)*rho(i,k)+qrtot.lt.0._r8) then
2427 :
2428 0 : if (-pre(k)+pracs(k)+mnuccr(k).ge.qsmall) then
2429 :
2430 : ratio = (qrtot/(dz(i,k)*rho(i,k))+(prc(k)+pra(k))*lcldm(i,k))/&
2431 0 : ((-pre(k)+pracs(k)+mnuccr(k))*cldmax(i,k))*omsm
2432 :
2433 0 : pre(k) = pre(k)*ratio
2434 0 : pracs(k) = pracs(k)*ratio
2435 0 : mnuccr(k) = mnuccr(k)*ratio
2436 : end if
2437 : end if
2438 :
2439 : ! conservation of nr
2440 : ! for now neglect evaporation of nr
2441 0 : nsubr(k)=0._r8
2442 :
2443 0 : if ((nprc(k)*lcldm(i,k)+(-nnuccr(k)+nsubr(k)-npracs(k)&
2444 : +nragg(k))*cldmax(i,k))*dz(i,k)*rho(i,k)+nrtot.lt.0._r8) then
2445 :
2446 0 : if (-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k).ge.qsmall) then
2447 :
2448 : ratio = (nrtot/(dz(i,k)*rho(i,k))+nprc(k)*lcldm(i,k))/&
2449 0 : ((-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k))*cldmax(i,k))*omsm
2450 :
2451 0 : nsubr(k) = nsubr(k)*ratio
2452 0 : npracs(k) = npracs(k)*ratio
2453 0 : nnuccr(k) = nnuccr(k)*ratio
2454 0 : nragg(k) = nragg(k)*ratio
2455 : end if
2456 : end if
2457 :
2458 : ! conservation of snow mix ratio
2459 :
2460 0 : if (((bergs(k)+psacws(k))*lcldm(i,k)+(prai(k)+prci(k))*icldm(i,k)+(pracs(k)+&
2461 : mnuccr(k)+prds(k))*cldmax(i,k))*dz(i,k)*rho(i,k)+qstot.lt.0._r8) then
2462 :
2463 0 : if (-prds(k).ge.qsmall) then
2464 :
2465 : ratio = (qstot/(dz(i,k)*rho(i,k))+(bergs(k)+psacws(k))*lcldm(i,k)+(prai(k)+prci(k))*icldm(i,k)+&
2466 0 : (pracs(k)+mnuccr(k))*cldmax(i,k))/(-prds(k)*cldmax(i,k))*omsm
2467 :
2468 0 : prds(k) = prds(k)*ratio
2469 : end if
2470 : end if
2471 :
2472 : ! conservation of ns
2473 :
2474 : ! calculate loss of number due to sublimation
2475 : ! for now neglect sublimation of ns
2476 0 : nsubs(k)=0._r8
2477 :
2478 0 : if ((nprci(k)*icldm(i,k)+(nnuccr(k)+nsubs(k)+nsagg(k))*cldmax(i,k))*&
2479 : dz(i,k)*rho(i,k)+nstot.lt.0._r8) then
2480 :
2481 0 : if (-nsubs(k)-nsagg(k).ge.qsmall) then
2482 :
2483 : ratio = (nstot/(dz(i,k)*rho(i,k))+nprci(k)*icldm(i,k)+&
2484 0 : nnuccr(k)*cldmax(i,k))/((-nsubs(k)-nsagg(k))*cldmax(i,k))*omsm
2485 :
2486 0 : nsubs(k) = nsubs(k)*ratio
2487 0 : nsagg(k) = nsagg(k)*ratio
2488 : end if
2489 : end if
2490 :
2491 : ! get tendencies due to microphysical conversion processes
2492 : ! note: tendencies are multiplied by appropaiate cloud/precip
2493 : ! fraction to get grid-scale values
2494 : ! note: cmei is already grid-average values
2495 :
2496 0 : qvlat(i,k) = qvlat(i,k)-(pre(k)+prds(k))*cldmax(i,k)-cmei(i,k)
2497 :
2498 : tlat(i,k) = tlat(i,k)+((pre(k)*cldmax(i,k)) &
2499 : *xxlv+(prds(k)*cldmax(i,k)+cmei(i,k))*xxls+ &
2500 : ((bergs(k)+psacws(k)+mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k)+(mnuccr(k)+ &
2501 0 : pracs(k))*cldmax(i,k)+berg(i,k))*xlf)
2502 :
2503 : qctend(i,k) = qctend(i,k)+ &
2504 : (-pra(k)-prc(k)-mnuccc(k)-mnucct(k)-msacwi(k)- &
2505 0 : psacws(k)-bergs(k))*lcldm(i,k)-berg(i,k)
2506 :
2507 0 : if (do_cldice) then
2508 :
2509 0 : frztmp = mnuccc(k) + mnucct(k) + msacwi(k)
2510 0 : if (use_hetfrz_classnuc) frztmp = mnuccc(k) + mnucct(k) + mnudep(k) + msacwi(k)
2511 : qitend(i,k) = qitend(i,k) + frztmp*lcldm(i,k) + &
2512 0 : (-prci(k)-prai(k))*icldm(i,k) + cmei(i,k) + berg(i,k)
2513 :
2514 : end if
2515 :
2516 : qrtend(i,k) = qrtend(i,k)+ &
2517 : (pra(k)+prc(k))*lcldm(i,k)+(pre(k)-pracs(k)- &
2518 0 : mnuccr(k))*cldmax(i,k)
2519 :
2520 : qnitend(i,k) = qnitend(i,k)+ &
2521 : (prai(k)+prci(k))*icldm(i,k)+(psacws(k)+bergs(k))*lcldm(i,k)+(prds(k)+ &
2522 0 : pracs(k)+mnuccr(k))*cldmax(i,k)
2523 :
2524 : ! add output for cmei (accumulate)
2525 0 : cmeiout(i,k) = cmeiout(i,k) + cmei(i,k)
2526 :
2527 : ! assign variables for trop_mozart, these are grid-average
2528 : ! evaporation/sublimation is stored here as positive term
2529 :
2530 0 : evapsnow(i,k) = evapsnow(i,k)-prds(k)*cldmax(i,k)
2531 0 : nevapr(i,k) = nevapr(i,k)-pre(k)*cldmax(i,k)
2532 0 : nevapr2(i,k) = nevapr2(i,k)-pre(k)*cldmax(i,k)
2533 :
2534 : ! change to make sure prain is positive: do not remove snow from
2535 : ! prain used for wet deposition
2536 : prain(i,k) = prain(i,k)+(pra(k)+prc(k))*lcldm(i,k)+(-pracs(k)- &
2537 0 : mnuccr(k))*cldmax(i,k)
2538 : prodsnow(i,k) = prodsnow(i,k)+(prai(k)+prci(k))*icldm(i,k)+(psacws(k)+bergs(k))*lcldm(i,k)+(&
2539 0 : pracs(k)+mnuccr(k))*cldmax(i,k)
2540 :
2541 : ! following are used to calculate 1st order conversion rate of cloud water
2542 : ! to rain and snow (1/s), for later use in aerosol wet removal routine
2543 : ! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc
2544 : ! used to calculate pra, prc, ... in this routine
2545 : ! qcsinksum_rate1ord = sum over iterations{ rate of direct transfer of cloud water to rain & snow }
2546 : ! (no cloud ice or bergeron terms)
2547 : ! qcsum_rate1ord = sum over iterations{ qc used in calculation of the transfer terms }
2548 :
2549 0 : qcsinksum_rate1ord(k) = qcsinksum_rate1ord(k) + (pra(k)+prc(k)+psacws(k))*lcldm(i,k)
2550 0 : qcsum_rate1ord(k) = qcsum_rate1ord(k) + qc(i,k)
2551 :
2552 : ! microphysics output, note this is grid-averaged
2553 0 : prao(i,k)=prao(i,k)+pra(k)*lcldm(i,k)
2554 0 : prco(i,k)=prco(i,k)+prc(k)*lcldm(i,k)
2555 0 : mnuccco(i,k)=mnuccco(i,k)+mnuccc(k)*lcldm(i,k)
2556 0 : mnuccto(i,k)=mnuccto(i,k)+mnucct(k)*lcldm(i,k)
2557 0 : mnuccdo(i,k)=mnuccdo(i,k)+mnuccd(k)*lcldm(i,k)
2558 0 : msacwio(i,k)=msacwio(i,k)+msacwi(k)*lcldm(i,k)
2559 0 : psacwso(i,k)=psacwso(i,k)+psacws(k)*lcldm(i,k)
2560 0 : bergso(i,k)=bergso(i,k)+bergs(k)*lcldm(i,k)
2561 0 : bergo(i,k)=bergo(i,k)+berg(i,k)
2562 0 : prcio(i,k)=prcio(i,k)+prci(k)*icldm(i,k)
2563 0 : praio(i,k)=praio(i,k)+prai(k)*icldm(i,k)
2564 0 : mnuccro(i,k)=mnuccro(i,k)+mnuccr(k)*cldmax(i,k)
2565 0 : pracso (i,k)=pracso (i,k)+pracs (k)*cldmax(i,k)
2566 :
2567 : ! multiply activation/nucleation by mtime to account for fast timescale
2568 :
2569 : nctend(i,k) = nctend(i,k)+ npccn(k)*mtime+&
2570 : (-nnuccc(k)-nnucct(k)-npsacws(k)+nsubc(k) &
2571 0 : -npra(k)-nprc1(k))*lcldm(i,k)
2572 :
2573 0 : if (do_cldice) then
2574 :
2575 0 : frztmp = nnucct(k) + nsacwi(k)
2576 0 : if (use_hetfrz_classnuc) frztmp = nnucct(k) + nnuccc(k) + nnudep(k) + nsacwi(k)
2577 : nitend(i,k) = nitend(i,k) + nnuccd(k)*mtime + &
2578 0 : frztmp*lcldm(i,k) + (nsubi(k)-nprci(k)-nprai(k))*icldm(i,k)
2579 :
2580 : end if
2581 :
2582 : nstend(i,k) = nstend(i,k)+(nsubs(k)+ &
2583 0 : nsagg(k)+nnuccr(k))*cldmax(i,k)+nprci(k)*icldm(i,k)
2584 :
2585 : nrtend(i,k) = nrtend(i,k)+ &
2586 : nprc(k)*lcldm(i,k)+(nsubr(k)-npracs(k)-nnuccr(k) &
2587 0 : +nragg(k))*cldmax(i,k)
2588 :
2589 : ! make sure that nc and ni at advanced time step do not exceed
2590 : ! maximum (existing N + source terms*dt), which is possible due to
2591 : ! fast nucleation timescale
2592 :
2593 0 : if (nctend(i,k).gt.0._r8.and.nc(i,k)+nctend(i,k)*deltat.gt.ncmax) then
2594 0 : nctend(i,k)=max(0._r8,(ncmax-nc(i,k))/deltat)
2595 : end if
2596 :
2597 0 : if (do_cldice .and. nitend(i,k).gt.0._r8.and.ni(i,k)+nitend(i,k)*deltat.gt.nimax) then
2598 0 : nitend(i,k)=max(0._r8,(nimax-ni(i,k))/deltat)
2599 : end if
2600 :
2601 : ! get final values for precipitation q and N, based on
2602 : ! flux of precip from above, source/sink term, and terminal fallspeed
2603 : ! see eq. 15-16 in MG2008
2604 :
2605 : ! rain
2606 :
2607 0 : if (qric(i,k).ge.qsmall) then
2608 0 : if (k.eq.top_lev) then
2609 0 : qric(i,k)=qrtend(i,k)*dz(i,k)/cldmax(i,k)/umr(k)
2610 0 : nric(i,k)=nrtend(i,k)*dz(i,k)/cldmax(i,k)/unr(k)
2611 : else
2612 0 : qric(i,k) = (rho(i,k-1)*umr(k-1)*qric(i,k-1)*cldmax(i,k-1)+ &
2613 0 : (rho(i,k)*dz(i,k)*qrtend(i,k)))/(umr(k)*rho(i,k)*cldmax(i,k))
2614 : nric(i,k) = (rho(i,k-1)*unr(k-1)*nric(i,k-1)*cldmax(i,k-1)+ &
2615 0 : (rho(i,k)*dz(i,k)*nrtend(i,k)))/(unr(k)*rho(i,k)*cldmax(i,k))
2616 :
2617 : end if
2618 : else
2619 0 : qric(i,k)=0._r8
2620 0 : nric(i,k)=0._r8
2621 : end if
2622 :
2623 : ! snow
2624 :
2625 0 : if (qniic(i,k).ge.qsmall) then
2626 0 : if (k.eq.top_lev) then
2627 0 : qniic(i,k)=qnitend(i,k)*dz(i,k)/cldmax(i,k)/ums(k)
2628 0 : nsic(i,k)=nstend(i,k)*dz(i,k)/cldmax(i,k)/uns(k)
2629 : else
2630 0 : qniic(i,k) = (rho(i,k-1)*ums(k-1)*qniic(i,k-1)*cldmax(i,k-1)+ &
2631 0 : (rho(i,k)*dz(i,k)*qnitend(i,k)))/(ums(k)*rho(i,k)*cldmax(i,k))
2632 : nsic(i,k) = (rho(i,k-1)*uns(k-1)*nsic(i,k-1)*cldmax(i,k-1)+ &
2633 0 : (rho(i,k)*dz(i,k)*nstend(i,k)))/(uns(k)*rho(i,k)*cldmax(i,k))
2634 : end if
2635 : else
2636 0 : qniic(i,k)=0._r8
2637 0 : nsic(i,k)=0._r8
2638 : end if
2639 :
2640 : ! calculate precipitation flux at surface
2641 : ! divide by density of water to get units of m/s
2642 :
2643 : prect(i) = prect(i)+(qrtend(i,k)*dz(i,k)*rho(i,k)+&
2644 0 : qnitend(i,k)*dz(i,k)*rho(i,k))/rhow
2645 0 : preci(i) = preci(i)+qnitend(i,k)*dz(i,k)*rho(i,k)/rhow
2646 :
2647 : ! convert rain rate from m/s to mm/hr
2648 :
2649 0 : rainrt(i,k)=qric(i,k)*rho(i,k)*umr(k)/rhow*3600._r8*1000._r8
2650 :
2651 : ! vertically-integrated precip source/sink terms (note: grid-averaged)
2652 :
2653 0 : qrtot = max(qrtot+qrtend(i,k)*dz(i,k)*rho(i,k),0._r8)
2654 0 : qstot = max(qstot+qnitend(i,k)*dz(i,k)*rho(i,k),0._r8)
2655 0 : nrtot = max(nrtot+nrtend(i,k)*dz(i,k)*rho(i,k),0._r8)
2656 0 : nstot = max(nstot+nstend(i,k)*dz(i,k)*rho(i,k),0._r8)
2657 :
2658 : ! calculate melting and freezing of precip
2659 :
2660 : ! melt snow at +2 C
2661 :
2662 0 : if (t(i,k)+tlat(i,k)/cpp*deltat > 275.15_r8) then
2663 0 : if (qstot > 0._r8) then
2664 :
2665 : ! make sure melting snow doesn't reduce temperature below threshold
2666 0 : dum = -xlf/cpp*qstot/(dz(i,k)*rho(i,k))
2667 0 : if (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt.275.15_r8) then
2668 0 : dum = (t(i,k)+tlat(i,k)/cpp*deltat-275.15_r8)*cpp/xlf
2669 0 : dum = dum/(xlf/cpp*qstot/(dz(i,k)*rho(i,k)))
2670 0 : dum = max(0._r8,dum)
2671 0 : dum = min(1._r8,dum)
2672 : else
2673 : dum = 1._r8
2674 : end if
2675 :
2676 0 : qric(i,k)=qric(i,k)+dum*qniic(i,k)
2677 0 : nric(i,k)=nric(i,k)+dum*nsic(i,k)
2678 0 : qniic(i,k)=(1._r8-dum)*qniic(i,k)
2679 0 : nsic(i,k)=(1._r8-dum)*nsic(i,k)
2680 : ! heating tendency
2681 0 : tmp=-xlf*dum*qstot/(dz(i,k)*rho(i,k))
2682 0 : meltsdt(i,k)=meltsdt(i,k) + tmp
2683 :
2684 0 : tlat(i,k)=tlat(i,k)+tmp
2685 0 : qrtot=qrtot+dum*qstot
2686 0 : nrtot=nrtot+dum*nstot
2687 0 : qstot=(1._r8-dum)*qstot
2688 0 : nstot=(1._r8-dum)*nstot
2689 0 : preci(i)=(1._r8-dum)*preci(i)
2690 : end if
2691 : end if
2692 :
2693 : ! freeze all rain at -5C for Arctic
2694 :
2695 0 : if (t(i,k)+tlat(i,k)/cpp*deltat < (tmelt - 5._r8)) then
2696 :
2697 0 : if (qrtot > 0._r8) then
2698 :
2699 : ! make sure freezing rain doesn't increase temperature above threshold
2700 0 : dum = xlf/cpp*qrtot/(dz(i,k)*rho(i,k))
2701 0 : if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.(tmelt - 5._r8)) then
2702 0 : dum = -(t(i,k)+tlat(i,k)/cpp*deltat-(tmelt-5._r8))*cpp/xlf
2703 0 : dum = dum/(xlf/cpp*qrtot/(dz(i,k)*rho(i,k)))
2704 0 : dum = max(0._r8,dum)
2705 0 : dum = min(1._r8,dum)
2706 : else
2707 : dum = 1._r8
2708 : end if
2709 :
2710 0 : qniic(i,k)=qniic(i,k)+dum*qric(i,k)
2711 0 : nsic(i,k)=nsic(i,k)+dum*nric(i,k)
2712 0 : qric(i,k)=(1._r8-dum)*qric(i,k)
2713 0 : nric(i,k)=(1._r8-dum)*nric(i,k)
2714 : ! heating tendency
2715 0 : tmp = xlf*dum*qrtot/(dz(i,k)*rho(i,k))
2716 0 : frzrdt(i,k)=frzrdt(i,k) + tmp
2717 :
2718 0 : tlat(i,k)=tlat(i,k)+tmp
2719 0 : qstot=qstot+dum*qrtot
2720 0 : qrtot=(1._r8-dum)*qrtot
2721 0 : nstot=nstot+dum*nrtot
2722 0 : nrtot=(1._r8-dum)*nrtot
2723 0 : preci(i)=preci(i)+dum*(prect(i)-preci(i))
2724 : end if
2725 : end if
2726 :
2727 : ! Precip Flux Calculation (Diagnostic)
2728 0 : rflx(i,k+1)=(prect(i)-preci(i)) * rhow
2729 0 : sflx(i,k+1)=preci(i) * rhow
2730 :
2731 : ! if rain/snow mix ratio is zero so should number concentration
2732 :
2733 0 : if (qniic(i,k).lt.qsmall) then
2734 0 : qniic(i,k)=0._r8
2735 0 : nsic(i,k)=0._r8
2736 : end if
2737 :
2738 0 : if (qric(i,k).lt.qsmall) then
2739 0 : qric(i,k)=0._r8
2740 0 : nric(i,k)=0._r8
2741 : end if
2742 :
2743 : ! make sure number concentration is a positive number to avoid
2744 : ! taking root of negative
2745 :
2746 0 : nric(i,k)=max(nric(i,k),0._r8)
2747 0 : nsic(i,k)=max(nsic(i,k),0._r8)
2748 :
2749 : !.......................................................................
2750 : ! get size distribution parameters for fallspeed calculations
2751 : !......................................................................
2752 : ! rain
2753 :
2754 0 : if (qric(i,k).ge.qsmall) then
2755 0 : lamr(k) = (pi*rhow*nric(i,k)/qric(i,k))**(1._r8/3._r8)
2756 0 : n0r(k) = nric(i,k)*lamr(k)
2757 :
2758 : ! check for slope
2759 : ! change lammax and lammin for rain and snow
2760 : ! adjust vars
2761 :
2762 0 : if (lamr(k).lt.lamminr) then
2763 :
2764 0 : lamr(k) = lamminr
2765 :
2766 0 : n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
2767 0 : nric(i,k) = n0r(k)/lamr(k)
2768 0 : else if (lamr(k).gt.lammaxr) then
2769 0 : lamr(k) = lammaxr
2770 0 : n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
2771 0 : nric(i,k) = n0r(k)/lamr(k)
2772 : end if
2773 :
2774 :
2775 : ! 'final' values of number and mass weighted mean fallspeed for rain (m/s)
2776 :
2777 0 : unr(k) = min(arn(i,k)*cons4/lamr(k)**br,9.1_r8*rhof(i,k))
2778 0 : umr(k) = min(arn(i,k)*cons5/(6._r8*lamr(k)**br),9.1_r8*rhof(i,k))
2779 :
2780 : else
2781 0 : lamr(k) = 0._r8
2782 0 : n0r(k) = 0._r8
2783 0 : umr(k)=0._r8
2784 0 : unr(k)=0._r8
2785 : end if
2786 :
2787 : !calculate mean size of combined rain and snow
2788 :
2789 0 : if (lamr(k).gt.0._r8) then
2790 0 : Artmp = n0r(k) * pi / (2._r8 * lamr(k)**3._r8)
2791 : else
2792 : Artmp = 0._r8
2793 : endif
2794 :
2795 0 : if (lamc(k).gt.0._r8) then
2796 0 : Actmp = cdist1(k) * pi * gamma(pgam(k)+3._r8)/(4._r8 * lamc(k)**2._r8)
2797 : else
2798 : Actmp = 0._r8
2799 : endif
2800 :
2801 0 : if (Actmp.gt.0_r8.or.Artmp.gt.0) then
2802 0 : rercld(i,k)=rercld(i,k) + 3._r8 *(qric(i,k) + qcic(i,k)) / (4._r8 * rhow * (Actmp + Artmp))
2803 0 : arcld(i,k)=arcld(i,k)+1._r8
2804 : endif
2805 :
2806 : !......................................................................
2807 : ! snow
2808 :
2809 0 : if (qniic(i,k).ge.qsmall) then
2810 : lams(k) = (cons6*cs*nsic(i,k)/ &
2811 0 : qniic(i,k))**(1._r8/ds)
2812 0 : n0s(k) = nsic(i,k)*lams(k)
2813 :
2814 : ! check for slope
2815 : ! adjust vars
2816 :
2817 0 : if (lams(k).lt.lammins) then
2818 0 : lams(k) = lammins
2819 0 : n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
2820 0 : nsic(i,k) = n0s(k)/lams(k)
2821 :
2822 0 : else if (lams(k).gt.lammaxs) then
2823 0 : lams(k) = lammaxs
2824 0 : n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
2825 0 : nsic(i,k) = n0s(k)/lams(k)
2826 : end if
2827 :
2828 : ! 'final' values of number and mass weighted mean fallspeed for snow (m/s)
2829 :
2830 0 : ums(k) = min(asn(i,k)*cons8/(6._r8*lams(k)**bs),1.2_r8*rhof(i,k))
2831 0 : uns(k) = min(asn(i,k)*cons7/lams(k)**bs,1.2_r8*rhof(i,k))
2832 :
2833 : else
2834 0 : lams(k) = 0._r8
2835 0 : n0s(k) = 0._r8
2836 0 : ums(k) = 0._r8
2837 0 : uns(k) = 0._r8
2838 : end if
2839 :
2840 : !c........................................................................
2841 : ! sum over sub-step for average process rates
2842 :
2843 : ! convert rain/snow q and N for output to history, note,
2844 : ! output is for gridbox average
2845 :
2846 0 : qrout(i,k)=qrout(i,k)+qric(i,k)*cldmax(i,k)
2847 0 : qsout(i,k)=qsout(i,k)+qniic(i,k)*cldmax(i,k)
2848 0 : nrout(i,k)=nrout(i,k)+nric(i,k)*rho(i,k)*cldmax(i,k)
2849 0 : nsout(i,k)=nsout(i,k)+nsic(i,k)*rho(i,k)*cldmax(i,k)
2850 :
2851 0 : tlat1(i,k)=tlat1(i,k)+tlat(i,k)
2852 0 : qvlat1(i,k)=qvlat1(i,k)+qvlat(i,k)
2853 0 : qctend1(i,k)=qctend1(i,k)+qctend(i,k)
2854 0 : qitend1(i,k)=qitend1(i,k)+qitend(i,k)
2855 0 : nctend1(i,k)=nctend1(i,k)+nctend(i,k)
2856 0 : nitend1(i,k)=nitend1(i,k)+nitend(i,k)
2857 :
2858 0 : t(i,k)=t(i,k)+tlat(i,k)*deltat/cpp
2859 0 : q(i,k)=q(i,k)+qvlat(i,k)*deltat
2860 0 : qc(i,k)=qc(i,k)+qctend(i,k)*deltat
2861 0 : qi(i,k)=qi(i,k)+qitend(i,k)*deltat
2862 0 : nc(i,k)=nc(i,k)+nctend(i,k)*deltat
2863 0 : ni(i,k)=ni(i,k)+nitend(i,k)*deltat
2864 :
2865 0 : rainrt1(i,k)=rainrt1(i,k)+rainrt(i,k)
2866 :
2867 : !divide rain radius over substeps for average
2868 0 : if (arcld(i,k) .gt. 0._r8) then
2869 0 : rercld(i,k)=rercld(i,k)/arcld(i,k)
2870 : end if
2871 :
2872 : !! add to summing sub-stepping variable
2873 0 : rflx1(i,k+1)=rflx1(i,k+1)+rflx(i,k+1)
2874 0 : sflx1(i,k+1)=sflx1(i,k+1)+sflx(i,k+1)
2875 :
2876 : !c........................................................................
2877 :
2878 : end do ! k loop
2879 :
2880 0 : prect1(i)=prect1(i)+prect(i)
2881 0 : preci1(i)=preci1(i)+preci(i)
2882 :
2883 : end do ! it loop, sub-step
2884 :
2885 0 : do k = top_lev, pver
2886 0 : rate1ord_cw2pr_st(i,k) = qcsinksum_rate1ord(k)/max(qcsum_rate1ord(k),1.0e-30_r8)
2887 : end do
2888 :
2889 0 : 300 continue ! continue if no cloud water
2890 : end do ! i loop
2891 :
2892 : ! convert dt from sub-step back to full time step
2893 0 : deltat=deltat*real(iter)
2894 :
2895 : !c.............................................................................
2896 :
2897 0 : do i=1,ncol
2898 :
2899 : ! skip all calculations if no cloud water
2900 0 : if (ltrue(i).eq.0) then
2901 :
2902 0 : do k=1,top_lev-1
2903 : ! assign zero values for effective radius above 1 mbar
2904 0 : effc(i,k)=0._r8
2905 0 : effi(i,k)=0._r8
2906 0 : effc_fn(i,k)=0._r8
2907 0 : lamcrad(i,k)=0._r8
2908 0 : pgamrad(i,k)=0._r8
2909 0 : deffi(i,k)=0._r8
2910 : end do
2911 :
2912 0 : do k=top_lev,pver
2913 : ! assign default values for effective radius
2914 0 : effc(i,k)=10._r8
2915 0 : effi(i,k)=25._r8
2916 0 : effc_fn(i,k)=10._r8
2917 0 : lamcrad(i,k)=0._r8
2918 0 : pgamrad(i,k)=0._r8
2919 0 : deffi(i,k)=0._r8
2920 : end do
2921 : goto 500
2922 : end if
2923 :
2924 : ! initialize nstep for sedimentation sub-steps
2925 0 : nstep = 1
2926 :
2927 : ! divide precip rate by number of sub-steps to get average over time step
2928 :
2929 0 : prect(i)=prect1(i)/real(iter)
2930 0 : preci(i)=preci1(i)/real(iter)
2931 :
2932 0 : do k=top_lev,pver
2933 :
2934 : ! assign variables back to start-of-timestep values before updating after sub-steps
2935 :
2936 0 : t(i,k)=t1(i,k)
2937 0 : q(i,k)=q1(i,k)
2938 0 : qc(i,k)=qc1(i,k)
2939 0 : qi(i,k)=qi1(i,k)
2940 0 : nc(i,k)=nc1(i,k)
2941 0 : ni(i,k)=ni1(i,k)
2942 :
2943 : ! divide microphysical tendencies by number of sub-steps to get average over time step
2944 :
2945 0 : tlat(i,k)=tlat1(i,k)/real(iter)
2946 0 : qvlat(i,k)=qvlat1(i,k)/real(iter)
2947 0 : qctend(i,k)=qctend1(i,k)/real(iter)
2948 0 : qitend(i,k)=qitend1(i,k)/real(iter)
2949 0 : nctend(i,k)=nctend1(i,k)/real(iter)
2950 0 : nitend(i,k)=nitend1(i,k)/real(iter)
2951 :
2952 0 : rainrt(i,k)=rainrt1(i,k)/real(iter)
2953 :
2954 : ! divide by number of sub-steps to find final values
2955 0 : rflx(i,k+1)=rflx1(i,k+1)/real(iter)
2956 0 : sflx(i,k+1)=sflx1(i,k+1)/real(iter)
2957 :
2958 : ! divide output precip q and N by number of sub-steps to get average over time step
2959 :
2960 0 : qrout(i,k)=qrout(i,k)/real(iter)
2961 0 : qsout(i,k)=qsout(i,k)/real(iter)
2962 0 : nrout(i,k)=nrout(i,k)/real(iter)
2963 0 : nsout(i,k)=nsout(i,k)/real(iter)
2964 :
2965 : ! divide trop_mozart variables by number of sub-steps to get average over time step
2966 :
2967 0 : nevapr(i,k) = nevapr(i,k)/real(iter)
2968 0 : nevapr2(i,k) = nevapr2(i,k)/real(iter)
2969 0 : evapsnow(i,k) = evapsnow(i,k)/real(iter)
2970 0 : prain(i,k) = prain(i,k)/real(iter)
2971 0 : prodsnow(i,k) = prodsnow(i,k)/real(iter)
2972 0 : cmeout(i,k) = cmeout(i,k)/real(iter)
2973 :
2974 0 : cmeiout(i,k) = cmeiout(i,k)/real(iter)
2975 0 : meltsdt(i,k) = meltsdt(i,k)/real(iter)
2976 0 : frzrdt (i,k) = frzrdt (i,k)/real(iter)
2977 :
2978 :
2979 : ! microphysics output
2980 0 : prao(i,k)=prao(i,k)/real(iter)
2981 0 : prco(i,k)=prco(i,k)/real(iter)
2982 0 : mnuccco(i,k)=mnuccco(i,k)/real(iter)
2983 0 : mnuccto(i,k)=mnuccto(i,k)/real(iter)
2984 0 : msacwio(i,k)=msacwio(i,k)/real(iter)
2985 0 : psacwso(i,k)=psacwso(i,k)/real(iter)
2986 0 : bergso(i,k)=bergso(i,k)/real(iter)
2987 0 : bergo(i,k)=bergo(i,k)/real(iter)
2988 0 : prcio(i,k)=prcio(i,k)/real(iter)
2989 0 : praio(i,k)=praio(i,k)/real(iter)
2990 :
2991 0 : mnuccro(i,k)=mnuccro(i,k)/real(iter)
2992 0 : pracso (i,k)=pracso (i,k)/real(iter)
2993 :
2994 0 : mnuccdo(i,k)=mnuccdo(i,k)/real(iter)
2995 :
2996 : ! modify to include snow. in prain & evap (diagnostic here: for wet dep)
2997 0 : nevapr(i,k) = nevapr(i,k) + evapsnow(i,k)
2998 0 : prer_evap(i,k) = nevapr2(i,k)
2999 0 : prain(i,k) = prain(i,k) + prodsnow(i,k)
3000 :
3001 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3002 : ! calculate sedimentation for cloud water and ice
3003 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3004 :
3005 : ! update in-cloud cloud mixing ratio and number concentration
3006 : ! with microphysical tendencies to calculate sedimentation, assign to dummy vars
3007 : ! note: these are in-cloud values***, hence we divide by cloud fraction
3008 :
3009 0 : dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k)
3010 0 : dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)/icldm(i,k)
3011 0 : dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k),0._r8)
3012 0 : dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)/icldm(i,k),0._r8)
3013 :
3014 0 : if (nccons) then
3015 0 : dumnc(i,k) = ncnst/rho(i,k)
3016 : end if
3017 0 : if (nicons) then
3018 0 : dumni(i,k) = ninst/rho(i,k)
3019 : end if
3020 :
3021 : ! obtain new slope parameter to avoid possible singularity
3022 :
3023 0 : if (dumi(i,k).ge.qsmall) then
3024 : ! add upper limit to in-cloud number concentration to prevent numerical error
3025 0 : dumni(i,k)=min(dumni(i,k),dumi(i,k)*1.e20_r8)
3026 :
3027 : lami(k) = (cons1*ci* &
3028 0 : dumni(i,k)/dumi(i,k))**(1._r8/di)
3029 0 : lami(k)=max(lami(k),lammini)
3030 0 : lami(k)=min(lami(k),lammaxi)
3031 : else
3032 0 : lami(k)=0._r8
3033 : end if
3034 :
3035 0 : if (dumc(i,k).ge.qsmall) then
3036 : ! add upper limit to in-cloud number concentration to prevent numerical error
3037 0 : dumnc(i,k)=min(dumnc(i,k),dumc(i,k)*1.e20_r8)
3038 : ! add lower limit to in-cloud number concentration
3039 0 : dumnc(i,k)=max(dumnc(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3
3040 0 : pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
3041 0 : pgam(k)=1._r8/(pgam(k)**2)-1._r8
3042 0 : pgam(k)=max(pgam(k),2._r8)
3043 0 : pgam(k)=min(pgam(k),15._r8)
3044 :
3045 : lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
3046 0 : (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
3047 0 : lammin = (pgam(k)+1._r8)/50.e-6_r8
3048 0 : lammax = (pgam(k)+1._r8)/2.e-6_r8
3049 0 : lamc(k)=max(lamc(k),lammin)
3050 0 : lamc(k)=min(lamc(k),lammax)
3051 : else
3052 0 : lamc(k)=0._r8
3053 : end if
3054 :
3055 : ! calculate number and mass weighted fall velocity for droplets
3056 : ! include effects of sub-grid distribution of cloud water
3057 :
3058 :
3059 0 : if (dumc(i,k).ge.qsmall) then
3060 0 : unc = acn(i,k)*gamma(1._r8+bc+pgam(k))/(lamc(k)**bc*gamma(pgam(k)+1._r8))
3061 0 : umc = acn(i,k)*gamma(4._r8+bc+pgam(k))/(lamc(k)**bc*gamma(pgam(k)+4._r8))
3062 : ! fallspeed for output
3063 0 : vtrmc(i,k)=umc
3064 : else
3065 : umc = 0._r8
3066 : unc = 0._r8
3067 : end if
3068 :
3069 : ! calculate number and mass weighted fall velocity for cloud ice
3070 :
3071 0 : if (dumi(i,k).ge.qsmall) then
3072 0 : uni = ain(i,k)*cons16/lami(k)**bi
3073 0 : umi = ain(i,k)*cons17/(6._r8*lami(k)**bi)
3074 0 : uni=min(uni,1.2_r8*rhof(i,k))
3075 0 : umi=min(umi,1.2_r8*rhof(i,k))
3076 :
3077 : ! fallspeed
3078 0 : vtrmi(i,k)=umi
3079 : else
3080 : umi = 0._r8
3081 : uni = 0._r8
3082 : end if
3083 :
3084 0 : fi(k) = g*rho(i,k)*umi
3085 0 : fni(k) = g*rho(i,k)*uni
3086 0 : fc(k) = g*rho(i,k)*umc
3087 0 : fnc(k) = g*rho(i,k)*unc
3088 :
3089 : ! calculate number of split time steps to ensure courant stability criteria
3090 : ! for sedimentation calculations
3091 :
3092 0 : rgvm = max(fi(k),fc(k),fni(k),fnc(k))
3093 0 : nstep = max(int(rgvm*deltat/pdel(i,k)+1._r8),nstep)
3094 :
3095 : ! redefine dummy variables - sedimentation is calculated over grid-scale
3096 : ! quantities to ensure conservation
3097 :
3098 0 : dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)
3099 0 : dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)
3100 0 : dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat),0._r8)
3101 0 : dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat),0._r8)
3102 :
3103 0 : if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
3104 0 : if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8
3105 :
3106 : end do !!! vertical loop
3107 0 : do n = 1,nstep !! loop over sub-time step to ensure stability
3108 :
3109 0 : do k = top_lev,pver
3110 0 : if (do_cldice) then
3111 0 : falouti(k) = fi(k)*dumi(i,k)
3112 0 : faloutni(k) = fni(k)*dumni(i,k)
3113 : else
3114 0 : falouti(k) = 0._r8
3115 0 : faloutni(k) = 0._r8
3116 : end if
3117 :
3118 0 : faloutc(k) = fc(k)*dumc(i,k)
3119 0 : faloutnc(k) = fnc(k)*dumnc(i,k)
3120 : end do
3121 :
3122 : ! top of model
3123 :
3124 0 : k = top_lev
3125 0 : faltndi = falouti(k)/pdel(i,k)
3126 0 : faltndni = faloutni(k)/pdel(i,k)
3127 0 : faltndc = faloutc(k)/pdel(i,k)
3128 0 : faltndnc = faloutnc(k)/pdel(i,k)
3129 :
3130 : ! add fallout terms to microphysical tendencies
3131 :
3132 0 : qitend(i,k) = qitend(i,k)-faltndi/nstep
3133 0 : nitend(i,k) = nitend(i,k)-faltndni/nstep
3134 0 : qctend(i,k) = qctend(i,k)-faltndc/nstep
3135 0 : nctend(i,k) = nctend(i,k)-faltndnc/nstep
3136 :
3137 : ! sedimentation tendencies for output
3138 0 : qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep
3139 0 : qisedten(i,k)=qisedten(i,k)-faltndi/nstep
3140 :
3141 0 : dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep
3142 0 : dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep
3143 0 : dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep
3144 0 : dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep
3145 :
3146 0 : do k = top_lev+1,pver
3147 :
3148 : ! for cloud liquid and ice, if cloud fraction increases with height
3149 : ! then add flux from above to both vapor and cloud water of current level
3150 : ! this means that flux entering clear portion of cell from above evaporates
3151 : ! instantly
3152 :
3153 0 : dum=lcldm(i,k)/lcldm(i,k-1)
3154 0 : dum=min(dum,1._r8)
3155 0 : dum1=icldm(i,k)/icldm(i,k-1)
3156 0 : dum1=min(dum1,1._r8)
3157 :
3158 0 : faltndqie=(falouti(k)-falouti(k-1))/pdel(i,k)
3159 0 : faltndi=(falouti(k)-dum1*falouti(k-1))/pdel(i,k)
3160 0 : faltndni=(faloutni(k)-dum1*faloutni(k-1))/pdel(i,k)
3161 0 : faltndqce=(faloutc(k)-faloutc(k-1))/pdel(i,k)
3162 0 : faltndc=(faloutc(k)-dum*faloutc(k-1))/pdel(i,k)
3163 0 : faltndnc=(faloutnc(k)-dum*faloutnc(k-1))/pdel(i,k)
3164 :
3165 : ! add fallout terms to eulerian tendencies
3166 :
3167 0 : qitend(i,k) = qitend(i,k)-faltndi/nstep
3168 0 : nitend(i,k) = nitend(i,k)-faltndni/nstep
3169 0 : qctend(i,k) = qctend(i,k)-faltndc/nstep
3170 0 : nctend(i,k) = nctend(i,k)-faltndnc/nstep
3171 :
3172 : ! sedimentation tendencies for output
3173 0 : qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep
3174 0 : qisedten(i,k)=qisedten(i,k)-faltndi/nstep
3175 :
3176 : ! add terms to to evap/sub of cloud water
3177 :
3178 0 : qvlat(i,k)=qvlat(i,k)-(faltndqie-faltndi)/nstep
3179 : ! for output
3180 0 : qisevap(i,k)=qisevap(i,k)-(faltndqie-faltndi)/nstep
3181 0 : qvlat(i,k)=qvlat(i,k)-(faltndqce-faltndc)/nstep
3182 : ! for output
3183 0 : qcsevap(i,k)=qcsevap(i,k)-(faltndqce-faltndc)/nstep
3184 :
3185 0 : tlat(i,k)=tlat(i,k)+(faltndqie-faltndi)*xxls/nstep
3186 0 : tlat(i,k)=tlat(i,k)+(faltndqce-faltndc)*xxlv/nstep
3187 :
3188 0 : dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep
3189 0 : dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep
3190 0 : dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep
3191 0 : dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep
3192 :
3193 0 : Fni(K)=MAX(Fni(K)/pdel(i,K),Fni(K-1)/pdel(i,K-1))*pdel(i,K)
3194 0 : FI(K)=MAX(FI(K)/pdel(i,K),FI(K-1)/pdel(i,K-1))*pdel(i,K)
3195 0 : fnc(k)=max(fnc(k)/pdel(i,k),fnc(k-1)/pdel(i,k-1))*pdel(i,k)
3196 0 : Fc(K)=MAX(Fc(K)/pdel(i,K),Fc(K-1)/pdel(i,K-1))*pdel(i,K)
3197 :
3198 : end do !! k loop
3199 :
3200 : ! units below are m/s
3201 : ! cloud water/ice sedimentation flux at surface
3202 : ! is added to precip flux at surface to get total precip (cloud + precip water)
3203 : ! rate
3204 :
3205 0 : prect(i) = prect(i)+(faloutc(pver)+falouti(pver))/g/nstep/1000._r8
3206 0 : preci(i) = preci(i)+(falouti(pver))/g/nstep/1000._r8
3207 :
3208 : ! Add fallout to Precip Flux: note unit change m/s *kg/m3 = kg/m2
3209 0 : do k = top_lev,pver
3210 0 : rflx(i,k+1)=rflx(i,k+1)+(faloutc(k))/g/nstep/1000._r8 * rhow
3211 0 : sflx(i,k+1)=sflx(i,k+1)+(falouti(k))/g/nstep/1000._r8 * rhow
3212 : end do
3213 :
3214 : end do !! nstep loop
3215 :
3216 : ! end sedimentation
3217 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3218 :
3219 : ! get new update for variables that includes sedimentation tendency
3220 : ! note : here dum variables are grid-average, NOT in-cloud
3221 :
3222 0 : do k=top_lev,pver
3223 :
3224 0 : dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)
3225 0 : dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)
3226 0 : dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)
3227 0 : dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)
3228 :
3229 0 : if (nccons) then
3230 0 : dumnc(i,k) = ncnst/rho(i,k)*lcldm(i,k)
3231 : end if
3232 0 : if (nicons) then
3233 0 : dumni(i,k) = ninst/rho(i,k)*icldm(i,k)
3234 : end if
3235 :
3236 0 : if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
3237 0 : if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8
3238 :
3239 : ! calculate instantaneous processes (melting, homogeneous freezing)
3240 0 : if (do_cldice) then
3241 :
3242 0 : if (t(i,k)+tlat(i,k)/cpp*deltat > tmelt) then
3243 0 : if (dumi(i,k) > 0._r8) then
3244 :
3245 : ! limit so that melting does not push temperature below freezing
3246 0 : dum = -dumi(i,k)*xlf/cpp
3247 0 : if (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt.tmelt) then
3248 0 : dum = (t(i,k)+tlat(i,k)/cpp*deltat-tmelt)*cpp/xlf
3249 0 : dum = dum/dumi(i,k)*xlf/cpp
3250 0 : dum = max(0._r8,dum)
3251 0 : dum = min(1._r8,dum)
3252 : else
3253 : dum = 1._r8
3254 : end if
3255 :
3256 0 : qctend(i,k)=qctend(i,k)+dum*dumi(i,k)/deltat
3257 :
3258 : ! for output
3259 0 : melto(i,k)=dum*dumi(i,k)/deltat
3260 :
3261 : ! assume melting ice produces droplet
3262 : ! mean volume radius of 8 micron
3263 :
3264 : nctend(i,k)=nctend(i,k)+3._r8*dum*dumi(i,k)/deltat/ &
3265 0 : (4._r8*pi*5.12e-16_r8*rhow)
3266 :
3267 0 : qitend(i,k)=((1._r8-dum)*dumi(i,k)-qi(i,k))/deltat
3268 0 : nitend(i,k)=((1._r8-dum)*dumni(i,k)-ni(i,k))/deltat
3269 0 : tlat(i,k)=tlat(i,k)-xlf*dum*dumi(i,k)/deltat
3270 : end if
3271 : end if
3272 :
3273 : ! homogeneously freeze droplets at -40 C
3274 :
3275 0 : if (t(i,k)+tlat(i,k)/cpp*deltat < 233.15_r8) then
3276 0 : if (dumc(i,k) > 0._r8) then
3277 :
3278 : ! limit so that freezing does not push temperature above threshold
3279 0 : dum = dumc(i,k)*xlf/cpp
3280 0 : if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.233.15_r8) then
3281 0 : dum = -(t(i,k)+tlat(i,k)/cpp*deltat-233.15_r8)*cpp/xlf
3282 0 : dum = dum/dumc(i,k)*xlf/cpp
3283 0 : dum = max(0._r8,dum)
3284 0 : dum = min(1._r8,dum)
3285 : else
3286 : dum = 1._r8
3287 : end if
3288 :
3289 0 : qitend(i,k)=qitend(i,k)+dum*dumc(i,k)/deltat
3290 : ! for output
3291 0 : homoo(i,k)=dum*dumc(i,k)/deltat
3292 :
3293 : ! assume 25 micron mean volume radius of homogeneously frozen droplets
3294 : ! consistent with size of detrained ice in stratiform.F90
3295 : nitend(i,k)=nitend(i,k)+dum*3._r8*dumc(i,k)/(4._r8*3.14_r8*1.563e-14_r8* &
3296 0 : 500._r8)/deltat
3297 0 : qctend(i,k)=((1._r8-dum)*dumc(i,k)-qc(i,k))/deltat
3298 0 : nctend(i,k)=((1._r8-dum)*dumnc(i,k)-nc(i,k))/deltat
3299 0 : tlat(i,k)=tlat(i,k)+xlf*dum*dumc(i,k)/deltat
3300 : end if
3301 : end if
3302 :
3303 : ! remove any excess over-saturation, which is possible due to non-linearity when adding
3304 : ! together all microphysical processes
3305 : ! follow code similar to old CAM scheme
3306 :
3307 0 : qtmp=q(i,k)+qvlat(i,k)*deltat
3308 0 : ttmp=t(i,k)+tlat(i,k)/cpp*deltat
3309 :
3310 0 : esn = svp_water(ttmp) ! use rhw to allow ice supersaturation
3311 0 : qsn = svp_to_qsat(esn, p(i,k))
3312 :
3313 0 : if (qtmp > qsn .and. qsn > 0) then
3314 : ! expression below is approximate since there may be ice deposition
3315 0 : dum = (qtmp-qsn)/(1._r8+cons27*qsn/(cpp*rv*ttmp**2))/deltat
3316 : ! add to output cme
3317 0 : cmeout(i,k) = cmeout(i,k)+dum
3318 : ! now add to tendencies, partition between liquid and ice based on temperature
3319 0 : if (ttmp > 268.15_r8) then
3320 : dum1=0.0_r8
3321 : ! now add to tendencies, partition between liquid and ice based on te
3322 0 : else if (ttmp < 238.15_r8) then
3323 : dum1=1.0_r8
3324 : else
3325 0 : dum1=(268.15_r8-ttmp)/30._r8
3326 : end if
3327 :
3328 : dum = (qtmp-qsn)/(1._r8+(xxls*dum1+xxlv*(1._r8-dum1))**2 &
3329 0 : *qsn/(cpp*rv*ttmp**2))/deltat
3330 0 : qctend(i,k)=qctend(i,k)+dum*(1._r8-dum1)
3331 : ! for output
3332 0 : qcreso(i,k)=dum*(1._r8-dum1)
3333 0 : qitend(i,k)=qitend(i,k)+dum*dum1
3334 0 : qireso(i,k)=dum*dum1
3335 0 : qvlat(i,k)=qvlat(i,k)-dum
3336 : ! for output
3337 0 : qvres(i,k)=-dum
3338 0 : tlat(i,k)=tlat(i,k)+dum*(1._r8-dum1)*xxlv+dum*dum1*xxls
3339 : end if
3340 : end if
3341 :
3342 : !...............................................................................
3343 : ! calculate effective radius for pass to radiation code
3344 : ! if no cloud water, default value is 10 micron for droplets,
3345 : ! 25 micron for cloud ice
3346 :
3347 : ! update cloud variables after instantaneous processes to get effective radius
3348 : ! variables are in-cloud to calculate size dist parameters
3349 :
3350 0 : dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)/lcldm(i,k)
3351 0 : dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)/icldm(i,k)
3352 0 : dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)/lcldm(i,k)
3353 0 : dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)/icldm(i,k)
3354 :
3355 0 : if (nccons) then
3356 0 : dumnc(i,k) = ncnst/rho(i,k)
3357 : end if
3358 0 : if (nicons) then
3359 0 : dumni(i,k) = ninst/rho(i,k)
3360 : end if
3361 :
3362 : ! limit in-cloud mixing ratio to reasonable value of 5 g kg-1
3363 :
3364 0 : dumc(i,k)=min(dumc(i,k),5.e-3_r8)
3365 0 : dumi(i,k)=min(dumi(i,k),5.e-3_r8)
3366 :
3367 : !...................
3368 : ! cloud ice effective radius
3369 :
3370 0 : if (dumi(i,k).ge.qsmall) then
3371 :
3372 0 : if (nicons) then
3373 : ! make sure ni is consistent with the constant N by adjusting
3374 : ! tendency, need to multiply by cloud fraction
3375 : ! note that nitend may be further adjusted below if mean crystal
3376 : ! size is out of bounds
3377 0 : nitend(i,k) = (ninst/rho(i,k)*icldm(i,k) - ni(i,k))/deltat
3378 : end if
3379 :
3380 : ! add upper limit to in-cloud number concentration to prevent numerical error
3381 0 : dumni(i,k)=min(dumni(i,k),dumi(i,k)*1.e20_r8)
3382 0 : lami(k) = (cons1*ci*dumni(i,k)/dumi(i,k))**(1._r8/di)
3383 :
3384 0 : if (lami(k).lt.lammini) then
3385 0 : lami(k) = lammini
3386 0 : n0i(k) = lami(k)**(di+1._r8)*dumi(i,k)/(ci*cons1)
3387 0 : niic(i,k) = n0i(k)/lami(k)
3388 : ! adjust number conc if needed to keep mean size in reasonable range
3389 0 : if (do_cldice) nitend(i,k)=(niic(i,k)*icldm(i,k)-ni(i,k))/deltat
3390 :
3391 0 : else if (lami(k).gt.lammaxi) then
3392 0 : lami(k) = lammaxi
3393 0 : n0i(k) = lami(k)**(di+1._r8)*dumi(i,k)/(ci*cons1)
3394 0 : niic(i,k) = n0i(k)/lami(k)
3395 : ! adjust number conc if needed to keep mean size in reasonable range
3396 0 : if (do_cldice) nitend(i,k)=(niic(i,k)*icldm(i,k)-ni(i,k))/deltat
3397 : end if
3398 0 : effi(i,k) = 1.5_r8/lami(k)*1.e6_r8
3399 :
3400 : else
3401 0 : effi(i,k) = 25._r8
3402 : end if
3403 :
3404 : ! NOTE: If CARMA is doing the ice microphysics, then the ice effective
3405 : ! radius has already been determined from the size distribution.
3406 0 : if (.not. do_cldice) then
3407 0 : effi(i,k) = re_ice(i,k) * 1e6_r8 ! m -> um
3408 : end if
3409 :
3410 : !...................
3411 : ! cloud droplet effective radius
3412 :
3413 0 : if (dumc(i,k).ge.qsmall) then
3414 :
3415 0 : if (nccons) then
3416 : ! make sure nc is consistent with the constant N by adjusting
3417 : ! tendency, need to multiply by cloud fraction
3418 : ! note that nctend may be further adjusted below if mean droplet
3419 : ! size is out of bounds
3420 0 : nctend(i,k) = (ncnst/rho(i,k)*lcldm(i,k) - nc(i,k))/deltat
3421 : end if
3422 :
3423 : ! add upper limit to in-cloud number concentration to prevent numerical error
3424 0 : dumnc(i,k)=min(dumnc(i,k),dumc(i,k)*1.e20_r8)
3425 :
3426 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3427 : ! set tendency to ensure minimum droplet concentration
3428 : ! after update by microphysics, except when lambda exceeds bounds on mean drop
3429 : ! size or if there is no cloud water
3430 0 : if (dumnc(i,k).lt.cdnl/rho(i,k)) then
3431 0 : nctend(i,k)=(cdnl/rho(i,k)*lcldm(i,k)-nc(i,k))/deltat
3432 : end if
3433 0 : dumnc(i,k)=max(dumnc(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3
3434 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3435 0 : pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
3436 0 : pgam(k)=1._r8/(pgam(k)**2)-1._r8
3437 0 : pgam(k)=max(pgam(k),2._r8)
3438 0 : pgam(k)=min(pgam(k),15._r8)
3439 :
3440 : lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
3441 0 : (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
3442 0 : lammin = (pgam(k)+1._r8)/50.e-6_r8
3443 : ! Multiply by omsm to fit within RRTMG's table.
3444 0 : lammax = (pgam(k)+1._r8)*omsm/2.e-6_r8
3445 0 : if (lamc(k).lt.lammin) then
3446 0 : lamc(k) = lammin
3447 : ncic(i,k) = 6._r8*lamc(k)**3*dumc(i,k)* &
3448 : gamma(pgam(k)+1._r8)/ &
3449 0 : (pi*rhow*gamma(pgam(k)+4._r8))
3450 : ! adjust number conc if needed to keep mean size in reasonable range
3451 0 : nctend(i,k)=(ncic(i,k)*lcldm(i,k)-nc(i,k))/deltat
3452 :
3453 0 : else if (lamc(k).gt.lammax) then
3454 0 : lamc(k) = lammax
3455 : ncic(i,k) = 6._r8*lamc(k)**3*dumc(i,k)* &
3456 : gamma(pgam(k)+1._r8)/ &
3457 0 : (pi*rhow*gamma(pgam(k)+4._r8))
3458 : ! adjust number conc if needed to keep mean size in reasonable range
3459 0 : nctend(i,k)=(ncic(i,k)*lcldm(i,k)-nc(i,k))/deltat
3460 : end if
3461 :
3462 : effc(i,k) = &
3463 : gamma(pgam(k)+4._r8)/ &
3464 0 : gamma(pgam(k)+3._r8)/lamc(k)/2._r8*1.e6_r8
3465 : !assign output fields for shape here
3466 0 : lamcrad(i,k)=lamc(k)
3467 0 : pgamrad(i,k)=pgam(k)
3468 :
3469 : else
3470 0 : effc(i,k) = 10._r8
3471 0 : lamcrad(i,k)=0._r8
3472 0 : pgamrad(i,k)=0._r8
3473 : end if
3474 :
3475 : ! ice effective diameter for david mitchell's optics
3476 0 : if (do_cldice) then
3477 0 : deffi(i,k)=effi(i,k)*rhoi/917._r8*2._r8
3478 : else
3479 0 : deffi(i,k)=effi(i,k) * 2._r8
3480 : end if
3481 :
3482 :
3483 : !!! recalculate effective radius for constant number, in order to separate
3484 : ! first and second indirect effects
3485 : ! assume constant number of 10^8 kg-1
3486 :
3487 0 : dumnc(i,k)=1.e8_r8
3488 :
3489 0 : if (dumc(i,k).ge.qsmall) then
3490 0 : pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
3491 0 : pgam(k)=1._r8/(pgam(k)**2)-1._r8
3492 0 : pgam(k)=max(pgam(k),2._r8)
3493 0 : pgam(k)=min(pgam(k),15._r8)
3494 :
3495 : lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
3496 0 : (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
3497 0 : lammin = (pgam(k)+1._r8)/50.e-6_r8
3498 0 : lammax = (pgam(k)+1._r8)/2.e-6_r8
3499 0 : if (lamc(k).lt.lammin) then
3500 0 : lamc(k) = lammin
3501 0 : else if (lamc(k).gt.lammax) then
3502 0 : lamc(k) = lammax
3503 : end if
3504 : effc_fn(i,k) = &
3505 : gamma(pgam(k)+4._r8)/ &
3506 0 : gamma(pgam(k)+3._r8)/lamc(k)/2._r8*1.e6_r8
3507 :
3508 : else
3509 0 : effc_fn(i,k) = 10._r8
3510 : end if
3511 :
3512 :
3513 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1!
3514 :
3515 : end do ! vertical k loop
3516 :
3517 : 500 continue
3518 :
3519 0 : do k=top_lev,pver
3520 : ! if updated q (after microphysics) is zero, then ensure updated n is also zero
3521 :
3522 0 : if (qc(i,k)+qctend(i,k)*deltat.lt.qsmall) nctend(i,k)=-nc(i,k)/deltat
3523 0 : if (do_cldice .and. qi(i,k)+qitend(i,k)*deltat.lt.qsmall) nitend(i,k)=-ni(i,k)/deltat
3524 : end do
3525 :
3526 : end do ! i loop
3527 :
3528 : ! add snow ouptut
3529 0 : do i = 1,ncol
3530 0 : do k=top_lev,pver
3531 0 : if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then
3532 0 : dsout(i,k)=3._r8*rhosn/917._r8*(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)
3533 : endif
3534 : end do
3535 : end do
3536 :
3537 : !calculate effective radius of rain and snow in microns for COSP using Eq. 9 of COSP v1.3 manual
3538 0 : do i = 1,ncol
3539 0 : do k=top_lev,pver
3540 : !! RAIN
3541 0 : if (qrout(i,k).gt.1.e-7_r8.and.nrout(i,k).gt.0._r8) then
3542 0 : reff_rain(i,k)=1.5_r8*(pi * rhow * nrout(i,k)/qrout(i,k))**(-1._r8/3._r8)*1.e6_r8
3543 : endif
3544 : !! SNOW
3545 0 : if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then
3546 0 : reff_snow(i,k)=1.5_r8*(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)*1.e6_r8
3547 : end if
3548 : end do
3549 : end do
3550 :
3551 : ! analytic radar reflectivity
3552 : ! formulas from Matthew Shupe, NOAA/CERES
3553 : ! *****note: radar reflectivity is local (in-precip average)
3554 : ! units of mm^6/m^3
3555 :
3556 0 : do i = 1,ncol
3557 0 : do k=top_lev,pver
3558 0 : if (qc(i,k)+qctend(i,k)*deltat.ge.qsmall .and. nc(i,k)+nctend(i,k)*deltat.gt.10._r8) then
3559 : dum=((qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)*1000._r8)**2 &
3560 0 : /(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)/1.e6_r8)*lcldm(i,k)/cldmax(i,k)
3561 : else
3562 : dum=0._r8
3563 : end if
3564 0 : if (qi(i,k)+qitend(i,k)*deltat.ge.qsmall) then
3565 0 : dum1=((qi(i,k)+qitend(i,k)*deltat)*rho(i,k)/icldm(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)*icldm(i,k)/cldmax(i,k)
3566 : else
3567 : dum1=0._r8
3568 : end if
3569 :
3570 0 : if (qsout(i,k).ge.qsmall) then
3571 0 : dum1=dum1+(qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)
3572 : end if
3573 :
3574 0 : refl(i,k)=dum+dum1
3575 :
3576 : ! add rain rate, but for 37 GHz formulation instead of 94 GHz
3577 : ! formula approximated from data of Matrasov (2007)
3578 : ! rainrt is the rain rate in mm/hr
3579 : ! reflectivity (dum) is in DBz
3580 :
3581 0 : if (rainrt(i,k).ge.0.001_r8) then
3582 0 : dum=log10(rainrt(i,k)**6._r8)+16._r8
3583 :
3584 : ! convert from DBz to mm^6/m^3
3585 :
3586 0 : dum = 10._r8**(dum/10._r8)
3587 : else
3588 : ! don't include rain rate in R calculation for values less than 0.001 mm/hr
3589 : dum=0._r8
3590 : end if
3591 :
3592 : ! add to refl
3593 :
3594 0 : refl(i,k)=refl(i,k)+dum
3595 :
3596 : !output reflectivity in Z.
3597 0 : areflz(i,k)=refl(i,k)
3598 :
3599 : ! convert back to DBz
3600 :
3601 0 : if (refl(i,k).gt.minrefl) then
3602 0 : refl(i,k)=10._r8*log10(refl(i,k))
3603 : else
3604 0 : refl(i,k)=-9999._r8
3605 : end if
3606 :
3607 : !set averaging flag
3608 0 : if (refl(i,k).gt.mindbz) then
3609 0 : arefl(i,k)=refl(i,k)
3610 0 : frefl(i,k)=1.0_r8
3611 : else
3612 0 : arefl(i,k)=0._r8
3613 0 : areflz(i,k)=0._r8
3614 0 : frefl(i,k)=0._r8
3615 : end if
3616 :
3617 : ! bound cloudsat reflectivity
3618 :
3619 0 : csrfl(i,k)=min(csmax,refl(i,k))
3620 :
3621 : !set averaging flag
3622 0 : if (csrfl(i,k).gt.csmin) then
3623 0 : acsrfl(i,k)=refl(i,k)
3624 0 : fcsrfl(i,k)=1.0_r8
3625 : else
3626 0 : acsrfl(i,k)=0._r8
3627 0 : fcsrfl(i,k)=0._r8
3628 : end if
3629 :
3630 : end do
3631 : end do
3632 :
3633 :
3634 : ! averaging for snow and rain number and diameter
3635 :
3636 0 : qrout2(:,:)=0._r8
3637 0 : qsout2(:,:)=0._r8
3638 0 : nrout2(:,:)=0._r8
3639 0 : nsout2(:,:)=0._r8
3640 0 : drout2(:,:)=0._r8
3641 0 : dsout2(:,:)=0._r8
3642 0 : freqs(:,:)=0._r8
3643 0 : freqr(:,:)=0._r8
3644 0 : do i = 1,ncol
3645 0 : do k=top_lev,pver
3646 0 : if (qrout(i,k).gt.1.e-7_r8.and.nrout(i,k).gt.0._r8) then
3647 0 : qrout2(i,k)=qrout(i,k)
3648 0 : nrout2(i,k)=nrout(i,k)
3649 0 : drout2(i,k)=(pi * rhow * nrout(i,k)/qrout(i,k))**(-1._r8/3._r8)
3650 0 : freqr(i,k)=1._r8
3651 : endif
3652 0 : if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then
3653 0 : qsout2(i,k)=qsout(i,k)
3654 0 : nsout2(i,k)=nsout(i,k)
3655 0 : dsout2(i,k)=(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)
3656 0 : freqs(i,k)=1._r8
3657 : endif
3658 : end do
3659 : end do
3660 :
3661 : ! output activated liquid and ice (convert from #/kg -> #/m3)
3662 0 : do i = 1,ncol
3663 0 : do k=top_lev,pver
3664 0 : ncai(i,k)=dum2i(i,k)*rho(i,k)
3665 0 : ncal(i,k)=dum2l(i,k)*rho(i,k)
3666 : end do
3667 : end do
3668 :
3669 :
3670 : !redefine fice here....
3671 0 : nfice(:,:)=0._r8
3672 0 : do k=top_lev,pver
3673 0 : do i=1,ncol
3674 0 : dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)
3675 0 : dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)
3676 0 : dumfice=qsout(i,k) + qrout(i,k) + dumc(i,k) + dumi(i,k)
3677 :
3678 0 : if (dumfice.gt.qsmall.and.(qsout(i,k)+dumi(i,k).gt.qsmall)) then
3679 0 : nfice(i,k)=(qsout(i,k) + dumi(i,k))/dumfice
3680 : endif
3681 :
3682 0 : if (nfice(i,k).gt.1._r8) then
3683 0 : nfice(i,k)=1._r8
3684 : endif
3685 :
3686 : enddo
3687 : enddo
3688 :
3689 :
3690 0 : end subroutine micro_mg_tend
3691 :
3692 : !========================================================================
3693 : !UTILITIES
3694 : !========================================================================
3695 :
3696 0 : pure subroutine micro_mg_get_cols(ncol, nlev, top_lev, qcn, qin, &
3697 : mgncol, mgcols)
3698 :
3699 : ! Determines which columns microphysics should operate over by
3700 : ! checking for non-zero cloud water/ice.
3701 :
3702 : integer, intent(in) :: ncol ! Number of columns with meaningful data
3703 : integer, intent(in) :: nlev ! Number of levels to use
3704 : integer, intent(in) :: top_lev ! Top level for microphysics
3705 :
3706 : real(r8), intent(in) :: qcn(:,:) ! cloud water mixing ratio (kg/kg)
3707 : real(r8), intent(in) :: qin(:,:) ! cloud ice mixing ratio (kg/kg)
3708 :
3709 : integer, intent(out) :: mgncol ! Number of columns MG will use
3710 : integer, allocatable, intent(out) :: mgcols(:) ! column indices
3711 :
3712 : integer :: lev_offset ! top_lev - 1 (defined here for consistency)
3713 0 : logical :: ltrue(ncol) ! store tests for each column
3714 :
3715 : integer :: i, ii ! column indices
3716 :
3717 0 : if (allocated(mgcols)) deallocate(mgcols)
3718 :
3719 0 : lev_offset = top_lev - 1
3720 :
3721 : ! Using "any" along dimension 2 collapses across levels, but
3722 : ! not columns, so we know if water is present at any level
3723 : ! in each column.
3724 :
3725 0 : ltrue = any(qcn(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2)
3726 0 : ltrue = ltrue .or. any(qin(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2)
3727 :
3728 : ! Scan for true values to get a usable list of indices.
3729 :
3730 0 : mgncol = count(ltrue)
3731 0 : allocate(mgcols(mgncol))
3732 0 : i = 0
3733 0 : do ii = 1,ncol
3734 0 : if (ltrue(ii)) then
3735 0 : i = i + 1
3736 0 : mgcols(i) = ii
3737 : end if
3738 : end do
3739 :
3740 0 : end subroutine micro_mg_get_cols
3741 :
3742 : end module micro_mg1_0
|
