# P-Star Coordinate and Tracer Initialization

date: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

## Background: z-tilde and p-star in Omega

Omega's vertical coordinate system is built on **pseudo-height**,
$\tilde{z} = -p / (\rho_0 g)$, where $p$ is sea gauge pressure, $\rho_0$ is a reference
seawater density, and $g$ is gravitational acceleration. This family of coordinates is
called the **z-tilde coordinate** in Omega.

Within the z-tilde framework, Omega supports a specific ALE variant called the **p-star
coordinate**. The p-star coordinate defines a set of reference pseudo-thicknesses
`RefPseudoThickness` — the pseudo-thicknesses the layers would have if the sea-surface
pressure were zero — together with a set of `VertCoordMovementWeights` that control how
changes in total column mass are distributed across layers. At runtime, Omega's
`VertCoord::computeTargetThickness()` adjusts each layer's pseudo-thickness from its
reference value proportionally to the total column pressure difference and the layer's
movement weight. With uniform weights (all ones, as currently used), every layer scales
proportionally and the coordinate behaves as a pseudo-compressible z-star: the total
column pseudo-thickness tracks the total column pressure.

Polaris initializes the p-star coordinate by computing `RefPseudoThickness` from a
target bathymetric depth and a given ocean state. The general z-tilde conversion
utilities (pressure/pseudo-height conversions, geometric height recovery) live in
`polaris/ocean/vertical/ztilde.py` and are reusable for any z-tilde–based coordinate.
The p-star–specific coordinate construction lives in `polaris/ocean/vertical/pstar.py`,
and the iterative initialization algorithm lives in `polaris/ocean/vertical/pstar_init.py`.

**MinLayerCell scope**: all p-star columns are assumed to start at the surface
(`minLevelCell = 0` in zero-based indexing). Support for non-zero `minLevelCell`
(e.g., ice-shelf cavities where the water column does not start at the free surface) is
explicitly **out of scope** for this design.

## Summary

Because ocean observations record geometric depth rather than pseudo-height, establishing
`RefPseudoThickness` requires knowing the thermohaline structure of the water column and
the equation of state (EOS). In turn, the thermohaline structure must be sampled on the
vertical grid that depends on `RefPseudoThickness`. This circular dependency means the
vertical coordinate and the tracer state must be determined jointly through iteration.

This design document describes:

1. A fixed-point iteration for determining `BottomPressure` (and therefore
   `RefPseudoThickness` and the p-star coordinate) such that the geometric seafloor depth
   recovered from the coordinate matches a target bathymetric depth within a configurable
   tolerance.
2. A pluggable interface for supplying conservative temperature (CT) and absolute salinity
   (SA) at each outer iteration step, supporting multiple initialization strategies
   without changes to the core algorithm.
3. The complete set of variables produced at convergence and their distribution across the
   separate output files established by the Polaris ocean framework.
4. The interaction between the proportional-ratio update and full or partial bottom cell
   snapping.

This capability is needed wherever the p-star coordinate must be initialized from a
known geometric seafloor depth. Its first concrete users are idealized test cases such as
`ocean/horiz_press_grad`, which already implements the algorithm but embeds it in
task-specific code, and the planned realistic initialization task. The design is
successful if the iteration logic is cleanly separated from the CT/SA initialization
strategy, and if the same base class and algorithm serve both idealized and realistic use
cases without modification.

The primary software challenge is providing a generalization pattern consistent with the
existing Polaris framework — specifically the convention of using overridable methods on
step subclasses, as exemplified by the convergence task framework — while keeping the
coupling between p-star coordinate construction, tracer initialization, and EOS
evaluation explicit and inspectable.

## Requirements

### Requirement: The p-star vertical coordinate can be initialized with a geometric seafloor depth matching target bathymetry within a configurable tolerance

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

Given a target geometric seafloor depth (equivalently, a target geometric water-column
thickness) and a means of obtaining conservative temperature and absolute salinity on the
resulting vertical grid, Polaris shall be able to initialize the p-star vertical
coordinate such that the geometric seafloor depth recovered from the converged coordinate
matches the target to within a configurable fractional tolerance.

