diff --git a/parcel/Jensen_et_al_2022.jl b/parcel/Jensen_et_al_2022.jl index 99f5097baa..75d15e169c 100644 --- a/parcel/Jensen_et_al_2022.jl +++ b/parcel/Jensen_et_al_2022.jl @@ -31,59 +31,27 @@ 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 @@ -91,7 +59,6 @@ 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) diff --git a/parcel/parcel.jl b/parcel/parcel.jl index 5e8b119174..b8d91d0039 100644 --- a/parcel/parcel.jl +++ b/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) diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl index 3e50719bee..249dd419e5 100644 --- a/src/IceNucleation.jl +++ b/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