# Realistic Ocean Initial Conditions

Creation date: 2026/03/22

Contributors: Xylar Asay-Davis, Codex, Claude

## Summary

This design document describes a new Polaris capability for generating realistic
ocean initial conditions in Python, starting from an ocean mesh and topography
that have already been culled by the `e3sm/init` workflow. This capability
belongs to the `ocean/spherical/realistic_global` framework, which collects task families focused
on observationally-constrained ocean initialization and evaluation, in contrast
to the idealized task families elsewhere in the ocean component. The intent is to
port the required functionality from the MPAS-Ocean init mode, especially
`mpas_ocn_init_global_ocean.F`, and from the Compass workflows in
`compass/ocean/tests/global_ocean/init` and
`compass/ocean/tests/utility/extrap_woa`.

Each task initializes the single ocean model selected by the `[ocean] model`
config option (resolved to `omega` or `mpas-ocean` during component setup).
The new workflow produces scientifically consistent initial conditions for the
configured model. In the first phase, the same capability supports MPAS-Ocean
tracers based on potential temperature and practical salinity and Omega tracers
based on conservative temperature and absolute salinity. The initial phase does
not need to support ice-shelf cavities, Haney-number vertical coordinates under
ice shelves, or SSH adjustment below ice shelves, but the design should preserve
a clean path to add those features later.

A key challenge specific to the Omega p-star coordinate is that the
geometric seafloor depth cannot be prescribed directly: `BottomPressure` (related
to the pseudo-height of the seafloor) must be found iteratively so that the geometric
depth recovered from the coordinate matches observed bathymetry. Because CT and
SA must be initialized on the p-star grid before specific volume can be
evaluated, the vertical coordinate and the tracer state are tightly coupled. The
algorithm for this joint initialization is described in the companion design
document [pstar_init.md](pstar_init.md). Its existence means the vertical
coordinate and tracer initial state cannot be produced in isolation from each
other; the decomposition into Polaris steps must reflect this coupling.

Because the p-star iteration is model-independent — the converged geometric
vertical grid depends only on bathymetry and hydrography, not on which ocean
model is the final target — running the task for different models on the same
mesh and hydrography produces the same geometric layer thicknesses. This makes
direct model comparison possible by running the task once per model: the
geometric grids at initialization will be identical, and differences in
simulated state can be attributed to model formulation rather than to different
initial geometries. The p-star iteration runs even when the configured model is
MPAS-Ocean; there is no separate z-star path.

