Subfilter available potential energy equation

The SubFilterAvailablePotentialEnergyEquation module provides diagnostics for the available potential energy carried by the scales that a low-pass spatial filter $\overline{(\,\cdot\,)}$ removes from the buoyancy field. It is the potential-energy counterpart of the Subfilter kinetic energy equation: where that module splits the kinetic energy across the filter scale, this one splits the local available potential energy $e_a$ of the available potential energy equation. The other half of the split, the APE of the filtered buoyancy and its dissipation, lives in the Filtered available potential energy equation and is re-exported here.

Both the full and the filtered buoyancy are measured against one shared reference profile $(b^\star, z^\star)$, ordinarily the sorted state of the full buoyancy field. The subfilter available potential energy is the filtered full APE minus the APE of the filtered buoyancy $\bar b$,

\[e_a^s = \overline{e_a(b, z)} - e_a(\bar b, z) , \qquad e_a(b, z) = \int_{z^\star(b)}^{z} \left[b^\star(\tilde z) - b\right] \mathrm{d}\tilde z ,\]

computed by SubFilterAvailablePotentialEnergy. The second term is the APE of the filtered buoyancy, $e_a^l = e_a(\bar b, z)$, computed by FilteredAvailablePotentialEnergy. Since the filter acts at fixed $z$, $e_a^s$ is also $\overline{e_a(b, z^\star(\bar b))}$: the average over the filter of the APE each parcel holds between its own reference height $z^\star(b)$ and the one the filtered buoyancy would assign it, $z^\star(\bar b)$. Because $e_a$ is convex in buoyancy, a filter that averages at fixed $z$ keeps $e_a^s \geq 0$ pointwise, by Jensen's inequality.

Looking the filtered buoyancy up in a profile it did not itself produce is what ProfileLookup was built for, so these diagnostics accept only that reference-height method: the default sorts the model's own buoyancy into a VerticalSort column on every compute!, a column you built yourself can be shared across diagnostics, and a profile given as plain arrays holds the reference state fixed in time (which also makes the diagnostics sort-free).

Horizontal filters only

The filter has to act only in the horizontal directions, and every diagnostic on this page throws an ArgumentError for a filter with a vertical component. Such a filter mixes heights as well as buoyancies and needs extra steps, not currently implemented, to guarantee a non-negative APE.

Deriving the subfilter available potential energy equation

Since $e_a^s = \overline{e_a} - e_a^l$, its budget is the budget of $\overline{e_a}$ minus that of $e_a^l$. Filtering the full-field APE budget (the filter commutes with derivatives) and splitting the filtered advective flux into a resolved and a subfilter part, $\overline{u_i e_a} = \bar u_i \, \overline{e_a} + \tau(u_i, e_a)$, gives the budget of $\overline{e_a}$,

\[\partial_t \overline{e_a} = -\partial_i(\bar u_i \, \overline{e_a}) - \overline{w b_r} - \partial_i\left[\tau(u_i, e_a) + \overline{\Upsilon q_i}\right] - \overline{\varepsilon_a} + \overline{R} ,\]

where $\tau(f, g) = \overline{f g} - \bar f \, \bar g$ is the subfilter covariance of two fields, and the buoyancy anomaly $b_r = b - b^\star(z)$, the displacement potential $\Upsilon = z^\star(b) - z$, the diffusive buoyancy flux $q_i$, the dissipation $\varepsilon_a$ and the reference tendency $R$ are those of the full-field budget. The filtered APE budget is

\[\partial_t e_a^l = -\partial_i(\bar u_i e_a^l) - \bar w \, \overline{b_r} - \Pi_a - \partial_i\left[\Upsilon^l \left(\tau(u_i, b) + \bar q_i\right)\right] - \varepsilon_a^l + R^l ,\]

and subtracting it from the budget of $\overline{e_a}$ term by term gives the subfilter APE budget

\[\partial_t e_a^s = \underbrace{-\partial_i(\bar u_i e_a^s)}_{\text{advection}} \underbrace{-\,\tau(w, b_r)}_{\text{APE to KE conversion}} + \underbrace{\Pi_a}_{\text{cross-scale flux}} \underbrace{-\,\partial_i F^s_i}_{\text{subfilter transport}} \underbrace{-\,\varepsilon_a^s}_{\text{dissipation}} + \underbrace{R^s}_{\text{reference tendency}} ,\]

