diff --git a/docs/src/API.md b/docs/src/API.md index d23f1db8d..e0e6450cb 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -10,6 +10,8 @@ CurrentModule = CloudMicrophysics MicrophysicsNonEq MicrophysicsNonEq.τ_relax MicrophysicsNonEq.conv_q_vap_to_q_lcl_icl +MicrophysicsNonEq.conv_q_vap_to_q_lcl_MM2015 +MicrophysicsNonEq.conv_q_vap_to_q_icl_MM2015 MicrophysicsNonEq.conv_q_vap_to_q_lcl_icl_MM2015 MicrophysicsNonEq.INP_limiter MicrophysicsNonEq.terminal_velocity diff --git a/docs/src/MicrophysicsNonEq.md b/docs/src/MicrophysicsNonEq.md index f627adc2a..90fac4271 100644 --- a/docs/src/MicrophysicsNonEq.md +++ b/docs/src/MicrophysicsNonEq.md @@ -99,7 +99,7 @@ To see this, it is necessary to use the definitions of ``\tau``, ``q_{sl}``, and ```math \begin{equation} - \tau = 4 \pi N_{tot} \bar{r}, \;\;\;\;\;\;\;\; + \tau = \left( 4 \pi N_{tot} D_v \bar{r} \right)^{-1}, \;\;\;\;\;\;\;\; q_{sl} = \frac{e_{sl}}{\rho R_v T}, \;\;\;\;\;\;\;\; D_v = \frac{K}{\rho c_p}. \end{equation} @@ -127,6 +127,41 @@ include("plots/NonEqCondEvapRate.jl") +### Liquid and ice relaxation timescales as functions of number concentration + +We also provide functionality to set the relaxation timescales of liquid and ice, ``\tau_l`` and ``\tau_i``, as functions of cloud droplet number concentrations ``N_l`` and ``N_i`` following [MorrisonGrabowski2008_supersat](@cite) and [MorrisonMilbrandt2015](@cite). First, we approximate average radius of droplets from mass and number: + +```math +\begin{equation} +\bar{r} = \left( \frac{3q_c}{4 π N(T) \rho_c} \right)^{1/3} +\end{equation} +``` + +where ``q_c`` is the mass of condensate and ``ρ_c`` is the density of the condensate (either liquid or water). Then, we calculate the relaxation timescale: + +```math +\begin{equation} +\tau = \left( 4 \pi N_c D_v \bar{r} \right)^{-1} +\end{equation} +``` + +where ``D_v`` is the thermal diffusivity and ``N_c`` is the droplet number concentration of the condensate, either ``N_l`` or ``N_i``. + +### Ice relaxation timescale as function of temperature through Frostenberg 2023 + +We can also calculate ``\tau_i`` as a function of temperature using the [Frostenberg2023](@cite) parameterization. Namely, we approximate the concentration of ice crystals as equal to the mean amount of ice nucleating particles for a given temperature using [Frostenberg2023](@cite): +```math +N_i(T) = ln(-(b \cdot T)^9 \times 10^{-9}) +``` + + ``\sigma^2`` is the variance, ``a`` and ``b`` are coefficients. The parameters defined in [Frostenberg2023](@cite) for marine data sets are ``\sigma=1.37``, ``a=1`` m``^3``, ``b=1``/C. + +Then, we use ``N_i`` to solve for ``\tau_i`` as demonstrated above. Here we show an example of how ``\tau_i`` changes with temperature for a set ``q_{icl}=1e-6``: +```@example +include("plots/NonEqFrostenberg.jl") +``` +![](taui_Frostenberg.svg) + ## Cloud condensate sedimentation diff --git a/docs/src/plots/NonEqCondEvapRate.jl b/docs/src/plots/NonEqCondEvapRate.jl index 308aa557a..98d743f51 100644 --- a/docs/src/plots/NonEqCondEvapRate.jl +++ b/docs/src/plots/NonEqCondEvapRate.jl @@ -50,10 +50,8 @@ function generate_cond_evap_rate(::HydrostaticBalance_qₗ_z) ρ = 1; title *= "ρ=$(ρ)kg/m³, " dt = 1; title *= "dt=$(dt)s, " τ_relax = 1; title *= "τ_relax=$(τ_relax)s, " - cm_params = CM.Parameters.CloudLiquid(FT) - cm_params = CM.Parameters.CloudLiquid(FT(τ_relax), cm_params.ρw, cm_params.r_eff, cm_params.N_0) # overwrite τ_relax - data = @. CMNe.conv_q_vap_to_q_lcl_icl_MM2015(cm_params, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T) + data = @. CMNe.conv_q_vap_to_q_lcl_MM2015(thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T, τ_relax) colorrange = extrema(data) @. data[iszero(data)] = NaN # set zero values to NaN, then use `nan_color=:gray` to show them as gray. These are clipped values. colormap = spliced_cmap(:blues, :reds, colorrange...; mid = 0, categorical = true, symmetrize_color_ranges = true) @@ -85,7 +83,7 @@ function generate_cond_evap_rate(::Range_qₗ_T) cm_params = CM.Parameters.CloudLiquid(FT) cm_params = CM.Parameters.CloudLiquid(FT(τ_relax), cm_params.ρw, cm_params.r_eff, cm_params.N_0) # overwrite τ_relax - data = @. CMNe.conv_q_vap_to_q_lcl_icl_MM2015(cm_params, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T') + data = @. CMNe.conv_q_vap_to_q_lcl_MM2015(thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T', τ_relax) colorrange = extrema(data) @. data[iszero(data)] = NaN # set zero values to NaN, then use `nan_color=:gray` to show them as gray. These are clipped values. colormap = spliced_cmap(:blues, :reds, colorrange...; mid = 0, categorical = true, symmetrize_color_ranges = true) diff --git a/docs/src/plots/NonEqFrostenberg.jl b/docs/src/plots/NonEqFrostenberg.jl new file mode 100644 index 000000000..c51536c1b --- /dev/null +++ b/docs/src/plots/NonEqFrostenberg.jl @@ -0,0 +1,35 @@ +import Plots as PL +import ClimaParams as CP +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.MicrophysicsNonEq as CMNe + +FT = Float32 + +# set constants +ip = CMP.Frostenberg2023(FT) +aps = CMP.AirProperties(FT) +q_icl = FT(1e-6) + +A = FT(1e-5) +override_file = Dict( + "sublimation_deposition_timescale" => + Dict("value" => A, "type" => "float"), +) +override_toml_dict = CP.create_toml_dict(FT; override_file) +ice = CMP.CloudIce(override_toml_dict) + +T_range = range(233, stop = 271, length = 500) + +τᵢ_values = [CMNe.τ_Frostenberg(ice, aps, ip, q_icl, T) for T in T_range] + +PL.plot( + T_range, + τᵢ_values, + xlabel = "T [K]", + ylabel = "τᵢ [s]", + label = "τᵢ(T)", + legend = :topright, + yaxis = :log, +) + +PL.savefig("τᵢ_Frostenberg.svg") diff --git a/parcel/Example_NonEq_Frostenberg.jl b/parcel/Example_NonEq_Frostenberg.jl new file mode 100644 index 000000000..b6bb19dd8 --- /dev/null +++ b/parcel/Example_NonEq_Frostenberg.jl @@ -0,0 +1,96 @@ +import OrdinaryDiffEq as ODE +import CairoMakie as MK + +import CloudMicrophysics as CM +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.ThermodynamicsInterface as TDI +import ClimaParams as CP + +FT = Float64 + +include(joinpath(pkgdir(CM), "parcel", "Parcel.jl")) + +τₗ = FT(10) + +@info("using τ value ", τ) +override_file = Dict( + "condensation_evaporation_timescale" => + Dict("value" => τ, "type" => "float"), + "sublimation_deposition_timescale" => + Dict("value" => τ, "type" => "float"), +) +override_toml_dict = CP.create_toml_dict(FT; override_file) +liquid = CMP.CloudLiquid(override_toml_dict) +ice = CMP.CloudIce(override_toml_dict) +@info("relaxations:", liquid.τ_relax, ice.τ_relax) + +# Get free parameters +tps = TDI.TD.Parameters.ThermodynamicsParameters(FT) +wps = CMP.WaterProperties(FT) +aps = CMP.AirProperties(FT) +ip = CMP.Frostenberg2023(FT) + +# Constants +ρₗ = wps.ρw +ρᵢ = wps.ρi +R_v = TDI.Rᵥ(tps) +R_d = TDI.Rd(tps) + +# Initial conditions +Nₐ = FT(0) +Nₗ = FT(200 * 1e6) +Nᵢ = FT(1e6) +r₀ₗ = FT(1e-6) +r₀ᵢ = FT(8e-6) +p₀ = FT(800 * 1e2) +T₀ = FT(243) +ln_INPC = FT(0) +e_sat = TDI.saturation_vapor_pressure_over_liquid(tps, T₀) +Sₗ = FT(1) +e = Sₗ * e_sat +md_v = (p₀ - e) / R_d / T₀ +mv_v = e / R_v / T₀ +ml_v = Nₗ * 4 / 3 * FT(π) * ρₗ * r₀ₗ^3 +mi_v = Nᵢ * 4 / 3 * FT(π) * ρᵢ * r₀ᵢ^3 +qᵥ = mv_v / (md_v + mv_v + ml_v + mi_v) +qₗ = ml_v / (md_v + mv_v + ml_v + mi_v) +qᵢ = mi_v / (md_v + mv_v + ml_v + mi_v) +IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, ln_INPC] + +# Simulation parameters passed into ODE solver +w = FT(1) # updraft speed +const_dt = FT(0.001) # model timestep +t_max = FT(20) +condensation_growth = "NonEq_Condensation" +deposition_growth = "NonEq_Deposition_Frostenberg" + +# Setup the plots +fig = MK.Figure(size = (800, 600)) +ax1 = MK.Axis(fig[1, 1], ylabel = "Liquid Supersaturation [%]") +ax2 = MK.Axis(fig[2, 1], ylabel = "Temperature [K]") +ax3 = MK.Axis(fig[3, 1], xlabel = "Time [s]", ylabel = "q_lcl [g/kg]") +ax4 = MK.Axis(fig[1, 2], ylabel = "Ice Supersaturation [%]") +ax5 = MK.Axis(fig[2, 2], ylabel = "q_vap [g/kg]") +ax6 = MK.Axis(fig[3, 2], xlabel = "Time [s]", ylabel = "q_icl [g/kg]") + +params = parcel_params{FT}( + condensation_growth = condensation_growth, + deposition_growth = deposition_growth, + const_dt = const_dt, + w = w, + liquid = liquid, + ice = ice, +) + +# solve ODE +sol = run_parcel(IC, FT(0), t_max, params) + +# Plot results +MK.lines!(ax1, sol.t, (sol[1, :] .- 1) * 100.0) +MK.lines!(ax2, sol.t, sol[3, :]) +MK.lines!(ax3, sol.t, sol[5, :] * 1e3) +MK.lines!(ax4, sol.t, (S_i.(tps, sol[3, :], sol[1, :]) .- 1) * 100.0) +MK.lines!(ax5, sol.t, sol[4, :] * 1e3) +MK.lines!(ax6, sol.t, sol[6, :] * 1e3) + +MK.save("noneq_frostenberg_parcel.svg", fig) diff --git a/parcel/ParcelModel.jl b/parcel/ParcelModel.jl index eb0e5c950..8e11aa2a8 100644 --- a/parcel/ParcelModel.jl +++ b/parcel/ParcelModel.jl @@ -328,6 +328,8 @@ function run_parcel(IC, t_0, t_end, pp) ds_params = NonEqDepParams_simple{FT}(pp.tps, pp.ice) elseif pp.deposition_growth == "NonEq_Deposition" ds_params = NonEqDepParams{FT}(pp.tps, pp.ice, pp.const_dt) + elseif pp.deposition_growth == "NonEq_Deposition_Frostenberg" + ds_params = NonEqDepFrostenbergParams{FT}(pp.aps, pp.tps, pp.ip, pp.ice, pp.const_dt) else throw("Unrecognized deposition growth mode") end diff --git a/parcel/ParcelParameters.jl b/parcel/ParcelParameters.jl index fb647affd..226e9f8d4 100644 --- a/parcel/ParcelParameters.jl +++ b/parcel/ParcelParameters.jl @@ -110,3 +110,11 @@ struct NonEqDepParams{FT} <: CMP.ParametersType ice::CMP.CloudIce{FT} dt::FT end + +struct NonEqDepFrostenbergParams{FT} <: CMP.ParametersType + aps::CMP.AirProperties{FT} + tps::TDP.ThermodynamicsParameters{FT} + ip::CMP.Frostenberg2023{FT} + ice::CMP.CloudIce{FT} + dt::FT +end diff --git a/parcel/ParcelTendencies.jl b/parcel/ParcelTendencies.jl index 8805f92a1..10f9a687f 100644 --- a/parcel/ParcelTendencies.jl +++ b/parcel/ParcelTendencies.jl @@ -253,7 +253,7 @@ function condensation(params::NonEqCondParams, PSD, state, ρ_air) qₜ = qᵥ + qₗ + qᵢ - cond_rate = MNE.conv_q_vap_to_q_lcl_icl_MM2015(liquid, tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T) + cond_rate = MNE.conv_q_vap_to_q_lcl_MM2015(tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T, liquid.τ_relax) # Using same limiter as ClimaAtmos for now # Not sure why, but without intermediate storing of the tendencies for the @@ -310,7 +310,31 @@ function deposition(params::NonEqDepParams, PSD, state, ρ_air) if qᵥ + qᵢ > FT(0) qₜ = qᵥ + qₗ + qᵢ - dep_rate = MNE.conv_q_vap_to_q_lcl_icl_MM2015(ice, tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T) + dep_rate = MNE.conv_q_vap_to_q_icl_MM2015(tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T, ice.τ_relax) + + # Using same limiter as ClimaAtmos for now + # Not sure why, but without intermediate storing of the tendencies for the + # if/else branch this code segfaults on julia v1.11 (works fine on v1.10) + dep_limit = min(dep_rate, limit(qᵥ, dt, 1)) + sub_limit = min(abs(dep_rate), limit(qᵢ, dt, 1)) + ret = ifelse(dep_rate > FT(0), dep_limit, -1 * sub_limit) + return ret + else + return FT(0) + end +end + +function deposition(params::NonEqDepFrostenbergParams, PSD, state, ρ_air) + FT = eltype(state) + (; T, qₗ, qᵥ, qᵢ) = state + + (; aps, tps, ip, ice, dt) = params + + if qᵥ + qᵢ > FT(0) + qₜ = qᵥ + qₗ + qᵢ + + τᵢ = MNE.τ_Frostenberg(ice, aps, ip, qᵢ, T) + dep_rate = MNE.conv_q_vap_to_q_icl_MM2015(tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T, τᵢ) # Using same limiter as ClimaAtmos for now # Not sure why, but without intermediate storing of the tendencies for the diff --git a/src/BulkMicrophysicsTendencies.jl b/src/BulkMicrophysicsTendencies.jl index 0e08b130e..39605828c 100644 --- a/src/BulkMicrophysicsTendencies.jl +++ b/src/BulkMicrophysicsTendencies.jl @@ -77,6 +77,20 @@ This unified scheme handles: """ struct Microphysics2Moment <: MicrophysicsScheme end +""" + NonEq Scheme + +Abstract type for NonEquilibrium microphysics parameterizations. Namely, +this scheme specifies how to set the relaxation timescales (τₗ, τᵢ) for +liquid condensation/evaporation and ice deposition/sublimation. +NonEq_Constant sets τₗ and τᵢ are determined by the τ_relax arguments in +lcl and icl, whereas NonEq_N calculates τₗ and τᵢ from liquid and ice +number concentration. +""" +abstract type NonEqScheme end +struct NonEq_Constant <: NonEqScheme end +struct NonEq_N <: NonEqScheme end + # --- Helper Functions --- """ @@ -100,8 +114,18 @@ end """ bulk_microphysics_tendencies( - ::Microphysics1Moment, mp, tps, - ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, + ::Microphysics1Moment, + noneq_scheme::NonEqScheme, + mp, + tps, + aps, + ρ, + T, + q_tot, + q_lcl, + q_icl, + q_rai, + q_sno, ) Compute all 1-moment microphysics tendencies in one fused call. @@ -113,6 +137,7 @@ This is a pure function of local thermodynamic state, suitable for: - Unit testing in isolation # Arguments +- `noneq_scheme` Type of NonEquilibrium scheme to use - `mp`: NamedTuple of microphysics parameters with fields: - `lcl`: CloudLiquid parameters - `icl`: CloudIce parameters @@ -148,8 +173,18 @@ This is a pure function of local thermodynamic state, suitable for: Limiters depend on timestep `dt` and should be applied by the caller after computing tendencies. """ @inline function bulk_microphysics_tendencies( - ::Microphysics1Moment, mp::CMP.Microphysics1MParams, tps, - ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, N_lcl = zero(ρ), + ::Microphysics1Moment, + noneq_scheme::NonEqScheme, + mp::CMP.Microphysics1MParams, + tps, + ρ, + T, + q_tot, + q_lcl, + q_icl, + q_rai, + q_sno, + N_lcl = zero(ρ), ) # Clamp negative inputs to zero (robustness against numerical errors) ρ = UT.clamp_to_nonneg(ρ) @@ -177,12 +212,22 @@ This is a pure function of local thermodynamic state, suitable for: # --- Cloud condensate formation (non-equilibrium) --- + if noneq_scheme isa NonEq_Constant + τₗ = lcl.τ_relax + τᵢ = icl.τ_relax + elseif noneq_scheme isa NonEq_N + τₗ = CMNonEq.τ_N(q_icl, N_lcl, lcl.ρw, aps.D_vapor) + τᵢ = CMNonEq.τ_Frostenberg(icl, aps, ip, q_icl, T) + end + # Condensation/evaporation of cloud liquid - S_lcl_cond = CMNonEq.conv_q_vap_to_q_lcl_icl_MM2015(lcl, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T) + S_lcl_cond = CMNonEq.conv_q_vap_to_q_lcl_MM2015(tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T, τₗ) dq_lcl_dt += S_lcl_cond - # Deposition/sublimation of cloud ice (INP limiter applied inside conv_q_vap_to_q_lcl_icl_MM2015) - S_icl_dep = CMNonEq.conv_q_vap_to_q_lcl_icl_MM2015(icl, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T) + # Deposition/sublimation of cloud ice + S_icl_dep = CMNonEq.conv_q_vap_to_q_icl_MM2015(tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T, τᵢ) + # No ice deposition above freezing (lack of INPs) - use ifelse to avoid GPU branching + S_icl_dep = ifelse(T > tps.T_freeze, min(S_icl_dep, zero(T)), S_icl_dep) dq_icl_dt += S_icl_dep # --- Autoconversion (cloud → precipitation) --- @@ -663,7 +708,7 @@ end # --- 2-Moment Microphysics Helper Functions --- """ - warm_rain_tendencies_2m(sb, q_lcl, q_rai, ρ, n_lcl, n_rai) + warm_rain_tendencies_2m(warm_rain, noneq_scheme, tps, T, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai) Internal helper function that computes core 2M warm rain processes: autoconversion, self-collection, accretion, and rain breakup. @@ -671,7 +716,10 @@ autoconversion, self-collection, accretion, and rain breakup. Used by both warm-only and warm+ice dispatch methods to reduce code duplication. # Arguments -- `sb`: SB2006 parameters +- `warm_rain` + - `sb`: SB2006 parameters + - `aps`: air property parameters +- `noneq_scheme` Type of NonEquilibrium scheme to use - `q_lcl`: Cloud liquid specific content (kg/kg) - `q_rai`: Rain specific content (kg/kg) - `ρ`: Air density (kg/m³) @@ -685,7 +733,7 @@ Used by both warm-only and warm+ice dispatch methods to reduce code duplication. - `dn_lcl_dt`: Cloud number tendency (1/kg/s) - `dn_rai_dt`: Rain number tendency (1/kg/s) """ -@inline function warm_rain_tendencies_2m(warm_rain, tps, T, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai) +@inline function warm_rain_tendencies_2m(warm_rain, noneq_scheme, tps, T, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai) # Unpack parameters sb = warm_rain.seifert_beheng @@ -704,7 +752,8 @@ Used by both warm-only and warm+ice dispatch methods to reduce code duplication. dn_rai_dt = zero(FT) # --- Condensation of vapor / evaporation of cloud liquid water --- - ∂ₜq_lcl_cond = CMNonEq.conv_q_vap_to_q_lcl_icl_MM2015(condevap, tps, q_tot, q_lcl, q_ice, q_rai, zero(q_ice), ρ, T) + ∂ₜq_lcl_cond = + CMNonEq.conv_q_vap_to_q_lcl_MM2015(tps, q_tot, q_lcl, q_ice, q_rai, zero(q_ice), ρ, T, condevap.τ_relax) ∂ₜn_lcl_cond = zero(∂ₜq_lcl_cond) # neglect number change from condensation/evaporation dq_lcl_dt += ∂ₜq_lcl_cond dn_lcl_dt += ∂ₜn_lcl_cond @@ -759,6 +808,7 @@ end """ bulk_microphysics_tendencies( ::Microphysics2Moment, + noneq_scheme::NonEqScheme, mp::Microphysics2MParams{WR, Nothing}, ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai, ) @@ -769,6 +819,7 @@ This method is type-stable and GPU-optimized for warm rain processes only. For warm rain + P3 ice, see the method that accepts `Microphysics2MParams{FT, WR, <:P3IceParams}`. # Arguments +- `noneq_scheme` Type of NonEquilibrium scheme to use - `mp`: Microphysics2MParams with `mp.ice == nothing` (warm rain only) - `tps`: Thermodynamics parameters - `ρ`: Air density (kg/m³) @@ -790,7 +841,7 @@ For warm rain + P3 ice, see the method that accepts `Microphysics2MParams{FT, WR - `db_rim_dt`: Rime volume tendency (always zero for warm-only) """ @inline function bulk_microphysics_tendencies( - ::Microphysics2Moment, mp::CMP.Microphysics2MParams{WR, Nothing}, tps, + ::Microphysics2Moment, noneq_scheme::NonEqScheme, mp::CMP.Microphysics2MParams{WR, Nothing}, tps, ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai, q_ice = zero(ρ), n_ice = zero(ρ), q_rim = zero(ρ), b_rim = zero(ρ), logλ = zero(ρ), ) where {WR} @@ -812,7 +863,7 @@ For warm rain + P3 ice, see the method that accepts `Microphysics2MParams{FT, WR db_rim_dt = zero(ρ) # --- Core Warm Rain Processes (shared helper) --- - warm = warm_rain_tendencies_2m(mp.warm_rain, tps, T, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai) + warm = warm_rain_tendencies_2m(mp.warm_rain, noneq_scheme, tps, T, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai) dq_lcl_dt = warm.dq_lcl_dt dn_lcl_dt = warm.dn_lcl_dt dq_rai_dt = warm.dq_rai_dt @@ -824,6 +875,7 @@ end """ bulk_microphysics_tendencies( ::Microphysics2Moment, + noneq_scheme::NonEqScheme, mp::Microphysics2MParams{FT, WR, <:P3IceParams}, ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai, q_ice, n_ice, q_rim, b_rim, logλ, @@ -836,6 +888,7 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch. # Arguments ## Required +- `noneq_scheme` Type of NonEquilibrium scheme to use - `mp`: Microphysics2MParams with P3 ice parameters present - `tps`: Thermodynamics parameters - `ρ`: Air density (kg/m³) @@ -864,7 +917,7 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch. - `db_rim_dt`: Rime volume tendency (m³/kg/s) """ @inline function bulk_microphysics_tendencies( - ::Microphysics2Moment, mp::CMP.Microphysics2MParams{WR, ICE}, tps, + ::Microphysics2Moment, noneq_scheme::NonEqScheme, mp::CMP.Microphysics2MParams{WR, ICE}, tps, ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai, q_ice = zero(ρ), n_ice = zero(ρ), q_rim = zero(ρ), b_rim = zero(ρ), logλ = zero(ρ), ) where {WR, ICE <: CMP.P3IceParams} @@ -891,7 +944,7 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch. db_rim_dt = zero(ρ) # --- Core Warm Rain Processes (shared helper) --- - warm = warm_rain_tendencies_2m(mp.warm_rain, tps, T, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai) + warm = warm_rain_tendencies_2m(mp.warm_rain, noneq_scheme, tps, T, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai) dq_lcl_dt = warm.dq_lcl_dt dn_lcl_dt = warm.dn_lcl_dt dq_rai_dt = warm.dq_rai_dt diff --git a/src/MicrophysicsNonEq.jl b/src/MicrophysicsNonEq.jl index 7ca4fbe8e..f961443ed 100644 --- a/src/MicrophysicsNonEq.jl +++ b/src/MicrophysicsNonEq.jl @@ -13,9 +13,12 @@ import ..Parameters as CMP import ..ThermodynamicsInterface as TDI import ..Common as CO import ..Utilities as UT +import ..HetIceNucleation as IN export τ_relax export conv_q_vap_to_q_lcl_icl +export conv_q_vap_to_q_lcl_MM2015 +export conv_q_vap_to_q_icl_MM2015 export conv_q_vap_to_q_lcl_icl_MM2015 export INP_limiter export dqcld_dT @@ -37,8 +40,34 @@ Returns the relaxation timescale for phase change processes. - Relaxation timescale [s] for condensation/evaporation (liquid) or deposition/sublimation (ice) """ -@inline τ_relax(p::CondEvap) = p.τ_relax -@inline τ_relax(p::SubDep) = p.τ_relax +@inline τ_relax(p::CMP.CloudLiquid) = p.τ_relax +@inline τ_relax(p::CMP.CloudIce) = p.τ_relax + +""" Calculate relaxation timescales from some number concentration N """ + +""" Includes options to calculate relaxation timescales from cloud droplet +number concentration for both liquid and ice. Also includes +an option to calculate ice relaxation timescale through +the Frostenberg et al., (2023). See DOI: 10.5194/acp-23-10883-2023 +parameterization that approximates ice droplet number from temperature. """ + +function τ_N(q_c, N_c, ρ_c, D_vapor) + r = ((3 * q_c) / (4 * pi * N_c * ρ_c))^(1 / 3) + τ = (4 * pi * D_vapor * N_c * r)^(-1) + return τ +end + +function τ_Frostenberg( + (; ρᵢ)::CMP.CloudIce, + (; D_vapor)::CMP.AirProperties, + ip::CMP.Frostenberg2023, + q_icl, + T, +) + N_ice = exp(IN.INP_concentration_mean(ip, T)) + τᵢ = τ_N(q_icl, N_ice, ρᵢ, D_vapor) + return τᵢ +end """ conv_q_vap_to_q_lcl_icl(params::CloudLiquid, q_sat_liq, q_lcl) @@ -121,7 +150,6 @@ This formulation includes a thermodynamic adjustment factor Γ that accounts for latent heat release modifying the saturation state. # Arguments -- `params` - cloud liquid or ice parameters struct containing `τ_relax` - `tps` - thermodynamics parameters struct - `q_tot` - total water specific content [kg/kg] - `q_lcl` - cloud liquid water specific content [kg/kg] @@ -130,6 +158,7 @@ accounts for latent heat release modifying the saturation state. - `q_sno` - snow specific content [kg/kg] - `ρ` - air density [kg/m³] - `T` - air temperature [K] +- `τ_relax` - relaxation timescale [s] # Returns - Cloud condensate tendency [kg/kg/s] @@ -138,9 +167,16 @@ accounts for latent heat release modifying the saturation state. This function does NOT apply limiters for small or negative specific humidities. Users should apply appropriate bounds checking when integrating in a model. """ -function conv_q_vap_to_q_lcl_icl_MM2015( - (; τ_relax)::CondEvap, tps::TDI.PS, - q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T, +function conv_q_vap_to_q_lcl_MM2015( + tps::TDI.PS, + q_tot, + q_lcl, + q_icl, + q_rai, + q_sno, + ρ, + T, + τ_relax, ) FT = eltype(tps) Rᵥ = TDI.Rᵥ(tps) @@ -163,9 +199,42 @@ function conv_q_vap_to_q_lcl_icl_MM2015( sat_excess / timescale, ) end + +# just to be able to still use old functionality in ClimaAtmos -- to be deleted function conv_q_vap_to_q_lcl_icl_MM2015( - (; τ_relax)::SubDep, tps::TDI.PS, - q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T, + liquid::CMP.CloudLiquid, + tps::TDI.PS, + q_tot, + q_lcl, + q_icl, + q_rai, + q_sno, + ρ, + T, +) + return conv_q_vap_to_q_lcl_MM2015( + tps, + q_tot, + q_lcl, + q_icl, + q_rai, + q_sno, + ρ, + T, + liquid.τ_relax, + ) +end + +function conv_q_vap_to_q_icl_MM2015( + tps::TDI.PS, + q_tot, + q_lcl, + q_icl, + q_rai, + q_sno, + ρ, + T, + τ_relax, ) FT = eltype(tps) Rᵥ = TDI.Rᵥ(tps) @@ -191,9 +260,30 @@ function conv_q_vap_to_q_lcl_icl_MM2015( return ifelse(limiter, zero(tendency), tendency) end -### -------------------------- ### -### 1-moment terminal velocity ### -### -------------------------- ### +# just to be able to still use old functionality -- to be deleted +function conv_q_vap_to_q_lcl_icl_MM2015( + ice::CMP.CloudIce, + tps::TDI.PS, + q_tot, + q_lcl, + q_icl, + q_rai, + q_sno, + ρ, + T, +) + return conv_q_vap_to_q_icl_MM2015( + tps, + q_tot, + q_lcl, + q_icl, + q_rai, + q_sno, + ρ, + T, + ice.τ_relax, + ) +end """ terminal_velocity(sediment::CloudLiquid, vel::StokesRegimeVelType, ρₐ, q) diff --git a/test/bulk_tendencies_tests.jl b/test/bulk_tendencies_tests.jl index d265a30d2..f0920c909 100644 --- a/test/bulk_tendencies_tests.jl +++ b/test/bulk_tendencies_tests.jl @@ -150,6 +150,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -172,6 +173,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -194,6 +196,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -218,6 +221,7 @@ function test_bulk_microphysics_1m_tendencies(FT) # Call fused function tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -248,6 +252,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -268,6 +273,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -290,6 +296,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -317,6 +324,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -343,6 +351,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -365,6 +374,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -390,6 +400,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -410,6 +421,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -431,6 +443,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies_cold = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ_cold, T_cold, q_tot_cold, FT(0), q_icl_cold, FT(0), q_sno_cold, ) @@ -446,6 +459,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies_warm = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ_warm, T_warm, q_tot_warm, q_lcl_warm, FT(0), q_rai_warm, FT(0), ) @@ -473,6 +487,7 @@ function test_bulk_microphysics_1m_tendencies(FT) # @inferred will throw if return type cannot be inferred tendencies = @inferred BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -493,6 +508,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -522,6 +538,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -548,6 +565,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, FT(0), FT(0), q_sno, ) @@ -589,6 +607,7 @@ function test_bulk_microphysics_1m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -843,6 +862,7 @@ function test_average_bulk_microphysics_1m_tendencies(FT) inst = BMT.bulk_microphysics_tendencies( BMT.Microphysics1Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, ) @@ -1053,6 +1073,7 @@ function test_bulk_microphysics_2m_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics2Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, @@ -1082,6 +1103,7 @@ function test_bulk_microphysics_2m_tendencies(FT) tend_cond = BMT.bulk_microphysics_tendencies( BMT.Microphysics2Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot_super, q_lcl, n_lcl, q_rai, n_rai, ) @test tend_cond.dq_lcl_dt > FT(0) # Cloud increases @@ -1092,6 +1114,7 @@ function test_bulk_microphysics_2m_tendencies(FT) tend_evap = BMT.bulk_microphysics_tendencies( BMT.Microphysics2Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, T, q_tot_sub, q_lcl, n_lcl, q_rai, n_rai, ) @test tend_evap.dq_lcl_dt < FT(0) # Cloud decreases @@ -1105,6 +1128,7 @@ function test_bulk_microphysics_2m_tendencies(FT) # so we only test that the return is a NamedTuple with correct keys tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics2Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, @@ -1163,6 +1187,7 @@ function test_bulk_microphysics_p3_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics2Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, @@ -1209,6 +1234,7 @@ function test_bulk_microphysics_p3_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics2Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, @@ -1255,6 +1281,7 @@ function test_bulk_microphysics_p3_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics2Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, @@ -1299,6 +1326,7 @@ function test_bulk_microphysics_p3_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics2Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, @@ -1346,6 +1374,7 @@ function test_bulk_microphysics_p3_tendencies(FT) tendencies = BMT.bulk_microphysics_tendencies( BMT.Microphysics2Moment(), + BMT.NonEq_Constant(), mp, tps, ρ, diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index 4d37ab7e9..3712433af 100644 --- a/test/gpu_tests.jl +++ b/test/gpu_tests.jl @@ -80,8 +80,8 @@ end i = @index(Global, Linear) FT = eltype(tps) q_lcl = q_icl = q_rai = q_sno = FT(0) # set to zero in this test - S_cond_MM2015 = CMN.conv_q_vap_to_q_lcl_icl_MM2015( - lcl, tps, qᵥ_sl[i], q_lcl, q_icl, q_rai, q_sno, ρ[i], T[i], + S_cond_MM2015 = CMN.conv_q_vap_to_q_lcl_MM2015( + tps, qᵥ_sl[i], q_lcl, q_icl, q_rai, q_sno, ρ[i], T[i], lcl.τ_relax, ) S_cond = CMN.conv_q_vap_to_q_lcl_icl(icl, qᵢ_s[i], qᵢ[i]) output[i] = (; S_cond_MM2015, S_cond) @@ -352,8 +352,9 @@ end ) i = @index(Global, Linear) CM1M = BMT.Microphysics1Moment() + noneq_scheme = BMT.NonEq_Constant() output[i] = BMT.bulk_microphysics_tendencies( - CM1M, mp, tps, ρ[i], T[i], q_tot[i], q_lcl[i], q_icl[i], q_rai[i], q_sno[i], + CM1M, noneq_scheme, mp, tps, ρ[i], T[i], q_tot[i], q_lcl[i], q_icl[i], q_rai[i], q_sno[i], ) end @@ -363,8 +364,9 @@ end ) i = @index(Global, Linear) CM2M = BMT.Microphysics2Moment() + noneq_scheme = BMT.NonEq_Constant() output[i] = BMT.bulk_microphysics_tendencies( - CM2M, mp, tps, ρ[i], T[i], q_tot[i], q_lcl[i], n_lcl[i], q_rai[i], n_rai[i], + CM2M, noneq_scheme, mp, tps, ρ[i], T[i], q_tot[i], q_lcl[i], n_lcl[i], q_rai[i], n_rai[i], ) end @@ -373,8 +375,9 @@ end ) i = @index(Global, Linear) CM2M = BMT.Microphysics2Moment() + noneq_scheme = BMT.NonEq_Constant() output[i] = BMT.bulk_microphysics_tendencies( - CM2M, mp, tps, ρ[i], T[i], q_tot[i], q_lcl[i], n_lcl[i], q_rai[i], n_rai[i], + CM2M, noneq_scheme, mp, tps, ρ[i], T[i], q_tot[i], q_lcl[i], n_lcl[i], q_rai[i], n_rai[i], q_ice[i], n_ice[i], q_rim[i], b_rim[i], ) end diff --git a/test/microphysics_noneq_tests.jl b/test/microphysics_noneq_tests.jl index 789259c29..56ad43cbc 100644 --- a/test/microphysics_noneq_tests.jl +++ b/test/microphysics_noneq_tests.jl @@ -57,25 +57,24 @@ function test_microphysics_noneq(FT) #! format: off # test sign - TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(liquid, tps, FT(0.5 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T) == FT(0) - TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(liquid, tps, FT(1.5 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T) > FT(0) - TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(liquid, tps, qᵥ_sl, FT(0), FT(0), FT(0), FT(0), ρ, T) ≈ FT(0) + TT.@test CMNe.conv_q_vap_to_q_lcl_MM2015(tps, FT(0.5 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T, liquid.τ_relax) == FT(0) + TT.@test CMNe.conv_q_vap_to_q_lcl_MM2015(tps, FT(1.5 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T, liquid.τ_relax) > FT(0) + TT.@test CMNe.conv_q_vap_to_q_lcl_MM2015(tps, qᵥ_sl, FT(0), FT(0), FT(0), FT(0), ρ, T, liquid.τ_relax) ≈ FT(0) - TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(ice, tps, FT(0.5 * qᵥ_si), FT(0), FT(0), FT(0), FT(0), ρ, T) == FT(0) - TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(ice, tps, FT(1.5 * qᵥ_si), FT(0), FT(0), FT(0), FT(0), ρ, T) > FT(0) - TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(ice, tps, qᵥ_si, FT(0), FT(0), FT(0), FT(0), ρ, T) ≈ FT(0) + TT.@test CMNe.conv_q_vap_to_q_icl_MM2015(tps, FT(0.5 * qᵥ_si), FT(0), FT(0), FT(0), FT(0), ρ, T, ice.τ_relax) == FT(0) + TT.@test CMNe.conv_q_vap_to_q_icl_MM2015(tps, FT(1.5 * qᵥ_si), FT(0), FT(0), FT(0), FT(0), ρ, T, ice.τ_relax) > FT(0) + TT.@test CMNe.conv_q_vap_to_q_icl_MM2015(tps, qᵥ_si, FT(0), FT(0), FT(0), FT(0), ρ, T, ice.τ_relax) ≈ FT(0) # smoke test for values - TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(liquid, tps, FT(1.2 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T) ≈ 3.763045798130144e-5 rtol = 1e-6 - TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(ice, tps, FT(1.2 * qᵥ_si), FT(0), FT(0), FT(0), FT(0), ρ, T) ≈ 3.235984203087906e-5 rtol = 1e-6 + TT.@test CMNe.conv_q_vap_to_q_lcl_MM2015(tps, FT(1.2 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T, liquid.τ_relax) ≈ 3.763045798130144e-5 rtol = 1e-6 + TT.@test CMNe.conv_q_vap_to_q_icl_MM2015(tps, FT(1.2 * qᵥ_si), FT(0), FT(0), FT(0), FT(0), ρ, T, ice.τ_relax) ≈ 3.235984203087906e-5 rtol = 1e-6 - TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(liquid, tps, FT(1.2 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T) ≈ 3.7630474f-5 rtol = 1e-6 - TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(ice, tps, FT(1.2 * qᵥ_si), FT(0), FT(0), FT(0), FT(0), ρ, T) ≈ 3.2359854f-5 rtol = 1e-6 + TT.@test CMNe.conv_q_vap_to_q_lcl_MM2015(tps, FT(1.2 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T, liquid.τ_relax) ≈ 3.7630474f-5 rtol = 1e-6 + TT.@test CMNe.conv_q_vap_to_q_icl_MM2015(tps, FT(1.2 * qᵥ_si), FT(0), FT(0), FT(0), FT(0), ρ, T, ice.τ_relax) ≈ 3.2359854f-5 rtol = 1e-6 # ice grows faster than liquid - TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(liquid, tps, FT(1.2 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T) < - CMNe.conv_q_vap_to_q_lcl_icl_MM2015(ice, tps, FT(1.2 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T) - + TT.@test CMNe.conv_q_vap_to_q_lcl_MM2015(tps, FT(1.2 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T, liquid.τ_relax) < + CMNe.conv_q_vap_to_q_icl_MM2015(tps, FT(1.2 * qᵥ_sl), FT(0), FT(0), FT(0), FT(0), ρ, T, ice.τ_relax) # --- INP limiter: above freezing, positive ice deposition should be zero --- T_warm = FT(280) # above freezing diff --git a/test/type_stability_tests.jl b/test/type_stability_tests.jl index 7261cd67a..02846ef74 100644 --- a/test/type_stability_tests.jl +++ b/test/type_stability_tests.jl @@ -71,6 +71,7 @@ function run_type_stability_tests() tendencies_1M = BMT.bulk_microphysics_tendencies.( Ref(BMT.Microphysics1Moment()), + Ref(BMT.NonEq_Constant()), Ref(mp1), Ref(tps), ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, @@ -91,6 +92,7 @@ function run_type_stability_tests() tendencies_2M_warm = BMT.bulk_microphysics_tendencies.( Ref(BMT.Microphysics2Moment()), + Ref(BMT.NonEq_Constant()), Ref(mp2_warm), Ref(tps), ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai, @@ -113,6 +115,7 @@ function run_type_stability_tests() tendencies_2M_p3 = BMT.bulk_microphysics_tendencies.( Ref(BMT.Microphysics2Moment()), + Ref(BMT.NonEq_Constant()), Ref(mp2_p3), Ref(tps), ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai,