splitting into 3 plots

parcel/P3_ice_nuc.jl

@@ -0,0 +1,114 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+
+# definition of the ODE problem for parcel model
+include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
+FT = Float32
+# get free parameters
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+aps = CMP.AirProperties(FT)
+wps = CMP.WaterProperties(FT)
+ip = CMP.IceNucleationParameters(FT)
+
+# Constants
+ρₗ = wps.ρw
+R_v = TD.Parameters.R_v(tps)
+R_d = TD.Parameters.R_d(tps)
+
+# Initial conditions
+Nₐ = FT(2000)
+Nₗ = FT(2000)
+Nᵢ = FT(0)
+rₗ = FT(1.25e-6)
+p₀ = FT(20000)
+T₀_dep = FT(230)
+T₀_het = FT(240)
+T₀_hom = FT(233.17)
+qᵥ = FT(3.3e-4)
+qₗ = FT(Nₗ * 4 / 3 * π * rₗ^3 * ρₗ / 1.2)
+qᵢ = FT(0)
+x_sulph = FT(0)
+
+# Moisture dependent initial conditions
+q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q)
+ρₐ = TD.air_density(tps, ts)
+Rₐ = TD.gas_constant_air(tps, q)
+eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
+e = qᵥ * p₀ * R_v / Rₐ
+
+Sₗ = FT(e / eₛ)
+IC_dep = [Sₗ, p₀, T₀_dep, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC_het = [Sₗ, p₀, T₀_het, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC_hom = [Sₗ, p₀, T₀_hom, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+
+ξ(T) =
+    TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
+    TD.saturation_vapor_pressure(tps, T, TD.Ice())
+S_i(T, S_liq) = ξ(T) * S_liq
+
+# Simulation parameters passed into ODE solver
+r_nuc = FT(1.25e-6)                     # assumed size of nucleated particles
+w = FT(3.5 * 1e-2)                      # updraft speed
+α_m = FT(0.5)                           # accomodation coefficient
+const_dt = FT(0.1)                      # model timestep
+t_max = FT(80)
+aerosol = []
+ice_nucleation_modes_list = [["P3_Deposition"], ["P3_het"], ["P3_hom"]]
+growth_modes = ["Deposition"]
+droplet_size_distribution_list = [["Monodisperse"]]
+
+# Plotting
+fig = MK.Figure(resolution = (800, 600))
+ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation [-]")
+ax2 = MK.Axis(fig[1, 2], ylabel = "ICNC [cm^-3]")
+ax3 = MK.Axis(fig[2, 1], ylabel = "Ice Saturation [-]")
+ax4 = MK.Axis(fig[2, 2], ylabel = "ICNC [cm^-3]")
+ax5 = MK.Axis(fig[3, 1], ylabel = "Ice Saturation [-]")
+ax6 = MK.Axis(fig[3, 2], ylabel = "ICNC [cm^-3]")
+
+for ice_nucleation_modes in ice_nucleation_modes_list
+    nuc_mode = ice_nucleation_modes[1]
+    droplet_size_distribution = droplet_size_distribution_list[1]
+    p = (;
+        wps,
+        aps,
+        tps,
+        ip,
+        const_dt,
+        r_nuc,
+        w,
+        α_m,
+        aerosol,
+        ice_nucleation_modes,
+        growth_modes,
+        droplet_size_distribution,
+    )
+    # solve ODE
+    #! formatt: off
+    if "P3_Deposition" in ice_nucleation_modes
+        sol = run_parcel(IC_dep, FT(0), t_max, p)
+        MK.lines!(ax1, sol.t, S_i.(sol[3, :], (sol[1, :])), label = nuc_mode) # saturation
+        MK.lines!(ax2, sol.t, sol[9, :] * 1e-6, label = nuc_mode)             # ICNC
+        MK.axislegend(ax1, framevisible = false, labelsize = 12, orientation = :horizontal, position = :rt)
+    elseif "P3_het" in ice_nucleation_modes
+        sol = run_parcel(IC_het, FT(0), t_max, p)
+        MK.lines!(ax3, sol.t, S_i.(sol[3, :], (sol[1, :])), label = nuc_mode, color = :red) # saturation
+        MK.lines!(ax4, sol.t, sol[9, :] * 1e-6, label = nuc_mode, color = :red)             # ICNC
+        MK.axislegend(ax3, framevisible = false, labelsize = 12, orientation = :horizontal, position = :lt)
+    elseif "P3_hom" in ice_nucleation_modes
+        sol = run_parcel(IC_hom, FT(0), t_max, p)
+        MK.lines!(ax5, sol.t, S_i.(sol[3, :], (sol[1, :])), label = nuc_mode, color = :orange) # saturation
+        MK.lines!(ax6, sol.t, sol[9, :] * 1e-6, label = nuc_mode, color = :orange)             # ICNC
+        MK.axislegend(ax5, framevisible = false, labelsize = 12, orientation = :horizontal, position = :lt)
+    end
+    #! formatt: on
+end
+
+
+
+#MK.save("P3_ice_nuc.svg", fig)
+fig

parcel/parcel.jl

@@ -107,6 +107,11 @@ function parcel_model(dY, Y, p, t)
             dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air
         end
     end
+    if "P3_Deposition" in ice_nucleation_modes
+        Nᵢ_depo = CMI_het.P3_deposition_N_i(T)
+        dN_act_dt_depo = max(FT(0), (Nᵢ_depo - N_ice) / const_dt)
+        dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air
+    end
 
     dN_act_dt_immersion = FT(0)
     dqi_dt_new_immers = FT(0)
@@ -126,6 +131,19 @@ function parcel_model(dY, Y, p, t)
         dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_l)
         dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air
     end
+    if "P3_het" in ice_nucleation_modes
+        if "Monodisperse" in droplet_size_distribution
+            r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air)
+        elseif "Gamma" in droplet_size_distribution
+            λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air)
+            #A = N_liq* λ^2
+            r_l = 2 / λ
+        end
+        V_l = FT(4 / 3 * π * r_l^3)
+        Nᵢ_het = CMI_het.P3_het_N_i(T, N_liq, V_l, t)
+        dN_act_dt_immersion = max(FT(0), (Nᵢ_het - N_ice) / const_dt)
+        dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air
+    end
 
     dN_act_dt_hom = FT(0)
     dqi_dt_new_hom = FT(0)