with

\[\tau(w, b_r) = \overline{w b_r} - \bar w \, \overline{b_r} , \qquad F^s_i = \tau(u_i, e_a) + \overline{\Upsilon q_i} - \Upsilon^l \left(\tau(u_i, b) + \bar q_i\right) , \qquad \varepsilon_a^s = \overline{\varepsilon_a} - \varepsilon_a^l , \qquad R^s = \overline{R} - R^l .\]

The cross-scale flux $\Pi_a = -\tau(u_i, b) \, \partial_i \Upsilon^l$ appears with the opposite sign in the filtered budget, so it moves APE between the filtered and the subfilter scales, from the former to the latter where $\Pi_a > 0$. The conversion $\tau(w, b_r)$ likewise appears with the opposite sign in the Subfilter kinetic energy equation, which makes it a reversible exchange between the subfilter APE and the subfilter kinetic energy. With the filtered conversion $\bar w \, \overline{b_r}$ (FilteredAvailablePotentialToKineticEnergyConversion) it adds up to $\overline{w b_r}$. For the horizontal filters accepted here $b^\star$ depends on $z$ alone, so $\overline{b_r} = \bar b - b^\star(z)$, the reference profile cancels between the two terms of $\tau(w, b_r)$, and the conversion is the subfilter buoyancy flux, $\tau(w, b_r) = \tau(w, b) = \overline{w b} - \bar w \, \bar b$.

The dissipation $\varepsilon_a^s$ is the filtered full-field dissipation $\varepsilon_a = -q_i \, \partial_i \Upsilon$ (AvailablePotentialEnergyDissipationRate) minus the dissipation of the filtered state, $\varepsilon_a^l = -\bar q_i \, \partial_i \Upsilon^l$ (FilteredAvailablePotentialEnergyDissipationRate), where $\bar q_i$ is the closure's diffusive buoyancy flux low-pass filtered and $\Upsilon^l = z^\star(\bar b) - z$ is the displacement potential of the filtered buoyancy. Filtering the flux, rather than recomputing it from $\bar b$, is the same choice the Filtered kinetic energy equation makes for the viscous flux; the Filtered available potential energy equation has the details.

The flux $F^s_i$ gathers the subfilter advection of $e_a$, the filtered diffusive transport of the full field and, with the opposite sign, the transport of the filtered budget. Its divergence redistributes $e_a^s$ in space and carries no APE across the filter scale. $R^s$ is the subfilter part of the reference tendency, computed by SubFilterAvailablePotentialEnergyReferenceTendency. It satisfies $\int R^s \, \mathrm{d}V = -\int R^l \, \mathrm{d}V$, as the Filtered available potential energy equation explains, and vanishes identically when the reference profile is held fixed in time.

Terms and diagnostics

Five of the quantities above have diagnostics: four are defined in this module, and the cross-scale flux is defined in the Filtered available potential energy equation and re-exported here. The advection and the subfilter transport have none. The Rayleigh-Taylor instability example closes the volume-integrated budget.

Quantity	Expression	Diagnostic
Subfilter available potential energy	$e_a^s = \overline{e_a(b, z)} - e_a(\bar b, z)$	`SubFilterAvailablePotentialEnergy`
Advection	$\partial_i(\bar u_i e_a^s)$	not implemented
APE to KE conversion	$\tau(w, b_r) = \overline{w b_r} - \bar w \, \overline{b_r}$	`SubFilterAvailablePotentialToKineticEnergyConversion`
Cross-scale flux	$\Pi_a = -\tau(u_i, b) \, \partial_i \Upsilon^l$	`AvailablePotentialEnergyCrossScaleFlux`
Subfilter transport	$\partial_i F^s_i$	not implemented
Dissipation	$\varepsilon_a^s = \overline{\varepsilon_a} - \varepsilon_a^l$	`SubFilterAvailablePotentialEnergyDissipationRate`
Reference tendency	$R^s = \overline{R} - R^l$	`SubFilterAvailablePotentialEnergyReferenceTendency`

$\varepsilon_a^s$ also answers to DissipationRate. Like its full-field and filtered namesakes, the alias is scoped to its module: using Oceanostics.SubFilterAvailablePotentialEnergyEquation brings it in, using Oceanostics does not.

Summary of $e_a^s$ equation terms

Oceanostics.SubFilterAvailablePotentialEnergyEquation.SubFilterAvailablePotentialEnergy — Type

SubFilterAvailablePotentialEnergy(
    model,
    filter;
    method,
    geopotential_height
)

Return the subfilter-scale (SFS) available potential energy eₐˢ, the available potential energy carried by the scales that a low-pass filter removes from the buoyancy field — the filtered full APE minus the APE of the filtered buoyancy b̄ = filter(b):

    eₐˢ = filter(eₐ(b, z)) - eₐ(b̄, z) ,   eₐ(b, z) = ∫_{z✶(b)}^{z} [b✶(z̃) - b] dz̃

where eₐ is the local available potential energy density (AvailablePotentialEnergy), eₐ(b̄, z) is the APE of the filtered buoyancy (FilteredAvailablePotentialEnergy), and both terms are measured against one shared reference profile (b✶, z✶): the filtered buoyancy is looked up in the same profile the full field is measured against, which is exactly what ProfileLookup was built for. It is the potential-energy counterpart of the subfilter kinetic energy eₖˢ (SubFilterKineticEnergy), just as FilteredAvailablePotentialEnergy is that of the filtered kinetic energy eₖˡ.

Because the two states have to share one profile, method must be a ProfileLookup:

ProfileLookup() (the default) sorts the profile from the model's own buoyancy — a VerticalSort column built internally and re-sorted on every compute!, so the reference state tracks the flow.
ProfileLookup(z✶_column) borrows a column you already built with reference_height(model, method=VerticalSort()), sharing its sort across several diagnostics.
ProfileLookup(b✶, z✶) with plain arrays holds the reference profile fixed in time, which also makes the diagnostic sort-free: each compute! is then a filter plus two binary-search lookups.

eₐ is convex in buoyancy (∂²eₐ/∂b² = ∂z✶/∂b ≥ 0 on a stable profile), so for a filter that averages at fixed z Jensen's inequality gives eₐˢ ≥ 0 pointwise, up to a marginal error where a filtered buoyancy falls between the profile's entries (the lookup then takes the nearest class). A filter with a vertical component mixes heights as well as buoyancies and needs extra steps, not currently implemented, to guarantee a non-negative eₐˢ, so filter has to act only in the horizontal and a filter with a vertical component throws an ArgumentError.

filter is any callable mapping a field to its low-pass-filtered counterpart, e.g. a reusable GaussianFilter or BoxFilter over dims = (1, 2). The filtered buoyancy and the filtered APE are materialized as Fields internally (so the separable filter takes its fast staged path), and the returned object is a lazy operation over them, ready for Field, Integral and OutputWriters. It lives at (Center, Center, Center), per unit mass (units m² s⁻²):

using Oceananigans, Oceanosticsgrid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b)filter = GaussianFilter(; dims=(1, 2), σ=0.1)SubFilterAvailablePotentialEnergy(model, filter)# outputSubFilterAvailablePotentialEnergy KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: subfilter_ape_ccc (generic function with 1 method)└── arguments: ("BinaryOperation",)└── computes: subfilter available potential energy  filter(eₐ) - eₐˡ

A convenience method SubFilterAvailablePotentialEnergy(model; σ, dims, boundary, N) builds the Gaussian filter for you along dims, (1, 2) by default, from a standard deviation σ (with σ = ℓ / (2√(2 ln 2)) for a FWHM ℓ). geopotential_height enters the buoyancy construction exactly as in reference_height.

source

Oceanostics.SubFilterAvailablePotentialEnergyEquation.SubFilterAvailablePotentialToKineticEnergyConversion — Type

SubFilterAvailablePotentialToKineticEnergyConversion(
    model,
    filter;
    method,
    geopotential_height
)

