Divide a prescribed energy flux by the Exner function by ewquon · Pull Request #1020 · NumericalEarth/Breeze.jl

ewquon · 2026-09-23T05:23:09Z

Closes #976.

The conversion of a prescribed energy flux stopped one step short. 𝒬 / cᵖᵐ is a temperature flux
Jᵀ, and it was written into the ρθ budget as though it were a potential-temperature flux. The
missing step is Jᶿ = Jᵀ / Π. The function is named 𝒬_to_Jᶿ; it returned Jᵀ.

The interior ρE forcing and the radiative flux divergence already take both steps, dividing by
cᵖᵐ Π. This makes the flux path agree with them.

Effect

Potential temperature tendency unchanged at 1000 hPa, +4.7% at 850 hPa, +10.7% at 700 hPa.
Any coupled EarthSystemModel run picks this up through the energy flux field its coupler writes,
so surface heating over elevated terrain changes by that much and existing coupled baselines will
not reproduce.

Breaking

Both wrapper structs gain a standard_pressure field and a type parameter. They are exported from
Breeze.BoundaryConditions, so positional construction at the old arity now throws. Ordinary usage
is unaffected — the public one-argument EnergyFluxBoundaryCondition(flux) is unchanged, and
ThetaFluxBoundaryConditionFunction's shorter convenience constructor still works.

I deliberately did not add a convenience constructor preserving the old EnergyFlux arity: unlike
density, which resolves at runtime through the #777 Nothing fallback, there is nothing for
standard_pressure to fall back to, so omitting it should fail loudly rather than silently produce
a boundary condition carrying nothing.

standard_pressure converts to the grid's float type at materialization. This is load-bearing, not
defensive: boundary conditions materialize from the stub dynamics, before materialize_dynamics
normalizes, so on a Float32 grid the constructor would otherwise see a Float64 value and throw.

Version: strictly this is 0.12.0, since an exported struct's shape changed. Worth weighing
against the cost — NumericalEarth pins Breeze = "0.11", so a 0.12.0 release strands this fix
behind a separate NumericalEarth compat PR, and the fix is what makes coupled surface-energy results
interpretable. 0.11.3 with the behavior change called out in the release notes is the pragmatic
alternative. Maintainers' call.

Tests

The existing conversion test read the model's own boundary condition, and at the default base
pressure (Π = 1.003) it would have caught this change — the conventions differ by 3.3e-3, about ten
times isapprox's Float32 tolerance. It now runs two base pressures anyway, because a single
column cannot distinguish a Π that tracks the column's pressure from a constant that happens to be
right there. 𝒬 / cᵖᵐ is identical at both, so the ratio between them tests that the conversion
varies correctly, and depends on no tolerance.

The first commit is unrelated to the fix and lands first so history never goes red. getbc coverage for all boundary faces asserts Δρθ != 0 after one step at Δt = 1e-6, where the increment is
about 9.9e-9 against ρθ ≈ 353, whose Float32 spacing is 3.05e-5 — three thousand times
larger. Measured Δρθ is 9.85e-9 in Float64 and exactly 0.0 in Float32. That assertion could
only ever pass in single precision on rounding noise, and this change perturbs the flux enough to
flip it. Raising Δt to 1 s makes it measure what it claims. This has been latent because
test_float_types() returns (Float64,) unless BREEZE_TEST_FLOAT32=true, which is set nowhere in
the repo.

Verified

forcing_and_boundary_conditions.jl and balance_adiabatically.jl green in both precisions;
quality_assurance.jl (Aqua, ExplicitImports) clean. An independent review run covered
boundary_conditions, wall_fluxes, turbulence_closures, atmosphere_model_construction and
diagnostics — 352/352 — and confirmed zero allocations and concrete return types, including a
mixed Float32-grid / Float64-constants case.

