CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
micro_mg_utils.F90
1
4
27
28!--------------------------------------------------------------------------
29!
30! List of required external functions that must be supplied:
31! gamma --> standard mathematical gamma function (if gamma is an
32! intrinsic, define HAVE_GAMMA_INTRINSICS)
33!
34!--------------------------------------------------------------------------
35!
36! Constants that must be specified in the "init" method (module variables):
37!
38! kind kind of reals (to verify correct linkage only) -
39! gravit acceleration due to gravity m s-2
40! rair dry air gas constant for air J kg-1 K-1
41! rh2o gas constant for water vapor J kg-1 K-1
42! cpair specific heat at constant pressure for dry air J kg-1 K-1
43! tmelt temperature of melting point for water K
44! latvap latent heat of vaporization J kg-1
45! latice latent heat of fusion J kg-1
46!
47!--------------------------------------------------------------------------
48
49! 8 byte real and integer
50use machine, only : r8 => kind_phys
51use machine, only : i8 => kind_phys
52implicit none
53private
54save
55
56public :: &
82!++ag
90! graupel_sublimate_evap
91!--ag
92
93
94public :: mghydrometeorprops
95
97 ! Density (kg/m^3)
98 real(r8) :: rho
99 ! Information for size calculations.
100 ! Basic calculation of mean size is:
101 ! lambda = (shape_coef*nic/qic)^(1/eff_dim)
102 ! Then lambda is constrained by bounds.
103 real(r8) :: eff_dim
104 real(r8) :: shape_coef
105 real(r8) :: lambda_bounds(2)
106 ! Minimum average particle mass (kg).
107 ! Limit is applied at the beginning of the size distribution calculations.
108 real(r8) :: min_mean_mass
109end type mghydrometeorprops
110
112 module procedure newmghydrometeorprops
113end interface
114
115type(mghydrometeorprops), public :: mg_liq_props
116type(mghydrometeorprops), public :: mg_ice_props
117type(mghydrometeorprops), public :: mg_rain_props
118type(mghydrometeorprops), public :: mg_snow_props
119!++ag
120type(mghydrometeorprops), public :: mg_graupel_props
121!--ag
122
124 module procedure size_dist_param_liq_vect
125 module procedure size_dist_param_liq_line
126end interface
128 module procedure size_dist_param_basic_vect
129 module procedure size_dist_param_basic_line
130end interface
131
133 module procedure size_dist_param_ice_vect
134 module procedure size_dist_param_ice_line
135end interface
136
137!=================================================
138! Public module parameters (mostly for MG itself)
139!=================================================
140
142real(r8), parameter, public :: pi = 3.14159265358979323846_r8
143
145!real(r8), parameter, public :: omsm = 1._r8 - 1.e-5_r8
146real(r8), parameter, public :: omsm = 1._r8 - 1.e-6_r8
147
149real(r8), parameter, public :: qsmall = 1.e-18_r8
150
152 real(r8), parameter, public :: mincld = 0.000001_r8
153!real(r8), parameter, public :: mincld = 0.0001_r8
154!real(r8), parameter, public :: mincld = 0.0_r8
155
156real(r8), parameter, public :: rhosn = 250._r8
157real(r8), parameter, public :: rhoi = 500._r8
158real(r8), parameter, public :: rhow = 1000._r8
159real(r8), parameter, public :: rhows = 917._r8
160
161!++ag
162!Hail and Graupel (set in MG3)
163real(r8), parameter, public :: rhog = 500._r8
164real(r8), parameter, public :: rhoh = 400._r8
165!--ag
166
167! fall speed parameters, V = aD^b (V is in m/s)
168! droplets
169real(r8), parameter, public :: ac = 3.e7_r8
170real(r8), parameter, public :: bc = 2._r8
171! snow
172real(r8), parameter, public :: as = 11.72_r8
173real(r8), parameter, public :: bs = 0.41_r8
174! cloud ice
175real(r8), parameter, public :: ai = 700._r8
176real(r8), parameter, public :: bi = 1._r8
177! small cloud ice (r< 10 um) - sphere, bulk density
178real(r8), parameter, public :: aj = ac*((rhoi/rhows)**(bc/3._r8))*rhows/rhow
179real(r8), parameter, public :: bj = bc
180! rain
181real(r8), parameter, public :: ar = 841.99667_r8
182real(r8), parameter, public :: br = 0.8_r8
183!++ag
184! graupel
185real(r8), parameter, public :: ag = 19.3_r8
186real(r8), parameter, public :: bg = 0.37_r8
187! hail
188real(r8), parameter, public :: ah = 114.5_r8
189real(r8), parameter, public :: bh = 0.5_r8
190!--ag
191
195real(r8), parameter, public :: mi0 = 4._r8/3._r8*pi*rhoi*(1.e-6_r8)**3
196
197!++ag
198! mass of new graupel particle (assume same as mi0 for now, may want to make bigger?)
199!real(r8), parameter, public :: mg0 = 4._r8/3._r8*pi*rhoi*(1.e-6_r8)**3
200!or set based on M2005:
201real(r8), parameter, public :: mg0 = 1.6e-10_r8
202! radius of contact nuclei
203real(r8), parameter, public :: mmult = 4._r8/3._r8*pi*rhoi*(5.e-6_r8)**3
204!--ag
205
206!=================================================
207! Private module parameters
208!=================================================
209
210! Signaling NaN bit pattern that represents a limiter that's turned off.
211integer(i8), parameter :: limiter_off = int(z'7FF1111111111111', i8)
212
213! alternate threshold used for some in-cloud mmr
214real(r8), parameter :: icsmall = 1.e-8_r8
215
216! particle mass-diameter relationship
217! currently we assume spherical particles for cloud ice/snow
218! m = cD^d
219! exponent
220real(r8), parameter :: dsph = 3._r8
221
222! Bounds for mean diameter for different constituents.
223real(r8), parameter :: lam_bnd_rain(2) = 1._r8/[500.e-6_r8, 20.e-6_r8]
224real(r8), parameter :: lam_bnd_snow(2) = 1._r8/[2000.e-6_r8, 10.e-6_r8]
225
226! Minimum average mass of particles.
227real(r8), parameter :: min_mean_mass_liq = 1.e-20_r8
228real(r8), parameter :: min_mean_mass_ice = 1.e-20_r8
229
230! ventilation parameters
231! for snow
232real(r8), parameter :: f1s = 0.86_r8
233real(r8), parameter :: f2s = 0.28_r8
234! for rain
235real(r8), parameter :: f1r = 0.78_r8
236real(r8), parameter :: f2r = 0.308_r8
237
238! collection efficiencies
239! aggregation of cloud ice and snow
240!real(r8), parameter :: eii = 0.5_r8
241!real(r8), parameter :: eii = 0.1_r8
242 real(r8), parameter :: eii = 0.2_r8
243!++ag
244! collection efficiency, ice-droplet collisions
245real(r8), parameter, public :: ecid = 0.7_r8
246! collection efficiency between droplets/rain and snow/rain
247real(r8), parameter, public :: ecr = 1.0_r8
248!--ag
249
250! immersion freezing parameters, bigg 1953
251real(r8), parameter :: bimm = 100._r8
252real(r8), parameter :: aimm = 0.66_r8
253
254! Mass of each raindrop created from autoconversion.
255real(r8), parameter :: droplet_mass_25um = 4._r8/3._r8*pi*rhow*(25.e-6_r8)**3
256real(r8), parameter :: droplet_mass_40um = 4._r8/3._r8*pi*rhow*(40.e-6_r8)**3, &
257 droplet_mass_40umi = 1._r8/droplet_mass_40um
258
259!=========================================================
260! Constants set in initialization
261!=========================================================
262
263! Set using arguments to micro_mg_init
264real(r8) :: rv ! water vapor gas constant
265real(r8) :: cpp ! specific heat of dry air
266real(r8) :: tmelt ! freezing point of water (K)
267
268real(r8) :: ra ! dry air gas constant
269
270! latent heats of:
271real(r8) :: xxlv ! vaporization
272real(r8) :: xlf ! freezing
273real(r8) :: xxls ! sublimation
274
275! additional constants to help speed up code
276real(r8) :: gamma_bs_plus3
277real(r8) :: gamma_half_br_plus5
278real(r8) :: gamma_half_bs_plus5
279!++ag
280real(r8) :: gamma_2bs_plus2
281!--ag
282!
283real(r8), parameter :: zero = 0._r8, one = 1._r8, two = 2._r8, three = 3._r8, &
284 four = 4._r8, five = 5._r8, six = 6._r8, pio6 = pi/six, &
285 pio3 = pi/three, half = 0.5_r8, oneo3 = one/three, &
286 twopi = pi + pi
287
288!=========================================================
289! Utilities that are cheaper if the compiler knows that
290! some argument is an integer.
291!=========================================================
292
295 module procedure rising_factorial_r8
296 module procedure rising_factorial_integer
297end interface rising_factorial
298
300interface var_coef
301 module procedure var_coef_r8
302 module procedure var_coef_integer
303end interface var_coef
304
305!==========================================================================
306contains
307!==========================================================================
308
311!
312! "kind" serves no purpose here except to check for unlikely linking
313! issues; always pass in the kind for a double precision real.
314!
315!
316! Check the list at the top of this module for descriptions of all other
317! arguments.
318subroutine micro_mg_utils_init( kind, rair, rh2o, cpair, tmelt_in, latvap, &
319 latice, dcs)
320! latice, dcs, errstring)
321
322 integer, intent(in) :: kind
323!++ag
324 real(r8), intent(in) :: rair
325!--ag
326 real(r8), intent(in) :: rh2o
327 real(r8), intent(in) :: cpair
328 real(r8), intent(in) :: tmelt_in
329 real(r8), intent(in) :: latvap
330 real(r8), intent(in) :: latice
331 real(r8), intent(in) :: dcs
332
333
334 ! Name this array to workaround an XLF bug (otherwise could just use the
335 ! expression that sets it).
336 real(r8) :: ice_lambda_bounds(2)
337
338 !-----------------------------------------------------------------------
339
340
341 ! declarations for MG code (transforms variable names)
342
343 rv = rh2o ! water vapor gas constant
344 cpp = cpair ! specific heat of dry air
345 tmelt = tmelt_in
346 !++ag
347 ra = rair ! dry air gas constant
348 !--ag
349
350 ! latent heats
351
352 xxlv = latvap ! latent heat vaporization
353 xlf = latice ! latent heat freezing
354 xxls = xxlv + xlf ! latent heat of sublimation
355
356 ! Define constants to help speed up code (this limits calls to gamma function)
357 gamma_bs_plus3 = gamma(three+bs)
358 gamma_half_br_plus5 = gamma((five+br)*half)
359 gamma_half_bs_plus5 = gamma((five+bs)*half)
360!++ag
361 gamma_2bs_plus2 = gamma(bs+bs+two)
362!--ag
363
364
365 ! Don't specify lambda bounds for cloud liquid, as they are determined by
366 ! pgam dynamically.
367 mg_liq_props = mghydrometeorprops(rhow, dsph, min_mean_mass=min_mean_mass_liq)
368
369 ! Mean ice diameter can not grow bigger than twice the autoconversion
370 ! threshold for snow.
371 ice_lambda_bounds = one/[two*dcs, 1.e-6_r8]
372
373 mg_ice_props = mghydrometeorprops(rhoi, dsph, &
374 ice_lambda_bounds, min_mean_mass_ice)
375
376 mg_rain_props = mghydrometeorprops(rhow, dsph, lam_bnd_rain)
377 mg_snow_props = mghydrometeorprops(rhosn, dsph, lam_bnd_snow)
378!++ag
379 mg_graupel_props = mghydrometeorprops(rhog, dsph, lam_bnd_snow)
380!--ag
381
382end subroutine micro_mg_utils_init
383
386function newmghydrometeorprops(rho, eff_dim, lambda_bounds, min_mean_mass) &
387 result(res)
388 real(r8), intent(in) :: rho, eff_dim
389 real(r8), intent(in), optional :: lambda_bounds(2), min_mean_mass
390 type(mghydrometeorprops) :: res
391
392 res%rho = rho
393 res%eff_dim = eff_dim
394 if (present(lambda_bounds)) then
395 res%lambda_bounds = lambda_bounds
396 else
397 res%lambda_bounds = no_limiter()
398 end if
399 if (present(min_mean_mass)) then
400 res%min_mean_mass = min_mean_mass
401 else
402 res%min_mean_mass = no_limiter()
403 end if
404
405 res%shape_coef = rho * pio6 * gamma(eff_dim+one)
406
407end function newmghydrometeorprops
408
409!========================================================================
410!FORMULAS
411!========================================================================
412
413! Use gamma function to implement rising factorial extended to the reals.
414pure function rising_factorial_r8(x, n) result(res)
415 real(r8), intent(in) :: x, n
416 real(r8) :: res
417
418 res = gamma(x+n) / gamma(x)
419
420end function rising_factorial_r8
421
422! Rising factorial can be performed much cheaper if n is a small integer.
423pure function rising_factorial_integer(x, n) result(res)
424 real(r8), intent(in) :: x
425 integer, intent(in) :: n
426 real(r8) :: res
427
428 integer :: i
429 real(r8) :: factor
430
431 res = one
432 factor = x
433
434 do i = 1, n
435 res = res * factor
436 factor = factor + one
437 end do
438
439end function rising_factorial_integer
440
441! Calculate correction due to latent heat for evaporation/sublimation
442elemental function calc_ab(t, qv, xxl) result(ab)
443 real(r8), intent(in) :: t ! Temperature
444 real(r8), intent(in) :: qv ! Saturation vapor pressure
445 real(r8), intent(in) :: xxl ! Latent heat
446
447 real(r8) :: ab
448
449 real(r8) :: dqsdt
450
451 dqsdt = xxl*qv / (rv*t*t)
452 ab = one + dqsdt*xxl/cpp
453
454end function calc_ab
455
458elemental subroutine size_dist_param_liq_line(props, qcic, ncic, rho, pgam, lamc)
459 type(mghydrometeorprops), intent(in) :: props
460 real(r8), intent(in) :: qcic
461 real(r8), intent(inout) :: ncic
462 real(r8), intent(in) :: rho
463
464 real(r8), intent(out) :: pgam
465 real(r8), intent(out) :: lamc
466 real(r8) :: xs
467
468 type(mghydrometeorprops) :: props_loc
469 logical, parameter :: liq_gmao=.true.
470
471 if (qcic > qsmall) then
472
473 ! Local copy of properties that can be modified.
474 ! (Elemental routines that operate on arrays can't modify scalar
475 ! arguments.)
476 props_loc = props
477
478 ! Get pgam from fit to observations of martin et al. 1994
479
480 if (liq_gmao) then
481 pgam = 0.0005714_r8*1.e-6_r8*ncic*rho + 0.2714_r8
482 ! Anning modified lamc
483 if ((ncic > 1.0e-3_r8) .and. (qcic > 1.0e-11_r8)) then
484 xs = 0.07_r8*(1000._r8*qcic/ncic) ** (-0.14_r8)
485 else
486 xs = 1.2_r8
487 end if
488
489 xs = max(min(xs, 1.7_r8), 1.1_r8)
490 xs = xs*xs*xs
491 xs = (xs + sqrt(xs+8.0_r8)*sqrt(xs) - 4.0_r8)/8.0_r8
492 pgam = sqrt(xs)
493 else
494
495 pgam = one - 0.7_r8 * exp(-0.008_r8*1.e-6_r8*ncic*rho)
496! pgam = 0.0005714_r8*1.e-6_r8*ncic*rho + 0.2714_r8
497 endif
498 pgam = one / (pgam*pgam) - one
499 pgam = max(pgam, two)
500
501 ! Set coefficient for use in size_dist_param_basic.
502 ! The 3D case is so common and optimizable that we specialize it:
503 if (props_loc%eff_dim == three) then
504 props_loc%shape_coef = pio6 * props_loc%rho * &
505 rising_factorial(pgam+one, 3)
506 else
507 props_loc%shape_coef = pio6 * props_loc%rho * &
508 rising_factorial(pgam+one, props_loc%eff_dim)
509 end if
510
511 ! Limit to between 2 and 50 microns mean size.
512 props_loc%lambda_bounds = (pgam+one) * one/[50.e-6_r8, 2.e-6_r8]
513
514 call size_dist_param_basic(props_loc, qcic, ncic, lamc)
515
516 else
517 ! pgam not calculated in this case, so set it to a value likely to
518 ! cause an error if it is accidentally used
519 ! (gamma function undefined for negative integers)
520 pgam = -100._r8
521 lamc = zero
522 end if
523
524end subroutine size_dist_param_liq_line
525
528subroutine size_dist_param_liq_vect(props, qcic, ncic, rho, pgam, lamc, mgncol)
529
530 type(mghydrometeorprops), intent(in) :: props
531 integer, intent(in) :: mgncol
532 real(r8), dimension(mgncol), intent(in) :: qcic
533 real(r8), dimension(mgncol), intent(inout) :: ncic
534 real(r8), dimension(mgncol), intent(in) :: rho
535 real(r8), dimension(mgncol), intent(out) :: pgam
536 real(r8), dimension(mgncol), intent(out) :: lamc
537 real(r8) :: xs
538 type(mghydrometeorprops) :: props_loc
539 logical, parameter :: liq_gmao=.true.
540 integer :: i
541
542 do i=1,mgncol
543 if (qcic(i) > qsmall) then
544 ! Local copy of properties that can be modified.
545 ! (Elemental routines that operate on arrays can't modify scalar
546 ! arguments.)
547 props_loc = props
548 ! Get pgam from fit to observations of martin et al. 1994
549
550 if (liq_gmao) then
551 pgam(i) = 0.0005714_r8*1.e-6_r8*ncic(i)*rho(i) + 0.2714_r8
552 if ((ncic(i) > 1.0e-3_r8) .and. (qcic(i) > 1.0e-11_r8)) then
553 xs = 0.07_r8*(1000._r8*qcic(i)/ncic(i)) **(-0.14_r8)
554 else
555 xs = 1.2_r8
556 end if
557
558 xs = max(min(xs, 1.7_r8), 1.1_r8)
559 xs = xs*xs*xs
560 xs = (xs + sqrt(xs+8.0_r8)*sqrt(xs) - 4.0_r8)/8.0_r8
561 pgam(i) = sqrt(xs)
562 else
563 pgam(i) = one - 0.7_r8 * exp(-0.008_r8*1.e-6_r8*ncic(i)*rho(i))
564! pgam(i) = 0.0005714_r8*1.e-6_r8*ncic(i)*rho(i) + 0.2714_r8
565 endif
566
567 pgam(i) = one/(pgam(i)*pgam(i)) - one
568 pgam(i) = max(pgam(i), two)
569 endif
570 enddo
571 do i=1,mgncol
572 if (qcic(i) > qsmall) then
573 ! Set coefficient for use in size_dist_param_basic.
574 ! The 3D case is so common and optimizable that we specialize
575 ! it:
576 if (props_loc%eff_dim == three) then
577 props_loc%shape_coef = pio6 * props_loc%rho * &
578 rising_factorial(pgam(i)+one, 3)
579 else
580 props_loc%shape_coef = pio6 * props_loc%rho * &
581 rising_factorial(pgam(i)+one, props_loc%eff_dim)
582 end if
583 ! Limit to between 2 and 50 microns mean size.
584 props_loc%lambda_bounds(1) = (pgam(i)+one) / 50.e-6_r8
585 props_loc%lambda_bounds(2) = (pgam(i)+one) / 2.e-6_r8
586 call size_dist_param_basic(props_loc, qcic(i), ncic(i), lamc(i))
587 endif
588 enddo
589 do i=1,mgncol
590 if (qcic(i) <= qsmall) then
591 ! pgam not calculated in this case, so set it to a value likely to
592 ! cause an error if it is accidentally used
593 ! (gamma function undefined for negative integers)
594 pgam(i) = -100._r8
595 lamc(i) = zero
596 end if
597 enddo
598
599end subroutine size_dist_param_liq_vect
600
603elemental subroutine size_dist_param_basic_line(props, qic, nic, lam, n0)
604 type(mghydrometeorprops), intent(in) :: props
605 real(r8), intent(in) :: qic
606 real(r8), intent(inout) :: nic
607
608 real(r8), intent(out) :: lam
609 real(r8), intent(out), optional :: n0
610
611 if (qic > qsmall) then
612
613 ! add upper limit to in-cloud number concentration to prevent
614 ! numerical error
615 if (limiter_is_on(props%min_mean_mass)) then
616 nic = min(nic, qic / props%min_mean_mass)
617 end if
618
619 ! lambda = (c n/q)^(1/d)
620 lam = (props%shape_coef * nic/qic)**(one/props%eff_dim)
621
622 ! check for slope
623 ! adjust vars
624 if (lam < props%lambda_bounds(1)) then
625 lam = props%lambda_bounds(1)
626 nic = lam**(props%eff_dim) * qic/props%shape_coef
627 else if (lam > props%lambda_bounds(2)) then
628 lam = props%lambda_bounds(2)
629 nic = lam**(props%eff_dim) * qic/props%shape_coef
630 end if
631
632 else
633 lam = zero
634 end if
635
636 if (present(n0)) n0 = nic * lam
637
638end subroutine size_dist_param_basic_line
639
642subroutine size_dist_param_basic_vect(props, qic, nic, lam, mgncol, n0)
643
644 type (mghydrometeorprops), intent(in) :: props
645 integer, intent(in) :: mgncol
646 real(r8), dimension(mgncol), intent(in) :: qic
647 real(r8), dimension(mgncol), intent(inout) :: nic
648 real(r8), dimension(mgncol), intent(out) :: lam
649 real(r8), dimension(mgncol), intent(out), optional :: n0
650 integer :: i
651 do i=1,mgncol
652
653 if (qic(i) > qsmall) then
654
655 ! add upper limit to in-cloud number concentration to prevent
656 ! numerical error
657 if (limiter_is_on(props%min_mean_mass)) then
658 nic(i) = min(nic(i), qic(i) / props%min_mean_mass)
659 end if
660
661 ! lambda = (c n/q)^(1/d)
662 lam(i) = (props%shape_coef * nic(i)/qic(i))**(one/props%eff_dim)
663
664 ! check for slope
665 ! adjust vars
666 if (lam(i) < props%lambda_bounds(1)) then
667 lam(i) = props%lambda_bounds(1)
668 nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef
669 else if (lam(i) > props%lambda_bounds(2)) then
670 lam(i) = props%lambda_bounds(2)
671 nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef
672 end if
673
674 else
675 lam(i) = zero
676 end if
677
678 enddo
679
680 if (present(n0)) n0 = nic * lam
681
682end subroutine size_dist_param_basic_vect
683
686elemental subroutine size_dist_param_ice_line(props, qic, nic, lam, n0)
687 type(mghydrometeorprops), intent(in) :: props
688 real(r8), intent(in) :: qic
689 real(r8), intent(inout) :: nic
690
691 real(r8), intent(out) :: lam
692 real(r8):: miu_ice,tx1,tx2, aux
693 real(r8), intent(out), optional :: n0
694 logical, parameter :: ice_sep=.true.
695
696 if (qic > qsmall) then
697
698 ! add upper limit to in-cloud number concentration to prevent
699 ! numerical error
700 if (limiter_is_on(props%min_mean_mass)) then
701 nic = min(nic, qic / props%min_mean_mass)
702 end if
703
704 ! lambda = (c n/q)^(1/d)
705 lam = (props%shape_coef * nic/qic)**(1._r8/props%eff_dim)
706 if (ice_sep) then
707 miu_ice = max(min(0.008_r8*(lam*0.01)**0.87_r8, 10.0_r8), 0.1_r8)
708 tx1 = 1.0_r8 + miu_ice
709 tx2 = 1.0_r8 / gamma(tx1)
710 aux = (gamma(tx1+3.0_r8)*tx2) ** (1.0_r8/3.0_r8)
711 lam = lam*aux
712 else
713 aux = 1.0_r8
714 tx1 = 1.0_r8
715 tx2 = 1.0_r8
716 end if
717 if (present(n0)) n0 = nic * lam**tx1*tx2
718
719 ! check for slope
720 ! adjust vars
721 if (lam < props%lambda_bounds(1)*aux) then
722 lam = props%lambda_bounds(1)
723 nic = lam**(props%eff_dim) * qic/props%shape_coef
724 if (present(n0)) n0 = nic * lam
725 else if (lam > props%lambda_bounds(2)*aux) then
726 lam = props%lambda_bounds(2)
727 nic = lam**(props%eff_dim) * qic/props%shape_coef
728 if (present(n0)) n0 = nic * lam
729 end if
730
731 else
732 lam = 0.0_r8
733 end if
734
735
736end subroutine size_dist_param_ice_line
737
740subroutine size_dist_param_ice_vect(props, qic, nic, lam, mgncol, n0)
741
742 type (mghydrometeorprops), intent(in) :: props
743 integer, intent(in) :: mgncol
744 real(r8), dimension(mgncol), intent(in) :: qic
745 real(r8), dimension(mgncol), intent(inout) :: nic
746 real(r8), dimension(mgncol), intent(out) :: lam
747 real(r8), dimension(mgncol), intent(out), optional :: n0
748 real(r8) :: miu_ice,tx1,tx2, aux
749 integer :: i
750 logical, parameter :: ice_sep=.true.
751 do i=1,mgncol
752
753 if (qic(i) > qsmall) then
754
755 ! add upper limit to in-cloud number concentration to prevent
756 ! numerical error
757 if (limiter_is_on(props%min_mean_mass)) then
758 nic(i) = min(nic(i), qic(i) / props%min_mean_mass)
759 end if
760
761 ! lambda = (c n/q)^(1/d)
762 lam(i) = (props%shape_coef * nic(i)/qic(i))**(1._r8/props%eff_dim)
763 if (ice_sep) then
764 miu_ice = max(min(0.008_r8*(lam(i)*0.01)**0.87_r8, 10.0_r8), 0.1_r8)
765 tx1 = 1.0_r8 + miu_ice
766 tx2 = 1.0_r8 / gamma(tx1)
767 aux = (gamma(tx1+3.0_r8)*tx2) ** (1.0_r8/3.0_r8)
768 lam(i) = lam(i)*aux
769 else
770 aux = 1.0_r8
771 tx1 = 1.0_r8
772 tx2 = 1.0_r8
773 end if
774 if (present(n0)) n0(i) = nic(i) * lam(i)**tx1*tx2
775
776 ! check for slope
777 ! adjust vars
778 if (lam(i) < props%lambda_bounds(1)*aux) then
779 lam(i) = props%lambda_bounds(1)
780 nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef
781 if (present(n0)) n0(i) = nic(i) * lam(i)
782 else if (lam(i) > props%lambda_bounds(2)*aux) then
783 lam(i) = props%lambda_bounds(2)
784 nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef
785 if (present(n0)) n0(i) = nic(i) * lam(i)
786 end if
787
788 else
789 lam(i) = 0.0_r8
790 end if
791
792 enddo
793
794end subroutine size_dist_param_ice_vect
795
800real(r8) elemental function avg_diameter(q, n, rho_air, rho_sub)
801 real(r8), intent(in) :: q
802 real(r8), intent(in) :: n
803 real(r8), intent(in) :: rho_air
804 real(r8), intent(in) :: rho_sub
805
806 avg_diameter = (pi * rho_sub * n/(q*rho_air))**(-oneo3)
807
808end function avg_diameter
809
813elemental function var_coef_r8(relvar, a) result(res)
814 real(r8), intent(in) :: relvar
815 real(r8), intent(in) :: a
816 real(r8) :: res
817
818 res = rising_factorial(relvar, a) / relvar**a
819
820end function var_coef_r8
821
825elemental function var_coef_integer(relvar, a) result(res)
826 real(r8), intent(in) :: relvar
827 integer, intent(in) :: a
828 real(r8) :: res
829
830 res = rising_factorial(relvar, a) / relvar**a
831
832end function var_coef_integer
833
834!========================================================================
835!MICROPHYSICAL PROCESS CALCULATIONS
836!========================================================================
837!========================================================================
842subroutine ice_deposition_sublimation(t, qv, qi, ni, &
843 icldm, rho, dv,qvl, qvi, &
844 berg, vap_dep, ice_sublim, mgncol)
845
846 !INPUT VARS:
847 !===============================================
848 integer, intent(in) :: mgncol
849 real(r8), dimension(mgncol), intent(in) :: t
850 real(r8), dimension(mgncol), intent(in) :: qv
851 real(r8), dimension(mgncol), intent(in) :: qi
852 real(r8), dimension(mgncol), intent(in) :: ni
853 real(r8), dimension(mgncol), intent(in) :: icldm
854 real(r8), dimension(mgncol), intent(in) :: rho
855 real(r8), dimension(mgncol), intent(in) :: dv
856 real(r8), dimension(mgncol), intent(in) :: qvl
857 real(r8), dimension(mgncol), intent(in) :: qvi
858
859 !OUTPUT VARS:
860 !===============================================
861 real(r8), dimension(mgncol), intent(out) :: vap_dep !ice deposition (cell-ave value)
862 real(r8), dimension(mgncol), intent(out) :: ice_sublim !ice sublimation (cell-ave value)
863 real(r8), dimension(mgncol), intent(out) :: berg !bergeron enhancement (cell-ave value)
864
865 !INTERNAL VARS:
866 !===============================================
867 real(r8) :: ab
868 real(r8) :: epsi
869 real(r8) :: qiic
870 real(r8) :: niic
871 real(r8) :: lami
872 real(r8) :: n0i
873 real(r8) :: tx1
874 integer :: i
875
876 do i=1,mgncol
877 if (qi(i)>=qsmall) then
878
879 !GET IN-CLOUD qi, ni
880 !===============================================
881 tx1 = one / icldm(i)
882 qiic = qi(i) * tx1
883 niic = ni(i) * tx1
884
885 !Compute linearized condensational heating correction
886 ab = calc_ab(t(i), qvi(i), xxls)
887 !Get slope and intercept of gamma distn for ice.
888! call size_dist_param_basic(mg_ice_props, qiic, niic, lami, n0i)
889 call size_dist_param_ice(mg_ice_props, qiic, niic, lami, n0i)
890 !Get depletion timescale=1/eps
891 epsi = twopi*n0i*rho(i)*dv(i)/(lami*lami)
892
893 !Compute deposition/sublimation
894 vap_dep(i) = epsi/ab*(qv(i) - qvi(i))
895
896 !Make this a grid-averaged quantity
897 vap_dep(i) = vap_dep(i)*icldm(i)
898
899 !Split into deposition or sublimation.
900 if (t(i) < tmelt .and. vap_dep(i) > zero) then
901 ice_sublim(i) = zero
902 else
903 ! make ice_sublim negative for consistency with other evap/sub processes
904 ice_sublim(i) = min(vap_dep(i), zero)
905 vap_dep(i) = zero
906 end if
907
908 !sublimation occurs @ any T. Not so for berg.
909 if (t(i) < tmelt) then
910
911 !Compute bergeron rate assuming cloud for whole step.
912 berg(i) = max(epsi/ab*(qvl(i) - qvi(i)), zero)
913 else !T>frz
914 berg(i) = zero
915 end if !T<frz
916
917 else !where qi<qsmall
918 berg(i) = zero
919 vap_dep(i) = zero
920 ice_sublim(i) = zero
921 end if !qi>qsmall
922 enddo
923end subroutine ice_deposition_sublimation
924
925!========================================================================
930subroutine kk2000_liq_autoconversion(microp_uniform, qcic, &
931 ncic, rho, relvar, prc, nprc, nprc1, mgncol)
932
933 integer, intent(in) :: mgncol
934 logical, intent(in) :: microp_uniform
935
936 real(r8), dimension(mgncol), intent(in) :: qcic
937 real(r8), dimension(mgncol), intent(in) :: ncic
938 real(r8), dimension(mgncol), intent(in) :: rho
939
940 real(r8), dimension(mgncol), intent(in) :: relvar
941
942 real(r8), dimension(mgncol), intent(out) :: prc
943 real(r8), dimension(mgncol), intent(out) :: nprc
944 real(r8), dimension(mgncol), intent(out) :: nprc1
945
946 real(r8), dimension(mgncol) :: prc_coef
947 integer :: i
948
949 ! Take variance into account, or use uniform value.
950 if (.not. microp_uniform) then
951 prc_coef(:) = var_coef(relvar(:), 2.47_r8)
952 else
953 prc_coef(:) = one
954 end if
955
956 do i=1,mgncol
957 if (qcic(i) >= icsmall) then
958
959 ! nprc is increase in rain number conc due to autoconversion
960 ! nprc1 is decrease in cloud droplet conc due to autoconversion
961
962 ! assume exponential sub-grid distribution of qc, resulting in additional
963 ! factor related to qcvar below
964 ! switch for sub-columns, don't include sub-grid qc
965
966 prc(i) = prc_coef(i) * &
967 0.01_r8 * 1350._r8 * qcic(i)**2.47_r8 * (ncic(i)*1.e-6_r8*rho(i))**(-1.1_r8)
968 nprc(i) = prc(i) * (one/droplet_mass_25um)
969 nprc1(i) = prc(i)*ncic(i)/qcic(i)
970
971 else
972 prc(i) = zero
973 nprc(i) = zero
974 nprc1(i) = zero
975 end if
976 enddo
977end subroutine kk2000_liq_autoconversion
978
979 !========================================================================
982subroutine sb2001v2_liq_autoconversion(pgam,qc,nc,qr,rho,relvar,au,nprc,nprc1,mgncol)
983 !
984 ! ---------------------------------------------------------------------
985 ! AUTO_SB: calculates the evolution of mass- and number mxg-ratio for
986 ! drizzle drops due to autoconversion. The autoconversion rate assumes
987 ! f(x)=A*x**(nu_c)*exp(-Bx) in drop MASS x.
988
989 ! Code from Hugh Morrison, Sept 2014
990
991 ! autoconversion
992 ! use simple lookup table of dnu values to get mass spectral shape parameter
993 ! equivalent to the size spectral shape parameter pgam
994
995 integer, intent(in) :: mgncol
996
997 real(r8), dimension(mgncol), intent (in) :: pgam
998 real(r8), dimension(mgncol), intent (in) :: qc ! = qc (cld water mixing ratio)
999 real(r8), dimension(mgncol), intent (in) :: nc ! = nc (cld water number conc /kg)
1000 real(r8), dimension(mgncol), intent (in) :: qr ! = qr (rain water mixing ratio)
1001 real(r8), dimension(mgncol), intent (in) :: rho ! = rho : density profile
1002 real(r8), dimension(mgncol), intent (in) :: relvar
1003
1004 real(r8), dimension(mgncol), intent (out) :: au ! = prc autoconversion rate
1005 real(r8), dimension(mgncol), intent (out) :: nprc1 ! = number tendency
1006 real(r8), dimension(mgncol), intent (out) :: nprc ! = number tendency fixed size for rain
1007
1008 ! parameters for droplet mass spectral shape,
1009 ! used by Seifert and Beheng (2001)
1010 ! warm rain scheme only (iparam = 1)
1011 real(r8), parameter :: dnu(16) = [0._r8,-0.557_r8,-0.430_r8,-0.307_r8, &
1012 -0.186_r8,-0.067_r8,0.050_r8,0.167_r8,0.282_r8,0.397_r8,0.512_r8, &
1013 0.626_r8,0.739_r8,0.853_r8,0.966_r8,0.966_r8]
1014
1015 ! parameters for Seifert and Beheng (2001) autoconversion/accretion
1016 real(r8), parameter :: kc = 9.44e9_r8
1017 real(r8), parameter :: kr = 5.78e3_r8
1018 real(r8), parameter :: auf = kc / (20._r8*2.6e-7_r8) * 1000._r8, &
1019 con_nprc1 = two/2.6e-7_r8*1000._r8
1020 real(r8) :: dum, dum1, nu, pra_coef, tx1, tx2, tx3, tx4
1021 integer :: dumi, i
1022
1023 do i=1,mgncol
1024
1025 pra_coef = var_coef(relvar(i), 2.47_r8)
1026 if (qc(i) > qsmall) then
1027 dumi = max(1, min(int(pgam(i)), 15))
1028 nu = dnu(dumi) + (dnu(dumi+1)-dnu(dumi))* (pgam(i)-dumi)
1029
1030 !Anning fixed a bug here for FV3GFS 10/13/2017
1031 dum = max(one-qc(i)/(qc(i)+qr(i)), zero)
1032 tx1 = dum**0.68_r8
1033 tx2 = one - tx1
1034 dum1 = 600._r8 * tx1 * tx2 * tx2 * tx2 ! Moorthi
1035! dum1 = 600._r8*dum**0.68_r8*(one-dum**0.68_r8)**3
1036
1037 tx1 = nu + one
1038 tx2 = 0.001_r8 * rho(i) * qc(i)
1039 tx3 = tx2 * tx2 / (rho(i)*nc(i)*1.e-6_r8)
1040 tx2 = tx3 * tx3
1041 tx3 = one - dum
1042 au(i) = auf * (nu+two) * (nu+four) * tx2 &
1043 * (one+dum1/(tx3*tx3)) / (tx1*tx1*rho(i))
1044
1045! au(i) = kc/(20._r8*2.6e-7_r8)* &
1046! (nu+2._r8)*(nu+4._r8)/(nu+1._r8)**2._r8* &
1047! (rho(i)*qc(i)/1000._r8)**4._r8/(rho(i)*nc(i)/1.e6_r8)**2._r8* &
1048! (1._r8+dum1/(1._r8-dum)**2)*1000._r8 / rho(i)
1049
1050! nprc1(i) = au(i) * two / 2.6e-7_r8 * 1000._r8
1051! nprc(i) = au(i) / droplet_mass_40um
1052 nprc1(i) = au(i) * con_nprc1
1053 nprc(i) = au(i) * droplet_mass_40umi
1054 else
1055 au(i) = zero
1056 nprc1(i) = zero
1057 nprc(i) = zero
1058 end if
1059
1060 enddo
1061
1062 end subroutine sb2001v2_liq_autoconversion
1063
1064!========================================================================
1067 subroutine liu_liq_autoconversion(pgam,qc,nc,qr,rho,relvar, &
1068 au,nprc,nprc1,mgncol)
1069
1070
1071
1072 integer, intent(in) :: mgncol
1073
1074 real(r8), dimension(mgncol), intent (in) :: pgam
1075 real(r8), dimension(mgncol), intent (in) :: qc
1076 real(r8), dimension(mgncol), intent (in) :: nc
1077 real(r8), dimension(mgncol), intent (in) :: qr
1078 real(r8), dimension(mgncol), intent (in) :: rho
1079 real(r8), dimension(mgncol), intent (in) :: relvar
1080
1081 real(r8), dimension(mgncol), intent (out) :: au
1082 real(r8), dimension(mgncol), intent (out) :: nprc1
1083 real(r8), dimension(mgncol), intent (out) :: nprc
1084 real(r8) :: xs,lw, nw, beta6
1085! real(r8), parameter :: dcrit=1.0e-6, miu_disp=1.
1086! real(r8), parameter :: dcrit=1.0e-3, miu_disp=1.
1087 real(r8), parameter :: dcrit = 2.0e-3, miu_disp = 0.8, &
1088 con_nprc1 = two/2.6e-7_r8*1000._r8
1089 integer :: i
1090
1091 do i=1,mgncol
1092 if (qc(i) > qsmall) then
1093 xs = one / (one+pgam(i))
1094 beta6 = (one+three*xs)*(one+four*xs)*(one+five*xs) &
1095 / ((one+xs)*(one+xs+xs))
1096 lw = 1.0e-3_r8 * qc(i) * rho(i)
1097 nw = nc(i) * rho(i) * 1.e-6_r8
1098
1099 xs = min(20.0_r8, 1.03e16_r8*(lw*lw)/(nw*sqrt(nw)))
1100 au(i) = 1.1e10_r8*beta6*lw*lw*lw &
1101 * (one-exp(-(xs**miu_disp))) / nw
1102 au(i) = au(i)*1.0e3_r8/rho(i)
1103 au(i) = au(i) * gamma(two+relvar(i)) &
1104 / (gamma(relvar(i))*(relvar(i)*relvar(i)))
1105
1106 au(i) = au(i) * dcrit
1107! nprc1(i)= au(i) * (two/2.6e-7_r8*1000._r8)
1108! nprc(i) = au(i) / droplet_mass_40um
1109 nprc1(i)= au(i) * con_nprc1
1110 nprc(i) = au(i) * droplet_mass_40umi
1111 else
1112 au(i) = zero
1113 nprc1(i) = zero
1114 nprc(i) = zero
1115 end if
1116 enddo
1117
1118 end subroutine liu_liq_autoconversion
1119
1120
1121!========================================================================
1122!SB2001 Accretion V2
1124subroutine sb2001v2_accre_cld_water_rain(qc,nc,qr,rho,relvar,pra,npra,mgncol)
1125 !
1126 ! ---------------------------------------------------------------------
1127 ! ACCR_SB calculates the evolution of mass mxng-ratio due to accretion
1128 ! and self collection following Seifert & Beheng (2001).
1129 !
1130
1131 integer, intent(in) :: mgncol
1132
1133 real(r8), dimension(mgncol), intent (in) :: qc ! = qc (cld water mixing ratio)
1134 real(r8), dimension(mgncol), intent (in) :: nc ! = nc (cld water number conc /kg)
1135 real(r8), dimension(mgncol), intent (in) :: qr ! = qr (rain water mixing ratio)
1136 real(r8), dimension(mgncol), intent (in) :: rho ! = rho : density profile
1137 real(r8), dimension(mgncol), intent (in) :: relvar
1138
1139 ! Output tendencies
1140 real(r8), dimension(mgncol), intent(out) :: pra ! MMR
1141 real(r8), dimension(mgncol), intent(out) :: npra ! Number
1142
1143 ! parameters for Seifert and Beheng (2001) autoconversion/accretion
1144 real(r8), parameter :: kc = 9.44e9_r8
1145 real(r8), parameter :: kr = 5.78e3_r8
1146
1147 real(r8) :: dum, dum1, tx1, tx2
1148 integer :: i
1149
1150 ! accretion
1151
1152 do i =1,mgncol
1153
1154 if (qc(i) > qsmall) then
1155 dum = one - qc(i)/(qc(i)+qr(i))
1156 tx1 = dum / (dum+5.e-4_r8)
1157 dum1 = tx1 * tx1
1158 dum1 = dum1 * dum1
1159 pra(i) = kr*rho(i)*0.001_r8*qc(i)*qr(i)*dum1
1160
1161 npra(i) = pra(i) * nc(i) / qc(i)
1162
1163! npra(i) = pra(i)*rho(i)*0.001_r8*(nc(i)*rho(i)*1.e-6_r8)/ &
1164! (qc(i)*rho(i)*0.001_r8)*1.e6_r8 / rho(i)
1165 else
1166 pra(i) = zero
1167 npra(i) = zero
1168 end if
1169
1170 enddo
1171
1172 end subroutine sb2001v2_accre_cld_water_rain
1173
1174!========================================================================
1175! Autoconversion of cloud ice to snow
1176! similar to Ferrier (1994)
1180subroutine ice_autoconversion(t, qiic, lami, n0i, dcs, ac_time, prci, nprci, mgncol)
1181
1182 integer, intent(in) :: mgncol
1183 real(r8), dimension(mgncol), intent(in) :: t
1184 real(r8), dimension(mgncol), intent(in) :: qiic
1185 real(r8), dimension(mgncol), intent(in) :: lami
1186 real(r8), dimension(mgncol), intent(in) :: n0i
1187 real(r8), intent(in) :: dcs
1188 real(r8), dimension(mgncol), intent(in) :: ac_time
1189
1190 real(r8), dimension(mgncol), intent(out) :: prci
1191 real(r8), dimension(mgncol), intent(out) :: nprci
1192
1193 ! Assume autoconversion timescale of 180 seconds.
1194
1195 ! Average mass of an ice particle.
1196 real(r8) :: m_ip
1197 ! Ratio of autoconversion diameter to average diameter.
1198 real(r8) :: d_rat
1199 integer :: i
1200
1201 do i=1,mgncol
1202 if (t(i) <= tmelt .and. qiic(i) >= qsmall) then
1203
1204 d_rat = lami(i)*dcs
1205
1206 ! Rate of ice particle conversion (number).
1207 nprci(i) = n0i(i)/(lami(i)*ac_time(i))*exp(-d_rat)
1208
1209 m_ip = rhoi * pio6 / (lami(i)*lami(i)*lami(i))
1210
1211! m_ip = (rhoi*pi/6._r8) / lami(i)**3
1212
1213 ! Rate of mass conversion.
1214 ! Note that this is:
1215 ! m n (d^3 + 3 d^2 + 6 d + 6)
1216 prci(i) = m_ip * nprci(i) * (((d_rat + three)*d_rat + six)*d_rat + six)
1217
1218 else
1219 prci(i) = zero
1220 nprci(i) = zero
1221 end if
1222 enddo
1223end subroutine ice_autoconversion
1224!===================================
1225! Anning Cheng 10/5/2017 added GMAO ice autoconversion
1228subroutine gmao_ice_autoconversion(t, qiic, niic, lami, n0i, &
1229 dcs, ac_time, prci, nprci, mgncol)
1230
1231 integer, intent(in) :: mgncol
1232 real(r8), dimension(mgncol), intent(in) :: t
1233 real(r8), dimension(mgncol), intent(in) :: qiic
1234 real(r8), dimension(mgncol), intent(in) :: niic
1235 real(r8), dimension(mgncol), intent(in) :: lami
1236 real(r8), dimension(mgncol), intent(in) :: n0i
1237 real(r8), dimension(mgncol), intent(in) :: ac_time
1238 real(r8), intent(in) :: dcs
1239
1240 real(r8), dimension(mgncol), intent(out) :: prci
1241 real(r8), dimension(mgncol), intent(out) :: nprci
1242
1243
1244 real(r8) :: m_ip, tx1, tx2
1245 integer :: i
1246 do i=1,mgncol
1247 if (t(i) <= tmelt .and. qiic(i) >= qsmall) then
1248 m_ip = max(min(0.008_r8*(lami(i)*0.01)**0.87_r8, &
1249 10.0_r8), 0.1_r8)
1250 tx1 = lami(i)*dcs
1251 tx2 = one / ac_time(i)
1252 nprci(i) = niic(i)*tx2 * (one - gamma_incomp(m_ip, tx1))
1253 prci(i) = qiic(i)*tx2 * (one - gamma_incomp(m_ip+three, tx1))
1254 else
1255 prci(i) = zero
1256 nprci(i) = zero
1257 end if
1258 enddo
1259end subroutine gmao_ice_autoconversion
1260!===================================
1261! immersion freezing (Bigg, 1953)
1262!===================================
1265subroutine immersion_freezing(microp_uniform, t, pgam, lamc, &
1266 qcic, ncic, relvar, mnuccc, nnuccc, mgncol)
1267
1268 integer, intent(in) :: mgncol
1269 logical, intent(in) :: microp_uniform
1270
1271 ! Temperature
1272 real(r8), dimension(mgncol), intent(in) :: t
1273
1274 ! Cloud droplet size distribution parameters
1275 real(r8), dimension(mgncol), intent(in) :: pgam
1276 real(r8), dimension(mgncol), intent(in) :: lamc
1277
1278 ! MMR and number concentration of in-cloud liquid water
1279 real(r8), dimension(mgncol), intent(in) :: qcic
1280 real(r8), dimension(mgncol), intent(in) :: ncic
1281
1282 ! Relative variance of cloud water
1283 real(r8), dimension(mgncol), intent(in) :: relvar
1284
1285 ! Output tendencies
1286 real(r8), dimension(mgncol), intent(out) :: mnuccc ! MMR
1287 real(r8), dimension(mgncol), intent(out) :: nnuccc ! Number
1288
1289 ! Coefficients that will be omitted for sub-columns
1290 real(r8), dimension(mgncol) :: dum
1291 real(r8) :: tx1
1292 integer :: i
1293
1294 if (.not. microp_uniform) then
1295 dum(:) = var_coef(relvar, 2)
1296 else
1297 dum(:) = one
1298 end if
1299 do i=1,mgncol
1300
1301 if (qcic(i) >= qsmall .and. t(i) < 269.15_r8) then
1302
1303 tx1 = one / (lamc(i) * lamc(i) * lamc(i))
1304 nnuccc(i) = pio6*ncic(i)*rising_factorial(pgam(i)+one, 3) * &
1305 bimm*(exp(aimm*(tmelt - t(i)))-one) * tx1
1306
1307 mnuccc(i) = dum(i) * nnuccc(i) * pio6 * rhow * &
1308 rising_factorial(pgam(i)+four, 3) * tx1
1309
1310 else
1311 mnuccc(i) = zero
1312 nnuccc(i) = zero
1313 end if ! qcic > qsmall and t < 4 deg C
1314 enddo
1315
1316end subroutine immersion_freezing
1317
1321subroutine contact_freezing (microp_uniform, t, p, rndst, nacon, &
1322 pgam, lamc, qcic, ncic, relvar, mnucct, nnucct, mgncol, mdust)
1323
1324 logical, intent(in) :: microp_uniform
1325
1326 integer, intent(in) :: mgncol
1327 integer, intent(in) :: mdust
1328
1329 real(r8), dimension(mgncol), intent(in) :: t ! Temperature
1330 real(r8), dimension(mgncol), intent(in) :: p ! Pressure
1331 real(r8), dimension(mgncol, mdust), intent(in) :: rndst ! Radius (for multiple dust bins)
1332 real(r8), dimension(mgncol, mdust), intent(in) :: nacon ! Number (for multiple dust bins)
1333
1334 ! Size distribution parameters for cloud droplets
1335 real(r8), dimension(mgncol), intent(in) :: pgam
1336 real(r8), dimension(mgncol), intent(in) :: lamc
1337
1338 ! MMR and number concentration of in-cloud liquid water
1339 real(r8), dimension(mgncol), intent(in) :: qcic
1340 real(r8), dimension(mgncol), intent(in) :: ncic
1341
1342 ! Relative cloud water variance
1343 real(r8), dimension(mgncol), intent(in) :: relvar
1344
1345 ! Output tendencies
1346 real(r8), dimension(mgncol), intent(out) :: mnucct ! MMR
1347 real(r8), dimension(mgncol), intent(out) :: nnucct ! Number
1348
1349 real(r8) :: tcnt ! scaled relative temperature
1350 real(r8) :: viscosity ! temperature-specific viscosity (kg/m/s)
1351 real(r8) :: mfp ! temperature-specific mean free path (m)
1352
1353 ! Dimension these according to number of dust bins, inferred from rndst size
1354 real(r8) :: nslip(size(rndst,2)) ! slip correction factors
1355 real(r8) :: ndfaer(size(rndst,2)) ! aerosol diffusivities (m^2/sec)
1356
1357 ! Coefficients not used for subcolumns
1358 real(r8) :: dum, dum1, tx1
1359
1360 ! Common factor between mass and number.
1361 real(r8) :: contact_factor
1362
1363 integer :: i
1364
1365 do i = 1,mgncol
1366
1367 if (qcic(i) >= qsmall .and. t(i) < 269.15_r8) then
1368
1369 if (.not. microp_uniform) then
1370 dum = var_coef(relvar(i), four/three)
1371 dum1 = var_coef(relvar(i), oneo3)
1372 else
1373 dum = one
1374 dum1 = one
1375 endif
1376
1377 tcnt=(270.16_r8-t(i))**1.3_r8
1378 viscosity = 1.8e-5_r8*(t(i)/298.0_r8)**0.85_r8 ! Viscosity (kg/m/s)
1379 mfp = two*viscosity/ & ! Mean free path (m)
1380 (p(i)*sqrt( 8.0_r8*28.96e-3_r8/(pi*8.314409_r8*t(i)) ))
1381
1382 ! Note that these two are vectors.
1383 nslip = one+(mfp/rndst(i,:))*(1.257_r8+(0.4_r8*exp(-(1.1_r8*rndst(i,:)/mfp))))! Slip correction factor
1384
1385 ndfaer = 1.381e-23_r8*t(i)*nslip/(6._r8*pi*viscosity*rndst(i,:)) ! aerosol diffusivity (m2/s)
1386
1387 tx1 = one / lamc(i)
1388 contact_factor = dot_product(ndfaer,nacon(i,:)*tcnt) * pi * &
1389 ncic(i) * (pgam(i) + one) * tx1
1390
1391 mnucct(i) = dum * contact_factor * &
1392 pio3*rhow*rising_factorial(pgam(i)+two, 3) * tx1 * tx1 *tx1
1393
1394 nnucct(i) = (dum1+dum1) * contact_factor
1395
1396 else
1397
1398 mnucct(i) = zero
1399 nnucct(i) = zero
1400
1401 end if ! qcic > qsmall and t < 4 deg C
1402 end do
1403
1404end subroutine contact_freezing
1405
1408!===================================================================
1409! this is hard-wired for bs = 0.4 for now
1410! ignore self-collection of cloud ice
1411subroutine snow_self_aggregation(t, rho, asn, rhosn, qsic, nsic, nsagg, mgncol)
1412
1413 integer, intent(in) :: mgncol
1414
1415 real(r8), dimension(mgncol), intent(in) :: t ! Temperature
1416 real(r8), dimension(mgncol), intent(in) :: rho ! Density
1417 real(r8), dimension(mgncol), intent(in) :: asn ! fall speed parameter for snow
1418 real(r8), intent(in) :: rhosn ! density of snow
1419
1420 ! In-cloud snow
1421 real(r8), dimension(mgncol), intent(in) :: qsic ! MMR
1422 real(r8), dimension(mgncol), intent(in) :: nsic ! Number
1423
1424 ! Output number tendency
1425 real(r8), dimension(mgncol), intent(out) :: nsagg
1426
1427 integer :: i
1428
1429 do i=1,mgncol
1430 if (qsic(i) >= qsmall .and. t(i) <= tmelt) then
1431 nsagg(i) = -1108._r8*eii/(four*720._r8*rhosn)*asn(i)*qsic(i)*nsic(i)*rho(i)*&
1432 ((qsic(i)/nsic(i))*(one/(rhosn*pi)))**((bs-one)*oneo3)
1433 else
1434 nsagg(i) = zero
1435 end if
1436 enddo
1437end subroutine snow_self_aggregation
1438
1441!===================================================================
1442! here use continuous collection equation with
1443! simple gravitational collection kernel
1444! ignore collisions between droplets/cloud ice
1445! since minimum size ice particle for accretion is 50 - 150 micron
1446subroutine accrete_cloud_water_snow(t, rho, asn, uns, mu, qcic, ncic, qsic, &
1447 pgam, lamc, lams, n0s, psacws, npsacws, mgncol)
1448
1449 integer, intent(in) :: mgncol
1450 real(r8), dimension(mgncol), intent(in) :: t ! Temperature
1451 real(r8), dimension(mgncol), intent(in) :: rho ! Density
1452 real(r8), dimension(mgncol), intent(in) :: asn ! Fallspeed parameter (snow)
1453 real(r8), dimension(mgncol), intent(in) :: uns ! Current fallspeed (snow)
1454 real(r8), dimension(mgncol), intent(in) :: mu ! Viscosity
1455
1456 ! In-cloud liquid water
1457 real(r8), dimension(mgncol), intent(in) :: qcic ! MMR
1458 real(r8), dimension(mgncol), intent(in) :: ncic ! Number
1459
1460 ! In-cloud snow
1461 real(r8), dimension(mgncol), intent(in) :: qsic ! MMR
1462
1463 ! Cloud droplet size parameters
1464 real(r8), dimension(mgncol), intent(in) :: pgam
1465 real(r8), dimension(mgncol), intent(in) :: lamc
1466
1467 ! Snow size parameters
1468 real(r8), dimension(mgncol), intent(in) :: lams
1469 real(r8), dimension(mgncol), intent(in) :: n0s
1470
1471 ! Output tendencies
1472 real(r8), dimension(mgncol), intent(out) :: psacws ! Mass mixing ratio
1473 real(r8), dimension(mgncol), intent(out) :: npsacws ! Number concentration
1474
1475 real(r8) :: dc0 ! Provisional mean droplet size
1476 real(r8) :: dum
1477 real(r8) :: eci ! collection efficiency for riming of snow by droplets
1478
1479 ! Fraction of cloud droplets accreted per second
1480 real(r8) :: accrete_rate
1481 integer :: i
1482
1483 ! ignore collision of snow with droplets above freezing
1484
1485 do i=1,mgncol
1486 if (qsic(i) >= qsmall .and. t(i) <= tmelt .and. qcic(i) >= qsmall) then
1487
1488 ! put in size dependent collection efficiency
1489 ! mean diameter of snow is area-weighted, since
1490 ! accretion is function of crystal geometric area
1491 ! collection efficiency is approximation based on stoke's law (Thompson et al. 2004)
1492
1493 dc0 = (pgam(i)+one)/lamc(i)
1494 dum = dc0*dc0*uns(i)*rhow*lams(i)/(9._r8*mu(i))
1495 eci = dum*dum / ((dum+0.4_r8)*(dum+0.4_r8))
1496
1497 eci = max(eci,zero)
1498 eci = min(eci,one)
1499
1500 ! no impact of sub-grid distribution of qc since psacws
1501 ! is linear in qc
1502 accrete_rate = (pi/four)*asn(i)*rho(i)*n0s(i)*eci*gamma_bs_plus3 / lams(i)**(bs+three)
1503 psacws(i) = accrete_rate*qcic(i)
1504 npsacws(i) = accrete_rate*ncic(i)
1505 else
1506 psacws(i) = zero
1507 npsacws(i) = zero
1508 end if
1509 enddo
1510end subroutine accrete_cloud_water_snow
1511
1514!===================================================================
1515! (Hallet-Mossop process) (from Cotton et al., 1986)
1516subroutine secondary_ice_production(t, psacws, msacwi, nsacwi, mgncol)
1517
1518 integer, intent(in) :: mgncol
1519 real(r8), dimension(mgncol), intent(in) :: t ! Temperature
1520
1521 ! Accretion of cloud water to snow tendencies
1522 real(r8), dimension(mgncol), intent(inout) :: psacws ! MMR
1523
1524 ! Output (ice) tendencies
1525 real(r8), dimension(mgncol), intent(out) :: msacwi ! MMR
1526 real(r8), dimension(mgncol), intent(out) :: nsacwi ! Number
1527 integer :: i
1528
1529 do i=1,mgncol
1530 if((t(i) < 270.16_r8) .and. (t(i) >= 268.16_r8)) then
1531 nsacwi(i) = 3.5e8_r8*(270.16_r8-t(i))/two*psacws(i)
1532 else if((t(i) < 268.16_r8) .and. (t(i) >= 265.16_r8)) then
1533 nsacwi(i) = 3.5e8_r8*(t(i)-265.16_r8)*oneo3*psacws(i)
1534 else
1535 nsacwi(i) = zero
1536 endif
1537 enddo
1538
1539 do i=1,mgncol
1540 msacwi(i) = min(nsacwi(i)*mi0, psacws(i))
1541 psacws(i) = psacws(i) - msacwi(i)
1542 enddo
1543end subroutine secondary_ice_production
1544
1547!===================================================================
1548! formula from ikawa and saito, 1991, used by reisner et al., 1998
1549subroutine accrete_rain_snow(t, rho, umr, ums, unr, uns, qric, qsic, &
1550 lamr, n0r, lams, n0s, pracs, npracs, mgncol)
1551
1552 integer, intent(in) :: mgncol
1553
1554 real(r8), dimension(mgncol), intent(in) :: t ! Temperature
1555 real(r8), dimension(mgncol), intent(in) :: rho ! Density
1556
1557 ! Fallspeeds
1558 ! mass-weighted
1559 real(r8), dimension(mgncol), intent(in) :: umr ! rain
1560 real(r8), dimension(mgncol), intent(in) :: ums ! snow
1561 ! number-weighted
1562 real(r8), dimension(mgncol), intent(in) :: unr ! rain
1563 real(r8), dimension(mgncol), intent(in) :: uns ! snow
1564
1565 ! In cloud MMRs
1566 real(r8), dimension(mgncol), intent(in) :: qric ! rain
1567 real(r8), dimension(mgncol), intent(in) :: qsic ! snow
1568
1569 ! Size distribution parameters
1570 ! rain
1571 real(r8), dimension(mgncol), intent(in) :: lamr
1572 real(r8), dimension(mgncol), intent(in) :: n0r
1573 ! snow
1574 real(r8), dimension(mgncol), intent(in) :: lams
1575 real(r8), dimension(mgncol), intent(in) :: n0s
1576
1577 ! Output tendencies
1578 real(r8), dimension(mgncol), intent(out) :: pracs ! MMR
1579 real(r8), dimension(mgncol), intent(out) :: npracs ! Number
1580
1581 ! Collection efficiency for accretion of rain by snow
1582 real(r8), parameter :: ecr = one
1583
1584 ! Ratio of average snow diameter to average rain diameter.
1585 real(r8) :: d_rat
1586 ! Common factor between mass and number expressions
1587 real(r8) :: common_factor
1588 real(r8) :: tx1, tx2
1589 integer :: i
1590
1591 do i=1,mgncol
1592 if (qric(i) >= icsmall .and. qsic(i) >= icsmall .and. t(i) <= tmelt) then
1593
1594 tx2 = lamr(i)*lamr(i)*lamr(i)
1595
1596 common_factor = pi*ecr*rho(i)*n0r(i)*n0s(i) / (tx2 * lams(i))
1597
1598 d_rat = lamr(i)/lams(i)
1599
1600 tx1 = 1.2_r8*umr(i)-0.95_r8*ums(i)
1601 pracs(i) = common_factor*pi*rhow* &
1602 sqrt(tx1*tx1 + 0.08_r8*ums(i)*umr(i)) * &
1603 ((half*d_rat + two)*d_rat + five) / tx2
1604
1605 tx1 = unr(i)-uns(i)
1606 npracs(i) = common_factor*half * &
1607 sqrt(1.7_r8*tx1*tx1 + 0.3_r8*unr(i)*uns(i)) * &
1608 ((d_rat + one)*d_rat + one)
1609
1610 else
1611 pracs(i) = zero
1612 npracs(i) = zero
1613 end if
1614 enddo
1615end subroutine accrete_rain_snow
1616
1619!===================================================================
1620! follows from Bigg (1953)
1621subroutine heterogeneous_rain_freezing(t, qric, nric, lamr, mnuccr, nnuccr, mgncol)
1622
1623 integer, intent(in) :: mgncol
1624 real(r8), dimension(mgncol), intent(in) :: t ! Temperature
1625
1626 ! In-cloud rain
1627 real(r8), dimension(mgncol), intent(in) :: qric ! MMR
1628 real(r8), dimension(mgncol), intent(in) :: nric ! Number
1629 real(r8), dimension(mgncol), intent(in) :: lamr ! size parameter
1630
1631 ! Output tendencies
1632 real(r8), dimension(mgncol), intent(out) :: mnuccr ! MMR
1633 real(r8), dimension(mgncol), intent(out) :: nnuccr ! Number
1634 real(r8) :: tx1
1635 integer :: i
1636
1637 do i=1,mgncol
1638
1639 if (t(i) < 269.15_r8 .and. qric(i) >= qsmall) then
1640 tx1 = pi / (lamr(i)*lamr(i)*lamr(i))
1641 nnuccr(i) = nric(i)*bimm* (exp(aimm*(tmelt - t(i)))-one) * tx1
1642
1643 mnuccr(i) = nnuccr(i) * 20._r8*rhow * tx1
1644
1645 else
1646 mnuccr(i) = zero
1647 nnuccr(i) = zero
1648 end if
1649 enddo
1650end subroutine heterogeneous_rain_freezing
1651
1655! gravitational collection kernel, droplet fall speed neglected
1656subroutine accrete_cloud_water_rain(microp_uniform, qric, qcic, &
1657 ncic, relvar, accre_enhan, pra, npra, mgncol)
1658
1659 logical, intent(in) :: microp_uniform
1660 integer, intent(in) :: mgncol
1661 ! In-cloud rain
1662 real(r8), dimension(mgncol), intent(in) :: qric ! MMR
1663
1664 ! Cloud droplets
1665 real(r8), dimension(mgncol), intent(in) :: qcic ! MMR
1666 real(r8), dimension(mgncol), intent(in) :: ncic ! Number
1667
1668 ! SGS variability
1669 real(r8), dimension(mgncol), intent(in) :: relvar
1670 real(r8), dimension(mgncol), intent(in) :: accre_enhan
1671
1672 ! Output tendencies
1673 real(r8), dimension(mgncol), intent(out) :: pra ! MMR
1674 real(r8), dimension(mgncol), intent(out) :: npra ! Number
1675
1676 ! Coefficient that varies for subcolumns
1677 real(r8), dimension(mgncol) :: pra_coef
1678
1679 integer :: i
1680
1681 if (.not. microp_uniform) then
1682 pra_coef(:) = accre_enhan * var_coef(relvar(:), 1.15_r8)
1683 else
1684 pra_coef(:) = one
1685 end if
1686
1687 do i=1,mgncol
1688
1689 if (qric(i) >= qsmall .and. qcic(i) >= qsmall) then
1690
1691 ! include sub-grid distribution of cloud water
1692 pra(i) = pra_coef(i) * 67._r8*(qcic(i)*qric(i))**1.15_r8
1693
1694 npra(i) = pra(i)*ncic(i)/qcic(i)
1695
1696 else
1697 pra(i) = zero
1698 npra(i) = zero
1699 end if
1700 end do
1701end subroutine accrete_cloud_water_rain
1702
1706subroutine self_collection_rain(rho, qric, nric, nragg, mgncol)
1707
1708 integer, intent(in) :: mgncol
1709 real(r8), dimension(mgncol), intent(in) :: rho ! Air density
1710
1711 ! Rain
1712 real(r8), dimension(mgncol), intent(in) :: qric ! MMR
1713 real(r8), dimension(mgncol), intent(in) :: nric ! Number
1714
1715 ! Output number tendency
1716 real(r8), dimension(mgncol), intent(out) :: nragg
1717
1718 integer :: i
1719
1720 do i=1,mgncol
1721 if (qric(i) >= qsmall) then
1722 nragg(i) = -8._r8*nric(i)*qric(i)*rho(i)
1723 else
1724 nragg(i) = zero
1725 end if
1726 enddo
1727end subroutine self_collection_rain
1728
1731!===================================================================
1732! For this calculation, it is assumed that the Vs >> Vi
1733! and Ds >> Di for continuous collection
1734subroutine accrete_cloud_ice_snow(t, rho, asn, qiic, niic, qsic, &
1735 lams, n0s, prai, nprai, mgncol)
1736
1737 integer, intent(in) :: mgncol
1738 real(r8), dimension(mgncol), intent(in) :: t ! Temperature
1739 real(r8), dimension(mgncol), intent(in) :: rho ! Density
1740
1741 real(r8), dimension(mgncol), intent(in) :: asn ! Snow fallspeed parameter
1742
1743 ! Cloud ice
1744 real(r8), dimension(mgncol), intent(in) :: qiic ! MMR
1745 real(r8), dimension(mgncol), intent(in) :: niic ! Number
1746
1747 real(r8), dimension(mgncol), intent(in) :: qsic ! Snow MMR
1748
1749 ! Snow size parameters
1750 real(r8), dimension(mgncol), intent(in) :: lams
1751 real(r8), dimension(mgncol), intent(in) :: n0s
1752
1753 ! Output tendencies
1754 real(r8), dimension(mgncol), intent(out) :: prai ! MMR
1755 real(r8), dimension(mgncol), intent(out) :: nprai ! Number
1756
1757 ! Fraction of cloud ice particles accreted per second
1758 real(r8) :: accrete_rate
1759
1760 integer :: i
1761
1762 do i=1,mgncol
1763 if (qsic(i) >= qsmall .and. qiic(i) >= qsmall .and. t(i) <= tmelt) then
1764
1765 accrete_rate = (pi/four) * eii * asn(i) * rho(i) * n0s(i) * gamma_bs_plus3 &
1766 / lams(i)**(bs+three)
1767
1768 prai(i) = accrete_rate * qiic(i)
1769 nprai(i) = accrete_rate * niic(i)
1770
1771 else
1772 prai(i) = zero
1773 nprai(i) = zero
1774 end if
1775 enddo
1776end subroutine accrete_cloud_ice_snow
1777
1780!===================================================================
1781! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
1782! in-cloud condensation/deposition of rain and snow is neglected
1783! except for transfer of cloud water to snow through bergeron process
1784subroutine evaporate_sublimate_precip(t, rho, dv, mu, sc, q, qvl, qvi, &
1785 lcldm, precip_frac, arn, asn, qcic, qiic, qric, qsic, lamr, n0r, lams, n0s, &
1786 pre, prds, am_evp_st, mgncol)
1787
1788 integer, intent(in) :: mgncol
1789
1790 real(r8), dimension(mgncol), intent(in) :: t ! temperature
1791 real(r8), dimension(mgncol), intent(in) :: rho ! air density
1792 real(r8), dimension(mgncol), intent(in) :: dv ! water vapor diffusivity
1793 real(r8), dimension(mgncol), intent(in) :: mu ! viscosity
1794 real(r8), dimension(mgncol), intent(in) :: sc ! schmidt number
1795 real(r8), dimension(mgncol), intent(in) :: q ! humidity
1796 real(r8), dimension(mgncol), intent(in) :: qvl ! saturation humidity (water)
1797 real(r8), dimension(mgncol), intent(in) :: qvi ! saturation humidity (ice)
1798 real(r8), dimension(mgncol), intent(in) :: lcldm ! liquid cloud fraction
1799 real(r8), dimension(mgncol), intent(in) :: precip_frac ! precipitation fraction (maximum overlap)
1800
1801 ! fallspeed parameters
1802 real(r8), dimension(mgncol), intent(in) :: arn ! rain
1803 real(r8), dimension(mgncol), intent(in) :: asn ! snow
1804
1805 ! In-cloud MMRs
1806 real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid
1807 real(r8), dimension(mgncol), intent(in) :: qiic ! cloud ice
1808 real(r8), dimension(mgncol), intent(in) :: qric ! rain
1809 real(r8), dimension(mgncol), intent(in) :: qsic ! snow
1810
1811 ! Size parameters
1812 ! rain
1813 real(r8), dimension(mgncol), intent(in) :: lamr
1814 real(r8), dimension(mgncol), intent(in) :: n0r
1815 ! snow
1816 real(r8), dimension(mgncol), intent(in) :: lams
1817 real(r8), dimension(mgncol), intent(in) :: n0s
1818
1819 ! Output tendencies
1820 real(r8), dimension(mgncol), intent(out) :: pre
1821 real(r8), dimension(mgncol), intent(out) :: prds
1822 real(r8), dimension(mgncol), intent(out) :: am_evp_st ! Fractional area where rain evaporates.
1823
1824 real(r8) :: qclr ! water vapor mixing ratio in clear air
1825 real(r8) :: ab ! correction to account for latent heat
1826 real(r8) :: eps ! 1/ sat relaxation timescale
1827 real(r8) :: tx1, tx2, tx3
1828
1829 real(r8), dimension(mgncol) :: dum
1830
1831 integer :: i
1832
1833 am_evp_st = zero
1834 ! set temporary cloud fraction to zero if cloud water + ice is very small
1835 ! this will ensure that evaporation/sublimation of precip occurs over
1836 ! entire grid cell, since min cloud fraction is specified otherwise
1837 do i=1,mgncol
1838 if (qcic(i)+qiic(i) < 1.e-6_r8) then
1839 dum(i) = zero
1840 else
1841 dum(i) = lcldm(i)
1842 end if
1843 enddo
1844 do i=1,mgncol
1845 ! only calculate if there is some precip fraction > cloud fraction
1846
1847 if (precip_frac(i) > dum(i)) then
1848
1849 if (qric(i) >= qsmall .or. qsic(i) >= qsmall) then
1850 am_evp_st(i) = precip_frac(i) - dum(i)
1851
1852 ! calculate q for out-of-cloud region
1853 qclr = (q(i)-dum(i)*qvl(i)) / (one-dum(i))
1854 end if
1855
1856 ! evaporation of rain
1857 if (qric(i) >= qsmall) then
1858
1859 ab = calc_ab(t(i), qvl(i), xxlv)
1860 eps = two*pi*n0r(i)*rho(i)*dv(i) * &
1861 (f1r/(lamr(i)*lamr(i)) + &
1862 f2r*sqrt(arn(i)*rho(i)/mu(i)) * &
1863 sc(i)**oneo3*gamma_half_br_plus5 &
1864 / (lamr(i)**((five+br)*half)))
1865
1866 pre(i) = eps*(qclr-qvl(i)) / ab
1867
1868 ! only evaporate in out-of-cloud region
1869 ! and distribute across precip_frac
1870 pre(i) = min(pre(i)*am_evp_st(i), zero)
1871 pre(i) = pre(i) / precip_frac(i)
1872 else
1873 pre(i) = zero
1874 end if
1875
1876 ! sublimation of snow
1877 if (qsic(i) >= qsmall) then
1878 ab = calc_ab(t(i), qvi(i), xxls)
1879 eps = two*pi*n0s(i)*rho(i)*dv(i) * &
1880 ( f1s/(lams(i)*lams(i)) &
1881 + f2s*sqrt(asn(i)*rho(i)/mu(i)) * &
1882 sc(i)**oneo3*gamma_half_bs_plus5 &
1883 / (lams(i)**((five+bs)*half)))
1884 prds(i) = eps*(qclr-qvi(i)) / ab
1885
1886 ! only sublimate in out-of-cloud region and distribute over precip_frac
1887 prds(i) = min(prds(i)*am_evp_st(i), zero)
1888 prds(i) = prds(i) / precip_frac(i)
1889 else
1890 prds(i) = zero
1891 end if
1892
1893 else
1894 prds(i) = zero
1895 pre(i) = zero
1896 end if
1897 enddo
1898
1899end subroutine evaporate_sublimate_precip
1900
1903!===================================================================
1904! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
1905! in-cloud condensation/deposition of rain and snow is neglected
1906! except for transfer of cloud water to snow through bergeron process
1907subroutine evaporate_sublimate_precip_graupel(t, rho, dv, mu, sc, q, qvl, qvi, &
1908 lcldm, precip_frac, arn, asn, agn, bg, qcic, qiic, qric, qsic, qgic, lamr, n0r, lams, n0s, lamg, n0g, &
1909 pre, prds, prdg, am_evp_st, mgncol)
1910
1911 integer, intent(in) :: mgncol
1912
1913 real(r8), dimension(mgncol), intent(in) :: t ! temperature
1914 real(r8), dimension(mgncol), intent(in) :: rho ! air density
1915 real(r8), dimension(mgncol), intent(in) :: dv ! water vapor diffusivity
1916 real(r8), dimension(mgncol), intent(in) :: mu ! viscosity
1917 real(r8), dimension(mgncol), intent(in) :: sc ! schmidt number
1918 real(r8), dimension(mgncol), intent(in) :: q ! humidity
1919 real(r8), dimension(mgncol), intent(in) :: qvl ! saturation humidity (water)
1920 real(r8), dimension(mgncol), intent(in) :: qvi ! saturation humidity (ice)
1921 real(r8), dimension(mgncol), intent(in) :: lcldm ! liquid cloud fraction
1922 real(r8), dimension(mgncol), intent(in) :: precip_frac ! precipitation fraction (maximum overlap)
1923
1924 ! fallspeed parameters
1925 real(r8), dimension(mgncol), intent(in) :: arn ! rain
1926 real(r8), dimension(mgncol), intent(in) :: asn ! snow
1927!++ag
1928 real(r8), dimension(mgncol), intent(in) :: agn ! graupel
1929 real(r8), intent(in) :: bg
1930!--ag
1931
1932 ! In-cloud MMRs
1933 real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid
1934 real(r8), dimension(mgncol), intent(in) :: qiic ! cloud ice
1935 real(r8), dimension(mgncol), intent(in) :: qric ! rain
1936 real(r8), dimension(mgncol), intent(in) :: qsic ! snow
1937 real(r8), dimension(mgncol), intent(in) :: qgic ! graupel
1938
1939 ! Size parameters
1940 ! rain
1941 real(r8), dimension(mgncol), intent(in) :: lamr
1942 real(r8), dimension(mgncol), intent(in) :: n0r
1943 ! snow
1944 real(r8), dimension(mgncol), intent(in) :: lams
1945 real(r8), dimension(mgncol), intent(in) :: n0s
1946!++ag
1947 ! graupel
1948 real(r8), dimension(mgncol), intent(in) :: lamg
1949 real(r8), dimension(mgncol), intent(in) :: n0g
1950!--ag
1951
1952 ! Output tendencies
1953 real(r8), dimension(mgncol), intent(out) :: pre
1954 real(r8), dimension(mgncol), intent(out) :: prds
1955!++ag
1956 real(r8), dimension(mgncol), intent(out) :: prdg
1957!--ag
1958 real(r8), dimension(mgncol), intent(out) :: am_evp_st ! Fractional area where rain evaporates.
1959
1960 real(r8) :: qclr ! water vapor mixing ratio in clear air
1961 real(r8) :: ab ! correction to account for latent heat
1962 real(r8) :: eps ! 1/ sat relaxation timescale
1963
1964 real(r8), dimension(mgncol) :: dum
1965
1966 integer :: i
1967
1968 ! set temporary cloud fraction to zero if cloud water + ice is very small
1969 ! this will ensure that evaporation/sublimation of precip occurs over
1970 ! entire grid cell, since min cloud fraction is specified otherwise
1971 do i=1,mgncol
1972 if (qcic(i)+qiic(i) < 1.e-6_r8) then
1973 dum(i) = zero
1974 else
1975 dum(i) = lcldm(i)
1976 end if
1977 enddo
1978 do i=1,mgncol
1979 ! only calculate if there is some precip fraction > cloud fraction
1980
1981 if (precip_frac(i) > dum(i)) then
1982
1983 if (qric(i) >= qsmall .or. qsic(i) >= qsmall .or. qgic(i) >= qsmall) then
1984 am_evp_st(i) = precip_frac(i) - dum(i)
1985
1986 ! calculate q for out-of-cloud region
1987 qclr = (q(i)-dum(i)*qvl(i)) / (one-dum(i))
1988 end if
1989
1990 ! evaporation of rain
1991 if (qric(i) >= qsmall) then
1992
1993 ab = calc_ab(t(i), qvl(i), xxlv)
1994 eps = twopi*n0r(i)*rho(i)*dv(i)* &
1995 ( f1r/(lamr(i)*lamr(i)) &
1996 + f2r*sqrt(arn(i)*rho(i)/mu(i)) &
1997 * sc(i)**oneo3*gamma_half_br_plus5 &
1998 / (lamr(i)**((five+br)*half)))
1999
2000 pre(i) = eps*(qclr-qvl(i))/ab
2001
2002 ! only evaporate in out-of-cloud region
2003 ! and distribute across precip_frac
2004 pre(i) = min(pre(i)*am_evp_st(i), zero)
2005 pre(i) = pre(i)/precip_frac(i)
2006 else
2007 pre(i) = zero
2008 end if
2009
2010 ! sublimation of snow
2011 if (qsic(i) >= qsmall) then
2012 ab = calc_ab(t(i), qvi(i), xxls)
2013 eps = twopi*n0s(i)*rho(i)*dv(i)* &
2014 ( f1s/(lams(i)*lams(i)) &
2015 + f2s*sqrt(asn(i)*rho(i)/mu(i)) &
2016 * sc(i)**oneo3*gamma_half_bs_plus5 &
2017 / (lams(i)**((five+bs)*half)))
2018 prds(i) = eps*(qclr-qvi(i))/ab
2019
2020 ! only sublimate in out-of-cloud region and distribute over precip_frac
2021 prds(i) = min(prds(i)*am_evp_st(i), zero)
2022 prds(i) = prds(i)/precip_frac(i)
2023 else
2024 prds(i) = zero
2025 end if
2026
2027!++AG ADD GRAUPEL, do Same with prdg.
2028
2029 if (qgic(i) >= qsmall) then
2030 ab = calc_ab(t(i), qvi(i), xxls)
2031
2032 eps = twopi*n0g(i)*rho(i)*dv(i)* &
2033 ( f1s/(lamg(i)*lamg(i)) &
2034 + f2s*sqrt(agn(i)*rho(i)/mu(i)) &
2035 * sc(i)**oneo3*gamma((five+bg)*half) &
2036 / (lamg(i)**((five+bs)*half)))
2037! / (lamg(i)**((five+bg)*half))) ! changing bs to bg - Moorthi
2038 prdg(i) = eps*(qclr-qvi(i))/ab
2039
2040 ! only sublimate in out-of-cloud region and distribute over precip_frac
2041 prdg(i) = min(prdg(i)*am_evp_st(i), zero)
2042 prdg(i) = prdg(i)/precip_frac(i)
2043 else
2044 prdg(i) = zero
2045 end if
2046
2047 else
2048 prds(i) = zero
2049 pre(i) = zero
2050!++ag
2051 prdg(i) = zero
2052!--ag
2053 end if
2054 enddo
2055
2057
2060subroutine bergeron_process_snow(t, rho, dv, mu, sc, qvl, qvi, asn, &
2061 qcic, qsic, lams, n0s, bergs, mgncol)
2062
2063 integer, intent(in) :: mgncol
2064
2065 real(r8), dimension(mgncol), intent(in) :: t ! temperature
2066 real(r8), dimension(mgncol), intent(in) :: rho ! air density
2067 real(r8), dimension(mgncol), intent(in) :: dv ! water vapor diffusivity
2068 real(r8), dimension(mgncol), intent(in) :: mu ! viscosity
2069 real(r8), dimension(mgncol), intent(in) :: sc ! schmidt number
2070 real(r8), dimension(mgncol), intent(in) :: qvl ! saturation humidity (water)
2071 real(r8), dimension(mgncol), intent(in) :: qvi ! saturation humidity (ice)
2072
2073 ! fallspeed parameter for snow
2074 real(r8), dimension(mgncol), intent(in) :: asn
2075
2076 ! In-cloud MMRs
2077 real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid mixing ratio
2078 real(r8), dimension(mgncol), intent(in) :: qsic ! snow mixing ratio
2079
2080 ! Size parameters for snow
2081 real(r8), dimension(mgncol), intent(in) :: lams
2082 real(r8), dimension(mgncol), intent(in) :: n0s
2083
2084 ! Output tendencies
2085 real(r8), dimension(mgncol), intent(out) :: bergs
2086
2087 real(r8) :: ab ! correction to account for latent heat
2088 real(r8) :: eps ! 1/ sat relaxation timescale
2089
2090 integer :: i
2091
2092 do i=1,mgncol
2093 if (qsic(i) >= qsmall.and. qcic(i) >= qsmall .and. t(i) < tmelt) then
2094 ab = calc_ab(t(i), qvi(i), xxls)
2095 eps = two*pi*n0s(i)*rho(i)*dv(i) * &
2096 (f1s/(lams(i)*lams(i)) + &
2097 f2s*sqrt(asn(i)*rho(i)/mu(i)) * &
2098 sc(i)**oneo3*gamma_half_bs_plus5 / &
2099 (lams(i)**((five+bs)*half)))
2100 bergs(i) = eps*(qvl(i)-qvi(i)) / ab
2101 else
2102 bergs(i) = zero
2103 end if
2104 enddo
2105end subroutine bergeron_process_snow
2106
2107!========================================================================
2110subroutine graupel_collecting_snow(qsic,qric,umr,ums,rho,lamr,n0r,lams,n0s, &
2111 psacr, mgncol)
2112
2113 integer, intent(in) :: mgncol
2114
2115 ! In-cloud MMRs
2116 real(r8), dimension(mgncol), intent(in) :: qsic ! snow
2117 real(r8), dimension(mgncol), intent(in) :: qric ! rain
2118
2119 ! mass-weighted fall speeds
2120 real(r8), dimension(mgncol), intent(in) :: umr ! rain
2121 real(r8), dimension(mgncol), intent(in) :: ums ! snow
2122
2123 real(r8), dimension(mgncol), intent(in) :: rho ! air density
2124
2125
2126 ! Size parameters for rain
2127 real(r8), dimension(mgncol), intent(in) :: lamr
2128 real(r8), dimension(mgncol), intent(in) :: n0r
2129
2130 ! Size parameters for snow
2131 real(r8), dimension(mgncol), intent(in) :: lams
2132 real(r8), dimension(mgncol), intent(in) :: n0s
2133
2134 real(r8), dimension(mgncol), intent(out) :: psacr ! conversion due to coll of snow by rain
2135
2136 real(r8) :: cons31, tx1, tx2, tx3, tx4, tx5
2137 integer :: i
2138
2139 cons31 = pi*pi*ecr*rhosn
2140
2141 do i=1,mgncol
2142
2143 if (qsic(i) >= 0.1e-3_r8 .and. qric(i) >= 0.1e-3_r8) then
2144 tx1 = 1.2_r8*umr(i) - 0.95_r8*ums(i)
2145 tx1 = sqrt(tx1*tx1+0.08_r8*ums(i)*umr(i))
2146 tx2 = one / lams(i)
2147 tx3 = one / lamr(i)
2148 tx4 = tx2 * tx2
2149 tx5 = tx4 * tx4 * tx3
2150
2151 psacr(i) = cons31 * tx1 * rho(i) * n0r(i) * n0s(i) * tx5 &
2152 * (5.0_r8*tx4+tx3*(tx2+tx2+0.5_r8*tx3))
2153
2154! psacr(i) = cons31*(((1.2_r8*umr(i)-0.95_r8*ums(i))**2+ &
2155! 0.08_r8*ums(i)*umr(i))**0.5_r8*rho(i)* &
2156! n0r(i)*n0s(i)/lams(i)**3* &
2157! (5._r8/(lams(i)**3*lamr(i))+ &
2158! 2._r8/(lams(i)**2*lamr(i)**2)+ &
2159! 0.5_r8/(lams(i)*lamr(i)**3)))
2160 else
2161 psacr(i) = zero
2162 end if
2163
2164 end do
2165
2166end subroutine graupel_collecting_snow
2167
2168!========================================================================
2171subroutine graupel_collecting_cld_water(qgic,qcic,ncic,rho,n0g,lamg,bg,agn, &
2172 psacwg, npsacwg, mgncol)
2173
2174 integer, intent(in) :: mgncol
2175
2176 ! In-cloud MMRs
2177 real(r8), dimension(mgncol), intent(in) :: qgic ! graupel
2178 real(r8), dimension(mgncol), intent(in) :: qcic ! cloud water
2179
2180 real(r8), dimension(mgncol), intent(in) :: ncic ! cloud water number conc
2181
2182 real(r8), dimension(mgncol), intent(in) :: rho ! air density
2183
2184 ! Size parameters for graupel
2185 real(r8), dimension(mgncol), intent(in) :: lamg
2186 real(r8), dimension(mgncol), intent(in) :: n0g
2187
2188 ! fallspeed parameters for graupel
2189 ! Set AGN and BG as input (in micro_mg3_0.F90) (select hail or graupel as appropriate)
2190 real(r8), intent(in) :: bg
2191 real(r8), dimension(mgncol), intent(in) :: agn
2192
2193 ! Output tendencies
2194 real(r8), dimension(mgncol), intent(out) :: psacwg
2195 real(r8), dimension(mgncol), intent(out) :: npsacwg
2196
2197 real(r8) :: cons, tx1
2198 integer :: i
2199
2200 cons = gamma(bg+three) * pi/four * ecid
2201
2202 do i=1,mgncol
2203
2204 if (qgic(i) >= 1.e-8_r8 .and. qcic(i) >= qsmall) then
2205
2206 tx1 = cons*agn(i)*rho(i)*n0g(i) / lamg(i)**(bg+three)
2207
2208 psacwg(i) = tx1 * qcic(i)
2209 npsacwg(i) = tx1 * ncic(i)
2210 else
2211 psacwg(i) = zero
2212 npsacwg(i) = zero
2213 end if
2214 enddo
2215end subroutine graupel_collecting_cld_water
2216
2217!========================================================================
2220subroutine graupel_riming_liquid_snow(psacws,qsic,qcic,nsic,rho,rhosn,rhog,asn,lams,n0s,dtime, &
2221 pgsacw,nscng,mgncol)
2222
2223 integer, intent(in) :: mgncol
2224
2225 ! Accretion of cloud water to snow tendency (modified)
2226 real(r8), dimension(mgncol), intent(inout) :: psacws
2227
2228 real(r8), dimension(mgncol), intent(in) :: qsic ! snow mixing ratio
2229 real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid mixing ratio
2230 real(r8), dimension(mgncol), intent(in) :: nsic ! snow number concentration
2231
2232 real(r8), dimension(mgncol), intent(in) :: rho ! air density
2233 real(r8), intent(in) :: rhosn ! snow density
2234 real(r8), intent(in) :: rhog ! graupel density
2235
2236 real(r8), dimension(mgncol), intent(in) :: asn ! fall speed parameter for snow
2237
2238 ! Size parameters for snow
2239 real(r8), dimension(mgncol), intent(in) :: lams
2240 real(r8), dimension(mgncol), intent(in) :: n0s
2241
2242 real(r8), intent(in) :: dtime
2243
2244 !Output process rates
2245 real(r8), dimension(mgncol), intent(out) :: pgsacw ! dQ graupel due to collection droplets by snow
2246 real(r8), dimension(mgncol), intent(out) :: nscng ! dN graupel due to collection droplets by snow
2247
2248 real(r8) :: cons
2249 real(r8) :: rhosu
2250 real(r8) :: dum
2251 integer :: i
2252
2253!........................................................................
2254!Input: PSACWS,qs,qc,n0s,rho,lams,rhos,rhog
2255!Output:PSACWS,PGSACW,NSCNG
2256
2257 rhosu = 85000._r8/(ra * tmelt) ! typical air density at 850 mb
2258
2259 do i=1,mgncol
2260
2261 cons=4._r8 *2._r8 *3._r8 *rhosu*pi*ecid*ecid*gamma_2bs_plus2/(8._r8*(rhog-rhosn))
2262
2263 if (psacws(i).gt.0._r8 .and. qsic(i).GE.0.1e-3_r8 .AND. qcic(i).GE.0.5e-3_r8) then
2264! Only allow conversion if qni > 0.1 and qc > 0.5 g/kg following Rutledge and Hobbs (1984)
2265 !if (qsic(i).GE.0.1e-3_r8 .AND. qcic(i).GE.0.5E-3_r8) then
2266
2267! portion of riming converted to graupel (Reisner et al. 1998, originally IS1991)
2268! dtime here is correct.
2269 pgsacw(i) = min(psacws(i), cons*dtime*n0s(i)*qcic(i)*qcic(i)* &
2270 asn(i)*asn(i)/ (rho(i)*lams(i)**(bs+bs+two)))
2271
2272! if (pgsacw(i).lt.0_r8) then
2273! write(iulog,*) "pgsacw,i,lams,cons",i,pgsacw(i),lams(i),cons
2274! end if
2275
2276! Mix rat converted into graupel as embryo (Reisner et al. 1998, orig M1990)
2277 dum= max(rhosn/(rhog-rhosn)*pgsacw(i), zero)
2278
2279! Number concentraiton of embryo graupel from riming of snow
2280 nscng(i) = dum/mg0*rho(i)
2281! Limit max number converted to snow number (dtime here correct? We think yes)
2282 nscng(i) = min(nscng(i),nsic(i)/dtime)
2283
2284! Portion of riming left for snow
2285 psacws(i) = psacws(i) - pgsacw(i)
2286 else
2287 pgsacw(i) = zero
2288 nscng(i) = zero
2289 end if
2290
2291 enddo
2292
2293end subroutine graupel_riming_liquid_snow
2294
2295!========================================================================
2298subroutine graupel_collecting_rain(qric,qgic,umg,umr,ung,unr,rho,n0r,lamr,n0g,lamg,&
2299 pracg,npracg,mgncol)
2300
2301 integer, intent(in) :: mgncol
2302
2303 !MMR
2304 real(r8), dimension(mgncol), intent(in) :: qric !rain mixing ratio
2305 real(r8), dimension(mgncol), intent(in) :: qgic !graupel mixing ratio
2306
2307 !Mass weighted Fall speeds
2308 real(r8), dimension(mgncol), intent(in) :: umg ! graupel fall speed
2309 real(r8), dimension(mgncol), intent(in) :: umr ! rain fall speed
2310
2311 !Number weighted fall speeds
2312 real(r8), dimension(mgncol), intent(in) :: ung ! graupel fall speed
2313 real(r8), dimension(mgncol), intent(in) :: unr ! rain fall speed
2314
2315 real(r8), dimension(mgncol), intent(in) :: rho ! air density
2316
2317 ! Size parameters for rain
2318 real(r8), dimension(mgncol), intent(in) :: n0r
2319 real(r8), dimension(mgncol), intent(in) :: lamr
2320
2321 ! Size parameters for graupel
2322 real(r8), dimension(mgncol), intent(in) :: n0g
2323 real(r8), dimension(mgncol), intent(in) :: lamg
2324
2325
2326 !Output process rates
2327 real(r8), dimension(mgncol), intent(out) :: pracg ! Q collection rain by graupel
2328 real(r8), dimension(mgncol), intent(out) :: npracg ! N collection rain by graupel
2329
2330! Add collection of graupel by rain above freezing
2331! assume all rain collection by graupel above freezing is shed
2332! assume shed drops are 1 mm in size
2333
2334 integer :: i
2335 real(r8) :: cons41
2336 real(r8) :: cons32
2337 real(r8) :: dum, tx1, tx2, tx3, tx4, tx5, tx6
2338
2339 cons41 = pi*pi*ecr*rhow
2340 cons32 = 0.5*pi*ecr
2341
2342 do i=1,mgncol
2343
2344 if (qric(i) >= 1.e-8_r8 .and. qgic(i) >= 1.e-8_r8) then
2345
2346! pracg is mixing ratio of rain per sec collected by graupel/hail
2347 tx1 = 1.2_r8*umr(i) - 0.95_r8*umg(i)
2348 tx1 = sqrt(tx1*tx1+0.08_r8*umg(i)*umr(i))
2349 tx2 = 1.0_r8 / lamr(i)
2350 tx3 = 1.0_r8 / lamg(i)
2351 tx4 = tx2 * tx2
2352 tx5 = tx4 * tx4 * tx3
2353 tx6 = rho(i) * n0r(i) * n0g(i)
2354
2355
2356 pracg(i) = cons41 * tx1 * tx6 * tx5 * (5.0*tx4+tx3*(tx2+tx2+0.5*tx3))
2357
2358
2359! pracg(i) = cons41*(((1.2_r8*umr(i)-0.95_r8*umg(i))**2._r8+ &
2360! 0.08_r8*umg(i)*umr(i))**0.5_r8*rho(i)* &
2361! n0r(i)*n0g(i)/lamr(i)**3._r8* &
2362! (5._r8/(lamr(i)**3._r8*lamg(i))+ &
2363! 2._r8/(lamr(i)**2._r8*lamg(i)**2._r8)+ &
2364! 0.5_r8/(lamr(i)*lamg(i)**3._r8)))
2365
2366! assume 1 mm drops are shed, get number shed per sec
2367
2368 dum = pracg(i) / 5.2e-7_r8
2369
2370 tx1 = unr(i) - ung(i)
2371 tx1 = sqrt(1.7_r8 * tx1 * tx1 + 0.3_r8*unr(i)*ung(i))
2372 tx4 = tx2 * tx3
2373
2374 npracg(i) = cons32 * tx1 * tx6 * tx4 * (tx2*(tx2+tx3)+tx3*tx3)
2375
2376! npracg(i) = cons32*rho(i)*(1.7_r8*(unr(i)-ung(i))**2._r8+ &
2377! 0.3_r8*unr(i)*ung(i))**0.5_r8*n0r(i)*n0g(i)* &
2378! (1._r8/(lamr(i)**3._r8*lamg(i))+ &
2379! 1._r8/(lamr(i)**2._r8*lamg(i)**2._r8)+ &
2380! 1._r8/(lamr(i)*lamg(i)**3._r8))
2381
2382! hm 7/15/13, remove limit so that the number of collected drops can smaller than
2383! number of shed drops
2384! NPRACG(K)=MAX(NPRACG(K)-DUM,0.)
2385 npracg(i) = npracg(i) - dum
2386 else
2387 npracg(i) = zero
2388 pracg(i) = zero
2389 end if
2390
2391 enddo
2392
2393end subroutine graupel_collecting_rain
2394
2395!========================================================================
2398!========================================================================
2399! Conversion of rimed rainwater onto snow converted to graupel
2400subroutine graupel_rain_riming_snow(pracs,npracs,psacr,qsic,qric,nric,nsic,n0s, &
2401 lams,n0r,lamr,dtime,pgracs,ngracs,mgncol)
2402
2403 integer, intent(in) :: mgncol
2404
2405 ! Accretion of rain by snow
2406 real(r8), dimension(mgncol), intent(inout) :: pracs
2407 real(r8), dimension(mgncol), intent(inout) :: npracs
2408 real(r8), dimension(mgncol), intent(inout) :: psacr ! conversion due to coll of snow by rain
2409
2410
2411 real(r8), dimension(mgncol), intent(in) :: qsic !snow mixing ratio
2412 real(r8), dimension(mgncol), intent(in) :: qric !rain mixing ratio
2413
2414 real(r8), dimension(mgncol), intent(in) :: nric ! rain number concentration
2415 real(r8), dimension(mgncol), intent(in) :: nsic ! snow number concentration
2416
2417 ! Size parameters for snow
2418 real(r8), dimension(mgncol), intent(in) :: n0s
2419 real(r8), dimension(mgncol), intent(in) :: lams
2420
2421 ! Size parameters for rain
2422 real(r8), dimension(mgncol), intent(in) :: n0r
2423 real(r8), dimension(mgncol), intent(in) :: lamr
2424
2425 real(r8), intent(in) :: dtime
2426
2427 !Output process rates
2428 real(r8), dimension(mgncol), intent(out) :: pgracs ! Q graupel due to collection rain by snow
2429 real(r8), dimension(mgncol), intent(out) :: ngracs ! N graupel due to collection rain by snow
2430
2431!Input: PRACS,NPRACS,PSACR,qni,qr,lams,lamr,nr,ns Note: No PSACR in MG2
2432!Output:PGRACS,NGRACS,PRACS,PSACR
2433
2434 integer :: i
2435 real(r8) :: cons18
2436 real(r8) :: cons19
2437 real(r8) :: dum,fmult,tx1,tx2
2438
2439 cons18 = rhosn*rhosn
2440 cons19 = rhow*rhow
2441
2442 do i=1,mgncol
2443
2444 fmult = zero
2445
2446 if (pracs(i) > zero .and. qsic(i) >= 0.1e-3_r8 .and. qric(i) >= 0.1e-3_r8) then
2447 ! only allow conversion if qs > 0.1 and qr > 0.1 g/kg following rutledge and hobbs (1984)
2448 !if (qsic(i).ge.0.1e-3_r8.and.qric(i).ge.0.1e-3_r8) then
2449 ! portion of collected rainwater converted to graupel (reisner et al. 1998)
2450 tx1 = four / lams(i)
2451 tx2 = four / lamr(i)
2452 tx1 = tx1 * tx1 * tx1
2453 tx2 = tx2 * tx2 * tx2
2454 dum = cons18 * tx1 * tx1
2455 dum = one - max(zero, min(one, dum / (dum + cons19 * tx2 * tx2)))
2456
2457! dum = cons18*(4._r8/lams(i))**3*(4._r8/lams(i))**3 &
2458! /(cons18*(4._r8/lams(i))**3*(4._r8/lams(i))**3+ &
2459! cons19*(4._r8/lamr(i))**3*(4._r8/lamr(i))**3)
2460! dum = min(dum,one)
2461! dum = max(dum, zero)
2462!
2463! pgracs(i) = (one-dum) * pracs(i)
2464! ngracs(i) = (one-dum) * npracs(i)
2465
2466 pgracs(i) = dum * pracs(i)
2467 ngracs(i) = dum * npracs(i)
2468 ! limit max number converted to min of either rain or snow number concentration
2469 ngracs(i) = min(ngracs(i),nric(i)/dtime)
2470 ngracs(i) = min(ngracs(i),nsic(i)/dtime)
2471
2472 ! amount left for snow production
2473 pracs(i) = pracs(i) - pgracs(i)
2474 npracs(i) = npracs(i) - ngracs(i)
2475
2476 ! conversion to graupel due to collection of snow by rain
2477! psacr(i) = psacr(i) * (one-dum)
2478 psacr(i) = psacr(i) * dum
2479 else
2480 pgracs(i) = zero
2481 ngracs(i) = zero
2482 end if
2483 enddo
2484
2485end subroutine graupel_rain_riming_snow
2486
2487!========================================================================
2488! Rime Splintering
2489!========================================================================
2492subroutine graupel_rime_splintering(t,qcic,qric,qgic,psacwg,pracg,&
2493 qmultg,nmultg,qmultrg,nmultrg,mgncol)
2494
2495 integer, intent(in) :: mgncol
2496
2497 real(r8), dimension(mgncol), intent(in) :: t !temperature
2498
2499 !MMR
2500 real(r8), dimension(mgncol), intent(in) :: qcic !liquid mixing ratio
2501 real(r8), dimension(mgncol), intent(in) :: qric !rain mixing ratio
2502 real(r8), dimension(mgncol), intent(in) :: qgic !graupel mixing ratio
2503
2504 ! Already calculated terms for collection
2505 real(r8), dimension(mgncol), intent(inout) :: psacwg ! collection droplets by graupel
2506 real(r8), dimension(mgncol), intent(inout) :: pracg ! collection rain by graupel
2507
2508 !Output process rates for splintering
2509 real(r8), dimension(mgncol), intent(out) :: qmultg ! Q ice mult of droplets/graupel
2510 real(r8), dimension(mgncol), intent(out) :: nmultg ! N ice mult of droplets/graupel
2511 real(r8), dimension(mgncol), intent(out) :: qmultrg ! Q ice mult of rain/graupel
2512 real(r8), dimension(mgncol), intent(out) :: nmultrg ! N ice mult of rain/graupel
2513
2514
2515!Input: qg,qc,qr, PSACWG,PRACG,T
2516!Output: NMULTG,QMULTG,NMULTRG,QMULTRG,PSACWG,PRACG
2517
2518 integer :: i
2519 real(r8) :: fmult
2520 real(r8) :: tm_3,tm_5,tm_8
2521
2522 tm_3 = tmelt - 3._r8
2523 tm_5 = tmelt - 5._r8
2524 tm_8 = tmelt - 8._r8
2525
2526
2527!nmultg,qmultg .
2528!========================================================================
2529 do i=1,mgncol
2530
2531 nmultrg(i) = zero
2532 qmultrg(i) = zero
2533 nmultg(i) = zero
2534 qmultg(i) = zero
2535
2536 if (qgic(i) >= 0.1e-3_r8) then
2537 if (qcic(i) >= 0.5e-3_r8 .or. qric(i) >= 0.1e-3_r8) then
2538 if (psacwg(i) > zero .or. pracg(i) > zero) then
2539 if (t(i) < tm_3 .and. t(i) > tm_8) then
2540 if (t(i) > tm_3) then
2541 fmult = zero
2542 else if (t(i) <= tm_3 .and. t(i) > tm_5) then
2543 fmult = (tm_3-t(i)) * 0.5
2544 else if (t(i) >= tm_8 .and. t(i) <= tm_5) then
2545 fmult = (t(i)-tm_8) * (one/three)
2546 else if (t(i) < tm_8) then
2547 fmult = zero
2548 end if
2549
2550! 1000 is to convert from kg to g (Check if needed for MG)
2551
2552! splintering from droplets accreted onto graupel
2553
2554 if (psacwg(i) > zero) then
2555 nmultg(i) = 35.e4_r8*psacwg(i)*fmult*1000._r8
2556 qmultg(i) = nmultg(i)*mmult
2557
2558! constrain so that transfer of mass from graupel to ice cannot be more mass
2559! than was rimed onto graupel
2560
2561 qmultg(i) = min(qmultg(i),psacwg(i))
2562 psacwg(i) = psacwg(i) - qmultg(i)
2563
2564 end if
2565
2566
2567!nmultrg,qmultrg .
2568!========================================================================
2569
2570! riming and splintering from accreted raindrops
2571
2572! Factor of 1000. again (Check)
2573
2574 if (pracg(i) > zero) then
2575 nmultrg(i) = 35.e4*pracg(i)*fmult*1000._r8
2576 qmultrg(i) = nmultrg(i)*mmult
2577
2578! constrain so that transfer of mass from graupel to ice cannot be more mass
2579! than was rimed onto graupel
2580
2581 qmultrg(i) = min(qmultrg(i),pracg(i))
2582 pracg(i) = pracg(i) - qmultrg(i)
2583
2584 end if
2585
2586 end if
2587 end if
2588 end if
2589 end if
2590 enddo
2591
2592end subroutine graupel_rime_splintering
2593
2594!========================================================================
2595! Evaporation and sublimation of graupel
2596!========================================================================
2597
2598!MERGE WITH RAIN AND SNOW EVAP
2599!
2600!subroutine graupel_sublimate_evap(t,q,qgic,rho,n0g,lamg,qvi,dv,mu,sc,bg,agn,&
2601! prdg,eprdg,mgncol)
2602!
2603! integer, intent(in) :: mgncol
2604!
2605! real(r8), dimension(mgncol), intent(in) :: t !temperature
2606! real(r8), dimension(mgncol), intent(in) :: q !specific humidity (mixing ratio)
2607!
2608! !MMR
2609! real(r8), dimension(mgncol), intent(in) :: qgic !graupel mixing ratio
2610!
2611! real(r8), dimension(mgncol), intent(in) :: rho ! air density
2612!
2613! ! Size parameters for graupel
2614! real(r8), dimension(mgncol), intent(in) :: n0g
2615! real(r8), dimension(mgncol), intent(in) :: lamg
2616!
2617! real(r8), dimension(mgncol), intent(in) :: qvi !saturation humidity (ice)
2618!
2619! real(r8), dimension(mgncol), intent(in) :: dv ! water vapor diffusivity
2620! real(r8), dimension(mgncol), intent(in) :: mu ! viscosity
2621! real(r8), dimension(mgncol), intent(in) :: sc ! schmidt number
2622!
2623! ! fallspeed parameters for graupel
2624! ! Set AGN and BG as input (in micro_mg3_0.F90) (select hail or graupel as appropriate)
2625! real(r8), intent(in) :: bg
2626! real(r8), dimension(mgncol), intent(in) :: agn
2627!
2628! ! Output tendencies (sublimation or evaporation of graupel)
2629! real(r8), dimension(mgncol), intent(out) :: prdg
2630! real(r8), dimension(mgncol), intent(out) :: eprdg
2631!
2632! real(r8) :: cons11,cons36
2633! real(r8) :: abi
2634! real(r8) :: epsg
2635! integer :: i
2636!
2637! cons11=gamma(2.5_r8+bg/2._r8) !bg will be different for graupel (bg) or hail (bh)
2638! cons36=(2.5_r8+bg/2._r8)
2639!
2640!
2641! do i=1,mgncol
2642!
2643! abi = calc_ab(t(i), qvi(i), xxls)
2644!
2645! if (qgic(i).ge.qsmall) then
2646! epsg = 2._r8*pi*n0g(i)*rho(i)*dv(i)* &
2647! (f1s/(lamg(i)*lamg(i))+ &
2648! f2s*(agn(i)*rho(i)/mu(i))**0.5_r8* &
2649! sc(i)**(1._r8/3._r8)*cons11/ &
2650! (lamg(i)**cons36))
2651! else
2652! epsg = 0.
2653! end if
2654!
2655!! vapor dpeosition on graupel
2656! prdg(i) = epsg*(q(i)-qvi(i))/abi
2657!
2658!! make sure not pushed into ice supersat/subsat
2659!! put this in main mg3 code ... check for it ...
2660!! formula from reisner 2 scheme
2661
2662!!
2663!! dum = (qv3d(k)-qvi(k))/dt
2664!!
2665!! fudgef = 0.9999
2666!! sum_dep = prd(k)+prds(k)+mnuccd(k)+prdg(k)
2667!!
2668!! if( (dum.gt.0. .and. sum_dep.gt.dum*fudgef) .or. &
2669!! (dum.lt.0. .and. sum_dep.lt.dum*fudgef) ) then
2670!! prdg(k) = fudgef*prdg(k)*dum/sum_dep
2671!! endif
2672!
2673!! if cloud ice/snow/graupel vap deposition is neg, then assign to sublimation processes
2674!
2675! eprdg(i)=0._r8
2676!
2677! if (prdg(i).lt.0._r8) then
2678! eprdg(i)=prdg(i)
2679! prdg(i)=0.
2680! end if
2681!
2682! enddo
2683!
2684!end subroutine graupel_sublimate_evap
2685
2686!========================================================================
2687!UTILITIES
2688!========================================================================
2689
2691pure function no_limiter()
2692 real(r8) :: no_limiter
2693
2694 no_limiter = transfer(limiter_off, no_limiter)
2695
2696end function no_limiter
2697
2699pure function limiter_is_on(lim)
2700 real(r8), intent(in) :: lim
2701 logical :: limiter_is_on
2702
2703 limiter_is_on = transfer(lim, limiter_off) /= limiter_off
2704
2705end function limiter_is_on
2706
2708FUNCTION gamma_incomp(muice, x)
2709
2710 real(r8) :: gamma_incomp
2711 REAL(r8), intent(in) :: muice, x
2712 REAL(r8) :: xog, kg, alfa, auxx
2713 alfa = min(max(muice+1._r8, 1._r8), 20._r8)
2714
2715 xog = log(alfa -0.3068_r8)
2716 kg = 1.44818_r8*(alfa**0.5357_r8)
2717 auxx = max(min(kg*(log(x)-xog), 30._r8), -30._r8)
2718 gamma_incomp = max(one/(one+exp(-auxx)), 1.0e-20)
2719
2720END FUNCTION gamma_incomp
2721
2722
2723end module micro_mg_utils
elemental real(r8) function var_coef_r8(relvar, a)
Finds a coefficient for process rates based on the relative variance of cloud water.
elemental real(r8) function var_coef_integer(relvar, a)
Finds a coefficient for process rates based on the relative variance of cloud water.
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 size_dist_param_liq_vect(props, qcic, ncic, rho, pgam, lamc, mgncol)
This subroutine gets cloud droplet size distribution parameters.
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.
pure logical function limiter_is_on(lim)
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 size_dist_param_basic_vect(props, qic, nic, lam, mgncol, n0)
This subroutine calculates.
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
pure real(r8) function no_limiter()
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
subroutine size_dist_param_ice_vect(props, qic, nic, lam, mgncol, n0)
This subroutine.
real(r8) function gamma_incomp(muice, x)
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.
elemental subroutine size_dist_param_ice_line(props, qic, nic, lam, n0)
ice routine for getting size distribution parameters.
elemental subroutine size_dist_param_basic_line(props, qic, nic, lam, n0)
Basic routine for getting size distribution parameters.
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)
elemental subroutine size_dist_param_liq_line(props, qcic, ncic, rho, pgam, lamc)
get cloud droplet size distribution parameters
type(mghydrometeorprops) function newmghydrometeorprops(rho, eff_dim, lambda_bounds, min_mean_mass)
Constructor for a constituent property object.