CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
micro_mg3_0.F90
1
4
5!---------------------------------------------------------------------------------
6! Purpose:
70!
71!---------------------------------------------------------------------------------
72!Version 3.O based on micro_mg2_0.F90 and WRF3.8.1 module_mp_morr_two_moment.F
73!---------------------------------------------------------------------------------
74! Based on micro_mg (restructuring of former cldwat2m_micro)
75! Author: Andrew Gettelman, Hugh Morrison.
76! Contributions from: Xiaohong Liu and Steve Ghan
77! December 2005-May 2010
78! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
79! Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)
80! for questions contact Hugh Morrison, Andrew Gettelman
81! e-mail: morrison@ucar.edu, andrew@ucar.edu
82!---------------------------------------------------------------------------------
83! Code comments added by HM, 093011
84! General code structure:
85!
86! Code is divided into two main subroutines:
87! subroutine micro_mg_init --> initializes microphysics routine, should be called
88! once at start of simulation
89! subroutine micro_mg_tend --> main microphysics routine to be called each time step
90! this also calls several smaller subroutines to calculate
91! microphysical processes and other utilities
92!
93! List of external functions:
94! qsat_water --> for calculating saturation vapor pressure with respect to liquid water
95! qsat_ice --> for calculating saturation vapor pressure with respect to ice
96! gamma --> standard mathematical gamma function
97! .........................................................................
98! List of inputs through use statement in fortran90:
99! Variable Name Description Units
100! .........................................................................
101! gravit acceleration due to gravity m s-2
102! rair dry air gas constant for air J kg-1 K-1
103! tmelt temperature of melting point for water K
104! cpair specific heat at constant pressure for dry air J kg-1 K-1
105! rh2o gas constant for water vapor J kg-1 K-1
106! latvap latent heat of vaporization J kg-1
107! latice latent heat of fusion J kg-1
108! qsat_water external function for calculating liquid water
109! saturation vapor pressure/humidity -
110! qsat_ice external function for calculating ice
111! saturation vapor pressure/humidity pa
112! rhmini relative humidity threshold parameter for
113! nucleating ice -
114! .........................................................................
115! NOTE: List of all inputs/outputs passed through the call/subroutine statement
116! for micro_mg_tend is given below at the start of subroutine micro_mg_tend.
117!---------------------------------------------------------------------------------
118
119! Procedures required:
120! 1) An implementation of the gamma function (if not intrinsic).
121! 2) saturation vapor pressure and specific humidity over water
122! 3) svp over ice
123
125
126use machine, only : r8 => kind_phys
127use funcphys, only : fpvsl, fpvsi
128
129!use wv_sat_methods, only: &
130! qsat_water => wv_sat_qsat_water, &
131! qsat_ice => wv_sat_qsat_ice
132
133! Parameters from the utilities module.
134use micro_mg_utils, only : pi, omsm, qsmall, mincld, rhosn, rhoi, &
135 rhow, rhows, ac, bc, ai, bi, &
136 aj, bj, ar, br, as, bs, &
137!++ag
138 ag, bg, ah, bh, rhog, rhoh, &
139!--ag
141
142implicit none
143private
144save
145
146public :: micro_mg_init, micro_mg_tend, qcvar
147
148! Switches for specification rather than prediction of droplet and crystal number
149! note: number will be adjusted as needed to keep mean size within bounds,
150! even when specified droplet or ice number is used
151!
152! If constant cloud ice number is set (nicons = .true.),
153! then all microphysical processes except mass transfer due to ice nucleation
154! (mnuccd) are based on the fixed cloud ice number. Calculation of
155! mnuccd follows from the prognosed ice crystal number ni.
156
157logical :: nccons
158logical :: nicons
159!++ag kt
160logical :: ngcons
161!--ag kt
162
163! specified ice and droplet number concentrations
164! note: these are local in-cloud values, not grid-mean
165real(r8) :: ncnst
166real(r8) :: ninst
167!++ag kt
168real(r8) :: ngnst
169!--ag kt
170
171!=========================================================
172! Private module parameters
173!=========================================================
174
176real(r8), parameter :: csmin = -30._r8
177real(r8), parameter :: csmax = 26._r8
178real(r8), parameter :: mindbz = -99._r8
179real(r8), parameter :: minrefl = 1.26e-10_r8 ! minrefl = 10._r8**(mindbz/10._r8)
180
181! autoconversion size threshold for cloud ice to snow (m)
182real(r8) :: dcs, ts_au, ts_au_min, qcvar
183
184! minimum mass of new crystal due to freezing of cloud droplets done
185! externally (kg)
186real(r8), parameter :: mi0l_min = 4._r8/3._r8*pi*rhow*(4.e-6_r8)**3
187
188! Ice number sublimation parameter. Assume some decrease in ice number with sublimation if non-zero. Else, no decrease in number with sublimation.
189real(r8), parameter :: sublim_factor = 0.0_r8 !number sublimation factor.
190
191real(r8), parameter :: zero=0.0_r8, one=1.0_r8, two=2.0_r8, three=3.0_r8, &
192 four=4.0_r8, five=5.0_r8, six=6._r8, half=0.5_r8, &
193 ten=10.0_r8, forty=40.0_r8, oneo6=one/six
194
195!=========================================================
196! Constants set in initialization
197!=========================================================
198
199! Set using arguments to micro_mg_init
200real(r8) :: g
201real(r8) :: r
202real(r8) :: rv
203real(r8) :: cpp
204real(r8) :: tmelt
205
206! latent heats of:
207real(r8) :: xxlv
208real(r8) :: xlf
209real(r8) :: xxls
210
211real(r8) :: rhmini
212
213! flags
214logical :: microp_uniform, do_cldice, use_hetfrz_classnuc, &
215!++ag
216 do_hail, do_graupel
217!--ag
218
219real(r8) :: rhosu
220
221real(r8) :: icenuct
222
223real(r8) :: snowmelt
224real(r8) :: rainfrze
225
226real(r8) :: rhogtmp
227real(r8) :: agtmp
228real(r8) :: bgtmp
229
230! additional constants to help speed up code
231real(r8) :: gamma_br_plus1, gamma_bs_plus1, gamma_bi_plus1, gamma_bj_plus1, gamma_bg_plus1
232real(r8) :: gamma_br_plus4, gamma_bs_plus4, gamma_bi_plus4, gamma_bj_plus4, gamma_bg_plus4
233real(r8) :: xxlv_squared, xxls_squared
234real(r8) :: omeps, epsqs
235
236character(len=16) :: micro_mg_precip_frac_method
237real(r8) :: micro_mg_berg_eff_factor
238
239logical :: allow_sed_supersat
240logical :: do_sb_physics
241logical :: do_ice_gmao
242logical :: do_liq_liu
243
244!===============================================================================
245contains
246!===============================================================================
247
251subroutine micro_mg_init( &
252 kind, gravit, rair, rh2o, cpair, eps, &
253 tmelt_in, latvap, latice, &
254 rhmini_in, micro_mg_dcs,ts_auto, mg_qcvar, &
255!++ag
256 micro_mg_do_hail_in, micro_mg_do_graupel_in, &
257!--ag
258 microp_uniform_in, do_cldice_in, use_hetfrz_classnuc_in, &
259 micro_mg_precip_frac_method_in, micro_mg_berg_eff_factor_in, &
260 allow_sed_supersat_in, do_sb_physics_in, &
261 do_ice_gmao_in, do_liq_liu_in, &
262 nccons_in, nicons_in, ncnst_in, ninst_in, ngcons_in, ngnst_in)
263! nccons_in, nicons_in, ncnst_in, ninst_in, ngcons_in, ngnst_in, errstring)
264
266 use wv_saturation, only : gestbl
267
268 !-----------------------------------------------------------------------
269 !
270 ! Purpose:
271 ! initialize constants for MG microphysics
272 !
273 ! Author: Andrew Gettelman Dec 2005
274 !
275 !-----------------------------------------------------------------------
276
277 integer, intent(in) :: kind ! Kind used for reals
278 real(r8), intent(in) :: gravit
279 real(r8), intent(in) :: rair
280 real(r8), intent(in) :: rh2o
281 real(r8), intent(in) :: cpair
282 real(r8), intent(in) :: eps
283! real(r8), intent(in) :: fv
284 real(r8), intent(in) :: tmelt_in ! Freezing point of water (K)
285 real(r8), intent(in) :: latvap
286 real(r8), intent(in) :: latice
287 real(r8), intent(in) :: rhmini_in ! Minimum rh for ice cloud fraction > 0.
288 real(r8), intent(in) :: micro_mg_dcs
289 real(r8), intent(in) :: ts_auto(2)
290 real(r8), intent(in) :: mg_qcvar
291
292!++ag
293!MG3 dense precipitating ice. Note, only 1 can be true, or both false.
294 logical, intent(in) :: micro_mg_do_graupel_in ! .true. = configure with graupel
295 ! .false. = no graupel (hail possible)
296 logical, intent(in) :: micro_mg_do_hail_in ! .true. = configure with hail
297 ! .false. = no hail (graupel possible)
298!--ag
299
300 logical, intent(in) :: microp_uniform_in ! .true. = configure uniform for sub-columns
301 ! .false. = use w/o sub-columns (standard)
302 logical, intent(in) :: do_cldice_in ! .true. = do all processes (standard)
303 ! .false. = skip all processes affecting cloud ice
304 logical, intent(in) :: use_hetfrz_classnuc_in ! use heterogeneous freezing
305
306 character(len=16),intent(in) :: micro_mg_precip_frac_method_in ! type of precipitation fraction method
307 real(r8), intent(in) :: micro_mg_berg_eff_factor_in ! berg efficiency factor
308 logical, intent(in) :: allow_sed_supersat_in ! allow supersaturated conditions after sedimentation loop
309 logical, intent(in) :: do_sb_physics_in ! do SB autoconversion and accretion physics
310 logical, intent(in) :: do_ice_gmao_in
311 logical, intent(in) :: do_liq_liu_in
312
313 logical, intent(in) :: nccons_in, nicons_in, ngcons_in
314 real(r8), intent(in) :: ncnst_in, ninst_in, ngnst_in
315 logical ip
316 real(r8):: tmn, tmx, trice
317
318
319! character(128), intent(out) :: errstring ! Output status (non-blank for error return)
320
321 !-----------------------------------------------------------------------
322
323 dcs = micro_mg_dcs * 1.0e-6_r8
324 ts_au_min = ts_auto(1)
325 ts_au = ts_auto(2)
326 qcvar = mg_qcvar
327
328 ! Initialize subordinate utilities module.
329 call micro_mg_utils_init(kind, rair, rh2o, cpair, tmelt_in, latvap, latice, &
330 dcs)
331! dcs, errstring)
332
333! if (trim(errstring) /= "") return
334
335 ! declarations for MG code (transforms variable names)
336
337 g = gravit ! gravity
338 r = rair ! dry air gas constant: note units(phys_constants are in J/K/kmol)
339! write(0,*)' in micro_mg_utils_init=',' r=',r,' rair=',rair,' rh2o=',rh2o
340 rv = rh2o ! water vapor gas constant
341 cpp = cpair ! specific heat of dry air
342 tmelt = tmelt_in
343 rhmini = rhmini_in
344 micro_mg_precip_frac_method = micro_mg_precip_frac_method_in
345 micro_mg_berg_eff_factor = micro_mg_berg_eff_factor_in
346 allow_sed_supersat = allow_sed_supersat_in
347 do_sb_physics = do_sb_physics_in
348 do_ice_gmao = do_ice_gmao_in
349 do_liq_liu = do_liq_liu_in
350
351 nccons = nccons_in
352 nicons = nicons_in
353 ncnst = ncnst_in
354 ninst = ninst_in
355!++ag
356 ngcons = ngcons_in
357 ngnst = ngnst_in
358!--ag
359
360 ! latent heats
361
362 xxlv = latvap ! latent heat vaporization
363 xlf = latice ! latent heat freezing
364 xxls = xxlv + xlf ! latent heat of sublimation
365
366 ! flags
367 microp_uniform = microp_uniform_in
368 do_cldice = do_cldice_in
369 use_hetfrz_classnuc = use_hetfrz_classnuc_in
370!++ag
371 do_hail = micro_mg_do_hail_in
372 do_graupel = micro_mg_do_graupel_in
373!
374 if (do_hail) then
375 agtmp = ah
376 bgtmp = bh
377 rhogtmp = rhoh
378 elseif (do_graupel) then
379 agtmp = ag
380 bgtmp = bg
381 rhogtmp = rhog
382 else
383 agtmp = zero
384 bgtmp = zero
385 endif
386!--ag
387
388 ! typical air density at 850 mb
389
390 rhosu = 85000._r8 / (rair * tmelt)
391
392 ! Maximum temperature at which snow is allowed to exist
393 snowmelt = tmelt + two
394 ! Minimum temperature at which rain is allowed to exist
395 rainfrze = tmelt - forty
396
397 ! Ice nucleation temperature
398 icenuct = tmelt - five
399
400 ! Define constants to help speed up code (this limits calls to gamma function)
401 gamma_br_plus1 = gamma(br+one)
402 gamma_br_plus4 = gamma(br+four)
403 gamma_bs_plus1 = gamma(bs+one)
404 gamma_bs_plus4 = gamma(bs+four)
405 gamma_bi_plus1 = gamma(bi+one)
406 gamma_bi_plus4 = gamma(bi+four)
407 gamma_bj_plus1 = gamma(bj+one)
408 gamma_bj_plus4 = gamma(bj+four)
409!
410 gamma_bg_plus1 = gamma(bgtmp+one)
411 gamma_bg_plus4 = gamma(bgtmp+four)
412
413 xxlv_squared = xxlv * xxlv
414 xxls_squared = xxls * xxls
415 epsqs = eps
416 omeps = one - epsqs
417 tmn = 173.16_r8
418 tmx = 375.16_r8
419 trice = 35.00_r8
420 ip = .true.
422 call gestbl(tmn ,tmx ,trice ,ip ,epsqs , latvap ,latice ,rh2o , &
423 cpair ,tmelt_in )
424
425
426
427end subroutine micro_mg_init
428
429!===============================================================================
430!microphysics routine for each timestep goes here...
431
437subroutine micro_mg_tend ( &
438 mgncol, nlev, deltatin, &
439 t, q, &
440 qcn, qin, &
441 ncn, nin, &
442 qrn, qsn, &
443 nrn, nsn, &
444!++ag
445 qgr, ngr, &
446!--ag
447 relvar, accre_enhan_i, &
448 p, pdel, &
449 cldn, liqcldf, icecldf, qsatfac, &
450 qcsinksum_rate1ord, &
451 naai, npccnin, &
452 rndst, nacon, &
453 tlat, qvlat, &
454 qctend, qitend, &
455 nctend, nitend, &
456 qrtend, qstend, &
457 nrtend, nstend, &
458!++ag
459 qgtend, ngtend, &
460!--ag
461 effc, effc_fn, effi, &
462 sadice, sadsnow, &
463 prect, preci, &
464 nevapr, evapsnow, &
465 am_evp_st, &
466 prain, prodsnow, &
467 cmeout, deffi, &
468 pgamrad, lamcrad, &
469 qsout, dsout, &
470!++ag
471 qgout, ngout, dgout, &
472!--ag
473 lflx, iflx, &
474!++ag
475 gflx, &
476!--ag
477 rflx, sflx, qrout, &
478!++ag
479 reff_rain, reff_snow, reff_grau, &
480!--ag
481
482 qcsevap, qisevap, qvres, &
483 cmeitot, vtrmc, vtrmi, &
484 umr, ums, &
485!++ag
486 umg, qgsedten, &
487!--ag
488 qcsedten, qisedten, &
489 qrsedten, qssedten, &
490 pratot, prctot, &
491 mnuccctot, mnuccttot, msacwitot, &
492 psacwstot, bergstot, bergtot, &
493 melttot, homotot, &
494 qcrestot, prcitot, praitot, &
495!++ag
496 qirestot, mnuccrtot, mnuccritot, pracstot, &
497!--ag
498 meltsdttot, frzrdttot, mnuccdtot, &
499!++ag
500 pracgtot, psacwgtot, pgsacwtot, &
501 pgracstot, prdgtot, &
502 qmultgtot, qmultrgtot, psacrtot, &
503 npracgtot, nscngtot, ngracstot, &
504 nmultgtot, nmultrgtot, npsacwgtot, &
505!--ag
506 nrout, nsout, &
507 refl, arefl, areflz, &
508 frefl, csrfl, acsrfl, &
509 fcsrfl, rercld, &
510 ncai, ncal, &
511 qrout2, qsout2, &
512 nrout2, nsout2, &
513 drout2, dsout2, &
514!++ag
515 qgout2, ngout2, dgout2, freqg, &
516!--ag
517 freqs, freqr, &
518 nfice, qcrat, &
519 prer_evap, xlat, xlon, lprnt, iccn, nlball)
520
521 ! Constituent properties.
522 use micro_mg_utils, only: mg_liq_props, &
523 mg_ice_props, &
524 mg_rain_props, &
525!++ag
526 mg_graupel_props,&
527!--ag
528 mg_snow_props
529
530 ! Size calculation functions.
534
535 ! Microphysical processes.
554!++ag
562! graupel_sublimate_evap
563!--ag
566
567 !Authors: Hugh Morrison, Andrew Gettelman, NCAR, Peter Caldwell, LLNL
568 ! e-mail: morrison@ucar.edu, andrew@ucar.edu
569
570 ! input arguments
571 integer, intent(in) :: mgncol
572 integer, intent(in) :: nlev
573 integer, intent(in) :: nlball(mgncol)
574 real(r8), intent(in) :: xlat,xlon
575 real(r8), intent(in) :: deltatin
576 real(r8), intent(in) :: t(mgncol,nlev)
577 real(r8), intent(in) :: q(mgncol,nlev)
578
579 ! note: all input cloud variables are grid-averaged
580 real(r8), intent(in) :: qcn(mgncol,nlev)
581 real(r8), intent(in) :: qin(mgncol,nlev)
582 real(r8), intent(in) :: ncn(mgncol,nlev)
583 real(r8), intent(in) :: nin(mgncol,nlev)
584
585 real(r8), intent(in) :: qrn(mgncol,nlev)
586 real(r8), intent(in) :: qsn(mgncol,nlev)
587 real(r8), intent(in) :: nrn(mgncol,nlev)
588 real(r8), intent(in) :: nsn(mgncol,nlev)
589!++ag
590 real(r8), intent(in) :: qgr(mgncol,nlev)
591 real(r8), intent(in) :: ngr(mgncol,nlev)
592!--ag
593
594 real(r8) :: relvar(mgncol,nlev)
595 real(r8) :: accre_enhan(mgncol,nlev)
596! real(r8), intent(in) :: relvar_i !< cloud water relative variance (-)
597 real(r8), intent(in) :: accre_enhan_i
599
600 real(r8), intent(in) :: p(mgncol,nlev)
601 real(r8), intent(in) :: pdel(mgncol,nlev)
602
603 real(r8), intent(in) :: cldn(mgncol,nlev)
604 real(r8), intent(in) :: liqcldf(mgncol,nlev)
605 real(r8), intent(in) :: icecldf(mgncol,nlev)
606 real(r8), intent(in) :: qsatfac(mgncol,nlev)
607 logical, intent(in) :: lprnt
608 integer, intent(in) :: iccn
609
610
611 ! used for scavenging
612 ! Inputs for aerosol activation
613 real(r8), intent(inout) :: naai(mgncol,nlev)
614 real(r8), intent(in) :: npccnin(mgncol,nlev)
615! real(r8), intent(in) :: npccn(mgncol,nlev) !< ccn activated number tendency (from microp_aero_ts) (1/kg*s)
616 real(r8) :: npccn(mgncol,nlev)
617
618 ! Note that for these variables, the dust bin is assumed to be the last index.
619 ! (For example, in CAM, the last dimension is always size 4.)
620 real(r8), intent(in) :: rndst(mgncol,nlev,10)
621 real(r8), intent(in) :: nacon(mgncol,nlev,10)
622
623 ! output arguments
624
625 real(r8), intent(out) :: qcsinksum_rate1ord(mgncol,nlev)
627 real(r8), intent(out) :: tlat(mgncol,nlev)
628 real(r8), intent(out) :: qvlat(mgncol,nlev)
629 real(r8), intent(out) :: qctend(mgncol,nlev)
630 real(r8), intent(out) :: qitend(mgncol,nlev)
631 real(r8), intent(out) :: nctend(mgncol,nlev)
632 real(r8), intent(out) :: nitend(mgncol,nlev)
633
634 real(r8), intent(out) :: qrtend(mgncol,nlev)
635 real(r8), intent(out) :: qstend(mgncol,nlev)
636 real(r8), intent(out) :: nrtend(mgncol,nlev)
637 real(r8), intent(out) :: nstend(mgncol,nlev)
638!++ag
639 real(r8), intent(out) :: qgtend(mgncol,nlev)
640 real(r8), intent(out) :: ngtend(mgncol,nlev)
641!--ag
642 real(r8), intent(out) :: effc(mgncol,nlev)
643 real(r8), intent(out) :: effc_fn(mgncol,nlev)
644 real(r8), intent(out) :: effi(mgncol,nlev)
645 real(r8), intent(out) :: sadice(mgncol,nlev)
646 real(r8), intent(out) :: sadsnow(mgncol,nlev)
647 real(r8), intent(out) :: prect(mgncol)
648 real(r8), intent(out) :: preci(mgncol)
649 real(r8), intent(out) :: nevapr(mgncol,nlev)
650 real(r8), intent(out) :: evapsnow(mgncol,nlev)
651 real(r8), intent(out) :: am_evp_st(mgncol,nlev)
652 real(r8), intent(out) :: prain(mgncol,nlev)
653 real(r8), intent(out) :: prodsnow(mgncol,nlev)
654 real(r8), intent(out) :: cmeout(mgncol,nlev)
655 real(r8), intent(out) :: deffi(mgncol,nlev)
656 real(r8), intent(out) :: pgamrad(mgncol,nlev)
657 real(r8), intent(out) :: lamcrad(mgncol,nlev)
658 real(r8), intent(out) :: qsout(mgncol,nlev)
659 real(r8), intent(out) :: dsout(mgncol,nlev)
660 real(r8), intent(out) :: lflx(mgncol,2:nlev+1)
661 real(r8), intent(out) :: iflx(mgncol,2:nlev+1)
662 real(r8), intent(out) :: rflx(mgncol,2:nlev+1)
663 real(r8), intent(out) :: sflx(mgncol,2:nlev+1)
664!++ag
665 real(r8), intent(out) :: gflx(mgncol,2:nlev+1)
666!--ag
667 real(r8), intent(out) :: qrout(mgncol,nlev)
668 real(r8), intent(out) :: reff_rain(mgncol,nlev)
669 real(r8), intent(out) :: reff_snow(mgncol,nlev)
670!++ag
671 real(r8), intent(out) :: reff_grau(mgncol,nlev)
672!--ag
673 real(r8), intent(out) :: qcsevap(mgncol,nlev)
674 real(r8), intent(out) :: qisevap(mgncol,nlev)
675 real(r8), intent(out) :: qvres(mgncol,nlev)
676 real(r8), intent(out) :: cmeitot(mgncol,nlev)
677 real(r8), intent(out) :: vtrmc(mgncol,nlev)
678 real(r8), intent(out) :: vtrmi(mgncol,nlev)
679 real(r8), intent(out) :: umr(mgncol,nlev)
680 real(r8), intent(out) :: ums(mgncol,nlev)
681!++ag
682 real(r8), intent(out) :: umg(mgncol,nlev)
683 real(r8), intent(out) :: qgsedten(mgncol,nlev)
684!--ag
685
686 real(r8), intent(out) :: qcsedten(mgncol,nlev)
687 real(r8), intent(out) :: qisedten(mgncol,nlev)
688 real(r8), intent(out) :: qrsedten(mgncol,nlev)
689 real(r8), intent(out) :: qssedten(mgncol,nlev)
690
691 ! microphysical process rates for output (mixing ratio tendencies) (all have units of 1/s)
692 real(r8), intent(out) :: pratot(mgncol,nlev)
693 real(r8), intent(out) :: prctot(mgncol,nlev)
694 real(r8), intent(out) :: mnuccctot(mgncol,nlev)
695 real(r8), intent(out) :: mnuccttot(mgncol,nlev)
696 real(r8), intent(out) :: msacwitot(mgncol,nlev)
697 real(r8), intent(out) :: psacwstot(mgncol,nlev)
698 real(r8), intent(out) :: bergstot(mgncol,nlev)
699 real(r8), intent(out) :: bergtot(mgncol,nlev)
700 real(r8), intent(out) :: melttot(mgncol,nlev)
701 real(r8), intent(out) :: homotot(mgncol,nlev)
702 real(r8), intent(out) :: qcrestot(mgncol,nlev)
703 real(r8), intent(out) :: prcitot(mgncol,nlev)
704 real(r8), intent(out) :: praitot(mgncol,nlev)
705 real(r8), intent(out) :: qirestot(mgncol,nlev)
706 real(r8), intent(out) :: mnuccrtot(mgncol,nlev)
707 real(r8), intent(out) :: mnuccritot(mgncol,nlev)
708 real(r8), intent(out) :: pracstot(mgncol,nlev)
709 real(r8), intent(out) :: meltsdttot(mgncol,nlev)
710 real(r8), intent(out) :: frzrdttot(mgncol,nlev)
711 real(r8), intent(out) :: mnuccdtot(mgncol,nlev)
712!++ag Hail/Graupel Tendencies
713 real(r8), intent(out) :: pracgtot(mgncol,nlev)
714 real(r8), intent(out) :: psacwgtot(mgncol,nlev)
715 real(r8), intent(out) :: pgsacwtot(mgncol,nlev)
716 real(r8), intent(out) :: pgracstot(mgncol,nlev)
717 real(r8), intent(out) :: prdgtot(mgncol,nlev)
718! real(r8), intent(out) :: eprdgtot(mgncol,nlev) !< sub of graupel (precipf)
719 real(r8), intent(out) :: qmultgtot(mgncol,nlev)
720 real(r8), intent(out) :: qmultrgtot(mgncol,nlev)
721 real(r8), intent(out) :: psacrtot(mgncol,nlev)
722 real(r8), intent(out) :: npracgtot(mgncol,nlev)
723 real(r8), intent(out) :: nscngtot(mgncol,nlev)
724 real(r8), intent(out) :: ngracstot(mgncol,nlev)
725 real(r8), intent(out) :: nmultgtot(mgncol,nlev)
726 real(r8), intent(out) :: nmultrgtot(mgncol,nlev)
727 real(r8), intent(out) :: npsacwgtot(mgncol,nlev)
728!--ag
729 real(r8), intent(out) :: nrout(mgncol,nlev)
730 real(r8), intent(out) :: nsout(mgncol,nlev)
731 real(r8), intent(out) :: refl(mgncol,nlev)
732 real(r8), intent(out) :: arefl(mgncol,nlev)
733 real(r8), intent(out) :: areflz(mgncol,nlev)
734 real(r8), intent(out) :: frefl(mgncol,nlev)
735 real(r8), intent(out) :: csrfl(mgncol,nlev)
736 real(r8), intent(out) :: acsrfl(mgncol,nlev)
737 real(r8), intent(out) :: fcsrfl(mgncol,nlev)
738 real(r8), intent(out) :: rercld(mgncol,nlev)
739 real(r8), intent(out) :: ncai(mgncol,nlev)
740 real(r8), intent(out) :: ncal(mgncol,nlev)
741 real(r8), intent(out) :: qrout2(mgncol,nlev)
742 real(r8), intent(out) :: qsout2(mgncol,nlev)
743 real(r8), intent(out) :: nrout2(mgncol,nlev)
744 real(r8), intent(out) :: nsout2(mgncol,nlev)
745 real(r8), intent(out) :: drout2(mgncol,nlev)
746 real(r8), intent(out) :: dsout2(mgncol,nlev)
747 real(r8), intent(out) :: freqs(mgncol,nlev)
748 real(r8), intent(out) :: freqr(mgncol,nlev)
749 real(r8), intent(out) :: nfice(mgncol,nlev)
750 real(r8), intent(out) :: qcrat(mgncol,nlev)
751!++ag
752 real(r8), intent(out) :: qgout(mgncol,nlev)
753 real(r8), intent(out) :: dgout(mgncol,nlev)
754 real(r8), intent(out) :: ngout(mgncol,nlev)
755!Not sure if these are needed since graupel/hail is prognostic?
756 real(r8), intent(out) :: qgout2(mgncol,nlev)
757 real(r8), intent(out) :: ngout2(mgncol,nlev)
758 real(r8), intent(out) :: dgout2(mgncol,nlev)
759 real(r8), intent(out) :: freqg(mgncol,nlev)
760
761!--ag
762
763 real(r8), intent(out) :: prer_evap(mgncol,nlev)
764
765
766 ! Tendencies calculated by external schemes that can replace MG's native
767 ! process tendencies.
768
769 ! Used with CARMA cirrus microphysics
770 ! (or similar external microphysics model)
771 ! real(r8), intent(in) :: tnd_qsnow(:,:) !< snow mass tendency (kg/kg/s)
772 ! real(r8), intent(in) :: tnd_nsnow(:,:) !< snow number tendency (#/kg/s)
773 ! real(r8), intent(in) :: re_ice(:,:) !< ice effective radius (m)
774
775 ! From external ice nucleation.
776 !real(r8), intent(in) :: frzimm(:,:) !< Number tendency due to immersion freezing (1/cm3)
777 !real(r8), intent(in) :: frzcnt(:,:) !< Number tendency due to contact freezing (1/cm3)
778 !real(r8), intent(in) :: frzdep(:,:) !< Number tendency due to deposition nucleation (1/cm3)
779
780 ! local workspace
781 ! all units mks unless otherwise stated
782
783 ! local copies of input variables
784 real(r8) :: qc(mgncol,nlev)
785 real(r8) :: qi(mgncol,nlev)
786 real(r8) :: nc(mgncol,nlev)
787 real(r8) :: ni(mgncol,nlev)
788 real(r8) :: qr(mgncol,nlev)
789 real(r8) :: qs(mgncol,nlev)
790 real(r8) :: nr(mgncol,nlev)
791 real(r8) :: ns(mgncol,nlev)
792!++ag
793 real(r8) :: qg(mgncol,nlev)
794 real(r8) :: ng(mgncol,nlev)
795! real(r8) :: rhogtmp !< hail or graupel density (kg m-3)
796
797!--ag
798
799 ! general purpose variables
800 real(r8) :: deltat
801 real(r8) :: oneodt
802 real(r8) :: mtime
803
804 ! physical properties of the air at a given point
805 real(r8) :: rho(mgncol,nlev) ! density (kg m-3)
806 real(r8) :: rhoinv(mgncol,nlev) ! one / density (kg m-3)
807 real(r8) :: dv(mgncol,nlev) ! diffusivity of water vapor
808 real(r8) :: mu(mgncol,nlev) ! viscosity
809 real(r8) :: sc(mgncol,nlev) ! schmidt number
810 real(r8) :: rhof(mgncol,nlev) ! density correction factor for fallspeed
811
812 ! cloud fractions
813 real(r8) :: precip_frac(mgncol,nlev)! precip fraction assuming maximum overlap
814 real(r8) :: cldm(mgncol,nlev) ! cloud fraction
815 real(r8) :: icldm(mgncol,nlev) ! ice cloud fraction
816 real(r8) :: lcldm(mgncol,nlev) ! liq cloud fraction
817 real(r8) :: qsfm(mgncol,nlev) ! subgrid cloud water saturation scaling factor
818
819 ! mass mixing ratios
820 real(r8) :: qcic(mgncol,nlev) ! in-cloud cloud liquid
821 real(r8) :: qiic(mgncol,nlev) ! in-cloud cloud ice
822 real(r8) :: qsic(mgncol,nlev) ! in-precip snow
823 real(r8) :: qric(mgncol,nlev) ! in-precip rain
824!++ag
825 real(r8) :: qgic(mgncol,nlev) ! in-precip graupel/hail
826!++ag
827
828
829 ! number concentrations
830 real(r8) :: ncic(mgncol,nlev) ! in-cloud droplet
831 real(r8) :: niic(mgncol,nlev) ! in-cloud cloud ice
832 real(r8) :: nsic(mgncol,nlev) ! in-precip snow
833 real(r8) :: nric(mgncol,nlev) ! in-precip rain
834!++ag
835 real(r8) :: ngic(mgncol,nlev) ! in-precip graupel/hail
836!++ag
837
838 ! maximum allowed ni value
839 real(r8) :: nimax(mgncol,nlev)
840
841 ! Size distribution parameters for:
842 ! cloud ice
843 real(r8) :: lami(mgncol,nlev) ! slope
844 real(r8) :: n0i(mgncol,nlev) ! intercept
845 ! cloud liquid
846 real(r8) :: lamc(mgncol,nlev) ! slope
847 real(r8) :: pgam(mgncol,nlev) ! spectral width parameter
848 ! snow
849 real(r8) :: lams(mgncol,nlev) ! slope
850 real(r8) :: n0s(mgncol,nlev) ! intercept
851 ! rain
852 real(r8) :: lamr(mgncol,nlev) ! slope
853 real(r8) :: n0r(mgncol,nlev) ! intercept
854!++ag
855 ! graupel/hail
856 real(r8) :: lamg(mgncol,nlev) ! slope
857 real(r8) :: n0g(mgncol,nlev) ! intercept
858! real(r8) :: bgtmp ! tmp fall speed parameter
859!--ag
860
861 ! Rates/tendencies due to:
862
863 ! Instantaneous snow melting
864 real(r8) :: minstsm(mgncol,nlev) ! mass mixing ratio
865 real(r8) :: ninstsm(mgncol,nlev) ! number concentration
866!++ag
867 ! Instantaneous graupel melting
868 real(r8) :: minstgm(mgncol,nlev) ! mass mixing ratio
869 real(r8) :: ninstgm(mgncol,nlev) ! number concentration
870!--ag
871
872 ! Instantaneous rain freezing
873 real(r8) :: minstrf(mgncol,nlev) ! mass mixing ratio
874 real(r8) :: ninstrf(mgncol,nlev) ! number concentration
875
876 ! deposition of cloud ice
877 real(r8) :: vap_dep(mgncol,nlev) ! deposition from vapor to ice PMC 12/3/12
878 ! sublimation of cloud ice
879 real(r8) :: ice_sublim(mgncol,nlev) ! sublimation from ice to vapor PMC 12/3/12
880 ! ice nucleation
881 real(r8) :: nnuccd(mgncol,nlev) ! number rate from deposition/cond.-freezing
882 real(r8) :: mnuccd(mgncol,nlev) ! mass mixing ratio
883 ! freezing of cloud water
884 real(r8) :: mnuccc(mgncol,nlev) ! mass mixing ratio
885 real(r8) :: nnuccc(mgncol,nlev) ! number concentration
886 ! contact freezing of cloud water
887 real(r8) :: mnucct(mgncol,nlev) ! mass mixing ratio
888 real(r8) :: nnucct(mgncol,nlev) ! number concentration
889 ! deposition nucleation in mixed-phase clouds (from external scheme)
890 real(r8) :: mnudep(mgncol,nlev) ! mass mixing ratio
891 real(r8) :: nnudep(mgncol,nlev) ! number concentration
892 ! ice multiplication
893 real(r8) :: msacwi(mgncol,nlev) ! mass mixing ratio
894 real(r8) :: nsacwi(mgncol,nlev) ! number concentration
895 ! autoconversion of cloud droplets
896 real(r8) :: prc(mgncol,nlev) ! mass mixing ratio
897 real(r8) :: nprc(mgncol,nlev) ! number concentration (rain)
898 real(r8) :: nprc1(mgncol,nlev) ! number concentration (cloud droplets)
899 ! self-aggregation of snow
900 real(r8) :: nsagg(mgncol,nlev) ! number concentration
901 ! self-collection of rain
902 real(r8) :: nragg(mgncol,nlev) ! number concentration
903 ! collection of droplets by snow
904 real(r8) :: psacws(mgncol,nlev) ! mass mixing ratio
905 real(r8) :: npsacws(mgncol,nlev) ! number concentration
906 ! collection of rain by snow
907 real(r8) :: pracs(mgncol,nlev) ! mass mixing ratio
908 real(r8) :: npracs(mgncol,nlev) ! number concentration
909 ! freezing of rain
910 real(r8) :: mnuccr(mgncol,nlev) ! mass mixing ratio
911 real(r8) :: nnuccr(mgncol,nlev) ! number concentration
912 ! freezing of rain to form ice (mg add 4/26/13)
913 real(r8) :: mnuccri(mgncol,nlev) ! mass mixing ratio
914 real(r8) :: nnuccri(mgncol,nlev) ! number concentration
915 ! accretion of droplets by rain
916 real(r8) :: pra(mgncol,nlev) ! mass mixing ratio
917 real(r8) :: npra(mgncol,nlev) ! number concentration
918 ! autoconversion of cloud ice to snow
919 real(r8) :: prci(mgncol,nlev) ! mass mixing ratio
920 real(r8) :: nprci(mgncol,nlev) ! number concentration
921 ! accretion of cloud ice by snow
922 real(r8) :: prai(mgncol,nlev) ! mass mixing ratio
923 real(r8) :: nprai(mgncol,nlev) ! number concentration
924 ! evaporation of rain
925 real(r8) :: pre(mgncol,nlev) ! mass mixing ratio
926 ! sublimation of snow
927 real(r8) :: prds(mgncol,nlev) ! mass mixing ratio
928 ! number evaporation
929 real(r8) :: nsubi(mgncol,nlev) ! cloud ice
930 real(r8) :: nsubc(mgncol,nlev) ! droplet
931 real(r8) :: nsubs(mgncol,nlev) ! snow
932 real(r8) :: nsubr(mgncol,nlev) ! rain
933 ! bergeron process
934 real(r8) :: berg(mgncol,nlev) ! mass mixing ratio (cloud ice)
935 real(r8) :: bergs(mgncol,nlev) ! mass mixing ratio (snow)
936
937!++ag
938 !graupel/hail processes
939 real(r8) :: npracg(mgncol,nlev) ! change n collection rain by graupel (precipf)
940 real(r8) :: nscng(mgncol,nlev) ! change n conversion to graupel due to collection droplets by snow (lcldm)
941 real(r8) :: ngracs(mgncol,nlev) ! change n conversion to graupel due to collection rain by snow (precipf)
942 real(r8) :: nmultg(mgncol,nlev) ! ice mult due to acc droplets by graupel (lcldm)
943 real(r8) :: nmultrg(mgncol,nlev) ! ice mult due to acc rain by graupel (precipf)
944 real(r8) :: npsacwg(mgncol,nlev) ! change n collection droplets by graupel (lcldm)
945
946 real(r8) :: psacr(mgncol,nlev) ! conversion due to coll of snow by rain (precipf)
947 real(r8) :: pracg(mgncol,nlev) ! change in q collection rain by graupel (precipf)
948 real(r8) :: psacwg(mgncol,nlev) ! change in q collection droplets by graupel (lcldm)
949 real(r8) :: pgsacw(mgncol,nlev) ! conversion q to graupel due to collection droplets by snow (lcldm)
950 real(r8) :: pgracs(mgncol,nlev) ! conversion q to graupel due to collection rain by snow (precipf)
951 real(r8) :: prdg(mgncol,nlev) ! dep of graupel (precipf)
952! real(r8) :: eprdg(mgncol,nlev) ! evap/sub of graupel (precipf)
953 real(r8) :: qmultg(mgncol,nlev) ! change q due to ice mult droplets/graupel (lcldm)
954 real(r8) :: qmultrg(mgncol,nlev) ! change q due to ice mult rain/graupel (precipf)
955!--ag
956
957
958 ! fallspeeds
959 ! number-weighted
960 real(r8) :: uns(mgncol,nlev) ! snow
961 real(r8) :: unr(mgncol,nlev) ! rain
962!++ag
963 real(r8) :: ung(mgncol,nlev) ! graupel/hail
964!--ag
965 ! air density corrected fallspeed parameters
966 real(r8) :: arn(mgncol,nlev) ! rain
967 real(r8) :: asn(mgncol,nlev) ! snow
968!++a
969 real(r8) :: agn(mgncol,nlev) ! graupel
970!--ag
971 real(r8) :: acn(mgncol,nlev) ! cloud droplet
972 real(r8) :: ain(mgncol,nlev) ! cloud ice
973 real(r8) :: ajn(mgncol,nlev) ! cloud small ice
974
975 ! Mass of liquid droplets used with external heterogeneous freezing.
976 real(r8) :: mi0l(mgncol)
977
978 ! saturation vapor pressures
979 real(r8) :: esl(mgncol,nlev) ! liquid
980 real(r8) :: esi(mgncol,nlev) ! ice
981 real(r8) :: esn ! checking for RH after rain evap
982
983 ! saturation vapor mixing ratios
984 real(r8) :: qvl(mgncol,nlev) ! liquid
985 real(r8) :: qvi(mgncol,nlev) ! ice
986 real(r8) :: qvn ! checking for RH after rain evap
987
988 ! relative humidity
989 real(r8) :: relhum(mgncol,nlev)
990
991 ! parameters for cloud water and cloud ice sedimentation calculations
992 real(r8) :: fc(mgncol,nlev)
993 real(r8) :: fnc(mgncol,nlev)
994 real(r8) :: fi(mgncol,nlev)
995 real(r8) :: fni(mgncol,nlev)
996
997!++ag
998 real(r8) :: fg(mgncol,nlev)
999 real(r8) :: fng(mgncol,nlev)
1000!--ag
1001
1002 real(r8) :: fr(mgncol,nlev)
1003 real(r8) :: fnr(mgncol,nlev)
1004 real(r8) :: fs(mgncol,nlev)
1005 real(r8) :: fns(mgncol,nlev)
1006
1007 real(r8) :: faloutc(nlev)
1008 real(r8) :: faloutnc(nlev)
1009 real(r8) :: falouti(nlev)
1010 real(r8) :: faloutni(nlev)
1011
1012 real(r8) :: faloutr(nlev)
1013 real(r8) :: faloutnr(nlev)
1014 real(r8) :: falouts(nlev)
1015 real(r8) :: faloutns(nlev)
1016
1017 real(r8) :: faltndc
1018 real(r8) :: faltndnc
1019 real(r8) :: faltndi
1020 real(r8) :: faltndni
1021 real(r8) :: faltndqie
1022 real(r8) :: faltndqce
1023
1024 real(r8) :: faltndr
1025 real(r8) :: faltndnr
1026 real(r8) :: faltnds
1027 real(r8) :: faltndns
1028
1029!++ag
1030 real(r8) :: faloutg(nlev)
1031 real(r8) :: faloutng(nlev)
1032 real(r8) :: faltndg
1033 real(r8) :: faltndng
1034!--ag
1035
1036 real(r8) :: rainrt(mgncol,nlev) ! rain rate for reflectivity calculation
1037
1038 ! dummy variables
1039 real(r8) :: dum
1040 real(r8) :: dum1
1041 real(r8) :: dum2
1042!++ag
1043 real(r8) :: dum3
1044!--ag
1045 real(r8) :: dumni0
1046 real(r8) :: dumns0
1047 real(r8) :: tx1, tx2, tx3, tx4, tx5, tx6, tx7, grho
1048 ! dummies for checking RH
1049 real(r8) :: qtmp
1050 real(r8) :: ttmp
1051 ! dummies for conservation check
1052 real(r8) :: ratio
1053 real(r8) :: tmpfrz
1054 ! dummies for in-cloud variables
1055 real(r8) :: dumc(mgncol,nlev) ! qc
1056 real(r8) :: dumnc(mgncol,nlev) ! nc
1057 real(r8) :: dumi(mgncol,nlev) ! qi
1058 real(r8) :: dumni(mgncol,nlev) ! ni
1059 real(r8) :: dumr(mgncol,nlev) ! rain mixing ratio
1060 real(r8) :: dumnr(mgncol,nlev) ! rain number concentration
1061 real(r8) :: dums(mgncol,nlev) ! snow mixing ratio
1062 real(r8) :: dumns(mgncol,nlev) ! snow number concentration
1063!++ag
1064 real(r8) :: dumg(mgncol,nlev) ! graupel mixing ratio
1065 real(r8) :: dumng(mgncol,nlev) ! graupel number concentration
1066!--ag
1067 ! Array dummy variable
1068! real(r8) :: dum_2D(mgncol,nlev)
1069 real(r8) :: pdel_inv(mgncol,nlev)
1070 real(r8) :: ts_au_loc(mgncol)
1071
1072 ! loop array variables
1073 ! "i" and "k" are column/level iterators for internal (MG) variables
1074 ! "n" is used for other looping (currently just sedimentation)
1075 integer i, k, n
1076
1077 ! number of sub-steps for loops over "n" (for sedimentation)
1078 integer nstep, mdust, nlb, nstep_def
1079
1080 ! Varaibles to scale fall velocity between small and regular ice regimes.
1081! real(r8) :: irad, ifrac, tsfac
1082 real(r8) :: irad, ifrac
1083! logical, parameter :: do_ice_gmao=.false., do_liq_liu=.false.
1084! logical, parameter :: do_ice_gmao=.false., do_liq_liu=.true.
1085! logical, parameter :: do_ice_gmao=.true., do_liq_liu=.false.
1086! real(r8), parameter :: qimax=0.010, qimin=0.001, qiinv=one/(qimax-qimin), &
1087! real(r8), parameter :: qimax=0.010, qimin=0.001, qiinv=one/(qimax-qimin), &
1088 real(r8), parameter :: qimax=0.010_r8, qimin=0.005_r8, qiinv=one/(qimax-qimin)
1089! ts_au_min=180.0
1090
1091 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1092
1093
1094 ! Process inputs
1095
1097 deltat = deltatin
1098 oneodt = one / deltat
1099! nstep_def = max(1, nint(deltat/20))
1100 nstep_def = max(1, nint(deltat/5))
1101! tsfac = log(ts_au/ts_au_min) * qiinv
1102
1104 do k=1,nlev
1105 do i=1,mgncol
1106 qc(i,k) = qcn(i,k)
1107 nc(i,k) = ncn(i,k)
1108 qi(i,k) = qin(i,k)
1109 ni(i,k) = nin(i,k)
1110 qr(i,k) = qrn(i,k)
1111 nr(i,k) = nrn(i,k)
1112 qs(i,k) = qsn(i,k)
1113 ns(i,k) = nsn(i,k)
1114!++ag
1115 qg(i,k) = qgr(i,k)
1116 ng(i,k) = ngr(i,k)
1117 enddo
1118 enddo
1119
1120 ! cldn: used to set cldm, unused for subcolumns
1121 ! liqcldf: used to set lcldm, unused for subcolumns
1122 ! icecldf: used to set icldm, unused for subcolumns
1124 if (microp_uniform) then
1125 ! subcolumns, set cloud fraction variables to one
1126 ! if cloud water or ice is present, if not present
1127 ! set to mincld (mincld used instead of zero, to prevent
1128 ! possible division by zero errors).
1129
1130 do k=1,nlev
1131 do i=1,mgncol
1132
1133 if (qc(i,k) >= qsmall) then
1134 lcldm(i,k) = one
1135 else
1136 lcldm(i,k) = mincld
1137 endif
1138
1139 if (qi(i,k) >= qsmall) then
1140 icldm(i,k) = one
1141 else
1142 icldm(i,k) = mincld
1143 endif
1144
1145 cldm(i,k) = max(icldm(i,k), lcldm(i,k))
1146! qsfm(i,k) = one
1147 qsfm(i,k) = qsatfac(i,k)
1148 enddo
1149 enddo
1150
1151 else ! get cloud fraction, check for minimum
1152 do k=1,nlev
1153 do i=1,mgncol
1154 cldm(i,k) = max(cldn(i,k), mincld)
1155 lcldm(i,k) = max(liqcldf(i,k), mincld)
1156 icldm(i,k) = max(icecldf(i,k), mincld)
1157 qsfm(i,k) = qsatfac(i,k)
1158 enddo
1159 enddo
1160 endif
1161
1162! if (lprnt) write(0,*)' cldm=',cldm(1,nlev-20:nlev)
1163! if (lprnt) write(0,*)' liqcldf=',liqcldf(1,nlev-20:nlev)
1164! if (lprnt) write(0,*)' lcldm=',lcldm(1,nlev-20:nlev)
1165! if (lprnt) write(0,*)' icecldf=',icecldf(1,nlev-20:nlev)
1166! if (lprnt) write(0,*)' icldm=',icldm(1,nlev-20:nlev)
1167! if (lprnt) write(0,*)' qsfm=',qsfm(1,nlev-20:nlev)
1168
1170
1171 ! local physical properties
1172
1173! write(0,*)' in mg2 T=',t(1,:)
1174! write(0,*)' in mg2 P=',p(1,:),' r=',r
1175 do k=1,nlev
1176 do i=1,mgncol
1177! rho(i,k) = p(i,k) / (r*t(i,k)*(one+fv*q(i,k)))
1178 rho(i,k) = p(i,k) / (r*t(i,k))
1179 rhoinv(i,k) = one / rho(i,k)
1180 dv(i,k) = 8.794e-5_r8 * t(i,k)**1.81_r8 / p(i,k)
1181 mu(i,k) = 1.496e-6_r8 * t(i,k)*sqrt(t(i,k)) / (t(i,k) + 120._r8)
1182 sc(i,k) = mu(i,k) / (rho(i,k)*dv(i,k))
1183
1184 ! air density adjustment for fallspeed parameters
1185 ! includes air density correction factor to the
1186 ! power of 0.54 following Heymsfield and Bansemer 2007
1187
1188 rhof(i,k) = (rhosu*rhoinv(i,k))**0.54_r8
1189
1190 arn(i,k) = ar * rhof(i,k)
1191 asn(i,k) = as * rhof(i,k)
1192!++ag if do hail then agn = ah *rhof else ag*rhof
1193 agn(i,k) = agtmp * rhof(i,k)
1194 acn(i,k) = g*rhow/(18._r8*mu(i,k))
1195 tx1 = (rhosu*rhoinv(i,k))**0.35_r8
1196 ain(i,k) = ai * tx1
1197 ajn(i,k) = aj * tx1
1198
1199 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1200 ! Get humidity and saturation vapor pressures
1201
1202! do k=1,nlev
1203! do i=1,mgncol
1204! relvar(i,k) = relvar_i
1205 accre_enhan(i,k) = accre_enhan_i
1206! call qsat_water(t(i,k), p(i,k), esl(i,k), qvl(i,k))
1207 esl(i,k) = min(fpvsl(t(i,k)), p(i,k))
1208 qvl(i,k) = epsqs*esl(i,k) / (p(i,k)-omeps*esl(i,k))
1209
1210
1211 ! make sure when above freezing that esi=esl, not active yet
1212 if (t(i,k) >= tmelt) then
1213 esi(i,k) = esl(i,k)
1214 qvi(i,k) = qvl(i,k)
1215 else
1216! call qsat_ice(t(i,k), p(i,k), esi(i,k), qvi(i,k))
1217 esi(i,k) = min(fpvsi(t(i,k)), p(i,k))
1218 qvi(i,k) = epsqs*esi(i,k) / (p(i,k)-omeps*esi(i,k))
1219 end if
1220
1221 ! Scale the water saturation values to reflect subgrid scale
1222 ! ice cloud fraction, where ice clouds begin forming at a
1223 ! gridbox average relative humidity of rhmini (not 1).
1224 !
1225 ! NOTE: For subcolumns and other non-subgrid clouds, qsfm will be 1.
1226 qvi(i,k) = qsfm(i,k) * qvi(i,k)
1227! esi(i,k) = qsfm(i,k) * esi(i,k)
1228 qvl(i,k) = qsfm(i,k) * qvl(i,k)
1229! esl(i,k) = qsfm(i,k) * esl(i,k)
1230
1231 relhum(i,k) = max(zero, min(q(i,k)/max(qvl(i,k), qsmall), two))
1232 enddo
1233 enddo
1234
1235 !===============================================
1236
1237 ! set mtime here to avoid answer-changing
1238 mtime = deltat
1239
1241 do k=1,nlev
1242 do i=1,mgncol
1243 qcsevap(i,k) = zero
1244 qisevap(i,k) = zero
1245 qvres(i,k) = zero
1246 cmeitot(i,k) = zero
1247 vtrmc(i,k) = zero
1248 vtrmi(i,k) = zero
1249 qcsedten(i,k) = zero
1250 qisedten(i,k) = zero
1251 qrsedten(i,k) = zero
1252 qssedten(i,k) = zero
1253!++ag
1254 qgsedten(i,k) = zero
1255!--ag
1256
1257
1258 pratot(i,k) = zero
1259 prctot(i,k) = zero
1260 mnuccctot(i,k) = zero
1261 mnuccttot(i,k) = zero
1262 msacwitot(i,k) = zero
1263 psacwstot(i,k) = zero
1264 bergstot(i,k) = zero
1265 bergtot(i,k) = zero
1266 melttot(i,k) = zero
1267 homotot(i,k) = zero
1268 qcrestot(i,k) = zero
1269 prcitot(i,k) = zero
1270 praitot(i,k) = zero
1271 qirestot(i,k) = zero
1272 mnuccrtot(i,k) = zero
1273!++ag
1274 mnuccritot(i,k) = zero
1275!--ag
1276
1277 pracstot(i,k) = zero
1278 meltsdttot(i,k) = zero
1279 frzrdttot(i,k) = zero
1280 mnuccdtot(i,k) = zero
1281
1282!++ag
1283 psacrtot(i,k) = zero
1284 pracgtot(i,k) = zero
1285 psacwgtot(i,k) = zero
1286 pgsacwtot(i,k) = zero
1287 pgracstot(i,k) = zero
1288 prdgtot(i,k) = zero
1289! eprdgtot(i,k) = zero
1290 qmultgtot(i,k) = zero
1291 qmultrgtot(i,k) = zero
1292 npracgtot(i,k) = zero
1293 nscngtot(i,k) = zero
1294 ngracstot(i,k) = zero
1295 nmultgtot(i,k) = zero
1296 nmultrgtot(i,k) = zero
1297 npsacwgtot(i,k) = zero
1298!need to zero these out to be totally switchable (for conservation)
1299 psacr(i,k) = zero
1300 pracg(i,k) = zero
1301 psacwg(i,k) = zero
1302 pgsacw(i,k) = zero
1303 pgracs(i,k) = zero
1304
1305 prdg(i,k) = zero
1306! eprdg(i,k) = zero
1307 qmultg(i,k) = zero
1308 qmultrg(i,k) = zero
1309 npracg(i,k) = zero
1310 nscng(i,k) = zero
1311 ngracs(i,k) = zero
1312 nmultg(i,k) = zero
1313 nmultrg(i,k) = zero
1314 npsacwg(i,k) = zero
1315!--ag
1316 rflx(i,k+1) = zero
1317 sflx(i,k+1) = zero
1318 lflx(i,k+1) = zero
1319 iflx(i,k+1) = zero
1320!++ag
1321 gflx(i,k+1) = zero
1322!--ag
1323
1325
1326 qrout(i,k) = zero
1327 qsout(i,k) = zero
1328 nrout(i,k) = zero
1329 nsout(i,k) = zero
1330!++ag
1331 qgout(i,k) = zero
1332 ngout(i,k) = zero
1333 dgout(i,k) = zero
1334!--ag
1335
1336 ! for refl calc
1337 rainrt(i,k) = zero
1338
1340 rercld(i,k) = zero
1341
1342 qcsinksum_rate1ord(i,k) = zero
1343
1345 nevapr(i,k) = zero
1346 prer_evap(i,k) = zero
1347 evapsnow(i,k) = zero
1348 am_evp_st(i,k) = zero
1349 prain(i,k) = zero
1350 prodsnow(i,k) = zero
1351 cmeout(i,k) = zero
1352
1353 precip_frac(i,k) = mincld
1354
1355 lamc(i,k) = zero
1356
1358
1359 tlat(i,k) = zero
1360 qvlat(i,k) = zero
1361 qctend(i,k) = zero
1362 qitend(i,k) = zero
1363 qstend(i,k) = zero
1364 qrtend(i,k) = zero
1365 nctend(i,k) = zero
1366 nitend(i,k) = zero
1367 nrtend(i,k) = zero
1368 nstend(i,k) = zero
1369!++ag
1370 qgtend(i,k) = zero
1371 ngtend(i,k) = zero
1372!--ag
1373
1375 qcic(i,k) = zero
1376 qiic(i,k) = zero
1377 qsic(i,k) = zero
1378 qric(i,k) = zero
1379!++ag
1380 qgic(i,k) = zero
1381!--ag
1382
1383
1384 ncic(i,k) = zero
1385 niic(i,k) = zero
1386 nsic(i,k) = zero
1387 nric(i,k) = zero
1388!++ag
1389 ngic(i,k) = zero
1390!--ag
1392 ums(i,k) = zero
1393 uns(i,k) = zero
1394 umr(i,k) = zero
1395 unr(i,k) = zero
1396!++ag
1397 umg(i,k) = zero
1398 ung(i,k) = zero
1399!--ag
1400
1402 qcrat(i,k) = one
1403
1404 ! Many outputs have to be initialized here at the top to work around
1405 ! ifort problems, even if they are always overwritten later.
1406 effc(i,k) = ten
1407 lamcrad(i,k) = zero
1408 pgamrad(i,k) = zero
1409 effc_fn(i,k) = ten
1410 effi(i,k) = 25._r8
1411 sadice(i,k) = zero
1412 sadsnow(i,k) = zero
1413 deffi(i,k) = 50._r8
1414
1415 qrout2(i,k) = zero
1416 nrout2(i,k) = zero
1417 drout2(i,k) = zero
1418 qsout2(i,k) = zero
1419 nsout2(i,k) = zero
1420 dsout(i,k) = zero
1421 dsout2(i,k) = zero
1422!++ag
1423 qgout2(i,k) = zero
1424 ngout2(i,k) = zero
1425 freqg(i,k) = zero
1426 dgout2(i,k) = zero
1427!--ag
1428
1429 freqr(i,k) = zero
1430 freqs(i,k) = zero
1431
1432 reff_rain(i,k) = zero
1433 reff_snow(i,k) = zero
1434!++ag
1435 reff_grau(i,k) = zero
1436 lamg(i,k) = zero
1437 n0g(i,k) = zero
1438!--ag
1439
1440 refl(i,k) = -9999._r8
1441 arefl(i,k) = zero
1442 areflz(i,k) = zero
1443 frefl(i,k) = zero
1444 csrfl(i,k) = zero
1445 acsrfl(i,k) = zero
1446 fcsrfl(i,k) = zero
1447
1448 ncal(i,k) = zero
1449 ncai(i,k) = zero
1450
1451 nfice(i,k) = zero
1452 npccn(i,k) = zero
1453 enddo
1454 enddo
1456 if (iccn == 1) then
1457 do k=1,nlev
1458 do i=1,mgncol
1459 npccn(i,k) = npccnin(i,k)
1460 enddo
1461 enddo
1462 else
1463 do k=1,nlev
1464 do i=1,mgncol
1465 npccn(i,k) = max((npccnin(i,k)*lcldm(i,k)-nc(i,k))*oneodt, zero)
1466 enddo
1467 enddo
1468 endif
1469
1471
1472 do i=1,mgncol
1473 prect(i) = zero
1474 preci(i) = zero
1475 enddo
1476
1477 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1478 ! droplet activation
1479 ! get provisional droplet number after activation. This is used for
1480 ! all microphysical process calculations, for consistency with update of
1481 ! droplet mass before microphysics
1482
1483 ! calculate potential for droplet activation if cloud water is present
1484 ! tendency from activation (npccn) is read in from companion routine
1485
1486 ! output activated liquid and ice (convert from #/kg -> #/m3)
1487 !--------------------------------------------------
1488! where (qc >= qsmall .and. lcldm > mincld)
1489! where (qc >= qsmall)
1490! npccn = max((npccnin*lcldm-nc)*oneodt, zero)
1491! nc = max(nc + npccn*deltat, zero)
1492! ncal = nc*rho/lcldm ! sghan minimum in #/cm3
1493! elsewhere
1494! ncal = zero
1495! end where
1496
1497! if (lprnt) write(0,*)' nc1=',nc(1,:)
1498 do k=1,nlev
1499 do i=1,mgncol
1500 if (qc(i,k) > qsmall .and. lcldm(i,k) >= mincld) then
1501 npccn(i,k) = max((npccnin(i,k)*lcldm(i,k)-nc(i,k))*oneodt, zero)
1502 nc(i,k) = max(nc(i,k) + npccn(i,k)*deltat, zero)
1503 ncal(i,k) = nc(i,k) * rho(i,k) / lcldm(i,k)
1504 else
1505 ncal(i,k) = 0.0
1506 endif
1507 enddo
1508 enddo
1509
1510 if (iccn == 1) then
1511 do k=1,nlev
1512 do i=1,mgncol
1513 if (t(i,k) < icenuct) then
1514 ncai(i,k) = 0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k))) * 1000._r8
1515! ncai(i,k) = min(ncai(i,k), 208.9e3_r8)
1516 ncai(i,k) = min(ncai(i,k), 355.0e3_r8)
1517 naai(i,k) = (ncai(i,k)*rhoinv(i,k) + naai(i,k)) * half
1518 ncai(i,k) = naai(i,k)*rho(i,k)
1519 else
1520 naai(i,k) = zero
1521 ncai(i,k) = zero
1522 endif
1523 enddo
1524 enddo
1525 elseif (iccn == 2) then
1526 do k=1,nlev
1527 do i=1,mgncol
1528 if (t(i,k) < icenuct) then
1529 ncai(i,k) = naai(i,k)*rho(i,k)
1530 ncai(i,k) = min(ncai(i,k), 710.0e3_r8)
1531 naai(i,k) = ncai(i,k)*rhoinv(i,k)
1532 else
1533 naai(i,k) = zero
1534 ncai(i,k) = zero
1535 endif
1536 enddo
1537 enddo
1538 else
1539 do k=1,nlev
1540 do i=1,mgncol
1541 if (t(i,k) < icenuct) then
1542 ncai(i,k) = 0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k))) * 1000._r8
1543 ncai(i,k) = min(ncai(i,k), 355.0e3_r8)
1544 naai(i,k) = ncai(i,k)*rhoinv(i,k)
1545 else
1546 naai(i,k) = zero
1547 ncai(i,k) = zero
1548 endif
1549 enddo
1550 enddo
1551
1552 endif
1553
1554
1555 !===============================================
1556
1557 ! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5%
1558 !
1559 ! NOTE: If using gridbox average values, condensation will not occur until rh=1,
1560 ! so the threshold seems like it should be 1.05 and not rhmini + 0.05. For subgrid
1561 ! clouds (using rhmini and qsfacm), the relhum has already been adjusted, and thus
1562 ! the nucleation threshold should also be 1.05 and not rhmini + 0.05.
1563
1564 !-------------------------------------------------------
1565
1566 if (do_cldice) then
1567 where (naai > zero .and. t < icenuct .and. relhum*esl/esi > 1.05_r8)
1568! where (naai > zero .and. t < icenuct .and. relhum*esl/esi > 1.05_r8 &
1569! .and. icldm > mincld )
1570
1571 !if NAAI > 0. then set numice = naai (as before)
1572 !note: this is gridbox averaged
1573 nnuccd = (naai-ni/icldm)/mtime*icldm
1574 nnuccd = max(nnuccd, zero)
1575 nimax = naai*icldm
1576
1577 !Calc mass of new particles using new crystal mass...
1578 !also this will be multiplied by mtime as nnuccd is...
1579
1580 mnuccd = nnuccd * mi0
1581
1582 elsewhere
1583 nnuccd = zero
1584 nimax = zero
1585 mnuccd = zero
1586 end where
1587
1588 endif
1589
1590
1591 !=============================================================================
1592 do k=1,nlev
1593
1594 do i=1,mgncol
1595
1596 ! calculate instantaneous precip processes (melting and homogeneous freezing)
1597
1598 ! melting of snow at +2 C
1599
1600 if (t(i,k) > snowmelt) then
1601 if (qs(i,k) > zero) then
1602
1603 ! make sure melting snow doesn't reduce temperature below threshold
1604 dum = -(xlf/cpp) * qs(i,k)
1605 if (t(i,k)+dum < snowmelt) then
1606 dum = min(one, max(zero, (cpp/xlf)*(t(i,k)-snowmelt)/qs(i,k)))
1607 else
1608 dum = one
1609 end if
1610
1611 minstsm(i,k) = dum*qs(i,k)
1612 ninstsm(i,k) = dum*ns(i,k)
1613
1614 dum1 = - minstsm(i,k) * (xlf*oneodt)
1615 tlat(i,k) = tlat(i,k) + dum1
1616 meltsdttot(i,k) = meltsdttot(i,k) + dum1
1617
1618! if (lprnt .and. k >=40) write(0,*)' tlats=',tlat(i,k),' dum1=',dum1,&
1619! ' minstsm=',minstsm(i,k),' qs=',qs(i,k),' xlf=',xlf,' oneodt=',oneodt, &
1620! ' snowmelt=',snowmelt,' t=',t(i,k),' dum=',dum,' k=',k
1621
1622 qs(i,k) = max(qs(i,k) - minstsm(i,k), zero)
1623 ns(i,k) = max(ns(i,k) - ninstsm(i,k), zero)
1624 qr(i,k) = max(qr(i,k) + minstsm(i,k), zero)
1625 nr(i,k) = max(nr(i,k) + ninstsm(i,k), zero)
1626 endif
1627 endif
1628
1629 enddo
1630 enddo
1631! if (lprnt) write(0,*)' tlat1=',tlat(1,:)*deltat
1632! if (lprnt) write(0,*)' qg1=',qg(1,:)
1633
1634!++ag
1635
1636 if (do_graupel .or. do_hail) then
1637! melting of graupel at +2 C
1638
1639 do k=1,nlev
1640 do i=1,mgncol
1641
1642 if (t(i,k) > snowmelt) then
1643 if (qg(i,k) > zero) then
1644
1645! make sure melting graupel doesn't reduce temperature below threshold
1646 dum = -(xlf/cpp) * qg(i,k)
1647 if (t(i,k)+dum < snowmelt) then
1648 dum = max(zero, min(one, (cpp/xlf)*(t(i,k)-snowmelt)/qg(i,k)))
1649 else
1650 dum = one
1651 end if
1652
1653 minstgm(i,k) = dum*qg(i,k)
1654 ninstgm(i,k) = dum*ng(i,k)
1655
1656 dum1 = - minstgm(i,k) * (xlf*oneodt)
1657 tlat(i,k) = dum1 + tlat(i,k)
1658 meltsdttot(i,k) = dum1 + meltsdttot(i,k)
1659
1660! if (lprnt .and. k >=40) write(0,*)' tlatg=',tlat(i,k),' dum1=',dum1,&
1661! ' minstgm=',minstgm(i,k),' qg=',qg(i,k),' xlf=',xlf,' oneodt=',oneodt, &
1662! ' snowmelt=',snowmelt,' t=',t(i,k),' k=',k,' cpp=',cpp
1663
1664 qg(i,k) = max(qg(i,k) - minstgm(i,k), zero)
1665 ng(i,k) = max(ng(i,k) - ninstgm(i,k), zero)
1666 qr(i,k) = max(qr(i,k) + minstgm(i,k), zero)
1667 nr(i,k) = max(nr(i,k) + ninstgm(i,k), zero)
1668 endif
1669 endif
1670
1671 enddo
1672 enddo
1673 endif
1674
1675! if (lprnt) write(0,*)' tlat1g=',tlat(1,:)*deltat
1676!--ag
1677
1678 do k=1,nlev
1679 do i=1,mgncol
1680 ! freezing of rain at -5 C
1681
1682 if (t(i,k) < rainfrze) then
1683
1684 if (qr(i,k) > zero) then
1685
1686 ! make sure freezing rain doesn't increase temperature above threshold
1687 dum = (xlf/cpp) * qr(i,k)
1688 if (t(i,k)+dum > rainfrze) then
1689 dum = -(t(i,k)-rainfrze) * (cpp/xlf)
1690 dum = min(one, max(zero, dum/qr(i,k)))
1691 else
1692 dum = one
1693 end if
1694
1695 minstrf(i,k) = dum*qr(i,k)
1696 ninstrf(i,k) = dum*nr(i,k)
1697
1698 ! heating tendency
1699 dum1 = minstrf(i,k) * (xlf*oneodt)
1700 tlat(i,k) = tlat(i,k) + dum1
1701 frzrdttot(i,k) = frzrdttot(i,k) + dum1
1702
1703 qr(i,k) = max(qr(i,k) - minstrf(i,k), zero)
1704 nr(i,k) = max(nr(i,k) - ninstrf(i,k), zero)
1705
1706!++ag
1707! freeze rain to graupel not snow.
1708 if(do_hail .or. do_graupel) then
1709 qg(i,k) = max(qg(i,k) + minstrf(i,k), zero)
1710 ng(i,k) = max(ng(i,k) + ninstrf(i,k), zero)
1711 else
1712 qs(i,k) = max(qs(i,k) + minstrf(i,k), zero)
1713 ns(i,k) = max(ns(i,k) + ninstrf(i,k), zero)
1714 end if
1715!--ag
1716
1717 endif
1718 endif
1719 enddo
1720 enddo
1721
1722! if (lprnt) then
1723! write(0,*)' tlat2=',tlat(1,:)*deltat
1724! write(0,*)' lcldm=',lcldm(1,:)
1725! write(0,*)' qc=',qc(1,:)
1726! write(0,*)' nc=',nc(1,:)
1727! write(0,*)' qg2=',qg(1,:)
1728! endif
1729
1730 do k=1,nlev
1731 do i=1,mgncol
1732 ! obtain in-cloud values of cloud water/ice mixing ratios and number concentrations
1733 !-------------------------------------------------------
1734 ! for microphysical process calculations
1735 ! units are kg/kg for mixing ratio, 1/kg for number conc
1736
1737! if (qc(i,k) >= qsmall .and. lcldm(i,k) > mincld) then
1738 if (qc(i,k) >= qsmall) then
1739 ! limit in-cloud values to 0.005 kg/kg
1740 dum = one / lcldm(i,k)
1741! qcic(i,k) = min(qc(i,k)*dum, 5.e-3_r8) ! limit in-cloud values to 0.005 kg/kg
1742 qcic(i,k) = min(qc(i,k)*dum, 0.05_r8) ! limit in-cloud values to 0.05 kg/kg
1743 ncic(i,k) = max(nc(i,k)*dum, zero)
1744
1745 ! specify droplet concentration
1746 if (nccons) then
1747 ncic(i,k) = ncnst * rhoinv(i,k)
1748 endif
1749 else
1750 qcic(i,k) = zero
1751 ncic(i,k) = zero
1752 endif
1753
1754! if (qi(i,k) >= qsmall .and. icldm(i,k) > mincld) then
1755 if (qi(i,k) >= qsmall) then
1756 dum = one / icldm(i,k)
1757! qiic(i,k) = min(qi(i,k)*dum, 5.e-3_r8) ! limit in-cloud values to 0.005 kg/kg
1758 qiic(i,k) = min(qi(i,k)*dum, 0.05_r8) ! limit in-cloud values to 0.05 kg/kg
1759 niic(i,k) = max(ni(i,k)*dum, zero)
1760
1761 ! switch for specification of cloud ice number
1762 if (nicons) then
1763 niic(i,k) = ninst * rhoinv(i,k)
1764 endif
1765 else
1766 qiic(i,k) = zero
1767 niic(i,k) = zero
1768 endif
1769
1770 enddo
1771 enddo
1772
1773 !========================================================================
1774
1775 ! for sub-columns cldm has already been set to 1 if cloud
1776 ! water or ice is present, so precip_frac will be correctly set below
1777 ! and nothing extra needs to be done here
1778
1779 precip_frac = cldm
1780
1781 micro_vert_loop: do k=1,nlev
1782
1783 if (trim(micro_mg_precip_frac_method) == 'in_cloud') then
1784
1785 if (k /= 1) then
1786 where (qc(:,k) < qsmall .and. qi(:,k) < qsmall)
1787 precip_frac(:,k) = precip_frac(:,k-1)
1788 end where
1789 endif
1790
1791 else if (trim(micro_mg_precip_frac_method) == 'max_overlap') then
1792
1793!++ag add graupel to precip frac?
1794 ! calculate precip fraction based on maximum overlap assumption
1795
1796 ! if rain or snow mix ratios are smaller than threshold,
1797 ! then leave precip_frac as cloud fraction at current level
1798 if (k /= 1) then
1799!++ag
1800! where (qr(:,k-1) >= qsmall .or. qs(:,k-1) >= qsmall .or. qg(:,k-1) >= qsmall)
1801!--ag
1802 where (qr(:,k-1) >= qsmall .or. qs(:,k-1) >= qsmall)
1803 precip_frac(:,k) = max(precip_frac(:,k-1), precip_frac(:,k))
1804 end where
1805 endif
1806
1807 endif
1808
1809
1810 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1811 ! get size distribution parameters based on in-cloud cloud water
1812 ! these calculations also ensure consistency between number and mixing ratio
1813 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1814
1815 ! cloud liquid
1816 !-------------------------------------------
1817
1818! if (lprnt .and. k>=60 .and. k<=65) then
1819! if (lprnt .and. k>=100) then
1820! if (lprnt) then
1821! write(0,*)' pgam=',pgam(1,k), ' qcic=',qcic(1,k),' ncic=',ncic(1,k),' rho=',rho(1,k),' k=',k
1822! endif
1823 call size_dist_param_liq(mg_liq_props, qcic(:,k), ncic(:,k), rho(:,k), &
1824 pgam(:,k), lamc(:,k), mgncol)
1825! if (lprnt .and. k>=60 .and. k<=65) then
1826! if (lprnt .and. k>=100) then
1827! if (lprnt) then
1828! write(0,*)' pgam=',pgam(1,k), ' lamc=',lamc(1,k),' k=',k
1829! endif
1830
1831
1832 !========================================================================
1833 ! autoconversion of cloud liquid water to rain
1834 ! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
1835 ! minimum qc of 1 x 10^-8 prevents floating point error
1836
1837 if (.not. do_sb_physics) then
1838 call kk2000_liq_autoconversion(microp_uniform, qcic(:,k), &
1839 ncic(:,k), rho(:,k), relvar(:,k), prc(:,k), nprc(:,k), nprc1(:,k), mgncol)
1840 endif
1841
1842 ! assign qric based on prognostic qr, using assumed precip fraction
1843 ! note: this could be moved above for consistency with qcic and qiic calculations
1844 do i=1,mgncol
1845 if (precip_frac(i,k) > mincld) then
1846 dum = one / precip_frac(i,k)
1847 else
1848 dum = zero
1849 endif
1850! qric(i,k) = min(qr(i,k)*dum, 0.01_r8) ! limit in-precip mixing ratios to 10 g/kg
1851 qric(i,k) = min(qr(i,k)*dum, 0.05_r8) ! limit in-precip mixing ratios to 50 g/kg
1852 nric(i,k) = nr(i,k) * dum
1853
1854
1855 ! add autoconversion to precip from above to get provisional rain mixing ratio
1856 ! and number concentration (qric and nric)
1857
1858 if(qric(i,k) < qsmall) then
1859 qric(i,k) = zero
1860 nric(i,k) = zero
1861 endif
1862
1863 ! make sure number concentration is a positive number to avoid
1864 ! taking root of negative later
1865
1866 nric(i,k) = max(nric(i,k),zero)
1867 enddo
1868 ! Get size distribution parameters for cloud ice
1869
1870 call size_dist_param_ice(mg_ice_props, qiic(:,k), niic(:,k), &
1871 lami(:,k), mgncol, n0=n0i(:,k))
1872
1873 ! Alternative autoconversion
1874 if (do_sb_physics) then
1875 if (do_liq_liu) then
1876 call liu_liq_autoconversion(pgam(:,k),qcic(:,k),ncic(:,k), &
1877 qric(:,k),rho(:,k),relvar(:,k),prc(:,k),nprc(:,k),nprc1(:,k),mgncol)
1878 else
1879 call sb2001v2_liq_autoconversion(pgam(:,k),qcic(:,k),ncic(:,k), &
1880 qric(:,k),rho(:,k),relvar(:,k),prc(:,k),nprc(:,k),nprc1(:,k), mgncol)
1881 endif
1882 endif
1883
1884 !.......................................................................
1885 ! Autoconversion of cloud ice to snow
1886 ! similar to Ferrier (1994)
1887
1888 if (do_cldice) then
1889 do i=1,mgncol
1890 if (qiic(i,k) >= qimax) then
1891! if (qi(i,k) >= qimax) then
1892 ts_au_loc(i) = ts_au_min
1893 elseif (qiic(i,k) <= qimin) then
1894! elseif (qi(i,k) <= qimin) then
1895 ts_au_loc(i) = ts_au
1896 else
1897! ts_au_loc(i) = (ts_au*(qimax-qi(i,k)) + ts_au_min*(qi(i,k)-qimin)) * qiinv
1898 ts_au_loc(i) = (ts_au*(qimax-qiic(i,k)) + ts_au_min*(qiic(i,k)-qimin)) * qiinv
1899! ts_au_loc(i) = ts_au * exp(-tsfac*(qiic(i,k)-qimin))
1900 endif
1901! if (ts_au_loc(i) > ts_au_min) ts_au_loc(i) = ts_au_loc(i)*min(five,sqrt(p(i,nlev)/p(i,k)))
1902 enddo
1903! if (lprnt) write(0,*)' ts_au_loc=',ts_au_loc(1),' k=',k, ' qiic=',qiic(1,k),&
1904! if (lprnt) write(0,*)' ts_au_loc=',ts_au_loc(1),' k=',k, ' qi=',qi(1,k),&
1905! ' ts_au=',ts_au,' ts_au_min=',ts_au_min,' qimin=',qimin,' qimax=',qimax
1906! ' ts_au=',ts_au,' ts_au_min=',ts_au_min,' tsfac=',tsfac,' qimin=',qimin,' qimax=',qimax
1907
1908 if(do_ice_gmao) then
1909 call gmao_ice_autoconversion(t(:,k), qiic(:,k), niic(:,k), lami(:,k), &
1910 n0i(:,k), dcs, ts_au_loc(:), prci(:,k), nprci(:,k), mgncol)
1911 else
1912 call ice_autoconversion(t(:,k), qiic(:,k), lami(:,k), n0i(:,k), &
1913 dcs, ts_au_loc(:), prci(:,k), nprci(:,k), mgncol)
1914 endif
1915 !else
1916 ! Add in the particles that we have already converted to snow, and
1917 ! don't do any further autoconversion of ice.
1918 !prci(:,k) = tnd_qsnow(:,k) / cldm(:,k)
1919 !nprci(:,k) = tnd_nsnow(:,k) / cldm(:,k)
1920 endif
1921
1922 ! note, currently we don't have this
1923 ! inside the do_cldice block, should be changed later
1924 ! assign qsic based on prognostic qs, using assumed precip fraction
1925 do i=1,mgncol
1926 if (precip_frac(i,k) > mincld) then
1927 dum = one / precip_frac(i,k)
1928 else
1929 dum = zero
1930 endif
1931! qsic(i,k) = min(qs(i,k)*dum, 0.01_r8) ! limit in-precip mixing ratios to 10 g/kg
1932 qsic(i,k) = min(qs(i,k)*dum, 0.05_r8) ! limit in-precip mixing ratios to 50 g/kg
1933 nsic(i,k) = ns(i,k) * dum
1934
1935 ! if precip mix ratio is zero so should number concentration
1936
1937 if(qsic(i,k) < qsmall) then
1938 qsic(i,k) = zero
1939 nsic(i,k) = zero
1940 endif
1941
1942 ! make sure number concentration is a positive number to avoid
1943 ! taking root of negative later
1944
1945 nsic(i,k) = max(nsic(i,k), zero)
1946
1947!++ also do this for graupel, which is assumed to be 'precip_frac'
1948 qgic(i,k) = min(qg(i,k)*dum, 0.01_r8) ! limit in-precip mixing ratios to 10 g/kg)
1949 ngic(i,k) = ng(i,k) * dum
1950
1951 ! if precip mix ratio is zero so should number concentration
1952 if (qgic(i,k) < qsmall) then
1953 qgic(i,k) = zero
1954 ngic(i,k) = zero
1955 endif
1956
1957 ! make sure number concentration is a positive number to avoid
1958 ! taking root of negative later
1959
1960 ngic(i,k) = max(ngic(i,k), zero)
1961!--ag
1962 enddo
1963
1964 !.......................................................................
1965 ! get size distribution parameters for precip
1966 !......................................................................
1967 ! rain
1968
1969 call size_dist_param_basic(mg_rain_props, qric(:,k), nric(:,k), &
1970 lamr(:,k), mgncol, n0=n0r(:,k))
1971
1972 do i=1,mgncol
1973 if (lamr(i,k) >= qsmall) then
1974 dum = arn(i,k) / lamr(i,k)**br
1975 dum1 = 9.1_r8*rhof(i,k)
1976
1977 ! provisional rain number and mass weighted mean fallspeed (m/s)
1978
1979 umr(i,k) = min(dum1, dum*gamma_br_plus4*oneo6)
1980 unr(i,k) = min(dum1, dum*gamma_br_plus1)
1981 else
1982
1983 umr(i,k) = zero
1984 unr(i,k) = zero
1985 endif
1986 enddo
1987
1988 !......................................................................
1989 ! snow
1990
1991 call size_dist_param_basic(mg_snow_props, qsic(:,k), nsic(:,k), &
1992 lams(:,k), mgncol, n0=n0s(:,k))
1993
1994 do i=1,mgncol
1995 if (lams(i,k) >= qsmall) then
1996
1997 ! provisional snow number and mass weighted mean fallspeed (m/s)
1998
1999 dum = asn(i,k) / lams(i,k)**bs
2000 dum1 = 1.2_r8*rhof(i,k)
2001 ums(i,k) = min(dum1, dum*gamma_bs_plus4*oneo6)
2002 uns(i,k) = min(dum1, dum*gamma_bs_plus1)
2003
2004 else
2005 ums(i,k) = zero
2006 uns(i,k) = zero
2007 endif
2008 enddo
2009
2010 if (do_graupel .or. do_hail) then
2011!++ag
2012!use correct bg or bh (bgtmp=bg or bh)
2013 !......................................................................
2014 ! graupel/hail
2015
2016!++AG SET rhog here and for mg_graupel_props?
2017! For now: rhog is constant. Set to same in micro_mg_utils.F90
2018! Ideally: find a method to set once. (Hail = 400, Graupel = 500 from M2005)
2019
2020!mg,snow_props or mg_graupel props?
2021
2022 call size_dist_param_basic(mg_graupel_props, qgic(:,k), ngic(:,k), &
2023 lamg(:,k), mgncol, n0=n0g(:,k))
2024
2025 do i=1,mgncol
2026 if (lamg(i,k) >= qsmall) then
2027
2028 ! provisional graupel/hail number and mass weighted mean fallspeed (m/s)
2029
2030 dum = agn(i,k) / lamg(i,k)**bgtmp
2031 dum1 = 20._r8*rhof(i,k)
2032 umg(i,k) = min(dum1, dum*gamma_bg_plus4*oneo6)
2033 ung(i,k) = min(dum1, dum*gamma_bg_plus1)
2034! umg(i,k) = min(dum1, dum*gamma(four+bgtmp)*oneo6)
2035! ung(i,k) = min(dum1, dum*gamma(one+bgtmp))
2036
2037 else
2038 umg(i,k) = zero
2039 ung(i,k) = zero
2040 endif
2041 enddo
2042!--ag
2043 endif
2044
2045 if (do_cldice) then
2046 if (.not. use_hetfrz_classnuc) then
2047
2048 ! heterogeneous freezing of cloud water
2049 !----------------------------------------------
2050
2051! if (lprnt .and. k>=60 .and. k<=65) then
2052! if (lprnt .and. k>=100) then
2053! if (lprnt) then
2054! write(0,*)' pgam=',pgam(1,k), ' lamc=',lamc(1,k),' qcic=',qcic(1,k),' ncic=',ncic(1,k),' t=',t(1,k),' k=',k,&
2055! ' relvar=',relvar(1,k)
2056! endif
2057
2058 call immersion_freezing(microp_uniform, t(:,k), pgam(:,k), lamc(:,k), &
2059 qcic(:,k), ncic(:,k), relvar(:,k), mnuccc(:,k), nnuccc(:,k), mgncol)
2060
2061! if (lprnt .and. k>=60 .and. k<=65) then
2062! if (lprnt .and. k>=100) then
2063! if (lprnt) then
2064! write(0,*)' mnuccca=',mnuccc(1,k),' lcldm=',lcldm(1,k),' nnuccc=',nnuccc(1,k),' k=',k
2065! endif
2066
2067 ! make sure number of droplets frozen does not exceed available ice nuclei concentration
2068 ! this prevents 'runaway' droplet freezing
2069
2070! where (qcic(:,k) >= qsmall .and. t(:,k) < 269.15_r8 .and. lcldm(:,k) > mincld)
2071 where (qcic(:,k) >= qsmall .and. t(:,k) < 269.15_r8)
2072 where (nnuccc(:,k)*lcldm(:,k) > nnuccd(:,k))
2073 ! scale mixing ratio of droplet freezing with limit
2074 mnuccc(:,k) = mnuccc(:,k)*(nnuccd(:,k)/(nnuccc(:,k)*lcldm(:,k)))
2075 nnuccc(:,k) = nnuccd(:,k)/lcldm(:,k)
2076 end where
2077 end where
2078
2079! if (lprnt .and. k >= 60 .and. k <=65) write(0,*)' mnuccc=',mnuccc(1,60:65)
2080! if (lprnt .and. k >= 100) write(0,*)' mnuccc=',mnuccc(1,k)
2081! if (lprnt) write(0,*)' mnuccc=',mnuccc(1,k)
2082
2083 mdust = size(rndst,3)
2084 call contact_freezing(microp_uniform, t(:,k), p(:,k), rndst(:,k,:), &
2085 nacon(:,k,:), pgam(:,k), lamc(:,k), qcic(:,k), ncic(:,k), &
2086 relvar(:,k), mnucct(:,k), nnucct(:,k), mgncol, mdust)
2087
2088! if (lprnt .and. k >= 60 .and. k <=65) write(0,*)' mnucct=',mnucct(1,:)
2089! if (lprnt .and. k >= 100 ) write(0,*)' mnucct=',mnucct(1,k)
2090! if (lprnt) write(0,*)' mnucct=',mnucct(1,k)
2091
2092 mnudep(:,k) = zero
2093 nnudep(:,k) = zero
2094
2095 !else
2096
2097 ! Mass of droplets frozen is the average droplet mass, except
2098 ! with two limiters: concentration must be at least 1/cm^3, and
2099 ! mass must be at least the minimum defined above.
2100 !mi0l = qcic(:,k)/max(ncic(:,k), 1.0e6_r8/rho(:,k))
2101 !mi0l = max(mi0l_min, mi0l)
2102
2103 !where (qcic(:,k) >= qsmall)
2104 !nnuccc(:,k) = frzimm(:,k)*1.0e6_r8*rhoinv(:,k)
2105 !mnuccc(:,k) = nnuccc(:,k)*mi0l
2106
2107 !nnucct(:,k) = frzcnt(:,k)*1.0e6_r8*rhoinv(:,k)
2108 !mnucct(:,k) = nnucct(:,k)*mi0l
2109
2110 !nnudep(:,k) = frzdep(:,k)*1.0e6_r8*rhoinv(:,k)
2111 !mnudep(:,k) = nnudep(:,k)*mi0
2112 !elsewhere
2113 !nnuccc(:,k) = zero
2114 !mnuccc(:,k) = zero
2115
2116 !nnucct(:,k) = zero
2117 !mnucct(:,k) = zero
2118
2119 !nnudep(:,k) = zero
2120 !mnudep(:,k) = zero
2121 !end where
2122
2123 endif
2124
2125 else
2126 do i=1,mgncol
2127 mnuccc(i,k) = zero
2128 nnuccc(i,k) = zero
2129 mnucct(i,k) = zero
2130 nnucct(i,k) = zero
2131 mnudep(i,k) = zero
2132 nnudep(i,k) = zero
2133 enddo
2134 endif
2135
2136 call snow_self_aggregation(t(:,k), rho(:,k), asn(:,k), rhosn, qsic(:,k), nsic(:,k), &
2137 nsagg(:,k), mgncol)
2138
2139 call accrete_cloud_water_snow(t(:,k), rho(:,k), asn(:,k), uns(:,k), mu(:,k), &
2140 qcic(:,k), ncic(:,k), qsic(:,k), pgam(:,k), lamc(:,k), lams(:,k), n0s(:,k), &
2141 psacws(:,k), npsacws(:,k), mgncol)
2142
2143 if (do_cldice) then
2144 call secondary_ice_production(t(:,k), psacws(:,k), msacwi(:,k), nsacwi(:,k), mgncol)
2145 else
2146 nsacwi(:,k) = zero
2147 msacwi(:,k) = zero
2148 endif
2149
2150 call accrete_rain_snow(t(:,k), rho(:,k), umr(:,k), ums(:,k), unr(:,k), uns(:,k), &
2151 qric(:,k), qsic(:,k), lamr(:,k), n0r(:,k), lams(:,k), n0s(:,k), &
2152 pracs(:,k), npracs(:,k), mgncol)
2153
2154 call heterogeneous_rain_freezing(t(:,k), qric(:,k), nric(:,k), lamr(:,k), &
2155 mnuccr(:,k), nnuccr(:,k), mgncol)
2156
2157 if (do_sb_physics) then
2158 call sb2001v2_accre_cld_water_rain(qcic(:,k), ncic(:,k), qric(:,k), &
2159 rho(:,k), relvar(:,k), pra(:,k), npra(:,k), mgncol)
2160 else
2161 call accrete_cloud_water_rain(microp_uniform, qric(:,k), qcic(:,k), &
2162 ncic(:,k), relvar(:,k), accre_enhan(:,k), pra(:,k), npra(:,k), mgncol)
2163 endif
2164
2165 call self_collection_rain(rho(:,k), qric(:,k), nric(:,k), nragg(:,k), mgncol)
2166
2167 if (do_cldice) then
2168 call accrete_cloud_ice_snow(t(:,k), rho(:,k), asn(:,k), qiic(:,k), niic(:,k), &
2169 qsic(:,k), lams(:,k), n0s(:,k), prai(:,k), nprai(:,k), mgncol)
2170 else
2171 prai(:,k) = zero
2172 nprai(:,k) = zero
2173 endif
2174
2175!++ag Moved below graupel conditional, now two different versions
2176! if (.not. (do_hail .or. do_graupel)) then
2177! call evaporate_sublimate_precip(t(:,k), rho(:,k), &
2178! dv(:,k), mu(:,k), sc(:,k), q(:,k), qvl(:,k), qvi(:,k), &
2179! lcldm(:,k), precip_frac(:,k), arn(:,k), asn(:,k), qcic(:,k), qiic(:,k), &
2180! qric(:,k), qsic(:,k), lamr(:,k), n0r(:,k), lams(:,k), n0s(:,k), &
2181! pre(:,k), prds(:,k), am_evp_st(:,k), mgncol)
2182! endif
2183!--ag
2184
2185 call bergeron_process_snow(t(:,k), rho(:,k), dv(:,k), mu(:,k), sc(:,k), &
2186 qvl(:,k), qvi(:,k), asn(:,k), qcic(:,k), qsic(:,k), lams(:,k), n0s(:,k), &
2187 bergs(:,k), mgncol)
2188! if(lprnt) write(0,*)' bergs1=',bergs(1,k),' k=',k,' micro_mg_berg_eff_factor=',micro_mg_berg_eff_factor
2189! if(lprnt) write(0,*)' t=',t(1,k),' rho=',rho(1,k),' dv=',dv(1,k),' mu=',mu(1,k),&
2190! 'qcic=',qcic(1,k),' qsic=',qsic(1,k),' qvl=',qvl(1,k),' qvi=',qvi(1,k), &
2191! ' mu=',mu(1,k),' sc=',sc(1,k),' asn=',asn(1,k),' lams=',lams(1,k),' n0s=',n0s(1,k),' ni=',ni(1,k)
2192
2193 bergs(:,k) = bergs(:,k) * micro_mg_berg_eff_factor
2194
2195 !+++PMC 12/3/12 - NEW VAPOR DEP/SUBLIMATION GOES HERE!!!
2196 if (do_cldice) then
2197
2198 call ice_deposition_sublimation(t(:,k), q(:,k), qi(:,k), ni(:,k), &
2199 icldm(:,k), rho(:,k), dv(:,k), qvl(:,k), qvi(:,k), &
2200 berg(:,k), vap_dep(:,k), ice_sublim(:,k), mgncol)
2201
2202! if(lprnt) write(0,*)' t=',t(1,k),' k=',k,' q=',q(1,k),' qi=',qi(1,k),&
2203! ' ni=',ni(1,k),' icldm=',icldm(1,k),' rho=',rho(1,k),' dv=',dv(1,k),&
2204! ' qvl=',qvl(1,k),' qvi=',qvi(1,k),' berg=',berg(1,k),' vap_dep=',&
2205! vap_dep(1,k),' ice_sublim=',ice_sublim(1,k)
2206! if(lprnt) write(0,*)' berg1=',berg(1,k),' k=',k,' micro_mg_berg_eff_factor=',micro_mg_berg_eff_factor
2207 do i=1,mgncol
2208! sublimation should not exceed available ice
2209 ice_sublim(i,k) = max(ice_sublim(i,k), -qi(i,k)*oneodt)
2210 berg(i,k) = berg(i,k) * micro_mg_berg_eff_factor
2211 if (ice_sublim(i,k) < zero .and. qi(i,k) > qsmall .and. icldm(i,k) > mincld) then
2212 nsubi(i,k) = sublim_factor * ice_sublim(i,k) * ni(i,k) / (qi(i,k) * icldm(i,k))
2213 else
2214 nsubi(i,k) = zero
2215 endif
2216
2217 ! bergeron process should not reduce nc unless
2218 ! all ql is removed (which is handled elsewhere)
2219 !in fact, nothing in this entire file makes nsubc nonzero.
2220 nsubc(i,k) = zero
2221 enddo
2222
2223 endif !do_cldice
2224 !---PMC 12/3/12
2225
2226!++ag Process rate calls for graupel here.
2227! (Should this be in do_cldice loop?)
2228!===================================================================
2229
2230 if(do_hail .or. do_graupel) then
2231 call graupel_collecting_snow(qsic(:,k),qric(:,k),umr(:,k),ums(:,k), &
2232 rho(:,k),lamr(:,k),n0r(:,k),lams(:,k),n0s(:,k), psacr(:,k), mgncol)
2233
2234 call graupel_collecting_cld_water(qgic(:,k),qcic(:,k),ncic(:,k),rho(:,k), &
2235 n0g(:,k),lamg(:,k),bgtmp,agn(:,k), psacwg(:,k), npsacwg(:,k), mgncol)
2236
2237 call graupel_riming_liquid_snow(psacws(:,k),qsic(:,k),qcic(:,k),nsic(:,k), &
2238 rho(:,k),rhosn,rhogtmp,asn(:,k),lams(:,k),n0s(:,k),deltat, &
2239 pgsacw(:,k),nscng(:,k),mgncol)
2240
2241! if(lprnt .and. k >=100) then
2242! if(lprnt) then
2243! write(0,*)' k=',k,' qric=',qric(1,k),' qgic=',qgic(1,k),' umg=',umg(1,k),' umr=',umr(1,k),&
2244! ' ung=',ung(1,k),' unr=',unr(1,k),' rho=',rho(1,k),' n0r=',n0r(1,k),' lamr=',lamr(1,k),&
2245! ' n0g=',n0g(1,k),' lamg=',lamg(1,k),' pracg=',pracg(1,k)
2246! endif
2247 call graupel_collecting_rain(qric(:,k),qgic(:,k),umg(:,k), &
2248 umr(:,k),ung(:,k),unr(:,k),rho(:,k),n0r(:,k),lamr(:,k),n0g(:,k), &
2249 lamg(:,k), pracg(:,k),npracg(:,k),mgncol)
2250! if(lprnt .and. k >=100) write(0,*)' k=',k,' pracg=',pracg(1,k),' npracg=',npracg(1,k)
2251
2252!AG note: Graupel rain riming snow changes
2253! pracs, npracs, (accretion of rain by snow) psacr (collection of snow by rain)
2254
2255! if (lprnt .and. abs(k-81) <5) &
2256! write(0,*)' k=',k,' pracs=',pracs(1,k),' npracs=',npracs(1,k),' psacr=',psacr(1,k),&
2257! ' qsic=',qsic(1,k),' qric=',qric(1,k),' nric=',nric(1,k),' nsic=',nsic(1,k), &
2258! ' n0s=',n0s(1,k),' lams=',lams(1,k),' n0r=',n0r(1,k),' lamr=',lamr(1,k), &
2259! ' pgracs=',pgracs(1,k),' ngracs=',ngracs(1,k)
2260
2261 call graupel_rain_riming_snow(pracs(:,k),npracs(:,k),psacr(:,k),qsic(:,k), &
2262 qric(:,k),nric(:,k),nsic(:,k),n0s(:,k),lams(:,k),n0r(:,k),lamr(:,k), &
2263 deltat,pgracs(:,k),ngracs(:,k),mgncol)
2264! if (lprnt .and. abs(k-81) <5) &
2265! write(0,*)' k=',k,' pracs=',pracs(1,k),' npracs=',npracs(1,k),' psacr=',psacr(1,k),&
2266! ' pgracs=',pgracs(1,k),' ngracs=',ngracs(1,k)
2267
2268 call graupel_rime_splintering(t(:,k),qcic(:,k),qric(:,k),qgic(:,k), &
2269 psacwg(:,k),pracg(:,k),qmultg(:,k),nmultg(:,k),qmultrg(:,k), &
2270 nmultrg(:,k),mgncol)
2271
2272! if(lprnt .and. k >=100) write(0,*)' k=',k,' pracg2=',pracg(1,k)
2273! if (lprnt .and. abs(k-81) <5) &
2274! write(0,*)' k=',k,' pracg2=',pracg(1,k)
2275
2276 call evaporate_sublimate_precip_graupel(t(:,k), rho(:,k), &
2277 dv(:,k), mu(:,k), sc(:,k), q(:,k), qvl(:,k), qvi(:,k), &
2278 lcldm(:,k), precip_frac(:,k), arn(:,k), asn(:,k), agn(:,k), bgtmp, &
2279 qcic(:,k), qiic(:,k), qric(:,k), qsic(:,k), qgic(:,k), &
2280 lamr(:,k), n0r(:,k), lams(:,k), n0s(:,k), lamg(:,k), n0g(:,k), &
2281 pre(:,k), prds(:,k), prdg(:,k), am_evp_st(:,k), mgncol)
2282
2283!!Not used: part of above
2284!! call graupel_sublimate_evap(t(:,k),q(:,k),qgic(:,k),rho(:,k),n0g(:,k), &
2285!! lamg(:,k),qvi(:,k),dv(:,k),mu(:,k),sc(:,k),bgtmp,agn(:,k), &
2286!! prdg(:,k),eprdg(:,k),mgncol)
2287
2288!Checks for Debugging
2289
2290! if (minval(qmultg(:,k)).lt.0._r8) &
2291! write(iulog,*) "OOPS, qmultg < 0 : min=",minval(qmultg(:,k))
2292!
2293! if (minval(qmultrg(:,k)).lt.0._r8) &
2294! write(iulog,*) "OOPS, qmultrg < 0 : min=",minval(qmultrg(:,k))
2295!
2296! if (minval(pgracs(:,k)).lt.0._r8) &
2297! write(iulog,*) "OOPS, pgracs < 0 : min=",minval(pgracs(:,k))
2298!
2299! if (minval(psacwg(:,k)).lt.0._r8) &
2300! write(iulog,*) "OOPS, psacwg < 0 : min=",minval(psacwg(:,k))
2301!
2302! if (minval(npsacwg(:,k)).lt.0._r8) &
2303! write(iulog,*) "OOPS, npsacwg < 0 : min=",minval(npsacwg(:,k))
2304!
2305! if (minval(pracg(:,k)).lt.0._r8) &
2306! write(iulog,*) "OOPS, pracg < 0 : min=",minval(pracg(:,k))
2307!
2308! if (maxval(prdg(:,k)).gt.0._r8) &
2309! write(iulog,*) "OOPS, prdg > 0 : max=",maxval(prdg(:,k))
2310!
2311! if (minval(nmultg(:,k)).lt.0._r8) &
2312! write(iulog,*) "OOPS, nmultg < 0 : min=",minval(nmultg(:,k))
2313!
2314! if (minval(nmultrg(:,k)).lt.0._r8) &
2315! write(iulog,*) "OOPS, nmultrg < 0 : min=",minval(nmultrg(:,k))
2316!
2317! if (minval(ngracs(:,k)).lt.0._r8) &
2318! write(iulog,*) "OOPS, ngracs < 0 : min=",minval(ngracs(:,k))
2319!
2320! if (minval(psacr(:,k)).lt.0._r8) &
2321! write(iulog,*) "OOPS, psacr < 0 : min=",minval(psacr(:,k))
2322!
2323! if (minval(nscng(:,k)).lt.0._r8) &
2324! write(iulog,*) "OOPS, nscng < 0 : min=",minval(nscng(:,k))
2325
2326 else
2327! Routine without Graupel (original)
2328 call evaporate_sublimate_precip(t(:,k), rho(:,k), &
2329 dv(:,k), mu(:,k), sc(:,k), q(:,k), qvl(:,k), qvi(:,k), &
2330 lcldm(:,k), precip_frac(:,k), arn(:,k), asn(:,k), qcic(:,k), qiic(:,k), &
2331 qric(:,k), qsic(:,k), lamr(:,k), n0r(:,k), lams(:,k), n0s(:,k), &
2332 pre(:,k), prds(:,k), am_evp_st(:,k), mgncol)
2333
2334
2335 endif ! end do_graupel/hail loop
2336!--ag
2337
2338 do i=1,mgncol
2339
2340 ! conservation to ensure no negative values of cloud water/precipitation
2341 ! in case microphysical process rates are large
2342 !===================================================================
2343
2344 ! note: for check on conservation, processes are multiplied by omsm
2345 ! to prevent problems due to round off error
2346
2347 ! conservation of qc
2348 !-------------------------------------------------------------------
2349
2350!++ag Add graupel tendencies for qc to equation ON
2351! dum = ((prc(i,k)+pra(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k)+ &
2352! psacws(i,k)+bergs(i,k))*lcldm(i,k)+berg(i,k))*deltat
2353 dum = ( (prc(i,k) + pra(i,k) + mnuccc(i,k) + mnucct(i,k) + msacwi(i,k) &
2354 + psacws(i,k) + bergs(i,k) + qmultg(i,k) + psacwg(i,k) + pgsacw(i,k))*lcldm(i,k) &
2355 + berg(i,k) ) * deltat
2356!--ag
2357
2358 if (dum > qc(i,k) .and. abs(dum) > qsmall) then
2359!++ag
2360! ratio = qc(i,k)/deltat/((prc(i,k)+pra(i,k)+mnuccc(i,k)+mnucct(i,k)+ &
2361! msacwi(i,k)+psacws(i,k)+bergs(i,k))*lcldm(i,k)+berg(i,k))*omsm
2362
2363 ratio = qc(i,k) / dum * omsm
2364
2365 qmultg(i,k) = ratio * qmultg(i,k)
2366 psacwg(i,k) = ratio * psacwg(i,k)
2367 pgsacw(i,k) = ratio * pgsacw(i,k)
2368!--ag
2369 prc(i,k) = ratio * prc(i,k)
2370 pra(i,k) = ratio * pra(i,k)
2371 mnuccc(i,k) = ratio * mnuccc(i,k)
2372 mnucct(i,k) = ratio * mnucct(i,k)
2373 msacwi(i,k) = ratio * msacwi(i,k)
2374 psacws(i,k) = ratio * psacws(i,k)
2375 bergs(i,k) = ratio * bergs(i,k)
2376 berg(i,k) = ratio * berg(i,k)
2377 qcrat(i,k) = ratio
2378 else
2379 qcrat(i,k) = one
2380 endif
2381
2382! if(lprnt) write(0,*)' bergs2=',bergs(1,k),' k=',k,' ratio=',ratio
2383
2384 !PMC 12/3/12: ratio is also frac of step w/ liquid.
2385 !thus we apply berg for "ratio" of timestep and vapor
2386 !deposition for the remaining frac of the timestep.
2387 if (qc(i,k) >= qsmall) then
2388 vap_dep(i,k) = vap_dep(i,k) * (one-qcrat(i,k))
2389 endif
2390
2391 enddo
2392
2393 do i=1,mgncol
2394
2395 !=================================================================
2396 ! apply limiter to ensure that ice/snow sublimation and rain evap
2397 ! don't push conditions into supersaturation, and ice deposition/nucleation don't
2398 ! push conditions into sub-saturation
2399 ! note this is done after qc conservation since we don't know how large
2400 ! vap_dep is before then
2401 ! estimates are only approximate since other process terms haven't been limited
2402 ! for conservation yet
2403
2404 ! first limit ice deposition/nucleation vap_dep + mnuccd
2405
2406 dum1 = vap_dep(i,k) + mnuccd(i,k)
2407 if (dum1 > 1.e-20_r8) then
2408 dum = (q(i,k)-qvi(i,k))/(one + xxls_squared*qvi(i,k)/(cpp*rv*t(i,k)*t(i,k)))*oneodt
2409 dum = max(dum, zero)
2410 if (dum1 > dum) then
2411 ! Allocate the limited "dum" tendency to mnuccd and vap_dep
2412 ! processes. Don't divide by cloud fraction; these are grid-
2413 ! mean rates.
2414 dum1 = mnuccd(i,k) / (vap_dep(i,k)+mnuccd(i,k))
2415 mnuccd(i,k) = dum*dum1
2416 vap_dep(i,k) = dum - mnuccd(i,k)
2417 endif
2418 endif
2419
2420 enddo
2421
2422 do i=1,mgncol
2423
2424 !===================================================================
2425 ! conservation of nc
2426 !-------------------------------------------------------------------
2427!++ag NEW ONE ON
2428! dum = (nprc1(i,k)+npra(i,k)+nnuccc(i,k)+nnucct(i,k)+ &
2429! npsacws(i,k)-nsubc(i,k))*lcldm(i,k)*deltat
2430 dum = (nprc1(i,k) + npra(i,k) + nnuccc(i,k) + nnucct(i,k) &
2431 + npsacws(i,k) - nsubc(i,k) + npsacwg(i,k))*lcldm(i,k)*deltat
2432!--ag
2433
2434 if (dum > nc(i,k) .and. abs(dum) > qsmall) then
2435 ratio = nc(i,k) / dum * omsm
2436!++ag
2437 npsacwg(i,k) = ratio * npsacwg(i,k)
2438!--ag
2439
2440 nprc1(i,k) = ratio * nprc1(i,k)
2441 npra(i,k) = ratio * npra(i,k)
2442 nnuccc(i,k) = ratio * nnuccc(i,k)
2443 nnucct(i,k) = ratio * nnucct(i,k)
2444 npsacws(i,k) = ratio * npsacws(i,k)
2445 nsubc(i,k) = ratio * nsubc(i,k)
2446 endif
2447
2448 mnuccri(i,k) = zero
2449 nnuccri(i,k) = zero
2450
2451 if (do_cldice) then
2452
2453 ! freezing of rain to produce ice if mean rain size is smaller than Dcs
2454 if (lamr(i,k) > qsmall) then
2455 if (one/lamr(i,k) < dcs) then
2456 mnuccri(i,k) = mnuccr(i,k)
2457 nnuccri(i,k) = nnuccr(i,k)
2458 mnuccr(i,k) = zero
2459 nnuccr(i,k) = zero
2460 endif
2461 endif
2462 endif
2463
2464 enddo
2465
2466 do i=1,mgncol
2467
2468 ! conservation of rain mixing ratio
2469 !-------------------------------------------------------------------
2470!++ag Implemented change for graupel
2471 dum1 = - pre(i,k) + pracs(i,k) + mnuccr(i,k) + mnuccri(i,k) &
2472 + qmultrg(i,k) + pracg(i,k) + pgracs(i,k)
2473 dum3 = dum1 * precip_frac(i,k)
2474 dum2 = (pra(i,k)+prc(i,k))*lcldm(i,k)
2475 dum = (dum3 - dum2) * deltat
2476!--ag
2477
2478 ! note that qrtend is included below because of instantaneous freezing/melt
2479 if (dum > qr(i,k) .and. dum1 >= qsmall .and. abs(dum3) > qsmall) then
2480 ratio = (qr(i,k)*oneodt + dum2) / dum3 * omsm
2481!++ag
2482 qmultrg(i,k) = ratio * qmultrg(i,k)
2483 pracg(i,k) = ratio * pracg(i,k)
2484 pgracs(i,k) = ratio * pgracs(i,k)
2485!--ag
2486 pre(i,k) = ratio * pre(i,k)
2487 pracs(i,k) = ratio * pracs(i,k)
2488 mnuccr(i,k) = ratio * mnuccr(i,k)
2489 mnuccri(i,k) = ratio * mnuccri(i,k)
2490 endif
2491
2492 enddo
2493
2494 do i=1,mgncol
2495
2496 ! conservation of rain number
2497 !-------------------------------------------------------------------
2498
2499 ! Add evaporation of rain number.
2500 if (pre(i,k) < zero) then
2501 dum = max(-one, pre(i,k)*deltat/qr(i,k))
2502 nsubr(i,k) = dum*nr(i,k) * oneodt
2503 else
2504 nsubr(i,k) = zero
2505 endif
2506
2507 enddo
2508
2509 do i=1,mgncol
2510
2511!++ag IMplemented change for graupel
2512! dum1 = (-nsubr(i,k)+npracs(i,k)+nnuccr(i,k)+nnuccri(i,k)-nragg(i,k))*precip_frac(i,k)
2513! nprc(i,k)*lcldm(i,k))*deltat
2514
2515 dum1 = (-nsubr(i,k)+npracs(i,k)+nnuccr(i,k)+nnuccri(i,k)-nragg(i,k) &
2516 +npracg(i,k)+ngracs(i,k))*precip_frac(i,k)
2517 dum2 = nprc(i,k)*lcldm(i,k)
2518 dum = (dum1 - dum2) * deltat
2519!--ag
2520
2521 if (dum > nr(i,k) .and. abs(dum1) > qsmall) then
2522 ratio = (nr(i,k)*oneodt + dum2) / dum1 * omsm
2523
2524!++ag
2525 npracg(i,k) = ratio * npracg(i,k)
2526 ngracs(i,k) = ratio * ngracs(i,k)
2527!--ag
2528 nragg(i,k) = ratio * nragg(i,k)
2529 npracs(i,k) = ratio * npracs(i,k)
2530 nnuccr(i,k) = ratio * nnuccr(i,k)
2531 nsubr(i,k) = ratio * nsubr(i,k)
2532 nnuccri(i,k) = ratio * nnuccri(i,k)
2533 endif
2534
2535 enddo
2536
2537 if (do_cldice) then
2538
2539 do i=1,mgncol
2540
2541 ! conservation of qi
2542 !-------------------------------------------------------------------
2543
2544!++ag
2545
2546 dum1 = (prci(i,k)+prai(i,k))*icldm(i,k)-ice_sublim(i,k)
2547! dum2 = vap_dep(i,k)+berg(i,k)+mnuccd(i,k) &
2548! + (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k))*lcldm(i,k) &
2549! + mnuccri(i,k)*precip_frac(i,k)
2550 dum2 = vap_dep(i,k)+berg(i,k)+mnuccd(i,k) &
2551 + (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k)+qmultg(i,k))*lcldm(i,k) &
2552 + (qmultrg(i,k)+mnuccri(i,k))*precip_frac(i,k)
2553 dum = (dum1 - dum2) * deltat
2554!-ag
2555
2556 if (dum > qi(i,k) .and. abs(dum1) > qsmall) then
2557 ratio = (qi(i,k)*oneodt + dum2) / dum1 * omsm
2558
2559!++ag
2560! Only sink terms are limited.
2561! qmultg(i,k) = ratio * qmultg(i,k)
2562! qmultrg(i,k) = ratio * qmultrg(i,k)
2563!--ag
2564 prci(i,k) = ratio * prci(i,k)
2565 prai(i,k) = ratio * prai(i,k)
2566 ice_sublim(i,k) = ratio * ice_sublim(i,k)
2567 endif
2568
2569 enddo
2570
2571 endif
2572
2573 if (do_cldice) then
2574
2575 do i=1,mgncol
2576
2577 ! conservation of ni
2578 !-------------------------------------------------------------------
2579 if (use_hetfrz_classnuc) then
2580 tmpfrz = nnuccc(i,k)
2581 else
2582 tmpfrz = zero
2583 endif
2584!++ag
2585 dum1 = (nprci(i,k)+nprai(i,k)-nsubi(i,k))*icldm(i,k)
2586! dum2 = nnuccd(i,k)+(nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k))*lcldm(i,k) &
2587! + nnuccri(i,k)*precip_frac(i,k)
2588 dum2 = nnuccd(i,k)+(nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k)+nmultg(i,k))*lcldm(i,k) &
2589 + (nmultrg(i,k)+nnuccri(i,k))*precip_frac(i,k)
2590!--ag
2591 dum = (dum1 - dum2) * deltat
2592
2593 if (dum > ni(i,k) .and. abs(dum1) > qsmall) then
2594 ratio = (ni(i,k)*oneodt + dum2) / dum1 * omsm
2595
2596 nprci(i,k) = ratio * nprci(i,k)
2597 nprai(i,k) = ratio * nprai(i,k)
2598 nsubi(i,k) = ratio * nsubi(i,k)
2599 endif
2600
2601 enddo
2602
2603 endif
2604
2605 do i=1,mgncol
2606
2607 ! conservation of snow mixing ratio
2608 !-------------------------------------------------------------------
2609!++ag
2610 if (do_hail .or. do_graupel) then
2611!NOTE: mnuccr is moved to graupel when active
2612!psacr is a positive value, but a loss for snow
2613!HM: psacr is positive in dum (two negatives)
2614
2615 dum1 = (psacr(i,k) - prds(i,k)) * precip_frac(i,k)
2616 dum2 = pracs(i,k)*precip_frac(i,k) &
2617 + (prai(i,k)+prci(i,k))*icldm(i,k) + (bergs(i,k)+psacws(i,k))*lcldm(i,k)
2618 dum = (dum1 - dum2) * deltat
2619 if (dum > qs(i,k) .and. psacr(i,k)-prds(i,k) >= qsmall) then
2620 ratio = (qs(i,k)*oneodt + dum2) / dum1 * omsm
2621 psacr(i,k) = ratio * psacr(i,k)
2622 prds(i,k) = ratio * prds(i,k)
2623 endif
2624 else
2625 dum1 = - prds(i,k) * precip_frac(i,k)
2626 dum2 = (pracs(i,k)+mnuccr(i,k))*precip_frac(i,k) &
2627 + (prai(i,k)+prci(i,k))*icldm(i,k) + (bergs(i,k)+psacws(i,k))*lcldm(i,k)
2628 dum = (dum1 - dum2) * deltat
2629 if (dum > qs(i,k) .and. -prds(i,k) >= qsmall) then
2630 ratio = (qs(i,k)*oneodt + dum2) / dum1 * omsm
2631 prds(i,k) = ratio * prds(i,k)
2632 endif
2633 endif
2634
2635!--ag
2636! dum1 = - prds(i,k) * precip_frac(i,k)
2637! dum2 = (pracs(i,k)+mnuccr(i,k))*precip_frac(i,k) &
2638! + (prai(i,k)+prci(i,k))*icldm(i,k) + (bergs(i,k)+psacws(i,k))*lcldm(i,k)
2639
2640! dum = (dum1 - dum2) * deltat
2641
2642! if (dum > qs(i,k) .and. -prds(i,k) >= qsmall) then
2643! ratio = (qs(i,k)*oneodt + dum2) / dum1 * omsm
2644
2645! prds(i,k) = ratio * prds(i,k)
2646! endif
2647
2648 enddo
2649
2650 do i=1,mgncol
2651
2652 ! conservation of snow number
2653 !-------------------------------------------------------------------
2654 ! calculate loss of number due to sublimation
2655 ! for now neglect sublimation of ns
2656 nsubs(i,k) = zero
2657
2658 ratio = one
2659!++ag Watch sign of nscng and ngracs. What is sign of nnuccr? Negative? Should be a source here.
2660
2661 if (do_hail .or. do_graupel) then
2662! dum1 = precip_frac(i,k)* (-nsubs(i,k)-nsagg(i,k)+ngracs(i,k))
2663! dum2 = nprci(i,k)*icldm(i,k) + nscng(i,k)*lcldm(i,k)
2664! dum = (dum1 - dum2) * deltat
2665! check here - this is slightly different from ag version - moorthi
2666
2667 dum1 = precip_frac(i,k)* (-nsubs(i,k)-nsagg(i,k)+ngracs(i,k)) &
2668 - nscng(i,k)*lcldm(i,k)
2669 dum2 = nprci(i,k)*icldm(i,k)
2670 dum = (dum1 - dum2) * deltat
2671
2672 if (dum > ns(i,k) .and. abs(dum1) > qsmall) then
2673 ratio = (ns(i,k)*oneodt + dum2) / dum1 * omsm
2674 nscng(i,k) = ratio * nscng(i,k)
2675 ngracs(i,k) = ratio * ngracs(i,k)
2676 endif
2677
2678 else
2679 dum1 = precip_frac(i,k)* (-nsubs(i,k)-nsagg(i,k))
2680 dum2 = nnuccr(i,k)*precip_frac(i,k) + nprci(i,k)*icldm(i,k)
2681 dum = (dum1 - dum2) * deltat
2682
2683 if (dum > ns(i,k) .and. abs(dum1) > qsmall) then
2684 ratio = (ns(i,k)*oneodt + dum2) / dum1 * omsm
2685 end if
2686 endif
2687 nsubs(i,k) = ratio * nsubs(i,k)
2688 nsagg(i,k) = ratio * nsagg(i,k)
2689
2690 enddo
2691
2692!++ag Graupel Conservation Checks
2693!-------------------------------------------------------------------
2694 if (do_hail .or. do_graupel) then
2695! conservation of graupel mass
2696!-------------------------------------------------------------------
2697 do i=1,mgncol
2698
2699 dum1 = -prdg(i,k) * precip_frac(i,k)
2700 dum2 = (pracg(i,k)+pgracs(i,k)+psacr(i,k)+mnuccr(i,k))*precip_frac(i,k) &
2701 + (psacwg(i,k)+pgsacw(i,k))*lcldm(i,k)
2702 dum = (dum1 - dum2) * deltat
2703
2704 if (dum > qg(i,k) .and. abs(dum1) > qsmall) then
2705
2706! hm added
2707! note: prdg is always negative (like prds), so it needs to be subtracted in ratio
2708 ratio = (qg(i,k)*oneodt + dum2) / dum1 * omsm
2709
2710 prdg(i,k) = ratio * prdg(i,k)
2711
2712 endif
2713
2714 enddo
2715
2716! conservation of graupel number: not needed, no sinks
2717!-------------------------------------------------------------------
2718 endif
2719!--ag
2720
2721
2722 do i=1,mgncol
2723
2724 ! next limit ice and snow sublimation and rain evaporation
2725 ! get estimate of q and t at end of time step
2726 ! don't include other microphysical processes since they haven't
2727 ! been limited via conservation checks yet
2728
2729!++ag need to add graupel sublimation/evap here too (prdg)? May not need eprdg?
2730!++ag
2731 tx1 = pre(i,k) * precip_frac(i,k)
2732 tx2 = prds(i,k) * precip_frac(i,k)
2733 tx6 = prdg(i,k) * precip_frac(i,k)
2734 tx5 = tx2 + tx6
2735 tx3 = tx1 + tx5 + ice_sublim(i,k)
2736
2737 if (tx3 < -1.e-20_r8) then
2738
2739 tx4 = tx5 + ice_sublim(i,k) + vap_dep(i,k) + mnuccd(i,k)
2740 qtmp = q(i,k) - (tx1 + tx4) * deltat
2741 ttmp = t(i,k) + (tx1*xxlv + tx4*xxls) * (deltat/cpp)
2742
2743 ! use rhw to allow ice supersaturation
2744 ! call qsat_water(ttmp, p(i,k), esn, qvn)
2745 esn = min(fpvsl(ttmp), p(i,k))
2746 qvn = epsqs*esn/(p(i,k)-omeps*esn) * qsfm(i,k)
2747! qvn = epsqs*esn/(p(i,k)-omeps*esn)
2748
2749 ! modify ice/precip evaporation rate if q > qsat
2750 if (qtmp > qvn) then
2751
2752 tx4 = one / tx3
2753 dum1 = tx1 * tx4
2754 dum2 = tx2 * tx4
2755!++ag
2756 dum3 = tx6 * tx4
2757!--ag
2758 ! recalculate q and t after vap_dep and mnuccd but without evap or sublim
2759 tx5 = (vap_dep(i,k)+mnuccd(i,k)) * deltat
2760 qtmp = q(i,k) - tx5
2761 ttmp = t(i,k) + tx5 * (xxls/cpp)
2762
2763 ! use rhw to allow ice supersaturation
2764 !call qsat_water(ttmp, p(i,k), esn, qvn)
2765 esn = min(fpvsl(ttmp), p(i,k))
2766 qvn = epsqs*esn / (p(i,k)-omeps*esn) * qsfm(i,k)
2767! qvn = epsqs*esn / (p(i,k)-omeps*esn)
2768
2769 dum = min(zero, (qtmp-qvn)/(one + xxlv_squared*qvn/(cpp*rv*ttmp*ttmp)))
2770
2771 ! modify rates if needed, divide by precip_frac to get local (in-precip) value
2772 if (precip_frac(i,k) > mincld) then
2773 tx4 = oneodt / precip_frac(i,k)
2774 else
2775 tx4 = zero
2776 endif
2777 pre(i,k) = dum*dum1*tx4
2778
2779 ! do separately using RHI for prds and ice_sublim
2780 !call qsat_ice(ttmp, p(i,k), esn, qvn)
2781 esn = min(fpvsi(ttmp), p(i,k))
2782 qvn = epsqs*esn / (p(i,k)-omeps*esn) * qsfm(i,k)
2783! qvn = epsqs*esn / (p(i,k)-omeps*esn)
2784
2785
2786 dum = min(zero, (qtmp-qvn)/(one + xxls_squared*qvn/(cpp*rv*ttmp*ttmp)))
2787
2788 ! modify rates if needed, divide by precip_frac to get local (in-precip) value
2789 prds(i,k) = dum*dum2*tx4
2790!++ag
2791 prdg(i,k) = dum*dum3*tx4
2792!--ag
2793!++ag
2794 ! don't divide ice_sublim by cloud fraction since it is grid-averaged
2795! dum1 = one - dum1 - dum2
2796 dum1 = one - dum1 - dum2 - dum3
2797!--ag
2798 ice_sublim(i,k) = dum*dum1*oneodt
2799 endif
2800 endif
2801
2802 enddo
2803
2804 ! Big "administration" loop enforces conservation, updates variables
2805 ! that accumulate over substeps, and sets output variables.
2806
2807 do i=1,mgncol
2808
2809 ! get tendencies due to microphysical conversion processes
2810 !==========================================================
2811 ! note: tendencies are multiplied by appropriate cloud/precip
2812 ! fraction to get grid-scale values
2813 ! note: vap_dep is already grid-average values
2814
2815 ! The net tendencies need to be added to rather than overwritten,
2816 ! because they may have a value already set for instantaneous
2817 ! melting/freezing.
2818
2819!++ag
2820! qvlat(i,k) = qvlat(i,k) - (pre(i,k)+prds(i,k))*precip_frac(i,k)-&
2821! vap_dep(i,k)-ice_sublim(i,k)-mnuccd(i,k)-mnudep(i,k)*lcldm(i,k)
2822
2823 qvlat(i,k) = qvlat(i,k)-(pre(i,k)+prds(i,k))*precip_frac(i,k) &
2824 -vap_dep(i,k)-ice_sublim(i,k)-mnuccd(i,k)-mnudep(i,k)*lcldm(i,k) &
2825 -prdg(i,k)*precip_frac(i,k)
2826
2827! tlat(i,k) = tlat(i,k) + ((pre(i,k)*precip_frac(i,k)) &
2828! *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+ &
2829! ((bergs(i,k)+psacws(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k))*lcldm(i,k)+(mnuccr(i,k)+ &
2830! pracs(i,k)+mnuccri(i,k))*precip_frac(i,k)+berg(i,k))*xlf)
2831
2832! if (lprnt .and. k >= 60 .and. k <=65) &
2833! if (lprnt .and. k >= 100 ) &
2834! if (lprnt .and. abs(k-81) <5) &
2835! if (lprnt .and. k >= 60 ) &
2836! write(0,*)' k=',k,' tlat=',tlat(i,k),' pre=',pre(i,k),' precip_frac=',precip_frac(i,k),&
2837! ' prds=',prds(i,k),' prdg=',prdg(i,k),' vap_dep=',vap_dep(i,k),' ice_sublim=',ice_sublim(i,k), &
2838! ' mnuccd=',mnuccd(i,k),' mnudep=',mnudep(i,k),' lcldm=',lcldm(i,k),' bergs=',bergs(i,k), &
2839! ' psacws=',psacws(i,k),' mnuccc=',mnuccc(i,k),' mnucct=',mnucct(i,k),' msacwi=',msacwi(i,k), &
2840! ' psacwg=',psacwg(i,k),' qmultg=',qmultg(i,k),' pgsacw=',pgsacw(i,k),' mnuccr=',mnuccr(i,k), &
2841! ' pracs=',pracs(i,k),' mnuccri=',mnuccri(i,k),' pracg=',pracg(i,k),' pgracs=',pgracs(i,k), &
2842! ' qmultrg=',qmultrg(i,k),' xlf=',xlf,' xxlv=',xxlv,' xxls=',xxls
2843
2844
2845 tlat(i,k) = tlat(i,k)+((pre(i,k)*precip_frac(i,k))*xxlv+ &
2846 ((prds(i,k)+prdg(i,k))*precip_frac(i,k)+vap_dep(i,k)+ice_sublim(i,k)+ &
2847 mnuccd(i,k)+mnudep(i,k)*lcldm(i,k))*xxls+ &
2848 ((bergs(i,k)+psacws(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k)+psacwg(i,k)+ &
2849 qmultg(i,k)+pgsacw(i,k))*lcldm(i,k)+ &
2850 (mnuccr(i,k)+pracs(i,k)+mnuccri(i,k)+pracg(i,k)+pgracs(i,k)+qmultrg(i,k))*precip_frac(i,k)+ &
2851 berg(i,k))*xlf)
2852
2853! if (lprnt .and. k >= 100 ) write(0,*)' k=',k,' tlat=',tlat(i,k)
2854! if (lprnt) write(0,*)' k=',k,' tlat=',tlat(i,k)
2855! if (lprnt .and. k >= 60) write(0,*)' k=',k,' tlat=',tlat(i,k)
2856
2857! qctend(i,k) = qctend(i,k) + (-pra(i,k)-prc(i,k)-mnuccc(i,k)-mnucct(i,k)-msacwi(i,k)- &
2858! psacws(i,k)-bergs(i,k))*lcldm(i,k)-berg(i,k)
2859
2860 qctend(i,k) = qctend(i,k) + &
2861 (-pra(i,k)-prc(i,k)-mnuccc(i,k)-mnucct(i,k)-msacwi(i,k) - &
2862 psacws(i,k)-bergs(i,k)-qmultg(i,k)-psacwg(i,k)-pgsacw(i,k))*lcldm(i,k)-berg(i,k)
2863
2864 if (do_cldice) then
2865! qitend(i,k) = qitend(i,k) + &
2866! (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k))*lcldm(i,k)+(-prci(i,k)- &
2867! prai(i,k))*icldm(i,k)+vap_dep(i,k)+berg(i,k)+ice_sublim(i,k)+ &
2868! mnuccd(i,k)+mnuccri(i,k)*precip_frac(i,k)
2869
2870 qitend(i,k) = qitend(i,k)+ &
2871 (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k)+qmultg(i,k)) * lcldm(i,k) &
2872 + (-prci(i,k)-prai(i,k)) * icldm(i,k) &
2873 + vap_dep(i,k)+berg(i,k)+ice_sublim(i,k)+mnuccd(i,k) &
2874 + (mnuccri(i,k)+qmultrg(i,k)) * precip_frac(i,k)
2875
2876 end if
2877
2878! qrtend(i,k) = qrtend(i,k) + (pra(i,k)+prc(i,k))*lcldm(i,k)+(pre(i,k)-pracs(i,k)- &
2879! mnuccr(i,k)-mnuccri(i,k))*precip_frac(i,k)
2880
2881 qrtend(i,k) = qrtend(i,k) + (pra(i,k)+prc(i,k))*lcldm(i,k)+(pre(i,k)-pracs(i,k)- &
2882 mnuccr(i,k)-mnuccri(i,k)-qmultrg(i,k)-pracg(i,k)-pgracs(i,k))*precip_frac(i,k)
2883
2884
2885! qstend(i,k) = qstend(i,k) + (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k) &
2886! + (prds(i,k)+pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)
2887
2888 if (do_hail.or.do_graupel) then
2889 qgtend(i,k) = qgtend(i,k) + (pracg(i,k)+pgracs(i,k)+prdg(i,k)+psacr(i,k)+mnuccr(i,k))*precip_frac(i,k) &
2890 + (psacwg(i,k)+pgsacw(i,k))*lcldm(i,k)
2891
2892 qstend(i,k) = qstend(i,k) + (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k) &
2893 + (prds(i,k)+pracs(i,k)-psacr(i,k))*precip_frac(i,k)
2894
2895 else
2896 !necessary since mnuccr moved to graupel
2897 qstend(i,k) = qstend(i,k) + (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k) &
2898 + (prds(i,k)+pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)
2899
2900 endif
2901!--ag
2902
2903
2904 cmeout(i,k) = vap_dep(i,k) + ice_sublim(i,k) + mnuccd(i,k)
2905
2906 ! add output for cmei (accumulate)
2907 cmeitot(i,k) = vap_dep(i,k) + ice_sublim(i,k) + mnuccd(i,k)
2908
2909 ! assign variables for trop_mozart, these are grid-average
2910 !-------------------------------------------------------------------
2911 ! evaporation/sublimation is stored here as positive term
2912
2913!++add evaporation/sublimation of graupel too? YES: After conservation checks.
2914
2915!++ag
2916!ADD GRAUPEL to evapsnow: prdg. (sign? same as prds: negative, so this is a positive number)
2917! evapsnow(i,k) = -prds(i,k) * precip_frac(i,k)
2918 evapsnow(i,k) = (-prds(i,k)-prdg(i,k)) * precip_frac(i,k)
2919!--ag
2920 nevapr(i,k) = -pre(i,k) * precip_frac(i,k)
2921 prer_evap(i,k) = -pre(i,k) * precip_frac(i,k)
2922
2923 ! change to make sure prain is positive: do not remove snow from
2924 ! prain used for wet deposition
2925
2926!++AG NEED TO MAKE CONSISTENT WITH BUDGETS
2927 prain(i,k) = (pra(i,k)+prc(i,k))*lcldm(i,k) &
2928 - (pracs(i,k)+mnuccr(i,k)+mnuccri(i,k))*precip_frac(i,k)
2929 if (do_hail .or. do_graupel) then
2930! Subtract PSACR here or not? Ask Hugh
2931 prodsnow(i,k) = (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+ &
2932 pracs(i,k)*precip_frac(i,k)
2933 else
2934 prodsnow(i,k) = (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+ &
2935 (pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)
2936 endif
2937
2938 ! following are used to calculate 1st order conversion rate of cloud water
2939 ! to rain and snow (1/s), for later use in aerosol wet removal routine
2940 ! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc
2941 ! used to calculate pra, prc, ... in this routine
2942 ! qcsinksum_rate1ord = { rate of direct transfer of cloud water to rain & snow }
2943 ! (no cloud ice or bergeron terms)
2944
2945!++AG NEED TO MAKE CONSITANT: PGSACW, PSACWG (check budgets)? More sink terms? Check. No. Just loss to precip.
2946!Ask Hugh
2947! qcsinksum_rate1ord(i,k) = (pra(i,k)+prc(i,k)+psacws(i,k))*lcldm(i,k)
2948 qcsinksum_rate1ord(i,k) = (pra(i,k)+prc(i,k)+psacws(i,k)+psacwg(i,k)+pgsacw(i,k))*lcldm(i,k)
2949!--ag
2950 ! Avoid zero/near-zero division.
2951 qcsinksum_rate1ord(i,k) = qcsinksum_rate1ord(i,k) / max(qc(i,k),1.0e-30_r8)
2952
2953
2954 ! microphysics output, note this is grid-averaged
2955 pratot(i,k) = pra(i,k) * lcldm(i,k)
2956 prctot(i,k) = prc(i,k) * lcldm(i,k)
2957 mnuccctot(i,k) = mnuccc(i,k) * lcldm(i,k)
2958 mnuccttot(i,k) = mnucct(i,k) * lcldm(i,k)
2959 msacwitot(i,k) = msacwi(i,k) * lcldm(i,k)
2960 psacwstot(i,k) = psacws(i,k) * lcldm(i,k)
2961 bergstot(i,k) = bergs(i,k) * lcldm(i,k)
2962 bergtot(i,k) = berg(i,k)
2963 prcitot(i,k) = prci(i,k) * icldm(i,k)
2964 praitot(i,k) = prai(i,k) * icldm(i,k)
2965 mnuccdtot(i,k) = mnuccd(i,k) * icldm(i,k)
2966
2967 pracstot(i,k) = pracs(i,k) * precip_frac(i,k)
2968 mnuccrtot(i,k) = mnuccr(i,k) * precip_frac(i,k)
2969!++ag
2970 mnuccritot(i,k) = mnuccri(i,k) * precip_frac(i,k)
2971!--ag
2972
2973!++ag Hail/Graupel tendencies for output
2974 psacrtot(i,k) = psacr(i,k) * precip_frac(i,k)
2975 pracgtot(i,k) = pracg(i,k) * precip_frac(i,k)
2976 psacwgtot(i,k) = psacwg(i,k) * lcldm(i,k)
2977 pgsacwtot(i,k) = pgsacw(i,k) * lcldm(i,k)
2978 pgracstot(i,k) = pgracs(i,k) * precip_frac(i,k)
2979 prdgtot(i,k) = prdg(i,k) * precip_frac(i,k)
2980 qmultgtot(i,k) = qmultg(i,k) * lcldm(i,k)
2981 qmultrgtot(i,k) = qmultrg(i,k) * precip_frac(i,k)
2982 npracgtot(i,k) = npracg(i,k) * precip_frac(i,k)
2983 nscngtot(i,k) = nscng(i,k) * lcldm(i,k)
2984 ngracstot(i,k) = ngracs(i,k) * precip_frac(i,k)
2985 nmultgtot(i,k) = nmultg(i,k) * lcldm(i,k)
2986 nmultrgtot(i,k) = nmultrg(i,k) * precip_frac(i,k)
2987 npsacwgtot(i,k) = npsacwg(i,k) * lcldm(i,k)
2988!--ag
2989
2990!++ag
2991! nctend(i,k) = nctend(i,k) + (-nnuccc(i,k)-nnucct(i,k)-npsacws(i,k)+nsubc(i,k) &
2992! - npra(i,k)-nprc1(i,k))*lcldm(i,k)
2993
2994 nctend(i,k) = nctend(i,k) + (-nnuccc(i,k)-nnucct(i,k)-npsacws(i,k)+nsubc(i,k) &
2995 -npra(i,k)-nprc1(i,k)-npsacwg(i,k))*lcldm(i,k)
2996
2997 if (do_cldice) then
2998 if (use_hetfrz_classnuc) then
2999 tmpfrz = nnuccc(i,k)
3000 else
3001 tmpfrz = zero
3002 end if
3003! nitend(i,k) = nitend(i,k) + nnuccd(i,k)+ &
3004! (nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k))*lcldm(i,k)+(nsubi(i,k)-nprci(i,k)- &
3005! nprai(i,k))*icldm(i,k)+nnuccri(i,k)*precip_frac(i,k)
3006
3007 nitend(i,k) = nitend(i,k) + nnuccd(i,k) &
3008 + (nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k)+nmultg(i,k))*lcldm(i,k) &
3009 + (nsubi(i,k)-nprci(i,k)-nprai(i,k))*icldm(i,k) &
3010 + (nnuccri(i,k)+nmultrg(i,k))*precip_frac(i,k)
3011 endif
3012
3013 if(do_graupel.or.do_hail) then
3014! nstend(i,k) = nstend(i,k) + (nsubs(i,k)+nsagg(i,k)+nnuccr(i,k))*precip_frac(i,k) &
3015! + nprci(i,k)*icldm(i,k)
3016
3017 nstend(i,k) = nstend(i,k) + (nsubs(i,k)+nsagg(i,k)-ngracs(i,k))*precip_frac(i,k) &
3018 + nprci(i,k)*icldm(i,k)-nscng(i,k)*lcldm(i,k)
3019
3020 ngtend(i,k) = ngtend(i,k) + nscng(i,k)*lcldm(i,k) &
3021 + (ngracs(i,k)+nnuccr(i,k))*precip_frac(i,k)
3022
3023 else
3024 !necessary since mnuccr moved to graupel
3025 nstend(i,k) = nstend(i,k) + (nsubs(i,k)+nsagg(i,k)+nnuccr(i,k))*precip_frac(i,k) &
3026 + nprci(i,k)*icldm(i,k)
3027
3028 endif
3029
3030! nrtend(i,k) = nrtend(i,k) + nprc(i,k)*lcldm(i,k)+(nsubr(i,k)-npracs(i,k)-nnuccr(i,k) &
3031! - nnuccri(i,k)+nragg(i,k))*precip_frac(i,k)
3032
3033 nrtend(i,k) = nrtend(i,k)+ nprc(i,k)*lcldm(i,k) &
3034 + (nsubr(i,k)-npracs(i,k)-nnuccr(i,k) &
3035 -nnuccri(i,k)+nragg(i,k)-npracg(i,k)-ngracs(i,k))*precip_frac(i,k)
3036!--ag
3037
3038 ! make sure that ni at advanced time step does not exceed
3039 ! maximum (existing N + source terms*dt), which is possible if mtime < deltat
3040 ! note that currently mtime = deltat
3041 !================================================================
3042
3043 if (do_cldice .and. nitend(i,k) > zero .and. ni(i,k)+nitend(i,k)*deltat > nimax(i,k)) then
3044 nitend(i,k) = max(zero, (nimax(i,k)-ni(i,k))*oneodt)
3045 endif
3046
3047 enddo
3048
3049 ! End of "administration" loop
3050
3051 end do micro_vert_loop ! end k loop
3052
3053! if (lprnt) write(0,*)' tlat3=',tlat(1,:)*deltat
3054 !-----------------------------------------------------
3055 ! convert rain/snow q and N for output to history, note,
3056 ! output is for gridbox average
3057
3058 do k=1,nlev
3059 do i=1,mgncol
3060 qrout(i,k) = qr(i,k)
3061 nrout(i,k) = nr(i,k) * rho(i,k)
3062 qsout(i,k) = qs(i,k)
3063 nsout(i,k) = ns(i,k) * rho(i,k)
3064!++ag
3065 qgout(i,k) = qg(i,k)
3066 ngout(i,k) = ng(i,k) * rho(i,k)
3067!--ag
3068 enddo
3069 enddo
3070
3071 ! calculate n0r and lamr from rain mass and number
3072 ! divide by precip fraction to get in-precip (local) values of
3073 ! rain mass and number, divide by rhow to get rain number in kg^-1
3074
3075 do k=1,nlev
3076
3077 call size_dist_param_basic(mg_rain_props, qric(:,k), nric(:,k), lamr(:,k), mgncol, n0=n0r(:,k))
3078
3079 enddo
3080 ! Calculate rercld
3081
3082 ! calculate mean size of combined rain and cloud water
3083
3084 call calc_rercld(lamr, n0r, lamc, pgam, qric, qcic, ncic, rercld, mgncol, nlev)
3085
3086
3087 ! Assign variables back to start-of-timestep values
3088 ! Some state variables are changed before the main microphysics loop
3089 ! to make "instantaneous" adjustments. Afterward, we must move those changes
3090 ! back into the tendencies.
3091 ! These processes:
3092 ! - Droplet activation (npccn, impacts nc)
3093 ! - Instantaneous snow melting (minstsm/ninstsm, impacts qr/qs/nr/ns)
3094 ! - Instantaneous rain freezing (minstfr/ninstrf, impacts qr/qs/nr/ns)
3095 !================================================================================
3096
3097 do k=1,nlev
3098 do i=1,mgncol
3099 ! Re-apply droplet activation tendency
3100 nc(i,k) = ncn(i,k)
3101 nctend(i,k) = nctend(i,k) + npccn(i,k)
3102
3103 ! Re-apply rain freezing and snow melting.
3104 qstend(i,k) = qstend(i,k) + (qs(i,k)-qsn(i,k)) * oneodt
3105 qs(i,k) = qsn(i,k)
3106
3107 nstend(i,k) = nstend(i,k) + (ns(i,k)-nsn(i,k)) * oneodt
3108 ns(i,k) = nsn(i,k)
3109
3110 qrtend(i,k) = qrtend(i,k) + (qr(i,k)-qrn(i,k)) * oneodt
3111 qr(i,k) = qrn(i,k)
3112
3113 nrtend(i,k) = nrtend(i,k) + (nr(i,k)-nrn(i,k)) * oneodt
3114 nr(i,k) = nrn(i,k)
3115
3116!++ag Re-apply graupel freezing/melting
3117 qgtend(i,k) = qgtend(i,k) + (qg(i,k)-qgr(i,k)) * oneodt
3118 qg(i,k) = qgr(i,k)
3119
3120! if (maxval(dum_2D-qg).gt.0._r8) &
3121! write(iulog,*) "OOPS, qg diff > 0 : max=",maxval(dum_2D-qg)
3122! if (minval(dum_2D-qg).lt.0._r8) &
3123! write(iulog,*) "OOPS, qg diff < 0 : min=",minval(dum_2D-qg)
3124!
3125! write(iulog,*) "Max qgtend: 1st = ",maxval(qgtend)
3126! write(iulog,*) "Min qgtend: 1st = ",minval(qgtend)
3127! write(iulog,*) "Max qvtend: 1st = ",maxval(qvlat)
3128! write(iulog,*) "Min qvtend: 1st = ",minval(qvlat)
3129
3130 ngtend(i,k) = ngtend(i,k) + (ng(i,k)-ngr(i,k)) * oneodt
3131 ng(i,k) = ngr(i,k)
3132!--ag
3133
3134 !.............................................................................
3135
3136 !================================================================================
3137
3138 ! modify to include snow. in prain & evap (diagnostic here: for wet dep)
3139 nevapr(i,k) = nevapr(i,k) + evapsnow(i,k)
3140 prain(i,k) = prain(i,k) + prodsnow(i,k)
3141
3142
3143 enddo
3144 enddo
3145
3146 do k=1,nlev
3147
3148 do i=1,mgncol
3149
3150 ! calculate sedimentation for cloud water and ice
3151!++ag ! and Graupel (mg3)
3152 !================================================================================
3153
3154 ! update in-cloud cloud mixing ratio and number concentration
3155 ! with microphysical tendencies to calculate sedimentation, assign to dummy vars
3156 ! note: these are in-cloud values***, hence we divide by cloud fraction
3157
3158 if (lcldm(i,k) > mincld) then
3159 tx1 = one / lcldm(i,k)
3160 dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat) * tx1
3161 dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)*tx1, zero)
3162 ! switch for specification of droplet and crystal number
3163 if (nccons) then
3164 dumnc(i,k) = ncnst*rhoinv(i,k)
3165 end if
3166 else
3167 dumc(i,k) = zero
3168 dumnc(i,k) = zero
3169 endif
3170 if (icldm(i,k) > mincld) then
3171 tx1 = one / icldm(i,k)
3172 dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat) * tx1
3173 dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)*tx1, zero)
3174 ! switch for specification of cloud ice number
3175 if (nicons) then
3176 dumni(i,k) = ninst*rhoinv(i,k)
3177 end if
3178 else
3179 dumi(i,k) = zero
3180 dumni(i,k) = zero
3181 endif
3182 if (precip_frac(i,k) > mincld) then
3183 tx1 = one / precip_frac(i,k)
3184 dumr(i,k) = (qr(i,k)+qrtend(i,k)*deltat) * tx1
3185 dums(i,k) = (qs(i,k)+qstend(i,k)*deltat) * tx1
3186
3187 dumnr(i,k) = max((nr(i,k)+nrtend(i,k)*deltat)*tx1, zero)
3188 dumns(i,k) = max((ns(i,k)+nstend(i,k)*deltat)*tx1, zero)
3189
3190!++ag Add graupel
3191 dumg(i,k) = (qg(i,k)+qgtend(i,k)*deltat) * tx1
3192! Moorthi testing
3193 if (dumg(i,k) > 0.01_r8) then
3194 tx2 = dumg(i,k) - 0.01_r8
3195 dumg(i,k) = 0.01_r8
3196 dums(i,k) = dums(i,k) + tx2
3197 qstend(i,k) = (dums(i,k)*precip_frac(i,k) - qs(i,k)) * oneodt
3198 qgtend(i,k) = (dumg(i,k)*precip_frac(i,k) - qg(i,k)) * oneodt
3199 endif
3200! Moorthi testing
3201
3202 dumng(i,k) = max((ng(i,k)+ngtend(i,k)*deltat)*tx1, zero)
3203 ! switch for specification of droplet and crystal number
3204 if (ngcons) then
3205 dumng(i,k) = ngnst*rhoinv(i,k)
3206 endif
3207!--ag
3208 else
3209 dumr(i,k) = zero
3210 dumr(i,k) = zero
3211 dums(i,k) = zero
3212 dumns(i,k) = zero
3213!++ag Add graupel
3214 dumg(i,k) = zero
3215 dumng(i,k) = zero
3216 endif
3217!--ag
3218 enddo
3219 enddo
3220
3221 do k=1,nlev
3222
3223! obtain new slope parameter to avoid possible singularity
3224
3225 call size_dist_param_ice(mg_ice_props, dumi(:,k), dumni(:,k), &
3226 lami(:,k), mgncol)
3227
3228 call size_dist_param_liq(mg_liq_props, dumc(:,k), dumnc(:,k), rho(:,k), &
3229 pgam(:,k), lamc(:,k), mgncol)
3230
3231! call size_dist_param_basic(mg_ice_props, dumi(:,k), dumni(:,k), &
3232! lami(:,k), mgncol)
3233! fallspeed for rain
3234
3235 call size_dist_param_basic(mg_rain_props, dumr(:,k), dumnr(:,k), &
3236 lamr(:,k), mgncol)
3237! fallspeed for snow
3238 call size_dist_param_basic(mg_snow_props, dums(:,k), dumns(:,k), &
3239 lams(:,k), mgncol)
3240! fallspeed for graupel/hail
3241 if (do_graupel .or. do_hail) then
3242 call size_dist_param_basic(mg_graupel_props, dumg(:,k), dumng(:,k), &
3243 lamg(:,k), mgncol)
3244 endif
3245 enddo
3246
3247 do k=1,nlev
3248 do i=1,mgncol
3249
3250 ! calculate number and mass weighted fall velocity for droplets and cloud ice
3251 !-------------------------------------------------------------------
3252
3253 grho = g*rho(i,k)
3254
3255 if (dumc(i,k) >= qsmall) then
3256
3257 tx1 = lamc(i,k)**bc
3258 vtrmc(i,k) = acn(i,k)*gamma(pgam(i,k)+four+bc) &
3259 / (tx1*gamma(pgam(i,k)+four))
3260
3261 fc(i,k) = grho * vtrmc(i,k)
3262 fnc(i,k) = grho * acn(i,k)*gamma(pgam(i,k)+one+bc) &
3263 / (tx1*gamma(pgam(i,k)+one))
3264 else
3265 fc(i,k) = zero
3266 fnc(i,k) = zero
3267 end if
3268
3269 ! calculate number and mass weighted fall velocity for cloud ice
3270
3271 if (dumi(i,k) >= qsmall) then
3272
3273 tx3 = one / lami(i,k)
3274 tx1 = ain(i,k) * tx3**bi
3275 tx2 = 1.2_r8*rhof(i,k)
3276 vtrmi(i,k) = min(tx1*gamma_bi_plus4*oneo6, tx2)
3277
3278 fi(i,k) = grho * vtrmi(i,k)
3279 fni(i,k) = grho * min(tx1*gamma_bi_plus1, tx2)
3280
3281 ! adjust the ice fall velocity for smaller (r < 20 um) ice
3282 ! particles (blend over 18-20 um)
3283 irad = (1.5_r8 * 1e6_r8) * tx3
3284 ifrac = min(one, max(zero, (irad-18._r8)*half))
3285
3286 if (ifrac < one) then
3287 tx1 = ajn(i,k) / lami(i,k)**bj
3288 vtrmi(i,k) = ifrac*vtrmi(i,k) + (one-ifrac) * min(tx1*gamma_bj_plus4*oneo6, tx2)
3289
3290 fi(i,k) = grho * vtrmi(i,k)
3291 fni(i,k) = ifrac * fni(i,k) + (one-ifrac) * grho * min(tx1*gamma_bj_plus1, tx2)
3292 end if
3293 else
3294 fi(i,k) = zero
3295 fni(i,k)= zero
3296 endif
3297
3298 ! fallspeed for rain
3299
3300! if (lamr(i,k) >= qsmall) then
3301 if (dumr(i,k) >= qsmall) then
3302
3303 ! 'final' values of number and mass weighted mean fallspeed for rain (m/s)
3304
3305 tx1 = arn(i,k) / lamr(i,k)**br
3306 tx2 = 9.1_r8*rhof(i,k)
3307 umr(i,k) = min(tx1*gamma_br_plus4*oneo6, tx2)
3308 unr(i,k) = min(tx1*gamma_br_plus1, tx2)
3309
3310 fr(i,k) = grho * umr(i,k)
3311 fnr(i,k) = grho * unr(i,k)
3312
3313 else
3314 fr(i,k) = zero
3315 fnr(i,k) = zero
3316 endif
3317
3318 ! fallspeed for snow
3319
3320! if (lams(i,k) >= qsmall) then
3321 if (dums(i,k) >= qsmall) then
3322
3323 ! 'final' values of number and mass weighted mean fallspeed for snow (m/s)
3324 tx1 = asn(i,k) / lams(i,k)**bs
3325 tx2 = 1.2_r8*rhof(i,k)
3326 ums(i,k) = min(tx1*gamma_bs_plus4*oneo6, tx2)
3327 uns(i,k) = min(tx1*gamma_bs_plus1, tx2)
3328
3329 fs(i,k) = grho * ums(i,k)
3330 fns(i,k) = grho * uns(i,k)
3331
3332 else
3333 fs(i,k) = zero
3334 fns(i,k) = zero
3335 endif
3336
3337 if (do_graupel .or. do_hail) then
3338!++ag
3339 ! fallspeed for graupel
3340
3341
3342! if (lamg(i,k) >= qsmall) then
3343 if (dumg(i,k) >= qsmall) then
3344
3345 ! 'final' values of number and mass weighted mean fallspeed for graupel (m/s)
3346 tx1 = agn(i,k) / lamg(i,k)**bgtmp
3347 tx2 = 20._r8 * rhof(i,k)
3348 umg(i,k) = min(tx1*gamma_bg_plus4*oneo6, tx2)
3349 ung(i,k) = min(tx1*gamma_bg_plus1, tx2)
3350
3351 fg(i,k) = grho * umg(i,k)
3352 fng(i,k) = grho * ung(i,k)
3353
3354 else
3355 fg(i,k) = zero
3356 fng(i,k) = zero
3357 endif
3358 endif
3359
3360 ! redefine dummy variables - sedimentation is calculated over grid-scale
3361 ! quantities to ensure conservation
3362
3363 dumc(i,k) = qc(i,k) + qctend(i,k)*deltat
3364 dumi(i,k) = qi(i,k) + qitend(i,k)*deltat
3365 dumr(i,k) = qr(i,k) + qrtend(i,k)*deltat
3366 dums(i,k) = qs(i,k) + qstend(i,k)*deltat
3367
3368 dumnc(i,k) = nc(i,k) + nctend(i,k)*deltat
3369 dumni(i,k) = ni(i,k) + nitend(i,k)*deltat
3370 dumnr(i,k) = nr(i,k) + nrtend(i,k)*deltat
3371 dumns(i,k) = ns(i,k) + nstend(i,k)*deltat
3372!++ag
3373 dumg(i,k) = qg(i,k) + qgtend(i,k)*deltat
3374 dumng(i,k) = ng(i,k) + ngtend(i,k)*deltat
3375!--ag
3376
3377 if (dumc(i,k) < qsmall) dumnc(i,k) = zero
3378 if (dumi(i,k) < qsmall) dumni(i,k) = zero
3379 if (dumr(i,k) < qsmall) dumnr(i,k) = zero
3380 if (dums(i,k) < qsmall) dumns(i,k) = zero
3381 if (dumg(i,k) < qsmall) dumng(i,k) = zero
3382
3383 enddo
3384 enddo !!! vertical loop
3385
3386 do k=1,nlev
3387 do i=1,mgncol
3388 pdel_inv(i,k) = one / pdel(i,k)
3389 enddo
3390 enddo
3391! if (lprnt) write(0,*)' bef sedimentation dumc=',dumc(i,nlev-10:nlev)
3392
3393 ! initialize nstep for sedimentation sub-steps
3394
3395 ! calculate number of split time steps to ensure courant stability criteria
3396 ! for sedimentation calculations
3397 !-------------------------------------------------------------------
3398 do i=1,mgncol
3399 nlb = nlball(i)
3400 nstep = 1 + nint(max( maxval( fi(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
3401 maxval(fni(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)
3402 nstep = min(nstep, nstep_def)
3403
3404
3405 ! loop over sedimentation sub-time step to ensure stability
3406 !==============================================================
3407 if (do_cldice) then
3408 tx2 = one / nstep
3409 tx1 = tx2 * deltat
3410 tx3 = tx2 / g
3411
3412 do n = 1,nstep
3413
3414 ! top of model
3415
3416 k = 1
3417
3418 ! add fallout terms to microphysical tendencies
3419
3420 tx5 = dumi(i,k)
3421 tx7 = pdel_inv(i,k) * tx1
3422 dumi(i,k) = tx5 / (one + fi(i,k)*tx7)
3423 tx6 = (dumi(i,k)-tx5) * oneodt
3424 qitend(i,k) = qitend(i,k) + tx6
3425 tx5 = dumni(i,k)
3426 dumni(i,k) = tx5 / (one + fni(i,k)*tx7)
3427 nitend(i,k) = nitend(i,k) + (dumni(i,k)-tx5) * oneodt
3428
3429 ! sedimentation tendency for output
3430 qisedten(i,k) = qisedten(i,k) + tx6
3431
3432 falouti(k) = fi(i,k) * dumi(i,k)
3433 faloutni(k) = fni(i,k) * dumni(i,k)
3434
3435 iflx(i,k+1) = iflx(i,k+1) + falouti(k) * tx3 ! Ice flux
3436
3437 do k = 2,nlev
3438
3439 ! for cloud liquid and ice, if cloud fraction increases with height
3440 ! then add flux from above to both vapor and cloud water of current level
3441 ! this means that flux entering clear portion of cell from above evaporates
3442 ! instantly
3443
3444 ! note: this is not an issue with precip, since we assume max overlap
3445
3446 if (icldm(i,k-1) > mincld) then
3447 dum1 = max(zero, min(one, icldm(i,k)/icldm(i,k-1)))
3448 else
3449 dum1 = one
3450 endif
3451
3452 tx5 = dumi(i,k)
3453 tx7 = pdel_inv(i,k) * tx1
3454 dum2 = tx7 * dum1
3455 dumi(i,k) = (tx5 + falouti(k-1)*dum2) / (one + fi(i,k)*tx7)
3456 tx6 = (dumi(i,k)-tx5) * oneodt
3457 ! add fallout terms to eulerian tendencies
3458 qitend(i,k) = qitend(i,k) + tx6
3459 tx5 = dumni(i,k)
3460 dumni(i,k) = (tx5 + faloutni(k-1)*dum2) / (one + fni(i,k)*tx7)
3461 nitend(i,k) = nitend(i,k) + (dumni(i,k)-tx5) * oneodt
3462
3463
3464 qisedten(i,k) = qisedten(i,k) + tx6 ! sedimentation tendency for output
3465
3466
3467 falouti(k) = fi(i,k) * dumi(i,k)
3468 faloutni(k) = fni(i,k) * dumni(i,k)
3469
3470 dum2 = (one-dum1) * falouti(k-1) * pdel_inv(i,k) * tx2
3471 qvlat(i,k) = qvlat(i,k) + dum2 ! add terms to evap/sub of cloud ice
3472 qisevap(i,k) = qisevap(i,k) + dum2 ! for output
3473
3474 tlat(i,k) = tlat(i,k) - dum2 * xxls
3475
3476 iflx(i,k+1) = iflx(i,k+1) + falouti(k) * tx3 ! Ice flux
3477 end do
3478
3479 ! units below are m/s
3480 ! sedimentation flux at surface is added to precip flux at surface
3481 ! to get total precip (cloud + precip water) rate
3482
3483 prect(i) = prect(i) + falouti(nlev) * (tx3*0.001_r8)
3484 preci(i) = preci(i) + falouti(nlev) * (tx3*0.001_r8)
3485
3486 enddo
3487 endif
3488
3489! if (lprnt) write(0,*)' tlat4=',tlat(1,:)*deltat
3490 ! calculate number of split time steps to ensure courant stability criteria
3491 ! for sedimentation calculations
3492 !-------------------------------------------------------------------
3493 nstep = 1 + nint(max( maxval( fc(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
3494 maxval(fnc(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)
3495 nstep = min(nstep, nstep_def)
3496
3497 ! loop over sedimentation sub-time step to ensure stability
3498 !==============================================================
3499 tx2 = one / nstep
3500 tx1 = tx2 * deltat
3501 tx3 = tx2 / g
3502
3503 do n = 1,nstep
3504
3505 ! top of model
3506 k = 1
3507
3508 tx5 = dumc(i,k)
3509 tx7 = pdel_inv(i,k) * tx1
3510 dumc(i,k) = tx5 / (one + fc(i,k)*tx7)
3511 tx6 = (dumc(i,k)-tx5) * oneodt
3512 qctend(i,k) = qctend(i,k) + tx6
3513 tx5 = dumnc(i,k)
3514 dumnc(i,k) = tx5 / (one + fnc(i,k)*tx7)
3515 nctend(i,k) = nctend(i,k) + (dumnc(i,k)-tx5) * oneodt
3516
3517
3518 ! sedimentation tendency for output
3519 qcsedten(i,k) = qcsedten(i,k) + tx6
3520
3521 faloutc(k) = fc(i,k) * dumc(i,k)
3522 faloutnc(k) = fnc(i,k) * dumnc(i,k)
3523
3524 lflx(i,k+1) = lflx(i,k+1) + faloutc(k) * tx3
3525 do k = 2,nlev
3526
3527 if (lcldm(i,k-1) > mincld) then
3528 dum1 = max(zero, min(one, lcldm(i,k)/lcldm(i,k-1)))
3529 else
3530 dum1 = one
3531 endif
3532
3533 tx5 = dumc(i,k)
3534 tx7 = pdel_inv(i,k) * tx1
3535 dum2 = tx7 * dum1
3536 dumc(i,k) = (tx5 + faloutc(k-1)*dum2) / (one + fc(i,k)*tx7)
3537 tx6 = (dumc(i,k)-tx5) * oneodt
3538 qctend(i,k) = qctend(i,k) + tx6
3539 tx5 = dumnc(i,k)
3540 dumnc(i,k) = (tx5 + faloutnc(k-1)*dum2) / (one + fnc(i,k)*tx7)
3541 nctend(i,k) = nctend(i,k) + (dumnc(i,k)-tx5) * oneodt
3542
3543
3544
3545 qcsedten(i,k) = qcsedten(i,k) + tx6 ! sedimentation tendency for output
3546
3547 faloutc(k) = fc(i,k) * dumc(i,k)
3548 faloutnc(k) = fnc(i,k) * dumnc(i,k)
3549
3550 dum2 = (one-dum1) * faloutc(k-1) * pdel_inv(i,k) * tx2
3551 qvlat(i,k) = qvlat(i,k) + dum2 ! add terms to to evap/sub of cloud water
3552 qcsevap(i,k) = qcsevap(i,k) + dum2 ! for output
3553
3554 tlat(i,k) = tlat(i,k) - dum2 * xxlv
3555
3556 lflx(i,k+1) = lflx(i,k+1) + faloutc(k) * tx3 ! Liquid condensate flux here
3557 enddo
3558
3559 prect(i) = prect(i) + faloutc(nlev) * (tx3*0.001_r8)
3560
3561 end do
3562! if (lprnt) write(0,*)' tlat5=',tlat(1,:)*deltat
3563! if (lprnt) write(0,*)' maxval=',maxval( fr(i,nlb:nlev)*pdel_inv(i,nlb:nlev))&
3564! ,maxval(fnr(i,nlb:nlev)*pdel_inv(i,nlb:nlev))
3565
3566 ! calculate number of split time steps to ensure courant stability criteria
3567 ! for sedimentation calculations
3568 !-------------------------------------------------------------------
3569 nstep = 1 + nint(max( maxval( fr(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
3570 maxval(fnr(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)
3571
3572 nstep = min(nstep, nstep_def)
3573
3574 ! loop over sedimentation sub-time step to ensure stability
3575 !==============================================================
3576 tx2 = one / nstep
3577 tx1 = tx2 * deltat
3578 tx3 = tx2 / g
3579
3580! if(lprnt) then
3581! write(0,*)' nstep=',nstep,' tx1=',tx1,' tx2=',tx2,' tx3=',tx3,' qsmall=',qsmall
3582! write(0,*)' fr=',fr(i,:)
3583! write(0,*)' dumr=',dumr(i,:)
3584! endif
3585
3586 do n = 1,nstep
3587
3588 ! top of model
3589 k = 1
3590
3591 ! add fallout terms to microphysical tendencies
3592
3593 tx5 = dumr(i,k)
3594 tx7 = pdel_inv(i,k) * tx1
3595 dumr(i,k) = tx5 / (one + fr(i,k)*tx7)
3596 tx6 = (dumr(i,k)-tx5) * oneodt
3597 qrtend(i,k) = qrtend(i,k) + tx6
3598 tx5 = dumnr(i,k)
3599 dumnr(i,k) = tx5 / (one + fnr(i,k)*tx7)
3600 nrtend(i,k) = nrtend(i,k) + (dumnr(i,k)-tx5) * oneodt
3601
3602 ! sedimentation tendency for output
3603 qrsedten(i,k) = qrsedten(i,k) + tx6
3604
3605 faloutr(k) = fr(i,k) * dumr(i,k)
3606 faloutnr(k) = fnr(i,k) * dumnr(i,k)
3607
3608 rflx(i,k+1) = rflx(i,k+1) + faloutr(k) * tx3
3609
3610 do k = 2,nlev
3611
3612 tx5 = dumr(i,k)
3613 tx7 = pdel_inv(i,k) * tx1
3614 dumr(i,k) = (tx5 + faloutr(k-1)*tx7) / (one + fr(i,k)*tx7)
3615 tx6 = (dumr(i,k)-tx5) * oneodt
3616 qrtend(i,k) = qrtend(i,k) + tx6
3617 tx5 = dumnr(i,k)
3618 dumnr(i,k) = (tx5 + faloutnr(k-1)*tx7) / (one + fnr(i,k)*tx7)
3619 nrtend(i,k) = nrtend(i,k) + (dumnr(i,k)-tx5) * oneodt
3620
3621 qrsedten(i,k) = qrsedten(i,k) + tx6 ! sedimentation tendency for output
3622
3623 faloutr(k) = fr(i,k) * dumr(i,k)
3624 faloutnr(k) = fnr(i,k) * dumnr(i,k)
3625
3626 rflx(i,k+1) = rflx(i,k+1) + faloutr(k) * tx3 ! Rain Flux
3627 enddo
3628
3629 prect(i) = prect(i) + faloutr(nlev) * (tx3*0.001_r8)
3630
3631 enddo
3632
3633! if (lprnt) write(0,*)' prectaftrain=',prect(i),' preci=',preci(i)
3634
3635 ! calculate number of split time steps to ensure courant stability criteria
3636 ! for sedimentation calculations
3637 !-------------------------------------------------------------------
3638 nstep = 1 + nint(max( maxval( fs(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
3639 maxval(fns(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)
3640 nstep = min(nstep, nstep_def)
3641
3642
3643 ! loop over sedimentation sub-time step to ensure stability
3644 !==============================================================
3645 tx2 = one / nstep
3646 tx1 = tx2 * deltat
3647 tx3 = tx2 / g
3648 do n = 1,nstep
3649
3650 ! top of model
3651 k = 1
3652
3653 ! add fallout terms to microphysical tendencies
3654
3655 tx5 = dums(i,k)
3656 tx7 = pdel_inv(i,k) * tx1
3657 dums(i,k) = tx5 / (one + fs(i,k)*tx7)
3658 tx6 = (dums(i,k)-tx5) * oneodt
3659 qstend(i,k) = qstend(i,k) + tx6
3660 tx5 = dumns(i,k)
3661 dumns(i,k) = tx5 / (one + fns(i,k)*tx7)
3662 nstend(i,k) = nstend(i,k) + (dumns(i,k)-tx5) * oneodt
3663
3664 ! sedimentation tendency for output
3665 qssedten(i,k) = qssedten(i,k) + tx6
3666
3667 falouts(k) = fs(i,k) * dums(i,k)
3668 faloutns(k) = fns(i,k) * dumns(i,k)
3669
3670 sflx(i,k+1) = sflx(i,k+1) + falouts(k) * tx3
3671
3672 do k = 2,nlev
3673
3674
3675 tx5 = dums(i,k)
3676 tx7 = pdel_inv(i,k) * tx1
3677 dums(i,k) = (tx5 + falouts(k-1)*tx7) / (one + fs(i,k)*tx7)
3678 tx6 = (dums(i,k)-tx5) * oneodt
3679 qstend(i,k) = qstend(i,k) + tx6
3680 tx5 = dumns(i,k)
3681 dumns(i,k) = (tx5 + faloutns(k-1)*tx7) / (one + fns(i,k)*tx7)
3682 nstend(i,k) = nstend(i,k) + (dumns(i,k)-tx5) * oneodt
3683
3684
3685 qssedten(i,k) = qssedten(i,k) + tx6 ! sedimentation tendency for output
3686
3687 falouts(k) = fs(i,k) * dums(i,k)
3688 faloutns(k) = fns(i,k) * dumns(i,k)
3689
3690 sflx(i,k+1) = sflx(i,k+1) + falouts(k) * tx3 ! Snow Flux
3691 end do !! k loop
3692
3693 prect(i) = prect(i) + falouts(nlev) * (tx3*0.001_r8)
3694 preci(i) = preci(i) + falouts(nlev) * (tx3*0.001_r8)
3695
3696 enddo !! nstep loop
3697
3698! if (lprnt) write(0,*)' prectaftssno=',prect(i),' preci=',preci(i)
3699! if (lprnt) write(0,*)' qgtnd1=',qgtend(1,:)
3700
3701 if (do_graupel .or. do_hail) then
3702!++ag Graupel Sedimentation
3703 ! calculate number of split time steps to ensure courant stability criteria
3704 ! for sedimentation calculations
3705 !-------------------------------------------------------------------
3706 nstep = 1 + nint(max( maxval( fg(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
3707 maxval(fng(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)
3708 nstep = min(nstep, nstep_def)
3709
3710 tx2 = one / nstep
3711 tx1 = tx2 * deltat
3712 tx3 = tx2 / g
3713 ! loop over sedimentation sub-time step to ensure stability
3714 !==============================================================
3715 do n = 1,nstep
3716
3717 ! top of model
3718 k = 1
3719
3720 ! add fallout terms to microphysical tendencies
3721
3722 tx5 = dumg(i,k)
3723 tx7 = pdel_inv(i,k) * tx1
3724 dumg(i,k) = tx5 / (one + fg(i,k)*tx7)
3725 tx6 = (dumg(i,k)-tx5) * oneodt
3726 qgtend(i,k) = qgtend(i,k) + tx6
3727 tx5 = dumng(i,k)
3728 dumng(i,k) = tx5 / (one + fng(i,k)*tx7)
3729 ngtend(i,k) = ngtend(i,k) + (dumng(i,k)-tx5) * oneodt
3730
3731 ! sedimentation tendency for output
3732 qgsedten(i,k) = qgsedten(i,k) + tx6
3733
3734 faloutg(k) = fg(i,k) * dumg(i,k)
3735 faloutng(k) = fng(i,k) * dumng(i,k)
3736
3737 gflx(i,k+1) = gflx(i,k+1) + faloutg(k) * tx3 ! Ice flux
3738
3739 do k = 2,nlev
3740
3741 tx5 = dumg(i,k)
3742 tx7 = pdel_inv(i,k) * tx1
3743 dumg(i,k) = (tx5 + faloutg(k-1)*tx7) / (one + fg(i,k)*tx7)
3744 tx6 = (dumg(i,k)-tx5) * oneodt
3745 ! add fallout terms to eulerian tendencies
3746 qgtend(i,k) = qgtend(i,k) + tx6
3747 tx5 = dumng(i,k)
3748 dumng(i,k) = (tx5 + faloutng(k-1)*tx7) / (one + fng(i,k)*tx7)
3749 ngtend(i,k) = ngtend(i,k) + (dumng(i,k)-tx5) * oneodt
3750
3751
3752 qgsedten(i,k) = qgsedten(i,k) + tx6 ! sedimentation tendency for output
3753
3754
3755 faloutg(k) = fg(i,k) * dumg(i,k)
3756 faloutng(k) = fng(i,k) * dumng(i,k)
3757
3758 gflx(i,k+1) = gflx(i,k+1) + faloutg(k) * tx3 ! Ice flux
3759 enddo
3760
3761 ! units below are m/s
3762 ! sedimentation flux at surface is added to precip flux at surface
3763 ! to get total precip (cloud + precip water) rate
3764
3765 prect(i) = prect(i) + faloutg(nlev) * (tx3*0.001_r8)
3766 preci(i) = preci(i) + faloutg(nlev) * (tx3*0.001_r8)
3767
3768 enddo !! nstep loop
3769 endif
3770! if (lprnt) write(0,*)' qgtnds=',qgtend(1,:)
3771!--ag
3772 enddo ! end of i loop
3773 ! end sedimentation
3774
3775! if (lprnt) write(0,*)' prectaftsed=',prect(i),' preci=',preci(i)
3776
3777 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3778
3779 ! get new update for variables that includes sedimentation tendency
3780 ! note : here dum variables are grid-average, NOT in-cloud
3781
3782 do k=1,nlev
3783 do i=1,mgncol
3784 dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat, zero)
3785 dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat, zero)
3786 dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat, zero)
3787 dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat, zero)
3788
3789 dumr(i,k) = max(qr(i,k)+qrtend(i,k)*deltat, zero)
3790 dumnr(i,k) = max(nr(i,k)+nrtend(i,k)*deltat, zero)
3791 dums(i,k) = max(qs(i,k)+qstend(i,k)*deltat, zero)
3792 dumns(i,k) = max(ns(i,k)+nstend(i,k)*deltat, zero)
3793
3794!++ag
3795 dumg(i,k) = max(qg(i,k)+qgtend(i,k)*deltat, zero)
3796! Moorthi testing
3797 if (dumg(i,k) > 0.01_r8) then
3798 tx2 = dumg(i,k) - 0.01_r8
3799 dumg(i,k) = 0.01_r8
3800 dums(i,k) = dums(i,k) + tx2
3801 qstend(i,k) = (dums(i,k) - qs(i,k)) * oneodt
3802 qgtend(i,k) = (dumg(i,k) - qg(i,k)) * oneodt
3803 endif
3804! Moorthi testing
3805 dumng(i,k) = max(ng(i,k)+ngtend(i,k)*deltat, zero)
3806!--ag
3807
3808 ! switch for specification of droplet and crystal number
3809 if (nccons) then
3810 dumnc(i,k) = ncnst*rhoinv(i,k)*lcldm(i,k)
3811 endif
3812
3813 ! switch for specification of cloud ice number
3814 if (nicons) then
3815 dumni(i,k) = ninst*rhoinv(i,k)*icldm(i,k)
3816 endif
3817
3818!++ag
3819 ! switch for specification of graupel number
3820 if (ngcons) then
3821 dumng(i,k) = ngnst*rhoinv(i,k)*precip_frac(i,k)
3822 endif
3823!--ag
3824
3825 if (dumc(i,k) < qsmall) dumnc(i,k) = zero
3826 if (dumi(i,k) < qsmall) dumni(i,k) = zero
3827 if (dumr(i,k) < qsmall) dumnr(i,k) = zero
3828 if (dums(i,k) < qsmall) dumns(i,k) = zero
3829!++ag
3830 if (dumg(i,k) < qsmall) dumng(i,k) = zero
3831!--ag
3832
3833 enddo
3834
3835 enddo
3836
3837 ! calculate instantaneous processes (melting, homogeneous freezing)
3838 !====================================================================
3839
3840 ! melting of snow at +2 C
3841 do k=1,nlev
3842
3843 do i=1,mgncol
3844
3845 tx1 = t(i,k) + tlat(i,k)*(deltat/cpp) - snowmelt
3846 if (tx1 > zero) then
3847 if (dums(i,k) > zero) then
3848
3849 ! make sure melting snow doesn't reduce temperature below threshold
3850 dum = -(xlf/cpp) * dums(i,k)
3851 if (tx1+dum < zero) then
3852 dum = min(one, max(zero, -tx1/dum))
3853 else
3854 dum = one
3855 end if
3856
3857 tx2 = dum * oneodt
3858 qstend(i,k) = qstend(i,k) - tx2*dums(i,k)
3859 nstend(i,k) = nstend(i,k) - tx2*dumns(i,k)
3860 qrtend(i,k) = qrtend(i,k) + tx2*dums(i,k)
3861 nrtend(i,k) = nrtend(i,k) + tx2*dumns(i,k)
3862
3863 dum1 = - xlf * tx2 * dums(i,k)
3864 tlat(i,k) = dum1 + tlat(i,k)
3865 meltsdttot(i,k) = dum1 + meltsdttot(i,k)
3866 endif
3867 endif
3868 enddo
3869 enddo
3870
3871 if (do_graupel .or. do_hail) then
3872!++ag
3873
3874 ! melting of graupel at +2 C
3875
3876 do k=1,nlev
3877
3878 do i=1,mgncol
3879
3880 tx1 = t(i,k) + tlat(i,k)*(deltat/cpp) - snowmelt
3881 if (tx1 > zero) then
3882 if (dumg(i,k) > zero) then
3883
3884 ! make sure melting graupel doesn't reduce temperature below threshold
3885 dum = -(xlf/cpp) * dumg(i,k)
3886 if (tx1+dum < zero) then
3887 dum = max(zero, min(one, -tx1/dum))
3888 else
3889 dum = one
3890 end if
3891
3892 tx2 = dum * oneodt
3893
3894 qgtend(i,k) = qgtend(i,k) - tx2*dumg(i,k)
3895 ngtend(i,k) = ngtend(i,k) - tx2*dumng(i,k)
3896 qrtend(i,k) = qrtend(i,k) + tx2*dumg(i,k)
3897 nrtend(i,k) = nrtend(i,k) + tx2*dumng(i,k)
3898
3899 dum1 = - xlf*tx2*dumg(i,k)
3900 tlat(i,k) = dum1 + tlat(i,k)
3901 meltsdttot(i,k) = dum1 + meltsdttot(i,k)
3902 endif
3903 endif
3904 enddo
3905 enddo
3906
3907!--ag
3908 endif
3909
3910 do k=1,nlev
3911 do i=1,mgncol
3912
3913 ! freezing of rain at -5 C
3914
3915 tx1 = t(i,k) + tlat(i,k) * (deltat/cpp) - rainfrze
3916 if (tx1 < zero) then
3917
3918 if (dumr(i,k) > zero) then
3919
3920 ! make sure freezing rain doesn't increase temperature above threshold
3921 dum = (xlf/cpp) * dumr(i,k)
3922 if (tx1+dum > zero) then
3923 dum = min(one, max(zero, -tx1/dum))
3924 else
3925 dum = one
3926 end if
3927 tx2 = dum * oneodt
3928 qrtend(i,k) = qrtend(i,k) - tx2 * dumr(i,k)
3929 nrtend(i,k) = nrtend(i,k) - tx2 * dumnr(i,k)
3930
3931 ! get mean size of rain = 1/lamr, add frozen rain to either snow or cloud ice
3932 ! depending on mean rain size
3933
3934 call size_dist_param_basic(mg_rain_props, dumr(i,k), dumnr(i,k), lamr(i,k))
3935
3936 if (lamr(i,k) < one/dcs) then
3937!++ag freeze rain to graupel
3938 if (do_hail .or. do_graupel) then
3939 qgtend(i,k) = qgtend(i,k) + tx2 * dumr(i,k)
3940 ngtend(i,k) = ngtend(i,k) + tx2 * dumnr(i,k)
3941 else
3942 qstend(i,k) = qstend(i,k) + tx2 * dumr(i,k)
3943 nstend(i,k) = nstend(i,k) + tx2 * dumnr(i,k)
3944 end if
3945!--ag
3946 else
3947 qitend(i,k) = qitend(i,k) + tx2 * dumr(i,k)
3948 nitend(i,k) = nitend(i,k) + tx2 * dumnr(i,k)
3949 end if
3950 ! heating tendency
3951 dum1 = xlf*dum*dumr(i,k)*oneodt
3952 frzrdttot(i,k) = dum1 + frzrdttot(i,k)
3953 tlat(i,k) = dum1 + tlat(i,k)
3954
3955 endif
3956 endif
3957
3958 enddo
3959 enddo
3960 if (do_cldice) then
3961 do k=1,nlev
3962 do i=1,mgncol
3963 tx1 = t(i,k) + tlat(i,k) * (deltat/cpp) - tmelt
3964 if (tx1 > zero) then
3965 if (dumi(i,k) > zero) then
3966
3967 ! limit so that melting does not push temperature below freezing
3968 !-----------------------------------------------------------------
3969 dum = -dumi(i,k)*xlf/cpp
3970 if (tx1+dum < zero) then
3971 dum = min(one, max(zero, -tx1/dum))
3972 else
3973 dum = one
3974 end if
3975
3976 tx2 = dum * oneodt
3977 qctend(i,k) = qctend(i,k) + tx2*dumi(i,k)
3978
3979 ! for output
3980 melttot(i,k) = tx2*dumi(i,k)
3981
3982 ! assume melting ice produces droplet
3983 ! mean volume radius of 8 micron
3984
3985 nctend(i,k) = nctend(i,k) + three*tx2*dumi(i,k)/(four*pi*5.12e-16_r8*rhow)
3986
3987 qitend(i,k) = ((one-dum)*dumi(i,k)-qi(i,k)) * oneodt
3988 nitend(i,k) = ((one-dum)*dumni(i,k)-ni(i,k)) * oneodt
3989 tlat(i,k) = tlat(i,k) - xlf*tx2*dumi(i,k)
3990 endif
3991 endif
3992 enddo
3993 enddo
3994
3995! if (lprnt) write(0,*)' tlat6=',tlat(1,:)*deltat
3996! if (lprnt) write(0,*)' qitend=',qitend(1,nlev-45:nlev)*deltat
3997! if (lprnt) write(0,*)' qctend=',qctend(1,nlev-45:nlev)*deltat
3998
3999 ! homogeneously freeze droplets at -40 C
4000 !-----------------------------------------------------------------
4001
4002 do k=1,nlev
4003 do i=1,mgncol
4004 tx1 = t(i,k) + tlat(i,k)*(deltat/cpp) - 233.15_r8
4005 if (tx1 < zero) then
4006 if (dumc(i,k) > zero) then
4007
4008 ! limit so that freezing does not push temperature above threshold
4009 dum = (xlf/cpp) * dumc(i,k)
4010 if (tx1+dum > zero) then
4011 dum = min(one, max(zero, -tx1/dum))
4012 else
4013 dum = one
4014 end if
4015
4016 tx2 = dum * oneodt * dumc(i,k)
4017 qitend(i,k) = tx2 + qitend(i,k)
4018 homotot(i,k) = tx2 ! for output
4019
4020 ! assume 25 micron mean volume radius of homogeneously frozen droplets
4021 ! consistent with size of detrained ice in stratiform.F90
4022
4023 nitend(i,k) = nitend(i,k) + tx2*(three/(four*pi*1.563e-14_r8* 500._r8))
4024 qctend(i,k) = ((one-dum)*dumc(i,k)-qc(i,k)) * oneodt
4025 nctend(i,k) = ((one-dum)*dumnc(i,k)-nc(i,k)) * oneodt
4026 tlat(i,k) = tlat(i,k) + xlf*tx2
4027 endif
4028 endif
4029 enddo
4030 enddo
4031 ! remove any excess over-saturation, which is possible due to non-linearity when adding
4032 ! together all microphysical processes
4033 !-----------------------------------------------------------------
4034 ! follow code similar to old CAM scheme
4035 do k=1,nlev
4036 do i=1,mgncol
4037
4038 qtmp = q(i,k) + qvlat(i,k) * deltat
4039 ttmp = t(i,k) + tlat(i,k) * (deltat/cpp)
4040
4041 ! use rhw to allow ice supersaturation
4042 !call qsat_water(ttmp, p(i,k), esn, qvn)
4043 esn = min(fpvsl(ttmp), p(i,k))
4044 qvn = epsqs*esn/(p(i,k)-omeps*esn) * qsfm(i,k)
4045! qvn = epsqs*esn/(p(i,k)-omeps*esn)
4046
4047
4048 if (qtmp > qvn .and. qvn > zero .and. allow_sed_supersat) then
4049 ! expression below is approximate since there may be ice deposition
4050 dum = (qtmp-qvn)/(one+xxlv_squared*qvn/(cpp*rv*ttmp*ttmp)) * oneodt
4051 ! add to output cme
4052 cmeout(i,k) = cmeout(i,k) + dum
4053 ! now add to tendencies, partition between liquid and ice based on temperature
4054 if (ttmp > 268.15_r8) then
4055 dum1 = zero
4056 ! now add to tendencies, partition between liquid and ice based on te
4057 !-------------------------------------------------------
4058 else if (ttmp < 238.15_r8) then
4059 dum1 = one
4060 else
4061 dum1 = (268.15_r8-ttmp)/30._r8
4062 end if
4063
4064 tx1 = xxls*dum1 + xxlv*(one-dum1)
4065 dum = (qtmp-qvn)/(one+tx1*tx1*qvn/(cpp*rv*ttmp*ttmp)) * oneodt
4066 tx2 = dum*(one-dum1)
4067 qctend(i,k) = qctend(i,k) + tx2
4068 qcrestot(i,k) = tx2 ! for output
4069 qitend(i,k) = qitend(i,k) + dum*dum1
4070 qirestot(i,k) = dum*dum1
4071 qvlat(i,k) = qvlat(i,k) - dum
4072 ! for output
4073 qvres(i,k) = -dum
4074 tlat(i,k) = tlat(i,k) + dum*tx1
4075 endif
4076 enddo
4077 enddo
4078 endif
4079
4080! if (lprnt) write(0,*)' tlat7=',tlat(1,:)*deltat
4081
4082 ! calculate effective radius for pass to radiation code
4083 !=========================================================
4084 ! if no cloud water, default value is 10 micron for droplets,
4085 ! 25 micron for cloud ice
4086
4087 ! update cloud variables after instantaneous processes to get effective radius
4088 ! variables are in-cloud to calculate size dist parameters
4089 do k=1,nlev
4090 do i=1,mgncol
4091 if (lcldm(i,k) > mincld) then
4092 tx1 = one / lcldm(i,k)
4093 else
4094 tx1 = zero
4095 endif
4096 if (icldm(i,k) > mincld) then
4097 tx2 = one / icldm(i,k)
4098 else
4099 tx2 = zero
4100 endif
4101 if (precip_frac(i,k) > mincld) then
4102 tx3 = one / precip_frac(i,k)
4103 else
4104 tx3 = zero
4105 endif
4106 dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat, zero) * tx1
4107 dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat, zero) * tx2
4108 dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat, zero) * tx1
4109 dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat, zero) * tx2
4110
4111 dumr(i,k) = max(qr(i,k)+qrtend(i,k)*deltat, zero) * tx3
4112 dumnr(i,k) = max(nr(i,k)+nrtend(i,k)*deltat, zero) * tx3
4113 dums(i,k) = max(qs(i,k)+qstend(i,k)*deltat, zero) * tx3
4114 dumns(i,k) = max(ns(i,k)+nstend(i,k)*deltat, zero) * tx3
4115
4116!++ag
4117 dumg(i,k) = max(qg(i,k)+qgtend(i,k)*deltat, zero) * tx3
4118 dumng(i,k) = max(ng(i,k)+ngtend(i,k)*deltat, zero) * tx3
4119!--ag
4120
4121 ! switch for specification of droplet and crystal number
4122 if (nccons) then
4123 dumnc(i,k) = ncnst * rhoinv(i,k)
4124 end if
4125
4126 ! switch for specification of cloud ice number
4127 if (nicons) then
4128 dumni(i,k) = ninst * rhoinv(i,k)
4129 end if
4130
4131!++ag
4132 ! switch for specification of graupel number
4133 if (ngcons) then
4134 dumng(i,k) = ngnst*rhoinv(i,k)*precip_frac(i,k)
4135 end if
4136!--ag
4137
4138 ! limit in-cloud mixing ratio to reasonable value of 5 g kg-1
4139! dumc(i,k) = min(dumc(i,k), 5.e-3_r8)
4140! dumi(i,k) = min(dumi(i,k), 5.e-3_r8)
4141 dumc(i,k) = min(dumc(i,k), 10.e-3_r8)
4142 dumi(i,k) = min(dumi(i,k), 10.e-3_r8)
4143 ! limit in-precip mixing ratios
4144 dumr(i,k) = min(dumr(i,k), 10.e-3_r8)
4145 dums(i,k) = min(dums(i,k), 10.e-3_r8)
4146!++ag
4147 dumg(i,k) = min(dumg(i,k), 10.e-3_r8)
4148!--ag
4149 enddo
4150 enddo
4151 ! cloud ice effective radius
4152 !-----------------------------------------------------------------
4153
4154 if (do_cldice) then
4155 do k=1,nlev
4156 do i=1,mgncol
4157 if (dumi(i,k) >= qsmall) then
4158
4159 tx1 = dumni(i,k)
4160 call size_dist_param_basic(mg_ice_props, dumi(i,k), dumni(i,k), &
4161 lami(i,k), dumni0)
4162
4163 if (dumni(i,k) /= tx1) then
4164 ! adjust number conc if needed to keep mean size in reasonable range
4165 nitend(i,k) = (dumni(i,k)*icldm(i,k)-ni(i,k)) * oneodt
4166 end if
4167
4168 tx1 = one / lami(i,k)
4169! effi(i,k) = (1.5_r8*1.e6_r8) * tx1
4170 effi(i,k) = (three*1.e6_r8) * tx1
4171 sadice(i,k) = two*pi*(tx1*tx1*tx1)*dumni0*rho(i,k)*1.e-2_r8 ! m2/m3 -> cm2/cm3
4172
4173 else
4174 effi(i,k) = 50._r8
4175 sadice(i,k) = zero
4176 end if
4177
4178 ! ice effective diameter for david mitchell's optics
4179 deffi(i,k) = effi(i,k) * (rhoi+rhoi)/rhows
4180 enddo
4181 enddo
4182 !else
4183 !do k=1,nlev
4184 !do i=1,mgncol
4185 ! NOTE: If CARMA is doing the ice microphysics, then the ice effective
4186 ! radius has already been determined from the size distribution.
4187 !effi(i,k) = re_ice(i,k) * 1.e6_r8 ! m -> um
4188 !deffi(i,k)=effi(i,k) * 2._r8
4189 !sadice(i,k) = 4._r8*pi*(effi(i,k)**2)*ni(i,k)*rho(i,k)*1e-2_r8
4190 !enddo
4191 !enddo
4192 end if
4193
4194 ! cloud droplet effective radius
4195 !-----------------------------------------------------------------
4196 do k=1,nlev
4197 do i=1,mgncol
4198 if (dumc(i,k) >= qsmall) then
4199
4200
4201 ! switch for specification of droplet and crystal number
4202 if (nccons) then
4203 ! make sure nc is consistence with the constant N by adjusting tendency, need
4204 ! to multiply by cloud fraction
4205 ! note that nctend may be further adjusted below if mean droplet size is
4206 ! out of bounds
4207
4208 nctend(i,k) = (ncnst*rhoinv(i,k)*lcldm(i,k)-nc(i,k)) * oneodt
4209
4210 end if
4211
4212 dum = dumnc(i,k)
4213
4214 call size_dist_param_liq(mg_liq_props, dumc(i,k), dumnc(i,k), rho(i,k), &
4215 pgam(i,k), lamc(i,k))
4216
4217 if (dum /= dumnc(i,k)) then
4218 ! adjust number conc if needed to keep mean size in reasonable range
4219 nctend(i,k) = (dumnc(i,k)*lcldm(i,k)-nc(i,k)) * oneodt
4220 end if
4221
4222 effc(i,k) = (half*1.e6_r8) * (pgam(i,k)+three) / lamc(i,k)
4223 !assign output fields for shape here
4224 lamcrad(i,k) = lamc(i,k)
4225 pgamrad(i,k) = pgam(i,k)
4226
4227
4228 ! recalculate effective radius for constant number, in order to separate
4229 ! first and second indirect effects
4230 !======================================
4231 ! assume constant number of 10^8 kg-1
4232
4233 dumnc(i,k) = 1.e8_r8
4234
4235 ! Pass in "false" adjust flag to prevent number from being changed within
4236 ! size distribution subroutine.
4237 call size_dist_param_liq(mg_liq_props, dumc(i,k), dumnc(i,k), rho(i,k), &
4238 pgam(i,k), lamc(i,k))
4239
4240 effc_fn(i,k) = (half*1.e6_r8) * (pgam(i,k)+three)/lamc(i,k)
4241
4242 else
4243 effc(i,k) = ten
4244 lamcrad(i,k) = zero
4245 pgamrad(i,k) = zero
4246 effc_fn(i,k) = ten
4247 endif
4248 enddo
4249 enddo
4250 ! recalculate 'final' rain size distribution parameters
4251 ! to ensure that rain size is in bounds, adjust rain number if needed
4252 do k=1,nlev
4253 do i=1,mgncol
4254
4255 if (dumr(i,k) >= qsmall) then
4256
4257 dum = dumnr(i,k)
4258
4259 call size_dist_param_basic(mg_rain_props, dumr(i,k), dumnr(i,k), lamr(i,k))
4260
4261 if (dum /= dumnr(i,k)) then
4262 ! adjust number conc if needed to keep mean size in reasonable range
4263 nrtend(i,k) = (dumnr(i,k)*precip_frac(i,k)-nr(i,k)) *oneodt
4264 endif
4265
4266 endif
4267 enddo
4268 enddo
4269 ! recalculate 'final' snow size distribution parameters
4270 ! to ensure that snow size is in bounds, adjust snow number if needed
4271 do k=1,nlev
4272 do i=1,mgncol
4273 if (dums(i,k) >= qsmall) then
4274
4275 dum = dumns(i,k)
4276
4277 call size_dist_param_basic(mg_snow_props, dums(i,k), dumns(i,k), &
4278 lams(i,k), n0=dumns0)
4279
4280 if (dum /= dumns(i,k)) then
4281 ! adjust number conc if needed to keep mean size in reasonable range
4282 nstend(i,k) = (dumns(i,k)*precip_frac(i,k)-ns(i,k)) * oneodt
4283 end if
4284
4285 tx1 = (two*pi*1.e-2_r8) / (lams(i,k)*lams(i,k)*lams(i,k))
4286 sadsnow(i,k) = tx1*dumns0*rho(i,k) ! m2/m3 -> cm2/cm3
4287
4288 endif
4289
4290
4291 enddo ! vertical k loop
4292 enddo
4293 do k=1,nlev
4294 do i=1,mgncol
4295 ! if updated q (after microphysics) is zero, then ensure updated n is also zero
4296 !=================================================================================
4297 if (qc(i,k)+qctend(i,k)*deltat < qsmall) nctend(i,k) = -nc(i,k) * oneodt
4298 if (do_cldice .and. qi(i,k)+qitend(i,k)*deltat < qsmall) nitend(i,k) = -ni(i,k) * oneodt
4299 if (qr(i,k)+qrtend(i,k)*deltat < qsmall) nrtend(i,k) = -nr(i,k) * oneodt
4300 if (qs(i,k)+qstend(i,k)*deltat < qsmall) nstend(i,k) = -ns(i,k) * oneodt
4301!++ag
4302 if (qg(i,k)+qgtend(i,k)*deltat < qsmall) ngtend(i,k) = -ng(i,k) * oneodt
4303!--ag
4304
4305 enddo
4306
4307 enddo
4308
4309 ! DO STUFF FOR OUTPUT:
4310 !==================================================
4311
4312 do k=1,nlev
4313 do i=1,mgncol
4314
4315 ! qc and qi are only used for output calculations past here,
4316 ! so add qctend and qitend back in one more time
4317 qc(i,k) = qc(i,k) + qctend(i,k)*deltat
4318 qi(i,k) = qi(i,k) + qitend(i,k)*deltat
4319
4320 ! averaging for snow and rain number and diameter
4321 !--------------------------------------------------
4322
4323 ! drout2/dsout2:
4324 ! diameter of rain and snow
4325 ! dsout:
4326 ! scaled diameter of snow (passed to radiation in CAM)
4327 ! reff_rain/reff_snow:
4328 ! calculate effective radius of rain and snow in microns for COSP using Eq. 9 of COSP v1.3 manual
4329
4330 if (qrout(i,k) > 1.e-7_r8 .and. nrout(i,k) > zero) then
4331 qrout2(i,k) = qrout(i,k) * precip_frac(i,k)
4332 nrout2(i,k) = nrout(i,k) * precip_frac(i,k)
4333 ! The avg_diameter call does the actual calculation; other diameter
4334 ! outputs are just drout2 times constants.
4335 drout2(i,k) = avg_diameter(qrout(i,k), nrout(i,k), rho(i,k), rhow)
4336 freqr(i,k) = precip_frac(i,k)
4337
4338 reff_rain(i,k) = (1.e6_r8*1.5_r8) * drout2(i,k)
4339 else
4340 qrout2(i,k) = zero
4341 nrout2(i,k) = zero
4342 drout2(i,k) = zero
4343 freqr(i,k) = zero
4344 reff_rain(i,k) = zero
4345 endif
4346
4347 if (qsout(i,k) > 1.e-7_r8 .and. nsout(i,k) > zero) then
4348 qsout2(i,k) = qsout(i,k) * precip_frac(i,k)
4349 nsout2(i,k) = nsout(i,k) * precip_frac(i,k)
4350 ! The avg_diameter call does the actual calculation; other diameter
4351 ! outputs are just dsout2 times constants.
4352 dsout2(i,k) = avg_diameter(qsout(i,k), nsout(i,k), rho(i,k), rhosn)
4353 freqs(i,k) = precip_frac(i,k)
4354
4355 dsout(i,k) = three*rhosn/rhows*dsout2(i,k)
4356
4357 reff_snow(i,k) = (1.e6_r8*three) * dsout2(i,k)
4358 else
4359 dsout(i,k) = zero
4360 qsout2(i,k) = zero
4361 nsout2(i,k) = zero
4362 dsout2(i,k) = zero
4363 freqs(i,k) = zero
4364 reff_snow(i,k) = zero
4365 endif
4366
4367 enddo
4368 enddo
4369
4370 ! analytic radar reflectivity
4371 !--------------------------------------------------
4372 ! formulas from Matthew Shupe, NOAA/CERES
4373 ! *****note: radar reflectivity is local (in-precip average)
4374 ! units of mm^6/m^3
4375
4376 do k=1,nlev
4377 do i = 1,mgncol
4378! if (qc(i,k) >= qsmall .and. (nc(i,k)+nctend(i,k)*deltat) > ten .and. lcldm(i,k) > mincld) then
4379 if (qc(i,k) >= qsmall .and. (nc(i,k)+nctend(i,k)*deltat) > ten) then
4380 tx1 = rho(i,k) / lcldm(i,k)
4381 tx2 = 1000._r8 * qc(i,k) * tx1
4382 dum = tx2 * tx2 * lcldm(i,k) &
4383 /(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)*tx1*1.e-6_r8*precip_frac(i,k))
4384! dum = (qc(i,k)/lcldm(i,k)*rho(i,k)*1000._r8)**2 &
4385! /(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)
4386 else
4387 dum = zero
4388 end if
4389! if (qi(i,k) >= qsmall .and. icldm(i,k) > mincld) then
4390 if (qi(i,k) >= qsmall) then
4391! 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)
4392 dum1 = (qi(i,k)*rho(i,k)/icldm(i,k)*10000._r8)**(one/0.63_r8)*icldm(i,k)/precip_frac(i,k)
4393 else
4394 dum1 = zero
4395 end if
4396
4397 if (qsout(i,k) >= qsmall) then
4398! dum1 = dum1 + (qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)**(one/0.63_r8)
4399 dum1 = dum1 + (qsout(i,k)*rho(i,k)*10000._r8)**(one/0.63_r8)
4400 end if
4401
4402 refl(i,k) = dum + dum1
4403
4404 ! add rain rate, but for 37 GHz formulation instead of 94 GHz
4405 ! formula approximated from data of Matrasov (2007)
4406 ! rainrt is the rain rate in mm/hr
4407 ! reflectivity (dum) is in DBz
4408
4409 if (rainrt(i,k) >= 0.001_r8) then
4410 dum = rainrt(i,k) * rainrt(i,k)
4411 dum = log10(dum*dum*dum) + 16._r8
4412
4413 ! convert from DBz to mm^6/m^3
4414
4415 dum = ten**(dum/ten)
4416 else
4417 ! don't include rain rate in R calculation for values less than 0.001 mm/hr
4418 dum = zero
4419 end if
4420
4421 ! add to refl
4422
4423 refl(i,k) = refl(i,k) + dum
4424
4425 !output reflectivity in Z.
4426 areflz(i,k) = refl(i,k) * precip_frac(i,k)
4427
4428 ! convert back to DBz
4429
4430 if (refl(i,k) > minrefl) then
4431 refl(i,k) = ten*log10(refl(i,k))
4432 else
4433 refl(i,k) = -9999._r8
4434 end if
4435
4436 !set averaging flag
4437 if (refl(i,k) > mindbz) then
4438 arefl(i,k) = refl(i,k) * precip_frac(i,k)
4439 frefl(i,k) = precip_frac(i,k)
4440 else
4441 arefl(i,k) = zero
4442 areflz(i,k) = zero
4443 frefl(i,k) = zero
4444 end if
4445
4446 ! bound cloudsat reflectivity
4447
4448 csrfl(i,k) = min(csmax,refl(i,k))
4449
4450 !set averaging flag
4451 if (csrfl(i,k) > csmin) then
4452 acsrfl(i,k) = refl(i,k) * precip_frac(i,k)
4453 fcsrfl(i,k) = precip_frac(i,k)
4454 else
4455 acsrfl(i,k) = zero
4456 fcsrfl(i,k) = zero
4457 endif
4458
4459 enddo
4460 enddo
4461
4462 do k=1,nlev
4463 do i = 1,mgncol
4464 !redefine fice here....
4465 tx2 = qsout(i,k) + qi(i,k)
4466 tx1 = tx2 + qrout(i,k) + qc(i,k)
4467 if ( tx2 > qsmall .and. tx1 > qsmall) then
4468 nfice(i,k) = min(tx2/tx1, one)
4469 else
4470 nfice(i,k) = zero
4471 endif
4472 enddo
4473 enddo
4474
4475end subroutine micro_mg_tend
4477
4478!========================================================================
4479!OUTPUT CALCULATIONS
4480!========================================================================
4481
4483subroutine calc_rercld(lamr, n0r, lamc, pgam, qric, qcic, ncic, rercld, mgncol,nlev)
4484 integer, intent(in) :: mgncol, nlev ! horizontal and vertical dimension
4485 real(r8), dimension(mgncol,nlev), intent(in) :: lamr ! rain size parameter (slope)
4486 real(r8), dimension(mgncol,nlev), intent(in) :: n0r ! rain size parameter (intercept)
4487 real(r8), dimension(mgncol,nlev), intent(in) :: lamc ! size distribution parameter (slope)
4488 real(r8), dimension(mgncol,nlev), intent(in) :: pgam ! droplet size parameter
4489 real(r8), dimension(mgncol,nlev), intent(in) :: qric ! in-cloud rain mass mixing ratio
4490 real(r8), dimension(mgncol,nlev), intent(in) :: qcic ! in-cloud cloud liquid
4491 real(r8), dimension(mgncol,nlev), intent(in) :: ncic ! in-cloud droplet number concentration
4492
4493 real(r8), dimension(mgncol,nlev), intent(inout) :: rercld ! effective radius calculation for rain + cloud
4494
4495 ! combined size of precip & cloud drops
4496 real(r8) :: Atmp
4497
4498 integer :: i, k
4499
4500 do k=1,nlev
4501 do i=1,mgncol
4502 ! Rain drops
4503 if (lamr(i,k) > zero) then
4504 atmp = n0r(i,k) * (half*pi) / (lamr(i,k)*lamr(i,k)*lamr(i,k))
4505 else
4506 atmp = zero
4507 end if
4508
4509 ! Add cloud drops
4510 if (lamc(i,k) > zero) then
4511 atmp = atmp + ncic(i,k) * pi * rising_factorial(pgam(i,k)+one, 2) &
4512 / (four*lamc(i,k)*lamc(i,k))
4513 end if
4514
4515 if (atmp > zero) then
4516 rercld(i,k) = rercld(i,k) + three *(qric(i,k) + qcic(i,k)) / (four * rhow * atmp)
4517 endif
4518 enddo
4519 enddo
4520end subroutine calc_rercld
4521
4522!========================================================================
4523
4524end module micro_mg3_0
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 graupel_rain_riming_snow(pracs, npracs, psacr, qsic, qric, nric, nsic, n0s, lams, n0r, lamr, dtime, pgracs, ngracs, mgncol)
Rain riming snow to graupel.
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 graupel_riming_liquid_snow(psacws, qsic, qcic, nsic, rho, rhosn, rhog, asn, lams, n0s, dtime, pgsacw, nscng, mgncol)
Conversion of rimed cloud water onto snow to graupel/hail.
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 graupel_rime_splintering(t, qcic, qric, qgic, psacwg, pracg, qmultg, nmultg, qmultrg, nmultrg, mgncol)
Rime splintering.
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 graupel_collecting_rain(qric, qgic, umg, umr, ung, unr, rho, n0r, lamr, n0g, lamg, pracg, npracg, mgncol)
CHANGE IN Q,N COLLECTION RAIN BY GRAUPEL.
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 evaporate_sublimate_precip_graupel(t, rho, dv, mu, sc, q, qvl, qvi, lcldm, precip_frac, arn, asn, agn, bg, qcic, qiic, qric, qsic, qgic, lamr, n0r, lams, n0s, lamg, n0g, pre, prds, prdg, am_evp_st, mgncol)
evaporation/sublimation of rain, snow and graupel
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 graupel_collecting_snow(qsic, qric, umr, ums, rho, lamr, n0r, lams, n0s, psacr, mgncol)
Collection of snow by rain to form graupel.
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 graupel_collecting_cld_water(qgic, qcic, ncic, rho, n0g, lamg, bg, agn, psacwg, npsacwg, mgncol)
Collection of cloud water by graupel.
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
MG microphysics version 3.0 - Update of MG microphysics with prognostic hail OR graupel.
Definition sflx.f:3