From ce0a52a0a8435913b48c53a6a265aca9ebc07979 Mon Sep 17 00:00:00 2001 From: amylu00 Date: Mon, 29 Jan 2024 17:16:07 -0800 Subject: [PATCH] adding suggestions --- docs/src/IceNucleationParcel0D.md | 8 ++-- parcel/Immersion_Freezing.jl | 6 +-- parcel/Jensen_et_al_2022.jl | 8 ++-- parcel/Liquid_only.jl | 2 +- parcel/Project.toml | 3 ++ parcel/parcel.jl | 69 ++++++++++++++----------------- 6 files changed, 47 insertions(+), 49 deletions(-) diff --git a/docs/src/IceNucleationParcel0D.md b/docs/src/IceNucleationParcel0D.md index db95c78b33..f5f0abe6fc 100644 --- a/docs/src/IceNucleationParcel0D.md +++ b/docs/src/IceNucleationParcel0D.md @@ -197,7 +197,7 @@ where: - ``N_{act}`` is the number of activated ice particles. ``N_{act}`` can be computed for example from - [activated fraction](https://clima.github.io/CloudMicrophysics.jl/previews/PR103/IceNucleation/#Activated-fraction-for-deposition-freezing-on-dust) ``f_i`` + [activated fraction](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/#Activated-fraction-for-deposition-freezing-on-dust) ``f_i`` ```math \begin{equation} N_{act} = N_{aer} f_i @@ -212,7 +212,7 @@ Following the water activity based immersion freezing model (ABIFM), the ABIFM d per second via immersion freezing can then be calculating using ```math \begin{equation} - P_{ice} = [\frac{dN_i}{dt}]_{immer} = J_{immer}A(N_{liq}) + P_{ice} = \left[ \frac{dN_i}{dt} \right]_{immer} = J_{immer}A(N_{liq}) \label{eq:ABIFM_P_ice} \end{equation} ``` @@ -221,12 +221,12 @@ where ``N_{liq}`` is total number of ice nuclei containing droplets and ## Homogeneous Freezing Homogeneous freezing follows the water-activity based model described in the - ``Ice Nucleation`` section which gives a nucleation rate coefficient of + [Ice Nucleation](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/) 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}) + P_{ice} = \left[ \frac{dN_i}{dt} \right]_{hom} = J_{hom}V(N_{liq}) \label{eq:hom_P_ice} \end{equation} ``` diff --git a/parcel/Immersion_Freezing.jl b/parcel/Immersion_Freezing.jl index 12a413f9ba..3c2630a5e4 100644 --- a/parcel/Immersion_Freezing.jl +++ b/parcel/Immersion_Freezing.jl @@ -26,7 +26,7 @@ r₀ = FT(1e-6) p₀ = FT(800 * 1e2) T₀ = FT(251) qᵥ = FT(8.1e-4) -qₗ = Nₗ * 4 / 3 * π * r₀^3 * ρₗ / 1.2 # 1.2 should be ρₐ +qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / 1.2 # 1.2 should be ρₐ qᵢ = FT(0) x_sulph = FT(0.01) @@ -102,10 +102,10 @@ for droplet_size_distribution in droplet_size_distribution_list sol_Nᵢ = sol[9, :] sol_qₗ = sol[5, :] if DSD == "Monodisperse" - rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * π) / ρₗ * ρₐ) + rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ) end if DSD == "Gamma" - λ = cbrt.(32 .* π .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ) + λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ) rₗ = 2 ./ λ end MK.lines!(ax5, sol.t / 60, sol_Nₗ) diff --git a/parcel/Jensen_et_al_2022.jl b/parcel/Jensen_et_al_2022.jl index 359f98fd57..ebdeb464e1 100644 --- a/parcel/Jensen_et_al_2022.jl +++ b/parcel/Jensen_et_al_2022.jl @@ -27,7 +27,7 @@ R_d = TD.Parameters.R_d(tps) ξ(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 +S_i(T, S_liq) = @. S_liq * ξ(T) # saturation ratio # Initial conditions Nₐ = FT(300 * 1e6) @@ -46,12 +46,12 @@ 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ₛ +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] +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 @@ -118,7 +118,7 @@ for droplet_size_distribution in droplet_size_distribution_list 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, 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 diff --git a/parcel/Liquid_only.jl b/parcel/Liquid_only.jl index ebd45ad76a..10cdd96718 100644 --- a/parcel/Liquid_only.jl +++ b/parcel/Liquid_only.jl @@ -36,7 +36,7 @@ e = eₛ # mass per volume for dry air, vapor and liquid md_v = (p₀ - e) / R_d / T₀ mv_v = e / R_v / T₀ -ml_v = Nₗ * 4 / 3 * π * ρₗ * r₀^3 +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) diff --git a/parcel/Project.toml b/parcel/Project.toml index 2a1dcf6247..3bbf3b9d6b 100644 --- a/parcel/Project.toml +++ b/parcel/Project.toml @@ -10,3 +10,6 @@ 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 38559ab666..8ed506fa24 100644 --- a/parcel/parcel.jl +++ b/parcel/parcel.jl @@ -5,7 +5,6 @@ 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 """ @@ -28,7 +27,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 @@ -54,7 +53,7 @@ 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) @@ -80,11 +79,11 @@ 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 - dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * π * r_nuc^3 * ρ_ice / ρ_air + dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air end dN_act_dt_immersion = FT(0) @@ -95,15 +94,15 @@ function parcel_model(dY, Y, p, 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 * π) / ρ_liq * ρ_air) + r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air) elseif "Gamma" in droplet_size_distribution - λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air) + λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air) #A = N_liq* λ^2 r_l = 2 / λ end - A_aer = 4 * π * r_l^2 + A_aer = 4 * FT(π) * 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 + dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air end dN_act_dt_hom = FT(0) @@ -113,33 +112,29 @@ function parcel_model(dY, Y, p, 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 + 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 * π * r_l^3 * ρ_ice / ρ_air + 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 * π * N_liq / q_liq * ρ_liq / ρ_air) + λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air) #A = N_liq* λ^2 r_l = 2 / λ - V_l = 4 / 3 * π * r_l^3 + 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 * π * r_l^3 * ρ_ice / ρ_air + 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 * π * exp((6*log(R_mode_liq) + 9 * σ^2)/(2))) + 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 * π * r_l^3 * ρ_ice / ρ_air + dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * 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 + Δ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 * π * exp((6*log(r_nuc) + 9 * σ^2)/(2))) + 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 * π * r_aero^3 * ρ_ice / ρ_air + dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_aero^3 * ρ_ice / ρ_air end end @@ -158,19 +153,19 @@ function parcel_model(dY, Y, p, t) 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 * π) / ρ_ice * ρ_air) + r_i = cbrt(q_ice / N_ice / (4 / 3 * FT(π)) / ρ_ice * ρ_air) C_i = r_i - dqi_dt_depo = 4 * π / ρ_air * (s_i - 1) * G_i * r_i * N_ice + 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 * π * N_ice / q_ice * ρ_ice / ρ_air) + λ = cbrt(32 * FT(π) * 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 + 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 * π) / ρ_ice * ρ_air) * exp(σ) - dqi_dt_depo = 4 * π / ρ_air * (s_i - 1) * G_i * r_i * N_ice + 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 @@ -180,14 +175,14 @@ function parcel_model(dY, Y, p, t) 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 * π) / ρ_liq * ρ_air) + 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 * π * N_liq / q_liq * ρ_liq / ρ_air) + λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air) #A = N_liq* λ^2 r_l = 2 / λ end - dql_dt_cond = 4 * π / ρ_air * (s_liq - 1) * G_l * r_l * N_liq + 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 @@ -207,12 +202,12 @@ function parcel_model(dY, Y, p, t) 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 + 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[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 @@ -234,13 +229,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