P3 Adding Liquid Fraction (v2)

docs/bibliography.bib

@@ -120,6 +120,16 @@ @article{China2017
   year = {2017}
 }
 
+@article{Choletteetal2019,
+  author = {Mélissa Cholette et al},
+  title = {Parameterization of the Bulk Liquid Fraction on Mixed-Phase Particles in the Predicted Particle Properties (P3) Scheme: Description and Idealized Simulations},
+  journal = {Journal of the Atmospheric Sciences},
+  year = {2019},
+  volume = {76},
+  doi = {10.1175/JAS-D-18-0278.1},
+  pages= {561--582}
+}
+
 @article{Desai2019,
   title = {Aerosol-Mediated Glaciation of Mixed-Phase Clouds: Steady-State Laboratory Measurements},
   author = {Desai, Neel and Chandrakar, KK and Kinney, G and Cantrell, W and Shaw, RA},

docs/src/plots/Choletteetal2019_fig1.jl

@@ -0,0 +1,130 @@
+import ClimaParams
+import CloudMicrophysics as CM
+import CloudMicrophysics.P3Scheme as P3
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.Common as CO
+import CairoMakie as Plt
+
+const PSP3 = CMP.ParametersP3
+
+FT = Float64
+
+p3 = CMP.ParametersP3(FT)
+
+# Testing terminal velocity with liquid fraction
+
+function get_values(
+    p3::PSP3,
+    Chen2022::CMP.Chen2022VelType,
+    F_liq::FT,
+    F_r::FT,
+    ρ_a::FT,
+    th,
+    res::Int,
+) where {FT <: Real}
+
+    # particle dimension from 0 to 1 cm
+    D_ps = range(FT(0), stop = FT(0.01), length = res)
+
+    # ρ_r = 900
+    ρ_r = FT(900)
+
+    V = zeros(res)
+
+    for i in 1:res
+        V[i] = P3.mixed_phase_particle_velocity(
+            p3,
+            D_ps[i],
+            Chen2022,
+            ρ_a,
+            F_r,
+            F_liq,
+            th,
+        )
+    end
+    return (D_ps, V)
+end
+
+function fig1()
+    Chen2022 = CMP.Chen2022VelType(FT)
+    ρ_a = FT(1.2)
+    ρ_r = FT(900)
+
+    F_liqs = (FT(0), FT(0.33), FT(0.67), FT(1))
+    F_rs = (FT(0), FT(1 - eps(FT)))
+
+    labels = ["F_liq = 0", "F_liq = 0.33", "F_liq = 0.67", "F_liq = 1"]
+    colors = [:black, :blue, :red, :green]
+    res = 50
+
+    fig = Plt.Figure(size = (1200, 400))
+
+    #F_r = 0
+    ax1 = Plt.Axis(
+        fig[1:7, 1:9],
+        title = "Fᵣ = 0",
+        xlabel = "Dₚ (m)",
+        ylabel = "V (m s⁻¹)",
+    )
+
+    for i in 1:4
+        D_ps, V_0 = get_values(
+            p3,
+            Chen2022,
+            F_liqs[i],
+            F_rs[1],
+            ρ_a,
+            P3.thresholds(p3, ρ_r, F_rs[1]),
+            res,
+        )
+        Plt.lines!(ax1, D_ps, V_0, color = colors[i], label = labels[i])
+    end
+
+    # F_r = 1
+    ax2 = Plt.Axis(
+        fig[1:7, 10:18],
+        title = "Fᵣ = 1",
+        xlabel = "Dₚ (m)",
+        ylabel = "V (m s⁻¹)",
+    )
+    for i in 1:4
+        D_ps, V_1 = get_values(
+            p3,
+            Chen2022,
+            F_liqs[i],
+            F_rs[2],
+            ρ_a,
+            P3.thresholds(p3, ρ_r, F_rs[2]),
+            res,
+        )
+        Plt.lines!(ax2, D_ps, V_1, color = colors[i], label = labels[i])
+    end
+
+    Plt.Legend(
+        fig[4, 19:22],
+        [
+            Plt.LineElement(color = :black),
+            Plt.LineElement(color = :blue),
+            Plt.LineElement(color = :red),
+            Plt.LineElement(color = :green),
+        ],
+        labels,
+        framevisible = false,
+    )
+
+    Plt.linkaxes!(ax1, ax2)
+    ax1.xticks = (
+        [0, 0.002, 0.004, 0.006, 0.008, 0.01],
+        string.([0, 0.002, 0.004, 0.006, 0.008, 0.01]),
+    )
+    ax2.xticks = (
+        [0, 0.002, 0.004, 0.006, 0.008, 0.01],
+        string.([0, 0.002, 0.004, 0.006, 0.008, 0.01]),
+    )
+
+
+    Plt.resize_to_layout!(fig)
+    Plt.save("Choletteetal2019_fig1.svg", fig)
+end
+
+fig1()

docs/src/plots/P3DvsFliq.jl

@@ -0,0 +1,127 @@
+import ClimaParams
+import CloudMicrophysics as CM
+import CloudMicrophysics.P3Scheme as P3
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.Common as CO
+import CairoMakie as Plt
+
+const PSP3 = CMP.ParametersP3
+
+FT = Float64
+
+p3 = CMP.ParametersP3(FT)
+
+function get_values(
+    p3::PSP3,
+    q::FT,
+    N::FT,
+    F_r::FT,
+    res::Int,
+) where {FT <: Real}
+
+    # particle dimension from 0 to 0.5 cm
+    F_liqs = range(FT(0), stop = FT(1), length = res)
+
+    # ρ_r = 900
+    ρ_r = FT(900)
+
+    D = zeros(res)
+
+    for i in 1:res
+        D[i] = P3.D_m(p3, q, N, ρ_r, F_r, F_liqs[i])
+    end
+    return (F_liqs, D)
+end
+
+function fig1()
+
+    # small, med, large D_m
+    qs = (FT(0.0008), FT(0.22), FT(0.7))
+    Ns = (FT(1e6), FT(1e6), FT(1e6))
+    F_rs = (FT(0), FT(0.5), FT(0.8))
+
+    # get values
+    res = 100
+
+    fig = Plt.Figure()
+
+    labels = ["Fᵣᵢₘ = 0", "Fᵣᵢₘ = 0.5", "Fᵣᵢₘ = 0.8"]
+    colors = [:lightskyblue, :deepskyblue2, :deepskyblue4]
+
+    fig = Plt.Figure(size = (1200, 400))
+
+    #small
+    ax1 = Plt.Axis(
+        fig[1:7, 1:9],
+        title = "Small particles",
+        xlabel = "F_liq",
+        ylabel = "Dₘ (m)",
+    )
+
+    for i in 1:3
+        F_liqs, D_1 = get_values(
+            p3,
+            qs[1],
+            Ns[1],
+            F_rs[i],
+            res,
+        )
+        Plt.lines!(ax1, F_liqs, D_1, color = colors[i], label = labels[i], linewidth = 2)
+    end
+
+    #med
+    ax2 = Plt.Axis(
+        fig[1:7, 10:18],
+        title = "Medium particles",
+        xlabel = "F_liq",
+        ylabel = "Dₘ (m)",
+    )
+
+    for i in 1:3
+        F_liqs, D_2 = get_values(
+            p3,
+            qs[2],
+            Ns[2],
+            F_rs[i],
+            res,
+        )
+        Plt.lines!(ax2, F_liqs, D_2, color = colors[i], label = labels[i], linewidth = 2)
+    end
+
+    #large
+    ax3 = Plt.Axis(
+        fig[1:7, 19:27],
+        title = "Large particles",
+        xlabel = "F_liq",
+        ylabel = "Dₘ (m)",
+    )
+
+    for i in 1:3
+        F_liqs, D_3 = get_values(
+            p3,
+            qs[3],
+            Ns[3],
+            F_rs[i],
+            res,
+        )
+        Plt.lines!(ax3, F_liqs, D_3, color = colors[i], label = labels[i], linewidth = 2)
+    end
+
+    Plt.Legend(
+        fig[4, 28:31],
+        [
+            Plt.LineElement(color = :lightskyblue),
+            Plt.LineElement(color = :deepskyblue2),
+            Plt.LineElement(color = :deepskyblue4),
+        ],
+        labels,
+        framevisible = false,
+    )
+
+    # Plt.linkaxes!(ax1, ax2, ax3)
+    Plt.resize_to_layout!(fig)
+    Plt.display(fig)
+    # Plt.save("P3DvsFliq.svg", fig)
+end
+
+fig1()

docs/src/plots/P3LambdaErrorPlots.jl

@@ -10,6 +10,7 @@ FT = Float64
 
 const PSP3 = CMP.ParametersP3
 p3 = CMP.ParametersP3(FT)
+F_liq = FT(0)
 
 function λ_diff(F_r::FT, ρ_r::FT, N::FT, λ_ex::FT, p3::PSP3) where {FT}
 
@@ -20,9 +21,10 @@ function λ_diff(F_r::FT, ρ_r::FT, N::FT, λ_ex::FT, p3::PSP3) where {FT}
     # Convert λ to ensure it remains positive
     x = log(λ_ex)
     # Compute mass density based on input shape parameters
-    q_calc = N * P3.q_over_N_gamma(p3, F_r, x, μ, th)
+    q_calc = N * P3.q_over_N_gamma(p3, F_liq, F_r, x, μ, th)
 
-    (λ_calculated,) = P3.distribution_parameter_solver(p3, q_calc, N, ρ_r, F_r)
+    (λ_calculated,) =
+        P3.distribution_parameter_solver(p3, q_calc, N, ρ_r, F_liq, F_r)
     return abs(λ_ex - λ_calculated)
 end
 
@@ -149,12 +151,12 @@ function μ_approximation_effects(F_r::FT, ρ_r::FT) where {FT}
     λ_solved = [FT(0) for λ in λs]
 
     for i in 1:numpts
-        q = P3.q_over_N_gamma(p3, F_r, log(λs[i]), μs[i], th)
+        q = P3.q_over_N_gamma(p3, F_liq, F_r, log(λs[i]), μs[i], th)
         qs[i] = q
         N = FT(1e6)
-        (L, N) = P3.distribution_parameter_solver(p3, q * N, N, ρ_r, F_r)
+        (L, N) = P3.distribution_parameter_solver(p3, q * N, N, ρ_r, F_liq, F_r)
         λ_solved[i] = L
-        μs_approx[i] = P3.DSD_μ_approx(p3, N * q, N, ρ_r, F_r)
+        μs_approx[i] = P3.DSD_μ_approx(p3, N * q, N, ρ_r, F_liq, F_r)
     end
 
     # Plot

docs/src/plots/P3SchemePlots.jl

@@ -8,6 +8,7 @@ FT = Float64
 
 const PSP3 = CMP.ParametersP3
 p3 = CMP.ParametersP3(FT)
+F_liq = FT(0)
 
 
 
@@ -51,13 +52,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.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.p3_mass(p3, D, 0.0, F_liq        ) 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, F_liq, 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, F_liq, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
     # a(D)
-    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)
+    fig2_0 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.0, F_liq        ) 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, F_liq, 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, F_liq, 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)
@@ -72,13 +73,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.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.p3_mass(p3, D, 0.95, F_liq, 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, F_liq, 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, F_liq, sol_8) for D in D_range], color = cl[3], linewidth = lw)
     # a(D)