Return the subfilter-scale (SFS) conversion of available potential energy into kinetic energy τ(w, bᵣ), the rate at which the scales a low-pass filter removes release their APE to the subfilter flow:

    τ(w, bᵣ) = filter(w bᵣ) - w̄ b̄ᵣ ,   bᵣ = b - b✶(z) ,   b̄ᵣ = b̄ - b✶(z)

It is the subfilter half of the split whose filtered half is FilteredAvailablePotentialToKineticEnergyConversion w̄ b̄ᵣ: the two sum to filter(w bᵣ), so the subfilter and filtered budgets between them exchange exactly what the full field converts. It enters this budget as -τ(w, bᵣ) and the subfilter kinetic energy budget (SubFilterKineticEnergy) as +τ(w, bᵣ), so it is a reversible exchange rather than a source or a sink.

It is the subfilter covariance of w and bᵣ. For the horizontal filters this diagnostic accepts, b✶ depends on z alone and passes through the filter unchanged, so the b̄ᵣ = b̄ - b✶(z) the diagnostic forms is filter(bᵣ), and the reference profile cancels between the two terms, which leaves the subfilter buoyancy flux τ(w, bᵣ) = filter(w b) - w̄ b̄.

method has to be a ProfileLookup, for the reason SubFilterAvailablePotentialEnergy gives, and both halves are built on the one profile it supplies. filter has to act only in the horizontal, since the current implementation cannot guarantee a non-negative APE with a filter that has a vertical component; such a filter throws an ArgumentError. filter is any callable mapping a field to its low-pass-filtered counterpart, e.g. a reusable GaussianFilter or BoxFilter over dims = (1, 2). The result lives at (Center, Center, Center), per unit mass (units m² s⁻³):

using Oceananigans, Oceanosticsgrid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b)filter = GaussianFilter(; dims=(1, 2), σ=0.1)SubFilterAvailablePotentialToKineticEnergyConversion(model, filter)# outputSubFilterAvailablePotentialToKineticEnergyConversion KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: subfilter_ape_to_ke_conversion_ccc (generic function with 1 method)└── arguments: ("BinaryOperation",)└── computes: subfilter APE to KE conversion  filter(wbᵣ) - w̄ b̄ᵣ

A convenience method SubFilterAvailablePotentialToKineticEnergyConversion(model; σ, dims, boundary, N) builds the Gaussian filter for you along dims, (1, 2) by default, from a standard deviation σ (with σ = ℓ / (2√(2 ln 2)) for a FWHM ℓ).

source

Oceanostics.FilteredAvailablePotentialEnergyEquation.AvailablePotentialEnergyCrossScaleFlux — Type

AvailablePotentialEnergyCrossScaleFlux(
    model,
    filter;
    dims,
    method,
    geopotential_height
)

Return the cross-scale (scale-to-scale) available-potential-energy flux Πₐ, the rate at which a low-pass filter transfers available potential energy from the filtered to the subfilter scales:

    Πₐ = -τ(uᵢ, b) ∂ᵢΥˡ ,   τ(uᵢ, b) = filter(b uᵢ) - b̄ ūᵢ ,   Υˡ = z✶(b̄) - z

where τ(uᵢ, b) is the subfilter buoyancy flux (the subfilter_covariance of uᵢ and b) and Υˡ is the FilteredAvailablePotentialEnergyDisplacementPotential. It is the APE analogue of KineticEnergyCrossScaleFlux, with the subfilter buoyancy flux in place of the subfilter stress and ∇Υˡ in place of the resolved strain; Πₐ > 0 is forward (downscale, filtered → subfilter) transfer. It appears as -Πₐ in the filtered APE budget and as +Πₐ in the subfilter one, which is what makes it a transfer rather than a source or a sink. The result lives at (Center, Center, Center), per unit mass (units m² s⁻³).

method has to be a ProfileLookup, for the reason FilteredAvailablePotentialEnergy gives: Υˡ measures the filtered buoyancy against a profile it did not itself produce, ordinarily the sorted state of the full buoyancy. The lookup also makes z✶ a function of buoyancy alone, which is what differentiating Υˡ needs. filter has to act only in the horizontal, since the current implementation cannot guarantee a non-negative APE with a filter that has a vertical component; such a filter throws an ArgumentError.

