testing stuff out

parcel/Jensen_et_al_2022.jl

@@ -0,0 +1,135 @@
+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 = Float64
+
+# Get free parameters
+tps = CMP.ThermodynamicsParameters(FT)
+wps = CMP.WaterProperties(FT)
+aps = CMP.AirProperties(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(0)
+Nₗ = FT(300 * 1e6)       # Jensen2022 starts with Nₐ = 300 * 1e6 and activates them within parcel
+Nᵢ = FT(0)
+r₀ = FT(25 * 1e-9)      # Value of dry aerosol r from Jensen2022
+p₀ = FT(800 * 1e2)
+T₀ = FT(190)
+x_sulph = FT(0)
+
+
+eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
+ξ(T) =
+TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
+TD.saturation_vapor_pressure(tps, T, TD.Ice())
+S_l(T, S_i) = @. (S_i + 1) / ξ(T)
+Sᵢ = FT(1.55)
+Sₗ = S_l(T₀, Sᵢ)
+e = eₛ * Sₗ
+ϵ = R_d / R_v
+
+# mass per volume for dry air, vapor and liquid
+md_v = (p₀ - e) / R_d / T₀
+mv_v = e / R_v / T₀
+ml_v = Nₗ * 4 / 3 * π * ρₗ * r₀^3
+
+qᵥ = mv_v / (md_v + mv_v + ml_v)
+qₗ = ml_v / (md_v + mv_v + ml_v)
+qᵢ = 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 = qᵥ * p₀ * R_v / Rₐ
+
+println("S_l initial = ", Sₗ )
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+
+# Simulation parameters passed into ODE solver
+r_nuc = FT(25 * 1e-9)                      # assumed size of nucleated particles
+w = FT(1)                                  # updraft speed
+α_m = FT(0.5)                              # accomodation coefficient
+const_dt = FT(0.5)                         # model timestep
+#t_max = FT(120)
+t_max = FT(30)
+aerosol = []
+ice_nucleation_modes = ["HomogeneousFreezing"]     # homogeneous freezing only
+growth_modes = []                                  # no growth
+droplet_size_distribution_list = [["Monodisperse"]]
+
+# Data from Jensen(2022) Figure 1
+# https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022JD036535
+#! format: off
+Jensen_t_sat = [0, 62.71, 70.52, 76.87, 82.4, 84.84, 88.1, 92, 96.07, 100.63, 105.35, 112.51, 119.83]
+Jensen_sat = [1.55, 1.694, 1.7107, 1.7208, 1.725, 1.726, 1.7259, 1.722, 1.715, 1.702, 1.686, 1.653, 1.6126]
+Jensen_t_T = [0, 120]
+Jensen_T = [190, 189]
+Jensen_t_ICNC = [0.217, 42.69, 50.02, 54.41, 58.97, 65.316, 72.477, 82.08, 92.658, 94.123, 95.5877, 119.84]
+Jensen_ICNC = [0, 0, 0.282, 0.789, 1.804, 4.1165, 7.218, 12.12, 16.35, 16.8, 16.97, 17.086]
+#! format: on
+
+fig = MK.Figure(resolution = (800, 600))
+ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation")
+ax2 = MK.Axis(fig[3, 1], xlabel = "Time [s]", ylabel = "Temperature [K]")
+ax3 = MK.Axis(fig[2, 1], ylabel = "q_vap [g/kg]")
+ax4 = MK.Axis(fig[2, 2], xlabel = "Time [s]", ylabel = "q_liq [g/kg]")
+ax5 = MK.Axis(fig[1, 2], ylabel = "ICNC [cm^-3]")
+
+MK.lines!(ax1, Jensen_t_sat, Jensen_sat, label = "Jensen 2022", color = :green)
+MK.lines!(ax2, Jensen_t_T, Jensen_T, color = :green)
+MK.lines!(ax5, Jensen_t_ICNC, Jensen_ICNC, color = :green)
+
+for droplet_size_distribution in droplet_size_distribution_list
+    p = (;
+        wps,
+        aps,
+        tps,
+        ip,
+        const_dt,
+        r_nuc,
+        w,
+        α_m,
+        aerosol,
+        ice_nucleation_modes,
+        growth_modes,
+        droplet_size_distribution,
+    )
+    # solve ODE
+    sol = run_parcel(IC, FT(0), t_max, p)
+
+    DSD = droplet_size_distribution[1]
+
+    # Plot results
+    MK.lines!(ax1, sol.t, S_i(sol[3, :], (sol[1, :])), label = DSD)
+    MK.lines!(ax2, sol.t, sol[3, :])
+    MK.lines!(ax3, sol.t, sol[4, :] * 1e3)
+    MK.lines!(ax4, sol.t, sol[5, :] * 1e3)
+    MK.lines!(ax5, sol.t, sol[9, :] * 1e6)
+end
+
+#fig[3, 2] = MK.Legend(fig, ax1, framevisible = false)
+MK.axislegend(
+    ax1,
+    framevisible = false,
+    labelsize = 12,
+    orientation = :horizontal,
+    nbanks = 2,
+    position = :rb,
+)
+
+MK.save("Jensen_et_al_2022.svg", fig)

parcel/parcel.jl

@@ -109,12 +109,24 @@ function parcel_model(dY, Y, p, t)
         Δa_w = T > FT(185) && T < FT(235) ?
             CMO.a_w_xT(H2SO4_prs, tps, x_sulph, T) - CMO.a_w_ice(tps, T) :
             CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T)
-        J_hom = CMI_hom.homogeneous_J(ip, Δa_w)
+        println("a_w_ice = ", CMO.a_w_ice(tps, T))
+        println("p_ice = ", TD.saturation_vapor_pressure(tps, T, TD.Ice()))
+        println("Δa_w = ", Δa_w)
+        # Δa_w = T > FT(185) && T < FT(235) ?
+        # CMO.a_w_xT(H2SO4_prs, tps, x_sulph, T) - (0.047/TD.saturation_vapor_pressure(tps, T, TD.Liquid())) :
+        # CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T)
+        # println("a_w_ice = ", (0.047/TD.saturation_vapor_pressure(tps, T, TD.Liquid())))
+        # println("Δa_w = ", Δa_w)
+        # println("T = ", T)
+        J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w)
         if "Monodisperse" in droplet_size_distribution
-            r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air)
+            r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air) / 10
+            println("radius of drople = ", r_l)
             V_l = 4 / 3 * π * r_l^3
             dN_act_dt_hom = max(FT(0), J_hom * N_liq * V_l)
             dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air
+            println("J_hom = ", J_hom)
+            println("dqi_dt_new_hom = ", dqi_dt_new_hom)
         end
         if "Gamma" in droplet_size_distribution
             λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air)
@@ -179,7 +191,7 @@ function parcel_model(dY, Y, p, t)
     dY[10] = FT(0)         # sulphuric acid concentration
 
     # TODO - add diagnostics output (radius, S, etc)
-
+    println("\n")
 end
 #! format: on