diff --git a/docs/src/IceNucleationParcel0D.md b/docs/src/IceNucleationParcel0D.md index 5aeb800d5..5f97d8903 100644 --- a/docs/src/IceNucleationParcel0D.md +++ b/docs/src/IceNucleationParcel0D.md @@ -205,7 +205,7 @@ There are multiple ways of running deposition nucleation in the parcel. above this value will result in nucleation in a different mode. `"ActivityBasedDeposition"` will trigger a water activity based approach from [Alpert2022](@cite). In this approach, ice production rate ``P_{ice, depo}`` - is calculated from + is calculated from ```math \begin{equation} P_{ice, depo} = \left[ \frac{dN_i}{dt} \right]_{depo} = J_{depo}\;A_{aero}\;N_{aero} @@ -252,14 +252,14 @@ Here we show various example simulation results from the adiabatic parcel and homogeneous freezing with deposition growth. We start with deposition freezing on dust. -The model is run three times using the `"MohlerAF_Deposition"` approach +The model is run three times using the `"MohlerAF_Deposition"` approach for 30 minutes simulation time, (shown by three different colors on the plot). Between each run the water vapor specific humidity is changed, while keeping all other state variables the same as at the last time step of the previous run. The prescribed vertical velocity is equal to 3.5 cm/s. Supersaturation is plotted for both liquid (solid lines) and ice (dashed lines). -The pale blue line uses the `"MohlerRate_Deposition"` approach. +The pale blue line uses the `"MohlerRate_Deposition"` approach. We only run it for the first GCM timestep because the rate approach requires the change in ice saturation over time. With the discontinuous jump in saturation, the parameterization is unable to determine a proper nucleation rate. When we force @@ -268,7 +268,7 @@ The pale blue line uses the `"MohlerRate_Deposition"` approach. in the `"MohlerAF_Deposition"` approach for the first GCM timestep. ```@example -include("../../parcel/Tully_et_al_2023.jl") +include("../../parcel/Example_Tully_et_al_2023.jl") ``` ![](cirrus_box.svg) @@ -280,7 +280,7 @@ The water activity based parameterization for deposition nucleation shows is no common aerosol type between the two parameterizations. ```@example -include("../../parcel/Deposition_Nucleation.jl") +include("../../parcel/Example_Deposition_Nucleation.jl") ``` ![](deposition_nucleation.svg) @@ -289,7 +289,7 @@ In the plots below, the parcel model is ran with only condensation (no ice or fr It is compared to [Rogers1975](@cite). ```@example -include("../../parcel/Liquid_only.jl") +include("../../parcel/Example_Liquid_only.jl") ``` ![](liquid_only_parcel.svg) @@ -299,7 +299,7 @@ The plots below are the results of the adiabatic parcel model yet been validated against literature. ```@example -include("../../parcel/Immersion_Freezing.jl") +include("../../parcel/Example_Immersion_Freezing.jl") ``` ![](immersion_freezing.svg) @@ -312,13 +312,13 @@ The following plots show the parcel model running with homogeneous freezing and for demonstrative purposes. ```@example -include("../../parcel/Jensen_et_al_2022.jl") +include("../../parcel/Example_Jensen_et_al_2022.jl") ``` ![](Jensen_et_al_2022.svg) ## P3 Ice Nucleation Parameterizations The parcel also includes ice nucleation parameterizations used in - the P3 scheme as described in [MorrisonMilbrandt2015](@cite). + the P3 scheme as described in [MorrisonMilbrandt2015](@cite). Deposition nucleation is based on the ice crystal number parameterization from Cooper (1986). The heterogeneous freezing parameterization, which follows Bigg(1953) with parameters from Barklie aand Gokhale (1959), is @@ -330,6 +330,6 @@ Shown below are three separate parcel simulations for deposition nucleation, INP number, while the other two freezing modes do. Updraft velocity is set to 0.5 m/s. ```@example -include("../../parcel/P3_ice_nuc.jl") +include("../../parcel/Example_P3_ice_nuc.jl") ``` ![](P3_ice_nuc.svg) diff --git a/parcel/Deposition_Nucleation.jl b/parcel/Example_Deposition_Nucleation.jl similarity index 62% rename from parcel/Deposition_Nucleation.jl rename to parcel/Example_Deposition_Nucleation.jl index 11cb530ba..9093dd3a9 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"]] -growth_modes = ["Deposition"] -droplet_size_distribution_list = [["Monodisperse"]] +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", "ABDINM"] +deposition_growth = "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") 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, + local params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + aerosol = aerosol, + deposition = deposition, + deposition_growth = deposition_growth, + size_distribution = size_distribution, ) - # solve ODE - sol = run_parcel(IC, FT(0), t_max, p) + local 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 == "ABDINM" 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, + local params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + aerosol = aerosol, + deposition = deposition, + deposition_growth = deposition_growth, + size_distribution = size_distribution, ) - # solve ODE - sol = run_parcel(IC, FT(0), t_max, p) + local 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 55% rename from parcel/Immersion_Freezing.jl rename to parcel/Example_Immersion_Freezing.jl index 2f279a812..f931b0f71 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,29 +25,26 @@ 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"] -growth_modes = ["Condensation", "Deposition"] -droplet_size_distribution_list = [["Monodisperse"], ["Gamma"]] +heterogeneous = "ABIFM" +condensation_growth = "Condensation" +deposition_growth = "Deposition" +size_distribution_list = ["Monodisperse", "Gamma"] # Plotting -fig = MK.Figure(resolution = (800, 600)) +fig = MK.Figure(size = (800, 600)) ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Supersaturation [-]") ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]") ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]") @@ -67,47 +58,31 @@ 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 + local params = parcel_params{FT}( + const_dt = const_dt, + w = w, + aerosol = aerosol, + heterogeneous = heterogeneous, + condensation_growth = condensation_growth, + deposition_growth = deposition_growth, + size_distribution = DSD, ) # solve ODE - sol = run_parcel(IC, FT(0), t_max, p) - - DSD = droplet_size_distribution[1] + local 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.(tps, 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 @@ -118,7 +93,7 @@ MK.axislegend( labelsize = 12, orientation = :horizontal, nbanks = 2, - position = :rb, + position = :rt, ) MK.save("immersion_freezing.svg", fig) diff --git a/parcel/Example_Jensen_et_al_2022.jl b/parcel/Example_Jensen_et_al_2022.jl new file mode 100644 index 000000000..27db1a266 --- /dev/null +++ b/parcel/Example_Jensen_et_al_2022.jl @@ -0,0 +1,149 @@ +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(300 * 1e6) +Nᵢ = FT(0) +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ᵥ₀) +# Compute qₗ assuming that initially droplets are lognormally distributed +# with N(r₀, σ). We are not keeping that size distribution assumption +# in the simulation +r₀ = FT(25 * 1e-9) +σ = FT(2) +qₗ = Nₗ * FT(4 / 3 * π) * exp((6 * log(r₀) + 9 * σ^2) / (2)) +qᵢ = FT(0) +Sᵢ = FT(1.55) +Sₗ = Sᵢ / ξ(tps, T₀) +e = Sₗ * eₛ +p₀ = e / cᵥ₀ +IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] + +# Simulation parameters passed into ODE solver +w = FT(1) # updraft speed +const_dt = FT(0.01) # model timestep +t_max = FT(120) # total time +homogeneous = "ABHOM" # homogeneous freezing only +deposition_growth = "Deposition" + +# 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}( + w = w, + const_dt = const_dt, + homogeneous = homogeneous, + deposition_growth = deposition_growth, +) + +# 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 1377752f8..ceb134543 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,22 @@ 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"] +condensation_growth = "Condensation" # Data from Rogers(1975) Figure 1 # https://www.tandfonline.com/doi/abs/10.1080/00046973.1975.9648397 @@ -71,35 +52,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 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 + local params = parcel_params{FT}( + size_distribution = DSD, + condensation_growth = condensation_growth, + const_dt = const_dt, + w = w, ) # solve ODE - sol = run_parcel(IC, FT(0), t_max, p) - - DSD = droplet_size_distribution[1] + local 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 +79,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, :] + local q = TD.PhasePartition.(sol_qᵥ + sol_qₗ + sol_qᵢ, sol_qₗ, sol_qᵢ) + local ts = TD.PhaseNonEquil_pTq.(tps, sol_p, sol_T, q) + local ρₐ = 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ᵢ, ρᵢ) + local 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 000000000..3a9f0701d --- /dev/null +++ b/parcel/Example_P3_ice_nuc.jl @@ -0,0 +1,105 @@ +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) +deposition_growth = "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] + local T₀ = T₀_list[it] + local ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q) + local ρₐ = TD.air_density(tps, ts) + local eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) + local Sₗ = FT(e / eₛ) + local IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] + + #! format: off + if mode == "P3_dep" + local params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + deposition = mode, + deposition_growth = deposition_growth, + size_distribution = size_distribution, + ) + elseif mode == "P3_het" + local params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + heterogeneous = mode, + deposition_growth = deposition_growth, + size_distribution = size_distribution, + ) + elseif mode == "P3_hom" + local params = parcel_params{FT}( + const_dt = const_dt, + r_nuc = r_nuc, + w = w, + homogeneous = mode, + deposition_growth = deposition_growth, + size_distribution = size_distribution, + ) + end + # solve ODE + local 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 73df491f0..91f4fe8de 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 + deposition_growth = "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, + deposition_growth = deposition_growth, + 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 e09ee88c1..000000000 --- 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 032a465e5..000000000 --- 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 000000000..0c0d3a566 --- /dev/null +++ b/parcel/Parcel.jl @@ -0,0 +1,6 @@ +import CloudMicrophysics as CM +include(joinpath(pkgdir(CM), "parcel", "ParcelParameters.jl")) +include(joinpath(pkgdir(CM), "parcel", "ParcelDistributions.jl")) +include(joinpath(pkgdir(CM), "parcel", "ParcelCommon.jl")) +include(joinpath(pkgdir(CM), "parcel", "ParcelTendencies.jl")) +include(joinpath(pkgdir(CM), "parcel", "ParcelModel.jl")) diff --git a/parcel/ParcelCommon.jl b/parcel/ParcelCommon.jl new file mode 100644 index 000000000..5ec587610 --- /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 000000000..f4cfc7c99 --- /dev/null +++ b/parcel/ParcelDistributions.jl @@ -0,0 +1,56 @@ +import CloudMicrophysics.Parameters as CMP + +struct Monodisperse{FT} <: CMP.ParametersType{FT} end + +struct Gamma{FT} <: CMP.ParametersType{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 diff --git a/parcel/ParcelModel.jl b/parcel/ParcelModel.jl new file mode 100644 index 000000000..61a1e3522 --- /dev/null +++ b/parcel/ParcelModel.jl @@ -0,0 +1,273 @@ +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" + condensation_growth = "None" + deposition_growth = "None" + size_distribution = "Monodisperse" + 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) + r_nuc = FT(0.5 * 1.e-4 * 1e-6) +end + +""" + If negative, clip the tracers to zero when computing tendencies +""" +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) = 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 - dNᵢ_dt_hom + # ... water mass budget + dqₗ_dt_v2l = dqₗ_dt_ce + dqᵢ_dt_l2i = dqᵢ_dt_imm + dqᵢ_dt_hom + dqᵢ_dt_v2i = dqᵢ_dt_dep + dqᵢ_dt_ds + + # 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"] + - condensation_growth - string with condensation/evaporation options for cloud droplets ["None" (default), "Condensation"] + - deposition_growth - string with deposition/sublimation options for ice crystals ["None" (default), "Deposition"] + - size_distribution - string with cloud droplet and ice crystal size disribution choice ["Monodisperse" (default), "Gamma"] + - 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 + - r_nuc - assumed size of nucleating ice crystals. Default value is 5e-11 [m] +""" +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}() + else + throw("Unrecognized size distribution") + end + + println("Aerosol :", chop(string(typeof(pp.aerosol)), head = 29, tail = 9)) + + 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 == "ABDINM" + dep_params = ABDINM{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 == "ABIFM" + 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 == "ABHOM" + hom_params = ABHOM{FT}(pp.tps, pp.ips) + elseif pp.homogeneous == "P3_hom" + hom_params = P3_hom{FT}(pp.const_dt) + else + throw("Unrecognized homogeneous mode") + end + + println("Condensation growth: ", pp.condensation_growth) + if pp.condensation_growth == "None" + ce_params = Empty{FT}() + elseif pp.condensation_growth == "Condensation" + ce_params = CondParams{FT}(pp.aps, pp.tps) + else + throw("Unrecognized condensation growth mode") + end + + println("Deposition growth: ", pp.deposition_growth) + if pp.deposition_growth == "None" + ds_params = Empty{FT}() + elseif pp.deposition_growth == "Deposition" + ds_params = DepParams{FT}(pp.aps, pp.tps) + else + throw("Unrecognized deposition growth mode") + 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, + ) + + 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 000000000..ce3abf879 --- /dev/null +++ b/parcel/ParcelParameters.jl @@ -0,0 +1,58 @@ +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 ABDINM{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.H2SO4SolutionParameters{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 ABHOM{FT} <: CMP.ParametersType{FT} + tps::TDP.ThermodynamicsParameters{FT} + ips::CMP.ParametersType{FT} +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 000000000..e215da7db --- /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::ABDINM, 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 - get rif of a_w_x option + Δ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) + + 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::ABHOM, PSD, state) + FT = eltype(state) + (; tps, ips) = 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 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 842e4b8d4..0c0168db2 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 56dc3b863..000000000 --- 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