The workflow shall support a configurable maximum number of iterations and shall terminate
early when the fractional change in the recovered geometric water-column thickness falls
below the tolerance. Both the convergence tolerance and the maximum iteration count shall
be configurable through the standard Polaris configuration mechanism.

### Requirement: CT and SA can be initialized on the p-star coordinate through a pluggable interface

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The workflow shall define a clear interface by which conservative temperature and absolute
salinity are supplied at each outer iteration step. This interface shall be general enough
to support multiple initialization strategies — including analytic profiles, piecewise
profiles defined by configuration values, and interpolation from an observational
hydrography dataset — without any change to the core iteration logic. The interface shall
be consistent in style with the generalization mechanisms already used in the Polaris
framework.

### Requirement: The initialized state includes all p-star coordinate variables and derived geometric quantities needed for subsequent steps

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The workflow shall produce at convergence all variables needed to configure the ocean
model, including p-star coordinate fields (pseudo-thickness, reference pseudo-thickness,
pseudo-height at midpoints and interfaces, minimum and maximum valid level indices, cell
mask, coordinate movement weights), converged tracer fields, specific volume, pressure,
`bottomDepth` consistent with the converged geometric water-column thickness, and
geometric height at layer midpoints and interfaces.

These outputs shall be distributed across model-specific files according to the Polaris
ocean framework's split-file convention: vertical coordinate variables shall be written to
a separate vertical-coordinate file for Omega (where they feed the `InitialVertCoord`
stream) or kept in the initial-state file for MPAS-Ocean. The outputs shall be sufficient
for subsequent forward-model steps without re-running the iteration.

### Requirement: Full and partial bottom cells are handled correctly within the iteration

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

When full or partial bottom cell behavior is requested, the iteration shall handle the
discrete snapping of the pseudo-column-thickness to layer boundaries in a manner that is
consistent with the converged state. The design shall specify how to detect and handle the
case where snapping prevents further convergence.

### Requirement: The capability is reusable across idealized and realistic initialization tasks

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The same iteration algorithm shall be applicable to both idealized test cases and
observationally-constrained realistic initialization tasks without duplication of the core
loop. Task-specific initialization of CT and SA shall be separable from the shared
iteration logic.

## Algorithm Design

### Algorithm Design: The p-star vertical coordinate can be initialized with a geometric seafloor depth matching target bathymetry within a configurable tolerance

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The core algorithm is a fixed-point iteration in pseudo-height space. Let:

- $\eta_0$ = prescribed sea-surface height (m); defaults to $-P_\text{surf}/(\rho_0 g)$,
  the resting surface-pressure depression for a reference-density fluid
- $H_\text{geo}^\star$ = target geometric water-column thickness $= \eta_0 - z_\text{bot}$
- $P_\text{surf}$ = prescribed sea-surface pressure (Pa; due to atmosphere,
  sea ice, ice shelves, etc.)
- $P_\text{bot}^{(k)}$ = `BottomPressure` at iteration $k$
- $\hat{P}_\text{bot}^{(k)}$ = post-snap `BottomPressure` after applying full or partial
  bottom-cell constraints inside the coordinate-construction step
- $H_\text{geo}^{(k)}$ = geometric water-column thickness recovered at iteration $k$

The three constraints ($z_\text{bot}$, $P_\text{surf}$, $\eta_0$) are fixed; the iteration
adjusts `BottomPressure` until all three are satisfied simultaneously.  The initial
`BottomPressure` is computed from the reference density $\rho_0$ and gravitational
acceleration $g$:

$$P_\text{bot}^{(0)} = P_\text{surf} + \rho_0\, g\, H_\text{geo}^\star$$

For the default $\eta_0 = -P_\text{surf}/(\rho_0 g)$ this simplifies to
$P_\text{bot}^{(0)} = \rho_0\, g\,(-z_\text{bot})$ because the $P_\text{surf}$ terms
cancel.  This keeps `pseudo_bottom_depth` $= P_\text{bot}^{(0)}/(\rho_0 g)$ equal to the
bathymetric depth $-z_\text{bot}$, safely within the reference grid regardless of
$P_\text{surf}$.

$P_\text{surf}$ and $\eta_0$ are optional inputs.  The diagnostic `ssh` at convergence
equals $\eta_0$ to good approximation (exactly for a constant-density fluid).

