P3 Terminal Velocity

Project.toml

@@ -7,6 +7,7 @@ version = "0.18.0"
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
@@ -29,6 +30,7 @@ DocStringExtensions = "0.8, 0.9"
 EnsembleKalmanProcesses = "1.1.5"
 ForwardDiff = "0.10"
 MLJ = "0.20"
+QuadGK = "2.9.4"
 RootSolvers = "0.3, 0.4"
 SpecialFunctions = "1, 2"
 Thermodynamics = "0.12.4"

docs/Project.toml

@@ -9,6 +9,7 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"

docs/src/API.md

@@ -67,7 +67,6 @@ MicrophysicsFlexible.weighted_vt
 ```@docs
 P3Scheme
 P3Scheme.thresholds
-P3Scheme.p3_mass
 P3Scheme.distribution_parameter_solver
 ```
 

docs/src/plots/P3SchemePlots.jl

@@ -9,6 +9,54 @@ 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)
 
@@ -43,17 +91,13 @@ function area(
     # Area regime:
     if P3.D_th_helper(p3) > D
         return A_s(D)                      # small spherical ice
-    end
-    if F_r == 0
+    elseif F_r == 0
         return A_ns(p3, D)                 # large nonspherical unrimed ice
-    end
-    if th.D_gr > D >= P3.D_th_helper(p3)
+    elseif th.D_gr > D >= P3.D_th_helper(p3)
         return A_ns(p3, D)                 # dense nonspherical ice
-    end
-    if th.D_cr > D >= th.D_gr
+    elseif th.D_cr > D >= th.D_gr
         return A_s(D)                      # graupel
-    end
-    if D >= th.D_cr
+    elseif D >= th.D_cr
         return A_r(p3, F_r, D)             # partially rimed ice
     end
 
@@ -76,6 +120,69 @@ function define_axis(fig, row_range, col_range, title, ylabel, yticks, aspect)
         aspect = aspect,
         limits = ((0.02, 10.0), nothing),
     )
