Materialize a time-varying pressure-level geopotential once per state update

src/Atmospheres/prescribed_atmosphere.jl

@@ -158,6 +158,10 @@ end
     for fts in ftses
         update_field_time_series!(fts, time)
     end
+
+    # A no-op unless the grid's vertical is a time-varying geopotential.
+    NumericalEarth.Grids.materialize_geopotential!(atmos.grid)
+
     return nothing
 end
 

src/DataWrangling/ERA5/ERA5_pressure_levels.jl

@@ -254,8 +254,9 @@ end
 Build a [`PressureLevelVerticalDiscretization`](@ref) whose per-column heights
 follow the ERA5 geopotential in time. The geopotential Φ(λ,φ,p,t) is loaded over
 the whole date window of `metadata` (all snapshots in memory) and wrapped in a
-`TimeSeriesInterpolation` bound to `clock`; `znode` reads Φ interpolated to
-`clock.time` and divides by `g`. Sub-surface levels are clipped to the local
+`TimeSeriesInterpolation` bound to `clock`, which the discretization materializes
+into a `Field` refreshed at each `update_state!(::PrescribedAtmosphere)`; `znode`
+reads that snapshot and divides by `g`. Sub-surface levels are clipped to the local
 surface geopotential so that columns stay monotonic. The single-level surface
 geopotential (orography × g, static in time) is downloaded from the matching
 `ERA5*SingleLevel` dataset and used as the clip source.

src/Grids/pressure_level_vertical_discretization.jl

@@ -34,9 +34,12 @@ Two `znodes` paths:
   locations, e.g. after `mean(...; dims=(1, 2))`); otherwise it returns
   a 3-D `Field` of per-cell heights at the field's own location.
 
-`geopotential` (units m²/s²) is a 3-D `Field` or a `TimeSeriesInterpolation`
-over a `FieldTimeSeries`. The former gives a static z-coordinate; the latter
-gives a time-evolving one driven by an attached `Clock`.
+`geopotential` (units m²/s²) is always a 3-D `Field`. Passing a `Field` to the
+constructor gives a static z-coordinate. Passing a `TimeSeriesInterpolation`
+over a `FieldTimeSeries` gives a time-evolving one driven by an attached
+`Clock`: the constructor materializes it into a `Field` that
+[`materialize_geopotential!`](@ref) refreshes from
+`update_state!(::PrescribedAtmosphere)`.
 
 The `LatitudeLongitudeGrid` constructor needs a value for `Lz`; we compute it
 as `extrema(geopotential) / g` inside `generate_coordinate`.
@@ -64,12 +67,18 @@ sub-surface *data* stays on those levels, so `column_fractional_z_index` also
 clamps interpolation to the first above-ground level (never sampling it). The
 clip source is retained on the discretization and exposed through
 [`surface_elevation`](@ref).
+
+A `TimeSeriesInterpolation` `geopotential` is clipped before it is materialized,
+so every snapshot the discretization goes on to compute inherits the clipped
+columns.
 """
 function PressureLevelVerticalDiscretization(geopotential;
                                               gravitational_acceleration,
                                               surface_geopotential = nothing)
     isnothing(surface_geopotential) || clip_subsurface!(geopotential, surface_geopotential)
-    return PressureLevelVerticalDiscretization(gravitational_acceleration, geopotential, surface_geopotential)
+    return PressureLevelVerticalDiscretization(gravitational_acceleration,
+                                               geopotential_field(geopotential),
+                                               surface_geopotential)
 end
 
 # Skip the generic validator (which would `length`-check the missing 1-D fields).
@@ -95,12 +104,26 @@ function generate_coordinate(FT, topo, sz, halo,
     return Lz, arch_discretization
 end
 
+"""
+    MaterializedGeopotential
+
+A geopotential `Field` computed from a `TimeSeriesInterpolation`: Φ at the time
+of the last [`materialize_geopotential!`](@ref), stored as plain data, with the
+source series still reachable through the field's `operand`.
+"""
+const MaterializedGeopotential = Field{<:Any, <:Any, <:Any, <:TimeSeriesInterpolation}
+
+geopotential_field(Φ::Field) = Φ
+geopotential_field(Φ::TimeSeriesInterpolation) = Field(Φ)
+
 geopotential_data_for_extrema(Φ::Field) = interior(Φ)
 # Use `interior(fts)` — not `parent(fts)` — so halo zeros don't dominate the
 # extrema / column mean. For a TSI this reads every time slice, so `Lz` spans the
 # heights reached at any time in the window — the extent a child must be able to
 # bracket over the whole run. One-time setup, not a hot path.
 geopotential_data_for_extrema(Φ::TimeSeriesInterpolation) = interior(Φ.time_series)
+# A materialized snapshot holds one time; read the source series behind it instead.
+geopotential_data_for_extrema(Φ::MaterializedGeopotential) = geopotential_data_for_extrema(Φ.operand)
 
 Adapt.adapt_structure(to, z::PressureLevelVerticalDiscretization) =
     PressureLevelVerticalDiscretization(z.gravitational_acceleration,
@@ -126,6 +149,21 @@ Type alias for any underlying grid whose vertical coordinate is a
 const PressureLevelGrid =
     AbstractUnderlyingGrid{<:Any, <:Any, <:Any, <:Any, <:PressureLevelVerticalDiscretization}
 
+"""
+    materialize_geopotential!(grid)
+
+Recompute the geopotential snapshot that a [`PressureLevelGrid`](@ref) with a
+time-varying vertical reads, at the current time of the `Clock` its
+`TimeSeriesInterpolation` is bound to. `znode` returns Φ at the time of the last
+refresh, not at the instant of access; `update_state!(::PrescribedAtmosphere)`
+refreshes.
+
+A no-op for any other grid, and for a static-`Field` geopotential.
+"""
+materialize_geopotential!(grid) = nothing
+materialize_geopotential!(grid::PressureLevelGrid) = materialize_geopotential!(grid.z.geopotential)
+materialize_geopotential!(Φ::MaterializedGeopotential) = (compute!(Φ); nothing)
+
 # Override the LLG show, which reads `grid.z.cᵃᵃᶠ` and crashes on PLVD.
 # We print the horizontal axes and report the z coordinate via the
 # `PressureLevelVerticalDiscretization` `show` defined above. Targeting the