dims selects which directions are summed over and which velocities are filtered: the default dims = (1, 2, 3) gives the full flux, while dims = (1, 3) gives the 2D x–z flux Πₐ = -(τ₁∂₁Υˡ + τ₃∂₃Υˡ). The filter's own directions are set independently inside filter, so the vertical term τ₃∂₃Υˡ stays in the sum although the filter acts only in the horizontal.

filter is any callable mapping a field to its low-pass-filtered counterpart, e.g. a reusable GaussianFilter or BoxFilter. The filtered fields and Υˡ are materialized internally, so the returned object is a lazy operation ready for Field, Integral and OutputWriters, re-filtering and re-sorting as the simulation evolves:

using Oceananigans, Oceanosticsgrid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b)filter = GaussianFilter(; dims=(1, 2), σ=0.1)AvailablePotentialEnergyCrossScaleFlux(model, filter)# outputAvailablePotentialEnergyCrossScaleFlux KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: cross_scale_ape_flux_ccc (generic function with 1 method)└── arguments: ("UnaryOperation",)└── computes: cross-scale available potential energy flux  -τ(uᵢ, b)∂ᵢΥˡ

A second method, AvailablePotentialEnergyCrossScaleFlux(model, filter, z✶ˡ; upsilon, dims), takes a reference height you built from the filtered buoyancy (see FilteredAvailablePotentialEnergy), so one lookup — and, through upsilon, one Υˡ (FilteredAvailablePotentialEnergyDisplacementPotential) — can serve this flux and the other filtered-state diagnostics. The filtered buoyancy is read off z✶ˡ itself; filter is still needed there, to build the subfilter buoyancy flux.

A convenience method AvailablePotentialEnergyCrossScaleFlux(model; σ, dims, boundary, N) builds the Gaussian filter for you from a standard deviation σ (with σ = ℓ / (2√(2 ln 2)) for a FWHM ℓ). There dims still selects the directions summed over, and the filter acts along the horizontal ones among them: the default filters along x and y and sums over all three directions, and dims = (1, 3) filters along x and sums over x and z.

source

Oceanostics.SubFilterAvailablePotentialEnergyEquation.SubFilterAvailablePotentialEnergyDissipationRate — Type

SubFilterAvailablePotentialEnergyDissipationRate(
    model,
    filter;
    method,
    geopotential_height
)

Return the subfilter-scale (SFS) available potential energy dissipation rate εₐˢ, the APE destruction by diffusion carried by the scales that a low-pass filter removes:

    εₐˢ = filter(εₐ) - εₐˡ ,   εₐˡ = -q̄ᵢ ∂ᵢΥˡ ,   q̄ᵢ = filter(qᵢ) ,   Υˡ = z✶(b̄) - z

where εₐ = -qᵢ∂ᵢΥ is the dissipation rate of the full field (AvailablePotentialEnergyDissipationRate), qᵢ is the closure's own diffusive buoyancy flux, and εₐˡ is the same contraction evaluated on the filtered state (FilteredAvailablePotentialEnergyDissipationRate): the filtered flux q̄ᵢ against the displacement potential Υˡ of the filtered buoyancy b̄ = filter(b). Filtering the flux rather than recomputing it from b̄ is what makes εₐˡ the sink of the filtered-state budget when κ varies in space; that docstring has the details. Both states are measured against one shared reference profile, exactly as in SubFilterAvailablePotentialEnergy, whose budget this is the diffusive sink of; it mirrors what SubFilterKineticEnergyDissipationRate is to the subfilter kinetic energy.

method has to be a ProfileLookup, for the reason SubFilterAvailablePotentialEnergy gives, and the lookup also makes each z✶ a function of buoyancy alone — the property that differentiating Υ and Υˡ needs (see AvailablePotentialEnergyDisplacementPotential). filter has to act only in the horizontal, since the current implementation cannot guarantee a non-negative APE with a filter that has a vertical component; such a filter throws an ArgumentError. Like AvailablePotentialEnergyDissipationRate, this diagnostic needs the buoyancy to be a tracer the closure diffuses (BuoyancyTracer only) and a closure that supplies a diffusive flux.

