From 045993c7b9aa907e1a6b23aec12f87f13dc52015 Mon Sep 17 00:00:00 2001 From: Anastasia Popova Date: Mon, 22 Apr 2024 15:51:15 -0700 Subject: [PATCH] Add P3 collision tendencies --- Project.toml | 2 + docs/src/P3Scheme.md | 28 ++ docs/src/plots/P3SchemePlots.jl | 116 +------- src/Microphysics2M.jl | 38 ++- src/P3Scheme.jl | 456 ++++++++++++++++++++++++++++++++ test/microphysics2M_tests.jl | 8 + test/p3_tests.jl | 93 +++++++ test/performance_tests.jl | 24 +- 8 files changed, 646 insertions(+), 119 deletions(-) diff --git a/Project.toml b/Project.toml index 88f64479b..3d078e32c 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.19.0" ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" @@ -27,6 +28,7 @@ Cloudy = "0.4" DataFrames = "1.6" DocStringExtensions = "0.8, 0.9" ForwardDiff = "0.10" +HCubature = "1.6" MLJ = "0.20" QuadGK = "2.9" RootSolvers = "0.3, 0.4" diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md index a197f080f..1fdf40d30 100644 --- a/docs/src/P3Scheme.md +++ b/docs/src/P3Scheme.md @@ -188,3 +188,31 @@ They can be compared with Figure 2 from [MorrisonMilbrandt2015](@cite). include("plots/P3TerminalVelocityPlots.jl") ``` ![](MorrisonandMilbrandtFig2.svg) + +## Collisions + +Collisions are calculated through the following equation: + +```math +\frac{dq_c}{dt} = \int_{0}^{\infty} \! \int_{0}^{\infty} \! \frac{E_ci}{\rho} N'_i(D_p) N'_c(D_c) A(D_i, D_c) m_{collected}(D_{collected}) |V_i(D_i) - V_c(D_c)| \mathrm{d}D_i \mathrm{d}D_c +``` + where the values are defined as follows: + - subscript ``i`` - ice particles (either cloud ice or precipitating ice) + - subscript ``c`` - corresponding colliding species + - subscript ``collected`` - whichever species is being collected by the collisions + - ``E_ci`` - collision efficiency between particles + - ``N'_x`` - distribution of particles of x type + - ``A`` - collision kernel between ice and colliding species + - ``m_x(D)`` - mass of given species x at diameter D + - ``V_x(D)`` - Chen terminal velocity of given species x at diameter D + + In our parametrization we take ``E_ci = 1`` and ``A(D_i, D_c) = \pi \frac{(D_i + D_c)^2}{4}``. We also assume that if the temperature is above freezing then ice particles get collected, otherwise the colliding species is collected by the ice. + + Furthermore, the distributions of non-ice species are taken to be of the following forms: + - cloud water: ``N_0 D^8 e^{-\lambda D^3}`` + - rain : ``N_0 e^{-\lambda D}`` + + ``N_0`` and ``\lambda`` can be solved for respectively in each scheme through setting: + - ``N_c = \int_{0}^{\infty} N'_c(D) \mathrm{d}D`` + - ``q_c = \int_{0}^{\infty} m_c(D) N'_c(D) \mathrm{d}D`` + where ``m_c(D)`` is the mass of particle of size D. Rain and cloud water particles are assumed to be spherical therefore, ``m_c(D) = \pi \rho_w \frac{D^3}{6}`` for both (``\rho_w`` = density of water). \ No newline at end of file diff --git a/docs/src/plots/P3SchemePlots.jl b/docs/src/plots/P3SchemePlots.jl index c106b8f99..b3bc2ecee 100644 --- a/docs/src/plots/P3SchemePlots.jl +++ b/docs/src/plots/P3SchemePlots.jl @@ -9,99 +9,7 @@ FT = Float64 const PSP3 = CMP.ParametersP3 p3 = CMP.ParametersP3(FT) -""" - mass_(p3, D, ρ, F_r) - - p3 - a struct with P3 scheme parameters - - D - maximum particle dimension [m] - - ρ - bulk ice density (ρ_i for small ice, ρ_g for graupel) [kg/m3] - - F_r - rime mass fraction [q_rim/q_i] - -Returns mass as a function of size for differen particle regimes [kg] -""" -# for spherical ice (small ice or completely rimed ice) -mass_s(D::FT, ρ::FT) where {FT <: Real} = FT(π) / 6 * ρ * D^3 -# for large nonspherical ice (used for unrimed and dense types) -mass_nl(p3::PSP3, D::FT) where {FT <: Real} = P3.α_va_si(p3) * D^p3.β_va -# for partially rimed ice -mass_r(p3::PSP3, D::FT, F_r::FT) where {FT <: Real} = - P3.α_va_si(p3) / (1 - F_r) * D^p3.β_va - -""" - p3_mass(p3, D, F_r, th) - - - p3 - a struct with P3 scheme parameters - - D - maximum particle dimension - - F_r - rime mass fraction (q_rim/q_i) - - th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d) - -Returns mass(D) regime, used to create figures for the docs page. -""" -function p3_mass( - p3::PSP3, - D::FT, - F_r::FT, - th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)), -) where {FT <: Real} - D_th = P3.D_th_helper(p3) - if D_th > D - return mass_s(D, p3.ρ_i) # small spherical ice - elseif F_r == 0 - return mass_nl(p3, D) # large nonspherical unrimed ice - elseif th.D_gr > D >= D_th - return mass_nl(p3, D) # dense nonspherical ice - elseif th.D_cr > D >= th.D_gr - return mass_s(D, th.ρ_g) # graupel - elseif D >= th.D_cr - return mass_r(p3, D, F_r) # partially rimed ice - end -end - -""" - A_(p3, D) - - - p3 - a struct with P3 scheme parameters - - D - maximum particle dimension - -Returns particle projected area as a function of size for different particle regimes -""" -# for spherical particles -A_s(D::FT) = FT(π) / 4 * D^2 -# for nonspherical particles -A_ns(p3::PSP3, D::FT) = p3.γ * D^p3.σ -# partially rimed ice -A_r(p3::PSP3, F_r::FT, D::FT) = F_r * A_s(D) + (1 - F_r) * A_ns(p3, D) - -""" - area(p3, D, F_r, th) - - - p3 - a struct with P3 scheme parameters - - D - maximum particle dimension - - F_r - rime mass fraction (q_rim/q_i) - - th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d) - -Returns area(D), used to create figures for the documentation. -""" -function area( - p3::PSP3, - D::FT, - F_r::FT, - th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)), -) where {FT <: Real} - # Area regime: - if P3.D_th_helper(p3) > D - return A_s(D) # small spherical ice - elseif F_r == 0 - return A_ns(p3, D) # large nonspherical unrimed ice - elseif th.D_gr > D >= P3.D_th_helper(p3) - return A_ns(p3, D) # dense nonspherical ice - elseif th.D_cr > D >= th.D_gr - return A_s(D) # graupel - elseif D >= th.D_cr - return A_r(p3, F_r, D) # partially rimed ice - end - -end function define_axis(fig, row_range, col_range, title, ylabel, yticks, aspect) return CMK.Axis( @@ -143,13 +51,13 @@ function p3_relations_plot() sol4_5 = P3.thresholds(p3, 400.0, 0.5) sol4_8 = P3.thresholds(p3, 400.0, 0.8) # m(D) - fig1_0 = CMK.lines!(ax1, D_range * 1e3, [p3_mass(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw) - fig1_5 = CMK.lines!(ax1, D_range * 1e3, [p3_mass(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw) - fig1_8 = CMK.lines!(ax1, D_range * 1e3, [p3_mass(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw) + fig1_0 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw) + fig1_5 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw) + fig1_8 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw) # a(D) - fig2_0 = CMK.lines!(ax2, D_range * 1e3, [area(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw) - fig2_5 = CMK.lines!(ax2, D_range * 1e3, [area(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw) - fig2_8 = CMK.lines!(ax2, D_range * 1e3, [area(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw) + fig2_0 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw) + fig2_5 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw) + fig2_8 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw) # plot verical lines for ax in [ax1, ax2] global d_tha = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw) @@ -164,13 +72,13 @@ function p3_relations_plot() sol_4 = P3.thresholds(p3, 400.0, 0.95) sol_8 = P3.thresholds(p3, 800.0, 0.95) # m(D) - fig3_200 = CMK.lines!(ax3, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw) - fig3_400 = CMK.lines!(ax3, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw) - fig3_800 = CMK.lines!(ax3, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw) + fig3_200 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw) + fig3_400 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw) + fig3_800 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw) # a(D) - fig3_200 = CMK.lines!(ax4, D_range * 1e3, [area(p3, D, 0.5, sol_2) for D in D_range], color = cl[1], linewidth = lw) - fig3_400 = CMK.lines!(ax4, D_range * 1e3, [area(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw) - fig3_800 = CMK.lines!(ax4, D_range * 1e3, [area(p3, D, 0.5, sol_8) for D in D_range], color = cl[3], linewidth = lw) + fig3_200 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_2) for D in D_range], color = cl[1], linewidth = lw) + fig3_400 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw) + fig3_800 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_8) for D in D_range], color = cl[3], linewidth = lw) # plot verical lines for ax in [ax3, ax4] global d_thb = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw) diff --git a/src/Microphysics2M.jl b/src/Microphysics2M.jl index 1e2ec3b3b..703ceb0fb 100644 --- a/src/Microphysics2M.jl +++ b/src/Microphysics2M.jl @@ -448,7 +448,7 @@ function rain_evaporation( end """ - radar_reflectivity(struct, q_liq, q_rai, N_liq, N_rai, ρ_air, ρ_w) + radar_reflectivity(struct, q_liq, q_rai, N_liq, N_rai, ρ_air, ρ_w, τ_q, τ_N) - `struct` - type for 2-moment rain autoconversion parameterization - `q_liq` - cloud water specific humidity @@ -457,6 +457,8 @@ end - `N_rai` - rain droplet number density - `ρ_air` - air density - `ρ_w` - water density + - `τ_q` - threshold for minimum specific humidity value + - `τ_N` - threshold for minimum number density value Returns logarithmic radar reflectivity from the assumed cloud and rain particle size distribuions normalized by the reflectivty of 1 millimiter drop in a volume of @@ -470,6 +472,8 @@ function radar_reflectivity( N_rai::FT, ρ_air::FT, ρ_w::FT, + τ_q::FT, + τ_N::FT, ) where {FT} # we assume a cloud droplets gamma distribution, @@ -482,14 +486,16 @@ function radar_reflectivity( νr = FT(-2 / 3) μr = FT(1 / 3) - #change of units for better accuracy - ρ_air *= FT(1e-3) - ρ_w *= FT(1e-3) + # change of units for N_liq and N_rai + # from m^-3 to mm^-3 N_liq *= FT(1e-9) N_rai *= FT(1e-9) - xc = (N_liq == 0) ? FT(0) : ((q_liq * ρ_air) / N_liq) - xr = (N_rai == 0) ? FT(0) : ((q_rai * ρ_air) / N_rai) + q_liq = (q_liq < τ_q) ? FT(0) : q_liq + q_rai = (q_rai < τ_q) ? FT(0) : q_rai + + xc = (N_liq < τ_N) ? FT(0) : ((q_liq * ρ_air) / N_liq) + xr = (N_rai < τ_N) ? FT(0) : ((q_rai * ρ_air) / N_rai) C = FT(4 / 3 * π * ρ_w) Z₀ = FT(1e-18) @@ -505,12 +511,12 @@ function radar_reflectivity( Zc = (Bc == 0) ? FT(0) : (FT(24) * Ac / (Bc^FT(5) * C^FT(2))) Zr = (Br == 0) ? FT(0) : (FT(2160) * Ar / (Br^FT(7) * C^FT(2))) - return ((Zc + Zr) == 0) ? FT(-200) : + return ((Zc + Zr) == 0) ? FT(-135) : FT(10) * (log10((Zc + Zr) / Z₀) + log10(FT(1e-9))) end """ - effective_radius(struct, q_liq, q_rai, N_liq, N_rai, ρ_air, ρ_w) + effective_radius(struct, q_liq, q_rai, N_liq, N_rai, ρ_air, ρ_w, τ_q, τ_N) - `struct` - type for 2-moment rain autoconversion parameterization - `q_liq` - cloud water specific humidity @@ -519,6 +525,8 @@ end - `N_rai` - rain droplet number density - `ρ_air` - air density - `ρ_w` - water density + - `τ_q` - threshold for minimum specific humidity value + - `τ_N` - threshold for minimum number density value Returns effective radius using the 2-moment scheme cloud and rain particle size distributions @@ -531,6 +539,8 @@ function effective_radius( N_rai::FT, ρ_air::FT, ρ_w::FT, + τ_q::FT, + τ_N::FT, ) where {FT} # we assume a cloud droplets gamma distribution, @@ -543,14 +553,16 @@ function effective_radius( νr = FT(-2 / 3) μr = FT(1 / 3) - #change of units for better accuracy - ρ_air *= FT(1e-3) - ρ_w *= FT(1e-3) + # change of units for N_liq and N_rai + # from m^-3 to mm^-3 N_liq *= FT(1e-9) N_rai *= FT(1e-9) - xc = (N_liq == 0) ? FT(0) : ((q_liq * ρ_air) / N_liq) - xr = (N_rai == 0) ? FT(0) : ((q_rai * ρ_air) / N_rai) + q_liq = (q_liq < τ_q) ? FT(0) : q_liq + q_rai = (q_rai < τ_q) ? FT(0) : q_rai + + xc = (N_liq < τ_N) ? FT(0) : ((q_liq * ρ_air) / N_liq) + xr = (N_rai < τ_N) ? FT(0) : ((q_rai * ρ_air) / N_rai) C = FT((4 / 3) * π * ρ_w) Bc = diff --git a/src/P3Scheme.jl b/src/P3Scheme.jl index 6992b0b82..db57dff9b 100644 --- a/src/P3Scheme.jl +++ b/src/P3Scheme.jl @@ -14,6 +14,7 @@ module P3Scheme import SpecialFunctions as SF import QuadGK as QGK import RootSolvers as RS +import HCubature as HC import ClimaParams as CP import CloudMicrophysics.Parameters as CMP @@ -156,6 +157,103 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT} end end +""" + mass_(p3, D, ρ, F_r) + + - p3 - a struct with P3 scheme parameters + - D - maximum particle dimension [m] + - ρ - bulk ice density (ρ_i for small ice, ρ_g for graupel) [kg/m3] + - F_r - rime mass fraction [q_rim/q_i] + +Returns mass as a function of size for differen particle regimes [kg] +""" +# for spherical ice (small ice or completely rimed ice) +mass_s(D::FT, ρ::FT) where {FT <: Real} = FT(π) / 6 * ρ * D^3 +# for large nonspherical ice (used for unrimed and dense types) +mass_nl(p3::PSP3, D::FT) where {FT <: Real} = α_va_si(p3) * D^p3.β_va +# for partially rimed ice +mass_r(p3::PSP3, D::FT, F_r::FT) where {FT <: Real} = + α_va_si(p3) / (1 - F_r) * D^p3.β_va + +""" + p3_mass(p3, D, F_r, th) + + - p3 - a struct with P3 scheme parameters + - D - maximum particle dimension + - F_r - rime mass fraction (q_rim/q_i) + - th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d) + +Returns mass(D) regime, used to create figures for the docs page. +""" +function p3_mass( + p3::PSP3, + D::FT, + F_r::FT, + th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)), +) where {FT <: Real} + D_th = D_th_helper(p3) + if D_th > D + return mass_s(D, p3.ρ_i) # small spherical ice + elseif F_r == 0 + return mass_nl(p3, D) # large nonspherical unrimed ice + elseif th.D_gr > D >= D_th + return mass_nl(p3, D) # dense nonspherical ice + elseif th.D_cr > D >= th.D_gr + return mass_s(D, th.ρ_g) # graupel + elseif D >= th.D_cr + return mass_r(p3, D, F_r) # partially rimed ice + end +end + +""" + A_(p3, D) + + - p3 - a struct with P3 scheme parameters + - D - maximum particle dimension + +Returns particle projected area as a function of size for different particle regimes +""" +# for spherical particles +A_s(D::FT) where {FT <: Real} = FT(π) / 4 * D^2 +# for nonspherical particles +A_ns(p3::PSP3, D::FT) where {FT <: Real} = p3.γ * D^p3.σ +# partially rimed ice +A_r(p3::PSP3, F_r::FT, D::FT) where {FT <: Real} = + F_r * A_s(D) + (1 - F_r) * A_ns(p3, D) + +""" + p3_area(p3, D, F_r, th) + + - p3 - a struct with P3 scheme parameters + - D - maximum particle dimension + - F_r - rime mass fraction (q_rim/q_i) + - th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d) + +Returns area(D), used to create figures for the documentation. +""" +function p3_area( + p3::PSP3, + D::FT, + F_r::FT, + th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)), +) where {FT <: Real} + # Area regime: + if D_th_helper(p3) > D + return A_s(D) # small spherical ice + elseif F_r == 0 + return A_ns(p3, D) # large nonspherical unrimed ice + elseif th.D_gr > D >= D_th_helper(p3) + return A_ns(p3, D) # dense nonspherical ice + elseif th.D_cr > D >= th.D_gr + return A_s(D) # graupel + elseif D >= th.D_cr + return A_r(p3, F_r, D) # partially rimed ice + else + println(D) + throw("D not in range") + end +end + # Some wrappers to cast types from SF.gamma # (which returns Float64 even when the input is Float32) Γ(a::FT, z::FT) where {FT <: Real} = FT(SF.gamma(a, z)) @@ -667,4 +765,362 @@ function D_m(p3::PSP3, q::FT, N::FT, ρ_r::FT, F_r::FT) where {FT} return n / q end +""" + get_rain_parameters(q, N) + + - q - mass mixing ratio of rain + - N - number mixing ratio of rain + + returns the parameters λ and N_0 where N' = N_0 * exp(-λ * D) +""" +function get_rain_parameters(q::FT, N::FT) where {FT} + ρ_w = FT(1000) + λ = (π * ρ_w * N / q)^(1 / 3) + N_0 = N * λ + return (λ, N_0) +end + +""" + get_cloud_parameters(q, N) + + - q - mass mixing ratio of rain + - N - number mixing ratio of rain + + returns the parameters λ and N_0 where N' = N_0 * D ^ 8 * exp(-λ * D ^ 3) +""" +function get_cloud_parameters(q::FT, N::FT) where {FT} + ρ_w = FT(1000) + λ = π * ρ_w * N / (2 * q) + N_0 = 3 / 2 * N * λ^3 + return (λ, N_0) +end + +""" + N′_rain(D, q, N) + + - D - diameter of particle + - q - mass mixing ratio + - N - number mixing ratio + + Returns the distribution of rain particles (assumed to be of the form + N'(D) = N0 * exp(-λD)) at given D +""" +function N′rain(D::FT, q::FT, N::FT) where {FT} + (λ, N_0) = get_rain_parameters(q, N) + return N_0 * exp(-λ * D) +end + +""" + N′_cloud(D, q, N) + + - D - diameter of particle + - q - mass mixing ratio + - N - number mixing ratio + + Returns the distribution of cloud particles (assumed to be of the form + N'(D) = N0 * D ^ 8 * exp(-λ D ^ 3)) at given D +""" +function N′cloud(D::FT, q::FT, N::FT) where {FT} + (λ, N_0) = get_cloud_parameters(q, N) + return N_0 * D^8 * exp(-λ * D^3) +end + +""" + N′_ice(D, p3, λ, N0) + + - D - diameter of particle + - p3 - a struct containing P3 scheme parameters + - λ - shape parameter of distribution + - N0 - shape parameter of distribution + + Returns the distribution of ice particles (assumed to be of the form + N'(D) = N0 * D ^ μ * exp(-λD)) at given D +""" +function N′ice(D::FT, p3::PSP3, λ::FT, N_0::FT) where {FT} + return N_0 * D^DSD_μ(p3, λ) * exp(-λ * D) +end + +""" + ϕ_ice(p3, D, F_r, th) + + - p3 - a struct containing P3 Scheme parameters + - D - maximum particle dimension + - F_r - rime mass fraction (q_rim/q_i) + - th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d) + + Returns the ϕ calculation for an ice particle in the P3 Scheme with given + maximum particle dimension +""" +function ϕ_ice( + D::FT, + p3::PSP3, + F_r::FT, + th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)), +) where {FT} + return 16 * p3.ρ_i^2 * (p3_area(p3, D, F_r, th))^3 / + (9 * π * (p3_mass(p3, D, F_r, th))^2) +end + +""" + velocity_chen(p3, Chen2022, ρ_a, D, F_r, th) + + - p3 - a struct containing P3 Scheme parameters + - Chen2022 - a struct with terminal velocity parameters as in Chen(2022) + - ρ_a - density of air + - D - maximum particle dimension + - F_r - rime mass fraction (q_rim/q_i) + - th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d) + + Returns the terminal velocity of ice at given particle dimension using P3 + distributions and Chen 2022 parametrizations +""" +function velocity_chen( + D::FT, + p3::PSP3, + Chen2022::CMP.Chen2022VelTypeSnowIce, + ρ_a::FT, + F_r::FT, + th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)), +) where {FT} + if D <= FT(0.000625) + (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) + else + (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) + end + + κ = FT(-1 / 6) + + v = 0 + for i in 1:2 + v += ai[i] * D^(bi[i]) * exp(-ci[i] * D) + end + + return ϕ_ice(D, p3, F_r, th)^κ * v +end + +""" + velocity_chen(p3, Chen2022, ρ_a, D) + + - Chen2022 - a struct with terminal velocity parameters as in Chen(2022) + - ρ_a - density of air + - D - maximum particle dimension + + Returns the terminal velocity of rain given Chen 2022 velocity parametrizations +""" +function velocity_chen( + D::FT, + Chen2022::CMP.Chen2022VelTypeRain, + ρ_a::FT, +) where {FT} + (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) + + v = 0 + for i in 1:3 + v += ai[i] * D^(bi[i]) * exp(-ci[i] * D) + end + + return v +end + +""" + integration_bounds(p3, tolerance, λ, N_0, Nᵢ, qᵣ, Nᵣ) + + - p3 - a struct containing P3 Scheme parameters + - tolerance - tolerance to which distributions need to be evaluated to + - λ - shape parameter of ice distribution + - N_0 - intercept size distribution of ice + - Nᵢ - number mixing ratio of ice + - q - mass mizing ratio of colliding species + - N - number mixing ratio of colliding species + + Returns the bounds over which to integrate rain, cloud, and ice distributions to ensure + coverage or more than (1 - tolerance) of each distribution +""" +function integration_bounds( + type::String, + p3::PSP3, + tolerance::FT, + λ::FT, + N_0::FT, + Nᵢ::FT, + q::FT, + N::FT, +) where {FT} + if type == "rain" + rain_λ, = get_rain_parameters(q, N) + colliding_x = -1 / rain_λ * log(tolerance) + elseif type == "cloud" + cloud_λ, = get_cloud_parameters(q, N) + cloud_problem(x) = + tolerance - + exp(-exp(x)^3 * cloud_λ) * + (1 + exp(x)^3 * cloud_λ + 1 / 2 * exp(x)^6 * cloud_λ^2) + guess = + log(0.5) + + (log(0.00025) - log(0.5)) / (log(1e12) - log(1e2)) * + (log(cloud_λ) - log(10^2)) + log_cloud_x = + RS.find_zero( + cloud_problem, + RS.NewtonsMethodAD(guess), + RS.CompactSolution(), + RS.RelativeSolutionTolerance(eps(FT)), + 5, + ).root + colliding_x = exp(log_cloud_x) + end + + ice_problem(x) = + tolerance - Γ(1 + DSD_μ(p3, λ), FT(exp(x)) * λ) / Γ(1 + DSD_μ(p3, λ)) + guess = log(19 / 6 * (DSD_μ(p3, λ) - 1) + 39) - log(λ) + log_ice_x = + RS.find_zero( + ice_problem, + RS.SecantMethod(guess - 1, guess), + RS.CompactSolution(), + RS.RelativeSolutionTolerance(eps(FT)), + 5, + ).root + + return (2 * colliding_x, 2 * exp(log_ice_x)) +end + +""" + a(D1, D2) + + - D1 - maximum dimension of first particle + - D2 - maximum dimension of second particle + + Returns the collision kernel (assumed to be of the form π(r1 + r2)^2) for the + two colliding particles +""" +function a(D1::FT, D2::FT) where {FT} + # TODO make this more accurate for non-spherical particles + return π * (D1 / 2 + D2 / 2)^2 +end + +""" + N′colliding(type, D, q, N) + + - type - defines what is colliding with ice ("rain" or "cloud") + - D - maximum dimension + - q - mass mixing ratio + - N - number mixing ratio + + Returns the corresponding distribution of particles depending on type +""" +function N′colliding(type::String, D::FT, q::FT, N::FT) where {FT} + if type == "rain" + return N′rain(D, q, N) + elseif type == "cloud" + return N′cloud(D, q, N) + else + throw("Undefined Type of Collisions") + end +end + +""" + vel_diff(type, D, Dᵢ, p3, Chen2022, ρ_a, F_r, th) + + - type - defines what is colliding with ice ("rain" or "cloud") + - D - maximum dimension of colliding particle + - Dᵢ - maximum dimension of ice particle + - p3 - a struct containing P3 parameters + - Chen2022 - a struct containing Chen 2022 velocity parameters + - ρ_a - density of air + - F_r - rime mass fraction (q_rim/ q_i) + - th - thresholds as calculated by thresholds() + + Returns the corresponding velocity difference of colliding particles depending on type +""" +function vel_diff( + type::String, + D::FT, + Dᵢ::FT, + p3::PSP3, + Chen2022::CMP.Chen2022VelType, + ρ_a::FT, + F_r::FT, + th, +) where {FT} + if type == "rain" + return abs( + velocity_chen(Dᵢ, p3, Chen2022.snow_ice, ρ_a, F_r, th) - + velocity_chen(D, Chen2022.rain, ρ_a), + ) + elseif type == "cloud" + return abs(velocity_chen(Dᵢ, p3, Chen2022.snow_ice, ρ_a, F_r, th)) + end +end + +""" + ice_collisions(type, p3, Chen2022, ρ_a, F_r, qᵣ, qᵢ, Nᵣ, Nᵢ, ρ, ρ_r, E_ri) + + - type - defines what is colliding with the ice ("cloud" or "rain") + - p3 - a struct with P3 scheme parameters + - Chen2022 - a struct with terminal velocity parameters as in Chen (2022) + - qᵢ - mass mixing ratio of ice + - Nᵢ - number mixing ratio of ice + - q_c - mass mixing ratio of colliding species + - N_c - number mixing ratio of colliding species + - ρ_a - density of air + - F_r - rime mass fraction (q_rim/ q_i) + - ρ_r - rime density (q_rim/B_rim) + - T - temperature (in K) + - E_ci - collision efficiency between ice and colliding species + + Returns the rate of collisions between cloud ice and rain + Equivalent to the measure of QRCOL in Morrison and Mildbrandt (2015) +""" +function ice_collisions( + type::String, + p3::PSP3, + Chen2022::CMP.Chen2022VelType, + qᵢ::FT, + Nᵢ::FT, + q_c::FT, + N_c::FT, + ρ_a::FT, + F_r::FT, + ρ_r::FT, + T::FT, + E_ci = FT(1), +) where {FT} + ρ_w = FT(1000) # TO DO: add this as density_liquid_water from ClimaParams + th = thresholds(p3, ρ_r, F_r) + (λ, N_0) = distribution_parameter_solver(p3, qᵢ, Nᵢ, ρ_r, F_r) + (colliding_bound, ice_bound) = + integration_bounds(type, p3, eps(FT), λ, N_0, Nᵢ, q_c, N_c) + + if T > 273.15 + f_warm(D_c, Dᵢ) = + E_ci / ρ_a * + N′colliding(type, D_c, q_c, N_c) * + N′ice(Dᵢ, p3, λ, N_0) * + a(D_c, Dᵢ) * + p3_mass(p3, Dᵢ, F_r, th) * + vel_diff(type, D_c, Dᵢ, p3, Chen2022, ρ_a, F_r, th) + (dqdt, error) = HC.hcubature( + d -> f_warm(d[1], d[2]), + (FT(0), FT(0)), + (colliding_bound, ice_bound), + ) + else + f_cold(D_c, Dᵢ) = + E_ci / ρ_a * + N′colliding(type, D_c, q_c, N_c) * + N′ice(Dᵢ, p3, λ, N_0) * + a(D_c, Dᵢ) * + mass_s(D_c, ρ_w) * + vel_diff(type, D_c, Dᵢ, p3, Chen2022, ρ_a, F_r, th) + (dqdt, error) = HC.hcubature( + d -> f_cold(d[1], d[2]), + (FT(0), FT(0)), + (colliding_bound, ice_bound), + ) + end + + return dqdt +end + end diff --git a/test/microphysics2M_tests.jl b/test/microphysics2M_tests.jl index 31ef002c6..9c3c3b718 100644 --- a/test/microphysics2M_tests.jl +++ b/test/microphysics2M_tests.jl @@ -495,6 +495,8 @@ function test_microphysics2M(FT) N_liq = FT(15053529) q_rai = FT(1.573e-4) N_rai = FT(510859) + τ_q = FT(1e-12) + τ_N = FT(1e-18) #action rr = CM2.radar_reflectivity( @@ -505,6 +507,8 @@ function test_microphysics2M(FT) N_rai, ρ_air, ρ_w, + τ_q, + τ_N, ) TT.@test rr ≈ FT(-13) atol = 2 @@ -519,6 +523,8 @@ function test_microphysics2M(FT) N_liq = FT(15053529) q_rai = FT(1.573e-4) N_rai = FT(510859) + τ_q = FT(1e-12) + τ_N = FT(1e-18) #action reff = CM2.effective_radius( @@ -529,6 +535,8 @@ function test_microphysics2M(FT) N_rai, ρ_air, ρ_w, + τ_q, + τ_N, ) #test diff --git a/test/p3_tests.jl b/test/p3_tests.jl index e2dab5a4f..f48f241f5 100644 --- a/test/p3_tests.jl +++ b/test/p3_tests.jl @@ -205,6 +205,98 @@ function test_velocities(FT) end end +function test_tendencies(FT) + + p3 = CMP.ParametersP3(FT) + Chen2022 = CMP.Chen2022VelType(FT) + + TT.@testset "Collision Tendencies Smoke Test" begin + N = FT(1e8) + ρ_a = FT(1.2) + ρ_r = FT(500) + F_r = FT(0.5) + T_warm = FT(300) + T_cold = FT(200) + + qs = range(0.001, stop = 0.005, length = 5) + q_const = FT(0.05) + + cloud_expected_warm = + [5.78e-5, 0.00019256, 0.00039239, 0.00065184, 0.00096698] + cloud_expected_cold = + [0.0016687, 0.0026921, 0.0035912, 0.0044255, 0.00522] + rain_expected_warm = [0.0003392, 0.000713, 0.001103, 0.001506, 0.00192] + rain_expected_cold = [0.2905, 0.2982, 0.3033, 0.3072, 0.3104] + + for i in axes(qs, 1) + cloud_warm = P3.ice_collisions( + "cloud", + p3, + Chen2022, + qs[i], + N, + q_const, + N, + ρ_a, + F_r, + ρ_r, + T_warm, + ) + cloud_cold = P3.ice_collisions( + "cloud", + p3, + Chen2022, + qs[i], + N, + q_const, + N, + ρ_a, + F_r, + ρ_r, + T_cold, + ) + + TT.@test cloud_warm >= 0 + TT.@test cloud_warm ≈ cloud_expected_warm[i] rtol = 1e-3 + TT.@test cloud_cold >= 0 + TT.@test cloud_cold ≈ cloud_expected_cold[i] rtol = 1e-3 + + rain_warm = P3.ice_collisions( + "rain", + p3, + Chen2022, + qs[i], + N, + q_const, + N, + ρ_a, + F_r, + ρ_r, + T_warm, + ) + rain_cold = P3.ice_collisions( + "rain", + p3, + Chen2022, + qs[i], + N, + q_const, + N, + ρ_a, + F_r, + ρ_r, + T_cold, + ) + + TT.@test rain_warm >= 0 + TT.@test rain_warm ≈ rain_expected_warm[i] rtol = 1e-3 + TT.@test rain_cold >= 0 + TT.@test rain_cold ≈ rain_expected_cold[i] rtol = 1e-3 + end + end + +end + println("Testing Float32") test_p3_thresholds(Float32) #TODO - only works for Float64 now. We should switch the units inside the solver @@ -215,3 +307,4 @@ println("Testing Float64") test_p3_thresholds(Float64) test_p3_shape_solver(Float64) test_velocities(Float64) +test_tendencies(Float64) diff --git a/test/performance_tests.jl b/test/performance_tests.jl index b50b18436..a1aa93a21 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -225,12 +225,32 @@ function benchmark_test(FT) ) bench_press( CM2.radar_reflectivity, - (sb2006.acnv, q_liq, q_rai, N_liq, N_rai, FT(1), FT(1000)), + ( + sb2006.acnv, + q_liq, + q_rai, + N_liq, + N_rai, + FT(1), + FT(1000), + FT(1e-12), + FT(1e-18), + ), 800, ) bench_press( CM2.effective_radius, - (sb2006.acnv, q_liq, q_rai, N_liq, N_rai, FT(1), FT(1000)), + ( + sb2006.acnv, + q_liq, + q_rai, + N_liq, + N_rai, + FT(1), + FT(1000), + FT(1e-12), + FT(1e-18), + ), 800, ) bench_press(