diff --git a/docs/bibliography.bib b/docs/bibliography.bib index 1a08af0377..0d2569d593 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -158,6 +158,20 @@ @article{Heymsfield2003 year = {2003} } +@article{Jensen2022, +author = {Jensen, E. J. and Diskin, G. S. and DiGangi, J. and Woods, S. and Lawson, R. P. and Bui, T. V.}, +title = {Homogeneous Freezing Events Sampled in the Tropical Tropopause Layer}, +journal = {Journal of Geophysical Research: Atmospheres}, +volume = {127}, +number = {17}, +pages = {e2022JD036535}, +keywords = {cirrus, nucleation, freezing}, +doi = {https://doi.org/10.1029/2022JD036535}, +url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022JD036535}, +note = {e2022JD036535 2022JD036535}, +year = {2022} +} + @article{Karcher2006, author = {Kärcher, B. and Hendricks, J. and Lohmann, U.}, title = {Physically based parameterization of cirrus cloud formation for use in global atmospheric models}, diff --git a/docs/src/IceNucleationParcel0D.md b/docs/src/IceNucleationParcel0D.md index 9405b9991d..7b688306d7 100644 --- a/docs/src/IceNucleationParcel0D.md +++ b/docs/src/IceNucleationParcel0D.md @@ -158,7 +158,7 @@ For a gamma distribution of droplets ``n(r) = A \; r \; exp(-\lambda r)``, ``\bar{r} = \frac{2}{\lambda}`` where ``\lambda = \left(\frac{32 \pi N_{tot} \rho_l}{q_l \rho_a}\right)^{1/3}``. -## Deposition on dust particles +## Deposition nucleation on dust particles Similarly, for a case of a spherical ice particle growing through water vapor deposition ```math @@ -212,10 +212,28 @@ Following the water activity based immersion freezing model (ABIFM), the ABIFM d where ``N_{liq}`` is total number of ice nuclei containing droplets and ``A`` is surface area of those droplets. +## Homogeneous Freezing +Homogeneous freezing follows the water-activity based model described in the + ``Ice Nucleation`` section which gives a nucleation rate coefficient of + units ``cm^{-3}s^{-1}``. +The ice production rate from homogeneous freezing can then be determined: +```math +\begin{equation} + P_{ice} = [\frac{dN_i}{dt}]_{hom} = J_{hom}V(N_{liq}) + \label{eq:hom_P_ice} +\end{equation} +``` +where ``N_{liq}`` is total number of ice nuclei containing droplets and + ``V`` is the volume of those droplets. + ## Example figures -Here we show example simulation results from the adiabatic parcel - model with deposition freezing on dust. +Here we show various example simulation results from the adiabatic parcel + model. This includes examples with deposition nucleation on dust, + liquid processes only, immersion freezing with condensation and deposition growth, + and homogeneous freezing with deposition growth. + +We start with deposition freezing on dust. The model is run three times for 30 minutes simulation time, (shown by three different colors on the plot). Between each run the water vapor specific humidity is changed, @@ -247,3 +265,16 @@ The plots below are the results of the adiabatic parcel model include("../../parcel/Immersion_Freezing.jl") ``` ![](immersion_freezing.svg) + +The following plots show the parcel model running with homogeneous freezing and + depositional growth assuming a lognormal distribution of aerosols. + It is compared against [Jensen2022](@cite). Note that running with the initial + conditions described in [Jensen2022](@cite) results in a ``\Delta a_w`` smaller + than the minimum valid value for the ``J_{hom}`` parameterization. We have forced + the ``\Delta a_w`` to be equal to the minimum valid value in this example only + for demonstrative purposes. + +```@example +include("../../parcel/Jensen_et_al_2022.jl") +``` +![](Jensen_et_al_2022.svg) diff --git a/parcel/Jensen_et_al_2022.jl b/parcel/Jensen_et_al_2022.jl new file mode 100644 index 0000000000..359f98fd57 --- /dev/null +++ b/parcel/Jensen_et_al_2022.jl @@ -0,0 +1,156 @@ +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 = Float64 + +# 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) + +# 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) +p₀ = FT(121 * 1e2) +T₀ = FT(190) +x_sulph = FT(0) + +# saturation and partial pressure +eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) +qᵥ = FT(5e-6) * FT(0.6) / FT(1.2) # ppmv converted to mass ppm +qₗ = Nₗ * 4 / 3 * π * r₀^3 * ρₗ / 1.2 +qᵢ = FT(0) +q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) +R_a = TD.gas_constant_air(tps, q) +e = qᵥ * p₀ * R_v / R_a +sₗ = e / eₛ + +## 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/Project.toml b/parcel/Project.toml index 5fbddd66e4..2a1dcf6247 100644 --- a/parcel/Project.toml +++ b/parcel/Project.toml @@ -4,11 +4,9 @@ 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.7.12" diff --git a/parcel/parcel.jl b/parcel/parcel.jl index 5b91f20a73..f2a8724832 100644 --- a/parcel/parcel.jl +++ b/parcel/parcel.jl @@ -3,7 +3,9 @@ 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 +import Random as RD #! format: off """ @@ -14,6 +16,11 @@ 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) @@ -21,7 +28,7 @@ function parcel_model(dY, Y, p, t) # 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 + 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 @@ -30,7 +37,7 @@ function parcel_model(dY, Y, p, t) 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 + x_sulph = Y[10] # percent mass sulphuric acid # Constants R_v = TD.Parameters.R_v(tps) @@ -47,25 +54,25 @@ function parcel_model(dY, Y, p, t) 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 + 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 + a2 = R_v * p_a / e_sl / R_a 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 - # TODO - add immersion, homogeneous, ... - #dN_act_dt_homogeneous = FT(0) dN_act_dt_depo = FT(0) dqi_dt_new_depo = FT(0) @@ -73,7 +80,7 @@ function parcel_model(dY, Y, p, t) AF = CMI_het.dust_activated_number_fraction( aerosol, ip.deposition, - S_i, + s_i, T, ) dN_act_dt_depo = max(FT(0), AF * N_aer - N_ice) / const_dt @@ -87,64 +94,129 @@ function parcel_model(dY, Y, p, t) 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 && "ImmersionFreezing" in ice_nucleation_modes + if "Monodisperse" in droplet_size_distribution r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air) - A_aer = 4 * π * r_l^2 - dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_aer) - dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air + elseif "Gamma" in droplet_size_distribution + λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air) + #A = N_liq* λ^2 + r_l = 2 / λ end - if "Gamma" in droplet_size_distribution && "ImmersionFreezing" in ice_nucleation_modes + A_aer = 4 * π * r_l^2 + dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_aer) + dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * π * 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 * π) / ρ_liq * ρ_air) + V_l = 4 /3 * π * r_l^3 + dN_act_dt_hom = max(FT(0), J_hom * N_liq * V_l) + dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air + elseif "Gamma" in droplet_size_distribution + J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w) λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air) #A = N_liq* λ^2 r_l = 2 / λ - A_aer = 4 * π * r_l^2 - dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_aer) - dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air + V_l = 4 / 3 * π * 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 * π * 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 * π * 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 * π * r_l^3 * ρ_ice / ρ_air + elseif "Jensen_Example" in droplet_size_distribution + if Δa_w < ip.homogeneous.Δa_w_min + Δa_w = ip.homogeneous.Δa_w_min + elseif Δa_w > ip.homogeneous.Δa_w_max + Δa_w = ip.homogeneous.Δa_w_max + end + J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w) + dN_act_dt_hom = max(FT(0), J_hom * N_aer * 4 / 3 * π * exp((6*log(r_nuc) + 9 * σ^2)/(2))) + r_aero = r_nuc * exp(σ) + dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * π * r_aero^3 * ρ_ice / ρ_air end end - dN_ice_dt = dN_act_dt_depo + dN_act_dt_immersion - dN_aer_dt = -dN_act_dt_depo - dN_liq_dt = -dN_act_dt_immersion + 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 - # Deposition on existing crystals (assuming all are the same...) G_i = CMO.G_func(aps, tps, T, TD.Ice()) - r_i = cbrt(q_ice / N_ice / (4 / 3 * π) / ρ_ice * ρ_air) - C_i = r_i - dqi_dt_depo = 4 * π / ρ_air * (S_i - 1) * G_i * r_i * N_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 * π) / ρ_ice * ρ_air) + C_i = r_i + dqi_dt_depo = 4 * π / ρ_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 * π * N_ice / q_ice * ρ_ice / ρ_air) + #A = N_ice* λ^2 + r_i = 2 / λ + C_i = r_i + dqi_dt_depo = 4 * π / ρ_air * (s_i - 1) * G_i * r_i * N_ice + elseif "Lognormal" in droplet_size_distribution || "Jensen_Example" in droplet_size_distribution + r_i = cbrt(q_ice / N_ice / (4 / 3 * π) / ρ_ice * ρ_air) * exp(σ) + dqi_dt_depo = 4 * π / ρ_air * (s_i - 1) * G_i * r_i * N_ice + end # TODO - check if r_i is non-zero if initial conditions say q_ice = 0 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 - G_l = CMO.G_func(aps, tps, T, TD.Liquid()) r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air) - dql_dt_cond = 4 * π / ρ_air * (S_liq - 1) * G_l * r_l * N_liq elseif "Gamma" in droplet_size_distribution # Condensation on existing droplets assuming n(r) = A r exp(-λr) - G_l = CMO.G_func(aps, tps, T, TD.Liquid()) λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air) #A = N_liq* λ^2 r_l = 2 / λ - dql_dt_cond = 4 * π / ρ_air * (S_liq - 1) * G_l * r_l * N_liq end + dql_dt_cond = 4 * π / ρ_air * (s_liq - 1) * G_l * r_l * N_liq + end + + dq_liq_dt_lv = dql_dt_cond # from liq-vap transitions + if "Jensen_Example" in droplet_size_distribution + # from liq-ice transitions + dq_liq_dt_li = -dqi_dt_new_immers + # from ice-liq transitions + dq_ice_dt_il = -dq_liq_dt_li + # from vap-ice transitions + dq_ice_dt_iv = dqi_dt_new_depo + dqi_dt_depo + dqi_dt_new_hom + else + # from liq-ice transitions + dq_liq_dt_li = -dqi_dt_new_immers - dqi_dt_new_hom + # from ice-liq transitions + dq_ice_dt_il = -dq_liq_dt_li + # from vap-ice transitions + dq_ice_dt_iv = dqi_dt_new_depo + dqi_dt_depo end # Update the tendecies - dq_ice_dt = dqi_dt_new_depo + dqi_dt_depo + dqi_dt_new_immers - dq_liq_dt = dql_dt_cond - dqi_dt_new_immers - dS_l_dt = a1 * w * S_liq - (a2 + a3) * S_liq * dq_liq_dt - (a2 + a4) * S_liq * dq_ice_dt + dq_vap_dt = -dq_ice_dt_iv - dq_liq_dt_lv + dq_ice_dt = dq_ice_dt_iv + dq_ice_dt_il + dq_liq_dt = dq_liq_dt_lv + dq_liq_dt_li + ds_l_dt = a1 * w * s_liq - (a2 + a3 * s_liq) * dq_liq_dt_lv - (a2 + a4 * s_liq) * dq_ice_dt_iv - a5 * s_liq * dq_ice_dt_il dp_a_dt = -p_a * grav / R_a / T * w - dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt + L_subl / cp_a * dq_ice_dt - dq_vap_dt = -dq_ice_dt - dq_liq_dt + dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt_lv + L_fus / cp_a * dq_ice_dt_il + L_subl / cp_a * dq_ice_dt_iv # Set tendencies - dY[1] = dS_l_dt # saturation ratio over liquid water + 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 @@ -156,7 +228,6 @@ function parcel_model(dY, Y, p, t) dY[10] = FT(0) # sulphuric acid concentration # TODO - add diagnostics output (radius, S, etc) - end #! format: on @@ -167,13 +238,13 @@ 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] + [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, + - s_l - saturation ratio over liquid water, - p_a - atmospheric pressure - T - temperature - q_vap - water vapor specific humidity @@ -199,9 +270,9 @@ The named tuple p should contain: - 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) = p + (; const_dt, ice_nucleation_modes, growth_modes, droplet_size_distribution) = p println(" ") println("Running parcel model with: ") @@ -210,24 +281,37 @@ function run_parcel(IC, t_0, t_end, p) print("Deposition on dust particles ") end if "ImmersionFreezing" in ice_nucleation_modes - print("Immersion freezing") + print("Immersion freezing ") + end + if "HomogeneousFreezing" in ice_nucleation_modes + print("Homogeneous freezing ") end print("\n") print("Growth modes: ") if "Condensation" in growth_modes - print( - "Condensation on liquid droplets with monodisperse size distribution", - ) + print("Condensation on liquid droplets",) end if "Condensation_DSD" in growth_modes - print("Condensation on liquid droplets with gamma distribution") + print("Condensation on liquid droplets") end if "Deposition" in growth_modes - print("Deposition on ice crystals with monodisperse size distribution") + print("Deposition on ice crystals") end print("\n") + print("Size distribution modes: ") + if "Monodisperse" in droplet_size_distribution + print("Monodisperse size distribution") + end + if "Gamma" in droplet_size_distribution + print("Gamma size distribution") + end + if "Lognormal" in droplet_size_distribution + print("Lognormal size distribution") + end println(" ") + #! format: on + problem = ODE.ODEProblem(parcel_model, IC, (FT(t_0), FT(t_end)), p) sol = ODE.solve( problem, diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl index 2c80312a4b..a6778ee0fd 100644 --- a/src/IceNucleation.jl +++ b/src/IceNucleation.jl @@ -85,8 +85,8 @@ Parameterization based on Koop 2000, DOI: 10.1038/35020537. """ function homogeneous_J(ip::CMP.Koop2000, Δa_w::FT) where {FT} - @assert Δa_w > ip.Δa_w_min - @assert Δa_w < ip.Δa_w_max + @assert Δa_w >= ip.Δa_w_min + @assert Δa_w <= ip.Δa_w_max logJ::FT = ip.c₁ + ip.c₂ * Δa_w - ip.c₃ * Δa_w^2 + ip.c₄ * Δa_w^3 diff --git a/test/homogeneous_ice_nucleation_tests.jl b/test/homogeneous_ice_nucleation_tests.jl index 773cf50435..0ef8b14e52 100644 --- a/test/homogeneous_ice_nucleation_tests.jl +++ b/test/homogeneous_ice_nucleation_tests.jl @@ -36,11 +36,11 @@ function test_homogeneous_J(FT) ) # If Δa_w out of range - TT.@test_throws AssertionError("Δa_w > ip.Δa_w_min") CMH.homogeneous_J( + TT.@test_throws AssertionError("Δa_w >= ip.Δa_w_min") CMH.homogeneous_J( ip.homogeneous, Δa_w_too_small, ) - TT.@test_throws AssertionError("Δa_w < ip.Δa_w_max") CMH.homogeneous_J( + TT.@test_throws AssertionError("Δa_w <= ip.Δa_w_max") CMH.homogeneous_J( ip.homogeneous, Δa_w_too_large, )