Add terminal velocities to P3 scheme

Author: anastasia-popova

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

@@ -68,7 +68,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
 
@@ -99,9 +143,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 +164,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)

docs/src/plots/P3TerminalVelocityPlots.jl

@@ -0,0 +1,255 @@
+import CLIMAParameters
+import CloudMicrophysics as CM
+import CloudMicrophysics.P3Scheme as P3
+import CloudMicrophysics.Parameters as CMP
+import CairoMakie as Plt
+
+const PSP3 = CMP.ParametersP3
+
+FT = Float64
+
+p3 = CMP.ParametersP3(FT)
+
+function get_values(
+    p3::PSP3,
+    Chen2022::CMP.Chen2022VelTypeSnowIce,
+    q::FT,
+    N::FT,
+    ρ_a::FT,
+    x_resolution::Int,
+    y_resolution::Int,
+) where {FT}
+    F_rs = range(FT(0), stop = FT(1 - eps(FT)), length = x_resolution)
+    ρ_rs = range(FT(25), stop = FT(975), length = y_resolution)
+
+    V_m = zeros(x_resolution, y_resolution)
+    D_m = zeros(x_resolution, y_resolution)
+
+    for i in 1:x_resolution
+        for j in 1:y_resolution
+            F_r = F_rs[i]
+            ρ_r = ρ_rs[j]
+
+            V_m[i, j] =
+                P3.terminal_velocity_mass(p3, Chen2022, q, N, ρ_r, F_r, ρ_a)
+            # get D_m in mm for plots 
+            D_m[i, j] = 1e3 * P3.D_m(p3, q, N, ρ_r, F_r)
+        end
+    end
+    return (; F_rs, ρ_rs, V_m, D_m)
+end
+
+function figure_2()
+    Chen2022 = CMP.Chen2022VelType(FT)
+    # density of air in kg/m^3
+    ρ_a = FT(1.2) #FT(1.293)
+
+    f = Plt.Figure()
+    xres = 100
+    yres = 100
+    min = FT(0)
+    max = FT(10)
+
+    # small D_m
+    q_s = FT(0.0008)
+    N_s = FT(1e6)
+
+    # medium D_m 
+    q_m = FT(0.22)
+    N_m = FT(1e6)
+
+    # large D_m
+    q_l = FT(0.7)
+    N_l = FT(1e6)
+
+    # get V_m and D_m
+    (F_rs, ρ_rs, V_ms, D_ms) =
+        get_values(p3, Chen2022.snow_ice, q_s, N_s, ρ_a, xres, yres)
+    (F_rm, ρ_rm, V_mm, D_mm) =
+        get_values(p3, Chen2022.snow_ice, q_m, N_m, ρ_a, xres, yres)
+    (F_rl, ρ_rl, V_ml, D_ml) =
+        get_values(p3, Chen2022.snow_ice, q_l, N_l, ρ_a, xres, yres)
+
+    ax1 = Plt.Axis(
+        f[1, 1],
+        xlabel = "F_r",
+        ylabel = "ρ_r",
+        title = "Small Dₘ",
+        width = 350,
+        height = 350,
+        limits = (0, 1.0, 25, 975),
+        xticks = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0],
+        yticks = [200, 400, 600, 800],
+    )
+
+    Plt.contourf!(
+        F_rs,
+        ρ_rs,
+        V_ms,
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+    )
+    Plt.contour!(F_rs, ρ_rs, D_ms, color = :black, labels = true, levels = 3)
+    Plt.Colorbar(
+        f[2, 1],
+        limits = (
+            minimum(v for v in V_ms if !isnan(v)),
+            maximum(v for v in V_ms if !isnan(v)),
+        ),
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+        flipaxis = true,
+        vertical = false,
+    )
+
+    ax2 = Plt.Axis(
+        f[1, 2],
+        xlabel = "F_r",
+        ylabel = "ρ_r",
+        title = "Medium Dₘ",
+        width = 350,
+        height = 350,
+        limits = (0, 1.0, 25, 975),
+        xticks = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0],
+        yticks = [200, 400, 600, 800],
+    )
+
+    Plt.contourf!(
+        F_rm,
+        ρ_rm,
+        V_mm,
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+    )
+    Plt.contour!(F_rm, ρ_rm, D_mm, color = :black, labels = true, levels = 3)
+    Plt.Colorbar(
+        f[2, 2],
+        limits = (
+            minimum(v for v in V_mm if !isnan(v)),
+            maximum(v for v in V_mm if !isnan(v)),
+        ),
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+        flipaxis = true,
+        vertical = false,
+    )
+
+    ax3 = Plt.Axis(
+        f[1, 3],
+        xlabel = "F_r",
+        ylabel = "ρ_r",
+        title = "Large Dₘ",
+        width = 350,
+        height = 350,
+        limits = (0, 1.0, 25, 975),
+        xticks = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0],
+        yticks = [200, 400, 600, 800],
+    )
+
+    Plt.contourf!(
+        F_rl,
+        ρ_rl,
+        V_ml,
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+    )
+    Plt.contour!(F_rl, ρ_rl, D_ml, color = :black, labels = true, levels = 3)
+    Plt.Colorbar(
+        f[2, 3],
+        limits = (
+            minimum(v for v in V_ml if !isnan(v)),
+            maximum(v for v in V_ml if !isnan(v)),
+        ),
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+        flipaxis = true,
+        vertical = false,
+    )
+
+    Plt.linkxaxes!(ax1, ax2)
+    Plt.linkxaxes!(ax2, ax3)
+    Plt.linkyaxes!(ax1, ax2)
+    Plt.linkyaxes!(ax2, ax3)
+
+    # Plot D_m as second row of comparisons
+
+    ax4 = Plt.Axis(
+        f[3, 1],
+        height = 350,
+        width = 350,
+        xlabel = "F_r",
+        ylabel = "ρ_r",
+        title = "Small Dₘ vs F_r and ρ_r",
+    )
+    Plt.contourf!(
+        F_rs,
+        ρ_rs,
+        D_ms,
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+    )
+    Plt.Colorbar(
+        f[4, 1],
+        limits = (
+            minimum(d for d in D_ms if !isnan(d)),
+            maximum(d for d in D_ms if !isnan(d)),
+        ),
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+        flipaxis = true,
+        vertical = false,
+    )
+    ax5 = Plt.Axis(
+        f[3, 2],
+        height = 350,
+        width = 350,
+        xlabel = "F_r",
+        ylabel = "ρ_r",
+        title = "Medium Dₘ vs F_r and ρ_r",
+    )
+    Plt.contourf!(
+        F_rm,
+        ρ_rm,
+        D_mm,
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+    )
+    Plt.Colorbar(
+        f[4, 2],
+        limits = (
+            minimum(d for d in D_mm if !isnan(d)),
+            maximum(d for d in D_mm if !isnan(d)),
+        ),
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+        flipaxis = true,
+        vertical = false,
+    )
+    ax6 = Plt.Axis(
+        f[3, 3],
+        height = 350,
+        width = 350,
+        xlabel = "F_r",
+        ylabel = "ρ_r",
+        title = "Large Dₘ vs F_r and ρ_r",
+    )
+    Plt.contourf!(
+        F_rl,
+        ρ_rl,
+        D_ml,
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+    )
+    Plt.Colorbar(
+        f[4, 3],
+        limits = (
+            minimum(d for d in D_ml if !isnan(d)),
+            maximum(d for d in D_ml if !isnan(d)),
+        ),
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+        flipaxis = true,
+        vertical = false,
+    )
+
+    Plt.linkxaxes!(ax1, ax4)
+    Plt.linkxaxes!(ax4, ax5)
+    Plt.linkxaxes!(ax5, ax6)
+    Plt.linkyaxes!(ax1, ax4)
+    Plt.linkyaxes!(ax4, ax5)
+    Plt.linkyaxes!(ax5, ax6)
+
+    Plt.resize_to_layout!(f)
+    Plt.save("MorrisonandMilbrandtFig2.svg", f)
+end
+
+# Terminal Velocity figure
+figure_2()