At each iteration $k$:

1. Construct the p-star coordinate from $P_\text{bot}^{(k)}$, obtaining
   `RefPseudoThickness`, `PseudoThickness`, `ZTildeMid`, `ZTildeInterface`, and related
   fields. The coordinate-construction step may snap `BottomPressure` to a layer boundary
   for full or partial bottom cells; denote the snapped value $\hat{P}_\text{bot}^{(k)}$.
2. Obtain CT and SA at the p-star layer midpoints through the CT/SA initialization
   interface (see the next algorithm design section).
3. Compute sea pressure at each layer midpoint from pseudo-height:
   $p_\text{mid} = -\rho_0 g \tilde{z}_\text{mid}$.
4. Evaluate the EOS to obtain specific volume:
   $\alpha = \text{EOS}(\text{CT},\, \text{SA},\, p_\text{mid})$.
5. Recover geometric layer thickness from pseudo-thickness and specific volume:
   $\delta z = \alpha\, \rho_0\, \delta\tilde{z}$.
6. Sum the geometric layer thicknesses upward from the seafloor to obtain the geometric
   water-column thickness $H_\text{geo}^{(k)}$.
7. Compute the proportional scaling factor:
   $s^{(k)} = H_\text{geo}^\star / H_\text{geo}^{(k)}$.
8. Check convergence: if $k > 0$ and
   $|H_\text{geo}^{(k)} - H_\text{geo}^{(k-1)}| / H_\text{geo}^{(k-1)} < \epsilon$
   (configurable tolerance), stop.
9. Check for full-cell stagnation: if
   $\hat{P}_\text{bot}^{(k)} = \hat{P}_\text{bot}^{(k-1)}$, stop early and log a warning
   (see the bottom-cell algorithm design section).
10. Update, scaling the **pseudo-column pressure** $\hat{P}_\text{bot}^{(k)} - P_\text{surf}$
    rather than $\hat{P}_\text{bot}^{(k)}$ directly:
    $$P_\text{bot}^{(k+1)} = P_\text{surf} + \bigl(\hat{P}_\text{bot}^{(k)} - P_\text{surf}\bigr)\cdot s^{(k)}$$

After convergence, `bottomDepth` is set to the actual recovered geometric water-column
thickness $H_\text{geo}^{(k)}$ (not the target $H_\text{geo}^\star$) for thermodynamic
self-consistency, even when bottom-cell snapping prevents an exact match to the target
bathymetry.  The diagnostic `ssh` is computed as the top interface height
$\tilde{z}^\text{top}_\text{min-level}$ after applying
`geom_height_from_pseudo_height`, which integrates upward from $z_\text{bot}$.

The iteration evaluates the EOS once per outer iteration (step 4 above). This is distinct
from the inner fixed-point iteration in
`pressure_and_spec_vol_from_state_at_geom_height` in
`polaris/ocean/vertical/ztilde.py`, which recovers pressure from geometric layer
thicknesses when pseudo-height is not known directly. In the present algorithm, pressure
is derived analytically from pseudo-height (step 3), so the inner iteration is not
invoked within the outer loop.

### Algorithm Design: CT and SA can be initialized on the p-star coordinate through a pluggable interface

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The CT/SA initialization interface is realized as an overridable method on a base step
class. The method receives the current p-star dataset — which contains `ZTildeMid`,
`ZTildeInterface`, `PseudoThickness`, `cellMask`, `minLevelCell`, `maxLevelCell`, and
associated coordinate information — and returns CT and SA as arrays with dimensions
`(Time, nCells, nVertLevels)`.

This follows the same generalization pattern used in the convergence task framework,
where concrete subclasses implement or override step methods that are invoked by shared
base-class logic.

Several concrete initialization strategies are envisioned:

1. **Piecewise pseudo-height profiles from config** (current `horiz_press_grad` approach):
   node values of CT and SA are read from the configuration at specified pseudo-height
   levels, and a monotone interpolant (e.g., PCHIP) maps them to each layer midpoint.
2. **Analytic profiles**: CT and SA are computed from a formula as a function of
   pseudo-height, geometric depth, or another column coordinate.
