From 5ea2106f88362b8cb1bac1f50c29d491a4bd4010 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Wed, 21 Feb 2024 00:12:05 -0800 Subject: [PATCH] Refactor parcel --- ...on.jl => Example_Deposition_Nucleation.jl} | 92 ++-- ...ezing.jl => Example_Immersion_Freezing.jl} | 66 +-- parcel/Example_Jensen_et_al_2022.jl | 150 +++++++ ...{Liquid_only.jl => Example_Liquid_only.jl} | 80 ++-- parcel/Example_P3_ice_nuc.jl | 106 +++++ ...al_2023.jl => Example_Tully_et_al_2023.jl} | 70 ++- parcel/Jensen_et_al_2022.jl | 159 ------- parcel/P3_ice_nuc.jl | 125 ------ parcel/Parcel.jl | 5 + parcel/ParcelCommon.jl | 10 + parcel/ParcelDistributions.jl | 79 ++++ parcel/ParcelModel.jl | 276 ++++++++++++ parcel/ParcelParameters.jl | 59 +++ parcel/ParcelTendencies.jl | 172 ++++++++ parcel/Project.toml | 6 - parcel/parcel.jl | 411 ------------------ 16 files changed, 960 insertions(+), 906 deletions(-) rename parcel/{Deposition_Nucleation.jl => Example_Deposition_Nucleation.jl} (63%) rename parcel/{Immersion_Freezing.jl => Example_Immersion_Freezing.jl} (60%) create mode 100644 parcel/Example_Jensen_et_al_2022.jl rename parcel/{Liquid_only.jl => Example_Liquid_only.jl} (59%) create mode 100644 parcel/Example_P3_ice_nuc.jl rename parcel/{Tully_et_al_2023.jl => Example_Tully_et_al_2023.jl} (76%) delete mode 100644 parcel/Jensen_et_al_2022.jl delete mode 100644 parcel/P3_ice_nuc.jl create mode 100644 parcel/Parcel.jl create mode 100644 parcel/ParcelCommon.jl create mode 100644 parcel/ParcelDistributions.jl create mode 100644 parcel/ParcelModel.jl create mode 100644 parcel/ParcelParameters.jl create mode 100644 parcel/ParcelTendencies.jl delete mode 100644 parcel/parcel.jl diff --git a/parcel/Deposition_Nucleation.jl b/parcel/Example_Deposition_Nucleation.jl similarity index 63% rename from parcel/Deposition_Nucleation.jl rename to parcel/Example_Deposition_Nucleation.jl index 11cb530bae..95e425d62b 100644 --- a/parcel/Deposition_Nucleation.jl +++ b/parcel/Example_Deposition_Nucleation.jl @@ -5,18 +5,11 @@ import CloudMicrophysics as CM import CLIMAParameters as CP # definition of the ODE problem for parcel model -include(joinpath(pkgdir(CM), "parcel", "parcel.jl")) +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 * 1e3) @@ -34,60 +27,45 @@ 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) +Rᵥ = TD.Parameters.R_v(tps) eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) -e = qᵥ * p₀ * R_v / Rₐ +e = eᵥ(qᵥ, p₀, Rₐ, Rᵥ) Sₗ = FT(e / eₛ) IC = [Sₗ, p₀, T₀, 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(100) -aerosol_1 = [CMP.DesertDust(FT), CMP.ArizonaTestDust(FT)] # aersol type for DustDeposition -aerosol_2 = [CMP.Feldspar(FT), CMP.Ferrihydrite(FT), CMP.Kaolinite(FT)] # aersol type for DepositionNucleation -ice_nucleation_modes_list = - [["MohlerRate_Deposition"], ["ActivityBasedDeposition"]] +aerosol_1 = [CMP.DesertDust(FT), CMP.ArizonaTestDust(FT)] # aersol type for DustDeposition +aerosol_2 = [CMP.Feldspar(FT), CMP.Ferrihydrite(FT), CMP.Kaolinite(FT)] # aersol type for DepositionNucleation +deposition_modes = ["MohlerRate", "ActivityBased"] growth_modes = ["Deposition"] -droplet_size_distribution_list = [["Monodisperse"]] +size_distribution = "Monodisperse" # Plotting -fig = MK.Figure(resolution = (800, 600)) +fig = MK.Figure(size = (800, 600)) ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation [-]") ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]") ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]", xlabel = "time") ax4 = MK.Axis(fig[2, 2], ylabel = "N_ice [m^-3]", xlabel = "time") -for ice_nucleation_modes in ice_nucleation_modes_list - nuc_mode = ice_nucleation_modes[1] - droplet_size_distribution = droplet_size_distribution_list[1] - - if nuc_mode == "MohlerRate_Deposition" +for deposition in deposition_modes + if deposition == "MohlerRate" for aerosol in aerosol_1 - p = (; - wps, - aps, - tps, - ip, - const_dt, - r_nuc, - w, - α_m, - aerosol, - ice_nucleation_modes, - growth_modes, - droplet_size_distribution, + params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + aerosol = aerosol, + deposition = deposition, + growth_modes = growth_modes, + size_distribution = size_distribution, ) - # solve ODE - sol = run_parcel(IC, FT(0), t_max, p) + sol = run_parcel(IC, FT(0), t_max, params) # Plot results if aerosol == CMP.DesertDust(FT) @@ -98,7 +76,7 @@ for ice_nucleation_modes in ice_nucleation_modes_list MK.lines!( # saturation ax1, sol.t, - S_i.(sol[3, :], sol[1, :]), + S_i.(tps, sol[3, :], sol[1, :]), label = aero_label, ) MK.lines!(ax2, sol.t, sol[3, :]) # temperature @@ -106,25 +84,19 @@ for ice_nucleation_modes in ice_nucleation_modes_list MK.lines!(ax4, sol.t, sol[9, :]) # N_ice end - elseif nuc_mode == "ActivityBasedDeposition" + elseif deposition == "ActivityBased" for aerosol in aerosol_2 - p = (; - wps, - aps, - tps, - ip, - const_dt, - r_nuc, - w, - α_m, - aerosol, - ice_nucleation_modes, - growth_modes, - droplet_size_distribution, + params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + aerosol = aerosol, + deposition = deposition, + growth_modes = growth_modes, + size_distribution = size_distribution, ) - # solve ODE - sol = run_parcel(IC, FT(0), t_max, p) + sol = run_parcel(IC, FT(0), t_max, params) # Plot results if aerosol == CMP.Feldspar(FT) @@ -137,7 +109,7 @@ for ice_nucleation_modes in ice_nucleation_modes_list MK.lines!( # saturation ax1, sol.t, - S_i.(sol[3, :], sol[1, :]), + S_i.(tps, sol[3, :], sol[1, :]), label = aero_label, linestyle = :dash, ) diff --git a/parcel/Immersion_Freezing.jl b/parcel/Example_Immersion_Freezing.jl similarity index 60% rename from parcel/Immersion_Freezing.jl rename to parcel/Example_Immersion_Freezing.jl index 2f279a812f..ab77084bee 100644 --- a/parcel/Immersion_Freezing.jl +++ b/parcel/Example_Immersion_Freezing.jl @@ -5,20 +5,14 @@ import CloudMicrophysics as CM import CLIMAParameters as CP # definition of the ODE problem for parcel model -include(joinpath(pkgdir(CM), "parcel", "parcel.jl")) +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 +ρₗ = wps.ρw Nₐ = FT(0) Nₗ = FT(500 * 1e3) Nᵢ = FT(0) @@ -31,26 +25,23 @@ qᵢ = FT(0) x_sulph = FT(0.01) # 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) +q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) +R_v = TD.Parameters.R_v(tps) Rₐ = TD.gas_constant_air(tps, q) eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) -e = qᵥ * p₀ * R_v / Rₐ - +e = eᵥ(qᵥ, p₀, Rₐ, R_v) Sₗ = FT(e / eₛ) IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] # Simulation parameters passed into ODE solver r_nuc = FT(0.5 * 1.e-4 * 1e-6) # assumed size of nucleated particles w = FT(0.7) # updraft speed -α_m = FT(0.5) # accomodation coefficient const_dt = FT(1) # model timestep t_max = FT(1200) aerosol = CMP.Illite(FT) -ice_nucleation_modes = ["ImmersionFreezing"] +heterogeneous = "ImmersionFreezing" growth_modes = ["Condensation", "Deposition"] -droplet_size_distribution_list = [["Monodisperse"], ["Gamma"]] +size_distribution_list = ["Monodisperse", "Gamma"] # Plotting fig = MK.Figure(resolution = (800, 600)) @@ -67,47 +58,26 @@ ax6 = MK.Axis( ) MK.ylims!(ax6, 1e-6, 1) -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, +for DSD in size_distribution_list + params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + aerosol = aerosol, + heterogeneous = heterogeneous, + growth_modes = growth_modes, + size_distribution = DSD, ) # solve ODE - sol = run_parcel(IC, FT(0), t_max, p) - - DSD = droplet_size_distribution[1] + sol = run_parcel(IC, FT(0), t_max, params) # Plot results - ξ(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 - 1 - - MK.lines!(ax1, sol.t / 60, S_i.(sol[3, :], sol[1, :]), label = DSD) + MK.lines!(ax1, sol.t / 60, S_i.(sol[3, :], sol[1, :]) .- 1, label = DSD) MK.lines!(ax2, sol.t / 60, sol[3, :]) MK.lines!(ax3, sol.t / 60, sol[6, :] * 1e3) MK.lines!(ax4, sol.t / 60, sol[5, :] * 1e3) - sol_Nₗ = sol[8, :] sol_Nᵢ = sol[9, :] - sol_qₗ = sol[5, :] - if DSD == "Monodisperse" - rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ) - end - if DSD == "Gamma" - λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ) - rₗ = 2 ./ λ - end MK.lines!(ax5, sol.t / 60, sol_Nₗ) MK.lines!(ax6, sol.t / 60, sol_Nᵢ ./ (sol_Nₗ .+ sol_Nᵢ)) end diff --git a/parcel/Example_Jensen_et_al_2022.jl b/parcel/Example_Jensen_et_al_2022.jl new file mode 100644 index 0000000000..730ed10604 --- /dev/null +++ b/parcel/Example_Jensen_et_al_2022.jl @@ -0,0 +1,150 @@ +import OrdinaryDiffEq as ODE +import CairoMakie as MK +import Thermodynamics as TD +import CloudMicrophysics as CM +import CLIMAParameters as CP +import Distributions as DS +import CloudMicrophysics.Parameters as CMP + +# 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) +wps = CMP.WaterProperties(FT) + +# Initial conditions +Nₐ = FT(300 * 1e6) +Nₗ = FT(0) +Nᵢ = FT(0) +r₀ = FT(25 * 1e-9) +T₀ = FT(190) +cᵥ₀ = FT(5 * 1e-6) +x_sulph = FT(0) + +# Constants +ρₗ = wps.ρw +R_v = TD.Parameters.R_v(tps) +R_d = TD.Parameters.R_d(tps) +ϵₘ = R_d / R_v +eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) +qᵥ = ϵₘ / (ϵₘ - 1 + 1 / cᵥ₀) +qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) +qᵢ = FT(0) +Sᵢ = FT(1.55) +Sₗ = Sᵢ / ξ(tps, T₀) +e = Sₗ * eₛ +p₀ = e / cᵥ₀ #TODO - ask +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 +const_dt = FT(0.01) # model timestep +t_max = FT(120) # total time +σ = FT(2) # lognormal distribution std deviation +homogeneous = "ActivityBased" # homogeneous freezing only +growth_modes = ["Deposition"] +size_distribution = "Lognormal" + +# 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(size = (1000, 1000)) +ax1 = MK.Axis(fig[1, 1], ylabel = "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]") +ax6 = MK.Axis(fig[3, 2], ylabel = "q_ice [g/kg]") + +MK.ylims!(ax2, 188.5, 190) +MK.xlims!(ax2, -5, 150) +MK.xlims!(ax3, -5, 150) +MK.xlims!(ax4, -5, 150) + +MK.lines!( + ax1, + Jensen_t_sat, + Jensen_sat, + label = "Jensen et al 2022", + color = :green, +) +MK.lines!( + ax2, + Jensen_t_T, + Jensen_T, + color = :green, + label = "Jensen et al 2022", +) +MK.lines!( + ax5, + Jensen_t_ICNC, + Jensen_ICNC, + color = :green, + label = "Jensen et al 2022", +) + +params = parcel_params{FT}( + r_nuc = r_nuc, + w = w, + const_dt = const_dt, + σ = σ, + homogeneous = homogeneous, + growth_modes = growth_modes, + size_distribution = size_distribution, +) + +# solve ODE +sol = run_parcel(IC, FT(0), t_max, params) + +# Plot results +MK.lines!( + ax1, + sol.t, + S_i.(tps, sol[3, :], (sol[1, :])), + label = "ice", + color = :blue, +) +MK.lines!(ax1, sol.t, (sol[1, :]), label = "liquid", color = :red) # liq saturation +MK.lines!(ax2, sol.t, sol[3, :], label = "CM.jl") # temperature +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, label = "CM.jl") # ICNC +MK.lines!(ax6, sol.t, sol[6, :] * 1e3) # q_ice + +MK.axislegend( + ax1, + framevisible = false, + labelsize = 12, + orientation = :horizontal, + nbanks = 2, + position = :lc, +) +MK.axislegend( + ax5, + framevisible = false, + labelsize = 12, + orientation = :horizontal, + nbanks = 2, + position = :lc, +) +MK.axislegend( + ax2, + framevisible = false, + labelsize = 12, + orientation = :horizontal, + nbanks = 2, + position = :lc, +) + +MK.save("Jensen_et_al_2022.svg", fig) diff --git a/parcel/Liquid_only.jl b/parcel/Example_Liquid_only.jl similarity index 59% rename from parcel/Liquid_only.jl rename to parcel/Example_Liquid_only.jl index 1377752f8a..85b6103014 100644 --- a/parcel/Liquid_only.jl +++ b/parcel/Example_Liquid_only.jl @@ -2,21 +2,19 @@ import OrdinaryDiffEq as ODE import CairoMakie as MK import Thermodynamics as TD import CloudMicrophysics as CM +import CloudMicrophysics.Parameters as CMP import CLIMAParameters as CP -# definition of the ODE problem for parcel model -include(joinpath(pkgdir(CM), "parcel", "parcel.jl")) +include(joinpath(pkgdir(CM), "parcel", "Parcel.jl")) FT = Float32 # Get free parameters tps = TD.Parameters.ThermodynamicsParameters(FT) wps = CMP.WaterProperties(FT) -aps = CMP.AirProperties(FT) -ip = CMP.IceNucleationParameters(FT) - # Constants ρₗ = wps.ρw +ρᵢ = wps.ρi R_v = TD.Parameters.R_v(tps) R_d = TD.Parameters.R_d(tps) @@ -28,39 +26,21 @@ r₀ = FT(8e-6) p₀ = FT(800 * 1e2) T₀ = FT(273.15 + 7.0) x_sulph = FT(0) - -eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) -e = eₛ -ϵ = R_d / R_v - -# mass per volume for dry air, vapor and liquid +e = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) +Sₗ = FT(1) md_v = (p₀ - e) / R_d / T₀ mv_v = e / R_v / T₀ ml_v = Nₗ * 4 / 3 * FT(π) * ρₗ * 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) -Sₗ = FT(e / eₛ) IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] # Simulation parameters passed into ODE solver -r_nuc = FT(0.5 * 1.e-4 * 1e-6) # assumed size of nucleated particles w = FT(10) # updraft speed -α_m = FT(0.5) # accomodation coefficient const_dt = FT(0.5) # model timestep t_max = FT(20) -aerosol = [] -ice_nucleation_modes = [] # no freezing -growth_modes = ["Condensation"] # switch on condensation -droplet_size_distribution_list = [["Monodisperse"], ["Gamma"]] +size_distribution_list = ["Monodisperse", "Gamma"] # Data from Rogers(1975) Figure 1 # https://www.tandfonline.com/doi/abs/10.1080/00046973.1975.9648397 @@ -71,35 +51,21 @@ Rogers_time_radius = [0.561, 2, 3.99, 10.7, 14.9, 19.9] Rogers_radius = [8.0, 8.08, 8.26, 8.91, 9.26, 9.68] #! format: on -fig = MK.Figure(resolution = (800, 600)) +# Setup the plots +fig = MK.Figure(size = (800, 600)) ax1 = MK.Axis(fig[1, 1], ylabel = "Supersaturation [%]") 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 = "radius [μm]") - MK.lines!(ax1, Rogers_time_supersat, Rogers_supersat, label = "Rogers_1975") MK.lines!(ax5, Rogers_time_radius, Rogers_radius) -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, - ) +for DSD in size_distribution_list + params = + parcel_params{FT}(size_distribution = DSD, const_dt = const_dt, w = w) # solve ODE - sol = run_parcel(IC, FT(0), t_max, p) - - DSD = droplet_size_distribution[1] + sol = run_parcel(IC, FT(0), t_max, params) # Plot results MK.lines!(ax1, sol.t, (sol[1, :] .- 1) * 100.0, label = DSD) @@ -108,18 +74,26 @@ for droplet_size_distribution in droplet_size_distribution_list MK.lines!(ax4, sol.t, sol[5, :] * 1e3) sol_Nₗ = sol[8, :] + sol_Nᵢ = sol[9, :] + # Compute the current air density + sol_T = sol[3, :] + sol_p = sol[2, :] + sol_qᵥ = sol[4, :] sol_qₗ = sol[5, :] - if DSD == "Monodisperse" - rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ) - end - if DSD == "Gamma" - λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ) - rₗ = 2 ./ λ + sol_qᵢ = sol[6, :] + q = TD.PhasePartition.(sol_qᵥ + sol_qₗ + sol_qᵢ, sol_qₗ, sol_qᵢ) + ts = TD.PhaseNonEquil_pTq.(tps, sol_p, sol_T, q) + ρₐ = TD.air_density.(tps, ts) + # Compute the mean particle size based on the distribution + distr = sol.prob.p.distr + moms = distribution_moments.(distr, sol_qₗ, sol_Nₗ, ρₗ, ρₐ, sol_qᵢ, Nᵢ, ρᵢ) + rₗ = similar(sol_T) + for it in range(1, length(sol_T)) + rₗ[it] = moms[it].rₗ end MK.lines!(ax5, sol.t, rₗ * 1e6) end -#fig[3, 2] = MK.Legend(fig, ax1, framevisible = false) MK.axislegend( ax1, framevisible = false, diff --git a/parcel/Example_P3_ice_nuc.jl b/parcel/Example_P3_ice_nuc.jl new file mode 100644 index 0000000000..db88c2b056 --- /dev/null +++ b/parcel/Example_P3_ice_nuc.jl @@ -0,0 +1,106 @@ +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) +wps = CM.Parameters.WaterProperties(FT) + +# Initial conditions +Nₐ = FT(2000) +Nₗ = FT(2000) +Nᵢ = FT(0) +rₗ = FT(1.25e-6) +p₀ = FT(20000) +qᵥ = FT(8.3e-4) +qₗ = FT(Nₗ * 4 / 3 * π * rₗ^3 * wps.ρw / 1.2) +qᵢ = FT(0) +x_sulph = FT(0) + +# Moisture dependent initial conditions +q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) +Rₐ = TD.gas_constant_air(tps, q) +Rᵥ = TD.Parameters.R_v(tps) +e = eᵥ(qᵥ, p₀, Rₐ, Rᵥ) + +# Simulation parameters passed into ODE solver +r_nuc = FT(1.25e-6) # assumed size of nucleated particles +w = FT(0.5) # updraft speed +const_dt = FT(0.1) # model timestep +t_max = FT(50) +growth_modes = ["Deposition"] +size_distribution = "Monodisperse" + +# Plotting +fig = MK.Figure(size = (900, 600)) +ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation [-]") +ax2 = MK.Axis(fig[1, 2], ylabel = "ICNC [cm^-3]") +ax7 = MK.Axis(fig[1, 3], ylabel = "Temperature [K]") +ax3 = MK.Axis(fig[2, 1], ylabel = "Ice Saturation [-]") +ax4 = MK.Axis(fig[2, 2], ylabel = "ICNC [cm^-3]") +ax8 = MK.Axis(fig[2, 3], ylabel = "Temperature [K]") +ax5 = MK.Axis(fig[3, 1], ylabel = "Ice Saturation [-]", xlabel = "Time [s]") +ax6 = MK.Axis(fig[3, 2], ylabel = "ICNC [cm^-3]", xlabel = "Time [s]") +ax9 = MK.Axis(fig[3, 3], ylabel = "Temperature [K]", xlabel = "Time [s]") + +ice_nucleation_modes_list = ["P3_dep", "P3_het", "P3_hom"] +T₀_list = [FT(235), FT(235), FT(233.2)] +color = [:blue, :red, :orange] +ax_row = [[ax1, ax2, ax7], [ax3, ax4, ax8], [ax5, ax6, ax9]] + +for it in [1, 2, 3] + mode = ice_nucleation_modes_list[it] + T₀ = T₀_list[it] + + ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q) + ρₐ = TD.air_density(tps, ts) + eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) + Sₗ = FT(e / eₛ) + IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] + + #! format: off + if mode == "P3_dep" + params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + deposition = mode, + growth_modes = growth_modes, + size_distribution = size_distribution, + ) + elseif mode == "P3_het" + params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + heterogeneous = mode, + growth_modes = growth_modes, + size_distribution = size_distribution, + ) + elseif mode == "P3_hom" + params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + homogeneous = mode, + growth_modes = growth_modes, + size_distribution = size_distribution, + ) + end + # solve ODE + sol = run_parcel(IC, FT(0), t_max, params) + + MK.lines!(ax_row[it][1], sol.t, S_i.(tps, sol[3, :], (sol[1, :])), label = mode, color = color[it]) # saturation + MK.lines!(ax_row[it][2], sol.t, sol[9, :] * 1e-6, color = color[it]) # ICNC + MK.lines!(ax_row[it][3], sol.t, sol[3, :], color = color[it]) # Temperature + MK.axislegend(ax_row[it][1], framevisible = false, labelsize = 12, orientation = :horizontal, position = :rt) + + #! format: on +end + +MK.save("P3_ice_nuc.svg", fig) diff --git a/parcel/Tully_et_al_2023.jl b/parcel/Example_Tully_et_al_2023.jl similarity index 76% rename from parcel/Tully_et_al_2023.jl rename to parcel/Example_Tully_et_al_2023.jl index 73df491f02..fec58ab21a 100644 --- a/parcel/Tully_et_al_2023.jl +++ b/parcel/Example_Tully_et_al_2023.jl @@ -6,7 +6,7 @@ import CloudMicrophysics as CM import CLIMAParameters as CP # definition of the ODE problem for parcel model -include(joinpath(pkgdir(CM), "parcel", "parcel.jl")) +include(joinpath(pkgdir(CM), "parcel", "Parcel.jl")) """ Wrapper for initial condition @@ -27,7 +27,7 @@ function get_initial_condition( R_a = TD.gas_constant_air(tps, q) R_v = TD.Parameters.R_v(tps) e_sl = TD.saturation_vapor_pressure(tps, T, TD.Liquid()) - e = q_vap * p_a * R_v / R_a + e = eᵥ(q_vap, p_a, R_a, R_v) S_liq = e / e_sl return [S_liq, p_a, T, q_vap, q_liq, q_ice, N_aer, N_liq, N_ice, x_sulph] @@ -45,9 +45,6 @@ function Tully_et_al_2023(FT) # get free parameters tps = TD.Parameters.ThermodynamicsParameters(FT) - aps = CMP.AirProperties(FT) - wps = CMP.WaterProperties(FT) - ip = CMP.IceNucleationParameters(FT) # Initial conditions for 1st period N_aerosol = FT(2000 * 1e3) N_droplets = FT(0) @@ -69,49 +66,34 @@ function Tully_et_al_2023(FT) t_max = 30 * 60 # Simulation parameters passed into ODE solver - r_nuc = FT(0.5 * 1.e-4 * 1e-6) # assumed size of nucleated particles - w = FT(3.5 * 1e-2) # updraft speed - α_m = FT(0.5) # accomodation coefficient - const_dt = 0.1 # model timestep - aerosol = CMP.DesertDust(FT) # aerosol type - ice_nucleation_modes_list = # ice nucleation modes - [["MohlerAF_Deposition"], ["MohlerRate_Deposition"]] - growth_modes = ["Deposition"] # switch on deposition growth - droplet_size_distribution = ["Monodisperse"] + r_nuc = FT(0.5 * 1.e-4 * 1e-6) # assumed size of nucleated particles + w = FT(3.5 * 1e-2) # updraft speed + const_dt = 0.1 # model timestep + aerosol = CMP.DesertDust(FT) # aerosol type + ice_nucleation_modes = ["MohlerAF", "MohlerRate"] # ice nucleation modes + growth_modes = ["Deposition"] # switch on deposition growth + size_distribution = "Monodisperse" # Plots - fig = MK.Figure(resolution = (800, 600)) + fig = MK.Figure(size = (800, 600)) ax1 = MK.Axis(fig[1, 1], ylabel = "Supersaturation [-]") ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]") ax3 = MK.Axis(fig[2, 1], ylabel = "N act [1/dm3]", yscale = log10) ax4 = MK.Axis(fig[2, 2], ylabel = "N areo [1/dm3]") ax5 = MK.Axis(fig[3, 1], ylabel = "q_vap [g/kg]", xlabel = "Height [m]") ax6 = MK.Axis(fig[3, 2], ylabel = "q_ice [g/kg]", xlabel = "Height [m]") - - #MK.ylims!(ax1, 1.0, 1.5) MK.ylims!(ax3, 3, 2e3) - ξ(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 - 1 - - for ice_nucleation_modes in ice_nucleation_modes_list - p = (; - wps, - aps, - tps, - ip, - const_dt, - r_nuc, - w, - α_m, - aerosol, - ice_nucleation_modes, - growth_modes, - droplet_size_distribution, + for mode in ice_nucleation_modes + params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + aerosol = aerosol, + deposition = mode, + growth_modes = growth_modes, + size_distribution = size_distribution, ) - # Simulation 1 IC1 = get_initial_condition( tps, @@ -125,8 +107,8 @@ function Tully_et_al_2023(FT) N_0, x_sulph, ) - sol1 = run_parcel(IC1, 0, t_max, p) - if ice_nucleation_modes == ["MohlerAF_Deposition"] + sol1 = run_parcel(IC1, 0, t_max, params) + if mode == "MohlerAF" # Simulation 2 # (alternatively set T and take q_vap from the previous simulation) #IC2 = get_initial_condition(sol1[2, end], sol1[3, end], T2, sol1[5, end], 0.0, sol1[6, end], sol1[7, end]) @@ -143,7 +125,7 @@ function Tully_et_al_2023(FT) sol1[9, end], x_sulph, ) - sol2 = run_parcel(IC2, sol1.t[end], sol1.t[end] + t_max, p) + sol2 = run_parcel(IC2, sol1.t[end], sol1.t[end] + t_max, params) # Simulation 3 # (alternatively set T and take q_vap from the previous simulation) @@ -161,7 +143,7 @@ function Tully_et_al_2023(FT) sol2[9, end], x_sulph, ) - sol3 = run_parcel(IC3, sol2.t[end], sol2.t[end] + t_max, p) + sol3 = run_parcel(IC3, sol2.t[end], sol2.t[end] + t_max, params) # Plot results sol = [sol1, sol2, sol3] @@ -177,14 +159,14 @@ function Tully_et_al_2023(FT) MK.lines!( ax1, sol[it].t * w, - S_i.(sol[it][3, :], sol[it][1, :]), + S_i.(tps, sol[it][3, :], sol[it][1, :]) .- 1, linestyle = :dash, color = clr[it], ) end #! format: on - elseif ice_nucleation_modes == ["MohlerRate_Deposition"] + elseif mode == "MohlerRate" MK.lines!(ax1, sol1.t * w, sol1[1, :] .- 1, color = :lightblue) MK.lines!(ax2, sol1.t * w, sol1[3, :], color = :lightblue) MK.lines!(ax3, sol1.t * w, sol1[9, :] * 1e-3, color = :lightblue) @@ -195,7 +177,7 @@ function Tully_et_al_2023(FT) MK.lines!( ax1, sol1.t * w, - S_i.(sol1[3, :], sol1[1, :]), + S_i.(tps, sol1[3, :], sol1[1, :]) .- 1, linestyle = :dash, color = :lightblue, ) diff --git a/parcel/Jensen_et_al_2022.jl b/parcel/Jensen_et_al_2022.jl deleted file mode 100644 index e09ee88c19..0000000000 --- a/parcel/Jensen_et_al_2022.jl +++ /dev/null @@ -1,159 +0,0 @@ -import OrdinaryDiffEq as ODE -import CairoMakie as MK -import Thermodynamics as TD -import CloudMicrophysics as CM -import CLIMAParameters as CP -import Distributions as DS -import CloudMicrophysics.Parameters as CMP - -# 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) -wps = CMP.WaterProperties(FT) -aps = CMP.AirProperties(FT) -ip = CMP.IceNucleationParameters(FT) -H2SO4_prs = CMP.H2SO4SolutionParameters(FT) - -# Constants -ρₗ = wps.ρw -R_v = TD.Parameters.R_v(tps) -R_d = TD.Parameters.R_d(tps) -ϵₘ = R_d / R_v - -# Saturation conversion equations -ξ(T) = - TD.saturation_vapor_pressure(tps, T, TD.Liquid()) / - TD.saturation_vapor_pressure(tps, T, TD.Ice()) -S_i(T, S_liq) = @. S_liq * ξ(T) # saturation ratio - -# Initial conditions -Nₐ = FT(300 * 1e6) -Nₗ = FT(0) -Nᵢ = FT(0) -r₀ = FT(25 * 1e-9) -T₀ = FT(190) -cᵥ₀ = FT(5 * 1e-6) -x_sulph = FT(0) - -# saturation and partial pressure -eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) -qᵥ = ϵₘ / (ϵₘ - 1 + 1 / cᵥ₀) -qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) -qᵢ = FT(0) -q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) -R_a = TD.gas_constant_air(tps, q) -Sᵢ = FT(1.55) -Sₗ = Sᵢ / ξ(T₀) -e = Sₗ * eₛ -p₀ = e / cᵥ₀ - -## Moisture dependent initial conditions -ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q) - -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(120) # total time -aerosol = [] # aerosol type. fill in if immersion or deposition freezing -R_mode_liq = r_nuc # lognormal distribution mode radius -σ = FT(2) # lognormal distribution std deviation -ice_nucleation_modes = ["HomogeneousFreezing"] # homogeneous freezing only -growth_modes = ["Deposition"] -droplet_size_distribution_list = [["Jensen_Example"]] - -# 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 = (1000, 1000)) -ax1 = MK.Axis(fig[1, 1], ylabel = "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]") -ax6 = MK.Axis(fig[3, 2], ylabel = "q_ice [g/kg]") - -MK.ylims!(ax2, 188.5, 190) -MK.xlims!(ax2, -5, 150) -MK.xlims!(ax3, -5, 150) -MK.xlims!(ax4, -5, 150) - -#! format: off -MK.lines!(ax1, Jensen_t_sat, Jensen_sat, label = "Jensen et al 2022", color = :green) -MK.lines!(ax2, Jensen_t_T, Jensen_T, color = :green, label = "Jensen et al 2022") -MK.lines!(ax5, Jensen_t_ICNC, Jensen_ICNC, color = :green, label = "Jensen et al 2022") - -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, - R_mode_liq, - σ, - ) - # 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 = "ice", color = :blue) - MK.lines!(ax1, sol.t, (sol[1, :]), label = "liquid", color = :red) # liq saturation - MK.lines!(ax2, sol.t, sol[3, :], label = "CM.jl") # temperature - 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, label = "CM.jl") # ICNC - MK.lines!(ax6, sol.t, sol[6, :] * 1e3) # q_ice - - MK.axislegend( - ax1, - framevisible = false, - labelsize = 12, - orientation = :horizontal, - nbanks = 2, - position = :lc, - ) - MK.axislegend( - ax5, - framevisible = false, - labelsize = 12, - orientation = :horizontal, - nbanks = 2, - position = :lc, - ) - MK.axislegend( - ax2, - framevisible = false, - labelsize = 12, - orientation = :horizontal, - nbanks = 2, - position = :lc, - ) -#! format: on -end - -MK.save("Jensen_et_al_2022.svg", fig) diff --git a/parcel/P3_ice_nuc.jl b/parcel/P3_ice_nuc.jl deleted file mode 100644 index 032a465e5d..0000000000 --- a/parcel/P3_ice_nuc.jl +++ /dev/null @@ -1,125 +0,0 @@ -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(235) -T₀_het = FT(235) -T₀_hom = FT(233.2) -qᵥ = FT(8.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_dep = TD.PhaseNonEquil_pTq(tps, p₀, T₀_dep, q) -ts_het = TD.PhaseNonEquil_pTq(tps, p₀, T₀_het, q) -ts_hom = TD.PhaseNonEquil_pTq(tps, p₀, T₀_hom, q) -ρₐ_dep = TD.air_density(tps, ts_dep) -ρₐ_het = TD.air_density(tps, ts_het) -ρₐ_hom = TD.air_density(tps, ts_hom) -Rₐ = TD.gas_constant_air(tps, q) -eₛ_dep = TD.saturation_vapor_pressure(tps, T₀_dep, TD.Liquid()) -eₛ_het = TD.saturation_vapor_pressure(tps, T₀_het, TD.Liquid()) -eₛ_hom = TD.saturation_vapor_pressure(tps, T₀_hom, TD.Liquid()) -e = qᵥ * p₀ * R_v / Rₐ - -Sₗ_dep = FT(e / eₛ_dep) -Sₗ_het = FT(e / eₛ_het) -Sₗ_hom = FT(e / eₛ_hom) -IC_dep = [Sₗ_dep, p₀, T₀_dep, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] -IC_het = [Sₗ_het, p₀, T₀_het, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] -IC_hom = [Sₗ_hom, 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(0.5) # updraft speed -α_m = FT(0.5) # accomodation coefficient -const_dt = FT(0.1) # model timestep -t_max = FT(50) -aerosol = [] -ice_nucleation_modes_list = [["P3_Deposition"], ["P3_het"], ["P3_hom"]] -growth_modes = ["Deposition"] -droplet_size_distribution_list = [["Monodisperse"]] - -# Plotting -fig = MK.Figure(resolution = (900, 600)) -ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation [-]") -ax2 = MK.Axis(fig[1, 2], ylabel = "ICNC [cm^-3]") -ax7 = MK.Axis(fig[1, 3], ylabel = "Temperature [K]") -ax3 = MK.Axis(fig[2, 1], ylabel = "Ice Saturation [-]") -ax4 = MK.Axis(fig[2, 2], ylabel = "ICNC [cm^-3]") -ax8 = MK.Axis(fig[2, 3], ylabel = "Temperature [K]") -ax5 = MK.Axis(fig[3, 1], ylabel = "Ice Saturation [-]", xlabel = "Time [s]") -ax6 = MK.Axis(fig[3, 2], ylabel = "ICNC [cm^-3]", xlabel = "Time [s]") -ax9 = MK.Axis(fig[3, 3], ylabel = "Temperature [K]", xlabel = "Time [s]") - -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 - #! format: 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) # ICNC - MK.lines!(ax7, sol.t, sol[3, :]) # Temperature - 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, color = :red) # ICNC - MK.lines!(ax8, sol.t, sol[3, :], color = :red) # Temperature - 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, color = :orange) # ICNC - MK.lines!(ax9, sol.t, sol[3, :], color = :orange) # Temperature - MK.axislegend(ax5, framevisible = false, labelsize = 12, orientation = :horizontal, position = :lt) - end - #! format: on -end - -MK.save("P3_ice_nuc.svg", fig) diff --git a/parcel/Parcel.jl b/parcel/Parcel.jl new file mode 100644 index 0000000000..0f312dcb8c --- /dev/null +++ b/parcel/Parcel.jl @@ -0,0 +1,5 @@ +include("ParcelParameters.jl") +include("ParcelDistributions.jl") +include("ParcelCommon.jl") +include("ParcelTendencies.jl") +include("ParcelModel.jl") diff --git a/parcel/ParcelCommon.jl b/parcel/ParcelCommon.jl new file mode 100644 index 0000000000..5ec5876100 --- /dev/null +++ b/parcel/ParcelCommon.jl @@ -0,0 +1,10 @@ +# Saturation ratio over ice +ξ(tps, T) = + TD.saturation_vapor_pressure(tps, T, TD.Liquid()) / + TD.saturation_vapor_pressure(tps, T, TD.Ice()) + +# Vapour partial pressure +eᵥ(qᵥ, p_air, R_air, Rᵥ) = qᵥ * p_air * Rᵥ / R_air + +# Saturation pressure over ice given saturation pressure over liquid +S_i(tps, T, S_liq) = ξ(tps, T) * S_liq diff --git a/parcel/ParcelDistributions.jl b/parcel/ParcelDistributions.jl new file mode 100644 index 0000000000..ed277ffd3b --- /dev/null +++ b/parcel/ParcelDistributions.jl @@ -0,0 +1,79 @@ +import CloudMicrophysics.Parameters as CMP + +struct Monodisperse{FT} <: CMP.ParametersType{FT} end + +struct Gamma{FT} <: CMP.ParametersType{FT} end + +struct Lognormal{FT} <: CMP.ParametersType{FT} + r::FT + σ::FT +end + +# Size distributiom moments +function distribution_moments(::Monodisperse, qₗ, Nₗ, ρₗ, ρ_air, qᵢ, Nᵢ, ρᵢ) + FT = typeof(qₗ) + # liquid droplet + if Nₗ == FT(0) + rₗ = FT(0) + Aₗ = FT(0) + Vₗ = FT(0) + else + rₗ = cbrt(qₗ / Nₗ / FT(4 / 3 * π) / ρₗ * ρ_air) + Aₗ = 4 * FT(π) * rₗ^2 + Vₗ = FT(4 / 3 * π) * rₗ^3 + end + # ice crystals + if Nᵢ == FT(0) + rᵢ = FT(0) + else + rᵢ = cbrt(qᵢ / Nᵢ / FT(4 / 3 * π) / ρᵢ * ρ_air) + end + return (; rₗ, Aₗ, Vₗ, rᵢ) +end + +function distribution_moments(::Gamma, qₗ, Nₗ, ρₗ, ρ_air, qᵢ, Nᵢ, ρᵢ) + #integral(n(r)) = 1/λ^2 + #integral(r*n(r)) = 2/λ^3 + # = 2/λ + #integral(r^2*n(r)) = 6/λ^4 + # = 6/λ^2 + #integral(r^3*n(r)) = 24/λ^5 + # = 24/λ^3 + # liquid droplets + if Nₗ == FT(0) + rₗ = FT(0) + Aₗ = FT(0) + Vₗ = FT(0) + else + λₗ = cbrt(32 * FT(π) * Nₗ / qₗ * ρₗ / ρ_air) + rₗ = 2 / λₗ + Aₗ = 4 * FT(π) * 6 / λₗ^2 + Vₗ = 4 / 3 * FT(π) * 24 / λₗ^3 + end + # ice crystals + if Nᵢ == FT(0) + rᵢ = FT(0) + else + λᵢ = cbrt(32 * FT(π) * Nᵢ / qᵢ * ρᵢ / ρ_air) + rᵢ = 2 / λᵢ + end + return (; rₗ, Aₗ, Vₗ, rᵢ) +end + +# Lognormal r = R_mode_liq, σ +# Jensen r = r_nuc, σ +function distribution_moments(params::Lognormal, qₗ, Nₗ, ρₗ, ρ_air, qᵢ, Nᵢ, ρᵢ) + (; r, σ) = params + # TODO - how does that connect with the modeled q_liq and N_liq? + # TODO - check if r_i is non-zero if initial conditions say q_ice = 0 + # liquid droplets + rₗ = r * exp(σ) + Vₗ = 4 / 3 * FT(π) * exp((6 * log(r) + 9 * σ^2) / (2)) + # ice crystals + if Nᵢ == FT(0) + rᵢ = FT(0) + else + rᵢ = cbrt(qᵢ / Nᵢ / (4 / 3 * FT(π)) / ρᵢ * ρ_air) * exp(σ) + end + return (; rₗ, Vₗ, rᵢ) +end diff --git a/parcel/ParcelModel.jl b/parcel/ParcelModel.jl new file mode 100644 index 0000000000..676084dd6e --- /dev/null +++ b/parcel/ParcelModel.jl @@ -0,0 +1,276 @@ +import Thermodynamics as TD +import CloudMicrophysics.Parameters as CMP + +include(joinpath(pkgdir(CM), "parcel", "ParcelParameters.jl")) + +""" + Parcel simulation parameters +""" +Base.@kwdef struct parcel_params{FT} <: CMP.ParametersType{FT} + deposition = "None" + heterogeneous = "None" + homogeneous = "None" + size_distribution = "Monodisperse" + growth_modes = ["Condensation"] + aerosol = Empty{FT}() + wps = CMP.WaterProperties(FT) + aps = CMP.AirProperties(FT) + tps = TD.Parameters.ThermodynamicsParameters(FT) + ips = CMP.IceNucleationParameters(FT) + H₂SO₄ps = CMP.H2SO4SolutionParameters(FT) + const_dt = 1 + w = FT(1) + Jensen = false + r_nuc = FT(0.5 * 1.e-4 * 1e-6) + σ = FT(2) +end +#TODO - should we have r and σ for aerosol, droplets and ice? + +clip!(x, lim) = max(x, lim) +clip!(x) = clip!(x, eltype(x)(0)) + +""" + ODE problem definitions +""" +function parcel_model(dY, Y, p, t) + # Numerical precision used in the simulation + FT = eltype(Y) + # Simulation parameters + (; wps, tps, r_nuc, w, Jensen) = p + (; distr, dep_params, imm_params, hom_params, ce_params, ds_params) = p + # Y values stored in a named tuple for ease of use + state = ( + Sₗ = Y[1], + p_air = Y[2], + T = Y[3], + qᵥ = clip!(Y[4]), + qₗ = clip!(Y[5]), + qᵢ = clip!(Y[6]), + Nₐ = clip!(Y[7]), + Nₗ = clip!(Y[8]), + Nᵢ = clip!(Y[9]), + xS = Y[10], + ) + + # Constants + Rᵥ = TD.Parameters.R_v(tps) + grav = TD.Parameters.grav(tps) + ρᵢ = wps.ρi + ρₗ = wps.ρw + + # Get the state values + (; Sₗ, p_air, T, qᵥ, qₗ, qᵢ, Nₗ, Nᵢ) = state + # Get thermodynamic parameters, phase partition and create thermo state. + q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) + ts = TD.PhaseNonEquil_pTq(tps, p_air, T, q) + + # Constants and variables that depend on the moisture content + R_air = TD.gas_constant_air(tps, q) + cp_air = TD.cp_m(tps, q) + L_subl = TD.latent_heat_sublim(tps, T) + L_fus = TD.latent_heat_fusion(tps, T) + L_vap = TD.latent_heat_vapor(tps, T) + ρ_air = TD.air_density(tps, ts) + + # Adiabatic parcel coefficients + a1 = L_vap * grav / cp_air / T^2 / Rᵥ - grav / R_air / T + a2 = 1 / qᵥ + a3 = L_vap^2 / Rᵥ / T^2 / cp_air + a4 = L_vap * L_subl / Rᵥ / T^2 / cp_air + a5 = L_vap * L_fus / Rᵥ / cp_air / (T^2) + + # Mean radius, area and volume of liquid droplets and ice crystals + PSD = distribution_moments(distr, qₗ, Nₗ, ρₗ, ρ_air, qᵢ, Nᵢ, ρᵢ) + + # Deposition ice nucleation + # (All deposition parameterizations assume monodisperse aerosol size distr) + dNᵢ_dt_dep = deposition_nucleation(dep_params, state, dY) + dqᵢ_dt_dep = dNᵢ_dt_dep * 4 / 3 * FT(π) * r_nuc^3 * ρᵢ / ρ_air + + # Heterogeneous ice nucleation + dNᵢ_dt_imm = immersion_freezing(imm_params, PSD, state) + dqᵢ_dt_imm = dNᵢ_dt_imm * PSD.Vₗ * ρᵢ / ρ_air + + # Homogeneous ice nucleation + dNᵢ_dt_hom = homogeneous_freezing(hom_params, PSD, state) + dqᵢ_dt_hom = dNᵢ_dt_hom * PSD.Vₗ * ρᵢ / ρ_air + + # Condensation/evaporation + dqₗ_dt_ce = condensation(ce_params, PSD, state, ρ_air) + # Deposition/sublimation + dqᵢ_dt_ds = deposition(ds_params, PSD, state, ρ_air) + + # number concentration and ... + dNᵢ_dt = dNᵢ_dt_dep + dNᵢ_dt_imm + dNᵢ_dt_hom + dNₐ_dt = -dNᵢ_dt_dep + dNₗ_dt = -dNᵢ_dt_imm + # ... water mass budget + dqₗ_dt_v2l = dqₗ_dt_ce + dqᵢ_dt_l2i = dqᵢ_dt_imm + dqᵢ_dt_v2i = dqᵢ_dt_dep + dqᵢ_dt_ds + # special case for the Jensen example + if Jensen + dNₐ_dt -= dNᵢ_dt_hom + dqᵢ_dt_v2i += dqᵢ_dt_hom + else + dNₗ_dt -= dNᵢ_dt_hom + dqᵢ_dt_l2i += dqᵢ_dt_hom + end + + # Update the tendecies + dqᵢ_dt = dqᵢ_dt_v2i + dqᵢ_dt_l2i + dqₗ_dt = dqₗ_dt_v2l - dqᵢ_dt_l2i + dqᵥ_dt = -dqᵢ_dt - dqₗ_dt + + dSₗ_dt = + a1 * w * Sₗ - (a2 + a3) * Sₗ * dqₗ_dt_v2l - + (a2 + a4) * Sₗ * dqᵢ_dt_v2i - a5 * Sₗ * dqᵢ_dt_l2i + + dp_air_dt = -p_air * grav / R_air / T * w + + dT_dt = + -grav / cp_air * w + + L_vap / cp_air * dqₗ_dt_v2l + + L_fus / cp_air * dqᵢ_dt_l2i + + L_subl / cp_air * dqᵢ_dt_v2i + + # Set tendencies + dY[1] = dSₗ_dt # saturation ratio over liquid water + dY[2] = dp_air_dt # pressure + dY[3] = dT_dt # temperature + dY[4] = dqᵥ_dt # vapor specific humidity + dY[5] = dqₗ_dt # liquid water specific humidity + dY[6] = dqᵢ_dt # ice specific humidity + dY[7] = dNₐ_dt # number concentration of interstitial aerosol + dY[8] = dNₗ_dt # mumber concentration of droplets + dY[9] = dNᵢ_dt # number concentration of activated particles + dY[10] = FT(0) # sulphuric acid concentration +end + +""" + run_parcel(IC, t_0, t_end, pp) + +Returns the solution of an ODE probelm defined by the parcel model. + +Inputs: + - IC - Vector with initial conditions: [Sₗ, p_air, T, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, xS] + - t_0 - simulation start time + - t_end - simulation end time + - pp - struct with parcel simulation parameters + +Parcel state vector (all variables are in base SI units): + - Sₗ - saturation ratio over liquid water + - p_air - air pressure + - T - temperature + - qᵥ - vapor specific humidity + - qₗ - liquid water specific humidity + - qᵢ - ice specific humidity + - Nₐ - number concentration of interstitial aerosol + - Nₗ - number concentration of existing water droplets + - Nᵢ - concentration of activated ice crystals + - xS - percent mass sulphuric acid + +The parcel parameters struct comes with default values that can be overwritten: + - deposition - string with deposition ice nucleation parameterization choice ["None" (default), "MohlerAF", "MohlerRate", "ActivityBased", "P3_dep"] + - heterogeneous - string with heterogeneous ice nucleation parameterization choice ["None" (default), "ImmersionFreezing", "P3_het"] + - homogeneous - string with homogeneous ice nucleation parameterization choice ["None" (default), "ActivityBased", "P3_hom"] + - size_distribution - string with cloud droplet and ice crystal size disribution choice ["Monodisperse" (default), "Gamma", "Lognormal"] + - growth_modes - a vector with particle growth modes. Supported options: "Condensation", "Deposition". Default: ["Condensation",] + - aerosol - a struct with aerosol parameters required by the nucleation parameterizations, see CloudMicrophysics documentation for all the options. The default is an Empty struct. + - wps, aps, tps, ips, H₂SO₄ps - structs with additional free parameters needed by the parameterizations. By default we use the values stored in ClimaParameters.jl. See CloudMicrophysics docs for more details. + - const_dt - model timestep [s]. Parcel model is using a simple Euler timestepper. Default value is 1 s + - w - parcel vertical velocity [m/s]. Default value is 1 m/s + - Jensen - boolean to flag the special case of running example Jensen_et_al_2022 + - r_nuc - assumed size of nucleating ice crystals. Default value is 5e-11 [m] + - σ - assumed standard deviation for the aerosol particle distribution. Defaul value is 2 +""" +function run_parcel(IC, t_0, t_end, pp) + + FT = eltype(IC) + + println(" ") + println("Size distribution: ", pp.size_distribution) + if pp.size_distribution == "Monodisperse" + distr = Monodisperse{FT}() + elseif pp.size_distribution == "Gamma" + distr = Gamma{FT}() + elseif pp.size_distribution == "Lognormal" + distr = Lognormal{FT}(pp.r_nuc, pp.σ) + else + throw("Unrecognized size_distribution") + end + + println("Deposition: ", pp.deposition) + if pp.deposition == "None" + dep_params = Empty{FT}() + elseif pp.deposition == "MohlerAF" + dep_params = MohlerAF{FT}(pp.ips, pp.aerosol, pp.tps, pp.const_dt) + elseif pp.deposition == "MohlerRate" + dep_params = MohlerRate{FT}(pp.ips, pp.aerosol, pp.tps) + elseif pp.deposition == "ActivityBased" + dep_params = ActivityBased{FT}(pp.tps, pp.aerosol, pp.r_nuc) + elseif pp.deposition == "P3_dep" + dep_params = P3_dep{FT}(pp.ips, pp.const_dt) + else + throw("Unrecognized deposition mode") + end + + println("Heterogeneous: ", pp.heterogeneous) + if pp.heterogeneous == "None" + imm_params = Empty{FT}() + elseif pp.heterogeneous == "ImmersionFreezing" + imm_params = ABIFM{FT}(pp.H₂SO₄ps, pp.tps, pp.aerosol) + elseif pp.heterogeneous == "P3_het" + imm_params = P3_het{FT}(pp.ips, pp.const_dt) + else + throw("Unrecognized heterogeneous mode") + end + + println("Homogeneous: ", pp.homogeneous) + if pp.homogeneous == "None" + hom_params = Empty{FT}() + elseif pp.homogeneous == "ActivityBased" + hom_params = WaterAct{FT}(pp.tps, pp.ips, pp.Jensen) + elseif pp.homogeneous == "P3_hom" + hom_params = P3_hom{FT}(pp.const_dt) + else + throw("Unrecognized homogeneous mode") + end + + println("Growth modes: ", pp.growth_modes) + if "Condensation" in pp.growth_modes + ce_params = CondParams{FT}(pp.aps, pp.tps) + else + ce_params = Empty{FT}() + end + if "Deposition" in pp.growth_modes + ds_params = DepParams{FT}(pp.aps, pp.tps) + else + ds_params = Empty{FT}() + end + println(" ") + + # Parameters for the ODE solver + p = ( + distr = distr, + dep_params = dep_params, + imm_params = imm_params, + hom_params = hom_params, + ce_params = ce_params, + ds_params = ds_params, + wps = pp.wps, + tps = pp.tps, + r_nuc = pp.r_nuc, + w = pp.w, + Jensen = pp.Jensen, + ) + + problem = ODE.ODEProblem(parcel_model, IC, (FT(t_0), FT(t_end)), p) + return ODE.solve( + problem, + ODE.Euler(), + dt = pp.const_dt, + reltol = 100 * eps(FT), + abstol = 100 * eps(FT), + ) +end diff --git a/parcel/ParcelParameters.jl b/parcel/ParcelParameters.jl new file mode 100644 index 0000000000..068c5f54a2 --- /dev/null +++ b/parcel/ParcelParameters.jl @@ -0,0 +1,59 @@ +import CloudMicrophysics.Parameters as CMP +import Thermodynamics.Parameters as TDP + +struct Empty{FT} <: CMP.ParametersType{FT} end + +struct MohlerAF{FT} <: CMP.ParametersType{FT} + ips::CMP.ParametersType{FT} + aerosol::CMP.AerosolType{FT} + tps::TDP.ThermodynamicsParameters{FT} + const_dt::FT +end + +struct MohlerRate{FT} <: CMP.ParametersType{FT} + ips::CMP.ParametersType{FT} + aerosol::CMP.AerosolType{FT} + tps::TDP.ThermodynamicsParameters{FT} +end + +struct ActivityBased{FT} <: CMP.ParametersType{FT} + tps::TDP.ThermodynamicsParameters{FT} + aerosol::CMP.AerosolType{FT} + r_nuc::FT +end + +struct P3_dep{FT} <: CMP.ParametersType{FT} + ips::CMP.ParametersType{FT} + const_dt::FT +end + +struct ABIFM{FT} <: CMP.ParametersType{FT} + H₂SO₄ps::CMP.ParametersType{FT} + tps::TDP.ThermodynamicsParameters{FT} + aerosol::CMP.AerosolType{FT} +end + +struct P3_het{FT} <: CMP.ParametersType{FT} + ips::CMP.ParametersType{FT} + const_dt::FT +end + +struct WaterAct{FT} <: CMP.ParametersType{FT} + tps::TDP.ThermodynamicsParameters{FT} + ips::CMP.ParametersType{FT} + Jensen::Bool +end + +struct P3_hom{FT} <: CMP.ParametersType{FT} + const_dt::FT +end + +struct CondParams{FT} <: CMP.ParametersType{FT} + aps::CMP.ParametersType{FT} + tps::TDP.ThermodynamicsParameters{FT} +end + +struct DepParams{FT} <: CMP.ParametersType{FT} + aps::CMP.ParametersType{FT} + tps::TDP.ThermodynamicsParameters{FT} +end diff --git a/parcel/ParcelTendencies.jl b/parcel/ParcelTendencies.jl new file mode 100644 index 0000000000..d619d66408 --- /dev/null +++ b/parcel/ParcelTendencies.jl @@ -0,0 +1,172 @@ +import Thermodynamics as TD +import CloudMicrophysics.Common as CMO +import CloudMicrophysics.HetIceNucleation as CMI_het +import CloudMicrophysics.HomIceNucleation as CMI_hom +import CloudMicrophysics.Parameters as CMP + +function deposition_nucleation(::Empty, state, dY) + FT = eltype(state) + return FT(0) +end + +function deposition_nucleation(params::MohlerAF, state, dY) + (; ips, aerosol, const_dt, tps) = params + (; Sₗ, T, Nₐ, Nᵢ) = state + Sᵢ = ξ(tps, T) * Sₗ + FT = eltype(state) + if Sᵢ >= ips.deposition.Sᵢ_max + @warn("Supersaturation exceeds Sᵢ_max. No dust will be activated.") + end + + AF = + Sᵢ >= ips.deposition.Sᵢ_max ? FT(0) : + AF = CMI_het.dust_activated_number_fraction( + aerosol, + ips.deposition, + Sᵢ, + T, + ) + return max(FT(0), AF * Nₐ - Nᵢ) / const_dt +end + +function deposition_nucleation(params::MohlerRate, state, dY) + (; ips, aerosol, tps) = params + (; T, Nₐ, Sₗ) = state + Sᵢ = ξ(tps, T) * Sₗ + dSᵢdt = ξ(tps, T) * dY[1] + + if Sᵢ >= ips.deposition.Sᵢ_max + @warn("Supersaturation exceeds Sᵢ_max. No dust will be activated.") + end + + return Sᵢ >= ips.deposition.Sᵢ_max ? FT(0) : + CMI_het.MohlerDepositionRate( + aerosol, + ips.deposition, + Sᵢ, + T, + dSᵢdt, + Nₐ, + ) +end + +function deposition_nucleation(params::ActivityBased, state, dY) + FT = eltype(state) + (; tps, aerosol, r_nuc) = params + (; T, p_air, qᵥ, qₗ, qᵢ) = state + + q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) + Rᵥ = TD.Parameters.R_v(tps) + R_air = TD.gas_constant_air(tps, q) + e = eᵥ(qᵥ, p_air, R_air, Rᵥ) + + Δa_w = CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T) + J = CMI_het.deposition_J(aerosol, Δa_w) + A = 4 * FT(π) * r_nuc^2 + return max(FT(0), J * Nₐ * A) +end + +function deposition_nucleation(params::P3_dep, state, dY) + FT = eltype(state) + (; ips, const_dt) = params + (; T, Nᵢ) = state + Nᵢ_depo = CMI_het.P3_deposition_N_i(ips.p3, T) + return max(FT(0), (Nᵢ_depo - Nᵢ) / const_dt) +end + +function immersion_freezing(::Empty, PSD, state) + FT = eltype(state) + return FT(0) +end + +function immersion_freezing(params::ABIFM, PSD, state) + (; T, xS, p_air, qᵥ, qₗ, qᵢ, Nₗ) = state + (; H₂SO₄ps, tps, aerosol) = params + + q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) + Rᵥ = TD.Parameters.R_v(tps) + R_air = TD.gas_constant_air(tps, q) + e = eᵥ(qᵥ, p_air, R_air, Rᵥ) + + # TODO - use CLIMAParameters + Δa_w = + T > FT(185) && T < FT(235) ? + CMO.a_w_xT(H₂SO₄ps, tps, xS, T) - CMO.a_w_ice(tps, T) : + CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T) + + if T < FT(185) || T > FT(235) + @warn("Clipping Δa_w for ABIFM immersion freezing") + end + + J = CMI_het.ABIFM_J(aerosol, Δa_w) + return max(FT(0), J * Nₗ * PSD.Aₗ) +end + +function immersion_freezing(params::P3_het, PSD, state) + FT = eltype(state) + (; const_dt, ips) = params + (; T, Nₗ, Nᵢ) = state + Nᵢ_het = CMI_het.P3_het_N_i(ips.p3, T, Nₗ, PSD.Vₗ, const_dt) + return max(FT(0), (Nᵢ_het - Nᵢ) / const_dt) +end + +function homogeneous_freezing(::Empty, PSD, state) + FT = eltype(state) + return FT(0) +end + +function homogeneous_freezing(params::WaterAct, PSD, state) + FT = eltype(state) + (; tps, ips, Jensen) = params + (; T, p_air, qᵥ, qₗ, qᵢ, Nₐ, Nₗ) = state + + q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) + Rᵥ = TD.Parameters.R_v(tps) + R_air = TD.gas_constant_air(tps, q) + e = eᵥ(qᵥ, p_air, R_air, Rᵥ) + + Δa_w = CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T) + + if Δa_w > ips.homogeneous.Δa_w_max || Δa_w < ips.homogeneous.Δa_w_min + @warn("Clipping Δa_w for Homogeneous freezing") + end + + Δa_w = max(min(ips.homogeneous.Δa_w_max, Δa_w), ips.homogeneous.Δa_w_min) + J = CMI_hom.homogeneous_J(ips.homogeneous, Δa_w) + + return Jensen ? max(FT(0), J * Nₐ * PSD.Vₗ) : max(FT(0), J * Nₗ * PSD.Vₗ) +end + +function homogeneous_freezing(params::P3_hom, PSD, state) + FT = eltype(state) + (; Nₗ, T) = state + (; const_dt) = params + #TODO - use CLIMAParameters + return T < FT(233.15) && Nₗ > FT(0) ? Nₗ / const_dt : FT(0) +end + +function condensation(::Empty, PSD, state, ρ_air) + FT = eltype(state) + return FT(0) +end + +function condensation(params::CondParams, PSD, state, ρ_air) + (; aps, tps) = params + (; Sₗ, T, Nₗ) = state + Gₗ = CMO.G_func(aps, tps, T, TD.Liquid()) + return 4 * FT(π) / ρ_air * (Sₗ - 1) * Gₗ * PSD.rₗ * Nₗ +end + +function deposition(::Empty, PSD, state, ρ_air) + FT = eltype(state) + return FT(0) +end + +function deposition(params::DepParams, PSD, state, ρ_air) + (; aps, tps) = params + (; T, Sₗ, Nᵢ) = state + FT = eltype(state) + Sᵢ = ξ(tps, T) * Sₗ + Gᵢ = CMO.G_func(aps, tps, T, TD.Ice()) + return 4 * FT(π) / ρ_air * (Sᵢ - 1) * Gᵢ * PSD.rᵢ * Nᵢ +end diff --git a/parcel/Project.toml b/parcel/Project.toml index 842e4b8d4c..0c0168db26 100644 --- a/parcel/Project.toml +++ b/parcel/Project.toml @@ -1,15 +1,9 @@ [deps] CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" -CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" -CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" -KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" - -[compat] -CLIMAParameters = "0.9" diff --git a/parcel/parcel.jl b/parcel/parcel.jl deleted file mode 100644 index 56dc3b8639..0000000000 --- a/parcel/parcel.jl +++ /dev/null @@ -1,411 +0,0 @@ -import Thermodynamics as TD -import CloudMicrophysics as CM - -import CloudMicrophysics.Common as CMO -import CloudMicrophysics.HetIceNucleation as CMI_het -import CloudMicrophysics.HomIceNucleation as CMI_hom -import CloudMicrophysics.Parameters as CMP - -#! format: off -""" - ODE problem definitions -""" -function parcel_model(dY, Y, p, t) - - # Get simulation parameters - (; aps, wps, tps, ip, const_dt, r_nuc, w, α_m, aerosol) = p - (; ice_nucleation_modes, growth_modes, droplet_size_distribution) = p - if "Jensen_Example" in droplet_size_distribution - (; σ) = p - elseif "Lognormal" in droplet_size_distribution - (; R_mode_liq, σ) = p - end - # Numerical precision used in the simulation - FT = eltype(Y) - - # Our state vector - # TODO - We are solving for both q_ and supersaturation. - # We should decide if we want to solve for S (as in the cirrus box model) - # or for q_ which is closer to what ClimaAtmos is doing. - S_liq = Y[1] # saturation ratio over liquid water - p_a = Y[2] # pressure - T = Y[3] # temperature - q_vap = Y[4] # vapor specific humidity - q_liq = Y[5] # liquid water specific humidity - q_ice = Y[6] # ice specific humidity - N_aer = Y[7] # number concentration of interstitial aerosol - N_liq = Y[8] # number concentration of existing water droplets - N_ice = Y[9] # number concentration of activated ice crystals - x_sulph = Y[10] # percent mass sulphuric acid - - # Constants - R_v = TD.Parameters.R_v(tps) - grav = TD.Parameters.grav(tps) - ρ_ice = wps.ρi - ρ_liq = wps.ρw - H2SO4_prs = CMP.H2SO4SolutionParameters(FT) - - # Get thermodynamic parameters, phase partition and create thermo state. - q = TD.PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice) - ts = TD.PhaseNonEquil_pTq(tps, p_a, T, q) - - # Saturation ratio over ice - e_si = TD.saturation_vapor_pressure(tps, T, TD.Ice()) - e_sl = TD.saturation_vapor_pressure(tps, T, TD.Liquid()) - ξ = e_sl / e_si - S_i = ξ * S_liq - - # Constants and variables that depend on the moisture content - R_a = TD.gas_constant_air(tps, q) - cp_a = TD.cp_m(tps, q) - L_subl = TD.latent_heat_sublim(tps, T) - L_fus = TD.latent_heat_fusion(tps, T) - L_vap = TD.latent_heat_vapor(tps, T) - ρ_air = TD.air_density(tps, ts) - e = q_vap * p_a * R_v / R_a - - # Adiabatic parcel coefficients - a1 = L_vap * grav / cp_a / T^2 / R_v - grav / R_a / T - a2 = 1 / q_vap - a3 = L_vap^2 / R_v / T^2 / cp_a - a4 = L_vap * L_subl / R_v / T^2 / cp_a - a5 = L_vap * L_fus / R_v / cp_a / (T^2) - - # TODO - we should zero out all tendencies and augemnt them - - dN_act_dt_depo = FT(0) - dqi_dt_new_depo = FT(0) - if "MohlerAF_Deposition" in ice_nucleation_modes - if S_i >= ip.deposition.Sᵢ_max - println("Supersaturation exceeds the allowed value. No dust will be activated.") - AF = FT(0) - else - AF = CMI_het.dust_activated_number_fraction( - aerosol, - ip.deposition, - S_i, - T, - ) - end - dN_act_dt_depo = max(FT(0), AF * N_aer - N_ice) / const_dt - dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air - elseif "MohlerRate_Deposition" in ice_nucleation_modes - if S_i >= ip.deposition.Sᵢ_max - println("Supersaturation exceeds the allowed value. No dust will be activated.") - dN_act_dt_depo = FT(0) - else - dN_act_dt_depo = CMI_het.MohlerDepositionRate(aerosol, ip.deposition, S_i, T, (dY[1] * ξ), N_aer) - end - dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air - end - if "ActivityBasedDeposition" in ice_nucleation_modes - Δa_w = CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T) - J_deposition = CMI_het.deposition_J(aerosol, Δa_w) - if "Monodisperse" in droplet_size_distribution - A_aer = 4 * FT(π) * r_nuc^2 - dN_act_dt_depo = max(FT(0), J_deposition * N_aer * A_aer) - 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(ip.p3, 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) - if "ImmersionFreezing" in ice_nucleation_modes - Δ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_immersion = CMI_het.ABIFM_J(aerosol, Δa_w) - 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 - A_l = 4 * FT(π) * r_l^2 - 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(ip.p3, T, N_liq, V_l, const_dt) - 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) - if "HomogeneousFreezing" in ice_nucleation_modes - a_w_ice = CMO.a_w_ice(tps, T) - Δa_w = CMO.a_w_eT(tps, e, T) - a_w_ice - if "Monodisperse" in droplet_size_distribution - J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w) - r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air) - V_l = 4 /3 * FT(π) * 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 * FT(π) * r_l^3 * ρ_ice / ρ_air - elseif "Gamma" in droplet_size_distribution - J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w) - λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air) - #A = N_liq* λ^2 - r_l = 2 / λ - V_l = 4 / 3 * FT(π) * r_l^3 - dN_act_dt_hom = max(FT(0), J_hom * N_aer * V_l) - dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air - elseif "Lognormal" in droplet_size_distribution - J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w) - dN_act_dt_hom = max(FT(0), J_hom * N_liq * 4 / 3 * FT(π) * exp((6*log(R_mode_liq) + 9 * σ^2)/(2))) - r_l = R_mode_liq * exp(σ) - dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air - elseif "Jensen_Example" in droplet_size_distribution - Δa_w = ifelse(Δa_w < ip.homogeneous.Δa_w_min, ip.homogeneous.Δa_w_min, ifelse(Δa_w > ip.homogeneous.Δa_w_max, ip.homogeneous.Δa_w_max, Δa_w)) - J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w) - dN_act_dt_hom = max(FT(0), J_hom * N_aer * 4 / 3 * FT(π) * exp((6*log(r_nuc) + 9 * σ^2)/(2))) - r_aero = r_nuc * exp(σ) - 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 - dN_aer_dt = -dN_act_dt_depo - dN_act_dt_hom - dN_liq_dt = -dN_act_dt_immersion - else - dN_aer_dt = -dN_act_dt_depo - dN_liq_dt = -dN_act_dt_immersion - dN_act_dt_hom - end - - # Growth - dqi_dt_depo = FT(0) - if "Deposition" in growth_modes && N_ice > 0 - G_i = CMO.G_func(aps, tps, T, TD.Ice()) - if "Monodisperse" in droplet_size_distribution - # Deposition on existing crystals (assuming all are the same...) - r_i = cbrt(q_ice / N_ice / (4 / 3 * FT(π)) / ρ_ice * ρ_air) - C_i = r_i - dqi_dt_depo = 4 * FT(π) / ρ_air * (S_i - 1) * G_i * r_i * N_ice - elseif "Gamma" in droplet_size_distribution - # Condensation on existing droplets assuming n(r) = A r exp(-λr) - λ = cbrt(32 * FT(π) * N_ice / q_ice * ρ_ice / ρ_air) - #A = N_ice* λ^2 - r_i = 2 / λ - C_i = r_i - dqi_dt_depo = 4 * FT(π) / ρ_air * (S_i - 1) * G_i * r_i * N_ice - elseif "Lognormal" in droplet_size_distribution || "Jensen_Example" in droplet_size_distribution - r_i = cbrt(q_ice / N_ice / (4 / 3 * FT(π)) / ρ_ice * ρ_air) * exp(σ) - dqi_dt_depo = 4 * FT(π) / ρ_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 - end - - dql_dt_cond = FT(0) - if "Condensation" in growth_modes && N_liq > 0 - G_l = CMO.G_func(aps, tps, T, TD.Liquid()) - if "Monodisperse" in droplet_size_distribution - # Condensation on existing droplets assuming all are the same - r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air) - elseif "Gamma" in droplet_size_distribution - # Condensation on existing droplets assuming n(r) = A r exp(-λr) - λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air) - #A = N_liq* λ^2 - r_l = 2 / λ - end - dql_dt_cond = 4 * FT(π) / ρ_air * (S_liq - 1) * G_l * r_l * N_liq - end - - dq_liq_dt_vap_to_liq = dql_dt_cond # change in q_liq from vap-liq transitions - if "Jensen_Example" in droplet_size_distribution - # change in q_ice from liq-ice transitions - dq_ice_dt_liq_to_ice = dqi_dt_new_immers - # change in q_ice from vap-ice transitions - dq_ice_dt_vap_to_ice = dqi_dt_new_depo + dqi_dt_depo + dqi_dt_new_hom - else - # change in q_ice from liq-ice transitions - dq_ice_dt_liq_to_ice = dqi_dt_new_immers + dqi_dt_new_hom - # change in q_ice from vap-ice transitions - dq_ice_dt_vap_to_ice = dqi_dt_new_depo + dqi_dt_depo - end - - # Update the tendecies - dq_ice_dt = dq_ice_dt_vap_to_ice + dq_ice_dt_liq_to_ice - dq_liq_dt = dq_liq_dt_vap_to_liq - dq_ice_dt_liq_to_ice - dq_vap_dt = -dq_ice_dt - dq_liq_dt - dS_l_dt = a1 * w * S_liq - (a2 + a3) * S_liq * dq_liq_dt_vap_to_liq - (a2 + a4) * S_liq * dq_ice_dt_vap_to_ice - a5 * S_liq * dq_ice_dt_liq_to_ice - dp_a_dt = -p_a * grav / R_a / T * w - dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt_vap_to_liq + L_fus / cp_a * dq_ice_dt_liq_to_ice + L_subl / cp_a * dq_ice_dt_vap_to_ice - - # Set tendencies - dY[1] = dS_l_dt # saturation ratio over liquid water - dY[2] = dp_a_dt # pressure - dY[3] = dT_dt # temperature - dY[4] = dq_vap_dt # vapor specific humidity - dY[5] = dq_liq_dt # liquid water specific humidity - dY[6] = dq_ice_dt # ice specific humidity - dY[7] = dN_aer_dt # number concentration of interstitial aerosol - dY[8] = dN_liq_dt # mumber concentration of droplets - dY[9] = dN_ice_dt # number concentration of activated particles - dY[10] = FT(0) # sulphuric acid concentration - - # TODO - add diagnostics output (radius, S, etc) -end -#! format: on - -""" - run_parcel(IC, t_0, t_end, p) - -Returns the solution of an ODE probelm defined by the parcel model. - -Inputs: - - IC - A vector with the initial conditions for - [S_l, p_a, T, q_vap, q_liq, q_ice, N_aer, N_liq, N_ice, x_sulph] - - t_0 - simulation start time - - t_end - simulation end time - - p - a named tuple with simulation parameters. - -Initial condition contains (all in base SI units): - - S_l - saturation ratio over liquid water, - - p_a - atmospheric pressure - - T - temperature - - q_vap - water vapor specific humidity - - q_liq - cloud liquid water specific humidity - - q_ice - cloud ice specific humidity - - N_aer - aerosol number concnetration - - N_liq - cloud droplet number concnentration - - N_ice - ice crystal number concentration - - x_sulph - sulphuric acid concentration - -The named tuple p should contain: - - wps - a struct with water parameters - - aps - a struct with air parameters - - tps - a struct with thermodynamics parameters - - aerosol - a struct with aerosol parameters - - ip - a struct with ice nucleation parameters - - const_dt - simulation timestep, - - r_nuc - assumed radius of newly nucleated ice crystals, - - w - vertical velocity, - - α_m - accomodation coefficient - - ice_nucleation_modes - a vector with enabled ice nucleation paths. Possible options: ("DustDeposition",) - - growth_modes - a vector with enabled growth modes. Possible options: ("Condensation", "Deposition") - - droplet_size_distribution - a vector with assumed droplet size distribution. Possible options: ("Monodisperse", "Gamma") -""" -function run_parcel(IC, t_0, t_end, p) - #! format: off - FT = eltype(IC) - (; const_dt, ice_nucleation_modes, growth_modes, droplet_size_distribution) = p - timestepper = ODE.Tsit5() - - println(" ") - println("Running parcel model with: ") - - print("Ice nucleation modes: ") - if "MohlerAF_Deposition" in ice_nucleation_modes - print("\n Deposition on dust particles using AF ") - print("(Note that this mode only runs for monodisperse size distributions.) ") - timestepper = ODE.Euler() - end - if "MohlerRate_Deposition" in ice_nucleation_modes - print("\n Deposition nucleation based on rate eqn in Mohler 2006 ") - print("(Note that this mode only runs for monodisperse size distributions.) ") - timestepper = ODE.Euler() - end - if "ActivityBasedDeposition" in ice_nucleation_modes - 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 (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 Water activity based homogeneous freezing ") - end - if "P3_hom" in ice_nucleation_modes - print("\n P3 homogeneous freezing ") - timestepper = ODE.Euler() - end - print("\n") - - print("Growth modes: ") - if "Condensation" in growth_modes - print("\n Condensation on liquid droplets",) - print("(Note that this mode only runs for monodisperse and gamma size distributions.) ") - end - if "Condensation_DSD" in growth_modes - print("\n Condensation on liquid droplets") - end - if "Deposition" in growth_modes - print("\n Deposition on ice crystals") - end - print("\n") - - print("Size distribution modes: ") - if "Monodisperse" in droplet_size_distribution - print("\n Monodisperse size distribution") - end - if "Gamma" in droplet_size_distribution - print("\n Gamma size distribution") - end - 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 || - "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(" ") - - #! format: on - - problem = ODE.ODEProblem(parcel_model, IC, (FT(t_0), FT(t_end)), p) - sol = ODE.solve( - problem, - timestepper, - dt = const_dt, - reltol = 100 * eps(FT), - abstol = 100 * eps(FT), - ) - return sol -end