polaris.tasks.ocean.horiz_press_grad.column.get_pchip_layer_mean

polaris.tasks.ocean.horiz_press_grad.column.get_pchip_layer_mean(z_tilde_nodes: ndarray, values_nodes: ndarray, name: str) → Callable[[ndarray, ndarray], ndarray]

Create a function giving the exact layer mean of a monotone PCHIP profile over a layer’s pseudo-height range.

Omega’s prognostic conservative temperature and absolute salinity are layer means, and the higher-order pressure gradient reconstructs them as mean-preserving polynomials in pressure, so the initial condition must supply layer means rather than samples at layer midpoints. Within a column pressure is exactly linear in pseudo-height, so a mean over a layer’s pseudo-height range is also its mean over the layer’s pressure range and there is no ambiguity about which variable to average in.

The mean is computed by two-point Gauss-Legendre quadrature over each profile-node interval the layer overlaps. PCHIP is piecewise cubic and a two-point Gauss rule is exact through cubics, so this is exact rather than approximate, and it is preferred over differencing the interpolant’s antiderivative – also exact in principle – because that differences two large numbers and so loses precision as the layers thin. Measured on the linear profile of hydrostatic_consistency_linear, the antiderivative form departs from the exact layer mean by an amount growing like 1/h (7e-14 at 256 m layers, 1.3e-12 at 32 m) while this form sits at 0.9 * eps of the salinity itself at every resolution. That matters because the departure limits how exactly the initial condition sits in the scheme’s exact set, and hence the floor of the machine-precision test.

Parameters:

• z_tilde_nodes (np.ndarray) – One-dimensional z-tilde node locations. Must be strictly monotonic (increasing or decreasing).

• values_nodes (np.ndarray) – One-dimensional values at z_tilde_nodes.

• name (str) – A descriptive name used in error messages.

Returns:

layer_mean (callable) – A function of (z_tilde_top, z_tilde_bot) – the pseudo-height of each layer’s upper and lower interface – returning the layer means. Layers must lie within the node range, but interfaces are allowed to fall a round-off distance outside it, which is unavoidable because a column’s outermost interfaces sit exactly on the outermost nodes in several configurations.