GPU: not run, but the path is precedented. compute_potential_temperature_tendency! already
builds a LiquidIcePotentialTemperatureState and calls exner_function on it inside a kernel
launched on arch (potential_temperature_tendency.jl:93-95), and the bulk-flux boundary
conditions already read dynamics_fields.p on-device. This adds no construct that is not already
exercised there — the added field is an isbits scalar, and Adapt.adapt_structure is updated for
both structs. What has not been run is this specific call site, since GPU jobs only trigger on
pull_request.

Coordination

Reversible leapfrog + digital-filter adiabatic balance (na_init) #834 reconstructs both structs with five positional arguments in without_microphysics and will
break against this. That is intended: it needs to carry standard_pressure through.
Fix moisture partitioning and unify moist-state conversions #997 rewrites the same using ..Thermodynamics: import block this touches.

Docs: total_energy_density_name's docstring and the boundary-conditions page both stated the
flux/forcing asymmetry as intended behavior; both now state the single conversion. Jᶿ joins the
notation table, since Jᵀ was there but Jᶿ was not and they differ by exactly this factor.

The last commit renames 𝒬_to_Jᶿ → 𝒬ᵀ_to_Jᶿ and Jᶿ_to_𝒬 → Jᶿ_to_𝒬ᵀ. A bare 𝒬 does not say
which energy flux, and a latent heat flux would be wrong here — that energy travels with the
moisture flux. The valid input is an enthalpy flux, which is 𝒬ᵀ as the table defines it, and the
table gives every 𝒬 a superscript. That is also why no bare-𝒬 row was added: with the rename,
none remains in the source. Both functions are internal to the one file and not exported, so the
rename is revertible in a pass if you prefer the bare spelling.

🤖 Generated with Claude Code

https://claude.ai/code/session_01D2fuRfqFoMFjvYyx3idQtQ

`getbc coverage for all boundary faces` asserts that one step with a prescribed `ρE` surface flux moves `ρθ`. At Δt = 1e-6 the increment 𝒬 Az Δt / (cᵖᵐ V) is about 9.9e-9 against ρθ ≈ 353, whose Float32 spacing is 3.05e-5 — three thousand times larger. Measured Δρθ is 9.85e-9 in Float64 and exactly 0.0 in Float32, so in single precision the assertion could only ever pass on rounding noise from the rest of the tendency. This has been latent because `test_float_types()` returns `(Float64,)` unless `BREEZE_TEST_FLOAT32=true`, which is set nowhere in the repo, so CI has never exercised it. Raising Δt to 1 s puts the increment at 9.9e-3, two orders above the spacing, and the test then measures what it claims to in both precisions. One step of an anelastic single-cell model has no stability constraint that 1 s approaches. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01D2fuRfqFoMFjvYyx3idQtQ

