works! enough at least

parcel/Jensen_et_al_2022.jl

@@ -31,67 +31,34 @@ s_i(T, s_liq) = @. s_liq * ξ(T) # saturation ratio
 
 # Initial conditions
 Nₐ = FT(300 * 1e6)
-Nₗ = FT(0)       # Jensen2022 starts with Nₐ = 300 * 1e6 and activates them within parcel
+Nₗ = FT(0)
 Nᵢ = FT(0)
-r₀ = FT(25 * 1e-9)         # Value of dry aerosol r from Jensen2022
-#r₀ = FT(7 * 1e-8)          # starting when Jensen freezes
+r₀ = FT(25 * 1e-9)
 p₀ = FT(121 * 1e2)
 T₀ = FT(190)
-#T₀ = FT(189.7)             # starting when Jensen freezes
 x_sulph = FT(0)
 
 # saturation and partial pressure
 eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
-# Sᵢ = FT(1.55)
-# Sₗ = S_l(T₀, Sᵢ)
-# e = eₛ * Sₗ
-# ϵ = R_d / R_v
-# starting when Jensen freezes
-#q_vap = FT(3.76e-7)
-q_vap = FT(5e-6) * FT(0.6) / FT(1.2)
-#q_vap = FT(4.1e-7) # this gives more valid ICNC and ice saturation curves
-q_liq = Nₗ * 4 / 3 * π * r₀^3 * ρₗ / 1.2
-q_ice = FT(0)
-#q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+qᵥ = FT(5e-6) * FT(0.6) / FT(1.2)
+qₗ = Nₗ * 4 / 3 * π * r₀^3 * ρₗ / 1.2
+qᵢ = FT(0)
 q = TD.PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice)
 R_a = TD.gas_constant_air(tps, q)
 e = q_vap * p₀ * R_v / R_a 
-sₗ = e / eₛ # saturation ratio
-println("sᵢ is ", s_i(T₀, sₗ))  # saturation ratio
-
-# 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)
+sₗ = e / eₛ 
 
 ## 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)
-
-#IC = [sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
-IC = [sₗ, p₀, T₀, q_vap, q_liq, q_ice, Nₐ, Nₗ, Nᵢ, x_sulph] # starting when Jensen freezes
 
-println("initial conditions:")
-println("   liquid saturation: ", sₗ)
-# println("   vapor mixing ratio: ", qᵥ)
-# println("   liquid mixing ratio: ", qₗ)
-# println("   ice mixing ratio: ", qᵢ)
-println("   vapor mixing ratio: ", q_vap)   # starting when Jensen freezes
-println("   liquid mixing ratio: ", q_liq)  # starting when Jensen freezes
+IC = [sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
 
 # Simulation parameters passed into ODE solver
 r_nuc = r₀                                 # assumed size of nucleated particles
 w = FT(1)                                  # updraft speed
 α_m = FT(1)                                # accomodation coefficient
 const_dt = FT(0.01)                        # model timestep
 t_max = FT(140)
-#t_max = FT(84)  # start when Jensen freezes
 aerosol = []
 ice_nucleation_modes = ["HomogeneousFreezing"]     # homogeneous freezing only
 growth_modes = ["Deposition"]
@@ -114,7 +81,7 @@ 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]")
-#ax6 = MK.Axis(fig[3, 2], ylabel = "r_liq [m]")
+ax6 = MK.Axis(fig[3, 2], ylabel = "q_ice [g/kg]")
 
 MK.ylims!(ax2, 188.5, 190)
 MK.xlims!(ax2, -5, 150)
@@ -152,6 +119,7 @@ for droplet_size_distribution in droplet_size_distribution_list
     MK.lines!(ax3, sol.t, sol[4, :] * 1e3) # q_vap
     MK.lines!(ax4, sol.t, sol[5, :] * 1e3) # q_liq
     MK.lines!(ax5, sol.t, sol[9, :] * 1e-6)# ICNC
+    MK.lines!(ax6, sol.t, sol[6, :] * 1e3) # q_ice
 end
 
 #MK.save("Jensen_et_al_2022.svg", fig)

parcel/parcel.jl

@@ -128,7 +128,7 @@ function parcel_model(dY, Y, p, t)
         elseif "Lognormal" in droplet_size_distribution
             #dN_act_dt_hom = max(FT(0), J_hom * N_aer * 4 / 3 * π * exp((6*log(R_mode_aero) + 9 * σ^2)/(2)))
             dN_act_dt_hom = max(FT(0), J_hom * N_aer * 4 / 3 * π * exp((6*log(r_nuc) + 36)/(2)))
-            r_aero = r_nuc * exp(2) # mean radius
+            r_aero = r_nuc * exp(2)  # mean radius, m. σ = 2
             dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * π * r_aero^3 * ρ_ice / ρ_air
         end
 
@@ -158,8 +158,8 @@ function parcel_model(dY, Y, p, t)
             r_i = 2 / λ
             C_i = r_i
             dqi_dt_depo = 4 * π / ρ_air * (s_i - 1) * G_i * r_i * N_ice
-        elseif "Lognormal" in droplet_size_distribution # TODO - monodisperse for now. need to figure out how to do this lognormally
-            r_i = cbrt(q_ice / N_ice / (4 / 3 * π) / ρ_ice * ρ_air)
+        elseif "Lognormal" in droplet_size_distribution 
+            r_i = cbrt(q_ice / N_ice / (4 / 3 * π) / ρ_ice * ρ_air) * exp(2)
             dqi_dt_depo = 4 * π / ρ_air * (s_i - 1) * G_i * r_i * N_ice
         end
         # TODO - check if r_i is non-zero if initial conditions say q_ice = 0
@@ -195,8 +195,7 @@ function parcel_model(dY, Y, p, t)
     dq_liq_dt = dq_liq_dt_lv + dq_liq_dt_li
     ds_l_dt = a1 * w * s_liq - (a2 + a3 * s_liq) * dq_liq_dt_lv - (a2 + a4 * s_liq) * dq_ice_dt_iv - a5 * s_liq * dq_ice_dt_il
     dp_a_dt = -p_a * grav / R_a / T * w
-    #dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt_lv + L_fus / cp_a * dq_ice_dt_il + L_subl / cp_a * dq_ice_dt_iv
-    dT_dt = FT(- 1 / 120)
+    dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt_lv + L_fus / cp_a * dq_ice_dt_il + L_subl / cp_a * dq_ice_dt_iv
     println("------------ dT/dt = ", dT_dt)
     println("------------ T = ", T)
     println("------------ q_vap = ", q_vap)

src/IceNucleation.jl

@@ -92,7 +92,7 @@ function homogeneous_J(ip::CMP.Koop2000, Δa_w::FT) where {FT}
     elseif Δa_w > ip.Δa_w_max
         Δa_w = ip.Δa_w_max
     end
-    
+
     logJ::FT = ip.c₁ + ip.c₂ * Δa_w - ip.c₃ * Δa_w^2 + ip.c₄ * Δa_w^3
 
     return FT(10)^(logJ) * 1e6