From d043a9e35edf014ac254259ae9fda7ba7f571389 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Thu, 22 Feb 2024 16:03:55 -0800 Subject: [PATCH] final version, if only Jensen was the same --- ...on.jl => Example_Deposition_Nucleation.jl} | 20 +- ...ezing.jl => Example_Immersion_Freezing.jl} | 38 +- ...l_2022.jl => Example_Jensen_et_al_2022.jl} | 30 +- ...{Liquid_only.jl => Example_Liquid_only.jl} | 53 +- parcel/Example_P3_ice_nuc.jl | 106 ++++ ...al_2023.jl => Example_Tully_et_al_2023.jl} | 17 +- parcel/P3_ice_nuc.jl | 131 ----- parcel/Parcel.jl | 5 + parcel/ParcelCommon.jl | 10 + parcel/ParcelDistributions.jl | 79 +++ parcel/ParcelModel.jl | 254 ++++++++ parcel/ParcelParameters.jl | 58 ++ parcel/ParcelTendencies.jl | 168 ++++++ parcel/parcel.jl | 552 ------------------ 14 files changed, 733 insertions(+), 788 deletions(-) rename parcel/{Deposition_Nucleation.jl => Example_Deposition_Nucleation.jl} (90%) rename parcel/{Immersion_Freezing.jl => Example_Immersion_Freezing.jl} (73%) rename parcel/{Jensen_et_al_2022.jl => Example_Jensen_et_al_2022.jl} (87%) rename parcel/{Liquid_only.jl => Example_Liquid_only.jl} (73%) create mode 100644 parcel/Example_P3_ice_nuc.jl rename parcel/{Tully_et_al_2023.jl => Example_Tully_et_al_2023.jl} (93%) 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 90% rename from parcel/Deposition_Nucleation.jl rename to parcel/Example_Deposition_Nucleation.jl index c35a5fffaa..95e425d62b 100644 --- a/parcel/Deposition_Nucleation.jl +++ b/parcel/Example_Deposition_Nucleation.jl @@ -5,16 +5,12 @@ 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) -# Constants -R_v = TD.Parameters.R_v(tps) -R_d = TD.Parameters.R_d(tps) - # Initial conditions Nₐ = FT(2000 * 1e3) Nₗ = FT(0) @@ -31,17 +27,13 @@ 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 @@ -54,7 +46,7 @@ growth_modes = ["Deposition"] 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") @@ -84,7 +76,7 @@ for deposition in deposition_modes 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 @@ -117,7 +109,7 @@ for deposition in deposition_modes 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 73% rename from parcel/Immersion_Freezing.jl rename to parcel/Example_Immersion_Freezing.jl index b2c3f09c6b..ab77084bee 100644 --- a/parcel/Immersion_Freezing.jl +++ b/parcel/Example_Immersion_Freezing.jl @@ -5,18 +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) wps = CMP.WaterProperties(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) @@ -29,13 +25,11 @@ 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] @@ -64,7 +58,7 @@ ax6 = MK.Axis( ) MK.ylims!(ax6, 1e-6, 1) -for size_distribution in size_distribution_list +for DSD in size_distribution_list params = parcel_params{FT}( const_dt = const_dt, r_nuc = r_nuc, @@ -72,34 +66,18 @@ for size_distribution in size_distribution_list aerosol = aerosol, heterogeneous = heterogeneous, growth_modes = growth_modes, - size_distribution = size_distribution, + size_distribution = DSD, ) # solve ODE sol = run_parcel(IC, FT(0), t_max, params) - DSD = size_distribution - # 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/Jensen_et_al_2022.jl b/parcel/Example_Jensen_et_al_2022.jl similarity index 87% rename from parcel/Jensen_et_al_2022.jl rename to parcel/Example_Jensen_et_al_2022.jl index fc8ca14ed1..1b1311332d 100644 --- a/parcel/Jensen_et_al_2022.jl +++ b/parcel/Example_Jensen_et_al_2022.jl @@ -7,26 +7,13 @@ import Distributions as DS import CloudMicrophysics.Parameters as CMP # 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) -# 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) @@ -36,18 +23,19 @@ T₀ = FT(190) cᵥ₀ = FT(5 * 1e-6) x_sulph = FT(0) -# saturation and partial pressure +# 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) -q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) Sᵢ = FT(1.55) -Sₗ = Sᵢ / ξ(T₀) +Sₗ = Sᵢ / ξ(tps, T₀) e = Sₗ * eₛ p₀ = e / cᵥ₀ #TODO - ask - -## Moisture dependent initial conditions IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] # Simulation parameters passed into ODE solver @@ -71,7 +59,7 @@ Jensen_t_ICNC = [0.217, 42.69, 50.02, 54.41, 58.97, 65.316, 72.477, 82.08, 92.65 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)) +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]") @@ -102,7 +90,7 @@ params = parcel_params{FT}( sol = run_parcel(IC, FT(0), t_max, params) # Plot results -MK.lines!(ax1, sol.t, S_i(sol[3, :], (sol[1, :])), label = "ice", color = :blue) +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 diff --git a/parcel/Liquid_only.jl b/parcel/Example_Liquid_only.jl similarity index 73% rename from parcel/Liquid_only.jl rename to parcel/Example_Liquid_only.jl index 0b01658a2c..41ee56ae5c 100644 --- a/parcel/Liquid_only.jl +++ b/parcel/Example_Liquid_only.jl @@ -2,19 +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) - # Constants ρₗ = wps.ρw +ρᵢ = wps.ρi R_v = TD.Parameters.R_v(tps) R_d = TD.Parameters.R_d(tps) @@ -26,25 +26,14 @@ r₀ = FT(8e-6) p₀ = FT(800 * 1e2) T₀ = FT(273.15 + 7.0) x_sulph = FT(0) - -# mass per volume for dry air, vapor and liquid -eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) -e = eₛ +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 @@ -62,27 +51,25 @@ 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 size_distribution in size_distribution_list +for DSD in size_distribution_list params = parcel_params{FT}( - size_distribution = size_distribution, + size_distribution = DSD, const_dt = const_dt, w = w, ) # solve ODE sol = run_parcel(IC, FT(0), t_max, params) - DSD = size_distribution - # Plot results MK.lines!(ax1, sol.t, (sol[1, :] .- 1) * 100.0, label = DSD) MK.lines!(ax2, sol.t, sol[3, :]) @@ -90,16 +77,26 @@ for size_distribution in 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, :] - rₗ = DSD == "Monodisperse" ? - cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ) : - DSD == "Gamma" ? - 2 ./ (cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)) : - FT(0) + 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..e8e1f6eeb1 --- /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 93% rename from parcel/Tully_et_al_2023.jl rename to parcel/Example_Tully_et_al_2023.jl index 8e5bc21a82..7a04664c9c 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] @@ -75,22 +75,15 @@ function Tully_et_al_2023(FT) 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 mode in ice_nucleation_modes params = parcel_params{FT}( const_dt = const_dt, @@ -166,7 +159,7 @@ 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], ) @@ -184,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/P3_ice_nuc.jl b/parcel/P3_ice_nuc.jl deleted file mode 100644 index 6e1dd4b161..0000000000 --- a/parcel/P3_ice_nuc.jl +++ /dev/null @@ -1,131 +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) -wps = CMP.WaterProperties(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 -const_dt = FT(0.1) # model timestep -t_max = FT(50) -ice_nucleation_modes_list = ["P3_dep", "P3_het", "P3_hom"] -growth_modes = ["Deposition"] -size_distribution = "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 mode in ice_nucleation_modes_list - #! 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, - ) - # solve ODE - sol = run_parcel(IC_dep, FT(0), t_max, params) - MK.lines!(ax1, sol.t, S_i.(sol[3, :], (sol[1, :])), label = 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 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, - ) - # solve ODE - sol = run_parcel(IC_het, FT(0), t_max, params) - MK.lines!(ax3, sol.t, S_i.(sol[3, :], (sol[1, :])), label = 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 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, - ) - # solve ODE - sol = run_parcel(IC_hom, FT(0), t_max, params) - MK.lines!(ax5, sol.t, S_i.(sol[3, :], (sol[1, :])), label = 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..a797b3b819 --- /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..ea429b3097 --- /dev/null +++ b/parcel/ParcelModel.jl @@ -0,0 +1,254 @@ +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..dea20b3136 --- /dev/null +++ b/parcel/ParcelParameters.jl @@ -0,0 +1,58 @@ +import CloudMicrophysics.Parameters as CMP + +struct Empty{FT} <: CMP.ParametersType{FT} end + +struct MohlerAF{FT} <: CMP.ParametersType{FT} + ips + aerosol + tps + const_dt::FT +end + +struct MohlerRate{FT} <: CMP.ParametersType{FT} + ips + aerosol + tps +end + +struct ActivityBased{FT} <: CMP.ParametersType{FT} + tps + aerosol + r_nuc +end + +struct P3_dep{FT} <: CMP.ParametersType{FT} + ips + const_dt +end + +struct ABIFM{FT} <: CMP.ParametersType{FT} + H₂SO₄ps + tps + aerosol +end + +struct P3_het{FT} <: CMP.ParametersType{FT} + ips + const_dt +end + +struct WaterAct{FT} <: CMP.ParametersType{FT} + tps + ips + Jensen +end + +struct P3_hom{FT} <: CMP.ParametersType{FT} + const_dt +end + +struct CondParams{FT} <: CMP.ParametersType{FT} + aps + tps +end + +struct DepParams{FT} <: CMP.ParametersType{FT} + aps + tps +end diff --git a/parcel/ParcelTendencies.jl b/parcel/ParcelTendencies.jl new file mode 100644 index 0000000000..504420b298 --- /dev/null +++ b/parcel/ParcelTendencies.jl @@ -0,0 +1,168 @@ +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 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/parcel.jl b/parcel/parcel.jl deleted file mode 100644 index 314baac988..0000000000 --- a/parcel/parcel.jl +++ /dev/null @@ -1,552 +0,0 @@ -import CLIMAParameters -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 - -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) - rₐ = FT(0.5 * 1.e-4 * 1e-6) - σ = FT(2) -end -#TODO - should we have r and σ for aerosol, droplets and ice? - -# Saturation ratio over ice -function ξ(tps, T) - e_si = TD.saturation_vapor_pressure(tps, T, TD.Ice()) - e_sl = TD.saturation_vapor_pressure(tps, T, TD.Liquid()) - return e_sl / e_si -end - -# Vapour partial pressure -function eᵥ(qᵥ, p_air, R_air, Rᵥ) - return qᵥ * p_air * Rᵥ / R_air -end - -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 - -struct Empty{FT} <: CMP.ParametersType{FT} end - -struct MohlerAF{FT} <: CMP.ParametersType{FT} - ips - aerosol - tps - const_dt::FT -end -struct MohlerRate{FT} <: CMP.ParametersType{FT} - ips - aerosol - tps -end -struct ActivityBased{FT} <: CMP.ParametersType{FT} - tps - aerosol - rₐ -end -struct P3_dep{FT} <: CMP.ParametersType{FT} - ips - const_dt -end -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 -# params = (; tps, aerosol, rₐ) -function deposition_nucleation(params::ActivityBased, state, dY) - FT = eltype(state) - (; tps, aerosol, rₐ) = 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ₐ^2 - return max(FT(0), J * Nₐ * Aₐ) -end -#params = (; ips, const_dt) -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 - -struct ABIFM{FT} <: CMP.ParametersType{FT} - H₂SO₄ps - tps - aerosol -end -struct P3_het{FT} <: CMP.ParametersType{FT} - ips - 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 - -struct WaterAct{FT} <: CMP.ParametersType{FT} - tps - ips - Jensen -end -struct P3_hom{FT} <: CMP.ParametersType{FT} - const_dt -end -function homogeneous_freezing(::Empty, PSD, state) - FT = eltype(state) - return FT(0) -end -# params = (; tps, ips) -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 FT(0) ? Nₗ / const_dt : FT(0) -end - -struct CondParams{FT} <: CMP.ParametersType{FT} - aps - tps -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 - -struct DepParams{FT} <: CMP.ParametersType{FT} - aps - tps -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 - -""" - ODE problem definitions - - Our state vector: - 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 -""" -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ᵥ = Y[4], qₗ = Y[5], - qᵢ = Y[6], Nₐ = Y[7], Nₗ = Y[8], Nᵢ = 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, 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 - - ips - a struct with ice nucleation parameters - - const_dt - simulation timestep, - - r_nuc - assumed radius of newly nucleated ice crystals, - - w - vertical velocity, - - 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, 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ₐ) - 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(" ") - - 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, rₐ=pp.rₐ, - w=pp.w, Jensen=pp.Jensen) - - problem = ODE.ODEProblem(parcel_model, IC, (FT(t_0), FT(t_end)), p) - sol = ODE.solve( - problem, - ODE.Euler(), - dt = pp.const_dt, - reltol = 100 * eps(FT), - abstol = 100 * eps(FT), - ) - return sol -end