An energy input reached the `ρθ` budget through two different conversions depending on how it arrived. The interior `ρE` forcing and the radiative flux divergence divide by `cᵖᵐ Π`; a prescribed flux under the same key divided by `cᵖᵐ` alone. Equating the two tendencies fixes the conversion: Oceananigans adds `getbc * Az / V` to `Gρθ`, while a volumetric source of the same strength enters as `FρE / (cᵖᵐ Π)`, so a flux must enter as `𝒬 / (cᵖᵐ Π)`. Physically, with `𝒬 = ρ cᵖᵐ ⟨w'T'⟩` and `T = Π θ + (ℒˡᵣ qˡ + ℒⁱᵣ qⁱ) / cᵖᵐ`, holding moisture fixed gives `δT = Π δθ`, so the condensate term cancels and the relation holds saturated as well as dry. Closes #976. Π comes from `dynamics_pressure`, reached through the `dynamics_fields` argument `getbc` already carried but these methods ignored. That is bit-identical to what the anelastic `diagnose_thermodynamic_state` reads, so flux and forcing divide by the same number. The compressible core instead diagnoses a density state whose pressure is ρRᵐT — which is what `dynamics.pressure` holds, written by the same inversion during `update_state!`, so the two agree up to its staleness inside an RK stage. Both wrapper structs gain `standard_pressure`, converted to the grid's float type at materialization: boundary conditions materialize from the stub dynamics, before `materialize_dynamics` normalizes, so on a Float32 grid the constructor would otherwise see a Float64 pˢᵗ and throw. Omitting the field is left a MethodError rather than given a convenience constructor — unlike `density`, which resolves at runtime through the #777 `Nothing` fallback, there is nothing to fall back to. `Jᶿ_to_𝒬` takes the same factor so the inverse stays the inverse; a flux supplied under `ρE` is unwrapped rather than reconverted, and round-trips exactly. Measured: the flux is unchanged at 1000 hPa, 4.7% larger at 850 hPa and 10.7% larger at 700 hPa. Every coupled NumericalEarth run picks this up through the `Jᵉ` field its coupler writes, so surface heating over elevated terrain changes by that much and existing baselines no longer reproduce. The test now runs two base pressures. A single column would have caught this change — at the default base pressure Π = 1.003, a relative difference of 3.3e-3, about ten times isapprox's Float32 tolerance — but it cannot distinguish a Π that tracks the column's pressure from a constant that happens to be right there. `𝒬 / cᵖᵐ` is identical at both pressures, so the ratio between them tests that the conversion varies correctly, and being a ratio it depends on no tolerance. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01D2fuRfqFoMFjvYyx3idQtQ

codecov · 2026-09-23T05:55:13Z

Codecov Report

❌ Patch coverage is 84.84848% with 5 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
...c/BoundaryConditions/thermodynamic_variable_bcs.jl	84.84%	5 Missing ⚠️

📢 Thoughts on this report? Let us know!

#974 wrote the flux/forcing asymmetry into the docs as intended behavior: `total_energy_density_name` said an energy input is divided by `cᵖᵐ` for fluxes and `cᵖᵐ Π` for forcings, and the boundary-condition page justified the split by noting the forcing "is applied to a potential temperature rather than a temperature" — which is equally true of the flux. Both now state the single conversion. The page's derivation also went from a dynamic flux to a kinematic one, dropping the ρ that `Jᶿ` carries, and asserted `θ = T / Π`, which holds only without condensate. The differential form is both correct and stronger. `Jᶿ` joins the notation table. It had `Jᵀ` but not `Jᶿ`, and they differ by exactly the Exner factor, so the distinction is now load-bearing. A bare `𝒬` needs no row: every `𝒬` in that table carries a superscript naming its flux, and `𝒬 = cᵖᵐ Π Jᶿ = cᵖᵐ Jᵀ` is the `𝒬ᵀ` already there. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01D2fuRfqFoMFjvYyx3idQtQ

ewquon · 2026-09-23T15:42:54Z

Another check to make sure this fix doesn't break anything else: what does the coupler actually put under ρE?

I.e., if it supplied cᵖᵐ Jᶿ rather than cᵖᵐ Jᵀ, then dividing by cᵖᵐ Π would yield Jᶿ / Π and
introduce an error of exactly the size this PR claims to remove, in every coupled run.

Traced through NumericalEarth (origin/main):

# InterfaceComputations/interface_states.jl:308
#   "Temperature increment including the ``lapse rate'' α = g / cᵖᵐ"
surface_atmosphere_temperature(Ψₐ, ℙₐ) = Tᵃᵗ + g * Δh / cᵃᵗ

# InterfaceComputations/compute_interface_state.jl:104
θᵃᵗ = surface_atmosphere_temperature(atmosphere_state, atmosphere_properties)
Δθ = θᵃᵗ - Tₛ

# InterfaceComputations/coefficient_based_turbulent_fluxes.jl:359
θ★ = Ch / sqrt(Cd) * Δθ