-    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)
+    fig3_200 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, F_liq, 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, F_liq, 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, F_liq, 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)

docs/src/plots/P3ShapeSolverPlots.jl

@@ -8,6 +8,7 @@ FT = Float64
 
 const PSP3 = CMP.ParametersP3
 p3 = CMP.ParametersP3(FT)
+F_liq = FT(0)
 
 function guess_value(λ::FT, p1::FT, p2::FT, q1::FT, q2::FT)
     return q1 * (λ / p1)^((log(q1) - log(q2)) / (log(p1) - log(p2)))
@@ -40,16 +41,16 @@ function lambda_guess_plot()
             λs = [10^logλ for logλ in logλs]
             th = P3.thresholds(p3, ρ_r, F_r)
             qs = [
-                P3.q_over_N_gamma(p3, F_r, log(λ), P3.DSD_μ(p3, λ), th) for
-                λ in λs
+                P3.q_over_N_gamma(p3, F_liq, F_r, log(λ), P3.DSD_μ(p3, λ), th) for λ in λs
             ]
             guesses = [FT(0) for λ in λs]
 
             for i in 1:length(λs)
                 (min,) = P3.get_bounds(
                     N,
                     qs[i] * N,
-                    P3.DSD_μ_approx(p3, qs[i] * N, N, ρ_r, F_r),
+                    P3.DSD_μ_approx(p3, qs[i] * N, N, ρ_r, F_r, F_liq),
+                    F_liq,
                     F_r,
                     p3,
                     th,