@@ -159,6 +177,23 @@ function parcel_model(dY, Y, p, t)
             dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_aero^3 * ρ_ice / ρ_air
         end
     end
+    if "P3_hom" in ice_nucleation_modes
+        if "Monodisperse" in droplet_size_distribution
+            r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air)
+        elseif "Gamma" in droplet_size_distribution
+            λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air)
+            #A = N_liq* λ^2
+            r_l = 2 / λ
+        elseif "Lognormal" in droplet_size_distribution
+            r_l = R_mode_liq * exp(σ)
+        end
+        if T < FT(233.15) && N_liq > FT(0)
+            dN_act_dt_hom = N_liq / const_dt
+        else
+            dN_act_dt_hom = FT(0)
+        end
+        dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air
+    end
 
     dN_ice_dt = dN_act_dt_depo + dN_act_dt_immersion + dN_act_dt_hom
     if "Jensen_Example" in droplet_size_distribution
@@ -306,12 +341,25 @@ function run_parcel(IC, t_0, t_end, p)
         print("\n    Water activity based deposition nucleation ")
         print("(Note that this mode only runs for monodisperse size distributions.) ")
     end
+    if "P3_Deposition" in ice_nucleation_modes
+        print("\n    P3 deposition nucleation ")
+        timestepper = ODE.Euler()
+    end
     if "ImmersionFreezing" in ice_nucleation_modes
-        print("\n    Immersion freezing ")
+        print("\n    Immersion freezing (ABIFM) ")
+        print("(Note that this mode only runs for monodisperse and gamma size distributions.) ")
+    end
+    if "P3_het" in ice_nucleation_modes
+        print("\n    P3 heterogeneous freezing ")
         print("(Note that this mode only runs for monodisperse and gamma size distributions.) ")
+        timestepper = ODE.Euler()
     end
     if "HomogeneousFreezing" in ice_nucleation_modes
-        print("\n    Homogeneous freezing ")
+        print("\n    Water activity based homogeneous freezing ")
+    end
+    if "P3_hom" in ice_nucleation_modes
+        print("\n    P3 homogeneous freezing ")
+        timestepper = ODE.Euler()
     end
     print("\n")
 
@@ -338,8 +386,14 @@ function run_parcel(IC, t_0, t_end, p)
     if "Lognormal" in droplet_size_distribution
         print("\n    Lognormal size distribution")
     end
-    if "MohlerAF_Deposition" in ice_nucleation_modes || "MohlerRate_Deposition" in ice_nucleation_modes || "ActivityBasedDeposition" in ice_nucleation_modes
-        print("\nWARNING: Multiple deposition nucleation parameterizations chosen to run in one simulation. Parcel will only run MohlerAF_Deposition for now.")
+    if "MohlerAF_Deposition" in ice_nucleation_modes ||
+        "MohlerRate_Deposition" in ice_nucleation_modes ||
+        "ActivityBasedDeposition" in ice_nucleation_modes ||
+        "P3_Deposition" in ice_nucleation_modes
+        print("\nWARNING: Multiple deposition nucleation parameterizations chosen to run in one simulation. Parcel can only properly handle one for now.")
+    end
+    if "P3_het" in ice_nucleation_modes
+        print("\nWARNING: Assuming P3_het parameterization is a parameterization for immeresion freezing.")
     end
     println(" ")