3. **Observational hydrography** (planned for realistic initialization): CT and SA are
   interpolated from a pre-processed hydrography product that has already been remapped
   to the MPAS horizontal mesh; at each outer iteration, vertical interpolation from the
   source depth levels to the current p-star midpoints is performed.

The interface constrains only what the method receives and returns; it does not constrain
how CT and SA are computed. The method may access `self.config`, `self.logger`, or any
other attributes of the step instance.

### Algorithm Design: The initialized state includes all p-star coordinate variables and derived geometric quantities needed for subsequent steps

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The Polaris ocean framework separates model inputs into three files. The horizontal mesh
file (`mesh.nc` or `culled_mesh.nc`) is not produced by the p-star initialization step;
it comes from an upstream mesh-construction or cull step. The p-star initialization step
produces the remaining two files.

**Vertical coordinate file** (`vert_coord.nc`, written by `write_vert_coord_dataset`):
For Omega this file feeds the `InitialVertCoord` stream. For MPAS-Ocean
`write_vert_coord_dataset` is a no-op and these variables remain in the initial-state
file instead. Variable names in parentheses are the Omega-native equivalents after the
framework's variable renaming.

| MPAS-Ocean variable | Omega variable | Description | Key dimensions |
|---|---|---|---|
| `minLevelCell` | `MinLayerCell` | First valid layer index (1-based) | nCells |
| `maxLevelCell` | `MaxLayerCell` | Last valid layer index (1-based) | nCells |
| `bottomDepth` | `BottomGeomDepth` | Actual geometric water-column thickness | nCells |
| `RefPseudoThickness` | `RefPseudoThickness` | Reference pseudo-thickness at zero SurfacePressure (no Time dim) | nCells, nVertLevels |
| `vertCoordMovementWeights` | `VertCoordMovementWeights` | Weights for coordinate movement (all ones) | nVertLevels |

**Initial-state file** (`init.nc`, written by `write_initial_state_dataset`): horizontal
mesh variables are stripped before writing; for Omega, the vertical coordinate variables
above are also stripped (they are in `vert_coord.nc` instead).

| MPAS-Ocean variable | Omega variable | Description | Key dimensions |
|---|---|---|---|
| `temperature` | `Temperature` | Conservative temperature (CT) at convergence | Time, nCells, nVertLevels |
| `salinity` | `Salinity` | Absolute salinity (SA) at convergence | Time, nCells, nVertLevels |
| `normalVelocity` | `NormalVelocity` | Initial velocity (typically zero) | Time, nEdges, nVertLevels |
| `PseudoThickness` | `PseudoThickness` | Pseudo-layer thickness | Time, nCells, nVertLevels |
| `ZTildeMid` | `ZTildeMid` | Pseudo-height at layer midpoints | Time, nCells, nVertLevels |
| `ZTildeInterface` | `ZTildeInterface` | Pseudo-height at layer interfaces | Time, nCells, nVertLevelsP1 |
| `SurfacePressure` | `SurfacePressure` | Surface pressure | nCells |
| `BottomPressure` | `BottomPressure` | Effective (post-snap) seafloor pressure | nCells |
| `cellMask` | `cellMask` | Boolean mask of valid layers | nCells, nVertLevels |
| `pressure` | `PressureMid` | Sea pressure at layer midpoints | Time, nCells, nVertLevels |
| `zMid` / `GeomZMid` | `GeomZMid` | Geometric height at layer midpoints | Time, nCells, nVertLevels |
| `zInterface` / `GeomZInterface` | `GeomZInterface` | Geometric height at layer interfaces | Time, nCells, nVertLevelsP1 |
| `ssh` | `SshCell` | Sea-surface geometric height (diagnostic) | nCells |

Additional task-specific diagnostic variables (e.g., `SpecVol`, `Density`, Montgomery
potential fields) may be included in `init.nc` by concrete subclasses.

### Algorithm Design: Full and partial bottom cells are handled correctly within the iteration

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The p-star coordinate-construction step (`init_pstar_vertical_coord` in
`polaris/ocean/vertical/pstar.py`) handles full and partial bottom cells by
snapping the pseudo-column-thickness to a discrete layer boundary before returning. Within
the outer iteration, the proportional-ratio update (step 10 of the algorithm) is applied
to the post-snap value $\hat{P}_\text{bot}^{(k)}$, not the raw value
$P_\text{bot}^{(k)}$. This ensures that each successive iterate moves relative to the
effective snapped position and prevents the iteration from swinging between two snapping
levels.

