From 5829cca34a27299454fc8de963a03f762b647e95 Mon Sep 17 00:00:00 2001 From: amylu00 Date: Tue, 31 Oct 2023 08:39:21 -0700 Subject: [PATCH] testing stuff out --- parcel/Jensen_et_al_2022.jl | 135 ++++++++++++++++++++++++++++++++++++ parcel/parcel.jl | 18 ++++- 2 files changed, 150 insertions(+), 3 deletions(-) create mode 100644 parcel/Jensen_et_al_2022.jl diff --git a/parcel/Jensen_et_al_2022.jl b/parcel/Jensen_et_al_2022.jl new file mode 100644 index 0000000000..3fd83411d6 --- /dev/null +++ b/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) diff --git a/parcel/parcel.jl b/parcel/parcel.jl index 0e24e446df..cfaeb32b83 100644 --- a/parcel/parcel.jl +++ b/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