CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
micro_mg2_0.F90
1
4
24!---------------------------------------------------------------------------------
25!
26! NOTE: Modified to allow other microphysics packages (e.g. CARMA) to do ice
27! microphysics in cooperation with the MG liquid microphysics. This is
28! controlled by the do_cldice variable.
29!
30! If do_cldice is false, then MG microphysics should not update CLDICE or
31! NUMICE; it is assumed that the other microphysics scheme will have updated
32! CLDICE and NUMICE. The other microphysics should handle the following
33! processes that would have been done by MG:
34! - Detrainment (liquid and ice)
35! - Homogeneous ice nucleation
36! - Heterogeneous ice nucleation
37! - Bergeron process
38! - Melting of ice
39! - Freezing of cloud drops
40! - Autoconversion (ice -> snow)
41! - Growth/Sublimation of ice
42! - Sedimentation of ice
43!
44! This option has not been updated since the introduction of prognostic
45! precipitation, and probably should be adjusted to cover snow as well.
46!
47!---------------------------------------------------------------------------------
48! Based on micro_mg (restructuring of former cldwat2m_micro)
49! Author: Andrew Gettelman, Hugh Morrison.
50! Contributions from: Xiaohong Liu and Steve Ghan
51! December 2005-May 2010
52! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
53! Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)
54! for questions contact Hugh Morrison, Andrew Gettelman
55! e-mail: morrison@ucar.edu, andrew@ucar.edu
56!---------------------------------------------------------------------------------
57! Code comments added by HM, 093011
58! General code structure:
59!
60! Code is divided into two main subroutines:
61! subroutine micro_mg_init --> initializes microphysics routine, should be called
62! once at start of simulation
63! subroutine micro_mg_tend --> main microphysics routine to be called each time step
64! this also calls several smaller subroutines to calculate
65! microphysical processes and other utilities
66!
67! List of external functions:
68! qsat_water --> for calculating saturation vapor pressure with respect to liquid water
69! qsat_ice --> for calculating saturation vapor pressure with respect to ice
70! gamma --> standard mathematical gamma function
71! .........................................................................
72! List of inputs through use statement in fortran90:
73! Variable Name Description Units
74! .........................................................................
75! gravit acceleration due to gravity m s-2
76! rair dry air gas constant for air J kg-1 K-1
77! tmelt temperature of melting point for water K
78! cpair specific heat at constant pressure for dry air J kg-1 K-1
79! rh2o gas constant for water vapor J kg-1 K-1
80! latvap latent heat of vaporization J kg-1
81! latice latent heat of fusion J kg-1
82! qsat_water external function for calculating liquid water
83! saturation vapor pressure/humidity -
84! qsat_ice external function for calculating ice
85! saturation vapor pressure/humidity pa
86! rhmini relative humidity threshold parameter for
87! nucleating ice -
88! .........................................................................
89! NOTE: List of all inputs/outputs passed through the call/subroutine statement
90! for micro_mg_tend is given below at the start of subroutine micro_mg_tend.
91!---------------------------------------------------------------------------------
92
93! Procedures required:
94! 1) An implementation of the gamma function (if not intrinsic).
95! 2) saturation vapor pressure and specific humidity over water
96! 3) svp over ice
97use machine, only : r8 => kind_phys
98use funcphys, only : fpvsl, fpvsi
99
100!use wv_sat_methods, only: &
101! qsat_water => wv_sat_qsat_water, &
102! qsat_ice => wv_sat_qsat_ice
103
104! Parameters from the utilities module.
105use micro_mg_utils, only : pi, omsm, qsmall, mincld, rhosn, rhoi, &
106 rhow, rhows, ac, bc, ai, bi, &
107 aj, bj, ar, br, as, bs, &
109
110implicit none
111private
112save
113
114public :: micro_mg_init, micro_mg_tend, qcvar
115
116! Switches for specification rather than prediction of droplet and crystal number
117! note: number will be adjusted as needed to keep mean size within bounds,
118! even when specified droplet or ice number is used
119!
120! If constant cloud ice number is set (nicons = .true.),
121! then all microphysical processes except mass transfer due to ice nucleation
122! (mnuccd) are based on the fixed cloud ice number. Calculation of
123! mnuccd follows from the prognosed ice crystal number ni.
124
125logical :: nccons ! nccons = .true. to specify constant cloud droplet number
126logical :: nicons ! nicons = .true. to specify constant cloud ice number
127
128! specified ice and droplet number concentrations
129! note: these are local in-cloud values, not grid-mean
130real(r8) :: ncnst ! droplet num concentration when nccons=.true. (m-3)
131real(r8) :: ninst ! ice num concentration when nicons=.true. (m-3)
132
133!=========================================================
134! Private module parameters
135!=========================================================
136
137!Range of cloudsat reflectivities (dBz) for analytic simulator
138real(r8), parameter :: csmin = -30._r8
139real(r8), parameter :: csmax = 26._r8
140real(r8), parameter :: mindbz = -99._r8
141real(r8), parameter :: minrefl = 1.26e-10_r8 ! minrefl = 10._r8**(mindbz/10._r8)
142
143! autoconversion size threshold for cloud ice to snow (m)
144real(r8) :: dcs, ts_au, ts_au_min, qcvar
145
146! minimum mass of new crystal due to freezing of cloud droplets done
147! externally (kg)
148real(r8), parameter :: mi0l_min = 4._r8/3._r8*pi*rhow*(4.e-6_r8)**3
149real(r8), parameter :: zero=0.0_r8, one=1.0_r8, two=2.0_r8, three=3.0_r8, &
150 four=4.0_r8, five=5.0_r8, six=6._r8, half=0.5_r8, &
151 ten=10.0_r8, forty=40.0_r8, oneo6=one/six
152
153!=========================================================
154! Constants set in initialization
155!=========================================================
156
157! Set using arguments to micro_mg_init
158real(r8) :: g ! gravity
159real(r8) :: r ! dry air gas constant
160real(r8) :: rv ! water vapor gas constant
161real(r8) :: cpp ! specific heat of dry air
162real(r8) :: tmelt ! freezing point of water (K)
163
164! latent heats of:
165real(r8) :: xxlv ! vaporization
166real(r8) :: xlf ! freezing
167real(r8) :: xxls ! sublimation
168
169real(r8) :: rhmini ! Minimum rh for ice cloud fraction > 0.
170
171! flags
172logical :: microp_uniform, do_cldice, use_hetfrz_classnuc
173
174real(r8) :: rhosu ! typical 850mn air density
175
176real(r8) :: icenuct ! ice nucleation temperature: currently -5 degrees C
177
178real(r8) :: snowmelt ! what temp to melt all snow: currently 2 degrees C
179real(r8) :: rainfrze ! what temp to freeze all rain: currently -5 degrees C
180
181! additional constants to help speed up code
182real(r8) :: gamma_br_plus1, gamma_bs_plus1, gamma_bi_plus1, gamma_bj_plus1
183real(r8) :: gamma_br_plus4, gamma_bs_plus4, gamma_bi_plus4, gamma_bj_plus4
184real(r8) :: xxlv_squared, xxls_squared
185real(r8) :: omeps, epsqs
186
187character(len=16) :: micro_mg_precip_frac_method ! type of precipitation fraction method
188real(r8) :: micro_mg_berg_eff_factor ! berg efficiency factor
189
190logical :: allow_sed_supersat ! Allow supersaturated conditions after sedimentation loop
191logical :: do_sb_physics ! do SB 2001 autoconversion or accretion physics
192logical :: do_ice_gmao
193logical :: do_liq_liu
194
195!===============================================================================
196contains
197!===============================================================================
198
201subroutine micro_mg_init( &
202 kind, gravit, rair, rh2o, cpair, eps, &
203 tmelt_in, latvap, latice, &
204 rhmini_in, micro_mg_dcs,ts_auto, mg_qcvar, &
205 microp_uniform_in, do_cldice_in, use_hetfrz_classnuc_in, &
206 micro_mg_precip_frac_method_in, micro_mg_berg_eff_factor_in, &
207 allow_sed_supersat_in, do_sb_physics_in, &
208 do_ice_gmao_in, do_liq_liu_in, &
209 nccons_in, nicons_in, ncnst_in, ninst_in)
210
212 use wv_saturation, only : gestbl
213
214 !-----------------------------------------------------------------------
215 !
216 ! Purpose:
217 ! initialize constants for MG microphysics
218 !
219 ! Author: Andrew Gettelman Dec 2005
220 !
221 !-----------------------------------------------------------------------
222
223 integer, intent(in) :: kind
224 real(r8), intent(in) :: gravit
225 real(r8), intent(in) :: rair
226 real(r8), intent(in) :: rh2o
227 real(r8), intent(in) :: cpair
228 real(r8), intent(in) :: eps
229! real(r8), intent(in) :: fv
230 real(r8), intent(in) :: tmelt_in
231 real(r8), intent(in) :: latvap
232 real(r8), intent(in) :: latice
233 real(r8), intent(in) :: rhmini_in
234 real(r8), intent(in) :: micro_mg_dcs
235 real(r8), intent(in) :: ts_auto(2)
236 real(r8), intent(in) :: mg_qcvar
237
238 logical, intent(in) :: microp_uniform_in
240 logical, intent(in) :: do_cldice_in
242 logical, intent(in) :: use_hetfrz_classnuc_in
243
244 character(len=16),intent(in) :: micro_mg_precip_frac_method_in
245 real(r8), intent(in) :: micro_mg_berg_eff_factor_in
246 logical, intent(in) :: allow_sed_supersat_in
247 logical, intent(in) :: do_sb_physics_in
248 logical, intent(in) :: do_ice_gmao_in
249 logical, intent(in) :: do_liq_liu_in
250
251 logical, intent(in) :: nccons_in, nicons_in
252 real(r8), intent(in) :: ncnst_in, ninst_in
253 logical ip
254 real(r8):: tmn, tmx, trice
255
256
257
258 !-----------------------------------------------------------------------
259
260 dcs = micro_mg_dcs * 1.0e-6
261 ts_au_min = ts_auto(1)
262 ts_au = ts_auto(2)
263 qcvar = mg_qcvar
264
265 ! Initialize subordinate utilities module.
266 call micro_mg_utils_init(kind, rair, rh2o, cpair, tmelt_in, latvap, latice, &
267 dcs)
268
269
270 ! declarations for MG code (transforms variable names)
271
272 g = gravit ! gravity
273 r = rair ! dry air gas constant: note units(phys_constants are in J/K/kmol)
274 rv = rh2o ! water vapor gas constant
275 cpp = cpair ! specific heat of dry air
276 tmelt = tmelt_in
277 rhmini = rhmini_in
278 micro_mg_precip_frac_method = micro_mg_precip_frac_method_in
279 micro_mg_berg_eff_factor = micro_mg_berg_eff_factor_in
280 allow_sed_supersat = allow_sed_supersat_in
281 do_sb_physics = do_sb_physics_in
282 do_ice_gmao = do_ice_gmao_in
283 do_liq_liu = do_liq_liu_in
284
285 nccons = nccons_in
286 nicons = nicons_in
287 ncnst = ncnst_in
288 ninst = ninst_in
289
290 ! latent heats
291
292 xxlv = latvap ! latent heat vaporization
293 xlf = latice ! latent heat freezing
294 xxls = xxlv + xlf ! latent heat of sublimation
295
296 ! flags
297 microp_uniform = microp_uniform_in
298 do_cldice = do_cldice_in
299 use_hetfrz_classnuc = use_hetfrz_classnuc_in
300
301 ! typical air density at 850 mb
302
303 rhosu = 85000._r8 / (rair * tmelt)
304
305 ! Maximum temperature at which snow is allowed to exist
306 snowmelt = tmelt + two
307 ! Minimum temperature at which rain is allowed to exist
308 rainfrze = tmelt - forty
309
310 ! Ice nucleation temperature
311 icenuct = tmelt - five
312
313 ! Define constants to help speed up code (this limits calls to gamma function)
314 gamma_br_plus1 = gamma(br+one)
315 gamma_br_plus4 = gamma(br+four)
316 gamma_bs_plus1 = gamma(bs+one)
317 gamma_bs_plus4 = gamma(bs+four)
318 gamma_bi_plus1 = gamma(bi+one)
319 gamma_bi_plus4 = gamma(bi+four)
320 gamma_bj_plus1 = gamma(bj+one)
321 gamma_bj_plus4 = gamma(bj+four)
322
323 xxlv_squared = xxlv * xxlv
324 xxls_squared = xxls * xxls
325 epsqs = eps
326 omeps = one - epsqs
327 tmn = 173.16_r8
328 tmx = 375.16_r8
329 trice = 35.00_r8
330 ip = .true.
331 call gestbl(tmn ,tmx ,trice ,ip ,epsqs , latvap ,latice ,rh2o , &
332 cpair ,tmelt_in )
333
334
335
336end subroutine micro_mg_init
337
338!===============================================================================
339!microphysics routine for each timestep goes here...
340
341!\ingroup mg2_0_mp
346subroutine micro_mg_tend ( &
347 mgncol, nlev, deltatin, &
348 t, q, &
349 qcn, qin, &
350 ncn, nin, &
351 qrn, qsn, &
352 nrn, nsn, &
353 relvar, accre_enhan_i, &
354 p, pdel, &
355 cldn, liqcldf, icecldf, qsatfac, &
356 qcsinksum_rate1ord, &
357 naai, npccnin, &
358 rndst, nacon, &
359 tlat, qvlat, &
360 qctend, qitend, &
361 nctend, nitend, &
362 qrtend, qstend, &
363 nrtend, nstend, &
364 effc, effc_fn, effi, &
365 sadice, sadsnow, &
366 prect, preci, &
367 nevapr, evapsnow, &
368 am_evp_st, &
369 prain, prodsnow, &
370 cmeout, deffi, &
371 pgamrad, lamcrad, &
372 qsout, dsout, &
373 lflx, iflx, &
374 rflx, sflx, qrout, &
375 reff_rain, reff_snow, &
376 qcsevap, qisevap, qvres, &
377 cmeitot, vtrmc, vtrmi, &
378 umr, ums, &
379 qcsedten, qisedten, &
380 qrsedten, qssedten, &
381 pratot, prctot, &
382 mnuccctot, mnuccttot, msacwitot, &
383 psacwstot, bergstot, bergtot, &
384 melttot, homotot, &
385 qcrestot, prcitot, praitot, &
386 qirestot, mnuccrtot, pracstot, &
387 meltsdttot, frzrdttot, mnuccdtot, &
388 nrout, nsout, &
389 refl, arefl, areflz, &
390 frefl, csrfl, acsrfl, &
391 fcsrfl, rercld, &
392 ncai, ncal, &
393 qrout2, qsout2, &
394 nrout2, nsout2, &
395 drout2, dsout2, &
396 freqs, freqr, &
397 nfice, qcrat, &
398 prer_evap, xlat, xlon, lprnt, iccn, nlball)
399
400 ! Constituent properties.
401 use micro_mg_utils, only: mg_liq_props, &
402 mg_ice_props, &
403 mg_rain_props, &
404 mg_snow_props
405
406 ! Size calculation functions.
410
411 ! Microphysical processes.
432
433 !Authors: Hugh Morrison, Andrew Gettelman, NCAR, Peter Caldwell, LLNL
434 ! e-mail: morrison@ucar.edu, andrew@ucar.edu
435
436 ! input arguments
437 integer, intent(in) :: mgncol ! number of microphysics columns
438 integer, intent(in) :: nlev ! number of layers
439 integer, intent(in) :: nlball(mgncol) ! sedimentation start level
440 real(r8), intent(in) :: xlat,xlon ! number of layers
441 real(r8), intent(in) :: deltatin ! time step (s)
442 real(r8), intent(in) :: t(mgncol,nlev) ! input temperature (K)
443 real(r8), intent(in) :: q(mgncol,nlev) ! input h20 vapor mixing ratio (kg/kg)
444
445 ! note: all input cloud variables are grid-averaged
446 real(r8), intent(in) :: qcn(mgncol,nlev) ! cloud water mixing ratio (kg/kg)
447 real(r8), intent(in) :: qin(mgncol,nlev) ! cloud ice mixing ratio (kg/kg)
448 real(r8), intent(in) :: ncn(mgncol,nlev) ! cloud water number conc (1/kg)
449 real(r8), intent(in) :: nin(mgncol,nlev) ! cloud ice number conc (1/kg)
450
451 real(r8), intent(in) :: qrn(mgncol,nlev) ! rain mixing ratio (kg/kg)
452 real(r8), intent(in) :: qsn(mgncol,nlev) ! snow mixing ratio (kg/kg)
453 real(r8), intent(in) :: nrn(mgncol,nlev) ! rain number conc (1/kg)
454 real(r8), intent(in) :: nsn(mgncol,nlev) ! snow number conc (1/kg)
455
456 real(r8) :: relvar(mgncol,nlev) ! cloud water relative variance (-)
457 real(r8) :: accre_enhan(mgncol,nlev)! optional accretion
458! real(r8), intent(in) :: relvar_i ! cloud water relative variance (-)
459 real(r8), intent(in) :: accre_enhan_i ! optional accretion
460 ! enhancement factor (-)
461
462 real(r8), intent(in) :: p(mgncol,nlev) ! air pressure (pa)
463 real(r8), intent(in) :: pdel(mgncol,nlev) ! pressure difference across level (pa)
464
465 real(r8), intent(in) :: cldn(mgncol,nlev) ! cloud fraction (no units)
466 real(r8), intent(in) :: liqcldf(mgncol,nlev) ! liquid cloud fraction (no units)
467 real(r8), intent(in) :: icecldf(mgncol,nlev) ! ice cloud fraction (no units)
468 real(r8), intent(in) :: qsatfac(mgncol,nlev) ! subgrid cloud water saturation scaling factor (no units)
469 logical, intent(in) :: lprnt
470 integer, intent(in) :: iccn
471
472
473 ! used for scavenging
474 ! Inputs for aerosol activation
475 real(r8), intent(inout) :: naai(mgncol,nlev) ! ice nucleation number (from microp_aero_ts) (1/kg)
476 real(r8), intent(in) :: npccnin(mgncol,nlev) ! ccn activated number tendency (from microp_aero_ts) (1/kg*s)
477! real(r8), intent(in) :: npccn(mgncol,nlev) ! ccn activated number tendency (from microp_aero_ts) (1/kg*s)
478 real(r8) :: npccn(mgncol,nlev) ! ccn activated number tendency (from microp_aero_ts) (1/kg*s)
479
480 ! Note that for these variables, the dust bin is assumed to be the last index.
481 ! (For example, in CAM, the last dimension is always size 4.)
482 real(r8), intent(in) :: rndst(mgncol,nlev,10) ! radius of each dust bin, for contact freezing (from microp_aero_ts) (m)
483 real(r8), intent(in) :: nacon(mgncol,nlev,10) ! number in each dust bin, for contact freezing (from microp_aero_ts) (1/m^3)
484
485 ! output arguments
486
487 real(r8), intent(out) :: qcsinksum_rate1ord(mgncol,nlev) ! 1st order rate for
488 ! direct cw to precip conversion
489 real(r8), intent(out) :: tlat(mgncol,nlev) ! latent heating rate (W/kg)
490 real(r8), intent(out) :: qvlat(mgncol,nlev) ! microphysical tendency qv (1/s)
491 real(r8), intent(out) :: qctend(mgncol,nlev) ! microphysical tendency qc (1/s)
492 real(r8), intent(out) :: qitend(mgncol,nlev) ! microphysical tendency qi (1/s)
493 real(r8), intent(out) :: nctend(mgncol,nlev) ! microphysical tendency nc (1/(kg*s))
494 real(r8), intent(out) :: nitend(mgncol,nlev) ! microphysical tendency ni (1/(kg*s))
495
496 real(r8), intent(out) :: qrtend(mgncol,nlev) ! microphysical tendency qr (1/s)
497 real(r8), intent(out) :: qstend(mgncol,nlev) ! microphysical tendency qs (1/s)
498 real(r8), intent(out) :: nrtend(mgncol,nlev) ! microphysical tendency nr (1/(kg*s))
499 real(r8), intent(out) :: nstend(mgncol,nlev) ! microphysical tendency ns (1/(kg*s))
500 real(r8), intent(out) :: effc(mgncol,nlev) ! droplet effective radius (micron)
501 real(r8), intent(out) :: effc_fn(mgncol,nlev) ! droplet effective radius, assuming nc = 1.e8 kg-1
502 real(r8), intent(out) :: effi(mgncol,nlev) ! cloud ice effective radius (micron)
503 real(r8), intent(out) :: sadice(mgncol,nlev) ! cloud ice surface area density (cm2/cm3)
504 real(r8), intent(out) :: sadsnow(mgncol,nlev) ! cloud snow surface area density (cm2/cm3)
505 real(r8), intent(out) :: prect(mgncol) ! surface precip rate (m/s)
506 real(r8), intent(out) :: preci(mgncol) ! cloud ice/snow precip rate (m/s)
507 real(r8), intent(out) :: nevapr(mgncol,nlev) ! evaporation rate of rain + snow (1/s)
508 real(r8), intent(out) :: evapsnow(mgncol,nlev) ! sublimation rate of snow (1/s)
509 real(r8), intent(out) :: am_evp_st(mgncol,nlev) ! stratiform evaporation area (frac)
510 real(r8), intent(out) :: prain(mgncol,nlev) ! production of rain + snow (1/s)
511 real(r8), intent(out) :: prodsnow(mgncol,nlev) ! production of snow (1/s)
512 real(r8), intent(out) :: cmeout(mgncol,nlev) ! evap/sub of cloud (1/s)
513 real(r8), intent(out) :: deffi(mgncol,nlev) ! ice effective diameter for optics (radiation) (micron)
514 real(r8), intent(out) :: pgamrad(mgncol,nlev) ! ice gamma parameter for optics (radiation) (no units)
515 real(r8), intent(out) :: lamcrad(mgncol,nlev) ! slope of droplet distribution for optics (radiation) (1/m)
516 real(r8), intent(out) :: qsout(mgncol,nlev) ! snow mixing ratio (kg/kg)
517 real(r8), intent(out) :: dsout(mgncol,nlev) ! snow diameter (m)
518 real(r8), intent(out) :: lflx(mgncol,2:nlev+1) ! grid-box average liquid condensate flux (kg m^-2 s^-1)
519 real(r8), intent(out) :: iflx(mgncol,2:nlev+1) ! grid-box average ice condensate flux (kg m^-2 s^-1)
520 real(r8), intent(out) :: rflx(mgncol,2:nlev+1) ! grid-box average rain flux (kg m^-2 s^-1)
521 real(r8), intent(out) :: sflx(mgncol,2:nlev+1) ! grid-box average snow flux (kg m^-2 s^-1)
522 real(r8), intent(out) :: qrout(mgncol,nlev) ! grid-box average rain mixing ratio (kg/kg)
523 real(r8), intent(out) :: reff_rain(mgncol,nlev) ! rain effective radius (micron)
524 real(r8), intent(out) :: reff_snow(mgncol,nlev) ! snow effective radius (micron)
525 real(r8), intent(out) :: qcsevap(mgncol,nlev) ! cloud water evaporation due to sedimentation (1/s)
526 real(r8), intent(out) :: qisevap(mgncol,nlev) ! cloud ice sublimation due to sedimentation (1/s)
527 real(r8), intent(out) :: qvres(mgncol,nlev) ! residual condensation term to ensure RH < 100% (1/s)
528 real(r8), intent(out) :: cmeitot(mgncol,nlev) ! grid-mean cloud ice sub/dep (1/s)
529 real(r8), intent(out) :: vtrmc(mgncol,nlev) ! mass-weighted cloud water fallspeed (m/s)
530 real(r8), intent(out) :: vtrmi(mgncol,nlev) ! mass-weighted cloud ice fallspeed (m/s)
531 real(r8), intent(out) :: umr(mgncol,nlev) ! mass weighted rain fallspeed (m/s)
532 real(r8), intent(out) :: ums(mgncol,nlev) ! mass weighted snow fallspeed (m/s)
533 real(r8), intent(out) :: qcsedten(mgncol,nlev) ! qc sedimentation tendency (1/s)
534 real(r8), intent(out) :: qisedten(mgncol,nlev) ! qi sedimentation tendency (1/s)
535 real(r8), intent(out) :: qrsedten(mgncol,nlev) ! qr sedimentation tendency (1/s)
536 real(r8), intent(out) :: qssedten(mgncol,nlev) ! qs sedimentation tendency (1/s)
537
538 ! microphysical process rates for output (mixing ratio tendencies) (all have units of 1/s)
539 real(r8), intent(out) :: pratot(mgncol,nlev) ! accretion of cloud by rain
540 real(r8), intent(out) :: prctot(mgncol,nlev) ! autoconversion of cloud to rain
541 real(r8), intent(out) :: mnuccctot(mgncol,nlev) ! mixing ratio tend due to immersion freezing
542 real(r8), intent(out) :: mnuccttot(mgncol,nlev) ! mixing ratio tend due to contact freezing
543 real(r8), intent(out) :: msacwitot(mgncol,nlev) ! mixing ratio tend due to H-M splintering
544 real(r8), intent(out) :: psacwstot(mgncol,nlev) ! collection of cloud water by snow
545 real(r8), intent(out) :: bergstot(mgncol,nlev) ! bergeron process on snow
546 real(r8), intent(out) :: bergtot(mgncol,nlev) ! bergeron process on cloud ice
547 real(r8), intent(out) :: melttot(mgncol,nlev) ! melting of cloud ice
548 real(r8), intent(out) :: homotot(mgncol,nlev) ! homogeneous freezing cloud water
549 real(r8), intent(out) :: qcrestot(mgncol,nlev) ! residual cloud condensation due to removal of excess supersat
550 real(r8), intent(out) :: prcitot(mgncol,nlev) ! autoconversion of cloud ice to snow
551 real(r8), intent(out) :: praitot(mgncol,nlev) ! accretion of cloud ice by snow
552 real(r8), intent(out) :: qirestot(mgncol,nlev) ! residual ice deposition due to removal of excess supersat
553 real(r8), intent(out) :: mnuccrtot(mgncol,nlev) ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s)
554 real(r8), intent(out) :: pracstot(mgncol,nlev) ! mixing ratio tendency due to accretion of rain by snow (1/s)
555 real(r8), intent(out) :: meltsdttot(mgncol,nlev)! latent heating rate due to melting of snow (W/kg)
556 real(r8), intent(out) :: frzrdttot(mgncol,nlev) ! latent heating rate due to homogeneous freezing of rain (W/kg)
557 real(r8), intent(out) :: mnuccdtot(mgncol,nlev) ! mass tendency from ice nucleation
558 real(r8), intent(out) :: nrout(mgncol,nlev) ! rain number concentration (1/m3)
559 real(r8), intent(out) :: nsout(mgncol,nlev) ! snow number concentration (1/m3)
560 real(r8), intent(out) :: refl(mgncol,nlev) ! analytic radar reflectivity
561 real(r8), intent(out) :: arefl(mgncol,nlev) ! average reflectivity will zero points outside valid range
562 real(r8), intent(out) :: areflz(mgncol,nlev) ! average reflectivity in z.
563 real(r8), intent(out) :: frefl(mgncol,nlev) ! fractional occurrence of radar reflectivity
564 real(r8), intent(out) :: csrfl(mgncol,nlev) ! cloudsat reflectivity
565 real(r8), intent(out) :: acsrfl(mgncol,nlev) ! cloudsat average
566 real(r8), intent(out) :: fcsrfl(mgncol,nlev) ! cloudsat fractional occurrence of radar reflectivity
567 real(r8), intent(out) :: rercld(mgncol,nlev) ! effective radius calculation for rain + cloud
568 real(r8), intent(out) :: ncai(mgncol,nlev) ! output number conc of ice nuclei available (1/m3)
569 real(r8), intent(out) :: ncal(mgncol,nlev) ! output number conc of CCN (1/m3)
570 real(r8), intent(out) :: qrout2(mgncol,nlev) ! copy of qrout as used to compute drout2
571 real(r8), intent(out) :: qsout2(mgncol,nlev) ! copy of qsout as used to compute dsout2
572 real(r8), intent(out) :: nrout2(mgncol,nlev) ! copy of nrout as used to compute drout2
573 real(r8), intent(out) :: nsout2(mgncol,nlev) ! copy of nsout as used to compute dsout2
574 real(r8), intent(out) :: drout2(mgncol,nlev) ! mean rain particle diameter (m)
575 real(r8), intent(out) :: dsout2(mgncol,nlev) ! mean snow particle diameter (m)
576 real(r8), intent(out) :: freqs(mgncol,nlev) ! fractional occurrence of snow
577 real(r8), intent(out) :: freqr(mgncol,nlev) ! fractional occurrence of rain
578 real(r8), intent(out) :: nfice(mgncol,nlev) ! fractional occurrence of ice
579 real(r8), intent(out) :: qcrat(mgncol,nlev) ! limiter for qc process rates (1=no limit --> 0. no qc)
580
581 real(r8), intent(out) :: prer_evap(mgncol,nlev)
582
583
584 ! Tendencies calculated by external schemes that can replace MG's native
585 ! process tendencies.
586
587 ! Used with CARMA cirrus microphysics
588 ! (or similar external microphysics model)
589 ! real(r8), intent(in) :: tnd_qsnow(:,:) ! snow mass tendency (kg/kg/s)
590 ! real(r8), intent(in) :: tnd_nsnow(:,:) ! snow number tendency (#/kg/s)
591 ! real(r8), intent(in) :: re_ice(:,:) ! ice effective radius (m)
592
593 ! From external ice nucleation.
594 !real(r8), intent(in) :: frzimm(:,:) ! Number tendency due to immersion freezing (1/cm3)
595 !real(r8), intent(in) :: frzcnt(:,:) ! Number tendency due to contact freezing (1/cm3)
596 !real(r8), intent(in) :: frzdep(:,:) ! Number tendency due to deposition nucleation (1/cm3)
597
598 ! local workspace
599 ! all units mks unless otherwise stated
600
601 ! local copies of input variables
602 real(r8) :: qc(mgncol,nlev) ! cloud liquid mixing ratio (kg/kg)
603 real(r8) :: qi(mgncol,nlev) ! cloud ice mixing ratio (kg/kg)
604 real(r8) :: nc(mgncol,nlev) ! cloud liquid number concentration (1/kg)
605 real(r8) :: ni(mgncol,nlev) ! cloud liquid number concentration (1/kg)
606 real(r8) :: qr(mgncol,nlev) ! rain mixing ratio (kg/kg)
607 real(r8) :: qs(mgncol,nlev) ! snow mixing ratio (kg/kg)
608 real(r8) :: nr(mgncol,nlev) ! rain number concentration (1/kg)
609 real(r8) :: ns(mgncol,nlev) ! snow number concentration (1/kg)
610
611 ! general purpose variables
612 real(r8) :: deltat ! sub-time step (s)
613 real(r8) :: oneodt ! one / deltat
614 real(r8) :: mtime ! the assumed ice nucleation timescale
615
616 ! physical properties of the air at a given point
617 real(r8) :: rho(mgncol,nlev) ! density (kg m-3)
618 real(r8) :: rhoinv(mgncol,nlev) ! one / density (kg m-3)
619 real(r8) :: dv(mgncol,nlev) ! diffusivity of water vapor
620 real(r8) :: mu(mgncol,nlev) ! viscosity
621 real(r8) :: sc(mgncol,nlev) ! schmidt number
622 real(r8) :: rhof(mgncol,nlev) ! density correction factor for fallspeed
623
624 ! cloud fractions
625 real(r8) :: precip_frac(mgncol,nlev)! precip fraction assuming maximum overlap
626 real(r8) :: cldm(mgncol,nlev) ! cloud fraction
627 real(r8) :: icldm(mgncol,nlev) ! ice cloud fraction
628 real(r8) :: lcldm(mgncol,nlev) ! liq cloud fraction
629 real(r8) :: qsfm(mgncol,nlev) ! subgrid cloud water saturation scaling factor
630
631 ! mass mixing ratios
632 real(r8) :: qcic(mgncol,nlev) ! in-cloud cloud liquid
633 real(r8) :: qiic(mgncol,nlev) ! in-cloud cloud ice
634 real(r8) :: qsic(mgncol,nlev) ! in-precip snow
635 real(r8) :: qric(mgncol,nlev) ! in-precip rain
636
637 ! number concentrations
638 real(r8) :: ncic(mgncol,nlev) ! in-cloud droplet
639 real(r8) :: niic(mgncol,nlev) ! in-cloud cloud ice
640 real(r8) :: nsic(mgncol,nlev) ! in-precip snow
641 real(r8) :: nric(mgncol,nlev) ! in-precip rain
642 ! maximum allowed ni value
643 real(r8) :: nimax(mgncol,nlev)
644
645 ! Size distribution parameters for:
646 ! cloud ice
647 real(r8) :: lami(mgncol,nlev) ! slope
648 real(r8) :: n0i(mgncol,nlev) ! intercept
649 ! cloud liquid
650 real(r8) :: lamc(mgncol,nlev) ! slope
651 real(r8) :: pgam(mgncol,nlev) ! spectral width parameter
652 ! snow
653 real(r8) :: lams(mgncol,nlev) ! slope
654 real(r8) :: n0s(mgncol,nlev) ! intercept
655 ! rain
656 real(r8) :: lamr(mgncol,nlev) ! slope
657 real(r8) :: n0r(mgncol,nlev) ! intercept
658
659 ! Rates/tendencies due to:
660
661 ! Instantaneous snow melting
662 real(r8) :: minstsm(mgncol,nlev) ! mass mixing ratio
663 real(r8) :: ninstsm(mgncol,nlev) ! number concentration
664 ! Instantaneous rain freezing
665 real(r8) :: minstrf(mgncol,nlev) ! mass mixing ratio
666 real(r8) :: ninstrf(mgncol,nlev) ! number concentration
667
668 ! deposition of cloud ice
669 real(r8) :: vap_dep(mgncol,nlev) ! deposition from vapor to ice PMC 12/3/12
670 ! sublimation of cloud ice
671 real(r8) :: ice_sublim(mgncol,nlev) ! sublimation from ice to vapor PMC 12/3/12
672 ! ice nucleation
673 real(r8) :: nnuccd(mgncol,nlev) ! number rate from deposition/cond.-freezing
674 real(r8) :: mnuccd(mgncol,nlev) ! mass mixing ratio
675 ! freezing of cloud water
676 real(r8) :: mnuccc(mgncol,nlev) ! mass mixing ratio
677 real(r8) :: nnuccc(mgncol,nlev) ! number concentration
678 ! contact freezing of cloud water
679 real(r8) :: mnucct(mgncol,nlev) ! mass mixing ratio
680 real(r8) :: nnucct(mgncol,nlev) ! number concentration
681 ! deposition nucleation in mixed-phase clouds (from external scheme)
682 real(r8) :: mnudep(mgncol,nlev) ! mass mixing ratio
683 real(r8) :: nnudep(mgncol,nlev) ! number concentration
684 ! ice multiplication
685 real(r8) :: msacwi(mgncol,nlev) ! mass mixing ratio
686 real(r8) :: nsacwi(mgncol,nlev) ! number concentration
687 ! autoconversion of cloud droplets
688 real(r8) :: prc(mgncol,nlev) ! mass mixing ratio
689 real(r8) :: nprc(mgncol,nlev) ! number concentration (rain)
690 real(r8) :: nprc1(mgncol,nlev) ! number concentration (cloud droplets)
691 ! self-aggregation of snow
692 real(r8) :: nsagg(mgncol,nlev) ! number concentration
693 ! self-collection of rain
694 real(r8) :: nragg(mgncol,nlev) ! number concentration
695 ! collection of droplets by snow
696 real(r8) :: psacws(mgncol,nlev) ! mass mixing ratio
697 real(r8) :: npsacws(mgncol,nlev) ! number concentration
698 ! collection of rain by snow
699 real(r8) :: pracs(mgncol,nlev) ! mass mixing ratio
700 real(r8) :: npracs(mgncol,nlev) ! number concentration
701 ! freezing of rain
702 real(r8) :: mnuccr(mgncol,nlev) ! mass mixing ratio
703 real(r8) :: nnuccr(mgncol,nlev) ! number concentration
704 ! freezing of rain to form ice (mg add 4/26/13)
705 real(r8) :: mnuccri(mgncol,nlev) ! mass mixing ratio
706 real(r8) :: nnuccri(mgncol,nlev) ! number concentration
707 ! accretion of droplets by rain
708 real(r8) :: pra(mgncol,nlev) ! mass mixing ratio
709 real(r8) :: npra(mgncol,nlev) ! number concentration
710 ! autoconversion of cloud ice to snow
711 real(r8) :: prci(mgncol,nlev) ! mass mixing ratio
712 real(r8) :: nprci(mgncol,nlev) ! number concentration
713 ! accretion of cloud ice by snow
714 real(r8) :: prai(mgncol,nlev) ! mass mixing ratio
715 real(r8) :: nprai(mgncol,nlev) ! number concentration
716 ! evaporation of rain
717 real(r8) :: pre(mgncol,nlev) ! mass mixing ratio
718 ! sublimation of snow
719 real(r8) :: prds(mgncol,nlev) ! mass mixing ratio
720 ! number evaporation
721 real(r8) :: nsubi(mgncol,nlev) ! cloud ice
722 real(r8) :: nsubc(mgncol,nlev) ! droplet
723 real(r8) :: nsubs(mgncol,nlev) ! snow
724 real(r8) :: nsubr(mgncol,nlev) ! rain
725 ! bergeron process
726 real(r8) :: berg(mgncol,nlev) ! mass mixing ratio (cloud ice)
727 real(r8) :: bergs(mgncol,nlev) ! mass mixing ratio (snow)
728
729
730 ! fallspeeds
731 ! number-weighted
732 real(r8) :: uns(mgncol,nlev) ! snow
733 real(r8) :: unr(mgncol,nlev) ! rain
734 ! air density corrected fallspeed parameters
735 real(r8) :: arn(mgncol,nlev) ! rain
736 real(r8) :: asn(mgncol,nlev) ! snow
737 real(r8) :: acn(mgncol,nlev) ! cloud droplet
738 real(r8) :: ain(mgncol,nlev) ! cloud ice
739 real(r8) :: ajn(mgncol,nlev) ! cloud small ice
740
741 ! Mass of liquid droplets used with external heterogeneous freezing.
742 real(r8) :: mi0l(mgncol)
743
744 ! saturation vapor pressures
745 real(r8) :: esl(mgncol,nlev) ! liquid
746 real(r8) :: esi(mgncol,nlev) ! ice
747 real(r8) :: esn ! checking for RH after rain evap
748
749 ! saturation vapor mixing ratios
750 real(r8) :: qvl(mgncol,nlev) ! liquid
751 real(r8) :: qvi(mgncol,nlev) ! ice
752 real(r8) :: qvn ! checking for RH after rain evap
753
754 ! relative humidity
755 real(r8) :: relhum(mgncol,nlev)
756
757 ! parameters for cloud water and cloud ice sedimentation calculations
758 real(r8) :: fc(mgncol,nlev)
759 real(r8) :: fnc(mgncol,nlev)
760 real(r8) :: fi(mgncol,nlev)
761 real(r8) :: fni(mgncol,nlev)
762
763 real(r8) :: fr(mgncol,nlev)
764 real(r8) :: fnr(mgncol,nlev)
765 real(r8) :: fs(mgncol,nlev)
766 real(r8) :: fns(mgncol,nlev)
767
768 real(r8) :: faloutc(nlev)
769 real(r8) :: faloutnc(nlev)
770 real(r8) :: falouti(nlev)
771 real(r8) :: faloutni(nlev)
772
773 real(r8) :: faloutr(nlev)
774 real(r8) :: faloutnr(nlev)
775 real(r8) :: falouts(nlev)
776 real(r8) :: faloutns(nlev)
777
778 real(r8) :: faltndc
779 real(r8) :: faltndnc
780 real(r8) :: faltndi
781 real(r8) :: faltndni
782 real(r8) :: faltndqie
783 real(r8) :: faltndqce
784
785 real(r8) :: faltndr
786 real(r8) :: faltndnr
787 real(r8) :: faltnds
788 real(r8) :: faltndns
789
790 real(r8) :: rainrt(mgncol,nlev) ! rain rate for reflectivity calculation
791
792 ! dummy variables
793 real(r8) :: dum
794 real(r8) :: dum1
795 real(r8) :: dum2
796 real(r8) :: dum3
797 real(r8) :: dumni0
798 real(r8) :: dumns0
799 real(r8) :: tx1, tx2, tx3, tx4, tx5, tx6, tx7, grho
800 ! dummies for checking RH
801 real(r8) :: qtmp
802 real(r8) :: ttmp
803 ! dummies for conservation check
804 real(r8) :: ratio
805 real(r8) :: tmpfrz
806 ! dummies for in-cloud variables
807 real(r8) :: dumc(mgncol,nlev) ! qc
808 real(r8) :: dumnc(mgncol,nlev) ! nc
809 real(r8) :: dumi(mgncol,nlev) ! qi
810 real(r8) :: dumni(mgncol,nlev) ! ni
811 real(r8) :: dumr(mgncol,nlev) ! rain mixing ratio
812 real(r8) :: dumnr(mgncol,nlev) ! rain number concentration
813 real(r8) :: dums(mgncol,nlev) ! snow mixing ratio
814 real(r8) :: dumns(mgncol,nlev) ! snow number concentration
815 ! Array dummy variable
816 !real(r8) :: dum_2D(mgncol,nlev)
817 real(r8) :: pdel_inv(mgncol,nlev)
818 real(r8) :: ts_au_loc(mgncol)
819
820 ! loop array variables
821 ! "i" and "k" are column/level iterators for internal (MG) variables
822 ! "n" is used for other looping (currently just sedimentation)
823 integer i, k, n
824
825 ! number of sub-steps for loops over "n" (for sedimentation)
826 integer nstep, mdust, nlb, nstep_def
827
828 ! Varaibles to scale fall velocity between small and regular ice regimes.
829 real(r8) :: irad, ifrac, tsfac
830! logical, parameter :: do_ice_gmao=.false., do_liq_liu=.false.
831! logical, parameter :: do_ice_gmao=.true., do_liq_liu=.true.
832! logical, parameter :: do_ice_gmao=.true., do_liq_liu=.false.
833 real(r8), parameter :: qimax=0.010, qimin=0.005, qiinv=one/(qimax-qimin)
834! ts_au_min=180.0
835
836 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
837
838
839 ! Process inputs
840
841 ! assign variable deltat to deltatin
842 deltat = deltatin
843 oneodt = one / deltat
844! nstep_def = max(1, nint(deltat/20))
845 nstep_def = max(1, nint(deltat/5))
846! tsfac = log(ts_au/ts_au_min) * qiinv
847
848 ! Copies of input concentrations that may be changed internally.
849 do k=1,nlev
850 do i=1,mgncol
851 qc(i,k) = qcn(i,k)
852 nc(i,k) = ncn(i,k)
853 qi(i,k) = qin(i,k)
854 ni(i,k) = nin(i,k)
855 qr(i,k) = qrn(i,k)
856 nr(i,k) = nrn(i,k)
857 qs(i,k) = qsn(i,k)
858 ns(i,k) = nsn(i,k)
859 enddo
860 enddo
861
862 ! cldn: used to set cldm, unused for subcolumns
863 ! liqcldf: used to set lcldm, unused for subcolumns
864 ! icecldf: used to set icldm, unused for subcolumns
865
866 if (microp_uniform) then
867 ! subcolumns, set cloud fraction variables to one
868 ! if cloud water or ice is present, if not present
869 ! set to mincld (mincld used instead of zero, to prevent
870 ! possible division by zero errors).
871
872 do k=1,nlev
873 do i=1,mgncol
874
875 if (qc(i,k) >= qsmall) then
876 lcldm(i,k) = one
877 else
878 lcldm(i,k) = mincld
879 endif
880
881 if (qi(i,k) >= qsmall) then
882 icldm(i,k) = one
883 else
884 icldm(i,k) = mincld
885 endif
886
887 cldm(i,k) = max(icldm(i,k), lcldm(i,k))
888! qsfm(i,k) = one
889 qsfm(i,k) = qsatfac(i,k)
890 enddo
891 enddo
892
893 else ! get cloud fraction, check for minimum
894 do k=1,nlev
895 do i=1,mgncol
896 cldm(i,k) = max(cldn(i,k), mincld)
897 lcldm(i,k) = max(liqcldf(i,k), mincld)
898 icldm(i,k) = max(icecldf(i,k), mincld)
899 qsfm(i,k) = qsatfac(i,k)
900 enddo
901 enddo
902 end if
903
904! if (lprnt) write(0,*)' cldm=',cldm(1,nlev-20:nlev)
905! if (lprnt) write(0,*)' liqcldf=',liqcldf(1,nlev-20:nlev)
906! if (lprnt) write(0,*)' lcldm=',lcldm(1,nlev-20:nlev)
907! if (lprnt) write(0,*)' icecldf=',icecldf(1,nlev-20:nlev)
908! if (lprnt) write(0,*)' icldm=',icldm(1,nlev-20:nlev)
909! if (lprnt) write(0,*)' qsfm=',qsfm(1,nlev-20:nlev)
910
911 ! Initialize local variables
912
913 ! local physical properties
914 do k=1,nlev
915 do i=1,mgncol
916! rho(i,k) = p(i,k) / (r*t(i,k)*(one+fv*q(i,k)))
917 rho(i,k) = p(i,k) / (r*t(i,k))
918 rhoinv(i,k) = one / rho(i,k)
919 dv(i,k) = 8.794e-5_r8 * t(i,k)**1.81_r8 / p(i,k)
920 mu(i,k) = 1.496e-6_r8 * t(i,k)*sqrt(t(i,k)) / (t(i,k) + 120._r8)
921 sc(i,k) = mu(i,k) / (rho(i,k)*dv(i,k))
922
923 ! air density adjustment for fallspeed parameters
924 ! includes air density correction factor to the
925 ! power of 0.54 following Heymsfield and Bansemer 2007
926
927 rhof(i,k) = (rhosu*rhoinv(i,k))**0.54_r8
928
929 arn(i,k) = ar*rhof(i,k)
930 asn(i,k) = as*rhof(i,k)
931 acn(i,k) = g*rhow/(18._r8*mu(i,k))
932 tx1 = (rhosu*rhoinv(i,k))**0.35_r8
933 ain(i,k) = ai*tx1
934 ajn(i,k) = aj*tx1
935
936 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
937 ! Get humidity and saturation vapor pressures
938
939! do k=1,nlev
940! do i=1,mgncol
941! relvar(i,k) = relvar_i
942 accre_enhan(i,k) = accre_enhan_i
943! call qsat_water(t(i,k), p(i,k), esl(i,k), qvl(i,k))
944 esl(i,k) = min(fpvsl(t(i,k)), p(i,k))
945 qvl(i,k) = epsqs*esl(i,k) / (p(i,k)-omeps*esl(i,k))
946
947
948 ! make sure when above freezing that esi=esl, not active yet
949 if (t(i,k) >= tmelt) then
950 esi(i,k) = esl(i,k)
951 qvi(i,k) = qvl(i,k)
952 else
953! call qsat_ice(t(i,k), p(i,k), esi(i,k), qvi(i,k))
954 esi(i,k) = min(fpvsi(t(i,k)), p(i,k))
955 qvi(i,k) = epsqs*esi(i,k) / (p(i,k)-omeps*esi(i,k))
956 end if
957
958 ! Scale the water saturation values to reflect subgrid scale
959 ! ice cloud fraction, where ice clouds begin forming at a
960 ! gridbox average relative humidity of rhmini (not 1).
961 !
962 ! NOTE: For subcolumns and other non-subgrid clouds, qsfm will be 1.
963 qvi(i,k) = qsfm(i,k) * qvi(i,k)
964! esi(i,k) = qsfm(i,k) * esi(i,k)
965 qvl(i,k) = qsfm(i,k) * qvl(i,k)
966! esl(i,k) = qsfm(i,k) * esl(i,k)
967
968 relhum(i,k) = max(zero, min(q(i,k)/max(qvl(i,k), qsmall), two))
969 end do
970 end do
971
972
973 !===============================================
974
975 ! set mtime here to avoid answer-changing
976 mtime = deltat
977
978 ! initialize microphysics output
979 do k=1,nlev
980 do i=1,mgncol
981 qcsevap(i,k) = zero
982 qisevap(i,k) = zero
983 qvres(i,k) = zero
984 cmeitot(i,k) = zero
985 vtrmc(i,k) = zero
986 vtrmi(i,k) = zero
987 qcsedten(i,k) = zero
988 qisedten(i,k) = zero
989 qrsedten(i,k) = zero
990 qssedten(i,k) = zero
991
992 pratot(i,k) = zero
993 prctot(i,k) = zero
994 mnuccctot(i,k) = zero
995 mnuccttot(i,k) = zero
996 msacwitot(i,k) = zero
997 psacwstot(i,k) = zero
998 bergstot(i,k) = zero
999 bergtot(i,k) = zero
1000 melttot(i,k) = zero
1001 homotot(i,k) = zero
1002 qcrestot(i,k) = zero
1003 prcitot(i,k) = zero
1004 praitot(i,k) = zero
1005 qirestot(i,k) = zero
1006 mnuccrtot(i,k) = zero
1007 pracstot(i,k) = zero
1008 meltsdttot(i,k) = zero
1009 frzrdttot(i,k) = zero
1010 mnuccdtot(i,k) = zero
1011
1012 rflx(i,k+1) = zero
1013 sflx(i,k+1) = zero
1014 lflx(i,k+1) = zero
1015 iflx(i,k+1) = zero
1016
1017 ! initialize precip output
1018
1019 qrout(i,k) = zero
1020 qsout(i,k) = zero
1021 nrout(i,k) = zero
1022 nsout(i,k) = zero
1023
1024 ! for refl calc
1025 rainrt(i,k) = zero
1026
1027 ! initialize rain size
1028 rercld(i,k) = zero
1029
1030 qcsinksum_rate1ord(i,k) = zero
1031
1032 ! initialize variables for trop_mozart
1033 nevapr(i,k) = zero
1034 prer_evap(i,k) = zero
1035 evapsnow(i,k) = zero
1036 am_evp_st(i,k) = zero
1037 prain(i,k) = zero
1038 prodsnow(i,k) = zero
1039 cmeout(i,k) = zero
1040
1041 precip_frac(i,k) = mincld
1042
1043 lamc(i,k) = zero
1044
1045 ! initialize microphysical tendencies
1046
1047 tlat(i,k) = zero
1048 qvlat(i,k) = zero
1049 qctend(i,k) = zero
1050 qitend(i,k) = zero
1051 qstend(i,k) = zero
1052 qrtend(i,k) = zero
1053 nctend(i,k) = zero
1054 nitend(i,k) = zero
1055 nrtend(i,k) = zero
1056 nstend(i,k) = zero
1057
1058 ! initialize in-cloud and in-precip quantities to zero
1059 qcic(i,k) = zero
1060 qiic(i,k) = zero
1061 qsic(i,k) = zero
1062 qric(i,k) = zero
1063
1064 ncic(i,k) = zero
1065 niic(i,k) = zero
1066 nsic(i,k) = zero
1067 nric(i,k) = zero
1068
1069 ! initialize precip fallspeeds to zero
1070 ums(i,k) = zero
1071 uns(i,k) = zero
1072 umr(i,k) = zero
1073 unr(i,k) = zero
1074
1075 ! initialize limiter for output
1076 qcrat(i,k) = one
1077
1078 ! Many outputs have to be initialized here at the top to work around
1079 ! ifort problems, even if they are always overwritten later.
1080 effc(i,k) = ten
1081 lamcrad(i,k) = zero
1082 pgamrad(i,k) = zero
1083 effc_fn(i,k) = ten
1084 effi(i,k) = 25._r8
1085 sadice(i,k) = zero
1086 sadsnow(i,k) = zero
1087 deffi(i,k) = 50._r8
1088
1089 qrout2(i,k) = zero
1090 nrout2(i,k) = zero
1091 drout2(i,k) = zero
1092 qsout2(i,k) = zero
1093 nsout2(i,k) = zero
1094 dsout(i,k) = zero
1095 dsout2(i,k) = zero
1096
1097 freqr(i,k) = zero
1098 freqs(i,k) = zero
1099
1100 reff_rain(i,k) = zero
1101 reff_snow(i,k) = zero
1102
1103 refl(i,k) = -9999._r8
1104 arefl(i,k) = zero
1105 areflz(i,k) = zero
1106 frefl(i,k) = zero
1107 csrfl(i,k) = zero
1108 acsrfl(i,k) = zero
1109 fcsrfl(i,k) = zero
1110
1111 ncal(i,k) = zero
1112 ncai(i,k) = zero
1113
1114 nfice(i,k) = zero
1115 npccn(i,k) = zero
1116 enddo
1117 enddo
1118
1119 if(iccn == 1) then
1120 do k=1,nlev
1121 do i=1,mgncol
1122 npccn(i,k) = npccnin(i,k)
1123 enddo
1124 enddo
1125 else
1126 do k=1,nlev
1127 do i=1,mgncol
1128 npccn(i,k) = max((npccnin(i,k)*lcldm(i,k)-nc(i,k))*oneodt, zero)
1129 enddo
1130 enddo
1131 endif
1132 ! initialize precip at surface
1133
1134 do i=1,mgncol
1135 prect(i) = zero
1136 preci(i) = zero
1137 enddo
1138
1139 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1140 ! droplet activation
1141 ! get provisional droplet number after activation. This is used for
1142 ! all microphysical process calculations, for consistency with update of
1143 ! droplet mass before microphysics
1144
1145 ! calculate potential for droplet activation if cloud water is present
1146 ! tendency from activation (npccn) is read in from companion routine
1147
1148 ! output activated liquid and ice (convert from #/kg -> #/m3)
1149 !--------------------------------------------------
1150 where (qc >= qsmall)
1151 npccn = max((npccnin*lcldm-nc)*oneodt, zero)
1152 nc = max(nc + npccn*deltat, zero)
1153 ncal = nc*rho/lcldm ! sghan minimum in #/cm3
1154 elsewhere
1155 ncal = zero
1156 end where
1157
1158 if (iccn == 1) then
1159 do k=1,nlev
1160 do i=1,mgncol
1161 if (t(i,k) < icenuct) then
1162 ncai(i,k) = 0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k))) * 1000._r8
1163! ncai(i,k) = min(ncai(i,k), 208.9e3_r8)
1164 ncai(i,k) = min(ncai(i,k), 355.0e3_r8)
1165 naai(i,k) = (ncai(i,k)*rhoinv(i,k) + naai(i,k)) * half
1166 ncai(i,k) = naai(i,k)*rho(i,k)
1167 else
1168 naai(i,k) = zero
1169 ncai(i,k) = zero
1170 endif
1171 enddo
1172 enddo
1173 elseif (iccn == 2) then
1174 do k=1,nlev
1175 do i=1,mgncol
1176 if (t(i,k) < icenuct) then
1177 ncai(i,k) = naai(i,k)*rho(i,k)
1178 else
1179 naai(i,k) = zero
1180 ncai(i,k) = zero
1181 endif
1182 enddo
1183 enddo
1184 else
1185 do k=1,nlev
1186 do i=1,mgncol
1187 if (t(i,k) < icenuct) then
1188 ncai(i,k) = 0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k))) * 1000._r8
1189 ncai(i,k) = min(ncai(i,k), 355.0e3_r8)
1190 naai(i,k) = ncai(i,k)*rhoinv(i,k)
1191 else
1192 naai(i,k) = zero
1193 ncai(i,k) = zero
1194 endif
1195 enddo
1196 enddo
1197 do k=1,nlev
1198 do i=1,mgncol
1199 naai(i,k) = zero
1200 ncai(i,k) = zero
1201 enddo
1202 enddo
1203 endif
1204
1205
1206 !===============================================
1207
1208 ! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5%
1209 !
1210 ! NOTE: If using gridbox average values, condensation will not occur until rh=1,
1211 ! so the threshold seems like it should be 1.05 and not rhmini + 0.05. For subgrid
1212 ! clouds (using rhmini and qsfacm), the relhum has already been adjusted, and thus
1213 ! the nucleation threshold should also be 1.05 and not rhmini + 0.05.
1214
1215 !-------------------------------------------------------
1216
1217 if (do_cldice) then
1218 where (naai > zero .and. t < icenuct .and. relhum*esl/esi > 1.05_r8)
1219
1220 !if NAAI > 0. then set numice = naai (as before)
1221 !note: this is gridbox averaged
1222 nnuccd = (naai-ni/icldm)/mtime*icldm
1223 nnuccd = max(nnuccd, zero)
1224 nimax = naai*icldm
1225
1226 !Calc mass of new particles using new crystal mass...
1227 !also this will be multiplied by mtime as nnuccd is...
1228
1229 mnuccd = nnuccd * mi0
1230
1231 elsewhere
1232 nnuccd = zero
1233 nimax = zero
1234 mnuccd = zero
1235 end where
1236
1237 end if
1238
1239
1240 !=============================================================================
1241 do k=1,nlev
1242
1243 do i=1,mgncol
1244
1245 ! calculate instantaneous precip processes (melting and homogeneous freezing)
1246
1247 ! melting of snow at +2 C
1248
1249 if (t(i,k) > snowmelt) then
1250 if (qs(i,k) > zero) then
1251
1252 ! make sure melting snow doesn't reduce temperature below threshold
1253 dum = -(xlf/cpp) * qs(i,k)
1254 if (t(i,k)+dum < snowmelt) then
1255 dum = min(one, max(zero, (cpp/xlf)*(t(i,k)-snowmelt)/qs(i,k)))
1256 else
1257 dum = one
1258 end if
1259
1260 minstsm(i,k) = dum*qs(i,k)
1261 ninstsm(i,k) = dum*ns(i,k)
1262
1263 dum1 = - minstsm(i,k) * (xlf*oneodt)
1264 tlat(i,k) = tlat(i,k) + dum1
1265 meltsdttot(i,k) = meltsdttot(i,k) + dum1
1266
1267 qs(i,k) = max(qs(i,k) - minstsm(i,k), zero)
1268 ns(i,k) = max(ns(i,k) - ninstsm(i,k), zero)
1269 qr(i,k) = max(qr(i,k) + minstsm(i,k), zero)
1270 nr(i,k) = max(nr(i,k) + ninstsm(i,k), zero)
1271 end if
1272 end if
1273
1274 end do
1275 end do
1276! if (lprnt) write(0,*)' tlat1=',tlat(1,:)*deltat
1277
1278 do k=1,nlev
1279 do i=1,mgncol
1280 ! freezing of rain at -5 C
1281
1282 if (t(i,k) < rainfrze) then
1283
1284 if (qr(i,k) > zero) then
1285
1286 ! make sure freezing rain doesn't increase temperature above threshold
1287 dum = (xlf/cpp) * qr(i,k)
1288 if (t(i,k)+dum > rainfrze) then
1289 dum = -(t(i,k)-rainfrze) * (cpp/xlf)
1290 dum = min(one, max(zero, dum/qr(i,k)))
1291 else
1292 dum = one
1293 end if
1294
1295 minstrf(i,k) = dum*qr(i,k)
1296 ninstrf(i,k) = dum*nr(i,k)
1297
1298 ! heating tendency
1299 dum1 = minstrf(i,k) * (xlf*oneodt)
1300 tlat(i,k) = tlat(i,k) + dum1
1301 frzrdttot(i,k) = frzrdttot(i,k) + dum1
1302
1303 qr(i,k) = max(qr(i,k) - minstrf(i,k), zero)
1304 nr(i,k) = max(nr(i,k) - ninstrf(i,k), zero)
1305 qs(i,k) = max(qs(i,k) + minstrf(i,k), zero)
1306 ns(i,k) = max(ns(i,k) + ninstrf(i,k), zero)
1307
1308 end if
1309 end if
1310 end do
1311 end do
1312
1313! if (lprnt) write(0,*)' tlat2=',tlat(1,:)*deltat
1314 do k=1,nlev
1315 do i=1,mgncol
1316 ! obtain in-cloud values of cloud water/ice mixing ratios and number concentrations
1317 !-------------------------------------------------------
1318 ! for microphysical process calculations
1319 ! units are kg/kg for mixing ratio, 1/kg for number conc
1320
1321 if (qc(i,k) >= qsmall) then
1322 dum = one / lcldm(i,k)
1323! qcic(i,k) = min(qc(i,k)*dum, 5.e-3_r8) ! limit in-cloud values to 0.005 kg/kg
1324 qcic(i,k) = min(qc(i,k)*dum, 0.05_r8) ! limit in-cloud values to 0.05 kg/kg
1325 ncic(i,k) = max(nc(i,k)*dum, zero)
1326
1327 ! specify droplet concentration
1328 if (nccons) then
1329 ncic(i,k) = ncnst * rhoinv(i,k)
1330 end if
1331 else
1332 qcic(i,k) = zero
1333 ncic(i,k) = zero
1334 end if
1335
1336 if (qi(i,k) >= qsmall) then
1337 dum = one / icldm(i,k)
1338! qiic(i,k) = min(qi(i,k)*dum, 5.e-3_r8) ! limit in-cloud values to 0.005 kg/kg
1339 qiic(i,k) = min(qi(i,k)*dum, 0.05_r8) ! limit in-cloud values to 0.05 kg/kg
1340 niic(i,k) = max(ni(i,k)*dum, zero)
1341
1342 ! switch for specification of cloud ice number
1343 if (nicons) then
1344 niic(i,k) = ninst * rhoinv(i,k)
1345 end if
1346 else
1347 qiic(i,k) = zero
1348 niic(i,k) = zero
1349 end if
1350
1351 end do
1352 end do
1353
1354 !========================================================================
1355
1356 ! for sub-columns cldm has already been set to 1 if cloud
1357 ! water or ice is present, so precip_frac will be correctly set below
1358 ! and nothing extra needs to be done here
1359
1360 precip_frac = cldm
1361
1362 micro_vert_loop: do k=1,nlev
1363
1364 if (trim(micro_mg_precip_frac_method) == 'in_cloud') then
1365
1366 if (k /= 1) then
1367 where (qc(:,k) < qsmall .and. qi(:,k) < qsmall)
1368 precip_frac(:,k) = precip_frac(:,k-1)
1369 end where
1370 endif
1371
1372 else if (trim(micro_mg_precip_frac_method) == 'max_overlap') then
1373
1374 ! calculate precip fraction based on maximum overlap assumption
1375
1376 ! if rain or snow mix ratios are smaller than threshold,
1377 ! then leave precip_frac as cloud fraction at current level
1378 if (k /= 1) then
1379 where (qr(:,k-1) >= qsmall .or. qs(:,k-1) >= qsmall)
1380 precip_frac(:,k) = max(precip_frac(:,k-1),precip_frac(:,k))
1381 end where
1382 end if
1383
1384 endif
1385
1386
1387 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1388 ! get size distribution parameters based on in-cloud cloud water
1389 ! these calculations also ensure consistency between number and mixing ratio
1390 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1391
1392 ! cloud liquid
1393 !-------------------------------------------
1394
1395 call size_dist_param_liq(mg_liq_props, qcic(:,k), ncic(:,k), rho(:,k), &
1396 pgam(:,k), lamc(:,k), mgncol)
1397
1398
1399 !========================================================================
1400 ! autoconversion of cloud liquid water to rain
1401 ! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
1402 ! minimum qc of 1 x 10^-8 prevents floating point error
1403
1404 if (.not. do_sb_physics) then
1405 call kk2000_liq_autoconversion(microp_uniform, qcic(:,k), &
1406 ncic(:,k), rho(:,k), relvar(:,k), prc(:,k), nprc(:,k), nprc1(:,k), mgncol)
1407 endif
1408
1409 ! assign qric based on prognostic qr, using assumed precip fraction
1410 ! note: this could be moved above for consistency with qcic and qiic calculations
1411 do i=1,mgncol
1412 if (precip_frac(i,k) > mincld) then
1413 dum = one / precip_frac(i,k)
1414 else
1415 dum = zero
1416 endif
1417! qric(i,k) = min(qr(i,k)*dum, 0.01_r8) ! limit in-precip mixing ratios to 10 g/kg
1418 qric(i,k) = min(qr(i,k)*dum, 0.05_r8) ! limit in-precip mixing ratios to 50 g/kg
1419 nric(i,k) = nr(i,k) * dum
1420
1421
1422 ! add autoconversion to precip from above to get provisional rain mixing ratio
1423 ! and number concentration (qric and nric)
1424
1425 if(qric(i,k) < qsmall) then
1426 qric(i,k) = zero
1427 nric(i,k) = zero
1428 endif
1429
1430 ! make sure number concentration is a positive number to avoid
1431 ! taking root of negative later
1432
1433 nric(i,k) = max(nric(i,k),zero)
1434 enddo
1435 ! Get size distribution parameters for cloud ice
1436
1437 call size_dist_param_ice(mg_ice_props, qiic(:,k), niic(:,k), &
1438 lami(:,k), mgncol, n0=n0i(:,k))
1439
1440! call size_dist_param_basic(mg_ice_props, qiic(:,k), niic(:,k), &
1441! lami(:,k), mgncol, n0=n0i(:,k))
1442
1443 ! Alternative autoconversion
1444 if (do_sb_physics) then
1445 if (do_liq_liu) then
1446 call liu_liq_autoconversion(pgam(:,k),qcic(:,k),ncic(:,k), &
1447 qric(:,k),rho(:,k),relvar(:,k),prc(:,k),nprc(:,k),nprc1(:,k),mgncol)
1448 else
1449 call sb2001v2_liq_autoconversion(pgam(:,k),qcic(:,k),ncic(:,k), &
1450 qric(:,k),rho(:,k),relvar(:,k),prc(:,k),nprc(:,k),nprc1(:,k), mgncol)
1451 endif
1452 endif
1453
1454 !.......................................................................
1455 ! Autoconversion of cloud ice to snow
1456 ! similar to Ferrier (1994)
1457
1458 if (do_cldice) then
1459 do i=1,mgncol
1460 if (qiic(i,k) >= qimax) then
1461! if (qi(i,k) >= qimax) then
1462 ts_au_loc(i) = ts_au_min
1463 elseif (qiic(i,k) <= qimin) then
1464! elseif (qi(i,k) <= qimin) then
1465 ts_au_loc(i) = ts_au
1466 else
1467! ts_au_loc(i) = (ts_au*(qimax-qi(i,k)) + ts_au_min*(qi(i,k)-qimin)) * qiinv
1468 ts_au_loc(i) = (ts_au*(qimax-qiic(i,k)) + ts_au_min*(qiic(i,k)-qimin)) * qiinv
1469! ts_au_loc(i) = ts_au * exp(-tsfac*(qiic(i,k)-qimin))
1470 endif
1471 enddo
1472
1473 if(do_ice_gmao) then
1474 call gmao_ice_autoconversion(t(:,k), qiic(:,k), niic(:,k), lami(:,k), &
1475 n0i(:,k), dcs, ts_au_loc(:), prci(:,k), nprci(:,k), mgncol)
1476 else
1477 call ice_autoconversion(t(:,k), qiic(:,k), lami(:,k), n0i(:,k), &
1478 dcs, ts_au_loc(:), prci(:,k), nprci(:,k), mgncol)
1479 end if
1480 !else
1481 ! Add in the particles that we have already converted to snow, and
1482 ! don't do any further autoconversion of ice.
1483 !prci(:,k) = tnd_qsnow(:,k) / cldm(:,k)
1484 !nprci(:,k) = tnd_nsnow(:,k) / cldm(:,k)
1485 end if
1486
1487 ! note, currently we don't have this
1488 ! inside the do_cldice block, should be changed later
1489 ! assign qsic based on prognostic qs, using assumed precip fraction
1490 do i=1,mgncol
1491 if (precip_frac(i,k) > mincld) then
1492 dum = one / precip_frac(i,k)
1493 else
1494 dum = zero
1495 endif
1496! qsic(i,k) = min(qs(i,k)*dum, 0.01_r8) ! limit in-precip mixing ratios to 10 g/kg
1497! qsic(i,k) = min(qs(i,k)*dum, 0.05_r8) ! limit in-precip mixing ratios to 50 g/kg
1498 qsic(i,k) = min(qs(i,k)*dum, 0.10_r8) ! limit in-precip mixing ratios to 50 g/kg
1499 nsic(i,k) = ns(i,k) * dum
1500
1501 ! if precip mix ratio is zero so should number concentration
1502
1503 if(qsic(i,k) < qsmall) then
1504 qsic(i,k) = zero
1505 nsic(i,k) = zero
1506 endif
1507
1508 ! make sure number concentration is a positive number to avoid
1509 ! taking root of negative later
1510
1511 nsic(i,k) = max(nsic(i,k), zero)
1512 enddo
1513
1514 !.......................................................................
1515 ! get size distribution parameters for precip
1516 !......................................................................
1517 ! rain
1518
1519 call size_dist_param_basic(mg_rain_props, qric(:,k), nric(:,k), &
1520 lamr(:,k), mgncol, n0=n0r(:,k))
1521
1522 do i=1,mgncol
1523 if (lamr(i,k) >= qsmall) then
1524 dum = arn(i,k) / lamr(i,k)**br
1525 dum1 = 9.1_r8*rhof(i,k)
1526
1527 ! provisional rain number and mass weighted mean fallspeed (m/s)
1528
1529 umr(i,k) = min(dum1, dum*gamma_br_plus4*oneo6)
1530 unr(i,k) = min(dum1, dum*gamma_br_plus1)
1531 else
1532
1533 umr(i,k) = zero
1534 unr(i,k) = zero
1535 endif
1536 enddo
1537
1538 !......................................................................
1539 ! snow
1540
1541 call size_dist_param_basic(mg_snow_props, qsic(:,k), nsic(:,k), &
1542 lams(:,k), mgncol, n0=n0s(:,k))
1543
1544 do i=1,mgncol
1545 if (lams(i,k) >= qsmall) then
1546
1547 ! provisional snow number and mass weighted mean fallspeed (m/s)
1548
1549 dum = asn(i,k) / lams(i,k)**bs
1550 dum1 = 1.2_r8*rhof(i,k)
1551 ums(i,k) = min(dum1, dum*gamma_bs_plus4*oneo6)
1552 uns(i,k) = min(dum1, dum*gamma_bs_plus1)
1553
1554 else
1555 ums(i,k) = zero
1556 uns(i,k) = zero
1557 endif
1558 enddo
1559
1560 if (do_cldice) then
1561 if (.not. use_hetfrz_classnuc) then
1562
1563 ! heterogeneous freezing of cloud water
1564 !----------------------------------------------
1565
1566 call immersion_freezing(microp_uniform, t(:,k), pgam(:,k), lamc(:,k), &
1567 qcic(:,k), ncic(:,k), relvar(:,k), mnuccc(:,k), nnuccc(:,k), mgncol)
1568
1569 ! make sure number of droplets frozen does not exceed available ice nuclei concentration
1570 ! this prevents 'runaway' droplet freezing
1571
1572 where (qcic(:,k) >= qsmall .and. t(:,k) < 269.15_r8)
1573 where (nnuccc(:,k)*lcldm(:,k) > nnuccd(:,k))
1574 ! scale mixing ratio of droplet freezing with limit
1575 mnuccc(:,k) = mnuccc(:,k)*(nnuccd(:,k)/(nnuccc(:,k)*lcldm(:,k)))
1576 nnuccc(:,k) = nnuccd(:,k)/lcldm(:,k)
1577 end where
1578 end where
1579
1580 mdust = size(rndst,3)
1581 call contact_freezing(microp_uniform, t(:,k), p(:,k), rndst(:,k,:), &
1582 nacon(:,k,:), pgam(:,k), lamc(:,k), qcic(:,k), ncic(:,k), &
1583 relvar(:,k), mnucct(:,k), nnucct(:,k), mgncol, mdust)
1584
1585 mnudep(:,k) = zero
1586 nnudep(:,k) = zero
1587
1588 !else
1589
1590 ! Mass of droplets frozen is the average droplet mass, except
1591 ! with two limiters: concentration must be at least 1/cm^3, and
1592 ! mass must be at least the minimum defined above.
1593 !mi0l = qcic(:,k)/max(ncic(:,k), 1.0e6_r8/rho(:,k))
1594 !mi0l = max(mi0l_min, mi0l)
1595
1596 !where (qcic(:,k) >= qsmall)
1597 !nnuccc(:,k) = frzimm(:,k)*1.0e6_r8/rho(:,k)
1598 !mnuccc(:,k) = nnuccc(:,k)*mi0l
1599
1600 !nnucct(:,k) = frzcnt(:,k)*1.0e6_r8/rho(:,k)
1601 !mnucct(:,k) = nnucct(:,k)*mi0l
1602
1603 !nnudep(:,k) = frzdep(:,k)*1.0e6_r8/rho(:,k)
1604 !mnudep(:,k) = nnudep(:,k)*mi0
1605 !elsewhere
1606 !nnuccc(:,k) = 0._r8
1607 !mnuccc(:,k) = 0._r8
1608
1609 !nnucct(:,k) = 0._r8
1610 !mnucct(:,k) = 0._r8
1611
1612 !nnudep(:,k) = 0._r8
1613 !mnudep(:,k) = 0._r8
1614 !end where
1615
1616 end if
1617
1618 else
1619 do i=1,mgncol
1620 mnuccc(i,k) = zero
1621 nnuccc(i,k) = zero
1622 mnucct(i,k) = zero
1623 nnucct(i,k) = zero
1624 mnudep(i,k) = zero
1625 nnudep(i,k) = zero
1626 enddo
1627 end if
1628
1629 call snow_self_aggregation(t(:,k), rho(:,k), asn(:,k), rhosn, qsic(:,k), nsic(:,k), &
1630 nsagg(:,k), mgncol)
1631
1632 call accrete_cloud_water_snow(t(:,k), rho(:,k), asn(:,k), uns(:,k), mu(:,k), &
1633 qcic(:,k), ncic(:,k), qsic(:,k), pgam(:,k), lamc(:,k), lams(:,k), n0s(:,k), &
1634 psacws(:,k), npsacws(:,k), mgncol)
1635
1636 if (do_cldice) then
1637 call secondary_ice_production(t(:,k), psacws(:,k), msacwi(:,k), nsacwi(:,k), mgncol)
1638 else
1639 nsacwi(:,k) = zero
1640 msacwi(:,k) = zero
1641 end if
1642
1643 call accrete_rain_snow(t(:,k), rho(:,k), umr(:,k), ums(:,k), unr(:,k), uns(:,k), &
1644 qric(:,k), qsic(:,k), lamr(:,k), n0r(:,k), lams(:,k), n0s(:,k), &
1645 pracs(:,k), npracs(:,k), mgncol)
1646
1647 call heterogeneous_rain_freezing(t(:,k), qric(:,k), nric(:,k), lamr(:,k), &
1648 mnuccr(:,k), nnuccr(:,k), mgncol)
1649
1650 if (do_sb_physics) then
1651 call sb2001v2_accre_cld_water_rain(qcic(:,k), ncic(:,k), qric(:,k), &
1652 rho(:,k), relvar(:,k), pra(:,k), npra(:,k), mgncol)
1653 else
1654 call accrete_cloud_water_rain(microp_uniform, qric(:,k), qcic(:,k), &
1655 ncic(:,k), relvar(:,k), accre_enhan(:,k), pra(:,k), npra(:,k), mgncol)
1656 endif
1657
1658 call self_collection_rain(rho(:,k), qric(:,k), nric(:,k), nragg(:,k), mgncol)
1659
1660 if (do_cldice) then
1661 call accrete_cloud_ice_snow(t(:,k), rho(:,k), asn(:,k), qiic(:,k), niic(:,k), &
1662 qsic(:,k), lams(:,k), n0s(:,k), prai(:,k), nprai(:,k), mgncol)
1663 else
1664 prai(:,k) = zero
1665 nprai(:,k) = zero
1666 end if
1667
1668 call evaporate_sublimate_precip(t(:,k), rho(:,k), &
1669 dv(:,k), mu(:,k), sc(:,k), q(:,k), qvl(:,k), qvi(:,k), &
1670 lcldm(:,k), precip_frac(:,k), arn(:,k), asn(:,k), qcic(:,k), qiic(:,k), &
1671 qric(:,k), qsic(:,k), lamr(:,k), n0r(:,k), lams(:,k), n0s(:,k), &
1672 pre(:,k), prds(:,k), am_evp_st(:,k), mgncol)
1673
1674 call bergeron_process_snow(t(:,k), rho(:,k), dv(:,k), mu(:,k), sc(:,k), &
1675 qvl(:,k), qvi(:,k), asn(:,k), qcic(:,k), qsic(:,k), lams(:,k), n0s(:,k), &
1676 bergs(:,k), mgncol)
1677
1678 bergs(:,k)=bergs(:,k)*micro_mg_berg_eff_factor
1679
1680 !+++PMC 12/3/12 - NEW VAPOR DEP/SUBLIMATION GOES HERE!!!
1681 if (do_cldice) then
1682
1683 call ice_deposition_sublimation(t(:,k), q(:,k), qi(:,k), ni(:,k), &
1684 icldm(:,k), rho(:,k), dv(:,k), qvl(:,k), qvi(:,k), &
1685 berg(:,k), vap_dep(:,k), ice_sublim(:,k), mgncol)
1686
1687 do i=1,mgncol
1688! sublimation should not exceed available ice
1689 ice_sublim(i,k) = max(ice_sublim(i,k), -qi(i,k)*oneodt)
1690
1691 berg(i,k) = berg(i,k)*micro_mg_berg_eff_factor
1692
1693 if (vap_dep(i,k) < zero .and. qi(i,k) > qsmall .and. icldm(i,k) > mincld) then
1694 nsubi(i,k) = vap_dep(i,k) * ni(i,k) / (qi(i,k) * icldm(i,k))
1695 else
1696 nsubi(i,k) = zero
1697 endif
1698
1699 ! bergeron process should not reduce nc unless
1700 ! all ql is removed (which is handled elsewhere)
1701 !in fact, nothing in this entire file makes nsubc nonzero.
1702
1703 nsubc(i,k) = zero
1704 enddo
1705
1706 end if !do_cldice
1707 !---PMC 12/3/12
1708
1709 do i=1,mgncol
1710
1711 ! conservation to ensure no negative values of cloud water/precipitation
1712 ! in case microphysical process rates are large
1713 !===================================================================
1714
1715 ! note: for check on conservation, processes are multiplied by omsm
1716 ! to prevent problems due to round off error
1717
1718 ! conservation of qc
1719 !-------------------------------------------------------------------
1720
1721 dum = ((prc(i,k)+pra(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k)+ &
1722 psacws(i,k)+bergs(i,k))*lcldm(i,k)+berg(i,k))*deltat
1723
1724 if (dum > qc(i,k)) then
1725 ratio = qc(i,k) / dum * omsm
1726
1727 prc(i,k) = ratio * prc(i,k)
1728 pra(i,k) = ratio * pra(i,k)
1729 mnuccc(i,k) = ratio * mnuccc(i,k)
1730 mnucct(i,k) = ratio * mnucct(i,k)
1731 msacwi(i,k) = ratio * msacwi(i,k)
1732 psacws(i,k) = ratio * psacws(i,k)
1733 bergs(i,k) = ratio * bergs(i,k)
1734 berg(i,k) = ratio * berg(i,k)
1735 qcrat(i,k) = ratio
1736 else
1737 qcrat(i,k) = one
1738 end if
1739
1740 !PMC 12/3/12: ratio is also frac of step w/ liquid.
1741 !thus we apply berg for "ratio" of timestep and vapor
1742 !deposition for the remaining frac of the timestep.
1743 if (qc(i,k) >= qsmall) then
1744 vap_dep(i,k) = vap_dep(i,k)*(1._r8-qcrat(i,k))
1745 end if
1746
1747 end do
1748
1749 do i=1,mgncol
1750
1751 !=================================================================
1752 ! apply limiter to ensure that ice/snow sublimation and rain evap
1753 ! don't push conditions into supersaturation, and ice deposition/nucleation don't
1754 ! push conditions into sub-saturation
1755 ! note this is done after qc conservation since we don't know how large
1756 ! vap_dep is before then
1757 ! estimates are only approximate since other process terms haven't been limited
1758 ! for conservation yet
1759
1760 ! first limit ice deposition/nucleation vap_dep + mnuccd
1761
1762 dum1 = vap_dep(i,k) + mnuccd(i,k)
1763 if (dum1 > 1.e-20_r8) then
1764 dum = (q(i,k)-qvi(i,k))/(one + xxls_squared*qvi(i,k)/(cpp*rv*t(i,k)*t(i,k)))*oneodt
1765 dum = max(dum, zero)
1766 if (dum1 > dum) then
1767 ! Allocate the limited "dum" tendency to mnuccd and vap_dep
1768 ! processes. Don't divide by cloud fraction; these are grid-
1769 ! mean rates.
1770 dum1 = mnuccd(i,k) / (vap_dep(i,k)+mnuccd(i,k))
1771 mnuccd(i,k) = dum*dum1
1772 vap_dep(i,k) = dum - mnuccd(i,k)
1773 end if
1774 end if
1775
1776 end do
1777
1778 do i=1,mgncol
1779
1780 !===================================================================
1781 ! conservation of nc
1782 !-------------------------------------------------------------------
1783 dum = (nprc1(i,k)+npra(i,k)+nnuccc(i,k)+nnucct(i,k)+ &
1784 npsacws(i,k)-nsubc(i,k))*lcldm(i,k) * deltat
1785
1786 if (dum > nc(i,k)) then
1787 ratio = nc(i,k) / dum * omsm
1788
1789 nprc1(i,k) = ratio * nprc1(i,k)
1790 npra(i,k) = ratio * npra(i,k)
1791 nnuccc(i,k) = ratio * nnuccc(i,k)
1792 nnucct(i,k) = ratio * nnucct(i,k)
1793 npsacws(i,k) = ratio * npsacws(i,k)
1794 nsubc(i,k) = ratio * nsubc(i,k)
1795 endif
1796
1797 mnuccri(i,k) = zero
1798 nnuccri(i,k) = zero
1799
1800 if (do_cldice) then
1801
1802 ! freezing of rain to produce ice if mean rain size is smaller than Dcs
1803 if (lamr(i,k) > qsmall) then
1804 if (one/lamr(i,k) < dcs) then
1805 mnuccri(i,k) = mnuccr(i,k)
1806 nnuccri(i,k) = nnuccr(i,k)
1807 mnuccr(i,k) = zero
1808 nnuccr(i,k) = zero
1809 endif
1810 endif
1811 endif
1812
1813 enddo
1814
1815 do i=1,mgncol
1816
1817 ! conservation of rain mixing ratio
1818 !-------------------------------------------------------------------
1819 dum1 = -pre(i,k) + pracs(i,k) + mnuccr(i,k) + mnuccri(i,k)
1820 dum3 = dum1 * precip_frac(i,k)
1821 dum2 = (pra(i,k)+prc(i,k))*lcldm(i,k)
1822 dum = (dum3 - dum2) * deltat
1823
1824 ! note that qrtend is included below because of instantaneous freezing/melt
1825 if (dum > qr(i,k) .and. dum1 >= qsmall) then
1826 ratio = (qr(i,k)*oneodt + dum2) / dum3 * omsm
1827 pre(i,k) = ratio * pre(i,k)
1828 pracs(i,k) = ratio * pracs(i,k)
1829 mnuccr(i,k) = ratio * mnuccr(i,k)
1830 mnuccri(i,k) = ratio * mnuccri(i,k)
1831 end if
1832
1833 end do
1834
1835 do i=1,mgncol
1836
1837 ! conservation of rain number
1838 !-------------------------------------------------------------------
1839
1840 ! Add evaporation of rain number.
1841 if (pre(i,k) < zero) then
1842 dum = max(-one, pre(i,k)*deltat/qr(i,k))
1843 nsubr(i,k) = dum*nr(i,k) * oneodt
1844 else
1845 nsubr(i,k) = zero
1846 end if
1847
1848 end do
1849
1850 do i=1,mgncol
1851
1852 dum1 = (-nsubr(i,k)+npracs(i,k)+nnuccr(i,k)+nnuccri(i,k)-nragg(i,k))*precip_frac(i,k)
1853 dum2 = nprc(i,k)*lcldm(i,k)
1854 dum = (dum1 - dum2) * deltat
1855
1856 if (dum > nr(i,k)) then
1857 ratio = (nr(i,k)*oneodt + dum2) / dum1 * omsm
1858
1859 nragg(i,k) = ratio * nragg(i,k)
1860 npracs(i,k) = ratio * npracs(i,k)
1861 nnuccr(i,k) = ratio * nnuccr(i,k)
1862 nsubr(i,k) = ratio * nsubr(i,k)
1863 nnuccri(i,k) = ratio * nnuccri(i,k)
1864 end if
1865
1866 end do
1867
1868 if (do_cldice) then
1869
1870 do i=1,mgncol
1871
1872 ! conservation of qi
1873 !-------------------------------------------------------------------
1874
1875 dum1 = (prci(i,k)+prai(i,k))*icldm(i,k)-ice_sublim(i,k)
1876 dum2 = vap_dep(i,k)+berg(i,k)+mnuccd(i,k) &
1877 + (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k))*lcldm(i,k) &
1878 + mnuccri(i,k)*precip_frac(i,k)
1879 dum = (dum1 - dum2) * deltat
1880
1881 if (dum > qi(i,k)) then
1882 ratio = (qi(i,k)*oneodt + dum2) / dum1 * omsm
1883
1884 prci(i,k) = ratio * prci(i,k)
1885 prai(i,k) = ratio * prai(i,k)
1886 ice_sublim(i,k) = ratio * ice_sublim(i,k)
1887 end if
1888
1889 end do
1890
1891 end if
1892
1893 if (do_cldice) then
1894
1895 do i=1,mgncol
1896
1897 ! conservation of ni
1898 !-------------------------------------------------------------------
1899 if (use_hetfrz_classnuc) then
1900 tmpfrz = nnuccc(i,k)
1901 else
1902 tmpfrz = zero
1903 end if
1904 dum1 = (nprci(i,k)+nprai(i,k)-nsubi(i,k))*icldm(i,k)
1905 dum2 = nnuccd(i,k)+(nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k))*lcldm(i,k) &
1906 + nnuccri(i,k)*precip_frac(i,k)
1907 dum = (dum1 - dum2) * deltat
1908
1909 if (dum > ni(i,k)) then
1910 ratio = (ni(i,k)*oneodt + dum2) / dum1 * omsm
1911
1912 nprci(i,k) = ratio * nprci(i,k)
1913 nprai(i,k) = ratio * nprai(i,k)
1914 nsubi(i,k) = ratio * nsubi(i,k)
1915 end if
1916
1917 end do
1918
1919 end if
1920
1921 do i=1,mgncol
1922
1923 ! conservation of snow mixing ratio
1924 !-------------------------------------------------------------------
1925 dum1 = - prds(i,k) * precip_frac(i,k)
1926 dum2 = (pracs(i,k)+mnuccr(i,k))*precip_frac(i,k) &
1927 + (prai(i,k)+prci(i,k))*icldm(i,k) + (bergs(i,k)+psacws(i,k))*lcldm(i,k)
1928 dum = (dum1 - dum2) * deltat
1929
1930 if (dum > qs(i,k) .and. -prds(i,k) >= qsmall) then
1931 ratio = (qs(i,k)*oneodt + dum2) / dum1 * omsm
1932
1933 prds(i,k) = ratio * prds(i,k)
1934 end if
1935
1936 end do
1937
1938 do i=1,mgncol
1939
1940 ! conservation of snow number
1941 !-------------------------------------------------------------------
1942 ! calculate loss of number due to sublimation
1943 ! for now neglect sublimation of ns
1944 nsubs(i,k) = zero
1945
1946 dum1 = precip_frac(i,k)* (-nsubs(i,k)-nsagg(i,k))
1947 dum2 = nnuccr(i,k)*precip_frac(i,k) + nprci(i,k)*icldm(i,k)
1948 dum = (dum1 - dum2) * deltat
1949
1950 if (dum > ns(i,k)) then
1951 ratio = (ns(i,k)*oneodt + dum2) / dum1 * omsm
1952
1953 nsubs(i,k) = ratio * nsubs(i,k)
1954 nsagg(i,k) = ratio * nsagg(i,k)
1955 end if
1956
1957 end do
1958
1959 do i=1,mgncol
1960
1961 ! next limit ice and snow sublimation and rain evaporation
1962 ! get estimate of q and t at end of time step
1963 ! don't include other microphysical processes since they haven't
1964 ! been limited via conservation checks yet
1965
1966 tx1 = pre(i,k) * precip_frac(i,k)
1967 tx2 = prds(i,k) * precip_frac(i,k)
1968 tx3 = tx1 + tx2 + ice_sublim(i,k)
1969 if (tx3 < -1.e-20_r8) then
1970
1971 tx4 = tx2 + ice_sublim(i,k) + vap_dep(i,k) + mnuccd(i,k)
1972 qtmp = q(i,k) - (tx1 + tx4) * deltat
1973 ttmp = t(i,k) + (tx1*xxlv + tx4*xxls) * (deltat/cpp)
1974
1975 ! use rhw to allow ice supersaturation
1976 ! call qsat_water(ttmp, p(i,k), esn, qvn)
1977 esn = min(fpvsl(ttmp), p(i,k))
1978 qvn = epsqs*esn/(p(i,k)-omeps*esn) * qsfm(i,k)
1979
1980 ! modify ice/precip evaporation rate if q > qsat
1981 if (qtmp > qvn) then
1982
1983 tx4 = one / tx3
1984 dum1 = tx1 * tx4
1985 dum2 = tx2 * tx4
1986 ! recalculate q and t after vap_dep and mnuccd but without evap or sublim
1987 tx5 = (vap_dep(i,k)+mnuccd(i,k)) * deltat
1988 qtmp = q(i,k) - tx5
1989 ttmp = t(i,k) + tx5 * (xxls/cpp)
1990
1991 ! use rhw to allow ice supersaturation
1992 !call qsat_water(ttmp, p(i,k), esn, qvn)
1993 esn = min(fpvsl(ttmp), p(i,k))
1994 qvn = epsqs*esn / (p(i,k)-omeps*esn) * qsfm(i,k)
1995
1996 dum = (qtmp-qvn) / (one + xxlv_squared*qvn/(cpp*rv*ttmp*ttmp))
1997 dum = min(dum, zero)
1998
1999 ! modify rates if needed, divide by precip_frac to get local (in-precip) value
2000 if (precip_frac(i,k) > mincld) then
2001 tx4 = oneodt / precip_frac(i,k)
2002 else
2003 tx4 = zero
2004 endif
2005 pre(i,k) = dum*dum1*tx4
2006
2007 ! do separately using RHI for prds and ice_sublim
2008 !call qsat_ice(ttmp, p(i,k), esn, qvn)
2009 esn = min(fpvsi(ttmp), p(i,k))
2010 qvn = epsqs*esn / (p(i,k)-omeps*esn) * qsfm(i,k)
2011
2012
2013 dum = (qtmp-qvn) / (one + xxls_squared*qvn/(cpp*rv*ttmp*ttmp))
2014 dum = min(dum, zero)
2015
2016 ! modify rates if needed, divide by precip_frac to get local (in-precip) value
2017 prds(i,k) = dum*dum2*tx4
2018
2019 ! don't divide ice_sublim by cloud fraction since it is grid-averaged
2020 dum1 = one - dum1 - dum2
2021 ice_sublim(i,k) = dum*dum1*oneodt
2022 end if
2023 end if
2024
2025 end do
2026
2027 ! Big "administration" loop enforces conservation, updates variables
2028 ! that accumulate over substeps, and sets output variables.
2029
2030 do i=1,mgncol
2031
2032 ! get tendencies due to microphysical conversion processes
2033 !==========================================================
2034 ! note: tendencies are multiplied by appropriate cloud/precip
2035 ! fraction to get grid-scale values
2036 ! note: vap_dep is already grid-average values
2037
2038 ! The net tendencies need to be added to rather than overwritten,
2039 ! because they may have a value already set for instantaneous
2040 ! melting/freezing.
2041
2042
2043 qvlat(i,k) = qvlat(i,k) - (pre(i,k)+prds(i,k))*precip_frac(i,k)-&
2044 vap_dep(i,k)-ice_sublim(i,k)-mnuccd(i,k)-mnudep(i,k)*lcldm(i,k)
2045
2046
2047! if (lprnt .and. k >= 60 ) &
2048! write(0,*)' k=',k,' tlat=',tlat(i,k),' pre=',pre(i,k),' precip_frac=',precip_frac(i,k),&
2049! ' prds=',prds(i,k),' vap_dep=',vap_dep(i,k),' ice_sublim=',ice_sublim(i,k), &
2050! ' mnuccd=',mnuccd(i,k),' mnudep=',mnudep(i,k),' lcldm=',lcldm(i,k),' bergs=',bergs(i,k), &
2051! ' psacws=',psacws(i,k),' mnuccc=',mnuccc(i,k),' mnucct=',mnucct(i,k),' msacwi=',msacwi(i,k), &
2052! ' mnuccr=',mnuccr(i,k), &
2053! ' pracs=',pracs(i,k),' mnuccri=',mnuccri(i,k),' xlf=',xlf,' xxlv=',xxlv,' xxls=',xxls
2054
2055 tlat(i,k) = tlat(i,k) + ((pre(i,k)*precip_frac(i,k)) &
2056 *xxlv+(prds(i,k)*precip_frac(i,k)+vap_dep(i,k)+ice_sublim(i,k)+mnuccd(i,k)+mnudep(i,k)*lcldm(i,k))*xxls+ &
2057 ((bergs(i,k)+psacws(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k))*lcldm(i,k)+(mnuccr(i,k)+ &
2058 pracs(i,k)+mnuccri(i,k))*precip_frac(i,k)+berg(i,k))*xlf)
2059
2060! if (lprnt .and. k >= 60) write(0,*)' k=',k,' tlat=',tlat(i,k)
2061
2062 qctend(i,k) = qctend(i,k) + (-pra(i,k)-prc(i,k)-mnuccc(i,k)-mnucct(i,k)-msacwi(i,k)- &
2063 psacws(i,k)-bergs(i,k))*lcldm(i,k)-berg(i,k)
2064
2065 if (do_cldice) then
2066 qitend(i,k) = qitend(i,k) + &
2067 (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k))*lcldm(i,k)+(-prci(i,k)- &
2068 prai(i,k))*icldm(i,k)+vap_dep(i,k)+berg(i,k)+ice_sublim(i,k)+ &
2069 mnuccd(i,k)+mnuccri(i,k)*precip_frac(i,k)
2070 end if
2071
2072 qrtend(i,k) = qrtend(i,k) + (pra(i,k)+prc(i,k))*lcldm(i,k)+(pre(i,k)-pracs(i,k)- &
2073 mnuccr(i,k)-mnuccri(i,k))*precip_frac(i,k)
2074
2075 qstend(i,k) = qstend(i,k) + (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k) &
2076 + (prds(i,k)+pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)
2077
2078
2079 cmeout(i,k) = vap_dep(i,k) + ice_sublim(i,k) + mnuccd(i,k)
2080
2081 ! add output for cmei (accumulate)
2082 cmeitot(i,k) = vap_dep(i,k) + ice_sublim(i,k) + mnuccd(i,k)
2083
2084 ! assign variables for trop_mozart, these are grid-average
2085 !-------------------------------------------------------------------
2086 ! evaporation/sublimation is stored here as positive term
2087
2088 evapsnow(i,k) = -prds(i,k) * precip_frac(i,k)
2089 nevapr(i,k) = -pre(i,k) * precip_frac(i,k)
2090 prer_evap(i,k) = -pre(i,k) * precip_frac(i,k)
2091
2092 ! change to make sure prain is positive: do not remove snow from
2093 ! prain used for wet deposition
2094 prain(i,k) = (pra(i,k)+prc(i,k))*lcldm(i,k)+(-pracs(i,k)- &
2095 mnuccr(i,k)-mnuccri(i,k))*precip_frac(i,k)
2096 prodsnow(i,k) = (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+(&
2097 pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)
2098
2099 ! following are used to calculate 1st order conversion rate of cloud water
2100 ! to rain and snow (1/s), for later use in aerosol wet removal routine
2101 ! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc
2102 ! used to calculate pra, prc, ... in this routine
2103 ! qcsinksum_rate1ord = { rate of direct transfer of cloud water to rain & snow }
2104 ! (no cloud ice or bergeron terms)
2105 qcsinksum_rate1ord(i,k) = (pra(i,k)+prc(i,k)+psacws(i,k))*lcldm(i,k)
2106 ! Avoid zero/near-zero division.
2107 qcsinksum_rate1ord(i,k) = qcsinksum_rate1ord(i,k) / max(qc(i,k),1.0e-30_r8)
2108
2109
2110 ! microphysics output, note this is grid-averaged
2111 pratot(i,k) = pra(i,k) * lcldm(i,k)
2112 prctot(i,k) = prc(i,k) * lcldm(i,k)
2113 mnuccctot(i,k) = mnuccc(i,k) * lcldm(i,k)
2114 mnuccttot(i,k) = mnucct(i,k) * lcldm(i,k)
2115 msacwitot(i,k) = msacwi(i,k) * lcldm(i,k)
2116 psacwstot(i,k) = psacws(i,k) * lcldm(i,k)
2117 bergstot(i,k) = bergs(i,k) * lcldm(i,k)
2118 bergtot(i,k) = berg(i,k)
2119 prcitot(i,k) = prci(i,k) * icldm(i,k)
2120 praitot(i,k) = prai(i,k) * icldm(i,k)
2121 mnuccdtot(i,k) = mnuccd(i,k) * icldm(i,k)
2122
2123 pracstot(i,k) = pracs(i,k) * precip_frac(i,k)
2124 mnuccrtot(i,k) = mnuccr(i,k) * precip_frac(i,k)
2125
2126
2127 nctend(i,k) = nctend(i,k) + (-nnuccc(i,k)-nnucct(i,k)-npsacws(i,k)+nsubc(i,k) &
2128 - npra(i,k)-nprc1(i,k))*lcldm(i,k)
2129
2130 if (do_cldice) then
2131 if (use_hetfrz_classnuc) then
2132 tmpfrz = nnuccc(i,k)
2133 else
2134 tmpfrz = zero
2135 end if
2136 nitend(i,k) = nitend(i,k) + nnuccd(i,k)+ &
2137 (nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k))*lcldm(i,k)+(nsubi(i,k)-nprci(i,k)- &
2138 nprai(i,k))*icldm(i,k)+nnuccri(i,k)*precip_frac(i,k)
2139 end if
2140
2141 nstend(i,k) = nstend(i,k) + (nsubs(i,k)+nsagg(i,k)+nnuccr(i,k))*precip_frac(i,k) &
2142 + nprci(i,k)*icldm(i,k)
2143
2144 nrtend(i,k) = nrtend(i,k) + nprc(i,k)*lcldm(i,k)+(nsubr(i,k)-npracs(i,k)-nnuccr(i,k) &
2145 - nnuccri(i,k)+nragg(i,k))*precip_frac(i,k)
2146
2147 ! make sure that ni at advanced time step does not exceed
2148 ! maximum (existing N + source terms*dt), which is possible if mtime < deltat
2149 ! note that currently mtime = deltat
2150 !================================================================
2151
2152 if (do_cldice .and. nitend(i,k) > zero .and. ni(i,k)+nitend(i,k)*deltat > nimax(i,k)) then
2153 nitend(i,k) = max(zero, (nimax(i,k)-ni(i,k))*oneodt)
2154 end if
2155
2156 end do
2157
2158 ! End of "administration" loop
2159
2160 end do micro_vert_loop ! end k loop
2161
2162! if (lprnt) write(0,*)' tlat3=',tlat(1,:)*deltat
2163 !-----------------------------------------------------
2164 ! convert rain/snow q and N for output to history, note,
2165 ! output is for gridbox average
2166
2167 do k=1,nlev
2168 do i=1,mgncol
2169 qrout(i,k) = qr(i,k)
2170 nrout(i,k) = nr(i,k) * rho(i,k)
2171 qsout(i,k) = qs(i,k)
2172 nsout(i,k) = ns(i,k) * rho(i,k)
2173 enddo
2174 enddo
2175
2176 ! calculate n0r and lamr from rain mass and number
2177 ! divide by precip fraction to get in-precip (local) values of
2178 ! rain mass and number, divide by rhow to get rain number in kg^-1
2179
2180 do k=1,nlev
2181
2182 call size_dist_param_basic(mg_rain_props, qric(:,k), nric(:,k), lamr(:,k), mgncol, n0=n0r(:,k))
2183
2184 enddo
2185 ! Calculate rercld
2186
2187 ! calculate mean size of combined rain and cloud water
2188
2189 call calc_rercld(lamr, n0r, lamc, pgam, qric, qcic, ncic, rercld, mgncol, nlev)
2190
2191
2192 ! Assign variables back to start-of-timestep values
2193 ! Some state variables are changed before the main microphysics loop
2194 ! to make "instantaneous" adjustments. Afterward, we must move those changes
2195 ! back into the tendencies.
2196 ! These processes:
2197 ! - Droplet activation (npccn, impacts nc)
2198 ! - Instantaneous snow melting (minstsm/ninstsm, impacts qr/qs/nr/ns)
2199 ! - Instantaneous rain freezing (minstfr/ninstrf, impacts qr/qs/nr/ns)
2200 !================================================================================
2201
2202 do k=1,nlev
2203 do i=1,mgncol
2204 ! Re-apply droplet activation tendency
2205 nc(i,k) = ncn(i,k)
2206 nctend(i,k) = nctend(i,k) + npccn(i,k)
2207
2208 ! Re-apply rain freezing and snow melting.
2209 qstend(i,k) = qstend(i,k) + (qs(i,k)-qsn(i,k)) * oneodt
2210 qs(i,k) = qsn(i,k)
2211
2212 nstend(i,k) = nstend(i,k) + (ns(i,k)-nsn(i,k)) * oneodt
2213 ns(i,k) = nsn(i,k)
2214
2215 qrtend(i,k) = qrtend(i,k) + (qr(i,k)-qrn(i,k)) * oneodt
2216 qr(i,k) = qrn(i,k)
2217
2218 nrtend(i,k) = nrtend(i,k) + (nr(i,k)-nrn(i,k)) * oneodt
2219 nr(i,k) = nrn(i,k)
2220
2221 !.............................................................................
2222
2223 !================================================================================
2224
2225 ! modify to include snow. in prain & evap (diagnostic here: for wet dep)
2226 nevapr(i,k) = nevapr(i,k) + evapsnow(i,k)
2227 prain(i,k) = prain(i,k) + prodsnow(i,k)
2228
2229 enddo
2230 enddo
2231
2232 do k=1,nlev
2233
2234 do i=1,mgncol
2235
2236 ! calculate sedimentation for cloud water and ice
2237 !================================================================================
2238
2239 ! update in-cloud cloud mixing ratio and number concentration
2240 ! with microphysical tendencies to calculate sedimentation, assign to dummy vars
2241 ! note: these are in-cloud values***, hence we divide by cloud fraction
2242
2243 if (lcldm(i,k) > mincld) then
2244 tx1 = one / lcldm(i,k)
2245 dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat) * tx1
2246 dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)*tx1, zero)
2247 else
2248 dumc(i,k) = zero
2249 dumnc(i,k) = zero
2250 endif
2251 if (icldm(i,k) > mincld) then
2252 tx1 = one / icldm(i,k)
2253 dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat) * tx1
2254 dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)*tx1, zero)
2255 else
2256 dumi(i,k) = zero
2257 dumni(i,k) = zero
2258 endif
2259 if (precip_frac(i,k) > mincld) then
2260 tx1 = one / precip_frac(i,k)
2261 dumr(i,k) = (qr(i,k)+qrtend(i,k)*deltat) * tx1
2262 dums(i,k) = (qs(i,k)+qstend(i,k)*deltat) * tx1
2263
2264 dumnr(i,k) = max((nr(i,k)+nrtend(i,k)*deltat)*tx1, zero)
2265 dumns(i,k) = max((ns(i,k)+nstend(i,k)*deltat)*tx1, zero)
2266 else
2267 dumr(i,k) = zero
2268 dumr(i,k) = zero
2269 dums(i,k) = zero
2270 dumns(i,k) = zero
2271 endif
2272
2273 ! switch for specification of droplet and crystal number
2274 if (nccons) then
2275 dumnc(i,k) = ncnst*rhoinv(i,k)
2276 end if
2277
2278 ! switch for specification of cloud ice number
2279 if (nicons) then
2280 dumni(i,k) = ninst*rhoinv(i,k)
2281 end if
2282 enddo
2283 enddo
2284
2285 do k=1,nlev
2286
2287! obtain new slope parameter to avoid possible singularity
2288
2289
2290 call size_dist_param_ice(mg_ice_props, dumi(:,k), dumni(:,k), &
2291 lami(:,k), mgncol)
2292
2293 call size_dist_param_liq(mg_liq_props, dumc(:,k), dumnc(:,k), rho(:,k), &
2294 pgam(:,k), lamc(:,k), mgncol)
2295
2296! call size_dist_param_basic(mg_ice_props, dumi(:,k), dumni(:,k), &
2297! lami(:,k), mgncol)
2298! fallspeed for rain
2299
2300 call size_dist_param_basic(mg_rain_props, dumr(:,k), dumnr(:,k), &
2301 lamr(:,k), mgncol)
2302! fallspeed for snow
2303 call size_dist_param_basic(mg_snow_props, dums(:,k), dumns(:,k), &
2304 lams(:,k), mgncol)
2305
2306 enddo
2307
2308 do k=1,nlev
2309 do i=1,mgncol
2310
2311 ! calculate number and mass weighted fall velocity for droplets and cloud ice
2312 !-------------------------------------------------------------------
2313
2314 grho = g*rho(i,k)
2315
2316 if (dumc(i,k) >= qsmall) then
2317
2318 tx1 = lamc(i,k)**bc
2319 vtrmc(i,k) = acn(i,k)*gamma(four+bc+pgam(i,k)) &
2320 / (tx1*gamma(pgam(i,k)+four))
2321
2322 fc(i,k) = grho*vtrmc(i,k)
2323
2324 fnc(i,k) = grho* acn(i,k)*gamma(pgam(i,k)+one+bc) &
2325 / (tx1*gamma(pgam(i,k)+one))
2326 else
2327 fc(i,k) = zero
2328 fnc(i,k) = zero
2329 end if
2330
2331 ! calculate number and mass weighted fall velocity for cloud ice
2332
2333 if (dumi(i,k) >= qsmall) then
2334
2335 tx3 = one / lami(i,k)
2336 tx1 = ain(i,k) * tx3**bi
2337 tx2 = 1.2_r8*rhof(i,k)
2338 vtrmi(i,k) = min(tx1*gamma_bi_plus4*oneo6, tx2)
2339
2340 fi(i,k) = grho * vtrmi(i,k)
2341 fni(i,k) = grho * min(tx1*gamma_bi_plus1, tx2)
2342
2343 ! adjust the ice fall velocity for smaller (r < 20 um) ice
2344 ! particles (blend over 18-20 um)
2345 irad = (1.5_r8 * 1e6_r8) * tx3
2346 ifrac = min(one, max(zero, (irad-18._r8)*half))
2347
2348 if (ifrac < one) then
2349 tx1 = ajn(i,k) / lami(i,k)**bj
2350 vtrmi(i,k) = ifrac*vtrmi(i,k) + (one-ifrac) * min(tx1*gamma_bj_plus4*oneo6, tx2)
2351
2352 fi(i,k) = grho*vtrmi(i,k)
2353 fni(i,k) = ifrac * fni(i,k) + (one-ifrac) * grho * min(tx1*gamma_bj_plus1, tx2)
2354 end if
2355 else
2356 fi(i,k) = zero
2357 fni(i,k)= zero
2358 end if
2359
2360
2361 ! fallspeed for rain
2362
2363! if (lamr(i,k) >= qsmall) then
2364 if (dumr(i,k) >= qsmall) then
2365
2366 ! 'final' values of number and mass weighted mean fallspeed for rain (m/s)
2367
2368 tx1 = arn(i,k) / lamr(i,k)**br
2369 tx2 = 9.1_r8*rhof(i,k)
2370 umr(i,k) = min(tx1*gamma_br_plus4*oneo6, tx2)
2371 unr(i,k) = min(tx1*gamma_br_plus1, tx2)
2372
2373 fr(i,k) = grho * umr(i,k)
2374 fnr(i,k) = grho * unr(i,k)
2375
2376 else
2377 fr(i,k) = zero
2378 fnr(i,k) = zero
2379 end if
2380
2381 ! fallspeed for snow
2382
2383
2384! if (lams(i,k) >= qsmall) then
2385 if (dums(i,k) >= qsmall) then
2386
2387 ! 'final' values of number and mass weighted mean fallspeed for snow (m/s)
2388 tx1 = asn(i,k) / lams(i,k)**bs
2389 tx2 = 1.2_r8*rhof(i,k)
2390 ums(i,k) = min(tx1*gamma_bs_plus4*oneo6, tx2)
2391 uns(i,k) = min(tx1*gamma_bs_plus1, tx2)
2392
2393 fs(i,k) = grho * ums(i,k)
2394 fns(i,k) = grho * uns(i,k)
2395
2396 else
2397 fs(i,k) = zero
2398 fns(i,k) = zero
2399 end if
2400
2401 ! redefine dummy variables - sedimentation is calculated over grid-scale
2402 ! quantities to ensure conservation
2403
2404 dumc(i,k) = qc(i,k) + qctend(i,k)*deltat
2405 dumi(i,k) = qi(i,k) + qitend(i,k)*deltat
2406 dumr(i,k) = qr(i,k) + qrtend(i,k)*deltat
2407 dums(i,k) = qs(i,k) + qstend(i,k)*deltat
2408
2409 dumnc(i,k) = nc(i,k) + nctend(i,k)*deltat
2410 dumni(i,k) = ni(i,k) + nitend(i,k)*deltat
2411 dumnr(i,k) = nr(i,k) + nrtend(i,k)*deltat
2412 dumns(i,k) = ns(i,k) + nstend(i,k)*deltat
2413
2414 if (dumc(i,k) < qsmall) dumnc(i,k) = zero
2415 if (dumi(i,k) < qsmall) dumni(i,k) = zero
2416 if (dumr(i,k) < qsmall) dumnr(i,k) = zero
2417 if (dums(i,k) < qsmall) dumns(i,k) = zero
2418
2419 enddo
2420 end do !!! vertical loop
2421
2422 do k=1,nlev
2423 do i=1,mgncol
2424 pdel_inv(i,k) = one / pdel(i,k)
2425 enddo
2426 enddo
2427! if (lprnt) write(0,*)' bef sedimentation dumc=',dumc(i,nlev-10:nlev)
2428
2429 ! initialize nstep for sedimentation sub-steps
2430
2431 ! calculate number of split time steps to ensure courant stability criteria
2432 ! for sedimentation calculations
2433 !-------------------------------------------------------------------
2434 do i=1,mgncol
2435 nlb = nlball(i)
2436 nstep = 1 + nint(max( maxval( fi(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
2437 maxval(fni(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)
2438
2439 nstep = min(nstep, nstep_def)
2440
2441
2442 ! loop over sedimentation sub-time step to ensure stability
2443 !==============================================================
2444 if (do_cldice) then
2445 tx2 = one / nstep
2446 tx1 = tx2 * deltat
2447 tx3 = tx2 / g
2448
2449 do n = 1,nstep
2450
2451 ! top of model
2452
2453 k = 1
2454
2455 ! add fallout terms to microphysical tendencies
2456
2457 tx5 = dumi(i,k)
2458 tx7 = pdel_inv(i,k) * tx1
2459 dumi(i,k) = tx5 / (one + fi(i,k)*tx7)
2460 tx6 = (dumi(i,k)-tx5) * oneodt
2461 qitend(i,k) = qitend(i,k) + tx6
2462 tx5 = dumni(i,k)
2463 dumni(i,k) = tx5 / (one + fni(i,k)*tx7)
2464 nitend(i,k) = nitend(i,k) + (dumni(i,k)-tx5) * oneodt
2465
2466 ! sedimentation tendency for output
2467 qisedten(i,k) = qisedten(i,k) + tx6
2468
2469 falouti(k) = fi(i,k) * dumi(i,k)
2470 faloutni(k) = fni(i,k) * dumni(i,k)
2471
2472 iflx(i,k+1) = iflx(i,k+1) + falouti(k) * tx3 ! Ice flux
2473
2474 do k = 2,nlev
2475
2476 ! for cloud liquid and ice, if cloud fraction increases with height
2477 ! then add flux from above to both vapor and cloud water of current level
2478 ! this means that flux entering clear portion of cell from above evaporates
2479 ! instantly
2480
2481 ! note: this is not an issue with precip, since we assume max overlap
2482
2483 if (icldm(i,k-1) > mincld) then
2484 dum1 = max(zero, min(one, icldm(i,k)/icldm(i,k-1)))
2485 else
2486 dum1 = one
2487 endif
2488
2489 tx5 = dumi(i,k)
2490 tx7 = pdel_inv(i,k) * tx1
2491 dum2 = tx7 * dum1
2492 dumi(i,k) = (tx5 + falouti(k-1)*dum2) / (one + fi(i,k)*tx7)
2493 tx6 = (dumi(i,k)-tx5) * oneodt
2494 ! add fallout terms to eulerian tendencies
2495 qitend(i,k) = qitend(i,k) + tx6
2496 tx5 = dumni(i,k)
2497 dumni(i,k) = (tx5 + faloutni(k-1)*dum2) / (one + fni(i,k)*tx7)
2498 nitend(i,k) = nitend(i,k) + (dumni(i,k)-tx5) * oneodt
2499
2500
2501 qisedten(i,k) = qisedten(i,k) + tx6 ! sedimentation tendency for output
2502
2503
2504 falouti(k) = fi(i,k) * dumi(i,k)
2505 faloutni(k) = fni(i,k) * dumni(i,k)
2506
2507 dum2 = (one-dum1) * falouti(k-1) * pdel_inv(i,k) * tx2
2508 qvlat(i,k) = qvlat(i,k) + dum2 ! add terms to evap/sub of cloud ice
2509 qisevap(i,k) = qisevap(i,k) + dum2 ! for output
2510
2511 tlat(i,k) = tlat(i,k) - dum2 * xxls
2512
2513 iflx(i,k+1) = iflx(i,k+1) + falouti(k) * tx3 ! Ice flux
2514 end do
2515
2516 ! units below are m/s
2517 ! sedimentation flux at surface is added to precip flux at surface
2518 ! to get total precip (cloud + precip water) rate
2519
2520 prect(i) = prect(i) + falouti(nlev) * (tx3*0.001_r8)
2521 preci(i) = preci(i) + falouti(nlev) * (tx3*0.001_r8)
2522
2523 end do
2524 end if
2525
2526! if (lprnt) write(0,*)' tlat4=',tlat(1,:)*deltat
2527 ! calculate number of split time steps to ensure courant stability criteria
2528 ! for sedimentation calculations
2529 !-------------------------------------------------------------------
2530 nstep = 1 + nint(max( maxval( fc(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
2531 maxval(fnc(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)
2532
2533 nstep = min(nstep, nstep_def)
2534
2535 ! loop over sedimentation sub-time step to ensure stability
2536 !==============================================================
2537 tx2 = one / nstep
2538 tx1 = tx2 * deltat
2539 tx3 = tx2 / g
2540
2541 do n = 1,nstep
2542
2543 ! top of model
2544 k = 1
2545
2546 tx5 = dumc(i,k)
2547 tx7 = pdel_inv(i,k) * tx1
2548 dumc(i,k) = tx5 / (one + fc(i,k)*tx7)
2549 tx6 = (dumc(i,k)-tx5) * oneodt
2550 qctend(i,k) = qctend(i,k) + tx6
2551 tx5 = dumnc(i,k)
2552 dumnc(i,k) = tx5 / (one + fnc(i,k)*tx7)
2553 nctend(i,k) = nctend(i,k) + (dumnc(i,k)-tx5) * oneodt
2554
2555
2556 ! sedimentation tendency for output
2557 qcsedten(i,k) = qcsedten(i,k) + tx6
2558
2559 faloutc(k) = fc(i,k) * dumc(i,k)
2560 faloutnc(k) = fnc(i,k) * dumnc(i,k)
2561
2562 lflx(i,k+1) = lflx(i,k+1) + faloutc(k) * tx3
2563 do k = 2,nlev
2564
2565 if (lcldm(i,k-1) > mincld) then
2566 dum1 = max(zero, min(one, lcldm(i,k)/lcldm(i,k-1)))
2567 else
2568 dum1 = one
2569 endif
2570
2571 tx5 = dumc(i,k)
2572 tx7 = pdel_inv(i,k) * tx1
2573 dum2 = tx7 * dum1
2574 dumc(i,k) = (tx5 + faloutc(k-1)*dum2) / (one + fc(i,k)*tx7)
2575 tx6 = (dumc(i,k)-tx5) * oneodt
2576 qctend(i,k) = qctend(i,k) + tx6
2577 tx5 = dumnc(i,k)
2578 dumnc(i,k) = (tx5 + faloutnc(k-1)*dum2) / (one + fnc(i,k)*tx7)
2579 nctend(i,k) = nctend(i,k) + (dumnc(i,k)-tx5) * oneodt
2580
2581
2582
2583 qcsedten(i,k) = qcsedten(i,k) + tx6 ! sedimentation tendency for output
2584
2585 faloutc(k) = fc(i,k) * dumc(i,k)
2586 faloutnc(k) = fnc(i,k) * dumnc(i,k)
2587
2588 dum2 = (one-dum1) * faloutc(k-1) * pdel_inv(i,k) * tx2
2589 qvlat(i,k) = qvlat(i,k) + dum2 ! add terms to to evap/sub of cloud water
2590 qcsevap(i,k) = qcsevap(i,k) + dum2 ! for output
2591
2592 tlat(i,k) = tlat(i,k) - dum2 * xxlv
2593
2594 lflx(i,k+1) = lflx(i,k+1) + faloutc(k) * tx3 ! Liquid condensate flux here
2595 end do
2596
2597 prect(i) = prect(i) + faloutc(nlev) * (tx3*0.001_r8)
2598
2599 end do
2600! if (lprnt) write(0,*)' tlat5=',tlat(1,:)*deltat
2601! if (lprnt) write(0,*)' maxval=',maxval( fr(i,nlb:nlev)*pdel_inv(i,nlb:nlev))&
2602! ,maxval(fnr(i,nlb:nlev)*pdel_inv(i,nlb:nlev))
2603
2604 ! calculate number of split time steps to ensure courant stability criteria
2605 ! for sedimentation calculations
2606 !-------------------------------------------------------------------
2607 nstep = 1 + nint(max( maxval( fr(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
2608 maxval(fnr(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)
2609
2610 nstep = min(nstep, nstep_def)
2611
2612 ! loop over sedimentation sub-time step to ensure stability
2613 !==============================================================
2614 tx2 = one / nstep
2615 tx1 = tx2 * deltat
2616 tx3 = tx2 / g
2617
2618! if(lprnt) then
2619! write(0,*)' nstep=',nstep,' tx1=',tx1,' tx2=',tx2,' tx3=',tx3,' qsmall=',qsmall
2620! write(0,*)' fr=',fr(i,:)
2621! write(0,*)' dumr=',dumr(i,:)
2622! endif
2623
2624 do n = 1,nstep
2625
2626 ! top of model
2627 k = 1
2628
2629 ! add fallout terms to microphysical tendencies
2630
2631 tx5 = dumr(i,k)
2632 tx7 = pdel_inv(i,k) * tx1
2633 dumr(i,k) = tx5 / (one + fr(i,k)*tx7)
2634 tx6 = (dumr(i,k)-tx5) * oneodt
2635 qrtend(i,k) = qrtend(i,k) + tx6
2636 tx5 = dumnr(i,k)
2637 dumnr(i,k) = tx5 / (one + fnr(i,k)*tx7)
2638 nrtend(i,k) = nrtend(i,k) + (dumnr(i,k)-tx5) * oneodt
2639
2640 ! sedimentation tendency for output
2641 qrsedten(i,k) = qrsedten(i,k) + tx6
2642
2643 faloutr(k) = fr(i,k) * dumr(i,k)
2644 faloutnr(k) = fnr(i,k) * dumnr(i,k)
2645
2646 rflx(i,k+1) = rflx(i,k+1) + faloutr(k) * tx3
2647
2648 do k = 2,nlev
2649
2650 tx5 = dumr(i,k)
2651 tx7 = pdel_inv(i,k) * tx1
2652 dumr(i,k) = (tx5 + faloutr(k-1)*tx7) / (one + fr(i,k)*tx7)
2653 tx6 = (dumr(i,k)-tx5) * oneodt
2654 qrtend(i,k) = qrtend(i,k) + tx6
2655 tx5 = dumnr(i,k)
2656 dumnr(i,k) = (tx5 + faloutnr(k-1)*tx7) / (one + fnr(i,k)*tx7)
2657 nrtend(i,k) = nrtend(i,k) + (dumnr(i,k)-tx5) * oneodt
2658
2659
2660 ! sedimentation tendency for output
2661 qrsedten(i,k) = qrsedten(i,k) + tx6 ! sedimentation tendency for output
2662
2663 faloutr(k) = fr(i,k) * dumr(i,k)
2664 faloutnr(k) = fnr(i,k) * dumnr(i,k)
2665
2666 rflx(i,k+1) = rflx(i,k+1) + faloutr(k) * tx3 ! Rain Flux
2667 end do
2668
2669 prect(i) = prect(i) + faloutr(nlev) * (tx3*0.001_r8)
2670
2671 end do
2672
2673! if (lprnt) write(0,*)' prectaftrain=',prect(i),' preci=',preci(i)
2674 ! calculate number of split time steps to ensure courant stability criteria
2675 ! for sedimentation calculations
2676 !-------------------------------------------------------------------
2677 nstep = 1 + nint(max( maxval( fs(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
2678 maxval(fns(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)
2679 nstep = min(nstep, nstep_def)
2680
2681 ! loop over sedimentation sub-time step to ensure stability
2682 !==============================================================
2683 tx2 = one / nstep
2684 tx1 = tx2 * deltat
2685 tx3 = tx2 / g
2686 do n = 1,nstep
2687
2688 ! top of model
2689 k = 1
2690
2691 ! add fallout terms to microphysical tendencies
2692
2693 tx5 = dums(i,k)
2694 tx7 = pdel_inv(i,k) * tx1
2695 dums(i,k) = tx5 / (one + fs(i,k)*tx7)
2696 tx6 = (dums(i,k)-tx5) * oneodt
2697 qstend(i,k) = qstend(i,k) + tx6
2698 tx5 = dumns(i,k)
2699 dumns(i,k) = tx5 / (one + fns(i,k)*tx7)
2700 nstend(i,k) = nstend(i,k) + (dumns(i,k)-tx5) * oneodt
2701
2702 ! sedimentation tendency for output
2703 qssedten(i,k) = qssedten(i,k) + tx6
2704
2705 falouts(k) = fs(i,k) * dums(i,k)
2706 faloutns(k) = fns(i,k) * dumns(i,k)
2707
2708 sflx(i,k+1) = sflx(i,k+1) + falouts(k) * tx3
2709
2710 do k = 2,nlev
2711
2712
2713 tx5 = dums(i,k)
2714 tx7 = pdel_inv(i,k) * tx1
2715 dums(i,k) = (tx5 + falouts(k-1)*tx7) / (one + fs(i,k)*tx7)
2716 tx6 = (dums(i,k)-tx5) * oneodt
2717 qstend(i,k) = qstend(i,k) + tx6
2718 tx5 = dumns(i,k)
2719 dumns(i,k) = (tx5 + faloutns(k-1)*tx7) / (one + fns(i,k)*tx7)
2720 nstend(i,k) = nstend(i,k) + (dumns(i,k)-tx5) * oneodt
2721
2722
2723 qssedten(i,k) = qssedten(i,k) + tx6 ! sedimentation tendency for output
2724
2725 falouts(k) = fs(i,k) * dums(i,k)
2726 faloutns(k) = fns(i,k) * dumns(i,k)
2727
2728 sflx(i,k+1) = sflx(i,k+1) + falouts(k) * tx3 ! Snow Flux
2729 end do !! k loop
2730
2731 prect(i) = prect(i) + falouts(nlev) * (tx3*0.001_r8)
2732 preci(i) = preci(i) + falouts(nlev) * (tx3*0.001_r8)
2733
2734 end do !! nstep loop
2735
2736 enddo ! end of i loop
2737 ! end sedimentation
2738
2739! if (lprnt) write(0,*)' prectaftsed=',prect(i),' preci=',preci(i)
2740
2741 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2742
2743 ! get new update for variables that includes sedimentation tendency
2744 ! note : here dum variables are grid-average, NOT in-cloud
2745
2746 do k=1,nlev
2747 do i=1,mgncol
2748 dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat, zero)
2749 dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat, zero)
2750 dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat, zero)
2751 dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat, zero)
2752
2753 dumr(i,k) = max(qr(i,k)+qrtend(i,k)*deltat, zero)
2754 dumnr(i,k) = max(nr(i,k)+nrtend(i,k)*deltat, zero)
2755 dums(i,k) = max(qs(i,k)+qstend(i,k)*deltat, zero)
2756 dumns(i,k) = max(ns(i,k)+nstend(i,k)*deltat, zero)
2757
2758 ! switch for specification of droplet and crystal number
2759 if (nccons) then
2760 dumnc(i,k) = ncnst*rhoinv(i,k)*lcldm(i,k)
2761 end if
2762
2763 ! switch for specification of cloud ice number
2764 if (nicons) then
2765 dumni(i,k) = ninst*rhoinv(i,k)*icldm(i,k)
2766 end if
2767
2768 if (dumc(i,k) < qsmall) dumnc(i,k) = zero
2769 if (dumi(i,k) < qsmall) dumni(i,k) = zero
2770 if (dumr(i,k) < qsmall) dumnr(i,k) = zero
2771 if (dums(i,k) < qsmall) dumns(i,k) = zero
2772
2773 enddo
2774
2775 enddo
2776
2777 ! calculate instantaneous processes (melting, homogeneous freezing)
2778 !====================================================================
2779
2780 ! melting of snow at +2 C
2781 do k=1,nlev
2782
2783 do i=1,mgncol
2784
2785 tx1 = t(i,k) + tlat(i,k)*(deltat/cpp) - snowmelt
2786 if (tx1 > zero) then
2787 if (dums(i,k) > zero) then
2788
2789 ! make sure melting snow doesn't reduce temperature below threshold
2790 dum = -(xlf/cpp) * dums(i,k)
2791 if (tx1+dum < zero) then
2792 dum = min(one, max(zero, -tx1/dum))
2793 else
2794 dum = one
2795 end if
2796
2797 tx1 = dum * oneodt
2798 qstend(i,k) = qstend(i,k) - tx1*dums(i,k)
2799 nstend(i,k) = nstend(i,k) - tx1*dumns(i,k)
2800 qrtend(i,k) = qrtend(i,k) + tx1*dums(i,k)
2801 nrtend(i,k) = nrtend(i,k) + tx1*dumns(i,k)
2802
2803 dum1 = - xlf * tx1 * dums(i,k)
2804 tlat(i,k) = tlat(i,k) + dum1
2805 meltsdttot(i,k) = meltsdttot(i,k) + dum1
2806 end if
2807 end if
2808 enddo
2809 enddo
2810 do k=1,nlev
2811 do i=1,mgncol
2812
2813 ! freezing of rain at -5 C
2814
2815 tx1 = t(i,k) + tlat(i,k) * (deltat/cpp) - rainfrze
2816 if (tx1 < zero) then
2817
2818 if (dumr(i,k) > zero) then
2819
2820 ! make sure freezing rain doesn't increase temperature above threshold
2821 dum = (xlf/cpp) * dumr(i,k)
2822 if (tx1+dum > zero) then
2823 dum = min(one, max(zero, -tx1/dum))
2824 else
2825 dum = one
2826 end if
2827 tx2 = dum * oneodt
2828 qrtend(i,k) = qrtend(i,k) - tx2 * dumr(i,k)
2829 nrtend(i,k) = nrtend(i,k) - tx2 * dumnr(i,k)
2830
2831 ! get mean size of rain = 1/lamr, add frozen rain to either snow or cloud ice
2832 ! depending on mean rain size
2833
2834 call size_dist_param_basic(mg_rain_props, dumr(i,k), dumnr(i,k), lamr(i,k))
2835
2836 if (lamr(i,k) < one/dcs) then
2837 qstend(i,k) = qstend(i,k) + tx2 * dumr(i,k)
2838 nstend(i,k) = nstend(i,k) + tx2 * dumnr(i,k)
2839 else
2840 qitend(i,k) = qitend(i,k) + tx2 * dumr(i,k)
2841 nitend(i,k) = nitend(i,k) + tx2 * dumnr(i,k)
2842 end if
2843 ! heating tendency
2844 dum1 = xlf*dum*dumr(i,k)*oneodt
2845 frzrdttot(i,k) = dum1 + frzrdttot(i,k)
2846 tlat(i,k) = dum1 + tlat(i,k)
2847
2848 end if
2849 end if
2850
2851 enddo
2852 enddo
2853 if (do_cldice) then
2854 do k=1,nlev
2855 do i=1,mgncol
2856 tx1 = t(i,k) + tlat(i,k) * (deltat/cpp) - tmelt
2857 if (tx1 > zero) then
2858 if (dumi(i,k) > zero) then
2859
2860 ! limit so that melting does not push temperature below freezing
2861 !-----------------------------------------------------------------
2862 dum = -dumi(i,k)*xlf/cpp
2863 if (tx1+dum < zero) then
2864 dum = min(one, max(zero, -tx1/dum))
2865 else
2866 dum = one
2867 end if
2868
2869 tx2 = dum * oneodt
2870 qctend(i,k) = qctend(i,k) + tx2*dumi(i,k)
2871
2872 ! for output
2873 melttot(i,k) = tx2*dumi(i,k)
2874
2875 ! assume melting ice produces droplet
2876 ! mean volume radius of 8 micron
2877
2878 nctend(i,k) = nctend(i,k) + three*tx2*dumi(i,k)/(four*pi*5.12e-16_r8*rhow)
2879
2880 qitend(i,k) = ((one-dum)*dumi(i,k)-qi(i,k)) * oneodt
2881 nitend(i,k) = ((one-dum)*dumni(i,k)-ni(i,k)) * oneodt
2882 tlat(i,k) = tlat(i,k) - xlf*tx2*dumi(i,k)
2883 end if
2884 end if
2885 enddo
2886 enddo
2887
2888! if (lprnt) write(0,*)' tlat6=',tlat(1,:)*deltat
2889! if (lprnt) write(0,*)' qitend=',qitend(1,nlev-10:nlev)*deltat
2890! if (lprnt) write(0,*)' qctend=',qctend(1,nlev-10:nlev)*deltat
2891 ! homogeneously freeze droplets at -40 C
2892 !-----------------------------------------------------------------
2893
2894 do k=1,nlev
2895 do i=1,mgncol
2896 tx1 = t(i,k) + tlat(i,k)*(deltat/cpp) - 233.15_r8
2897 if (tx1 < zero) then
2898 if (dumc(i,k) > zero) then
2899
2900 ! limit so that freezing does not push temperature above threshold
2901 dum = (xlf/cpp) * dumc(i,k)
2902 if (tx1+dum > zero) then
2903 dum = min(one, max(zero, -tx1/dum))
2904 else
2905 dum = one
2906 end if
2907
2908 tx2 = dum * oneodt * dumc(i,k)
2909 qitend(i,k) = tx2 + qitend(i,k)
2910 homotot(i,k) = tx2 ! for output
2911
2912 ! assume 25 micron mean volume radius of homogeneously frozen droplets
2913 ! consistent with size of detrained ice in stratiform.F90
2914
2915 nitend(i,k) = nitend(i,k) + tx2*(three/(four*pi*1.563e-14_r8* 500._r8))
2916 qctend(i,k) = ((one-dum)*dumc(i,k)-qc(i,k)) * oneodt
2917 nctend(i,k) = ((one-dum)*dumnc(i,k)-nc(i,k)) * oneodt
2918 tlat(i,k) = tlat(i,k) + xlf*tx2
2919 end if
2920 end if
2921 enddo
2922 enddo
2923 ! remove any excess over-saturation, which is possible due to non-linearity when adding
2924 ! together all microphysical processes
2925 !-----------------------------------------------------------------
2926 ! follow code similar to old CAM scheme
2927 do k=1,nlev
2928 do i=1,mgncol
2929
2930 qtmp = q(i,k) + qvlat(i,k) * deltat
2931 ttmp = t(i,k) + tlat(i,k) * (deltat/cpp)
2932
2933 ! use rhw to allow ice supersaturation
2934 !call qsat_water(ttmp, p(i,k), esn, qvn)
2935 esn = min(fpvsl(ttmp), p(i,k))
2936 qvn = epsqs*esn/(p(i,k)-omeps*esn) * qsfm(i,k)
2937
2938
2939 if (qtmp > qvn .and. qvn > 0 .and. allow_sed_supersat) then
2940 ! expression below is approximate since there may be ice deposition
2941 dum = (qtmp-qvn)/(one+xxlv_squared*qvn/(cpp*rv*ttmp*ttmp)) * oneodt
2942 ! add to output cme
2943 cmeout(i,k) = cmeout(i,k) + dum
2944 ! now add to tendencies, partition between liquid and ice based on temperature
2945 if (ttmp > 268.15_r8) then
2946 dum1 = zero
2947 ! now add to tendencies, partition between liquid and ice based on te
2948 !-------------------------------------------------------
2949 else if (ttmp < 238.15_r8) then
2950 dum1 = one
2951 else
2952 dum1 = (268.15_r8-ttmp)/30._r8
2953 end if
2954
2955 tx1 = xxls*dum1 + xxlv*(one-dum1)
2956 dum = (qtmp-qvn)/(one+tx1*tx1*qvn/(cpp*rv*ttmp*ttmp)) * oneodt
2957 tx2 = dum*(one-dum1)
2958 qctend(i,k) = qctend(i,k) + tx2
2959 qcrestot(i,k) = tx2 ! for output
2960 qitend(i,k) = qitend(i,k) + dum*dum1
2961 qirestot(i,k) = dum*dum1
2962 qvlat(i,k) = qvlat(i,k) - dum
2963 ! for output
2964 qvres(i,k) = -dum
2965 tlat(i,k) = tlat(i,k) + dum*tx1
2966 end if
2967 enddo
2968 enddo
2969 end if
2970
2971! if (lprnt) write(0,*)' tlat7=',tlat(1,:)*deltat
2972 ! calculate effective radius for pass to radiation code
2973 !=========================================================
2974 ! if no cloud water, default value is 10 micron for droplets,
2975 ! 25 micron for cloud ice
2976
2977 ! update cloud variables after instantaneous processes to get effective radius
2978 ! variables are in-cloud to calculate size dist parameters
2979 do k=1,nlev
2980 do i=1,mgncol
2981 if (lcldm(i,k) > mincld) then
2982 tx1 = one / lcldm(i,k)
2983 else
2984 tx1 = zero
2985 endif
2986 if (icldm(i,k) > mincld) then
2987 tx2 = one / icldm(i,k)
2988 else
2989 tx2 = zero
2990 endif
2991 if (precip_frac(i,k) > mincld) then
2992 tx3 = one / precip_frac(i,k)
2993 else
2994 tx3 = zero
2995 endif
2996 dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat, zero) * tx1
2997 dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat, zero) * tx2
2998 dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat, zero) * tx1
2999 dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat, zero) * tx2
3000
3001 dumr(i,k) = max(qr(i,k)+qrtend(i,k)*deltat, zero) * tx3
3002 dumnr(i,k) = max(nr(i,k)+nrtend(i,k)*deltat, zero) * tx3
3003 dums(i,k) = max(qs(i,k)+qstend(i,k)*deltat, zero) * tx3
3004 dumns(i,k) = max(ns(i,k)+nstend(i,k)*deltat, zero) * tx3
3005
3006 ! switch for specification of droplet and crystal number
3007 if (nccons) then
3008 dumnc(i,k) = ncnst * rhoinv(i,k)
3009 end if
3010
3011 ! switch for specification of cloud ice number
3012 if (nicons) then
3013 dumni(i,k) = ninst * rhoinv(i,k)
3014 end if
3015
3016 ! limit in-cloud mixing ratio to reasonable value of 5 g kg-1
3017! dumc(i,k) = min(dumc(i,k), 5.e-3_r8)
3018! dumi(i,k) = min(dumi(i,k), 5.e-3_r8)
3019 dumc(i,k) = min(dumc(i,k), 10.e-3_r8)
3020 dumi(i,k) = min(dumi(i,k), 10.e-3_r8)
3021 ! limit in-precip mixing ratios
3022 dumr(i,k) = min(dumr(i,k), 10.e-3_r8)
3023 dums(i,k) = min(dums(i,k), 10.e-3_r8)
3024 enddo
3025 enddo
3026 ! cloud ice effective radius
3027 !-----------------------------------------------------------------
3028
3029 if (do_cldice) then
3030 do k=1,nlev
3031 do i=1,mgncol
3032 if (dumi(i,k) >= qsmall) then
3033
3034 tx1 = dumni(i,k)
3035 call size_dist_param_basic(mg_ice_props, dumi(i,k), dumni(i,k), &
3036 lami(i,k), dumni0)
3037
3038 if (dumni(i,k) /= tx1) then
3039 ! adjust number conc if needed to keep mean size in reasonable range
3040 nitend(i,k) = (dumni(i,k)*icldm(i,k)-ni(i,k)) * oneodt
3041 end if
3042
3043 tx1 = one / lami(i,k)
3044! effi(i,k) = (1.5_r8*1.e6_r8) * tx1
3045 effi(i,k) = (three*1.e6_r8) * tx1
3046 sadice(i,k) = two*pi*(tx1*tx1*tx1)*dumni0*rho(i,k)*1.e-2_r8 ! m2/m3 -> cm2/cm3
3047
3048 else
3049! effi(i,k) = 25._r8
3050 effi(i,k) = 50._r8
3051 sadice(i,k) = zero
3052 end if
3053
3054 ! ice effective diameter for david mitchell's optics
3055 deffi(i,k) = effi(i,k) * (rhoi+rhoi)/rhows
3056 enddo
3057 enddo
3058 !else
3059 !do k=1,nlev
3060 !do i=1,mgncol
3061 ! NOTE: If CARMA is doing the ice microphysics, then the ice effective
3062 ! radius has already been determined from the size distribution.
3063 !effi(i,k) = re_ice(i,k) * 1.e6_r8 ! m -> um
3064 !deffi(i,k)=effi(i,k) * 2._r8
3065 !sadice(i,k) = 4._r8*pi*(effi(i,k)**2)*ni(i,k)*rho(i,k)*1e-2_r8
3066 !enddo
3067 !enddo
3068 end if
3069
3070 ! cloud droplet effective radius
3071 !-----------------------------------------------------------------
3072 do k=1,nlev
3073 do i=1,mgncol
3074 if (dumc(i,k) >= qsmall) then
3075
3076
3077 ! switch for specification of droplet and crystal number
3078 if (nccons) then
3079 ! make sure nc is consistence with the constant N by adjusting tendency, need
3080 ! to multiply by cloud fraction
3081 ! note that nctend may be further adjusted below if mean droplet size is
3082 ! out of bounds
3083
3084 nctend(i,k) = (ncnst*rhoinv(i,k)*lcldm(i,k)-nc(i,k)) * oneodt
3085
3086 end if
3087
3088 dum = dumnc(i,k)
3089
3090 call size_dist_param_liq(mg_liq_props, dumc(i,k), dumnc(i,k), rho(i,k), &
3091 pgam(i,k), lamc(i,k))
3092
3093 if (dum /= dumnc(i,k)) then
3094 ! adjust number conc if needed to keep mean size in reasonable range
3095 nctend(i,k) = (dumnc(i,k)*lcldm(i,k)-nc(i,k)) * oneodt
3096 end if
3097
3098 effc(i,k) = (half*1.e6_r8) * (pgam(i,k)+three) / lamc(i,k)
3099 !assign output fields for shape here
3100 lamcrad(i,k) = lamc(i,k)
3101 pgamrad(i,k) = pgam(i,k)
3102
3103
3104 ! recalculate effective radius for constant number, in order to separate
3105 ! first and second indirect effects
3106 !======================================
3107 ! assume constant number of 10^8 kg-1
3108
3109 dumnc(i,k) = 1.e8_r8
3110
3111 ! Pass in "false" adjust flag to prevent number from being changed within
3112 ! size distribution subroutine.
3113 call size_dist_param_liq(mg_liq_props, dumc(i,k), dumnc(i,k), rho(i,k), &
3114 pgam(i,k), lamc(i,k))
3115
3116 effc_fn(i,k) = (half*1.e6_r8) * (pgam(i,k)+three)/lamc(i,k)
3117
3118 else
3119 effc(i,k) = ten
3120 lamcrad(i,k) = zero
3121 pgamrad(i,k) = zero
3122 effc_fn(i,k) = ten
3123 end if
3124 enddo
3125 enddo
3126 ! recalculate 'final' rain size distribution parameters
3127 ! to ensure that rain size is in bounds, adjust rain number if needed
3128 do k=1,nlev
3129 do i=1,mgncol
3130
3131 if (dumr(i,k) >= qsmall) then
3132
3133 dum = dumnr(i,k)
3134
3135 call size_dist_param_basic(mg_rain_props, dumr(i,k), dumnr(i,k), lamr(i,k))
3136
3137 if (dum /= dumnr(i,k)) then
3138 ! adjust number conc if needed to keep mean size in reasonable range
3139 nrtend(i,k) = (dumnr(i,k)*precip_frac(i,k)-nr(i,k)) *oneodt
3140 end if
3141
3142 end if
3143 enddo
3144 enddo
3145 ! recalculate 'final' snow size distribution parameters
3146 ! to ensure that snow size is in bounds, adjust snow number if needed
3147 do k=1,nlev
3148 do i=1,mgncol
3149 if (dums(i,k) >= qsmall) then
3150
3151 dum = dumns(i,k)
3152
3153 call size_dist_param_basic(mg_snow_props, dums(i,k), dumns(i,k), &
3154 lams(i,k), n0=dumns0)
3155
3156 if (dum /= dumns(i,k)) then
3157 ! adjust number conc if needed to keep mean size in reasonable range
3158 nstend(i,k) = (dumns(i,k)*precip_frac(i,k)-ns(i,k)) * oneodt
3159 end if
3160
3161 tx1 = (two*pi*1.e-2_r8) / (lams(i,k)*lams(i,k)*lams(i,k))
3162 sadsnow(i,k) = tx1*dumns0*rho(i,k) ! m2/m3 -> cm2/cm3
3163
3164 end if
3165
3166
3167 end do ! vertical k loop
3168 enddo
3169 do k=1,nlev
3170 do i=1,mgncol
3171 ! if updated q (after microphysics) is zero, then ensure updated n is also zero
3172 !=================================================================================
3173 if (qc(i,k)+qctend(i,k)*deltat < qsmall) nctend(i,k) = -nc(i,k) * oneodt
3174 if (do_cldice .and. qi(i,k)+qitend(i,k)*deltat < qsmall) nitend(i,k) = -ni(i,k) * oneodt
3175 if (qr(i,k)+qrtend(i,k)*deltat < qsmall) nrtend(i,k) = -nr(i,k) * oneodt
3176 if (qs(i,k)+qstend(i,k)*deltat < qsmall) nstend(i,k) = -ns(i,k) * oneodt
3177
3178 end do
3179
3180 end do
3181
3182 ! DO STUFF FOR OUTPUT:
3183 !==================================================
3184
3185 do k=1,nlev
3186 do i=1,mgncol
3187
3188 ! qc and qi are only used for output calculations past here,
3189 ! so add qctend and qitend back in one more time
3190 qc(i,k) = qc(i,k) + qctend(i,k)*deltat
3191 qi(i,k) = qi(i,k) + qitend(i,k)*deltat
3192
3193 ! averaging for snow and rain number and diameter
3194 !--------------------------------------------------
3195
3196 ! drout2/dsout2:
3197 ! diameter of rain and snow
3198 ! dsout:
3199 ! scaled diameter of snow (passed to radiation in CAM)
3200 ! reff_rain/reff_snow:
3201 ! calculate effective radius of rain and snow in microns for COSP using Eq. 9 of COSP v1.3 manual
3202
3203 if (qrout(i,k) > 1.e-7_r8 .and. nrout(i,k) > zero) then
3204 qrout2(i,k) = qrout(i,k) * precip_frac(i,k)
3205 nrout2(i,k) = nrout(i,k) * precip_frac(i,k)
3206 ! The avg_diameter call does the actual calculation; other diameter
3207 ! outputs are just drout2 times constants.
3208 drout2(i,k) = avg_diameter(qrout(i,k), nrout(i,k), rho(i,k), rhow)
3209 freqr(i,k) = precip_frac(i,k)
3210
3211 reff_rain(i,k) = (1.e6_r8*three) * drout2(i,k)
3212 else
3213 qrout2(i,k) = zero
3214 nrout2(i,k) = zero
3215 drout2(i,k) = zero
3216 freqr(i,k) = zero
3217 reff_rain(i,k) = zero
3218 endif
3219
3220 if (qsout(i,k) > 1.e-7_r8 .and. nsout(i,k) > zero) then
3221 qsout2(i,k) = qsout(i,k) * precip_frac(i,k)
3222 nsout2(i,k) = nsout(i,k) * precip_frac(i,k)
3223 ! The avg_diameter call does the actual calculation; other diameter
3224 ! outputs are just dsout2 times constants.
3225 dsout2(i,k) = avg_diameter(qsout(i,k), nsout(i,k), rho(i,k), rhosn)
3226 freqs(i,k) = precip_frac(i,k)
3227
3228 dsout(i,k) = three*rhosn/rhows*dsout2(i,k)
3229
3230 reff_snow(i,k) = (1.e6_r8*three) * dsout2(i,k)
3231 else
3232 dsout(i,k) = zero
3233 qsout2(i,k) = zero
3234 nsout2(i,k) = zero
3235 dsout2(i,k) = zero
3236 freqs(i,k) = zero
3237 reff_snow(i,k) = zero
3238 endif
3239
3240 enddo
3241 enddo
3242
3243 ! analytic radar reflectivity
3244 !--------------------------------------------------
3245 ! formulas from Matthew Shupe, NOAA/CERES
3246 ! *****note: radar reflectivity is local (in-precip average)
3247 ! units of mm^6/m^3
3248
3249 do k=1,nlev
3250 do i = 1,mgncol
3251 if (qc(i,k) >= qsmall .and. (nc(i,k)+nctend(i,k)*deltat) > ten) then
3252 tx1 = rho(i,k) / lcldm(i,k)
3253 tx2 = 1000._r8 * qc(i,k) * tx1
3254 dum = tx2 * tx2 * lcldm(i,k) &
3255 /(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)*tx1*1.e-6_r8*precip_frac(i,k))
3256! dum = (qc(i,k)/lcldm(i,k)*rho(i,k)*1000._r8)**2 &
3257! /(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)/1.e6_r8)*lcldm(i,k)/precip_frac(i,k)
3258 else
3259 dum = zero
3260 end if
3261 if (qi(i,k) >= qsmall) then
3262! dum1 = (qi(i,k)*rho(i,k)/icldm(i,k)*1000._r8/0.1_r8)**(one/0.63_r8)*icldm(i,k)/precip_frac(i,k)
3263 dum1 = (qi(i,k)*rho(i,k)/icldm(i,k)*10000._r8)**(one/0.63_r8)*icldm(i,k)/precip_frac(i,k)
3264 else
3265 dum1 = zero
3266 end if
3267
3268 if (qsout(i,k) >= qsmall) then
3269! dum1 = dum1 + (qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)**(one/0.63_r8)
3270 dum1 = dum1 + (qsout(i,k)*rho(i,k)*10000._r8)**(one/0.63_r8)
3271 end if
3272
3273 refl(i,k) = dum + dum1
3274
3275 ! add rain rate, but for 37 GHz formulation instead of 94 GHz
3276 ! formula approximated from data of Matrasov (2007)
3277 ! rainrt is the rain rate in mm/hr
3278 ! reflectivity (dum) is in DBz
3279
3280 if (rainrt(i,k) >= 0.001_r8) then
3281 dum = rainrt(i,k) * rainrt(i,k)
3282 dum = log10(dum*dum*dum) + 16._r8
3283
3284 ! convert from DBz to mm^6/m^3
3285
3286 dum = ten**(dum/ten)
3287 else
3288 ! don't include rain rate in R calculation for values less than 0.001 mm/hr
3289 dum = zero
3290 end if
3291
3292 ! add to refl
3293
3294 refl(i,k) = refl(i,k) + dum
3295
3296 !output reflectivity in Z.
3297 areflz(i,k) = refl(i,k) * precip_frac(i,k)
3298
3299 ! convert back to DBz
3300
3301 if (refl(i,k) > minrefl) then
3302 refl(i,k) = ten*log10(refl(i,k))
3303 else
3304 refl(i,k) = -9999._r8
3305 end if
3306
3307 !set averaging flag
3308 if (refl(i,k) > mindbz) then
3309 arefl(i,k) = refl(i,k) * precip_frac(i,k)
3310 frefl(i,k) = precip_frac(i,k)
3311 else
3312 arefl(i,k) = zero
3313 areflz(i,k) = zero
3314 frefl(i,k) = zero
3315 end if
3316
3317 ! bound cloudsat reflectivity
3318
3319 csrfl(i,k) = min(csmax,refl(i,k))
3320
3321 !set averaging flag
3322 if (csrfl(i,k) > csmin) then
3323 acsrfl(i,k) = refl(i,k) * precip_frac(i,k)
3324 fcsrfl(i,k) = precip_frac(i,k)
3325 else
3326 acsrfl(i,k) = zero
3327 fcsrfl(i,k) = zero
3328 end if
3329
3330 end do
3331 end do
3332
3333 do k=1,nlev
3334 do i = 1,mgncol
3335 !redefine fice here....
3336 tx2 = qsout(i,k) + qi(i,k)
3337 tx1 = tx2 + qrout(i,k) + qc(i,k)
3338 if ( tx2 > qsmall .and. tx1 > qsmall) then
3339 nfice(i,k) = min(tx2/tx1, one)
3340 else
3341 nfice(i,k) = zero
3342 endif
3343 enddo
3344 enddo
3345
3346end subroutine micro_mg_tend
3347
3348!========================================================================
3349!OUTPUT CALCULATIONS
3350!========================================================================
3351
3354subroutine calc_rercld(lamr, n0r, lamc, pgam, qric, qcic, ncic, rercld, mgncol,nlev)
3355 integer, intent(in) :: mgncol, nlev
3356 real(r8), dimension(mgncol,nlev), intent(in) :: lamr ! rain size parameter (slope)
3357 real(r8), dimension(mgncol,nlev), intent(in) :: n0r ! rain size parameter (intercept)
3358 real(r8), dimension(mgncol,nlev), intent(in) :: lamc ! size distribution parameter (slope)
3359 real(r8), dimension(mgncol,nlev), intent(in) :: pgam ! droplet size parameter
3360 real(r8), dimension(mgncol,nlev), intent(in) :: qric ! in-cloud rain mass mixing ratio
3361 real(r8), dimension(mgncol,nlev), intent(in) :: qcic ! in-cloud cloud liquid
3362 real(r8), dimension(mgncol,nlev), intent(in) :: ncic ! in-cloud droplet number concentration
3363
3364 real(r8), dimension(mgncol,nlev), intent(inout) :: rercld ! effective radius calculation for rain + cloud
3365
3366 ! combined size of precip & cloud drops
3367 real(r8) :: Atmp
3368
3369 integer :: i, k
3370
3371 do k=1,nlev
3372 do i=1,mgncol
3373 ! Rain drops
3374 if (lamr(i,k) > zero) then
3375 atmp = n0r(i,k) * (half*pi) / (lamr(i,k)*lamr(i,k)*lamr(i,k))
3376 else
3377 atmp = zero
3378 end if
3379
3380 ! Add cloud drops
3381 if (lamc(i,k) > zero) then
3382 atmp = atmp + ncic(i,k) * pi * rising_factorial(pgam(i,k)+one, 2) &
3383 / (four*lamc(i,k)*lamc(i,k))
3384 end if
3385
3386 if (atmp > zero) then
3387 rercld(i,k) = rercld(i,k) + three *(qric(i,k) + qcic(i,k)) / (four * rhow * atmp)
3388 end if
3389 enddo
3390 enddo
3391end subroutine calc_rercld
3392
3393!========================================================================
3394
3395end module micro_mg2_0
subroutine calc_rercld(lamr, n0r, lamc, pgam, qric, qcic, ncic, rercld, mgncol, nlev)
This subroutine.
subroutine, public micro_mg_init(kind, gravit, rair, rh2o, cpair, eps, tmelt_in, latvap, latice, rhmini_in, micro_mg_dcs, ts_auto, mg_qcvar, microp_uniform_in, do_cldice_in, use_hetfrz_classnuc_in, micro_mg_precip_frac_method_in, micro_mg_berg_eff_factor_in, allow_sed_supersat_in, do_sb_physics_in, do_ice_gmao_in, do_liq_liu_in, nccons_in, nicons_in, ncnst_in, ninst_in)
This subroutine calculates.
subroutine, public ice_autoconversion(t, qiic, lami, n0i, dcs, ac_time, prci, nprci, mgncol)
Autoconversion of cloud ice to snow similar to Ferrier (1994)
subroutine, public liu_liq_autoconversion(pgam, qc, nc, qr, rho, relvar, au, nprc, nprc1, mgncol)
Anning Cheng 10/5/2017 add Liu et al. autoconversion.
subroutine, public bergeron_process_snow(t, rho, dv, mu, sc, qvl, qvi, asn, qcic, qsic, lams, n0s, bergs, mgncol)
bergeron process - evaporation of droplets and deposition onto snow
subroutine, public secondary_ice_production(t, psacws, msacwi, nsacwi, mgncol)
add secondary ice production due to accretion of droplets by snow
subroutine, public sb2001v2_accre_cld_water_rain(qc, nc, qr, rho, relvar, pra, npra, mgncol)
subroutine, public evaporate_sublimate_precip(t, rho, dv, mu, sc, q, qvl, qvi, lcldm, precip_frac, arn, asn, qcic, qiic, qric, qsic, lamr, n0r, lams, n0s, pre, prds, am_evp_st, mgncol)
calculate evaporation/sublimation of rain and snow
subroutine, public contact_freezing(microp_uniform, t, p, rndst, nacon, pgam, lamc, qcic, ncic, relvar, mnucct, nnucct, mgncol, mdust)
contact freezing (-40<T<-3 C) (Young, 1974) with hooks into simulated dust dust size and number in mu...
subroutine, public kk2000_liq_autoconversion(microp_uniform, qcic, ncic, rho, relvar, prc, nprc, nprc1, mgncol)
autoconversion of cloud liquid water to rain formula from Khrouditnov and Kogan (2000),...
subroutine, public self_collection_rain(rho, qric, nric, nragg, mgncol)
Self-collection of rain drops from Beheng(1994)
subroutine, public micro_mg_utils_init(kind, rair, rh2o, cpair, tmelt_in, latvap, latice, dcs)
Initialize module variables.
subroutine, public accrete_cloud_ice_snow(t, rho, asn, qiic, niic, qsic, lams, n0s, prai, nprai, mgncol)
Accretion of cloud ice by snow.
subroutine, public ice_deposition_sublimation(t, qv, qi, ni, icldm, rho, dv, qvl, qvi, berg, vap_dep, ice_sublim, mgncol)
Initial ice deposition and sublimation loop. Run before the main loop This subroutine written by Pete...
subroutine, public accrete_cloud_water_snow(t, rho, asn, uns, mu, qcic, ncic, qsic, pgam, lamc, lams, n0s, psacws, npsacws, mgncol)
accretion of cloud droplets onto snow/graupel
subroutine, public heterogeneous_rain_freezing(t, qric, nric, lamr, mnuccr, nnuccr, mgncol)
heterogeneous freezing of rain drops
subroutine, public snow_self_aggregation(t, rho, asn, rhosn, qsic, nsic, nsagg, mgncol)
snow self-aggregation from passarelli, 1978, used by reisner, 1998
real(r8) elemental function, public avg_diameter(q, n, rho_air, rho_sub)
Finds the average diameter of particles given their density, and mass/number concentrations in the ai...
subroutine, public sb2001v2_liq_autoconversion(pgam, qc, nc, qr, rho, relvar, au, nprc, nprc1, mgncol)
This subroutine.
subroutine, public gmao_ice_autoconversion(t, qiic, niic, lami, n0i, dcs, ac_time, prci, nprci, mgncol)
GMAO ice autoconversion.
subroutine, public accrete_cloud_water_rain(microp_uniform, qric, qcic, ncic, relvar, accre_enhan, pra, npra, mgncol)
accretion of cloud liquid water by rain formula from Khrouditnov and Kogan (2000)
subroutine, public accrete_rain_snow(t, rho, umr, ums, unr, uns, qric, qsic, lamr, n0r, lams, n0s, pracs, npracs, mgncol)
accretion of rain water by snow
subroutine, public immersion_freezing(microp_uniform, t, pgam, lamc, qcic, ncic, relvar, mnuccc, nnuccc, mgncol)
immersion freezing (Bigg, 1953)
subroutine, public gestbl(tmn, tmx, trice, ip, epsil, latvap, latice, rh2o, cpair, tmeltx)
This subroutine builds saturation vapor pressure table for later procedure.
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
Definition sflx.f:3