For full-cell configurations, `BottomPressure` is forced to the bottom of a reference
layer regardless of the incoming value. Once the iteration converges on the nearest full
cell, $\hat{P}_\text{bot}^{(k)}$ no longer changes between iterations even though the
scaling factor $s^{(k)}$ may differ from 1. The algorithm detects stagnation by comparing
$\hat{P}_\text{bot}^{(k)}$ to $\hat{P}_\text{bot}^{(k-1)}$ and stops early, logging a
warning that reports the residual mismatch between the recovered and target geometric
water-column thicknesses. The residual is bounded by the reference pseudo-layer spacing at
the seafloor.

For partial-cell configurations, `BottomPressure` is allowed to vary continuously between
the partial-cell minimum threshold and the full-cell boundary. Convergence within this
window is generally achievable, and the stagnation check triggers only if the partial-cell
constraints force a discrete jump.

### Algorithm Design: The capability is reusable across idealized and realistic initialization tasks

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The iteration algorithm is encapsulated in a base step class (see the implementation
section) that provides the outer fixed-point loop, output dataset assembly, convergence
checking, full-cell stagnation detection, and iteration logging. Concrete subclasses
supply only the CT/SA initialization method, making it straightforward to add new
initialization strategies without touching the core algorithm.

The `horiz_press_grad.Init` step is the first concrete implementation. The planned
realistic initialization step will be the second. No changes to the base class are
anticipated when adding new CT/SA initialization strategies.

## Implementation

### Implementation: The p-star vertical coordinate can be initialized with a geometric seafloor depth matching target bathymetry within a configurable tolerance

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The outer fixed-point loop is encapsulated in a `run_pstar_init` method on the
`PStarInitStep` base class in `polaris/ocean/vertical/pstar_init.py`. The method takes
`ds_mesh`, `geom_z_bot`, an optional `surface_pressure` (defaulting to zero), and an
optional `sea_surface_height` as arguments. `sea_surface_height` is the prescribed SSH;
it defaults to `-surface_pressure / (RhoSw * g)`. The target water-column thickness is
set to `sea_surface_height - geom_z_bot`; the iteration adjusts `BottomPressure` until
this target is met, at which point the diagnostic `ssh` equals `sea_surface_height`. For
the default `sea_surface_height`, the initial `BottomPressure` simplifies to
`RhoSw * g * (-geom_z_bot)` because the `surface_pressure` terms cancel. The
`horiz_press_grad.Init` class inherits from `PStarInitStep`, implements the `init_tracers`
abstract method, and delegates the entire iteration to `run_pstar_init`.

The convergence threshold and maximum iteration count are read from the `vertical_grid`
configuration section. `pseudothickness_iter_count` was already in that section;
`water_col_adjust_frac_change_threshold` is placed there as well so that it is available
to all users of the base class without task-specific configuration. The default value
(`1e-12`) is set in `polaris/ocean/ocean.cfg`.

### Implementation: CT and SA can be initialized on the p-star coordinate through a pluggable interface

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The CT/SA initialization interface is declared using Python's `abc.abstractmethod`
decorator on the `PStarInitStep` base class. Because the `Step` base class does not
include `ABCMeta`, `PStarInitStep` inherits from both `Step` and `abc.ABC` to
activate abstract-method enforcement. The required method signature is:

```python
@abstractmethod
def init_tracers(
    self, ds: xr.Dataset
) -> tuple[xr.DataArray, xr.DataArray]:
    """
    Initialize CT and SA at p-star layer midpoints for the current iteration.

    Parameters
    ----------
    ds : xarray.Dataset
        Current p-star dataset, including ZTildeMid, ZTildeInterface,
        PseudoThickness, cellMask, minLevelCell, and maxLevelCell.

    Returns
    -------
    conservative_temperature : xarray.DataArray
        CT with dimensions (Time, nCells, nVertLevels).
    absolute_salinity : xarray.DataArray
        SA with dimensions (Time, nCells, nVertLevels).
    """
```