# InterfaceComputations/atmosphere_interface_kernels.jl:56
interface_fluxes.sensible_heat[i, j, 1] = - ρᵃᵗ * cᵖᵐ * u★ * θ★

and that value reaches the ρE bottom-flux field via _assemble_net_atmosphere_fluxes!
(ext/NumericalEarthBreezeExt/breeze_atmosphere_interface.jl), for both the atmosphere–ocean and
atmosphere–land interfaces.

θᵃᵗ is a potential temperature — lifting dry-adiabatically is what makes one, and it is
equivalently dry static energy over cᵖᵐ, (cᵖᵐ T + g z) / cᵖᵐ. So the question is not
temperature versus potential temperature. It is the reference level.

Δh = zᵃᵗ is the MOST reference height: the lowest cell-centre elevation above ground,
½·Δz(i,j,1). The z = 0 in the surface_atmosphere_temperature(Ψₐ, ℙₐ) comment is the
surface-layer coordinate's origin, not sea level — worth being explicit about, since lifting
from a sea-level height while differencing against a skin temperature on 2.5 km terrain would be a
~24 K mismatch at the dry adiabatic lapse rate.

So θᵃᵗ is lifted to the local ground surface and differenced against Tₛ there.
A potential temperature referenced to the ground is numerically the temperature
there, so Δθ is a temperature difference at the surface, θ★ a temperature scale, and
− ρ cᵖᵐ u★ θ★ is ρ cᵖᵐ ⟨w′T′⟩ = 𝒬ᵀ.

Breeze's Π = (p / pˢᵗ)^κ references θ to pˢᵗ, not to the ground. These are two different
potential temperatures, and Π is exactly the conversion between them:

NE produces a flux of θ referenced to the surface, i.e. 𝒬ᵀ. Dividing by cᵖᵐ Π gives the
flux of θ(pˢᵗ), which is what ρθ carries. ✓
Had NE referenced to pˢᵗ, its flux would already be cᵖᵐ Jᶿ and dividing again would
double-count. ✗

Only the reference level separates those two cases.

Therefore 𝒬 / (cᵖᵐ Π) = Jᵀ / Π = Jᶿ, which is what the ρθ budget wants.

A naming question this raises. The function is 𝒬_to_Jᶿ, and a bare 𝒬 does not say which
energy flux — handing it a latent heat flux 𝒬ᵛ would be wrong.

🤖 Generated with Claude Code

https://claude.ai/code/session_01D2fuRfqFoMFjvYyx3idQtQ

`𝒬_to_Jᶿ` took a bare `𝒬`, which does not say which energy flux. Handing it a latent heat flux would be wrong — that energy is carried by the moisture flux, as the coupler that writes this field notes — so the valid input is specifically an enthalpy flux, `cᵖᵐ` times a temperature flux. That is `𝒬ᵀ` as the notation table defines it, and the table gives every `𝒬` a superscript naming its flux: `𝒬ᵀ = cᵖᵐ Jᵀ`, `𝒬ᵛ = ℒˡ Jᵛ`. What arrives under `ρE` from a coupled model is `𝒬ᵀ`. NumericalEarth's similarity theory forms `θᵃᵗ` by lifting the lowest-level air dry-adiabatically through the MOST reference height — the cell-centre elevation above ground, per column on a terrain-following grid — and differences it against the skin temperature at that same ground. A potential temperature referenced to the ground is numerically the temperature there, so `- ρ cᵖᵐ u★ θ★` is `ρ cᵖᵐ ⟨w'T'⟩`. Had it instead been referenced to `pˢᵗ`, the flux would already be `cᵖᵐ Jᶿ` and the preceding commit's division by `Π` would double-count. Same class of correction as that commit, applied to the other end of the signature: the function returned `Jᵀ` while claiming `Jᶿ`; it accepted `𝒬ᵀ` while claiming an unspecified `𝒬`. Mechanical and confined — `𝒬ᵀ_to_Jᶿ` and `Jᶿ_to_𝒬ᵀ` are internal to this file and not exported, so a maintainer preferring the bare spelling can revert it in one pass. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01D2fuRfqFoMFjvYyx3idQtQ