The Polaris ocean framework separates model inputs across staged files: a
horizontal mesh file (`mesh.nc`), a vertical coordinate file (`vert_coord.nc`;
Omega only), an initial-state file (`init.nc`) and — added after this document
was first written — a surface-forcing file (`forcing.nc`), described under
[Wind forcing](#wind-forcing). The realistic init workflow writes all of them,
through the framework's `write_horiz_mesh_dataset`, `write_vert_coord_dataset`,
`write_initial_state_dataset` and `write_forcing_dataset` helpers. The
horizontal geometry is entirely the upstream `e3sm/init` cull workflow's;
`mesh.nc` is that culled mesh with the fields the model's mesh stream needs but
the cull step does not produce. For MPAS-Ocean, `write_vert_coord_dataset` is a
no-op and the vertical coordinate variables remain in `init.nc`.

The primary software challenges are to replace a large, monolithic Fortran
workflow with an inspectable sequence of Polaris steps; to use a model-neutral
p-star intermediate step that produces the same geometric grid regardless of the
target model; and to keep the workflow performant on very large meshes through
chunked xarray/dask operations and materialized intermediate files. This design
is successful if it produces reusable, well-tested Polaris steps that scale to
production global meshes, initialize the configured model from a model-neutral
geometric vertical grid, and can eventually replace the corresponding legacy
Compass and MPAS-Ocean init workflows for the supported open-ocean use case.

## Requirements

### Requirement: Open-ocean global initial conditions can be created from a culled mesh

Date last modified: 2026/03/22

Contributors: Xylar Asay-Davis, Codex

Polaris shall provide a workflow that begins from an ocean mesh and
corresponding culled topography produced by the `e3sm/init` cull workflow and
produces a global ocean initial condition suitable for use by a standalone
ocean model.

The initial supported use case is open ocean without ice-shelf cavities. The
workflow shall therefore assume that the culled mesh does not require cavity
geometry, land-ice pressure initialization, or SSH adjustment beneath floating
ice.

Any topography smoothing, minimum-depth logic, or related preparation of the
ocean topography needed for initialization shall be handled upstream in the
`e3sm/init` topography workflow rather than inside the ocean-initialization
workflow. This supports the goal of a unified MPAS mesh and a common
topographic description across E3SM components.

### Requirement: A reusable global hydrography product is available from WOA

Date last modified: 2026/04/14

Contributors: Xylar Asay-Davis, Codex

Polaris shall provide a reusable hydrography product derived from the World
Ocean Atlas at 0.25-degree resolution. This product shall contain canonical
conservative temperature and absolute salinity fields needed to initialize the
ocean state over the full globe, including values in regions where the original
product is missing and extrapolation is required for later remapping. Because
WOA23 supplies in-situ temperature and practical salinity rather than
conservative temperature and absolute salinity, the workflow shall derive the
canonical tracers during preprocessing.

The hydrography product shall be suitable for caching because it is expected to
be computationally expensive and because it should be reusable across many
meshes and tasks.

### Requirement: The workflow produces consistent tracer initial conditions for the configured model

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Codex, Claude

Given a common hydrographic source, Polaris shall produce initial conditions
that are thermodynamically consistent with the tracer conventions of the
configured ocean model.

For MPAS-Ocean, the workflow shall produce temperature and salinity fields in
the form expected by MPAS-Ocean (potential temperature and practical salinity).
For Omega, the workflow shall produce temperature and salinity fields in the
form expected by Omega (conservative temperature and absolute salinity). The
design shall minimize duplicated logic between the two model targets by sharing
the p-star iteration and WOA preprocessing and branching only where tracer
conversion genuinely differs.

### Requirement: The p-star iteration yields a model-independent geometric vertical grid

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Codex, Claude

The p-star initialization step shall produce a geometric vertical grid that
depends only on bathymetry and hydrography, not on which ocean model is the
final target. As a consequence, running the task for MPAS-Ocean and for Omega
on the same horizontal mesh and hydrography shall produce the same initialized
geometric layer thicknesses, enabling direct geometric comparison between the
two models.

This shared-grid property is a natural outcome of the p-star algorithm: the
converged geometric layer thicknesses emerge from a model-neutral iteration and
are then consumed by the model-specific `initial_state` step. The p-star
iteration runs even when the configured model is MPAS-Ocean; there is no
separate z-star resting-thickness path.

The workflow shall compute any density, pressure, or specific-volume fields
needed to translate hydrography defined with respect to geometric depth or
height into the model's native vertical representation.

Haney-number coordinates are out of scope for the initial implementation but
should remain feasible to add later.

### Requirement: The capability is decomposed into inspectable Polaris steps

Date last modified: 2026/03/22

Contributors: Xylar Asay-Davis, Codex

The workflow shall be organized into multiple Polaris steps with clear,
meaningful outputs that users and developers can inspect independently for
sanity checking, visualization, and debugging.

Where appropriate, steps shall be reusable or cacheable so that expensive work
does not need to be repeated across tasks that use identical inputs.

### Requirement: The workflow is practical for very large global meshes

Date last modified: 2026/03/22

Contributors: Xylar Asay-Davis, Codex

The workflow shall be implementable in a way that is practical for large
production meshes, including meshes on the order of uniform 5-km global
resolution.

The implementation shall avoid assumptions that require the full global state to
fit comfortably in memory on a single node. It shall support sensible
chunk-based execution and the use of intermediate files where needed to keep
the workflow robust.

### Requirement: The design preserves a path to future cavity support

Date last modified: 2026/03/22

Contributors: Xylar Asay-Davis, Codex

Although the first implementation does not need to support ice-shelf cavities,
the workflow shall not hard-code assumptions that make later support for
cavities, land-ice pressure consistency, or SSH adjustment unnecessarily
difficult.

## Algorithm Design

### Algorithm Design: Open-ocean global initial conditions can be created from a culled mesh

Date last modified: 2026/03/22

Contributors: Xylar Asay-Davis, Codex

The workflow should adopt the culled ocean mesh and associated topography from
the `e3sm/init` cull workflow as its starting point rather than reproducing the
mesh-culling logic inside the ocean task. This keeps responsibilities separated:
`e3sm/init` produces physically consistent meshes and masks, while the new ocean
initial-condition workflow consumes those products to build the 3D ocean state.

For the initial phase, the relevant topographic fields should be restricted to
open-ocean quantities needed for bathymetry, land masks, and any mesh-derived
auxiliary fields required by the target model.

This design also assumes that any topography smoothing or minimum-depth
enforcement needed for ocean initialization has already been applied in the
`e3sm/init` topography workflow, so the ocean task treats the incoming
topography as authoritative.

### Algorithm Design: A reusable global hydrography product is available from WOA

Date last modified: 2026/04/14

Contributors: Xylar Asay-Davis, Codex

The workflow should create a canonical hydrography product on the native WOA23
latitude-longitude grid before any remapping to an MPAS mesh. The current
implementation does this with a concrete sequence:

1. Create combined topography on the native 0.25-degree latitude-longitude
   grid.
2. Read WOA23 January and annual temperature and salinity climatologies and
   combine them into a single source product, using annual values at depths
   where the monthly fields are not available, then derive conservative
   temperature and absolute salinity from WOA23 in-situ temperature and
   practical salinity.
3. Build a 3D ocean mask on the WOA grid from the combined topography product.
4. Extrapolate horizontally and then vertically within the ocean mask.
5. Extrapolate horizontally and then vertically into land and grounded-ice
   regions so later remapping will not sample missing values.
6. Optionally produce visualization products for sanity checking.

This is conceptually similar to the existing Compass `extrap_woa` utility but
the Polaris implementation should favor xarray-based operations and clear
intermediate products over multiprocessing code that is tightly coupled to the
legacy workflow.

### Algorithm Design: The workflow produces consistent tracer initial conditions for the configured model

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Codex, Claude

A common intermediate tracer representation should be used to avoid maintaining
two largely separate initialization pipelines. The leading candidate is a
canonical hydrographic state based on conservative temperature and absolute
salinity because those align naturally with Omega and with the desired WOA-based
source product.

Model-specific tracer fields are derived late in the workflow, in the
`initial_state` step:

1. Interpolate the canonical state onto the target mesh columns.
2. Convert to the tracer convention required by the configured model.
3. Write the resulting fields using the variable names and metadata expected by
   the model.

The exact thermodynamic conversions should be finalized during implementation,
but the conversion logic should be isolated so it can be tested independently of
horizontal and vertical interpolation.

### Algorithm Design: The p-star iteration yields a model-independent geometric vertical grid

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Codex, Claude

The model-neutral geometric grid is determined by the `pstar_init` step.
The high-level outcome is:

1. Horizontally remap the source hydrography to each MPAS column on a source
   vertical grid.
2. For each water column, run the coupled p-star and tracer initialization to
   convergence (see [pstar_init.md](pstar_init.md)), producing the p-star
   coordinate, CT and SA at layer midpoints, specific volume, and the converged
   geometric layer thicknesses.
3. Use the resulting initialized geometric layer thicknesses in the configured
   model: as part of the Omega initial condition, or as `restingThickness` in
   MPAS-Ocean.
4. The `initial_state` step derives any remaining model-specific vertical state,
   including ALE layer-thickness fields for MPAS-Ocean.

Because the Omega p-star coordinate is defined in pseudo-height rather than
geometric height, the geometric seafloor depth cannot be prescribed directly.
The pseudo-height of the seafloor (`BottomPressure`) must be found iteratively
so that the geometric depth recovered from the coordinate matches observed
bathymetry. This requires initializing CT and SA on the p-star grid at each
iteration step so that specific volume can be evaluated, making the vertical
coordinate and the tracer state tightly coupled. Full and partial bottom cells
add further complexity by introducing discrete jumps in the pseudo-bottom depth.

The detailed fixed-point algorithm for this joint p-star and tracer
initialization — including the proportional-ratio update for `BottomPressure`,
the CT/SA initialization interface, the handling of full and partial bottom
cells, and the complete list of output variables produced — is specified in the
companion design document [pstar_init.md](pstar_init.md).

The initialized geometric grid emerges from the coupled iteration rather than
from a geometric grid constructed in isolation; the workflow steps must reflect
this coupling.

### Algorithm Design: The workflow is practical for very large global meshes

Date last modified: 2026/03/22

Contributors: Xylar Asay-Davis, Codex

The workflow should treat chunking and materialization as part of the algorithm
design rather than as a later optimization. Likely design principles are:

1. Chunk horizontally over cells and only keep the full vertical column when
   needed.
2. Use explicit intermediate files between major phases so dask graphs stay
   bounded in size.
3. Reuse interpolation weights or cached preprocessed products where practical.
4. Keep expensive EOS calculations localized to the phases that need them.

## Implementation

### Implementation: Open-ocean global initial conditions can be created from a culled mesh

Date last modified: 2026/03/22

Contributors: Xylar Asay-Davis, Codex

The new task should live under the ocean framework and consume outputs from a
culled-mesh step in `e3sm/init`, through explicit step inputs rather
than implicit filesystem assumptions. The initial task variant should target
open-ocean meshes without cavities and should accept a culled mesh,
corresponding graph file, and culled topography as inputs.

The first draft of the workflow should port only the parts of the current
Compass and init-mode logic that are required for this open-ocean use case.
Features related only to cavities, SSH adjustment, or ecosystem forcing can be
left out until a later requirement is added for them.

This implementation should not perform its own topography smoothing or
minimum-depth preparation. Those operations should be implemented in
`e3sm/init` so the same prepared topography can be shared consistently across
ocean, sea-ice, land, and river workflows.

### Implementation: A reusable global hydrography product is available from WOA

Date last modified: 2026/08/11

Contributors: Xylar Asay-Davis, Codex, Claude

This part of the design has landed. `e3sm/init` creates combined topography on
a native 0.25-degree latitude-longitude grid, and the WOA23 steps under
`ocean/spherical/realistic_global/hydrography/woa23` build the corresponding
reusable product on that grid:

1. `combine_topo`, the shared `e3sm/init` `CombineStep` for the 0.25-degree
   lat-lon grid, cached by default so the expensive blending is not repeated
2. `combine`, which merges January and annual WOA23 climatologies into
   `woa_combined.nc`, using annual values at depths where the monthly fields
   are not available, and derives conservative temperature and absolute
   salinity from WOA23's in-situ temperature and practical salinity with `gsw`
3. `extrapolate`, which builds a 3D ocean mask from the combined topography and
   produces the reusable product `woa23_decav_0.25_jan_extrap.nc`
4. `viz`, a diagnostics step for maps and Antarctic transects, created as a
   shared step but not run by default

The steps are assembled by `get_woa23_steps(component, include_viz)`, which
returns them keyed by symlink name along with their shared config. A consumer
that needs only the hydrography product — the mesh-specific init steps, for
instance — gets the same step instances rather than a second copy of the
preprocessing.

The rename worked out differently from what this section anticipated. The work
directories did move under the `realistic_global` family, but the Python
package did not gain a `spherical` level: the code is
`polaris.tasks.ocean.realistic_global.hydrography.woa23`. Work-directory paths
and module paths are separate concerns in Polaris and only the first is what a
user navigates, so there was no reason to make the module path longer. The same
split holds throughout this workflow: mesh-specific steps run under
`ocean/spherical/realistic_global/{mesh_name}/init` and live in
`polaris.tasks.ocean.realistic_global.init`.

### Implementation: The workflow produces consistent tracer initial conditions for the configured model

Date last modified: 2026/08/11

Contributors: Xylar Asay-Davis, Codex, Claude

The workflow defines a model-neutral intermediate dataset (the `pstar_init`
step output) and leaves the model-specific translation to the write. That
separation is what lets the p-star iteration, WOA23 remapping and chunking
logic be shared regardless of the configured target model.

The translator layer ended up in the framework rather than in `initial_state`:
`OceanIOStep.write_initial_state_dataset` takes the convention the dataset's
tracers are in and converts them to the one the configured model expects,
defaulting to the convention implied by `eos_type`. The task step supplies only
what the framework cannot know — the per-cell longitude and latitude the
conversion needs, which `pstar_init.nc` does not carry. Putting it there means
every task that writes an initial state converts tracers the same way, and a
step cannot forget to.

Thermodynamic calculations that are independent of the task live in
`polaris.ocean.eos` and `polaris.ocean.init_state`, as this section intended.

### Implementation: The p-star iteration yields a model-independent geometric vertical grid

Date last modified: 2026/08/11

Contributors: Xylar Asay-Davis, Codex, Claude

The coupled p-star and tracer initialization described in
[pstar_init.md](pstar_init.md) shall be implemented as a concrete subclass of
`PStarInitStep` (defined in `polaris/ocean/vertical/pstar_init.py`). That
subclass, `RealisticPStarInitStep`, is defined in
`polaris/tasks/ocean/realistic_global/init/pstar_init.py`.
It implements the `init_tracers` method by interpolating CT and SA from the
pre-processed WOA hydrography product remapped to the MPAS horizontal mesh,
using the p-star layer midpoints from the current iteration as the vertical
interpolation target.

The reference 1D vertical grid that establishes the initial layer structure is
configured by setting `grid_type` in the `[vertical_grid]` section of
`realistic_global_init.cfg`. The `pstar_init` step passes this reference grid
to `generate_1d_grid` (in `polaris/ocean/vertical/grid_1d/__init__.py`),
which returns the interface depths used to dimension the p-star iteration. The
default for the realistic global init task is the pre-defined `80layerE3SMv1`
grid. The `grid_type` option must be present in the config; if it is absent,
`generate_1d_grid` raises `ValueError: Unexpected grid type: None`.

The vertical grid turned out to be one of the things that has to vary from mesh
to mesh, which this document did not anticipate. It is set per mesh through
`polaris.tasks.ocean.realistic_global.mesh_configs`: an optional
`<mesh_name>.cfg` whose options are added after the task's own config file and
therefore override it. The three 240 km meshes replace the 80-layer default
with a 16-level `tanh_dz` grid over a 3000 m bottom depth, because they exist
for fast smoke-testing rather than for realistic simulation. The same mechanism
carries the ocean-culled cell count used to size MPI tasks, and is where
per-mesh forward-run options belong as well. Anything that describes what the
*ocean* does on a mesh belongs there rather than in the mesh component's own
per-mesh config, which describes the mesh itself.

Three behaviors of the iteration were settled during implementation rather than
in [pstar_init.md](pstar_init.md), which specifies the fixed-point algorithm
itself. The column is anchored at the prescribed sea surface, so `ssh` matches
its prescribed value to machine precision and any residual left by partial-cell
snapping adjusts the diagnosed `bottomDepth` — the representable bathymetry —
rather than `ssh`, mirroring z-star partial cells. Columns too shallow to hold
`min_vert_levels` layers, or shallower than `min_bottom_depth`, are clamped.
And isolated bathymetry holes — cells deeper than every ocean neighbor — are
capped at their deepest neighbor's level and re-solved, through
`polaris.ocean.vertical.bathymetry_holes.fill_max_level_holes`.

The output of the coupled initialization step — converged geometric layer
thicknesses, CT and SA, specific volume, and associated coordinate fields — shall
be exposed as an inspectable intermediate product (`pstar_init`). That
intermediate is a single combined dataset written in neutral (MPAS-Ocean) naming
and includes at least the full set of variables listed in the output table of
[pstar_init.md](pstar_init.md). The model-specific split into separate output
files happens in the downstream `initial_state` step rather than in
`pstar_init`, keeping the intermediate inspectable without model-specific
naming or variable filtering.

The `initial_state` step consumes the `pstar_init` intermediate dataset,
performs any model-specific tracer conversions, and writes the final split
output:

- `vert_coord.nc` via `write_vert_coord_dataset` (Omega only; no-op for
  MPAS-Ocean). For Omega this converts `restingThickness` to `RefPseudoThickness`
  and writes the five `InitialVertCoord` fields (`MinLayerCell`, `MaxLayerCell`,
  `BottomGeomDepth`, `RefPseudoThickness`, `VertCoordMovementWeights`).
- `init.nc` via `write_initial_state_dataset`, which strips horizontal mesh
  variables and (for Omega) vertical coordinate variables before writing.

Surface forcing is *not* written by `initial_state`. Wind stress is produced by
a separate `forcing` step and written to its own `forcing.nc`, described under
[Wind forcing](#wind-forcing) below. Restoring fields remain future work.

The mesh file (`mesh.nc`) is written by `initial_state` rather than linked from
upstream. The culled mesh from `e3sm/init` is not by itself a complete ocean
mesh stream: `add_coriolis_to_dataset` adds `fCell`, `fEdge` and `fVertex`,
and for Omega `write_horiz_mesh_dataset` merges in the cell-centered
vector-reconstruction weights, which have to be the *culled* mesh's since that
is the mesh the initial condition is built on. The horizontal geometry itself
is still entirely the cull step's; nothing here re-derives it.

The p-star iteration shall run even when the configured model is MPAS-Ocean,
since the converged geometric grid is always defined through the p-star
iteration.

The MPAS-Ocean path derives `restingThickness` from the converged geometric
layer thicknesses rather than falling back to the legacy z-star resting-thickness
construction.

### Implementation: The capability is decomposed into inspectable Polaris steps

Date last modified: 2026/09/24

Contributors: Xylar Asay-Davis, Codex, Claude

The decomposition landed with more steps than this section first sketched,
because remapping turned out to be worth separating from the mapping files it
uses, and because wind forcing arrived as a second remapped source. The chain
is assembled by `get_realistic_init_steps(component, mesh_name, include_viz)`,
which builds every step with `Component.get_or_create_shared_step` and returns
them keyed by symlink name along with the per-mesh shared config. Steps run
under `ocean/spherical/realistic_global/{mesh_name}/init`, except the
mesh-independent preprocessing, which is shared across meshes:

1. the upstream `e3sm/init` cull-topography steps, from
   `get_cull_topo_steps`, which supply the culled ocean mesh, its graph and
   the remapped topography
2. the WOA23 hydrography steps described above, from `get_woa23_steps`
3. the JRA55-do wind-stress steps, from `get_jra55_steps` (see
   [Wind forcing](#wind-forcing))
4. `cull_topo`: reindex the remapped topography from the base mesh to the
   culled ocean mesh using `ocean_map_culled_to_base.nc`, producing
   `topography_culled.nc`, which is validated against a baseline when one is
   provided
5. `woa23_map`: build the bilinear mapping file from the 0.25-degree WOA23
   lat-lon grid to the culled MPAS mesh. This is the one MPI step of the
   hydrography chain, and its task count is sized from the approximate culled
   ocean cell count
6. `remap_woa23`: apply those weights with `ncremap`, producing
   `woa23_on_mesh.nc`
7. `jra55_map` and `remap_jra55`: the same two steps for the wind stress
8. `pstar_init`: perform the coupled p-star and hydrography initialization
   described in [pstar_init.md](pstar_init.md), producing the converged
   geometric grid, CT and SA on that grid, and all associated p-star
   coordinate fields; this step writes a single combined intermediate NetCDF
   file (`pstar_init.nc`) in neutral naming rather than the final split model
   files, so all outputs remain inspectable regardless of the target model
9. `initial_state`: consume the `pstar_init` intermediate dataset, apply
   model-specific tracer conversions, populate remaining required fields, and
   write `mesh.nc`, `vert_coord.nc` (Omega only) and `init.nc`
10. `forcing`: write the model-specific `forcing.nc` from `jra55_on_mesh.nc`
11. `viz`: plots and sanity checks of the initial condition, vertical
    coordinate and forcing

Each remap is two steps rather than one because the original single step both
reimplemented what `polaris.remap.MappingFileStep` already provides and ran an
MPI operation (`mbtempest`) and a serial one (`ncremap`) together, which does
not schedule well under task parallelism. Splitting them also lets the mapping
step size its own task count from the estimated ocean cell count while the
`ncremap` call stays serial.

The step this section called `diagnostics` landed as `viz`, matching the name
used elsewhere in Polaris. It is created as a shared step but returned — and so
run — only when `include_viz=True`, which the standalone `RealisticGlobalInit`
task passes and consumers that reuse the init outputs as dependencies do not.
Workflows that build on the initial condition, such as forward runs and
dynamic adjustment, are meant to call `get_realistic_init_steps` the same way
and get the same step instances rather than a second copy of a chain that
costs hours.

The `viz` step makes Omega's geometric `layerThickness` itself. Omega's initial
state carries `PseudoThickness` and no `SpecVol`, and `open_model_dataset` no
longer derives equation-of-state fields (#817), so it only turns a
pseudo-thickness into a geometric one for a file that has both. The step
computes the specific volume with `spec_vol_from_ds()` instead.

Within this decomposition, `initial_state` is what ports the pieces of legacy
init-mode functionality that are neither part of `pstar_init` nor already
handled upstream in `e3sm/init`. As implemented, it:

1. Consumes the outputs of `pstar_init`, including the converged geometric
   layer thicknesses and the model-agnostic hydrographic state on the target
   mesh.
2. Converts that hydrographic state into the tracer conventions required by the
   configured model — CT and SA kept for Omega, converted to potential
   temperature and practical salinity for MPAS-Ocean. The conversion itself
   belongs to the framework's initial-state helpers; this step supplies the
   per-cell longitude and latitude they need, since `pstar_init.nc` carries no
   horizontal mesh fields.
3. Populates quiescent dynamical initial conditions, setting `normalVelocity`
   to zero.
4. Computes the remaining derived fields, notably in-situ density from the
   specific volume the p-star iteration produced.
5. Writes `mesh.nc` through `write_horiz_mesh_dataset`, `vert_coord.nc`
   through `write_vert_coord_dataset` (Omega only; a no-op for MPAS-Ocean,
   which keeps vertical coordinate variables in `init.nc`) and `init.nc`
   through `write_initial_state_dataset`.

Restoring fields did not land. Surface and interior restoring values and their
piston velocities are not part of the supported workflow yet, so nothing
computes them; the forcing file that replaces the legacy
`init_mode_forcing_data.nc` carries wind stress alone and is written by the
separate `forcing` step rather than by `initial_state`.

The first implementation should explicitly exclude or defer several categories
of legacy init-mode behavior that are present in
`mpas_ocn_init_global_ocean.F` but are not part of the current design scope:

1. Ice-shelf and land-ice-specific logic, including cavity temperature
   modification, land-ice pressure consistency, and SSH adjustment.
2. Ecosystem tracers, ecosystem forcing fields, and related pH or sediment
   initialization.
3. Shortwave-absorption auxiliary fields needed for the `ohlmann00`
   parameterization, such as chlorophyll, zenith angle, and clear-sky
   radiation.
4. Debug tracers and other optional tracer packages initialized only to simple
   default values in legacy init mode.
5. Any inland-sea culling or other mesh/topography correction that is more
   appropriately handled upstream in `e3sm/init`.

The legacy diagnostics this section left open were decided as follows.
Haney-number fields are not written at all — the `viz` step's port of Compass'
`plot_initial_state` drops those panels rather than reproducing them — since
Haney-number coordinates are out of scope. In-situ density is written, because
the p-star iteration already produces the specific volume it comes from and
because it is what the Omega stratification check reads.

(wind-forcing)=
### Implementation: Wind forcing

Date last modified: 2026/08/11

Contributors: Xylar Asay-Davis, Claude

Dynamic adjustment needs a realistic, constant-in-time wind stress so the ocean
can spin down fast waves against a sensible circulation. The forcing does not
need to be a defensible climate product, only realistic enough that the
adjusted state is dynamically sensible. Surface restoring, thermal and
freshwater fluxes, atmospheric and sea-ice pressure, time-varying forcing, and
all forward-model settings remain out of scope.

**Source.** JRA55-do v1.5.0 (`MRI-JRA55-do-1-5-0`), variables `uas` and `vas`:
10-m winds on the TL319 Gaussian grid (640 x 320), 3-hourly, distributed as
input4MIPs on ESGF and anonymously downloadable from the CEDA node. E3SM forces
G-cases with JRA55, so the adjusted state is adjusted against something close to
what a coupled run will apply, and the input4MIPs record pins version and
checksums. This supersedes the earlier intent to reuse the Compass NCEP
climatology, whose provenance is a 2015 filename and nothing more. The time
window is January 1958, the first month of the record and of an interannually
forced G-case, configurable via `year` and `month`.

**Winds to stress.** JRA55-do supplies winds, so a bulk formula is needed. Two
choices matter. First, the stress is averaged rather than the winds: computing
`tau(mean wind)` discards the gust contribution and underestimates stress in the
storm tracks, so stress is computed at each 3-hourly step and then averaged --
which is why the 3-hourly data is needed rather than a published monthly mean.
Second, the Large and Yeager (2004, 2009) neutral 10-m drag law is used with
`rho_air = 1.22` and the wind speed clamped below at 0.5 m/s, with the stability
correction, the current-relative wind, and sea-ice drag deliberately omitted as
second-order for this purpose.

**Remapping.** Bilinear, in Cartesian components, with no padding of the
source grid, and a check that every cell is covered. Each of these was settled
by measurement or review, and each has a failure mode that is easy to
reintroduce:

- Bilinear rather than conservative, because the ocean responds to wind stress
  *curl*; first-order conservative gives a piecewise-constant stress whose curl
  is grid-scale noise, and pyremap's moab path hard-codes `--order 1`.
- Cartesian rather than zonal/meridional components, because the local
  east/north basis rotates with longitude, fastest near the poles. The stress
  is converted to its three Cartesian components on the source grid, each is
  remapped as a scalar, and the result is rotated back to east/north at each
  cell center, as E3SM's coupler does for vectors. Remapping zonal and
  meridional components directly worked under mbtempest (a uniform vector of
  magnitude 1.0 measured a flat 0.998 near the North Pole) but not under ESMF,
  whose default pole handling builds its pole point from the zonal average of
  the source's outermost row: the same vector measured 0.42 by 89.82 degrees N,
  tending to zero at the pole. In Cartesian components that average is
  meaningful, so the result no longer depends on `map_tool`.
- No padding, because bilinear is center-based for ESMF but corner-based for
  mbtempest, and padding a lat-lon source so either its corners or its centers
  reach the pole aborts mbtempest. Duplicating a longitude column breaks both
  tools and is unnecessary, since TL319's longitude corners already span exactly
  360 degrees.
- Every cell of the culled mesh must be covered, and the mapping files for
  `u-oi30-lr10` and `u-oi6to18-lr6to10` cover all of them, the cells nearest
  the North Pole included. `ncremap` writes zero rather than a missing value
  to a cell the weights do not reach, so a gap would otherwise reach the
  forcing file as a silent zero stress. The step instead reads each cell's
  coverage (`frac_b`) from the mapping file and fails if any cell is short of
  it.

**Steps.** Deriving the global product is mesh-independent, and although the
reduction itself takes seconds it needs several GiB of raw winds, so the
product is cached and the work lives outside the `init` workflow, mirroring the
WOA23 hydrography product:

| subdir | step | output |
| --- | --- | --- |
| `forcing/jra55/stress` | `Jra55StressStep` | `jra55_stress.nc` |
| `forcing/jra55/viz` | `Jra55VizStep` | diagnostic plots |
| `{mesh}/init/jra55_map` | `Jra55MapStep` | bilinear weights |
| `{mesh}/init/remap_jra55` | `RemapJra55Step` | `jra55_on_mesh.nc` |
| `{mesh}/init/forcing` | `ForcingStep` | `forcing.nc` |

**Output.** Both models take zonal and meridional components at cell centers and
project onto edges themselves, so no vector rotation is needed. The field is
authored in MPAS-Ocean names and renamed on write through the existing
`mpaso_to_omega.yaml` mapping. The framework gains a fourth staged file
alongside the mesh, vertical coordinate and initial state: `forcing_filename` in
`[ocean_staged_files]`, a `forcing_variables` list, and a
`write_forcing_dataset` helper.

The two models differ on the time dimension. Omega's `SfcStressForcingVars`
registers 1-D fields on `NCells`, while MPAS-Ocean's Registry declares
`dimensions="nCells Time"`, so `write_forcing_dataset` adds `Time=1` for
MPAS-Ocean only.

`Jra55StressStep` is `default_cached`, with `jra55_stress.nc` in the cache
database, so the multi-GiB download happens only when the standalone `jra55`
task regenerates the product.

The read side -- staging `forcing.nc` as a model *input*, and the associated
namelist and config settings -- belongs to the forward-model work rather than
to this design.

### Implementation: The workflow is practical for very large global meshes

Date last modified: 2026/10/08

Contributors: Xylar Asay-Davis, Codex, Claude

The implementation should use xarray lazily by default for large fields and
should choose chunking with the target access pattern in mind. A reasonable
starting assumption is that horizontal chunking over `nCells` will be more
robust than attempting to process the full mesh at once.

Major stages should write intermediate NetCDF outputs rather than carrying a
single dask graph across the entire workflow. This is likely to be especially
important after WOA preprocessing, after `pstar_init`, and before final
model-specific output writing.

Only the second half landed. Every stage writes its own NetCDF output, but the
steps hold whole fields in memory rather than chunking them with dask, which
has been enough for `u-oi6to18-lr6to10` (4.0 million ocean cells) on a single
node. The one place the work is chunked over `nCells` is the vertical
interpolation of the WOA23 tracers in `pstar_init`, which runs at every
fixed-point iteration and handles a chunk of cells at a time rather than one
column at a time. Meshes approaching uniform 5-km resolution will need more
than this.

### Implementation: The design preserves a path to future cavity support

Date last modified: 2026/03/22

Contributors: Xylar Asay-Davis, Codex

Even though cavities are out of scope initially, the data model should not
assume that the free surface always coincides with sea level or that land-ice
pressure fields are absent forever. Interfaces between topography handling,
vertical-coordinate construction, and final state assembly should remain clean
enough that cavity-aware logic can later be inserted without rewriting the full
task.

## Testing

### Testing and Validation: Open-ocean global initial conditions can be created from a culled mesh

Date last modified: 2026/09/24

Contributors: Xylar Asay-Davis, Codex, Claude

The first level of testing should verify that the workflow can consume outputs
from an `e3sm/init` culled-mesh step and produce complete initial-condition
files for at least one supported global mesh. A small or moderate mesh should
be used for routine regression tests, with larger meshes reserved for
integration testing.

That split landed differently than "small mesh in a regression suite" suggests.
The whole chain — culling, remapping and the p-star iteration — costs hours on
any mesh worth initializing, so no task in this family belongs in a PR suite;
routine regression testing is the unit tests described below, and the workflow
itself is exercised by running it on a mesh. The 240 km meshes exist for
exactly that, which is why they override the vertical grid down to 16 levels.
Initial conditions built this way have also been carried through dynamic
adjustment on `u-oi30-lr10` and `u-oi6to18-lr6to10` in downstream work, which
is where artifacts in the WOA23 source data were found.

### Testing and Validation: A reusable global hydrography product is available from WOA

Date last modified: 2026/03/22

Contributors: Xylar Asay-Davis, Codex

Testing should verify that the preprocessed WOA product contains no missing
values in the regions needed for later remapping, that metadata and dimensions
are consistent, and that rerunning the step with the same inputs gives
reproducible results.

### Testing and Validation: The workflow produces consistent tracer initial conditions for the configured model

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Codex, Claude

Unit tests should cover any tracer-conversion utilities independently of the
full task. Task-level validation should verify that model-specific outputs
contain the expected tracer variable names and are thermodynamically consistent
(e.g., potential temperature within physical bounds for MPAS-Ocean; conservative
temperature for Omega).

Where practical, outputs on a small mesh should also be compared against the
legacy Compass plus init-mode workflow, recognizing that exact bit-for-bit
agreement may not be the right success criterion if the new implementation uses
cleaner thermodynamic conversions.

### Testing and Validation: The p-star iteration yields a model-independent geometric vertical grid

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Codex, Claude

Testing should verify that the `pstar_init.nc` intermediate file contains
consistent geometric and pseudo-height fields (layer thicknesses, interface
depths, valid-level masks) that respect bathymetry and minimum-thickness
constraints.

Running the task for both MPAS-Ocean and Omega configurations on the same mesh
and comparing the resulting `pstar_init.nc` intermediates (or the
`restingThickness` / `RefPseudoThickness` fields in the final outputs) should
confirm that the geometric vertical grid at initialization is identical.

For Omega, tests should also verify the convergence and stability of the
p-star iteration used to relate geometric height, pressure, and pseudo-height.

For MPAS-Ocean, tests should verify that the initialized layer-thickness state
is derived from the converged geometric grid from the p-star iteration rather
than silently falling back to the legacy z-star construction.

### Testing and Validation: The capability is decomposed into inspectable Polaris steps

Date last modified: 2026/10/15

Contributors: Xylar Asay-Davis, Codex, Claude

Each major step has unit tests under `tests/ocean/realistic_global/`, mirroring
the package layout: `hydrography/woa23/test_combine.py` and
`hydrography/woa23/test_extrapolate.py` for the WOA23 preprocessing,
`forcing/jra55/test_stress.py` for the wind-stress product,
`init/test_woa23_map.py` and `init/test_jra55_map.py` for
the mapping steps and how their task counts scale, `init/test_remap_woa23.py`
and `init/test_remap_jra55.py` for what the remapped products contain and for
the coverage check, `init/test_pstar_init.py` for the column helpers,
`init/test_initial_state.py` for which upstream files the step reads, and
`init/test_forcing.py` for the per-model forcing file. `test_tasks.py` checks
that one task is registered per mesh and that its step list is what the
decomposition says, and `test_mesh_configs.py` that the per-mesh overrides
reach the right options without displacing the rest.

The `viz` step is covered by `init/test_viz_config.py` rather than by a smoke
test: it checks that the step reads only config sections that exist and that no
colormap is chosen by a literal in the plotting code. Rendering the plots needs
a real initial condition, so that part is exercised by running the task rather
than in CI.

### Testing and Validation: The workflow is practical for very large global meshes

Date last modified: 2026/03/22

Contributors: Xylar Asay-Davis, Codex

Performance testing should include at least one larger mesh case to assess
memory usage, wall-clock behavior, and whether intermediate-file boundaries are
sufficient to keep dask execution stable. These tests do not need to run in the
routine regression suite but should be documented and rerun periodically.

### Testing and Validation: The design preserves a path to future cavity support

Date last modified: 2026/03/22

Contributors: Xylar Asay-Davis, Codex

This requirement is largely architectural in the first phase. Review-based
validation is appropriate initially: the workflow should be checked for
hard-coded open-ocean assumptions at module boundaries and in intermediate
dataset schemas. Later, this requirement should graduate to task-level tests
once cavity support is added.