+    sol_5 = P3.thresholds(p3, 400.0, 0.5)
+    sol_8 = P3.thresholds(p3, 400.0, 0.8)
+
+    #! format: off
+    fig1_a_0 = CMK.lines!(ax1_a, D_range * 1e3, [p3_mass(p3, D, 0.0       ) for D in D_range], color = cl[1], linewidth = lw)
+    fig1_a_5 = CMK.lines!(ax1_a, D_range * 1e3, [p3_mass(p3, D, 0.5, sol_5) for D in D_range], color = cl[2], linewidth = lw)
+    fig1_a_8 = CMK.lines!(ax1_a, D_range * 1e3, [p3_mass(p3, D, 0.8, sol_8) for D in D_range], color = cl[3], linewidth = lw)
+
+    d_tha  = CMK.vlines!(ax1_a, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
+    d_cr_5 = CMK.vlines!(ax1_a, sol_5[1]           * 1e3, linestyle = :dot,  color = cl[2], linewidth = lw)
+    d_cr_8 = CMK.vlines!(ax1_a, sol_8[1]           * 1e3, linestyle = :dot,  color = cl[3], linewidth = lw)
+    d_gr_5 = CMK.vlines!(ax1_a, sol_5[2]           * 1e3, linestyle = :dash, color = cl[2], linewidth = lw)
+    d_gr_8 = CMK.vlines!(ax1_a, sol_8[2]           * 1e3, linestyle = :dash, color = cl[3], linewidth = lw)
+
+    leg1_a = CMK.Legend(fig1_a[8:9, 1], [fig1_a_0, fig1_a_5, fig1_a_8], [CMK.L"$F_{r} = 0.0$", CMK.L"$F_{r} = 0.5$", CMK.L"$F_{r} = 0.8$"], framevisible = false)
+    leg1_a_dth = CMK.Legend(fig1_a[8:9, 3], [d_tha], [CMK.L"$D_{th}$"], framevisible = false)
+    leg1_a_dcr = CMK.Legend(fig1_a[8:9, 7], [d_cr_5, d_cr_8], [CMK.L"$D_{cr}$ for $F_{r} = 0.5$", CMK.L"$D_{cr}$ for $F_{r} = 0.8$"], framevisible = false)
+    leg1_a_dgr = CMK.Legend(fig1_a[8:9, 5], [d_gr_5, d_gr_8], [CMK.L"$D_{gr}$ for $F_{r} = 0.5$", CMK.L"$D_{gr}$ for $F_{r} = 0.8$"], framevisible = false)
+
+    #! format: on
+    CMK.save("MorrisonandMilbrandtFig1a.svg", fig1_a)
+
+    fig1_b = CMK.Figure()
+    ax1_b = CMK.Axis(
+        fig1_b[1:10, 1:11],
+        title = CMK.L"m(D) regime for $F_r = 0.95$",
+        xlabel = CMK.L"$D$ (mm)",
+        ylabel = CMK.L"$m$ (kg)",
+        xscale = CMK.log10,
+        yscale = CMK.log10,
+        yminorticksvisible = true,
+        yminorticks = CMK.IntervalsBetween(3),
+        xminorticksvisible = true,
+        xminorticks = CMK.IntervalsBetween(5),
+        xticks = [0.01, 0.1, 1, 10],
+        aspect = 1.67,
+        limits = ((0.02, 10.0), nothing),
+    )
+
+    sol_2 = P3.thresholds(p3, 200.0, 0.95)
+    sol_4 = P3.thresholds(p3, 400.0, 0.95)
+    sol_8 = P3.thresholds(p3, 800.0, 0.95)
+
+    #! format: off
+    fig1_b200 = CMK.lines!(ax1_b, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw)
+    fig1_b400 = CMK.lines!(ax1_b, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw)
+    fig1_b800 = CMK.lines!(ax1_b, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw)
+
+    d_thb    = CMK.vlines!(ax1_b, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
+    d_cr_200 = CMK.vlines!(ax1_b, sol_2[1] * 1e3,           linestyle = :dot,  color = cl[1], linewidth = lw)
+    d_cr_400 = CMK.vlines!(ax1_b, sol_4[1] * 1e3,           linestyle = :dot,  color = cl[2], linewidth = lw)
+    d_cr_800 = CMK.vlines!(ax1_b, sol_8[1] * 1e3,           linestyle = :dot,  color = cl[3], linewidth = lw)
+    d_gr_200 = CMK.vlines!(ax1_b, sol_2[2] * 1e3,           linestyle = :dash, color = cl[1], linewidth = lw)
+    d_gr_400 = CMK.vlines!(ax1_b, sol_4[2] * 1e3,           linestyle = :dash, color = cl[2], linewidth = lw)
+    d_gr_800 = CMK.vlines!(ax1_b, sol_8[2] * 1e3,           linestyle = :dash, color = cl[3], linewidth = lw)
+
+    leg1_b     = CMK.Legend(fig1_b[11:12, 4], [fig1_b200, fig1_b400, fig1_b800], [CMK.L"$\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false)
+    leg1_b_dth = CMK.Legend(fig1_b[11:12, 5], [d_thb],                           [CMK.L"$D_{th}$"], framevisible = false)
+    leg1_b_dcr = CMK.Legend(fig1_b[11:12, 8], [d_cr_200, d_cr_400, d_cr_800],    [CMK.L"$D_{cr}$ for $\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$D_{cr}$ for $\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$D_{cr}$ for $\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false)
+    leg1_b_dgr = CMK.Legend(fig1_b[11:12, 6], [d_gr_200, d_gr_400, d_gr_800],    [CMK.L"$D_{gr}$ for $\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$D_{gr}$ for $\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$D_{gr}$ for $\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false)
+
+    #! format: on
+    CMK.save("MorrisonandMilbrandtFig1b.svg", fig1_b)
 end
 
 #! format: off
@@ -99,9 +206,9 @@ 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.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)
+    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)
     # 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)
@@ -120,9 +227,9 @@ 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.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)
+    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)
     # 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)