The former `horiz_press_grad.Init._interpolate_t_s` private method is the concrete
`init_tracers` implementation in the refactored `horiz_press_grad.Init` subclass; the
column-position array `x` is stored on `self` before `run_pstar_init` is called so that
`init_tracers` can access it via the step instance.

### Implementation: The initialized state includes all p-star coordinate variables and derived geometric quantities needed for subsequent steps

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The base class `run_pstar_init` method assembles and returns the output dataset from
the quantities accumulated during the outer loop. The concrete subclass `run` method
then writes the split output files using the framework helpers:

```python
self.write_vert_coord_dataset(ds, 'vert_coord.nc', config)
self.write_initial_state_dataset(ds, 'init.nc', config)
```

`write_vert_coord_dataset` is a no-op for MPAS-Ocean; `write_initial_state_dataset`
strips horizontal mesh variables and (for Omega) also strips vertical coordinate variables
before writing. Since `horiz_press_grad` is Omega-only, `vert_coord.nc` is registered
unconditionally as a step output file.

Any model-specific or task-specific fields not part of the p-star base output (e.g.,
`normalVelocity`, pressure-gradient diagnostics) are added to `ds` by the concrete
subclass before calling the write helpers. The horizontal mesh file is written separately
and is not produced by `PStarInitStep`.

### Implementation: Full and partial bottom cells are handled correctly within the iteration

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The stagnation-detection logic (comparing the post-snap `BottomPressure` between
consecutive iterations) is preserved in the base class. The convergence check (step 8 of
the algorithm) is evaluated before the stagnation check (step 9), so a perfect initial
guess exits via convergence rather than triggering a misleading stagnation warning. The
warning message reports the iteration number to aid debugging on meshes where full-cell
behaviour is unexpected.

### Implementation: The capability is reusable across idealized and realistic initialization tasks

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

`horiz_press_grad.Init` is the first concrete subclass of `PStarInitStep`. It uses the
three-file split (`culled_mesh.nc`, `vert_coord.nc`, `init.nc`) and demonstrates that the
base class is consistent with that pattern.

Because `horiz_press_grad` allows `z_tilde_bot` to vary per cell via a configurable
gradient, `horiz_press_grad.Init` also overrides the `_build_pstar_coord_ds(self,
ds_mesh, bottom_pressure, surface_pressure)` method. The base-class default calls
`init_pstar_vertical_coord` once on the full mesh; `horiz_press_grad.Init` overrides
it with a per-cell loop that sets a per-cell reference depth in a local config copy before
constructing each column's coordinate. This override pattern is consistent in style with
the `init_tracers` interface and keeps the per-cell logic in the subclass.

The planned realistic initialization step (`ocean/realistic/init`) will subclass
`PStarInitStep` and implement `init_tracers` to read from the WOA hydrography product.
In that context, CT/SA initialization at each outer iteration involves sampling a
pre-computed hydrography product that has been remapped to the MPAS horizontal mesh, then
interpolating vertically from the source depth levels to the current p-star layer
midpoints. The mesh file for realistic initialization comes from an upstream `e3sm/init`
cull step rather than being constructed by the init step itself, so only `vert_coord.nc`
and `init.nc` are produced by `PStarInitStep` in that context.

## Testing

### Testing and Validation: The p-star vertical coordinate can be initialized with a geometric seafloor depth matching target bathymetry within a configurable tolerance

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

Unit tests in `tests/ocean/vertical/test_pstar_init.py` cover:

- `test_constant_density_converges_cleanly`: single-column case with
  $\rho = \rho_0$ everywhere ($s^{(k)} = 1$ at every iteration). The loop exits via the
  convergence check on iteration 1 with `bottomDepth` matching the target to within
  floating-point precision and no stagnation warning.
- `test_wrong_reference_density_requires_multiple_iterations`: single-column case with
  a 1% density offset from $\rho_0$, requiring 3 outer iterations before the fractional
  change in water-column thickness falls below the threshold. The final `bottomDepth`
  matches the target to within a relative tolerance well below $10^{-10}$.
