From acad9b42d1c3d488def53b9bc81aee0916123cea Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Mon, 23 Jun 2025 15:01:10 -0700 Subject: [PATCH 1/3] implement ARG aerosol activation with prescribed aerosol densities for 2M microphysics --- .buildkite/Manifest-v1.11.toml | 6 +- Project.toml | 6 +- .../aquaplanet_nonequil_allsky_gw_res_2M.yml | 2 +- docs/src/equations.md | 11 ++ reproducibility_tests/ref_counter.jl | 5 +- .../microphysics/cloud_condensate.jl | 85 +++++++-- .../microphysics/microphysics_wrappers.jl | 178 +++++++++++------- src/parameters/create_parameters.jl | 22 +++ 8 files changed, 223 insertions(+), 92 deletions(-) diff --git a/.buildkite/Manifest-v1.11.toml b/.buildkite/Manifest-v1.11.toml index 891d84bf96..e0759210e8 100644 --- a/.buildkite/Manifest-v1.11.toml +++ b/.buildkite/Manifest-v1.11.toml @@ -461,10 +461,10 @@ uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9" version = "0.1.13" [[deps.CloudMicrophysics]] -deps = ["ClimaParams", "DocStringExtensions", "ForwardDiff", "LazyArtifacts", "LogExpFunctions", "QuadGK", "RootSolvers", "SpecialFunctions", "Thermodynamics"] -git-tree-sha1 = "cfc9decc96d616ab3b8d911f843eeb2aafd7d2ad" +deps = ["ClimaParams", "DocStringExtensions", "ForwardDiff", "LazyArtifacts", "LogExpFunctions", "QuadGK", "RootSolvers", "SpecialFunctions", "StaticArrays", "Thermodynamics"] +git-tree-sha1 = "8f7137d8bb3098d6c46b2c3f859fb873d5b28abd" uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" -version = "0.25.0" +version = "0.26.1" [deps.CloudMicrophysics.extensions] EmulatorModelsExt = ["DataFrames", "MLJ"] diff --git a/Project.toml b/Project.toml index bdefca8f3f..4611b5f871 100644 --- a/Project.toml +++ b/Project.toml @@ -49,7 +49,7 @@ ClimaInterpolations = "0.1.0" ClimaParams = "0.10.35" ClimaTimeSteppers = "0.8.2" ClimaUtilities = "0.1.22" -CloudMicrophysics = "0.25.0" +CloudMicrophysics = "0.26.1" Dates = "1" ForwardDiff = "0.10.38, 1" Insolation = "0.9.2" @@ -64,10 +64,10 @@ NullBroadcasts = "0.1" RRTMGP = "0.21.3" Random = "1" SciMLBase = "2.12" -StaticArrays = "1.7" +StaticArrays = "1.9" Statistics = "1" SurfaceFluxes = "0.11, 0.12" -Thermodynamics = "0.12.13" +Thermodynamics = "0.12.14" UnrolledUtilities = "0.1.9" YAML = "0.4" julia = "1.9" diff --git a/config/model_configs/aquaplanet_nonequil_allsky_gw_res_2M.yml b/config/model_configs/aquaplanet_nonequil_allsky_gw_res_2M.yml index 2b577ec244..e8f83c1e07 100644 --- a/config/model_configs/aquaplanet_nonequil_allsky_gw_res_2M.yml +++ b/config/model_configs/aquaplanet_nonequil_allsky_gw_res_2M.yml @@ -20,7 +20,7 @@ orographic_gravity_wave: "gfdl_restart" surface_setup: "DefaultMoninObukhov" prescribe_ozone: true reproducibility_test: true -prescribed_aerosols: ["SO4", "CB1", "OC1", "DST01", "SSLT01"] +prescribed_aerosols: ["SO4", "CB1", "OC1", "DST01", "SSLT01", "SSLT02", "SSLT03", "SSLT04", "SSLT05"] diagnostics: - short_name: [edt, evu, mmrso4, mmrbcpo, mmrocpo, mmrdust, mmrss, loadss, o3, od550aer, odsc550aer] reduction_time: average diff --git a/docs/src/equations.md b/docs/src/equations.md index 2358177ce0..4b3cfc447b 100644 --- a/docs/src/equations.md +++ b/docs/src/equations.md @@ -414,3 +414,14 @@ The source terms functions treat negative specific humidities as zeros, so the simulations should be stable even with small negative numbers. We do not apply hyperdiffusion for precipitation tracers. + +### Aerosol Activation for 2-Moment Microphysics + +Aerosol activation uses functions from the [CloudMicrophysics.jl](https://github.com/CliMA/CloudMicrophysics.jl) library, based on the Abdul-Razzak and Ghan (ARG) parameterization. ARG predicts the number of activated cloud droplets assuming a parcel of clear air rising adiabatically. This formulation is traditionally applied only at cloud base, where the maximum supersaturation typically occurs. + +To enable ARG to be used locally (i.e., without explicitly identifying cloud base), CloudMicrophysics.jl implements a modified equation for the maximum supersaturation that accounts for the presence of pre-existing liquid and ice particles. This allows activation to be applied inside clouds. To ensure that activation occurs only where physically appropriate, we apply additional clipping logic: + +- If the predicted maximum supersaturation ``S_max`` is less than the local supersaturation (i.e., supersaturation is decreasing), aerosol activation is not applied. +- If the predicted number of activated droplets is less than the existing local cloud droplet number concentration, activation is also suppressed. + +This ensures that droplet activation occurs only in physically meaningful regions—typically near cloud base—even though the activation routine can be applied throughout the domain. diff --git a/reproducibility_tests/ref_counter.jl b/reproducibility_tests/ref_counter.jl index 6491e2a62e..20b6cba58a 100644 --- a/reproducibility_tests/ref_counter.jl +++ b/reproducibility_tests/ref_counter.jl @@ -1,4 +1,4 @@ -248 +249 # **README** # @@ -20,6 +20,9 @@ #= +249 +- Add aerosol activation for 2M microphysics + 248 - update deps: climacore 0.14.35 diff --git a/src/parameterized_tendencies/microphysics/cloud_condensate.jl b/src/parameterized_tendencies/microphysics/cloud_condensate.jl index fd0fd6e385..3fa4cfde95 100644 --- a/src/parameterized_tendencies/microphysics/cloud_condensate.jl +++ b/src/parameterized_tendencies/microphysics/cloud_condensate.jl @@ -74,13 +74,13 @@ function cloud_condensate_tendency!( (; params, dt) = p FT = eltype(params) thp = CAP.thermodynamics_params(params) - cmc = CAP.microphysics_cloud_params(params) + cmp = CAP.microphysics_2m_params(params) Tₐ = @. lazy(TD.air_temperature(thp, ᶜts)) @. Yₜ.c.ρq_liq += Y.c.ρ * cloud_sources( - cmc.liquid, + cmp.liquid, thp, specific(Y.c.ρq_tot, Y.c.ρ), specific(Y.c.ρq_liq, Y.c.ρ), @@ -93,7 +93,7 @@ function cloud_condensate_tendency!( ) @. Yₜ.c.ρq_ice += Y.c.ρ * cloud_sources( - cmc.ice, + cmp.ice, thp, specific(Y.c.ρq_tot, Y.c.ρ), specific(Y.c.ρq_liq, Y.c.ρ), @@ -105,19 +105,68 @@ function cloud_condensate_tendency!( dt, ) - @. Yₜ.c.ρn_liq += - Y.c.ρ * aerosol_activation_sources( - cmc.liquid, - thp, - Y.c.ρ, - Tₐ, - specific(Y.c.ρq_tot, Y.c.ρ), - specific(Y.c.ρq_liq, Y.c.ρ), - specific(Y.c.ρq_ice, Y.c.ρ), - specific(Y.c.ρq_rai, Y.c.ρ), - specific(Y.c.ρq_sno, Y.c.ρ), - specific(Y.c.ρn_liq, Y.c.ρ), - cmc.N_cloud_liquid_droplets, - dt, - ) + # Aerosol activation using prescribed aerosol (Sea salt and sulfate) + if !(:prescribed_aerosols_field in propertynames(p.tracers)) + return + end + seasalt_num = p.scratch.ᶜtemp_scalar + seasalt_mean_radius = p.scratch.ᶜtemp_scalar_2 + sulfate_num = p.scratch.ᶜtemp_scalar_3 + @. seasalt_num = 0 + @. seasalt_mean_radius = 0 + @. sulfate_num = 0 + # Get aerosol concentrations if available + seasalt_names = [:SSLT01, :SSLT02, :SSLT03, :SSLT04, :SSLT05] + sulfate_names = [:SO4] + aerosol_field = p.tracers.prescribed_aerosols_field + for aerosol_name in propertynames(aerosol_field) + if aerosol_name in seasalt_names + seasalt_particle_radius = getproperty( + cmp.aerosol, + Symbol(string(aerosol_name) * "_radius"), + ) + seasalt_particle_mass = + FT(4 / 3 * pi) * + seasalt_particle_radius^3 * + cmp.aerosol.seasalt_density + seasalt_mass = getproperty(aerosol_field, aerosol_name) + @. seasalt_num += seasalt_mass / seasalt_particle_mass + @. seasalt_mean_radius += + seasalt_mass / seasalt_particle_mass * + log(seasalt_particle_radius) + elseif aerosol_name in sulfate_names + sulfate_particle_mass = + FT(4 / 3 * pi) * + cmp.aerosol.sulfate_radius^3 * + cmp.aerosol.sulfate_density + sulfate_mass = getproperty(aerosol_field, aerosol_name) + @. sulfate_num += sulfate_mass / sulfate_particle_mass + end + end + # Compute geometric mean radius of the log-normal distribution: + # exp(weighted average of log(radius)) + @. seasalt_mean_radius = + ifelse(seasalt_num == 0, 0, exp(seasalt_mean_radius / seasalt_num)) + + # Compute aerosol activation (ARG 2000) + (; ᶜu) = p.precomputed + ᶜw = p.scratch.ᶜtemp_scalar_4 + Snₗ = p.scratch.ᶜtemp_scalar_5 + @. ᶜw = max(0, w_component.(Geometry.WVector.(ᶜu))) + @. Snₗ = aerosol_activation_sources( + seasalt_num, + seasalt_mean_radius, + sulfate_num, + specific(Y.c.ρq_tot, Y.c.ρ), + specific(Y.c.ρq_liq + Y.c.ρq_rai, Y.c.ρ), + specific(Y.c.ρq_ice + Y.c.ρq_sno, Y.c.ρ), + specific(Y.c.ρn_liq + Y.c.ρn_rai, Y.c.ρ), + Y.c.ρ, + ᶜw, + (cmp,), + thp, + ᶜts, + dt, + ) + @. Yₜ.c.ρn_liq += Y.c.ρ * Snₗ end diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl index 84a6022f16..1998f29cb9 100644 --- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl +++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl @@ -1,14 +1,18 @@ # A set of wrappers for using CloudMicrophysics.jl functions inside EDMFX loops import Thermodynamics as TD +import CloudMicrophysics.ThermodynamicsInterface as CMTDI import CloudMicrophysics.Microphysics0M as CM0 import CloudMicrophysics.Microphysics1M as CM1 import CloudMicrophysics.Microphysics2M as CM2 import CloudMicrophysics.MicrophysicsNonEq as CMNe +import CloudMicrophysics.AerosolModel as CMAM +import CloudMicrophysics.AerosolActivation as CMAA import CloudMicrophysics.Parameters as CMP # Define some aliases and functions to make the code more readable const Tₐ = TD.air_temperature +const Pₐ = TD.air_pressure const PP = TD.PhasePartition const qᵥ = TD.vapor_specific_humidity @@ -25,33 +29,6 @@ function limit(q, dt, n::Int) return max(FT(0), q) / float(dt) / n end -""" - ml_N_cloud_liquid_droplets(cmc, c_dust, c_seasalt, c_SO4, q_liq) - - - cmc - a struct with cloud and aerosol parameters - - c_dust, c_seasalt, c_SO4 - dust, seasalt and ammonium sulfate mass concentrations [kg/kg] - - q_liq - liquid water specific humidity - -Returns the liquid cloud droplet number concentration diagnosed based on the -aerosol loading and cloud liquid water. -""" -function ml_N_cloud_liquid_droplets(cmc, c_dust, c_seasalt, c_SO4, q_liq) - # We can also add w, T, RH, w' ... - # Also consider lookind only at around cloud base height - (; α_dust, α_seasalt, α_SO4, α_q_liq) = cmc.aml - (; c₀_dust, c₀_seasalt, c₀_SO4, q₀_liq) = cmc.aml - N₀ = cmc.N_cloud_liquid_droplets - - FT = eltype(N₀) - return N₀ * ( - FT(1) + - α_dust * (log(max(c_dust, eps(FT))) - log(c₀_dust)) + - α_seasalt * (log(max(c_seasalt, eps(FT))) - log(c₀_seasalt)) + - α_SO4 * (log(max(c_SO4, eps(FT))) - log(c₀_SO4)) + - α_q_liq * (log(max(q_liq, eps(FT))) - log(q₀_liq)) - ) -end - """ cloud_sources(cm_params, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Tₐ, dt) @@ -386,60 +363,129 @@ end ##### """ - aerosol_activation_sources(cm_params, thp, ρ, Tₐ, qₜ, qₗ, qᵢ, qᵣ, qₛ, n_dp, n_dp_prescribed, dt) - - - cm_params - CloudMicrophysics parameters struct for cloud water or ice condensate - - thp - Thermodynamics parameters struct - - ρ - air density - - Tₐ - air temperature - - qₜ - total specific humidity - - qₗ - liquid specific humidity - - qᵢ - ice specific humidity - - qᵣ - rain specific humidity - - qₛ - snow specific humidity - - n_dp - number concentration droplets (liquid or ice) per mass - _ n_dp_prescribed - prescribed number concentration of droplets (liquid or ice) per mass - - dt - model time step + aerosol_activation_sources( + seasalt_num, + seasalt_mean_radius, + sulfate_num, + qₜ, + qₗ, + qᵢ, + nₗ, + ρ, + w, + cmp, + thermo_params, + ts, + dt, + ) -Returns the activation rate. #TODO This function temporarily computes activation rate -based on mass rates and a prescribed droplet mass (no activation parameterization yet). +Computes the source term for cloud droplet number concentration per mass due to aerosol activation, +based on the Abdul-Razzak and Ghan (2000) parameterization. + +This function estimates the number of aerosols activated into cloud droplets per mass of air per second +from a bi-modal aerosol distribution (sea salt and sulfate), given local supersaturation and vertical +velocity. The result is returned as a tendency (per second) of liquid droplet number concentration. + +# Arguments +- `seasalt_num`: Number concentration per mass of sea salt aerosols [kg⁻¹]. +- `seasalt_mean_radius`: Mean dry radius of sea salt aerosol mode [m]. +- `sulfate_num`: Number concentration per mass of sulfate aerosols [kg⁻¹]. +- `qₜ`, `qₗ`, `qᵢ` - total water, liquid (cloud liquid and rain) and ice (cloud ice and snow) specific humidity +- `nₗ` - lqiuid (cloud liquid and rain) number concentration per mass [kg⁻¹] +- `ρ`: Air density [kg/m³]. +- `w`: Vertical velocity [m/s]. +- `cmp`: Microphysics parameters +- `thermo_params`: Thermodynamic parameters for computing saturation, pressure, temperature, etc. +- `ts`: Thermodynamic state (e.g., prognostic variables) used for evaluating the phase partition. +- `dt`: Time step (s) over which the activation tendency is applied. + +# Returns +- Tendency of cloud liquid droplet number concentration per mass of air due to aerosol activation [m⁻³/s]. """ function aerosol_activation_sources( - cm_params::CMP.CloudLiquid{FT}, - thp, - ρ, - Tₐ, + seasalt_num, + seasalt_mean_radius, + sulfate_num, qₜ, qₗ, qᵢ, - qᵣ, - qₛ, - n_dp, - n_dp_prescribed, + nₗ, + ρ, + w, + cmp, + thermo_params, + ts, dt, -) where {FT} - r_dp = FT(2e-6) # 2 μm - m_dp = 4 / 3 * FT(π) * r_dp^3 * cm_params.ρw - Sn = cloud_sources(cm_params, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Tₐ, dt) / m_dp +) + + FT = eltype(nₗ) + air_params = cmp.aps + arg_params = cmp.arg + aerosol_params = cmp.aerosol + T = Tₐ(thermo_params, ts) + p = Pₐ(thermo_params, ts) + S = CMTDI.supersaturation_over_liquid(thermo_params, qₜ, qₗ, qᵢ, ρ, T) + n_aer = seasalt_num + sulfate_num + if S < 0 || n_aer < eps(FT) + return FT(0) + end + + seasalt_mode = CMAM.Mode_κ( + seasalt_mean_radius, + aerosol_params.seasalt_std, + seasalt_num * ρ, + (FT(1),), + (FT(1),), + (FT(0),), + (aerosol_params.seasalt_kappa,), + ) + sulfate_mode = CMAM.Mode_κ( + aerosol_params.sulfate_radius, + aerosol_params.sulfate_std, + sulfate_num * ρ, + (FT(1),), + (FT(1),), + (FT(0),), + (aerosol_params.sulfate_kappa,), + ) + + aerosol_dist = CMAM.AerosolDistribution((seasalt_mode, sulfate_mode)) + + args = ( + arg_params, + aerosol_dist, + air_params, + thermo_params, + T, + p, + w, + qₜ, + qₗ, + qᵢ, + ρ * nₗ, + FT(0), + ) #assuming no ice particles because we don't track n_ice for now + S_max = CMAA.max_supersaturation(args...) + n_act = CMAA.total_N_activated(args...) / ρ return ifelse( - Sn > FT(0), - triangle_inequality_limiter(Sn, limit((n_dp_prescribed - n_dp), dt, 2)), - -triangle_inequality_limiter(abs(Sn), limit(n_dp, dt, 2)), + S_max < S || isnan(n_act) || n_act < nₗ, + FT(0), + (n_act - nₗ) / dt, ) end """ compute_warm_precipitation_sources_2M!(Sᵖ, S₂ᵖ, Snₗᵖ, Snᵣᵖ, Sqₗᵖ, Sqᵣᵖ, ρ, nₗ, nᵣ, qₜ, qₗ, qᵢ, qᵣ, qₛ, ts, dt, sb, thp) - - Sᵖ, S₂ᵖ - temporary containters to help compute precipitation source terms - - Snₗᵖ, Snᵣᵖ, Sqₗᵖ, Sqᵣᵖ - cached storage for precipitation source terms - - ρ - air density - - nₗ, nᵣ - cloud liquid and rain number concentration per mass [1 / kg of moist air] - - qₜ, qₗ, qᵢ, qᵣ, qₛ - total water, cloud liquid, cloud ice, rain and snow specific humidity - - ts - thermodynamic state (see td package for details) - - dt - model time step - - thp, mp - structs with thermodynamic and microphysics parameters + - `Sᵖ`, `S₂ᵖ` - temporary containters to help compute precipitation source terms + - `Snₗᵖ`, `Snᵣᵖ`, `Sqₗᵖ`, `Sqᵣᵖ` - cached storage for precipitation source terms + - `ρ` - air density + - `nₗ`, `nᵣ` - cloud liquid and rain number concentration per mass [1 / kg of moist air] + - `qₜ`, `qₗ`, `qᵢ`, `qᵣ`, `qₛ` - total water, cloud liquid, cloud ice, rain and snow specific humidity + - `ts` - thermodynamic state (see td package for details) + - `dt` - model time step + - `thp`, `mp` - structs with thermodynamic and microphysics parameters Computes precipitation number and mass sources due to warm precipitation processes based on the 2-moment [Seifert and Beheng (2006) scheme](https://clima.github.io/CloudMicrophysics.jl/dev/Microphysics2M/). diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl index 117bef9636..b855da66a5 100644 --- a/src/parameters/create_parameters.jl +++ b/src/parameters/create_parameters.jl @@ -163,6 +163,10 @@ microphys_2m_parameters(toml_dict::CP.AbstractTOMLDict) = (; sb = CM.Parameters.SB2006(toml_dict), aps = CM.Parameters.AirProperties(toml_dict), tv = CM.Parameters.SB2006VelType(toml_dict), + liquid = CM.Parameters.CloudLiquid(toml_dict), + ice = CM.Parameters.CloudIce(toml_dict), + arg = CM.Parameters.AerosolActivationParameters(toml_dict), + aerosol = prescribed_aerosol_parameters(toml_dict), ) function vert_diff_parameters(toml_dict) @@ -195,6 +199,24 @@ function aerosol_ml_parameters(toml_dict) return CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos") end +function prescribed_aerosol_parameters(toml_dict) + name_map = (; + :MERRA2_seasalt_aerosol_bin01_radius => :SSLT01_radius, + :MERRA2_seasalt_aerosol_bin02_radius => :SSLT02_radius, + :MERRA2_seasalt_aerosol_bin03_radius => :SSLT03_radius, + :MERRA2_seasalt_aerosol_bin04_radius => :SSLT04_radius, + :MERRA2_seasalt_aerosol_bin05_radius => :SSLT05_radius, + :seasalt_aerosol_kappa => :seasalt_kappa, + :seasalt_aerosol_density => :seasalt_density, + :mam3_stdev_coarse => :seasalt_std, + :MERRA2_sulfate_aerosol_radius => :sulfate_radius, + :sulfate_aerosol_kappa => :sulfate_kappa, + :sulfate_aerosol_density => :sulfate_density, + :mam3_stdev_accum => :sulfate_std, + ) + return CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos") +end + to_svec(x::AbstractArray) = SA.SVector{length(x)}(x) to_svec(x) = x to_svec(x::NamedTuple) = map(x -> to_svec(x), x) From f450a5a4963a9272f2560507a35e785b88f14996 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Wed, 23 Jul 2025 16:44:41 -0700 Subject: [PATCH 2/3] minor --- docs/src/equations.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/equations.md b/docs/src/equations.md index 4b3cfc447b..b9f9c3ffad 100644 --- a/docs/src/equations.md +++ b/docs/src/equations.md @@ -421,7 +421,7 @@ Aerosol activation uses functions from the [CloudMicrophysics.jl](https://github To enable ARG to be used locally (i.e., without explicitly identifying cloud base), CloudMicrophysics.jl implements a modified equation for the maximum supersaturation that accounts for the presence of pre-existing liquid and ice particles. This allows activation to be applied inside clouds. To ensure that activation occurs only where physically appropriate, we apply additional clipping logic: -- If the predicted maximum supersaturation ``S_max`` is less than the local supersaturation (i.e., supersaturation is decreasing), aerosol activation is not applied. +- If the predicted maximum supersaturation is less than the local supersaturation (i.e., supersaturation is decreasing), aerosol activation is not applied. - If the predicted number of activated droplets is less than the existing local cloud droplet number concentration, activation is also suppressed. This ensures that droplet activation occurs only in physically meaningful regions—typically near cloud base—even though the activation routine can be applied throughout the domain. From 8aef1db977015eea72b0b7109da0c02800dfb665 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 24 Jul 2025 14:35:44 -0700 Subject: [PATCH 3/3] bugfix for gpu compatibility --- .../microphysics/microphysics_wrappers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl index 1998f29cb9..f457d1932b 100644 --- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl +++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl @@ -471,7 +471,7 @@ function aerosol_activation_sources( return ifelse( S_max < S || isnan(n_act) || n_act < nₗ, FT(0), - (n_act - nₗ) / dt, + (n_act - nₗ) / float(dt), ) end