The six `getbc` methods differ in exactly one thing: which boundary cell they hand to the conversion — `(i,j,1)` for bottom, `(i,j,Nz)` for top, `(1,j,k)` for west, and so on. Nothing asserted that. `getbc coverage for all boundary faces` covers bottom and west, and asserts only that `ρθ` moved, which a conversion reading the wrong cell also satisfies. Adding `Π` raised the cost of getting that index wrong. It previously reached only `qᵛ` and density, which vary weakly, so a misread cell was a percent. `Π` follows pressure: over the column used here it runs ≈0.97 at the lowest cell to ≈0.53 at the highest, so the same mistake is now most of a factor of two. So this checks the value rather than the motion: for each of the six faces, that the boundary condition returns `𝒬ᵀ / (cᵖᵐ Π)` with `Π` taken at that face's own cell, on a fully bounded 15 km column where the faces genuinely disagree. That premise is asserted first rather than last: if `Π` barely varied, the six would pass under any index, so it belongs before them, not after. No `θ` is set. Nothing the conversion reads depends on it — the density and pressure are the anelastic reference fields, built at construction, and `exner_function` ignores the state's `θ`, which is the property the whole conversion rests on. Setting it would imply otherwise. That holds only while the default microphysics is `nothing`: with one, `set!` would split `qᵗ` into condensate using temperature and the vapor-only expectation here would stop matching. What it does not catch: the anelastic reference pressure is a column, so `Π` has no horizontal structure and a pure i↔j or 1↔N mix-up between two horizontal faces is invisible. Errors in the vertical index, and errors confusing a horizontal index with it, do show. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01D2fuRfqFoMFjvYyx3idQtQ

glwagner · 2026-09-23T19:47:43Z

+        Jᶿ = Array(interior(Field(BoundaryConditionOperation(ρθ, side, model))))
+        expected = [𝒬ᵀ / (cᵖᵐ * Π_cell(c...)) for c in getproperty(cells, side)]
+        @test vec(Jᶿ) ≈ vec(expected)


may be possible to use more of the built in Field infra here, eg avoiding Array(interior()) chain. Typical claudism.

glwagner · 2026-09-23T19:48:47Z

+    p = Array(interior(dynamics_pressure(model.dynamics)))
+    Nx, Ny, Nz = size(grid)
+
+    # The anelastic reference pressure is a column, so `Π` varies in `k` only. That is enough to
+    # catch any error in the vertical index, and any error that confuses a horizontal index with
+    # it; a pure i↔j or 1↔N mix-up between two horizontal faces would not show here.
+    p_cell(i, j, k) = FT(p[min(i, size(p, 1)), min(j, size(p, 2)), k])
+    Π_cell(i, j, k) = exner_function(LiquidIcePotentialTemperatureState(zero(FT), q, pˢᵗ, p_cell(i, j, k)),
+                                     constants)
+
+    # Premise for everything below: `Π` must vary enough across the column that reading the wrong
+    # cell is visible. Without it the six assertions would pass under any index.
+    @test Π_cell(1, 1, 1) / Π_cell(1, 1, Nz) > 1.5


it may be possible to redesign this stuff to use more built-in Field infra

however it doesn't matter. we can probably do a sweep later to clean this (and many other) up.

glwagner · 2026-09-23T19:55:18Z

    q = grid_moisture_fractions(i, j, k, grid, ef.microphysics, ρ, qᵛ, fields)
    cᵖᵐ = mixture_heat_capacity(q, ef.thermodynamic_constants)
-    return 𝒬 / cᵖᵐ
+    Π = near_wall_exner_function(i, j, k, ef, q, dynamics_fields)


might it be possible to use a more generic exner_function(i, j, k, grid, ef, q, dynamics_field) ?