- `test_nonzero_surface_pressure_shifts_pressures`: single-column case with non-zero
  $P_\text{surf}$ and the default $\eta_0 = -P_\text{surf}/(\rho_0 g)$. Confirms that
  the converged `ssh` equals $-P_\text{surf}/(\rho_0 g)$ and `bottomDepth` equals
  the actual water-column thickness $\eta_0 - z_\text{bot}$, and that
  `BottomPressure > SurfacePressure`.
- `test_surface_pressure_sets_ztilde_surface_correctly`: verifies that
  `ZTildeInterface.isel(nVertLevelsP1=0)` equals $-P_\text{surf}/(\rho_0 g)$, confirming
  the top of the pseudo-height coordinate is correctly anchored at the surface pseudo-height.
- `test_surface_pressure_total_pseudo_thickness`: verifies that the sum of
  `PseudoThickness` over valid layers equals
  $(P_\text{bot} - P_\text{surf})/(\rho_0 g)$, confirming the z-star-style scaling
  $(P_\text{bot} - P_\text{surf})/P_\text{bot}$ is correctly applied.

The `horiz_press_grad` regression tests (all three task variants: `salinity_gradient`,
`temperature_gradient`, `ztilde_gradient`) validate that the refactored base class
produces the same outputs as the original embedded loop. All three variants passed the
`omega_pr` suite on Chrysalis (2026-06-01), with `salinity_gradient` producing results
identical to the `main`-branch baseline.

### Testing and Validation: CT and SA can be initialized on the p-star coordinate through a pluggable interface

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

Unit tests in `tests/ocean/vertical/test_pstar_init.py` cover:

- `test_abstract_class_not_instantiable`: confirms `init_tracers` is in
  `PStarInitStep.__abstractmethods__`, verifying that the abstract-method mechanism
  is active.
- `test_subclass_without_init_tracers_not_instantiable`: confirms that a subclass that
  does not override `init_tracers` also carries the abstract method.
- `test_minimal_subclass_returns_complete_dataset`: a minimal concrete subclass with a
  constant `init_tracers` completes the outer loop without error and the returned dataset
  contains all 17 required base-class output variables with correct dimensions.

The surface-pressure tests (`test_nonzero_surface_pressure_shifts_pressures`,
`test_surface_pressure_sets_ztilde_surface_correctly`,
`test_surface_pressure_total_pseudo_thickness`) also exercise `init_tracers` with a
shifted `ZTildeMid` range, confirming that the coordinate geometry produced under
non-zero $P_\text{surf}$ is compatible with the pluggable interface.

### Testing and Validation: The initialized state includes all p-star coordinate variables and derived geometric quantities needed for subsequent steps

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

The `test_minimal_subclass_returns_complete_dataset` unit test verifies that the output
dataset contains all variables listed in the algorithm design section tables. Regression
tests for `horiz_press_grad` will verify correct dimensions and non-NaN values in valid
cells across the three task variants.

### Testing and Validation: Full and partial bottom cells are handled correctly within the iteration

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

Unit test `test_full_cell_stagnation_warning_and_early_exit` in
`tests/ocean/vertical/test_pstar_init.py` overrides `_build_pstar_coord_ds` to always
return a fixed `BottomPressure` regardless of input, simulating full-cell snapping. The
test confirms that a warning containing `'full-cell snap'` is emitted and that the loop
terminates without raising an exception, and that the returned dataset still contains
`bottomDepth`.

The `horiz_press_grad` regression tests (all three task variants) passed the `omega_pr`
suite on Chrysalis (2026-06-01). The `horiz_press_grad` tasks all use
`partial_cell_type = partial`, exercising the partial-cell path through the iteration.

### Testing and Validation: The capability is reusable across idealized and realistic initialization tasks

Date last modified: 2026/06/13

Contributors: Xylar Asay-Davis, Claude

`horiz_press_grad.Init` is refactored onto `PStarInitStep` and its existing regression
tests serve as the first validation that the base class interface is correct and complete.
All three `horiz_press_grad` task variants passed the `omega_pr` suite on Chrysalis
(2026-06-01), confirming that the outer loop converges and the output dataset is complete.
Once the realistic initialization task (`ocean/realistic/init`) is implemented, its
regression tests will exercise the shared base class with WOA-derived CT/SA on at least
one small or moderate-resolution global mesh, confirming that no modification to the base
class is needed for realistic initialization.