filter is any callable mapping a field to its low-pass-filtered counterpart, e.g. a reusable GaussianFilter or BoxFilter. The filtered fluxes, the filtered buoyancy, Υˡ and the filtered full-field dissipation are materialized as Fields internally, so the returned object is a lazy operation ready for Field, Integral and OutputWriters, recomputing (re-sorting included) as the simulation evolves. It lives at (Center, Center, Center), per unit mass (units m² s⁻³):

using Oceananigans, Oceanosticsgrid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b, closure=ScalarDiffusivity(κ=1e-4))filter = GaussianFilter(; dims=(1, 2), σ=0.1)SubFilterAvailablePotentialEnergyDissipationRate(model, filter)# outputSubFilterAvailablePotentialEnergyDissipationRate KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: subfilter_ape_dissipation_rate_ccc (generic function with 1 method)└── arguments: ("BinaryOperation",)└── computes: subfilter available potential energy dissipation rate  filter(εₐ) - εₐˡ

A convenience method SubFilterAvailablePotentialEnergyDissipationRate(model; σ, dims, boundary, N) builds the Gaussian filter for you along dims, (1, 2) by default, from a standard deviation σ (with σ = ℓ / (2√(2 ln 2)) for a FWHM ℓ).

source

Oceanostics.SubFilterAvailablePotentialEnergyEquation.SubFilterAvailablePotentialEnergyReferenceTendency — Type

SubFilterAvailablePotentialEnergyReferenceTendency(
    model,
    filter,
    ∂ₜb✶;
    method,
    geopotential_height
)

Return the subfilter-scale (SFS) reference tendency of the available potential energy Rˢ, the share of the evolving reference profile's effect on the APE that falls on the scales a low-pass filter removes:

    Rˢ = filter(R) - Rˡ ,   R = ∫_{z✶(b)}^{z} ∂ₜb✶(z̃) dz̃ ,   Rˡ = ∫_{z✶(b̄)}^{z} ∂ₜb✶(z̃) dz̃

where R is the reference tendency of the full field (AvailablePotentialEnergyReferenceTendency) and Rˡ that of the filtered buoyancy (FilteredAvailablePotentialEnergyReferenceTendency), both built on one time derivative ∂ₜb✶ of the shared reference profile; the full-field docstring shows how to build it and keep it current. Rˢ enters the budget of SubFilterAvailablePotentialEnergy as +Rˢ. R integrates to zero over a closed or periodic domain and Rˡ does not, so for a filter that preserves volume integrals ∫Rˢ dV = -∫Rˡ dV: an evolving reference profile moves APE across the filter scale.

method has to be a ProfileLookup, for the reason SubFilterAvailablePotentialEnergy gives. Passing ProfileLookup(z✶_column), with the column ∂ₜb✶ differentiates, shares that column's sort. filter has to act only in the horizontal, since the current implementation cannot guarantee a non-negative APE with a filter that has a vertical component; such a filter throws an ArgumentError. The result lives at (Center, Center, Center), per unit mass (units m² s⁻³):

using Oceananigans, Oceanosticsgrid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b)z✶_column = reference_height(model, method=VerticalSort())∂ₜb✶ = TimeDerivative(reference_buoyancy(z✶_column))filter = GaussianFilter(; dims=(1, 2), σ=0.1)SubFilterAvailablePotentialEnergyReferenceTendency(model, filter, ∂ₜb✶; method=ProfileLookup(z✶_column))# outputSubFilterAvailablePotentialEnergyReferenceTendency KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: subfilter_ape_reference_tendency_ccc (generic function with 1 method)└── arguments: ("BinaryOperation",)└── computes: subfilter available potential energy reference tendency  filter(R) - Rˡ

A convenience method SubFilterAvailablePotentialEnergyReferenceTendency(model, ∂ₜb✶; σ, dims, boundary, N) builds the Gaussian filter for you along dims, (1, 2) by default, from a standard deviation σ (with σ = ℓ / (2√(2 ln 2)) for a FWHM ℓ).

source