also note to include the grid argument so this is compatible with KernelFunctionOperation

glwagner · 2026-09-23T19:55:51Z

+end
+
+# Convert energy flux to potential temperature flux: Jᶿ = 𝒬ᵀ / (cᵖᵐ Π)
+@inline function 𝒬ᵀ_to_Jᶿ(i, j, k, grid, ef, 𝒬ᵀ, fields, dynamics_fields)


could be possible to share a utility with the energy forcing conversion for tendencies. but maybe not.

@glwagner

Review points from @glwagner: the helper was called `near_wall_exner_function` but is valid at any `i, j, k`, and a grid-point Exner evaluation should carry `grid` so it is usable from a `KernelFunctionOperation`. Both resolve the same way. It is not a boundary-conditions helper that happens to compute Π — it is Π at a grid point, so it becomes a method of `exner_function` itself with the signature `exner_function(i, j, k, grid, ef, q, dynamics_fields)`. Both call sites already had `grid` in scope. The `near_wall_` prefix was wrong for the reason given. Its neighbours earn it — `near_wall_velocity` reads `fields.u[i, j, 1]`, fixed to the first cell — while this takes an arbitrary index and the six `getbc` methods are what choose a boundary cell. The prefix described the callers, not the function. Extending a name imported through `using ... :` is not allowed, so `exner_function` moves to its own `import` line. Not piracy: `ef` is a type this module owns. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01D2fuRfqFoMFjvYyx3idQtQ

ewquon and others added 2 commits September 22, 2026 00:10

ewquon force-pushed the eq/energy-flux-exner branch from 49964c5 to 059bcee Compare September 23, 2026 05:56

ewquon and others added 2 commits September 23, 2026 10:48

ewquon force-pushed the eq/energy-flux-exner branch from b7607ba to a10cb79 Compare September 23, 2026 18:26

Merge branch 'main' into eq/energy-flux-exner

6881962

ewquon requested a review from glwagner September 23, 2026 19:45

glwagner approved these changes Sep 23, 2026

View reviewed changes

glwagner reviewed Sep 23, 2026

View reviewed changes

Comment thread src/BoundaryConditions/thermodynamic_variable_bcs.jl Outdated

glwagner reviewed Sep 23, 2026

View reviewed changes

ewquon and others added 2 commits September 23, 2026 14:24

Merge branch 'main' into eq/energy-flux-exner

f31ab9e

ewquon merged commit b6f6ff5 into main Sep 23, 2026
23 of 24 checks passed

ewquon deleted the eq/energy-flux-exner branch September 23, 2026 21:54

Navigation Menu

Divide a prescribed energy flux by the Exner function - #1020

Divide a prescribed energy flux by the Exner function#1020
ewquon merged 8 commits into
mainfrom
eq/energy-flux-exner

ewquon commented Sep 23, 2026 •

edited

Loading

Uh oh!

codecov Bot commented Sep 23, 2026 •

edited

Loading

Uh oh!

ewquon commented Sep 23, 2026 •

edited

Loading

Uh oh!

glwagner Sep 23, 2026

Uh oh!

glwagner Sep 23, 2026

Uh oh!

Uh oh!

glwagner Sep 23, 2026

Uh oh!

glwagner Sep 23, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

ewquon commented Sep 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Effect

Breaking

Tests

Verified

Coordination

Uh oh!

codecov Bot commented Sep 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

ewquon commented Sep 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

glwagner Sep 23, 2026

Choose a reason for hiding this comment

Uh oh!

glwagner Sep 23, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

glwagner Sep 23, 2026

Choose a reason for hiding this comment

Uh oh!

glwagner Sep 23, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

ewquon commented Sep 23, 2026 •

edited

Loading

codecov Bot commented Sep 23, 2026 •

edited

Loading

ewquon commented Sep 23, 2026 •

edited

Loading