!>\file micro_mg_utils.F90 !! This file contains process rates and utility functions used by the !! MG microphysics. !>\ingroup mg2mg3 !>\defgroup micro_mg_utils_mod Morrison-Gettelman MP utils Module !! This module contains process rates and utility functions used by the MG !! microphysics. !! !! Original MG authors: Andrew Gettelman, Hugh Morrison !! Contributions from: Peter Caldwell, Xiaohong Liu and Steve Ghan !! !! Separated from MG 1.5 by B. Eaton. !! !! Separated module switched to MG 2.0 and further changes by S. Santos. !! !! Anning Cheng changed for FV3GFS 9/29/2017 !! added ac_time as an input !! !! S. Moorthi - Feb 2018 : code optimization !! !! This version: https://svn-ccsm-models.cgd.ucar.edu/cam1/branch_tags/mg3_tags/mg3_33_cam5_4_153/ !! !! for questions contact Hugh Morrison, Andrew Gettelman !! e-mail: morrison@ucar.edu, andrew@ucar.edu module micro_mg_utils !-------------------------------------------------------------------------- ! ! List of required external functions that must be supplied: ! gamma --> standard mathematical gamma function (if gamma is an ! intrinsic, define HAVE_GAMMA_INTRINSICS) ! !-------------------------------------------------------------------------- ! ! Constants that must be specified in the "init" method (module variables): ! ! kind kind of reals (to verify correct linkage only) - ! gravit acceleration due to gravity m s-2 ! rair dry air gas constant for air J kg-1 K-1 ! rh2o gas constant for water vapor J kg-1 K-1 ! cpair specific heat at constant pressure for dry air J kg-1 K-1 ! tmelt temperature of melting point for water K ! latvap latent heat of vaporization J kg-1 ! latice latent heat of fusion J kg-1 ! !-------------------------------------------------------------------------- ! 8 byte real and integer use machine, only : r8 => kind_phys use machine, only : i8 => kind_phys implicit none private save public :: & micro_mg_utils_init, & size_dist_param_liq, & size_dist_param_basic, & avg_diameter, & rising_factorial, & ice_deposition_sublimation, & sb2001v2_liq_autoconversion, & sb2001v2_accre_cld_water_rain, & kk2000_liq_autoconversion, & ice_autoconversion, & immersion_freezing, & contact_freezing, & snow_self_aggregation, & accrete_cloud_water_snow, & secondary_ice_production, & accrete_rain_snow, & heterogeneous_rain_freezing, & accrete_cloud_water_rain, & self_collection_rain, & accrete_cloud_ice_snow, & evaporate_sublimate_precip, & bergeron_process_snow, & liu_liq_autoconversion, & gmao_ice_autoconversion, & size_dist_param_ice, & !++ag graupel_collecting_snow, & graupel_collecting_rain, & graupel_collecting_cld_water, & graupel_riming_liquid_snow, & graupel_rain_riming_snow, & graupel_rime_splintering, & evaporate_sublimate_precip_graupel ! graupel_sublimate_evap !--ag public :: MGHydrometeorProps type :: MGHydrometeorProps ! Density (kg/m^3) real(r8) :: rho ! Information for size calculations. ! Basic calculation of mean size is: ! lambda = (shape_coef*nic/qic)^(1/eff_dim) ! Then lambda is constrained by bounds. real(r8) :: eff_dim real(r8) :: shape_coef real(r8) :: lambda_bounds(2) ! Minimum average particle mass (kg). ! Limit is applied at the beginning of the size distribution calculations. real(r8) :: min_mean_mass end type MGHydrometeorProps interface MGHydrometeorProps module procedure NewMGHydrometeorProps end interface type(MGHydrometeorProps), public :: mg_liq_props type(MGHydrometeorProps), public :: mg_ice_props type(MGHydrometeorProps), public :: mg_rain_props type(MGHydrometeorProps), public :: mg_snow_props !++ag type(MGHydrometeorProps), public :: mg_graupel_props !--ag interface size_dist_param_liq module procedure size_dist_param_liq_vect module procedure size_dist_param_liq_line end interface interface size_dist_param_basic module procedure size_dist_param_basic_vect module procedure size_dist_param_basic_line end interface interface size_dist_param_ice module procedure size_dist_param_ice_vect module procedure size_dist_param_ice_line end interface !================================================= ! Public module parameters (mostly for MG itself) !================================================= !> Pi to 20 digits; more than enough to reach the limit of double precision. real(r8), parameter, public :: pi = 3.14159265358979323846_r8 !> "One minus small number": number near unity for round-off issues. !real(r8), parameter, public :: omsm = 1._r8 - 1.e-5_r8 real(r8), parameter, public :: omsm = 1._r8 - 1.e-6_r8 !> Smallest mixing ratio considered in microphysics. real(r8), parameter, public :: qsmall = 1.e-18_r8 !> minimum allowed cloud fraction real(r8), parameter, public :: mincld = 0.000001_r8 !real(r8), parameter, public :: mincld = 0.0001_r8 !real(r8), parameter, public :: mincld = 0.0_r8 real(r8), parameter, public :: rhosn = 250._r8 !< bulk density snow real(r8), parameter, public :: rhoi = 500._r8 !< bulk density ice real(r8), parameter, public :: rhow = 1000._r8 !< bulk density liquid real(r8), parameter, public :: rhows = 917._r8 !< bulk density water solid !++ag !Hail and Graupel (set in MG3) real(r8), parameter, public :: rhog = 500._r8 real(r8), parameter, public :: rhoh = 400._r8 !--ag ! fall speed parameters, V = aD^b (V is in m/s) ! droplets real(r8), parameter, public :: ac = 3.e7_r8 real(r8), parameter, public :: bc = 2._r8 ! snow real(r8), parameter, public :: as = 11.72_r8 real(r8), parameter, public :: bs = 0.41_r8 ! cloud ice real(r8), parameter, public :: ai = 700._r8 real(r8), parameter, public :: bi = 1._r8 ! small cloud ice (r< 10 um) - sphere, bulk density real(r8), parameter, public :: aj = ac*((rhoi/rhows)**(bc/3._r8))*rhows/rhow real(r8), parameter, public :: bj = bc ! rain real(r8), parameter, public :: ar = 841.99667_r8 real(r8), parameter, public :: br = 0.8_r8 !++ag ! graupel real(r8), parameter, public :: ag = 19.3_r8 real(r8), parameter, public :: bg = 0.37_r8 ! hail real(r8), parameter, public :: ah = 114.5_r8 real(r8), parameter, public :: bh = 0.5_r8 !--ag !> mass of new crystal due to aerosol freezing and growth (kg) !! Make this consistent with the lower bound, to support UTLS and !! stratospheric ice, and the smaller ice size limit. real(r8), parameter, public :: mi0 = 4._r8/3._r8*pi*rhoi*(1.e-6_r8)**3 !++ag ! mass of new graupel particle (assume same as mi0 for now, may want to make bigger?) !real(r8), parameter, public :: mg0 = 4._r8/3._r8*pi*rhoi*(1.e-6_r8)**3 !or set based on M2005: real(r8), parameter, public :: mg0 = 1.6e-10_r8 ! radius of contact nuclei real(r8), parameter, public :: mmult = 4._r8/3._r8*pi*rhoi*(5.e-6_r8)**3 !--ag !================================================= ! Private module parameters !================================================= ! Signaling NaN bit pattern that represents a limiter that's turned off. integer(i8), parameter :: limiter_off = int(Z'7FF1111111111111', i8) ! alternate threshold used for some in-cloud mmr real(r8), parameter :: icsmall = 1.e-8_r8 ! particle mass-diameter relationship ! currently we assume spherical particles for cloud ice/snow ! m = cD^d ! exponent real(r8), parameter :: dsph = 3._r8 ! Bounds for mean diameter for different constituents. real(r8), parameter :: lam_bnd_rain(2) = 1._r8/[500.e-6_r8, 20.e-6_r8] real(r8), parameter :: lam_bnd_snow(2) = 1._r8/[2000.e-6_r8, 10.e-6_r8] ! Minimum average mass of particles. real(r8), parameter :: min_mean_mass_liq = 1.e-20_r8 real(r8), parameter :: min_mean_mass_ice = 1.e-20_r8 ! ventilation parameters ! for snow real(r8), parameter :: f1s = 0.86_r8 real(r8), parameter :: f2s = 0.28_r8 ! for rain real(r8), parameter :: f1r = 0.78_r8 real(r8), parameter :: f2r = 0.308_r8 ! collection efficiencies ! aggregation of cloud ice and snow !real(r8), parameter :: eii = 0.5_r8 !real(r8), parameter :: eii = 0.1_r8 real(r8), parameter :: eii = 0.2_r8 !++ag ! collection efficiency, ice-droplet collisions real(r8), parameter, public :: ecid = 0.7_r8 ! collection efficiency between droplets/rain and snow/rain real(r8), parameter, public :: ecr = 1.0_r8 !--ag ! immersion freezing parameters, bigg 1953 real(r8), parameter :: bimm = 100._r8 real(r8), parameter :: aimm = 0.66_r8 ! Mass of each raindrop created from autoconversion. real(r8), parameter :: droplet_mass_25um = 4._r8/3._r8*pi*rhow*(25.e-6_r8)**3 real(r8), parameter :: droplet_mass_40um = 4._r8/3._r8*pi*rhow*(40.e-6_r8)**3, & droplet_mass_40umi = 1._r8/droplet_mass_40um !========================================================= ! Constants set in initialization !========================================================= ! Set using arguments to micro_mg_init real(r8) :: rv ! water vapor gas constant real(r8) :: cpp ! specific heat of dry air real(r8) :: tmelt ! freezing point of water (K) real(r8) :: ra ! dry air gas constant ! latent heats of: real(r8) :: xxlv ! vaporization real(r8) :: xlf ! freezing real(r8) :: xxls ! sublimation ! additional constants to help speed up code real(r8) :: gamma_bs_plus3 real(r8) :: gamma_half_br_plus5 real(r8) :: gamma_half_bs_plus5 !++ag real(r8) :: gamma_2bs_plus2 !--ag ! real(r8), parameter :: zero = 0._r8, one = 1._r8, two = 2._r8, three = 3._r8, & four = 4._r8, five = 5._r8, six = 6._r8, pio6 = pi/six, & pio3 = pi/three, half = 0.5_r8, oneo3 = one/three, & twopi = pi + pi !========================================================= ! Utilities that are cheaper if the compiler knows that ! some argument is an integer. !========================================================= !>\ingroup micro_mg_utils_mod interface rising_factorial module procedure rising_factorial_r8 module procedure rising_factorial_integer end interface rising_factorial !>\ingroup micro_mg_utils_mod interface var_coef module procedure var_coef_r8 module procedure var_coef_integer end interface var_coef !========================================================================== contains !========================================================================== !>\ingroup micro_mg_utils_mod !! Initialize module variables. ! ! "kind" serves no purpose here except to check for unlikely linking ! issues; always pass in the kind for a double precision real. ! ! ! Check the list at the top of this module for descriptions of all other ! arguments. subroutine micro_mg_utils_init( kind, rair, rh2o, cpair, tmelt_in, latvap, & latice, dcs) ! latice, dcs, errstring) integer, intent(in) :: kind !++ag real(r8), intent(in) :: rair !--ag real(r8), intent(in) :: rh2o real(r8), intent(in) :: cpair real(r8), intent(in) :: tmelt_in real(r8), intent(in) :: latvap real(r8), intent(in) :: latice real(r8), intent(in) :: dcs ! Name this array to workaround an XLF bug (otherwise could just use the ! expression that sets it). real(r8) :: ice_lambda_bounds(2) !----------------------------------------------------------------------- ! declarations for MG code (transforms variable names) rv = rh2o ! water vapor gas constant cpp = cpair ! specific heat of dry air tmelt = tmelt_in !++ag ra = rair ! dry air gas constant !--ag ! latent heats xxlv = latvap ! latent heat vaporization xlf = latice ! latent heat freezing xxls = xxlv + xlf ! latent heat of sublimation ! Define constants to help speed up code (this limits calls to gamma function) gamma_bs_plus3 = gamma(three+bs) gamma_half_br_plus5 = gamma((five+br)*half) gamma_half_bs_plus5 = gamma((five+bs)*half) !++ag gamma_2bs_plus2 = gamma(bs+bs+two) !--ag ! Don't specify lambda bounds for cloud liquid, as they are determined by ! pgam dynamically. mg_liq_props = MGHydrometeorProps(rhow, dsph, min_mean_mass=min_mean_mass_liq) ! Mean ice diameter can not grow bigger than twice the autoconversion ! threshold for snow. ice_lambda_bounds = one/[two*dcs, 1.e-6_r8] mg_ice_props = MGHydrometeorProps(rhoi, dsph, & ice_lambda_bounds, min_mean_mass_ice) mg_rain_props = MGHydrometeorProps(rhow, dsph, lam_bnd_rain) mg_snow_props = MGHydrometeorProps(rhosn, dsph, lam_bnd_snow) !++ag mg_graupel_props = MGHydrometeorProps(rhog, dsph, lam_bnd_snow) !--ag end subroutine micro_mg_utils_init !>\ingroup micro_mg_utils_mod !! Constructor for a constituent property object. function NewMGHydrometeorProps(rho, eff_dim, lambda_bounds, min_mean_mass) & result(res) real(r8), intent(in) :: rho, eff_dim real(r8), intent(in), optional :: lambda_bounds(2), min_mean_mass type(MGHydrometeorProps) :: res res%rho = rho res%eff_dim = eff_dim if (present(lambda_bounds)) then res%lambda_bounds = lambda_bounds else res%lambda_bounds = no_limiter() end if if (present(min_mean_mass)) then res%min_mean_mass = min_mean_mass else res%min_mean_mass = no_limiter() end if res%shape_coef = rho * pio6 * gamma(eff_dim+one) end function NewMGHydrometeorProps !======================================================================== !FORMULAS !======================================================================== ! Use gamma function to implement rising factorial extended to the reals. pure function rising_factorial_r8(x, n) result(res) real(r8), intent(in) :: x, n real(r8) :: res res = gamma(x+n) / gamma(x) end function rising_factorial_r8 ! Rising factorial can be performed much cheaper if n is a small integer. pure function rising_factorial_integer(x, n) result(res) real(r8), intent(in) :: x integer, intent(in) :: n real(r8) :: res integer :: i real(r8) :: factor res = one factor = x do i = 1, n res = res * factor factor = factor + one end do end function rising_factorial_integer ! Calculate correction due to latent heat for evaporation/sublimation elemental function calc_ab(t, qv, xxl) result(ab) real(r8), intent(in) :: t ! Temperature real(r8), intent(in) :: qv ! Saturation vapor pressure real(r8), intent(in) :: xxl ! Latent heat real(r8) :: ab real(r8) :: dqsdt dqsdt = xxl*qv / (rv*t*t) ab = one + dqsdt*xxl/cpp end function calc_ab !>\ingroup micro_mg_utils_mod !! get cloud droplet size distribution parameters elemental subroutine size_dist_param_liq_line(props, qcic, ncic, rho, pgam, lamc) type(MGHydrometeorProps), intent(in) :: props real(r8), intent(in) :: qcic real(r8), intent(inout) :: ncic real(r8), intent(in) :: rho real(r8), intent(out) :: pgam real(r8), intent(out) :: lamc real(r8) :: xs type(MGHydrometeorProps) :: props_loc logical, parameter :: liq_gmao=.true. if (qcic > qsmall) then ! Local copy of properties that can be modified. ! (Elemental routines that operate on arrays can't modify scalar ! arguments.) props_loc = props ! Get pgam from fit to observations of martin et al. 1994 if (liq_gmao) then pgam = 0.0005714_r8*1.e-6_r8*ncic*rho + 0.2714_r8 ! Anning modified lamc if ((ncic > 1.0e-3_r8) .and. (qcic > 1.0e-11_r8)) then xs = 0.07_r8*(1000._r8*qcic/ncic) ** (-0.14_r8) else xs = 1.2_r8 end if xs = max(min(xs, 1.7_r8), 1.1_r8) xs = xs*xs*xs xs = (xs + sqrt(xs+8.0_r8)*sqrt(xs) - 4.0_r8)/8.0_r8 pgam = sqrt(xs) else pgam = one - 0.7_r8 * exp(-0.008_r8*1.e-6_r8*ncic*rho) ! pgam = 0.0005714_r8*1.e-6_r8*ncic*rho + 0.2714_r8 endif pgam = one / (pgam*pgam) - one pgam = max(pgam, two) ! Set coefficient for use in size_dist_param_basic. ! The 3D case is so common and optimizable that we specialize it: if (props_loc%eff_dim == three) then props_loc%shape_coef = pio6 * props_loc%rho * & rising_factorial(pgam+one, 3) else props_loc%shape_coef = pio6 * props_loc%rho * & rising_factorial(pgam+one, props_loc%eff_dim) end if ! Limit to between 2 and 50 microns mean size. props_loc%lambda_bounds = (pgam+one) * one/[50.e-6_r8, 2.e-6_r8] call size_dist_param_basic(props_loc, qcic, ncic, lamc) else ! pgam not calculated in this case, so set it to a value likely to ! cause an error if it is accidentally used ! (gamma function undefined for negative integers) pgam = -100._r8 lamc = zero end if end subroutine size_dist_param_liq_line !>\ingroup micro_mg_utils_mod !! This subroutine gets cloud droplet size distribution parameters subroutine size_dist_param_liq_vect(props, qcic, ncic, rho, pgam, lamc, mgncol) type(mghydrometeorprops), intent(in) :: props integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: qcic real(r8), dimension(mgncol), intent(inout) :: ncic real(r8), dimension(mgncol), intent(in) :: rho real(r8), dimension(mgncol), intent(out) :: pgam real(r8), dimension(mgncol), intent(out) :: lamc real(r8) :: xs type(mghydrometeorprops) :: props_loc logical, parameter :: liq_gmao=.true. integer :: i do i=1,mgncol if (qcic(i) > qsmall) then ! Local copy of properties that can be modified. ! (Elemental routines that operate on arrays can't modify scalar ! arguments.) props_loc = props ! Get pgam from fit to observations of martin et al. 1994 if (liq_gmao) then pgam(i) = 0.0005714_r8*1.e-6_r8*ncic(i)*rho(i) + 0.2714_r8 if ((ncic(i) > 1.0e-3_r8) .and. (qcic(i) > 1.0e-11_r8)) then xs = 0.07_r8*(1000._r8*qcic(i)/ncic(i)) **(-0.14_r8) else xs = 1.2_r8 end if xs = max(min(xs, 1.7_r8), 1.1_r8) xs = xs*xs*xs xs = (xs + sqrt(xs+8.0_r8)*sqrt(xs) - 4.0_r8)/8.0_r8 pgam(i) = sqrt(xs) else pgam(i) = one - 0.7_r8 * exp(-0.008_r8*1.e-6_r8*ncic(i)*rho(i)) ! pgam(i) = 0.0005714_r8*1.e-6_r8*ncic(i)*rho(i) + 0.2714_r8 endif pgam(i) = one/(pgam(i)*pgam(i)) - one pgam(i) = max(pgam(i), two) endif enddo do i=1,mgncol if (qcic(i) > qsmall) then ! Set coefficient for use in size_dist_param_basic. ! The 3D case is so common and optimizable that we specialize ! it: if (props_loc%eff_dim == three) then props_loc%shape_coef = pio6 * props_loc%rho * & rising_factorial(pgam(i)+one, 3) else props_loc%shape_coef = pio6 * props_loc%rho * & rising_factorial(pgam(i)+one, props_loc%eff_dim) end if ! Limit to between 2 and 50 microns mean size. props_loc%lambda_bounds(1) = (pgam(i)+one) / 50.e-6_r8 props_loc%lambda_bounds(2) = (pgam(i)+one) / 2.e-6_r8 call size_dist_param_basic(props_loc, qcic(i), ncic(i), lamc(i)) endif enddo do i=1,mgncol if (qcic(i) <= qsmall) then ! pgam not calculated in this case, so set it to a value likely to ! cause an error if it is accidentally used ! (gamma function undefined for negative integers) pgam(i) = -100._r8 lamc(i) = zero end if enddo end subroutine size_dist_param_liq_vect !>\ingroup micro_mg_utils_mod !! Basic routine for getting size distribution parameters. elemental subroutine size_dist_param_basic_line(props, qic, nic, lam, n0) type(MGHydrometeorProps), intent(in) :: props real(r8), intent(in) :: qic real(r8), intent(inout) :: nic real(r8), intent(out) :: lam real(r8), intent(out), optional :: n0 if (qic > qsmall) then ! add upper limit to in-cloud number concentration to prevent ! numerical error if (limiter_is_on(props%min_mean_mass)) then nic = min(nic, qic / props%min_mean_mass) end if ! lambda = (c n/q)^(1/d) lam = (props%shape_coef * nic/qic)**(one/props%eff_dim) ! check for slope ! adjust vars if (lam < props%lambda_bounds(1)) then lam = props%lambda_bounds(1) nic = lam**(props%eff_dim) * qic/props%shape_coef else if (lam > props%lambda_bounds(2)) then lam = props%lambda_bounds(2) nic = lam**(props%eff_dim) * qic/props%shape_coef end if else lam = zero end if if (present(n0)) n0 = nic * lam end subroutine size_dist_param_basic_line !>\ingroup micro_mg_utils_mod !! This subroutine calculates subroutine size_dist_param_basic_vect(props, qic, nic, lam, mgncol, n0) type (mghydrometeorprops), intent(in) :: props integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: qic real(r8), dimension(mgncol), intent(inout) :: nic real(r8), dimension(mgncol), intent(out) :: lam real(r8), dimension(mgncol), intent(out), optional :: n0 integer :: i do i=1,mgncol if (qic(i) > qsmall) then ! add upper limit to in-cloud number concentration to prevent ! numerical error if (limiter_is_on(props%min_mean_mass)) then nic(i) = min(nic(i), qic(i) / props%min_mean_mass) end if ! lambda = (c n/q)^(1/d) lam(i) = (props%shape_coef * nic(i)/qic(i))**(one/props%eff_dim) ! check for slope ! adjust vars if (lam(i) < props%lambda_bounds(1)) then lam(i) = props%lambda_bounds(1) nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef else if (lam(i) > props%lambda_bounds(2)) then lam(i) = props%lambda_bounds(2) nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef end if else lam(i) = zero end if enddo if (present(n0)) n0 = nic * lam end subroutine size_dist_param_basic_vect !>\ingroup micro_mg_utils_mod !! ice routine for getting size distribution parameters. elemental subroutine size_dist_param_ice_line(props, qic, nic, lam, n0) type(MGHydrometeorProps), intent(in) :: props real(r8), intent(in) :: qic real(r8), intent(inout) :: nic real(r8), intent(out) :: lam real(r8):: miu_ice,tx1,tx2, aux real(r8), intent(out), optional :: n0 logical, parameter :: ice_sep=.true. if (qic > qsmall) then ! add upper limit to in-cloud number concentration to prevent ! numerical error if (limiter_is_on(props%min_mean_mass)) then nic = min(nic, qic / props%min_mean_mass) end if ! lambda = (c n/q)^(1/d) lam = (props%shape_coef * nic/qic)**(1._r8/props%eff_dim) if (ice_sep) then miu_ice = max(min(0.008_r8*(lam*0.01)**0.87_r8, 10.0_r8), 0.1_r8) tx1 = 1.0_r8 + miu_ice tx2 = 1.0_r8 / gamma(tx1) aux = (gamma(tx1+3.0_r8)*tx2) ** (1.0_r8/3.0_r8) lam = lam*aux else aux = 1.0_r8 tx1 = 1.0_r8 tx2 = 1.0_r8 end if if (present(n0)) n0 = nic * lam**tx1*tx2 ! check for slope ! adjust vars if (lam < props%lambda_bounds(1)*aux) then lam = props%lambda_bounds(1) nic = lam**(props%eff_dim) * qic/props%shape_coef if (present(n0)) n0 = nic * lam else if (lam > props%lambda_bounds(2)*aux) then lam = props%lambda_bounds(2) nic = lam**(props%eff_dim) * qic/props%shape_coef if (present(n0)) n0 = nic * lam end if else lam = 0.0_r8 end if end subroutine size_dist_param_ice_line !>\ingroup micro_mg_utils_mod !! This subroutine subroutine size_dist_param_ice_vect(props, qic, nic, lam, mgncol, n0) type (mghydrometeorprops), intent(in) :: props integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: qic real(r8), dimension(mgncol), intent(inout) :: nic real(r8), dimension(mgncol), intent(out) :: lam real(r8), dimension(mgncol), intent(out), optional :: n0 real(r8) :: miu_ice,tx1,tx2, aux integer :: i logical, parameter :: ice_sep=.true. do i=1,mgncol if (qic(i) > qsmall) then ! add upper limit to in-cloud number concentration to prevent ! numerical error if (limiter_is_on(props%min_mean_mass)) then nic(i) = min(nic(i), qic(i) / props%min_mean_mass) end if ! lambda = (c n/q)^(1/d) lam(i) = (props%shape_coef * nic(i)/qic(i))**(1._r8/props%eff_dim) if (ice_sep) then miu_ice = max(min(0.008_r8*(lam(i)*0.01)**0.87_r8, 10.0_r8), 0.1_r8) tx1 = 1.0_r8 + miu_ice tx2 = 1.0_r8 / gamma(tx1) aux = (gamma(tx1+3.0_r8)*tx2) ** (1.0_r8/3.0_r8) lam(i) = lam(i)*aux else aux = 1.0_r8 tx1 = 1.0_r8 tx2 = 1.0_r8 end if if (present(n0)) n0(i) = nic(i) * lam(i)**tx1*tx2 ! check for slope ! adjust vars if (lam(i) < props%lambda_bounds(1)*aux) then lam(i) = props%lambda_bounds(1) nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef if (present(n0)) n0(i) = nic(i) * lam(i) else if (lam(i) > props%lambda_bounds(2)*aux) then lam(i) = props%lambda_bounds(2) nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef if (present(n0)) n0(i) = nic(i) * lam(i) end if else lam(i) = 0.0_r8 end if enddo end subroutine size_dist_param_ice_vect !>\ingroup micro_mg_utils_mod !> Finds the average diameter of particles given their density, and !! mass/number concentrations in the air. !! Assumes that diameter follows an exponential distribution. real(r8) elemental function avg_diameter(q, n, rho_air, rho_sub) real(r8), intent(in) :: q !< mass mixing ratio real(r8), intent(in) :: n !< number concentration (per volume) real(r8), intent(in) :: rho_air !< local density of the air real(r8), intent(in) :: rho_sub !< density of the particle substance avg_diameter = (pi * rho_sub * n/(q*rho_air))**(-oneo3) end function avg_diameter !>\ingroup mg2mg3 !> Finds a coefficient for process rates based on the relative variance !! of cloud water. elemental function var_coef_r8(relvar, a) result(res) real(r8), intent(in) :: relvar real(r8), intent(in) :: a real(r8) :: res res = rising_factorial(relvar, a) / relvar**a end function var_coef_r8 !>\ingroup mg2mg3 !> Finds a coefficient for process rates based on the relative variance !! of cloud water. elemental function var_coef_integer(relvar, a) result(res) real(r8), intent(in) :: relvar integer, intent(in) :: a real(r8) :: res res = rising_factorial(relvar, a) / relvar**a end function var_coef_integer !======================================================================== !MICROPHYSICAL PROCESS CALCULATIONS !======================================================================== !======================================================================== !>\ingroup micro_mg_utils_mod !! Initial ice deposition and sublimation loop. !! Run before the main loop !! This subroutine written by Peter Caldwell subroutine ice_deposition_sublimation(t, qv, qi, ni, & icldm, rho, dv,qvl, qvi, & berg, vap_dep, ice_sublim, mgncol) !INPUT VARS: !=============================================== integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t real(r8), dimension(mgncol), intent(in) :: qv real(r8), dimension(mgncol), intent(in) :: qi real(r8), dimension(mgncol), intent(in) :: ni real(r8), dimension(mgncol), intent(in) :: icldm real(r8), dimension(mgncol), intent(in) :: rho real(r8), dimension(mgncol), intent(in) :: dv real(r8), dimension(mgncol), intent(in) :: qvl real(r8), dimension(mgncol), intent(in) :: qvi !OUTPUT VARS: !=============================================== real(r8), dimension(mgncol), intent(out) :: vap_dep !ice deposition (cell-ave value) real(r8), dimension(mgncol), intent(out) :: ice_sublim !ice sublimation (cell-ave value) real(r8), dimension(mgncol), intent(out) :: berg !bergeron enhancement (cell-ave value) !INTERNAL VARS: !=============================================== real(r8) :: ab real(r8) :: epsi real(r8) :: qiic real(r8) :: niic real(r8) :: lami real(r8) :: n0i real(r8) :: tx1 integer :: i do i=1,mgncol if (qi(i)>=qsmall) then !GET IN-CLOUD qi, ni !=============================================== tx1 = one / icldm(i) qiic = qi(i) * tx1 niic = ni(i) * tx1 !Compute linearized condensational heating correction ab = calc_ab(t(i), qvi(i), xxls) !Get slope and intercept of gamma distn for ice. ! call size_dist_param_basic(mg_ice_props, qiic, niic, lami, n0i) call size_dist_param_ice(mg_ice_props, qiic, niic, lami, n0i) !Get depletion timescale=1/eps epsi = twopi*n0i*rho(i)*Dv(i)/(lami*lami) !Compute deposition/sublimation vap_dep(i) = epsi/ab*(qv(i) - qvi(i)) !Make this a grid-averaged quantity vap_dep(i) = vap_dep(i)*icldm(i) !Split into deposition or sublimation. if (t(i) < tmelt .and. vap_dep(i) > zero) then ice_sublim(i) = zero else ! make ice_sublim negative for consistency with other evap/sub processes ice_sublim(i) = min(vap_dep(i), zero) vap_dep(i) = zero end if !sublimation occurs @ any T. Not so for berg. if (t(i) < tmelt) then !Compute bergeron rate assuming cloud for whole step. berg(i) = max(epsi/ab*(qvl(i) - qvi(i)), zero) else !T>frz berg(i) = zero end if !Tqsmall enddo end subroutine ice_deposition_sublimation !======================================================================== !>\ingroup micro_mg_utils_mod !! autoconversion of cloud liquid water to rain !! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc !! minimum qc of 1 x 10^-8 prevents floating point error subroutine kk2000_liq_autoconversion(microp_uniform, qcic, & ncic, rho, relvar, prc, nprc, nprc1, mgncol) integer, intent(in) :: mgncol logical, intent(in) :: microp_uniform real(r8), dimension(mgncol), intent(in) :: qcic real(r8), dimension(mgncol), intent(in) :: ncic real(r8), dimension(mgncol), intent(in) :: rho real(r8), dimension(mgncol), intent(in) :: relvar real(r8), dimension(mgncol), intent(out) :: prc real(r8), dimension(mgncol), intent(out) :: nprc real(r8), dimension(mgncol), intent(out) :: nprc1 real(r8), dimension(mgncol) :: prc_coef integer :: i ! Take variance into account, or use uniform value. if (.not. microp_uniform) then prc_coef(:) = var_coef(relvar(:), 2.47_r8) else prc_coef(:) = one end if do i=1,mgncol if (qcic(i) >= icsmall) then ! nprc is increase in rain number conc due to autoconversion ! nprc1 is decrease in cloud droplet conc due to autoconversion ! assume exponential sub-grid distribution of qc, resulting in additional ! factor related to qcvar below ! switch for sub-columns, don't include sub-grid qc prc(i) = prc_coef(i) * & 0.01_r8 * 1350._r8 * qcic(i)**2.47_r8 * (ncic(i)*1.e-6_r8*rho(i))**(-1.1_r8) nprc(i) = prc(i) * (one/droplet_mass_25um) nprc1(i) = prc(i)*ncic(i)/qcic(i) else prc(i) = zero nprc(i) = zero nprc1(i) = zero end if enddo end subroutine kk2000_liq_autoconversion !======================================================================== !>\ingroup micro_mg_utils_mod !! This subroutine subroutine sb2001v2_liq_autoconversion(pgam,qc,nc,qr,rho,relvar,au,nprc,nprc1,mgncol) ! ! --------------------------------------------------------------------- ! AUTO_SB: calculates the evolution of mass- and number mxg-ratio for ! drizzle drops due to autoconversion. The autoconversion rate assumes ! f(x)=A*x**(nu_c)*exp(-Bx) in drop MASS x. ! Code from Hugh Morrison, Sept 2014 ! autoconversion ! use simple lookup table of dnu values to get mass spectral shape parameter ! equivalent to the size spectral shape parameter pgam integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent (in) :: pgam real(r8), dimension(mgncol), intent (in) :: qc ! = qc (cld water mixing ratio) real(r8), dimension(mgncol), intent (in) :: nc ! = nc (cld water number conc /kg) real(r8), dimension(mgncol), intent (in) :: qr ! = qr (rain water mixing ratio) real(r8), dimension(mgncol), intent (in) :: rho ! = rho : density profile real(r8), dimension(mgncol), intent (in) :: relvar real(r8), dimension(mgncol), intent (out) :: au ! = prc autoconversion rate real(r8), dimension(mgncol), intent (out) :: nprc1 ! = number tendency real(r8), dimension(mgncol), intent (out) :: nprc ! = number tendency fixed size for rain ! parameters for droplet mass spectral shape, ! used by Seifert and Beheng (2001) ! warm rain scheme only (iparam = 1) real(r8), parameter :: dnu(16) = [0._r8,-0.557_r8,-0.430_r8,-0.307_r8, & -0.186_r8,-0.067_r8,0.050_r8,0.167_r8,0.282_r8,0.397_r8,0.512_r8, & 0.626_r8,0.739_r8,0.853_r8,0.966_r8,0.966_r8] ! parameters for Seifert and Beheng (2001) autoconversion/accretion real(r8), parameter :: kc = 9.44e9_r8 real(r8), parameter :: kr = 5.78e3_r8 real(r8), parameter :: auf = kc / (20._r8*2.6e-7_r8) * 1000._r8, & con_nprc1 = two/2.6e-7_r8*1000._r8 real(r8) :: dum, dum1, nu, pra_coef, tx1, tx2, tx3, tx4 integer :: dumi, i do i=1,mgncol pra_coef = var_coef(relvar(i), 2.47_r8) if (qc(i) > qsmall) then dumi = max(1, min(int(pgam(i)), 15)) nu = dnu(dumi) + (dnu(dumi+1)-dnu(dumi))* (pgam(i)-dumi) !Anning fixed a bug here for FV3GFS 10/13/2017 dum = max(one-qc(i)/(qc(i)+qr(i)), zero) tx1 = dum**0.68_r8 tx2 = one - tx1 dum1 = 600._r8 * tx1 * tx2 * tx2 * tx2 ! Moorthi ! dum1 = 600._r8*dum**0.68_r8*(one-dum**0.68_r8)**3 tx1 = nu + one tx2 = 0.001_r8 * rho(i) * qc(i) tx3 = tx2 * tx2 / (rho(i)*nc(i)*1.e-6_r8) tx2 = tx3 * tx3 tx3 = one - dum au(i) = auf * (nu+two) * (nu+four) * tx2 & * (one+dum1/(tx3*tx3)) / (tx1*tx1*rho(i)) ! au(i) = kc/(20._r8*2.6e-7_r8)* & ! (nu+2._r8)*(nu+4._r8)/(nu+1._r8)**2._r8* & ! (rho(i)*qc(i)/1000._r8)**4._r8/(rho(i)*nc(i)/1.e6_r8)**2._r8* & ! (1._r8+dum1/(1._r8-dum)**2)*1000._r8 / rho(i) ! nprc1(i) = au(i) * two / 2.6e-7_r8 * 1000._r8 ! nprc(i) = au(i) / droplet_mass_40um nprc1(i) = au(i) * con_nprc1 nprc(i) = au(i) * droplet_mass_40umi else au(i) = zero nprc1(i) = zero nprc(i) = zero end if enddo end subroutine sb2001v2_liq_autoconversion !======================================================================== !>\ingroup micro_mg_utils_mod !! Anning Cheng 10/5/2017 add Liu et al. autoconversion subroutine liu_liq_autoconversion(pgam,qc,nc,qr,rho,relvar, & au,nprc,nprc1,mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent (in) :: pgam real(r8), dimension(mgncol), intent (in) :: qc real(r8), dimension(mgncol), intent (in) :: nc real(r8), dimension(mgncol), intent (in) :: qr real(r8), dimension(mgncol), intent (in) :: rho real(r8), dimension(mgncol), intent (in) :: relvar real(r8), dimension(mgncol), intent (out) :: au real(r8), dimension(mgncol), intent (out) :: nprc1 real(r8), dimension(mgncol), intent (out) :: nprc real(r8) :: xs,lw, nw, beta6 ! real(r8), parameter :: dcrit=1.0e-6, miu_disp=1. ! real(r8), parameter :: dcrit=1.0e-3, miu_disp=1. real(r8), parameter :: dcrit = 2.0e-3, miu_disp = 0.8, & con_nprc1 = two/2.6e-7_r8*1000._r8 integer :: i do i=1,mgncol if (qc(i) > qsmall) then xs = one / (one+pgam(i)) beta6 = (one+three*xs)*(one+four*xs)*(one+five*xs) & / ((one+xs)*(one+xs+xs)) LW = 1.0e-3_r8 * qc(i) * rho(i) NW = nc(i) * rho(i) * 1.e-6_r8 xs = min(20.0_r8, 1.03e16_r8*(LW*LW)/(NW*SQRT(NW))) au(i) = 1.1e10_r8*beta6*LW*LW*LW & * (one-exp(-(xs**miu_disp))) / NW au(i) = au(i)*1.0e3_r8/rho(i) au(i) = au(i) * gamma(two+relvar(i)) & / (gamma(relvar(i))*(relvar(i)*relvar(i))) au(i) = au(i) * dcrit ! nprc1(i)= au(i) * (two/2.6e-7_r8*1000._r8) ! nprc(i) = au(i) / droplet_mass_40um nprc1(i)= au(i) * con_nprc1 nprc(i) = au(i) * droplet_mass_40umi else au(i) = zero nprc1(i) = zero nprc(i) = zero end if enddo end subroutine liu_liq_autoconversion !======================================================================== !SB2001 Accretion V2 !>\ingroup micro_mg_utils_mod subroutine sb2001v2_accre_cld_water_rain(qc,nc,qr,rho,relvar,pra,npra,mgncol) ! ! --------------------------------------------------------------------- ! ACCR_SB calculates the evolution of mass mxng-ratio due to accretion ! and self collection following Seifert & Beheng (2001). ! integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent (in) :: qc ! = qc (cld water mixing ratio) real(r8), dimension(mgncol), intent (in) :: nc ! = nc (cld water number conc /kg) real(r8), dimension(mgncol), intent (in) :: qr ! = qr (rain water mixing ratio) real(r8), dimension(mgncol), intent (in) :: rho ! = rho : density profile real(r8), dimension(mgncol), intent (in) :: relvar ! Output tendencies real(r8), dimension(mgncol), intent(out) :: pra ! MMR real(r8), dimension(mgncol), intent(out) :: npra ! Number ! parameters for Seifert and Beheng (2001) autoconversion/accretion real(r8), parameter :: kc = 9.44e9_r8 real(r8), parameter :: kr = 5.78e3_r8 real(r8) :: dum, dum1, tx1, tx2 integer :: i ! accretion do i =1,mgncol if (qc(i) > qsmall) then dum = one - qc(i)/(qc(i)+qr(i)) tx1 = dum / (dum+5.e-4_r8) dum1 = tx1 * tx1 dum1 = dum1 * dum1 pra(i) = kr*rho(i)*0.001_r8*qc(i)*qr(i)*dum1 npra(i) = pra(i) * nc(i) / qc(i) ! npra(i) = pra(i)*rho(i)*0.001_r8*(nc(i)*rho(i)*1.e-6_r8)/ & ! (qc(i)*rho(i)*0.001_r8)*1.e6_r8 / rho(i) else pra(i) = zero npra(i) = zero end if enddo end subroutine sb2001v2_accre_cld_water_rain !======================================================================== ! Autoconversion of cloud ice to snow ! similar to Ferrier (1994) !>\ingroup micro_mg_utils_mod !! Autoconversion of cloud ice to snow !! similar to Ferrier (1994) subroutine ice_autoconversion(t, qiic, lami, n0i, dcs, ac_time, prci, nprci, mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t real(r8), dimension(mgncol), intent(in) :: qiic real(r8), dimension(mgncol), intent(in) :: lami real(r8), dimension(mgncol), intent(in) :: n0i real(r8), intent(in) :: dcs real(r8), dimension(mgncol), intent(in) :: ac_time real(r8), dimension(mgncol), intent(out) :: prci real(r8), dimension(mgncol), intent(out) :: nprci ! Assume autoconversion timescale of 180 seconds. ! Average mass of an ice particle. real(r8) :: m_ip ! Ratio of autoconversion diameter to average diameter. real(r8) :: d_rat integer :: i do i=1,mgncol if (t(i) <= tmelt .and. qiic(i) >= qsmall) then d_rat = lami(i)*dcs ! Rate of ice particle conversion (number). nprci(i) = n0i(i)/(lami(i)*ac_time(i))*exp(-d_rat) m_ip = rhoi * pio6 / (lami(i)*lami(i)*lami(i)) ! m_ip = (rhoi*pi/6._r8) / lami(i)**3 ! Rate of mass conversion. ! Note that this is: ! m n (d^3 + 3 d^2 + 6 d + 6) prci(i) = m_ip * nprci(i) * (((d_rat + three)*d_rat + six)*d_rat + six) else prci(i) = zero nprci(i) = zero end if enddo end subroutine ice_autoconversion !=================================== ! Anning Cheng 10/5/2017 added GMAO ice autoconversion !>\ingroup micro_mg_utils_mod !! GMAO ice autoconversion subroutine gmao_ice_autoconversion(t, qiic, niic, lami, n0i, & dcs, ac_time, prci, nprci, mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t real(r8), dimension(mgncol), intent(in) :: qiic real(r8), dimension(mgncol), intent(in) :: niic real(r8), dimension(mgncol), intent(in) :: lami real(r8), dimension(mgncol), intent(in) :: n0i real(r8), dimension(mgncol), intent(in) :: ac_time real(r8), intent(in) :: dcs real(r8), dimension(mgncol), intent(out) :: prci real(r8), dimension(mgncol), intent(out) :: nprci real(r8) :: m_ip, tx1, tx2 integer :: i do i=1,mgncol if (t(i) <= tmelt .and. qiic(i) >= qsmall) then m_ip = max(min(0.008_r8*(lami(i)*0.01)**0.87_r8, & 10.0_r8), 0.1_r8) tx1 = lami(i)*dcs tx2 = one / ac_time(i) nprci(i) = niic(i)*tx2 * (one - gamma_incomp(m_ip, tx1)) prci(i) = qiic(i)*tx2 * (one - gamma_incomp(m_ip+three, tx1)) else prci(i) = zero nprci(i) = zero end if enddo end subroutine gmao_ice_autoconversion !=================================== ! immersion freezing (Bigg, 1953) !=================================== !>\ingroup micro_mg_utils_mod !! immersion freezing (Bigg, 1953) subroutine immersion_freezing(microp_uniform, t, pgam, lamc, & qcic, ncic, relvar, mnuccc, nnuccc, mgncol) integer, intent(in) :: mgncol logical, intent(in) :: microp_uniform ! Temperature real(r8), dimension(mgncol), intent(in) :: t ! Cloud droplet size distribution parameters real(r8), dimension(mgncol), intent(in) :: pgam real(r8), dimension(mgncol), intent(in) :: lamc ! MMR and number concentration of in-cloud liquid water real(r8), dimension(mgncol), intent(in) :: qcic real(r8), dimension(mgncol), intent(in) :: ncic ! Relative variance of cloud water real(r8), dimension(mgncol), intent(in) :: relvar ! Output tendencies real(r8), dimension(mgncol), intent(out) :: mnuccc ! MMR real(r8), dimension(mgncol), intent(out) :: nnuccc ! Number ! Coefficients that will be omitted for sub-columns real(r8), dimension(mgncol) :: dum real(r8) :: tx1 integer :: i if (.not. microp_uniform) then dum(:) = var_coef(relvar, 2) else dum(:) = one end if do i=1,mgncol if (qcic(i) >= qsmall .and. t(i) < 269.15_r8) then tx1 = one / (lamc(i) * lamc(i) * lamc(i)) nnuccc(i) = pio6*ncic(i)*rising_factorial(pgam(i)+one, 3) * & bimm*(exp(aimm*(tmelt - t(i)))-one) * tx1 mnuccc(i) = dum(i) * nnuccc(i) * pio6 * rhow * & rising_factorial(pgam(i)+four, 3) * tx1 else mnuccc(i) = zero nnuccc(i) = zero end if ! qcic > qsmall and t < 4 deg C enddo end subroutine immersion_freezing !>\ingroup micro_mg_utils_mod !! contact freezing (-40= qsmall .and. t(i) < 269.15_r8) then if (.not. microp_uniform) then dum = var_coef(relvar(i), four/three) dum1 = var_coef(relvar(i), oneo3) else dum = one dum1 = one endif tcnt=(270.16_r8-t(i))**1.3_r8 viscosity = 1.8e-5_r8*(t(i)/298.0_r8)**0.85_r8 ! Viscosity (kg/m/s) mfp = two*viscosity/ & ! Mean free path (m) (p(i)*sqrt( 8.0_r8*28.96e-3_r8/(pi*8.314409_r8*t(i)) )) ! Note that these two are vectors. nslip = one+(mfp/rndst(i,:))*(1.257_r8+(0.4_r8*exp(-(1.1_r8*rndst(i,:)/mfp))))! Slip correction factor ndfaer = 1.381e-23_r8*t(i)*nslip/(6._r8*pi*viscosity*rndst(i,:)) ! aerosol diffusivity (m2/s) tx1 = one / lamc(i) contact_factor = dot_product(ndfaer,nacon(i,:)*tcnt) * pi * & ncic(i) * (pgam(i) + one) * tx1 mnucct(i) = dum * contact_factor * & pio3*rhow*rising_factorial(pgam(i)+two, 3) * tx1 * tx1 *tx1 nnucct(i) = (dum1+dum1) * contact_factor else mnucct(i) = zero nnucct(i) = zero end if ! qcic > qsmall and t < 4 deg C end do end subroutine contact_freezing !>\ingroup micro_mg_utils_mod !! snow self-aggregation from passarelli, 1978, used by reisner, 1998 !=================================================================== ! this is hard-wired for bs = 0.4 for now ! ignore self-collection of cloud ice subroutine snow_self_aggregation(t, rho, asn, rhosn, qsic, nsic, nsagg, mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t ! Temperature real(r8), dimension(mgncol), intent(in) :: rho ! Density real(r8), dimension(mgncol), intent(in) :: asn ! fall speed parameter for snow real(r8), intent(in) :: rhosn ! density of snow ! In-cloud snow real(r8), dimension(mgncol), intent(in) :: qsic ! MMR real(r8), dimension(mgncol), intent(in) :: nsic ! Number ! Output number tendency real(r8), dimension(mgncol), intent(out) :: nsagg integer :: i do i=1,mgncol if (qsic(i) >= qsmall .and. t(i) <= tmelt) then nsagg(i) = -1108._r8*eii/(four*720._r8*rhosn)*asn(i)*qsic(i)*nsic(i)*rho(i)*& ((qsic(i)/nsic(i))*(one/(rhosn*pi)))**((bs-one)*oneo3) else nsagg(i) = zero end if enddo end subroutine snow_self_aggregation !>\ingroup micro_mg_utils_mod !! accretion of cloud droplets onto snow/graupel !=================================================================== ! here use continuous collection equation with ! simple gravitational collection kernel ! ignore collisions between droplets/cloud ice ! since minimum size ice particle for accretion is 50 - 150 micron subroutine accrete_cloud_water_snow(t, rho, asn, uns, mu, qcic, ncic, qsic, & pgam, lamc, lams, n0s, psacws, npsacws, mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t ! Temperature real(r8), dimension(mgncol), intent(in) :: rho ! Density real(r8), dimension(mgncol), intent(in) :: asn ! Fallspeed parameter (snow) real(r8), dimension(mgncol), intent(in) :: uns ! Current fallspeed (snow) real(r8), dimension(mgncol), intent(in) :: mu ! Viscosity ! In-cloud liquid water real(r8), dimension(mgncol), intent(in) :: qcic ! MMR real(r8), dimension(mgncol), intent(in) :: ncic ! Number ! In-cloud snow real(r8), dimension(mgncol), intent(in) :: qsic ! MMR ! Cloud droplet size parameters real(r8), dimension(mgncol), intent(in) :: pgam real(r8), dimension(mgncol), intent(in) :: lamc ! Snow size parameters real(r8), dimension(mgncol), intent(in) :: lams real(r8), dimension(mgncol), intent(in) :: n0s ! Output tendencies real(r8), dimension(mgncol), intent(out) :: psacws ! Mass mixing ratio real(r8), dimension(mgncol), intent(out) :: npsacws ! Number concentration real(r8) :: dc0 ! Provisional mean droplet size real(r8) :: dum real(r8) :: eci ! collection efficiency for riming of snow by droplets ! Fraction of cloud droplets accreted per second real(r8) :: accrete_rate integer :: i ! ignore collision of snow with droplets above freezing do i=1,mgncol if (qsic(i) >= qsmall .and. t(i) <= tmelt .and. qcic(i) >= qsmall) then ! put in size dependent collection efficiency ! mean diameter of snow is area-weighted, since ! accretion is function of crystal geometric area ! collection efficiency is approximation based on stoke's law (Thompson et al. 2004) dc0 = (pgam(i)+one)/lamc(i) dum = dc0*dc0*uns(i)*rhow*lams(i)/(9._r8*mu(i)) eci = dum*dum / ((dum+0.4_r8)*(dum+0.4_r8)) eci = max(eci,zero) eci = min(eci,one) ! no impact of sub-grid distribution of qc since psacws ! is linear in qc accrete_rate = (pi/four)*asn(i)*rho(i)*n0s(i)*eci*gamma_bs_plus3 / lams(i)**(bs+three) psacws(i) = accrete_rate*qcic(i) npsacws(i) = accrete_rate*ncic(i) else psacws(i) = zero npsacws(i) = zero end if enddo end subroutine accrete_cloud_water_snow !>\ingroup micro_mg_utils_mod !! add secondary ice production due to accretion of droplets by snow !=================================================================== ! (Hallet-Mossop process) (from Cotton et al., 1986) subroutine secondary_ice_production(t, psacws, msacwi, nsacwi, mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t ! Temperature ! Accretion of cloud water to snow tendencies real(r8), dimension(mgncol), intent(inout) :: psacws ! MMR ! Output (ice) tendencies real(r8), dimension(mgncol), intent(out) :: msacwi ! MMR real(r8), dimension(mgncol), intent(out) :: nsacwi ! Number integer :: i do i=1,mgncol if((t(i) < 270.16_r8) .and. (t(i) >= 268.16_r8)) then nsacwi(i) = 3.5e8_r8*(270.16_r8-t(i))/two*psacws(i) else if((t(i) < 268.16_r8) .and. (t(i) >= 265.16_r8)) then nsacwi(i) = 3.5e8_r8*(t(i)-265.16_r8)*oneo3*psacws(i) else nsacwi(i) = zero endif enddo do i=1,mgncol msacwi(i) = min(nsacwi(i)*mi0, psacws(i)) psacws(i) = psacws(i) - msacwi(i) enddo end subroutine secondary_ice_production !>\ingroup micro_mg_utils_mod !! accretion of rain water by snow !=================================================================== ! formula from ikawa and saito, 1991, used by reisner et al., 1998 subroutine accrete_rain_snow(t, rho, umr, ums, unr, uns, qric, qsic, & lamr, n0r, lams, n0s, pracs, npracs, mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t ! Temperature real(r8), dimension(mgncol), intent(in) :: rho ! Density ! Fallspeeds ! mass-weighted real(r8), dimension(mgncol), intent(in) :: umr ! rain real(r8), dimension(mgncol), intent(in) :: ums ! snow ! number-weighted real(r8), dimension(mgncol), intent(in) :: unr ! rain real(r8), dimension(mgncol), intent(in) :: uns ! snow ! In cloud MMRs real(r8), dimension(mgncol), intent(in) :: qric ! rain real(r8), dimension(mgncol), intent(in) :: qsic ! snow ! Size distribution parameters ! rain real(r8), dimension(mgncol), intent(in) :: lamr real(r8), dimension(mgncol), intent(in) :: n0r ! snow real(r8), dimension(mgncol), intent(in) :: lams real(r8), dimension(mgncol), intent(in) :: n0s ! Output tendencies real(r8), dimension(mgncol), intent(out) :: pracs ! MMR real(r8), dimension(mgncol), intent(out) :: npracs ! Number ! Collection efficiency for accretion of rain by snow real(r8), parameter :: ecr = one ! Ratio of average snow diameter to average rain diameter. real(r8) :: d_rat ! Common factor between mass and number expressions real(r8) :: common_factor real(r8) :: tx1, tx2 integer :: i do i=1,mgncol if (qric(i) >= icsmall .and. qsic(i) >= icsmall .and. t(i) <= tmelt) then tx2 = lamr(i)*lamr(i)*lamr(i) common_factor = pi*ecr*rho(i)*n0r(i)*n0s(i) / (tx2 * lams(i)) d_rat = lamr(i)/lams(i) tx1 = 1.2_r8*umr(i)-0.95_r8*ums(i) pracs(i) = common_factor*pi*rhow* & sqrt(tx1*tx1 + 0.08_r8*ums(i)*umr(i)) * & ((half*d_rat + two)*d_rat + five) / tx2 tx1 = unr(i)-uns(i) npracs(i) = common_factor*half * & sqrt(1.7_r8*tx1*tx1 + 0.3_r8*unr(i)*uns(i)) * & ((d_rat + one)*d_rat + one) else pracs(i) = zero npracs(i) = zero end if enddo end subroutine accrete_rain_snow !>\ingroup micro_mg_utils_mod !! heterogeneous freezing of rain drops !=================================================================== ! follows from Bigg (1953) subroutine heterogeneous_rain_freezing(t, qric, nric, lamr, mnuccr, nnuccr, mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t ! Temperature ! In-cloud rain real(r8), dimension(mgncol), intent(in) :: qric ! MMR real(r8), dimension(mgncol), intent(in) :: nric ! Number real(r8), dimension(mgncol), intent(in) :: lamr ! size parameter ! Output tendencies real(r8), dimension(mgncol), intent(out) :: mnuccr ! MMR real(r8), dimension(mgncol), intent(out) :: nnuccr ! Number real(r8) :: tx1 integer :: i do i=1,mgncol if (t(i) < 269.15_r8 .and. qric(i) >= qsmall) then tx1 = pi / (lamr(i)*lamr(i)*lamr(i)) nnuccr(i) = nric(i)*bimm* (exp(aimm*(tmelt - t(i)))-one) * tx1 mnuccr(i) = nnuccr(i) * 20._r8*rhow * tx1 else mnuccr(i) = zero nnuccr(i) = zero end if enddo end subroutine heterogeneous_rain_freezing !>\ingroup micro_mg_utils_mod !! accretion of cloud liquid water by rain !! formula from Khrouditnov and Kogan (2000) ! gravitational collection kernel, droplet fall speed neglected subroutine accrete_cloud_water_rain(microp_uniform, qric, qcic, & ncic, relvar, accre_enhan, pra, npra, mgncol) logical, intent(in) :: microp_uniform integer, intent(in) :: mgncol ! In-cloud rain real(r8), dimension(mgncol), intent(in) :: qric ! MMR ! Cloud droplets real(r8), dimension(mgncol), intent(in) :: qcic ! MMR real(r8), dimension(mgncol), intent(in) :: ncic ! Number ! SGS variability real(r8), dimension(mgncol), intent(in) :: relvar real(r8), dimension(mgncol), intent(in) :: accre_enhan ! Output tendencies real(r8), dimension(mgncol), intent(out) :: pra ! MMR real(r8), dimension(mgncol), intent(out) :: npra ! Number ! Coefficient that varies for subcolumns real(r8), dimension(mgncol) :: pra_coef integer :: i if (.not. microp_uniform) then pra_coef(:) = accre_enhan * var_coef(relvar(:), 1.15_r8) else pra_coef(:) = one end if do i=1,mgncol if (qric(i) >= qsmall .and. qcic(i) >= qsmall) then ! include sub-grid distribution of cloud water pra(i) = pra_coef(i) * 67._r8*(qcic(i)*qric(i))**1.15_r8 npra(i) = pra(i)*ncic(i)/qcic(i) else pra(i) = zero npra(i) = zero end if end do end subroutine accrete_cloud_water_rain !>\ingroup micro_mg_utils_mod !! Self-collection of rain drops !! from Beheng(1994) subroutine self_collection_rain(rho, qric, nric, nragg, mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: rho ! Air density ! Rain real(r8), dimension(mgncol), intent(in) :: qric ! MMR real(r8), dimension(mgncol), intent(in) :: nric ! Number ! Output number tendency real(r8), dimension(mgncol), intent(out) :: nragg integer :: i do i=1,mgncol if (qric(i) >= qsmall) then nragg(i) = -8._r8*nric(i)*qric(i)*rho(i) else nragg(i) = zero end if enddo end subroutine self_collection_rain !>\ingroup micro_mg_utils_mod !! Accretion of cloud ice by snow !=================================================================== ! For this calculation, it is assumed that the Vs >> Vi ! and Ds >> Di for continuous collection subroutine accrete_cloud_ice_snow(t, rho, asn, qiic, niic, qsic, & lams, n0s, prai, nprai, mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t ! Temperature real(r8), dimension(mgncol), intent(in) :: rho ! Density real(r8), dimension(mgncol), intent(in) :: asn ! Snow fallspeed parameter ! Cloud ice real(r8), dimension(mgncol), intent(in) :: qiic ! MMR real(r8), dimension(mgncol), intent(in) :: niic ! Number real(r8), dimension(mgncol), intent(in) :: qsic ! Snow MMR ! Snow size parameters real(r8), dimension(mgncol), intent(in) :: lams real(r8), dimension(mgncol), intent(in) :: n0s ! Output tendencies real(r8), dimension(mgncol), intent(out) :: prai ! MMR real(r8), dimension(mgncol), intent(out) :: nprai ! Number ! Fraction of cloud ice particles accreted per second real(r8) :: accrete_rate integer :: i do i=1,mgncol if (qsic(i) >= qsmall .and. qiic(i) >= qsmall .and. t(i) <= tmelt) then accrete_rate = (pi/four) * eii * asn(i) * rho(i) * n0s(i) * gamma_bs_plus3 & / lams(i)**(bs+three) prai(i) = accrete_rate * qiic(i) nprai(i) = accrete_rate * niic(i) else prai(i) = zero nprai(i) = zero end if enddo end subroutine accrete_cloud_ice_snow !>\ingroup micro_mg_utils_mod !! calculate evaporation/sublimation of rain and snow !=================================================================== ! note: evaporation/sublimation occurs only in cloud-free portion of grid cell ! in-cloud condensation/deposition of rain and snow is neglected ! except for transfer of cloud water to snow through bergeron process subroutine 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) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t ! temperature real(r8), dimension(mgncol), intent(in) :: rho ! air density real(r8), dimension(mgncol), intent(in) :: dv ! water vapor diffusivity real(r8), dimension(mgncol), intent(in) :: mu ! viscosity real(r8), dimension(mgncol), intent(in) :: sc ! schmidt number real(r8), dimension(mgncol), intent(in) :: q ! humidity real(r8), dimension(mgncol), intent(in) :: qvl ! saturation humidity (water) real(r8), dimension(mgncol), intent(in) :: qvi ! saturation humidity (ice) real(r8), dimension(mgncol), intent(in) :: lcldm ! liquid cloud fraction real(r8), dimension(mgncol), intent(in) :: precip_frac ! precipitation fraction (maximum overlap) ! fallspeed parameters real(r8), dimension(mgncol), intent(in) :: arn ! rain real(r8), dimension(mgncol), intent(in) :: asn ! snow ! In-cloud MMRs real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid real(r8), dimension(mgncol), intent(in) :: qiic ! cloud ice real(r8), dimension(mgncol), intent(in) :: qric ! rain real(r8), dimension(mgncol), intent(in) :: qsic ! snow ! Size parameters ! rain real(r8), dimension(mgncol), intent(in) :: lamr real(r8), dimension(mgncol), intent(in) :: n0r ! snow real(r8), dimension(mgncol), intent(in) :: lams real(r8), dimension(mgncol), intent(in) :: n0s ! Output tendencies real(r8), dimension(mgncol), intent(out) :: pre real(r8), dimension(mgncol), intent(out) :: prds real(r8), dimension(mgncol), intent(out) :: am_evp_st ! Fractional area where rain evaporates. real(r8) :: qclr ! water vapor mixing ratio in clear air real(r8) :: ab ! correction to account for latent heat real(r8) :: eps ! 1/ sat relaxation timescale real(r8) :: tx1, tx2, tx3 real(r8), dimension(mgncol) :: dum integer :: i am_evp_st = zero ! set temporary cloud fraction to zero if cloud water + ice is very small ! this will ensure that evaporation/sublimation of precip occurs over ! entire grid cell, since min cloud fraction is specified otherwise do i=1,mgncol if (qcic(i)+qiic(i) < 1.e-6_r8) then dum(i) = zero else dum(i) = lcldm(i) end if enddo do i=1,mgncol ! only calculate if there is some precip fraction > cloud fraction if (precip_frac(i) > dum(i)) then if (qric(i) >= qsmall .or. qsic(i) >= qsmall) then am_evp_st(i) = precip_frac(i) - dum(i) ! calculate q for out-of-cloud region qclr = (q(i)-dum(i)*qvl(i)) / (one-dum(i)) end if ! evaporation of rain if (qric(i) >= qsmall) then ab = calc_ab(t(i), qvl(i), xxlv) eps = two*pi*n0r(i)*rho(i)*Dv(i) * & (f1r/(lamr(i)*lamr(i)) + & f2r*sqrt(arn(i)*rho(i)/mu(i)) * & sc(i)**oneo3*gamma_half_br_plus5 & / (lamr(i)**((five+br)*half))) pre(i) = eps*(qclr-qvl(i)) / ab ! only evaporate in out-of-cloud region ! and distribute across precip_frac pre(i) = min(pre(i)*am_evp_st(i), zero) pre(i) = pre(i) / precip_frac(i) else pre(i) = zero end if ! sublimation of snow if (qsic(i) >= qsmall) then ab = calc_ab(t(i), qvi(i), xxls) eps = two*pi*n0s(i)*rho(i)*Dv(i) * & ( f1s/(lams(i)*lams(i)) & + f2s*sqrt(asn(i)*rho(i)/mu(i)) * & sc(i)**oneo3*gamma_half_bs_plus5 & / (lams(i)**((five+bs)*half))) prds(i) = eps*(qclr-qvi(i)) / ab ! only sublimate in out-of-cloud region and distribute over precip_frac prds(i) = min(prds(i)*am_evp_st(i), zero) prds(i) = prds(i) / precip_frac(i) else prds(i) = zero end if else prds(i) = zero pre(i) = zero end if enddo end subroutine evaporate_sublimate_precip !>\ingroup micro_mg_utils_mod !! evaporation/sublimation of rain, snow and graupel !=================================================================== ! note: evaporation/sublimation occurs only in cloud-free portion of grid cell ! in-cloud condensation/deposition of rain and snow is neglected ! except for transfer of cloud water to snow through bergeron process subroutine 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) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t ! temperature real(r8), dimension(mgncol), intent(in) :: rho ! air density real(r8), dimension(mgncol), intent(in) :: dv ! water vapor diffusivity real(r8), dimension(mgncol), intent(in) :: mu ! viscosity real(r8), dimension(mgncol), intent(in) :: sc ! schmidt number real(r8), dimension(mgncol), intent(in) :: q ! humidity real(r8), dimension(mgncol), intent(in) :: qvl ! saturation humidity (water) real(r8), dimension(mgncol), intent(in) :: qvi ! saturation humidity (ice) real(r8), dimension(mgncol), intent(in) :: lcldm ! liquid cloud fraction real(r8), dimension(mgncol), intent(in) :: precip_frac ! precipitation fraction (maximum overlap) ! fallspeed parameters real(r8), dimension(mgncol), intent(in) :: arn ! rain real(r8), dimension(mgncol), intent(in) :: asn ! snow !++ag real(r8), dimension(mgncol), intent(in) :: agn ! graupel real(r8), intent(in) :: bg !--ag ! In-cloud MMRs real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid real(r8), dimension(mgncol), intent(in) :: qiic ! cloud ice real(r8), dimension(mgncol), intent(in) :: qric ! rain real(r8), dimension(mgncol), intent(in) :: qsic ! snow real(r8), dimension(mgncol), intent(in) :: qgic ! graupel ! Size parameters ! rain real(r8), dimension(mgncol), intent(in) :: lamr real(r8), dimension(mgncol), intent(in) :: n0r ! snow real(r8), dimension(mgncol), intent(in) :: lams real(r8), dimension(mgncol), intent(in) :: n0s !++ag ! graupel real(r8), dimension(mgncol), intent(in) :: lamg real(r8), dimension(mgncol), intent(in) :: n0g !--ag ! Output tendencies real(r8), dimension(mgncol), intent(out) :: pre real(r8), dimension(mgncol), intent(out) :: prds !++ag real(r8), dimension(mgncol), intent(out) :: prdg !--ag real(r8), dimension(mgncol), intent(out) :: am_evp_st ! Fractional area where rain evaporates. real(r8) :: qclr ! water vapor mixing ratio in clear air real(r8) :: ab ! correction to account for latent heat real(r8) :: eps ! 1/ sat relaxation timescale real(r8), dimension(mgncol) :: dum integer :: i ! set temporary cloud fraction to zero if cloud water + ice is very small ! this will ensure that evaporation/sublimation of precip occurs over ! entire grid cell, since min cloud fraction is specified otherwise do i=1,mgncol if (qcic(i)+qiic(i) < 1.e-6_r8) then dum(i) = zero else dum(i) = lcldm(i) end if enddo do i=1,mgncol ! only calculate if there is some precip fraction > cloud fraction if (precip_frac(i) > dum(i)) then if (qric(i) >= qsmall .or. qsic(i) >= qsmall .or. qgic(i) >= qsmall) then am_evp_st(i) = precip_frac(i) - dum(i) ! calculate q for out-of-cloud region qclr = (q(i)-dum(i)*qvl(i)) / (one-dum(i)) end if ! evaporation of rain if (qric(i) >= qsmall) then ab = calc_ab(t(i), qvl(i), xxlv) eps = twopi*n0r(i)*rho(i)*Dv(i)* & ( f1r/(lamr(i)*lamr(i)) & + f2r*sqrt(arn(i)*rho(i)/mu(i)) & * sc(i)**oneo3*gamma_half_br_plus5 & / (lamr(i)**((five+br)*half))) pre(i) = eps*(qclr-qvl(i))/ab ! only evaporate in out-of-cloud region ! and distribute across precip_frac pre(i) = min(pre(i)*am_evp_st(i), zero) pre(i) = pre(i)/precip_frac(i) else pre(i) = zero end if ! sublimation of snow if (qsic(i) >= qsmall) then ab = calc_ab(t(i), qvi(i), xxls) eps = twopi*n0s(i)*rho(i)*Dv(i)* & ( f1s/(lams(i)*lams(i)) & + f2s*sqrt(asn(i)*rho(i)/mu(i)) & * sc(i)**oneo3*gamma_half_bs_plus5 & / (lams(i)**((five+bs)*half))) prds(i) = eps*(qclr-qvi(i))/ab ! only sublimate in out-of-cloud region and distribute over precip_frac prds(i) = min(prds(i)*am_evp_st(i), zero) prds(i) = prds(i)/precip_frac(i) else prds(i) = zero end if !++AG ADD GRAUPEL, do Same with prdg. if (qgic(i) >= qsmall) then ab = calc_ab(t(i), qvi(i), xxls) eps = twopi*n0g(i)*rho(i)*Dv(i)* & ( f1s/(lamg(i)*lamg(i)) & + f2s*sqrt(agn(i)*rho(i)/mu(i)) & * sc(i)**oneo3*gamma((five+bg)*half) & / (lamg(i)**((five+bs)*half))) ! / (lamg(i)**((five+bg)*half))) ! changing bs to bg - Moorthi prdg(i) = eps*(qclr-qvi(i))/ab ! only sublimate in out-of-cloud region and distribute over precip_frac prdg(i) = min(prdg(i)*am_evp_st(i), zero) prdg(i) = prdg(i)/precip_frac(i) else prdg(i) = zero end if else prds(i) = zero pre(i) = zero !++ag prdg(i) = zero !--ag end if enddo end subroutine evaporate_sublimate_precip_graupel !>\ingroup micro_mg_utils_mod !! bergeron process - evaporation of droplets and deposition onto snow subroutine bergeron_process_snow(t, rho, dv, mu, sc, qvl, qvi, asn, & qcic, qsic, lams, n0s, bergs, mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t ! temperature real(r8), dimension(mgncol), intent(in) :: rho ! air density real(r8), dimension(mgncol), intent(in) :: dv ! water vapor diffusivity real(r8), dimension(mgncol), intent(in) :: mu ! viscosity real(r8), dimension(mgncol), intent(in) :: sc ! schmidt number real(r8), dimension(mgncol), intent(in) :: qvl ! saturation humidity (water) real(r8), dimension(mgncol), intent(in) :: qvi ! saturation humidity (ice) ! fallspeed parameter for snow real(r8), dimension(mgncol), intent(in) :: asn ! In-cloud MMRs real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid mixing ratio real(r8), dimension(mgncol), intent(in) :: qsic ! snow mixing ratio ! Size parameters for snow real(r8), dimension(mgncol), intent(in) :: lams real(r8), dimension(mgncol), intent(in) :: n0s ! Output tendencies real(r8), dimension(mgncol), intent(out) :: bergs real(r8) :: ab ! correction to account for latent heat real(r8) :: eps ! 1/ sat relaxation timescale integer :: i do i=1,mgncol if (qsic(i) >= qsmall.and. qcic(i) >= qsmall .and. t(i) < tmelt) then ab = calc_ab(t(i), qvi(i), xxls) eps = two*pi*n0s(i)*rho(i)*Dv(i) * & (f1s/(lams(i)*lams(i)) + & f2s*sqrt(asn(i)*rho(i)/mu(i)) * & sc(i)**oneo3*gamma_half_bs_plus5 / & (lams(i)**((five+bs)*half))) bergs(i) = eps*(qvl(i)-qvi(i)) / ab else bergs(i) = zero end if enddo end subroutine bergeron_process_snow !======================================================================== !>\ingroup micro_mg_utils_mod !! Collection of snow by rain to form graupel subroutine graupel_collecting_snow(qsic,qric,umr,ums,rho,lamr,n0r,lams,n0s, & psacr, mgncol) integer, intent(in) :: mgncol ! In-cloud MMRs real(r8), dimension(mgncol), intent(in) :: qsic ! snow real(r8), dimension(mgncol), intent(in) :: qric ! rain ! mass-weighted fall speeds real(r8), dimension(mgncol), intent(in) :: umr ! rain real(r8), dimension(mgncol), intent(in) :: ums ! snow real(r8), dimension(mgncol), intent(in) :: rho ! air density ! Size parameters for rain real(r8), dimension(mgncol), intent(in) :: lamr real(r8), dimension(mgncol), intent(in) :: n0r ! Size parameters for snow real(r8), dimension(mgncol), intent(in) :: lams real(r8), dimension(mgncol), intent(in) :: n0s real(r8), dimension(mgncol), intent(out) :: psacr ! conversion due to coll of snow by rain real(r8) :: cons31, tx1, tx2, tx3, tx4, tx5 integer :: i cons31 = pi*pi*ecr*rhosn do i=1,mgncol if (qsic(i) >= 0.1e-3_r8 .and. qric(i) >= 0.1e-3_r8) then tx1 = 1.2_r8*umr(i) - 0.95_r8*ums(i) tx1 = sqrt(tx1*tx1+0.08_r8*ums(i)*umr(i)) tx2 = one / lams(i) tx3 = one / lamr(i) tx4 = tx2 * tx2 tx5 = tx4 * tx4 * tx3 psacr(i) = cons31 * tx1 * rho(i) * n0r(i) * n0s(i) * tx5 & * (5.0_r8*tx4+tx3*(tx2+tx2+0.5_r8*tx3)) ! psacr(i) = cons31*(((1.2_r8*umr(i)-0.95_r8*ums(i))**2+ & ! 0.08_r8*ums(i)*umr(i))**0.5_r8*rho(i)* & ! n0r(i)*n0s(i)/lams(i)**3* & ! (5._r8/(lams(i)**3*lamr(i))+ & ! 2._r8/(lams(i)**2*lamr(i)**2)+ & ! 0.5_r8/(lams(i)*lamr(i)**3))) else psacr(i) = zero end if end do end subroutine graupel_collecting_snow !======================================================================== !>\ingroup micro_mg_utils_mod !! Collection of cloud water by graupel subroutine graupel_collecting_cld_water(qgic,qcic,ncic,rho,n0g,lamg,bg,agn, & psacwg, npsacwg, mgncol) integer, intent(in) :: mgncol ! In-cloud MMRs real(r8), dimension(mgncol), intent(in) :: qgic ! graupel real(r8), dimension(mgncol), intent(in) :: qcic ! cloud water real(r8), dimension(mgncol), intent(in) :: ncic ! cloud water number conc real(r8), dimension(mgncol), intent(in) :: rho ! air density ! Size parameters for graupel real(r8), dimension(mgncol), intent(in) :: lamg real(r8), dimension(mgncol), intent(in) :: n0g ! fallspeed parameters for graupel ! Set AGN and BG as input (in micro_mg3_0.F90) (select hail or graupel as appropriate) real(r8), intent(in) :: bg real(r8), dimension(mgncol), intent(in) :: agn ! Output tendencies real(r8), dimension(mgncol), intent(out) :: psacwg real(r8), dimension(mgncol), intent(out) :: npsacwg real(r8) :: cons, tx1 integer :: i cons = gamma(bg+three) * pi/four * ecid do i=1,mgncol if (qgic(i) >= 1.e-8_r8 .and. qcic(i) >= qsmall) then tx1 = cons*agn(i)*rho(i)*n0g(i) / lamg(i)**(bg+three) psacwg(i) = tx1 * qcic(i) npsacwg(i) = tx1 * ncic(i) else psacwg(i) = zero npsacwg(i) = zero end if enddo end subroutine graupel_collecting_cld_water !======================================================================== !>\ingroup micro_mg_utils_mod !! Conversion of rimed cloud water onto snow to graupel/hail subroutine graupel_riming_liquid_snow(psacws,qsic,qcic,nsic,rho,rhosn,rhog,asn,lams,n0s,dtime, & pgsacw,nscng,mgncol) integer, intent(in) :: mgncol ! Accretion of cloud water to snow tendency (modified) real(r8), dimension(mgncol), intent(inout) :: psacws real(r8), dimension(mgncol), intent(in) :: qsic ! snow mixing ratio real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid mixing ratio real(r8), dimension(mgncol), intent(in) :: nsic ! snow number concentration real(r8), dimension(mgncol), intent(in) :: rho ! air density real(r8), intent(in) :: rhosn ! snow density real(r8), intent(in) :: rhog ! graupel density real(r8), dimension(mgncol), intent(in) :: asn ! fall speed parameter for snow ! Size parameters for snow real(r8), dimension(mgncol), intent(in) :: lams real(r8), dimension(mgncol), intent(in) :: n0s real(r8), intent(in) :: dtime !Output process rates real(r8), dimension(mgncol), intent(out) :: pgsacw ! dQ graupel due to collection droplets by snow real(r8), dimension(mgncol), intent(out) :: nscng ! dN graupel due to collection droplets by snow real(r8) :: cons real(r8) :: rhosu real(r8) :: dum integer :: i !........................................................................ !Input: PSACWS,qs,qc,n0s,rho,lams,rhos,rhog !Output:PSACWS,PGSACW,NSCNG rhosu = 85000._r8/(ra * tmelt) ! typical air density at 850 mb do i=1,mgncol cons=4._r8 *2._r8 *3._r8 *rhosu*pi*ecid*ecid*gamma_2bs_plus2/(8._r8*(rhog-rhosn)) if (psacws(i).gt.0._r8 .and. qsic(i).GE.0.1e-3_r8 .AND. qcic(i).GE.0.5E-3_r8) then ! Only allow conversion if qni > 0.1 and qc > 0.5 g/kg following Rutledge and Hobbs (1984) !if (qsic(i).GE.0.1e-3_r8 .AND. qcic(i).GE.0.5E-3_r8) then ! portion of riming converted to graupel (Reisner et al. 1998, originally IS1991) ! dtime here is correct. pgsacw(i) = min(psacws(i), cons*dtime*n0s(i)*qcic(i)*qcic(i)* & asn(i)*asn(i)/ (rho(i)*lams(i)**(bs+bs+two))) ! if (pgsacw(i).lt.0_r8) then ! write(iulog,*) "pgsacw,i,lams,cons",i,pgsacw(i),lams(i),cons ! end if ! Mix rat converted into graupel as embryo (Reisner et al. 1998, orig M1990) dum= max(rhosn/(rhog-rhosn)*pgsacw(i), zero) ! Number concentraiton of embryo graupel from riming of snow nscng(i) = dum/mg0*rho(i) ! Limit max number converted to snow number (dtime here correct? We think yes) nscng(i) = min(nscng(i),nsic(i)/dtime) ! Portion of riming left for snow psacws(i) = psacws(i) - pgsacw(i) else pgsacw(i) = zero nscng(i) = zero end if enddo end subroutine graupel_riming_liquid_snow !======================================================================== !>\ingroup micro_mg_utils_mod !!CHANGE IN Q,N COLLECTION RAIN BY GRAUPEL subroutine graupel_collecting_rain(qric,qgic,umg,umr,ung,unr,rho,n0r,lamr,n0g,lamg,& pracg,npracg,mgncol) integer, intent(in) :: mgncol !MMR real(r8), dimension(mgncol), intent(in) :: qric !rain mixing ratio real(r8), dimension(mgncol), intent(in) :: qgic !graupel mixing ratio !Mass weighted Fall speeds real(r8), dimension(mgncol), intent(in) :: umg ! graupel fall speed real(r8), dimension(mgncol), intent(in) :: umr ! rain fall speed !Number weighted fall speeds real(r8), dimension(mgncol), intent(in) :: ung ! graupel fall speed real(r8), dimension(mgncol), intent(in) :: unr ! rain fall speed real(r8), dimension(mgncol), intent(in) :: rho ! air density ! Size parameters for rain real(r8), dimension(mgncol), intent(in) :: n0r real(r8), dimension(mgncol), intent(in) :: lamr ! Size parameters for graupel real(r8), dimension(mgncol), intent(in) :: n0g real(r8), dimension(mgncol), intent(in) :: lamg !Output process rates real(r8), dimension(mgncol), intent(out) :: pracg ! Q collection rain by graupel real(r8), dimension(mgncol), intent(out) :: npracg ! N collection rain by graupel ! Add collection of graupel by rain above freezing ! assume all rain collection by graupel above freezing is shed ! assume shed drops are 1 mm in size integer :: i real(r8) :: cons41 real(r8) :: cons32 real(r8) :: dum, tx1, tx2, tx3, tx4, tx5, tx6 cons41 = pi*pi*ecr*rhow cons32 = 0.5*pi*ecr do i=1,mgncol if (qric(i) >= 1.e-8_r8 .and. qgic(i) >= 1.e-8_r8) then ! pracg is mixing ratio of rain per sec collected by graupel/hail tx1 = 1.2_r8*umr(i) - 0.95_r8*umg(i) tx1 = sqrt(tx1*tx1+0.08_r8*umg(i)*umr(i)) tx2 = 1.0_r8 / lamr(i) tx3 = 1.0_r8 / lamg(i) tx4 = tx2 * tx2 tx5 = tx4 * tx4 * tx3 tx6 = rho(i) * n0r(i) * n0g(i) pracg(i) = cons41 * tx1 * tx6 * tx5 * (5.0*tx4+tx3*(tx2+tx2+0.5*tx3)) ! pracg(i) = cons41*(((1.2_r8*umr(i)-0.95_r8*umg(i))**2._r8+ & ! 0.08_r8*umg(i)*umr(i))**0.5_r8*rho(i)* & ! n0r(i)*n0g(i)/lamr(i)**3._r8* & ! (5._r8/(lamr(i)**3._r8*lamg(i))+ & ! 2._r8/(lamr(i)**2._r8*lamg(i)**2._r8)+ & ! 0.5_r8/(lamr(i)*lamg(i)**3._r8))) ! assume 1 mm drops are shed, get number shed per sec dum = pracg(i) / 5.2e-7_r8 tx1 = unr(i) - ung(i) tx1 = sqrt(1.7_r8 * tx1 * tx1 + 0.3_r8*unr(i)*ung(i)) tx4 = tx2 * tx3 npracg(i) = cons32 * tx1 * tx6 * tx4 * (tx2*(tx2+tx3)+tx3*tx3) ! npracg(i) = cons32*rho(i)*(1.7_r8*(unr(i)-ung(i))**2._r8+ & ! 0.3_r8*unr(i)*ung(i))**0.5_r8*n0r(i)*n0g(i)* & ! (1._r8/(lamr(i)**3._r8*lamg(i))+ & ! 1._r8/(lamr(i)**2._r8*lamg(i)**2._r8)+ & ! 1._r8/(lamr(i)*lamg(i)**3._r8)) ! hm 7/15/13, remove limit so that the number of collected drops can smaller than ! number of shed drops ! NPRACG(K)=MAX(NPRACG(K)-DUM,0.) npracg(i) = npracg(i) - dum else npracg(i) = zero pracg(i) = zero end if enddo end subroutine graupel_collecting_rain !======================================================================== !>\ingroup micro_mg_utils_mod !! Rain riming snow to graupel !======================================================================== ! Conversion of rimed rainwater onto snow converted to graupel subroutine graupel_rain_riming_snow(pracs,npracs,psacr,qsic,qric,nric,nsic,n0s, & lams,n0r,lamr,dtime,pgracs,ngracs,mgncol) integer, intent(in) :: mgncol ! Accretion of rain by snow real(r8), dimension(mgncol), intent(inout) :: pracs real(r8), dimension(mgncol), intent(inout) :: npracs real(r8), dimension(mgncol), intent(inout) :: psacr ! conversion due to coll of snow by rain real(r8), dimension(mgncol), intent(in) :: qsic !snow mixing ratio real(r8), dimension(mgncol), intent(in) :: qric !rain mixing ratio real(r8), dimension(mgncol), intent(in) :: nric ! rain number concentration real(r8), dimension(mgncol), intent(in) :: nsic ! snow number concentration ! Size parameters for snow real(r8), dimension(mgncol), intent(in) :: n0s real(r8), dimension(mgncol), intent(in) :: lams ! Size parameters for rain real(r8), dimension(mgncol), intent(in) :: n0r real(r8), dimension(mgncol), intent(in) :: lamr real(r8), intent(in) :: dtime !Output process rates real(r8), dimension(mgncol), intent(out) :: pgracs ! Q graupel due to collection rain by snow real(r8), dimension(mgncol), intent(out) :: ngracs ! N graupel due to collection rain by snow !Input: PRACS,NPRACS,PSACR,qni,qr,lams,lamr,nr,ns Note: No PSACR in MG2 !Output:PGRACS,NGRACS,PRACS,PSACR integer :: i real(r8) :: cons18 real(r8) :: cons19 real(r8) :: dum,fmult,tx1,tx2 cons18 = rhosn*rhosn cons19 = rhow*rhow do i=1,mgncol fmult = zero if (pracs(i) > zero .and. qsic(i) >= 0.1e-3_r8 .and. qric(i) >= 0.1e-3_r8) then ! only allow conversion if qs > 0.1 and qr > 0.1 g/kg following rutledge and hobbs (1984) !if (qsic(i).ge.0.1e-3_r8.and.qric(i).ge.0.1e-3_r8) then ! portion of collected rainwater converted to graupel (reisner et al. 1998) tx1 = four / lams(i) tx2 = four / lamr(i) tx1 = tx1 * tx1 * tx1 tx2 = tx2 * tx2 * tx2 dum = cons18 * tx1 * tx1 dum = one - max(zero, min(one, dum / (dum + cons19 * tx2 * tx2))) ! dum = cons18*(4._r8/lams(i))**3*(4._r8/lams(i))**3 & ! /(cons18*(4._r8/lams(i))**3*(4._r8/lams(i))**3+ & ! cons19*(4._r8/lamr(i))**3*(4._r8/lamr(i))**3) ! dum = min(dum,one) ! dum = max(dum, zero) ! ! pgracs(i) = (one-dum) * pracs(i) ! ngracs(i) = (one-dum) * npracs(i) pgracs(i) = dum * pracs(i) ngracs(i) = dum * npracs(i) ! limit max number converted to min of either rain or snow number concentration ngracs(i) = min(ngracs(i),nric(i)/dtime) ngracs(i) = min(ngracs(i),nsic(i)/dtime) ! amount left for snow production pracs(i) = pracs(i) - pgracs(i) npracs(i) = npracs(i) - ngracs(i) ! conversion to graupel due to collection of snow by rain ! psacr(i) = psacr(i) * (one-dum) psacr(i) = psacr(i) * dum else pgracs(i) = zero ngracs(i) = zero end if enddo end subroutine graupel_rain_riming_snow !======================================================================== ! Rime Splintering !======================================================================== !>\ingroup micro_mg_utils_mod !! Rime splintering subroutine graupel_rime_splintering(t,qcic,qric,qgic,psacwg,pracg,& qmultg,nmultg,qmultrg,nmultrg,mgncol) integer, intent(in) :: mgncol real(r8), dimension(mgncol), intent(in) :: t !temperature !MMR real(r8), dimension(mgncol), intent(in) :: qcic !liquid mixing ratio real(r8), dimension(mgncol), intent(in) :: qric !rain mixing ratio real(r8), dimension(mgncol), intent(in) :: qgic !graupel mixing ratio ! Already calculated terms for collection real(r8), dimension(mgncol), intent(inout) :: psacwg ! collection droplets by graupel real(r8), dimension(mgncol), intent(inout) :: pracg ! collection rain by graupel !Output process rates for splintering real(r8), dimension(mgncol), intent(out) :: qmultg ! Q ice mult of droplets/graupel real(r8), dimension(mgncol), intent(out) :: nmultg ! N ice mult of droplets/graupel real(r8), dimension(mgncol), intent(out) :: qmultrg ! Q ice mult of rain/graupel real(r8), dimension(mgncol), intent(out) :: nmultrg ! N ice mult of rain/graupel !Input: qg,qc,qr, PSACWG,PRACG,T !Output: NMULTG,QMULTG,NMULTRG,QMULTRG,PSACWG,PRACG integer :: i real(r8) :: fmult real(r8) :: tm_3,tm_5,tm_8 tm_3 = tmelt - 3._r8 tm_5 = tmelt - 5._r8 tm_8 = tmelt - 8._r8 !nmultg,qmultg . !======================================================================== do i=1,mgncol nmultrg(i) = zero qmultrg(i) = zero nmultg(i) = zero qmultg(i) = zero if (qgic(i) >= 0.1e-3_r8) then if (qcic(i) >= 0.5e-3_r8 .or. qric(i) >= 0.1e-3_r8) then if (psacwg(i) > zero .or. pracg(i) > zero) then if (t(i) < tm_3 .and. t(i) > tm_8) then if (t(i) > tm_3) then fmult = zero else if (t(i) <= tm_3 .and. t(i) > tm_5) then fmult = (tm_3-t(i)) * 0.5 else if (t(i) >= tm_8 .and. t(i) <= tm_5) then fmult = (t(i)-tm_8) * (one/three) else if (t(i) < tm_8) then fmult = zero end if ! 1000 is to convert from kg to g (Check if needed for MG) ! splintering from droplets accreted onto graupel if (psacwg(i) > zero) then nmultg(i) = 35.e4_r8*psacwg(i)*fmult*1000._r8 qmultg(i) = nmultg(i)*mmult ! constrain so that transfer of mass from graupel to ice cannot be more mass ! than was rimed onto graupel qmultg(i) = min(qmultg(i),psacwg(i)) psacwg(i) = psacwg(i) - qmultg(i) end if !nmultrg,qmultrg . !======================================================================== ! riming and splintering from accreted raindrops ! Factor of 1000. again (Check) if (pracg(i) > zero) then nmultrg(i) = 35.e4*pracg(i)*fmult*1000._r8 qmultrg(i) = nmultrg(i)*mmult ! constrain so that transfer of mass from graupel to ice cannot be more mass ! than was rimed onto graupel qmultrg(i) = min(qmultrg(i),pracg(i)) pracg(i) = pracg(i) - qmultrg(i) end if end if end if end if end if enddo end subroutine graupel_rime_splintering !======================================================================== ! Evaporation and sublimation of graupel !======================================================================== !MERGE WITH RAIN AND SNOW EVAP ! !subroutine graupel_sublimate_evap(t,q,qgic,rho,n0g,lamg,qvi,dv,mu,sc,bg,agn,& ! prdg,eprdg,mgncol) ! ! integer, intent(in) :: mgncol ! ! real(r8), dimension(mgncol), intent(in) :: t !temperature ! real(r8), dimension(mgncol), intent(in) :: q !specific humidity (mixing ratio) ! ! !MMR ! real(r8), dimension(mgncol), intent(in) :: qgic !graupel mixing ratio ! ! real(r8), dimension(mgncol), intent(in) :: rho ! air density ! ! ! Size parameters for graupel ! real(r8), dimension(mgncol), intent(in) :: n0g ! real(r8), dimension(mgncol), intent(in) :: lamg ! ! real(r8), dimension(mgncol), intent(in) :: qvi !saturation humidity (ice) ! ! real(r8), dimension(mgncol), intent(in) :: dv ! water vapor diffusivity ! real(r8), dimension(mgncol), intent(in) :: mu ! viscosity ! real(r8), dimension(mgncol), intent(in) :: sc ! schmidt number ! ! ! fallspeed parameters for graupel ! ! Set AGN and BG as input (in micro_mg3_0.F90) (select hail or graupel as appropriate) ! real(r8), intent(in) :: bg ! real(r8), dimension(mgncol), intent(in) :: agn ! ! ! Output tendencies (sublimation or evaporation of graupel) ! real(r8), dimension(mgncol), intent(out) :: prdg ! real(r8), dimension(mgncol), intent(out) :: eprdg ! ! real(r8) :: cons11,cons36 ! real(r8) :: abi ! real(r8) :: epsg ! integer :: i ! ! cons11=gamma(2.5_r8+bg/2._r8) !bg will be different for graupel (bg) or hail (bh) ! cons36=(2.5_r8+bg/2._r8) ! ! ! do i=1,mgncol ! ! abi = calc_ab(t(i), qvi(i), xxls) ! ! if (qgic(i).ge.qsmall) then ! epsg = 2._r8*pi*n0g(i)*rho(i)*dv(i)* & ! (f1s/(lamg(i)*lamg(i))+ & ! f2s*(agn(i)*rho(i)/mu(i))**0.5_r8* & ! sc(i)**(1._r8/3._r8)*cons11/ & ! (lamg(i)**cons36)) ! else ! epsg = 0. ! end if ! !! vapor dpeosition on graupel ! prdg(i) = epsg*(q(i)-qvi(i))/abi ! !! make sure not pushed into ice supersat/subsat !! put this in main mg3 code ... check for it ... !! formula from reisner 2 scheme !! !! dum = (qv3d(k)-qvi(k))/dt !! !! fudgef = 0.9999 !! sum_dep = prd(k)+prds(k)+mnuccd(k)+prdg(k) !! !! if( (dum.gt.0. .and. sum_dep.gt.dum*fudgef) .or. & !! (dum.lt.0. .and. sum_dep.lt.dum*fudgef) ) then !! prdg(k) = fudgef*prdg(k)*dum/sum_dep !! endif ! !! if cloud ice/snow/graupel vap deposition is neg, then assign to sublimation processes ! ! eprdg(i)=0._r8 ! ! if (prdg(i).lt.0._r8) then ! eprdg(i)=prdg(i) ! prdg(i)=0. ! end if ! ! enddo ! !end subroutine graupel_sublimate_evap !======================================================================== !UTILITIES !======================================================================== !>\ingroup micro_mg_utils_mod pure function no_limiter() real(r8) :: no_limiter no_limiter = transfer(limiter_off, no_limiter) end function no_limiter !>\ingroup micro_mg_utils_mod pure function limiter_is_on(lim) real(r8), intent(in) :: lim logical :: limiter_is_on limiter_is_on = transfer(lim, limiter_off) /= limiter_off end function limiter_is_on !>\ingroup micro_mg_utils_mod FUNCTION gamma_incomp(muice, x) real(r8) :: gamma_incomp REAL(r8), intent(in) :: muice, x REAL(r8) :: xog, kg, alfa, auxx alfa = min(max(muice+1._r8, 1._r8), 20._r8) xog = log(alfa -0.3068_r8) kg = 1.44818_r8*(alfa**0.5357_r8) auxx = max(min(kg*(log(x)-xog), 30._r8), -30._r8) gamma_incomp = max(one/(one+exp(-auxx)), 1.0e-20) END FUNCTION gamma_incomp end module micro_mg_utils