Support non-zero `SurfacePressure` by xylar · Pull Request #433 · E3SM-Project/Omega

xylar · 2026-06-11T17:47:50Z

Register SurfacePressure as an input/output field owned by the vertical coordinate so it is read from the initial-condition and restart files, and define all Omega pressure fields as relative (gauge) pressure. SurfacePressure was previously a zeroed placeholder in VertCoord; it is now allocated, registered as a field, and added to the State and Restart field groups next to the other prognostic variables, with a post-read step to exchange halos and update its host mirror. Because VertCoord is the sole consumer of SurfacePressure when computing 3D pressure, it owns the array directly. SurfacePressure is registered as an optional read, so it is read from the initial condition or restart file when present but defaults to zero when absent rather than causing the read to fail. Supporting this adds a general per-field optional-read capability to the IO layer.

Changes

Register `SurfacePressure` as a `VertCoord` IO field

Allocate SurfacePressure/SurfacePressureH in VertCoord and remove the placeholder TODO zero-fill
Register SurfacePressure as a 1D (NCells) field and add it to the State and Restart field groups so it is included in the initial-condition and restart streams
Keep ownership in VertCoord, its only consumer, so no array needs to be shared between classes to compute 3D pressure
Because VertCoord initializes before OceanState and both contribute to the State group, create the State/Restart groups from whichever initializes first; OceanState guards its State-group creation accordingly
Add VertCoord::initSurfacePressure(Halo*) to exchange halos and sync device→host after the initial-state/restart stream read, and call it from OceanInit after the read so the host mirror and halo are ready before the model advances

SurfacePressure is marked as an optional read, so it is not required in the Omega initial condition or restart file. When it is absent, the stream read leaves it at its fill value and initSurfacePressure defaults it to zero (0 Pa), which is the correct value for standard atmospheric forcing.

Add an optional-read capability to IO fields

Fields registered in a read stream are normally required to be present in the input file, which is why an owner like VertCoord could not simply default SurfacePressure when it lives in the shared State/Restart read. This adds a general way to make a field optional.

Add a per-field OptionalRead flag to Field (setOptionalRead/isOptionalRead, required by default) so a field may be absent from an input file
Update IOStream::readFieldData to skip a missing optional field — log, return success, and leave the array at the fill value set by attachData — instead of failing the whole stream read; required fields are unaffected
Mark VertCoord::SurfacePressure as an optional read and default it to zero in initSurfacePressure when its owned cells still hold the fill value after the read, so a missing SurfacePressure no longer aborts initialization
Add an IOStreamTest case that reads a stream containing an optional field absent from the file and verifies the read succeeds with the array left at its fill value
Document the optional-read behavior in the Field and IOStreams developer guides and the fill-values design doc

Define all pressure fields as relative (gauge) pressure

Omega's pressure fields (SurfacePressure, PressureInterface, PressureMid) are defined as relative pressure — absolute pressure minus the standard atmosphere (101325 Pa). Motivations:

TEOS-10 expects relative (gauge) pressure; Eos.h already passes Pressure * Pa2Db with no offset, consistent with this definition.
z-tilde pseudo-height ẑ = p / (ρ₀ g) equals zero at the sea surface under standard atmospheric forcing only when p is relative pressure. Under absolute pressure the surface pseudo-height would already be ≈ −10 m, which is conceptually strange.
The CF standard name sea_water_pressure is defined for absolute pressure, so it is not used for these fields; no CF standard name exists for relative pressure.

Specific changes:

VertCoord.cpp: create SurfacePressure with units Pa, long name "relative pressure at the top of the ocean column", and no CF standard name; remove the sea_water_pressure CF standard name from PressureInterface and PressureMid, update their long names to "relative pressure (gauge pressure)", and set their minimum valid value to −AtmRefP (≈ −101325 Pa, since absolute pressure ≥ 0 implies relative pressure ≥ −standard atmosphere)
Eos.h: document that Pressure inputs to the TEOS-10 specific-volume, freezing-temperature, and Brunt-Väisälä functors are relative pressure
VertCoord.h/VertCoord.cpp: clarify that computePressure computes relative pressure given relative SurfacePressure
User guide (VertCoord.md, OceanState.md): describe the pressure fields as relative, document that SurfacePressure is owned by VertCoord and is an optional read that defaults to zero, and note that it is 0 Pa for standard atmospheric forcing and represents the atmospheric anomaly in coupled runs
Developer guide (VertCoord.md, OceanState.md): document that SurfacePressure is a VertCoord field, is relative pressure, and is an optional read that defaults to zero

Checklist

xylar · 2026-06-11T17:48:08Z

This branch is based off of #428

xylar · 2026-06-11T17:55:29Z

Testing

Testing was with E3SM-Project/polaris#610

Polaris `omega_pr` suite

Baseline workdir: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260610/omega-pr-fix-fill-values
Baseline build: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260610/omega-pr-fix-fill-values/build
PR build: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260610/omega-pr-surface-pressure-init2/build
PR workdir: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260610/omega-pr-surface-pressure-init2
Machine: chrysalis
Partition: compute
Compiler: oneapi-ifx
Build type: Release
Log: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260610/omega-pr-surface-pressure-init2/polaris_omega_pr.o1231363
Result: All tests passed

A new realistic global initial condition is needed with SurfacePressure included before CTests can be expected to pass. I have not yet done this, and am hoping that the utility in E3SM-Project/polaris#562 might be used.

xylar · 2026-06-11T17:56:08Z

Documentation hasn't yet been updated.

xylar · 2026-06-12T16:31:42Z

We need to add to this that we subtract standard_atmosphere from pressure when computing SpecVol.

xylar · 2026-06-13T11:07:24Z

I will instead define pressure fields clearly as relative pressure and drop the CF standard name sea_water_pressure, see E3SM-Project/polaris#610 (comment)

xylar · 2026-06-13T12:27:21Z

Updated Testing

Again with E3SM-Project/polaris#610

Polaris `omega_pr` suite

Baseline workdir: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260611/omega-pr-baseline/
Baseline build: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260611/omega-pr-baseline/build
PR build: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260613/omega-pr-surface-pressure-init/build
PR workdir: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260613/omega-pr-surface-pressure-init
Machine: chrysalis
Partition: compute
Compiler: oneapi-ifx
Build type: Release
Log: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260613/omega-pr-surface-pressure-init/polaris_omega_pr.o1235245
Result: All tests passed

xylar · 2026-06-13T12:30:39Z

@cbegeman, @andrewdnolan and @alicebarthel, sorry for the see-sawing on relative vs. absolute pressure. I have settled on relative pressure and that is documented here now.

Please review and modify as needed (e.g. rebasing after #428 goes in). Please get @sbrus89 or someone else with the privileges to merge when you're happy, and once #428 has been merged.

cbegeman

I agree with the code changes you made here. I will circle back for testing.

Adds `surfacePressure` to Omega initial conditions ... ...but does not add it to yaml files. Thus, `surfacePressure` is not used by the model. This relies on E3SM-Project/Omega#433. Adding `surfacePressure` to the yaml files for Omega-supported tests (that do not use `State`) will be done in a further PR. This PR also adds `surfacePressure` --> `SurfacePressure` mapping from MPAS-Ocean --> Omega. We need to add `SurfacePressure` explicitly to to manufactured solution initial state, since it does not use the `State` field group.

sbrus89 · 2026-07-01T15:25:03Z

Would the Forcing class be an option for where SurfacePressure could live?

Compute the Laplacian of relative vorticity (Del2RelVortVertex) over the full vertex valid range [MinLayerVertexTop, MaxLayerVertexBot] and clamp each edge contribution to that edge's valid range [MinLayerEdgeTop, MaxLayerEdgeBot], matching VorticityAuxVars::computeVarsOnVertex. Previously Del2RelVortVertex was computed only over the narrower [MinLayerVertexBot, MaxLayerVertexTop] range (layers where every surrounding cell is active) and summed Del2Edge with no per-edge clamp. The biharmonic velocity tendency (VelocityHyperDiffOnEdge) reads Del2RelVortVertex over the edge range up to MaxLayerEdgeTop, which for a deep edge sharing a vertex with a shallower cell exceeds MaxLayerVertexTop. Those boundary-vertex layers were never computed, so on partial-bottom meshes the biharmonic term was under-computed there. This is a latent correctness bug on its own; on the fill-values branch the uncomputed layers hold FillValueReal, so the term instead blew up to ~1e45, corrupting NormalVelocity and, through the flux divergence and vertical advection, PseudoThickness and the tracers (surfaced by the new state validation in DRIVER_TEST). The wider range gives boundary vertices a valid, generally non-zero value from their active edges; inactive edges contribute zero via the per-edge clamp and Del2Edge's zeroed boundary band, so no fill value is read. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

This merge adds 5 constexpr fill value constants (FillValueI4, FillValueI8, FillValueR4, FillValueR8, FillValueReal) in namespace OMEGA, wrapping PIO_FILL_* from SCORPIO.

This merge adds private fillWithValue<T>() template helper (using if constexpr on T::value_type) and a call to it at the end of attachData<T>(). Arrays are now auto-filled with the declared fill value at attach time.

This merge removes all local fill value declarations and replaces references with the centralized constants.

We don't want fill values at every boundary edge becausee these are also valid. Instead, we want zeros except for edges between two inactive cells (e.g. both adjacent cells are below bathymetry)

Rather than assuming that the NormalVelocity read from an initial condition has correct masking, we ensure that it has fill values for inactive layers, zeros at boundary edges (both adjacent to land and to bathymetry), and leave its values unchanged for active, non-boundary layers.

We need to ensure that it is zero at active boundary edges.

Co-authored-by: Maciej Waruszewski <mwarusz@igf.fuw.edu.pl>

…variables

Add Test 5 and Test 6 to FillValueTest.cpp covering the cell and vertex layer-mask methods added in 27c02f5. Each test fills a synthetic field with a sentinel value, applies the mask, and verifies the resulting zones: cells expect FillValueReal in inactive layers and the sentinel in active layers; vertices expect the 3-zone pattern (fill / 0 / sentinel) matching the existing edge-mask test. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Document applyCellLayerMask and applyVertexLayerMask alongside the edge helpers in Section 4.5, and describe the OceanState::applyLayerMasks driver called from ocnInit. Correct stale facts in Section 5: the test lives at test/ocn/FillValueTest.cpp, uses direct comparison with Error accumulation, and now includes cell and vertex mask cases (5.5, 5.6); remove the unimplemented NetCDF attribute test. Update Requirement 2.7 for the expanded cell/vertex coverage. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Add a Layer masking subsection to the VertCoord devGuide describing zeroEdgeField, applyEdgeLayerMask, applyCellLayerMask, and applyVertexLayerMask, their zones, the inclusive active-layer convention, and the applyLayerMasks driver. Add a userGuide note explaining that inactive layers carry the fill value and boundary layers carry 0 in output. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

SurfacePressure is the relative (gauge) pressure at the top of the ocean column and the top boundary condition of the pressure field, consumed only by VertCoord::computePressure. Give VertCoord ownership of the array (device and host mirror) instead of passing it a borrowed handle from OceanState: - VertCoord allocates and registers the SurfacePressure field and adds it to the State and Restart groups (creating them if VertCoord initializes before OceanState) so it is still read from the initial-state and restart files. - VertCoord::initSurfacePressure() exchanges the halo and copies to host after the stream read; OceanInit calls it in place of the old OceanState hook. - OceanState no longer owns SurfacePressure and guards its State-group creation since VertCoord may create the group first. Also clarify in Eos that the pressure argument is relative (gauge) pressure. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Retarget the SurfacePressure size and read checks in StateTest to the VertCoord-owned array and add a host-mirror size check to VertCoordTest. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Move the SurfacePressure documentation from the OceanState user and developer guides to the VertCoord guides, describing that VertCoord owns it, that it is read from the initial-state and restart files, and that initSurfacePressure() finalizes it after the read. Note it is relative (gauge) pressure. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

xylar · 2026-07-02T15:29:09Z

@sbrus89, I have refactored the PR to leave SurfacePressure in the vertical coordinate. I think that's we're it belongs. There is still a little dance that has to happen because we read the static-in-time part of the vertical coordinate earlier, and surface pressure later in the initial state stream. But as you said, there's not a requirement about where the field lives.

I need to do some retesting to make sure the refactor actually works as expected, I just wanted to give you a heads up about the outcome.

I looked into it but I don't think the forcing is a good place for surface pressure to live. It can be updated based on coupling fields (atmospheric pressure, sea ice pressure, etc.) and still live in the vertical coordinate.

Add a per-field "optional read" flag to Field so that a field may be absent from an input file. When a read stream does not contain an optional field's variable, IOStream::readFieldData now logs and returns success, leaving the attached array at the fill value set by attachData, instead of failing the whole stream read. Fields are required by default, so existing behavior is unchanged. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Add a case to IOStreamTest that registers a field absent from the input file, marks it as an optional read, and adds it to the InitialState stream. Verify the stream read succeeds and the field array retains the fill value on all owned cells, confirming a missing optional field is skipped rather than failing the read. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Mark the VertCoord SurfacePressure field as an optional read so it is no longer required in the initial-condition or restart file. SurfacePressure stays in the State and Restart groups, so it is still written to history and restart files and read from the initial-state and restart files when present. When it is absent, initSurfacePressure detects the leftover fill value on the owned cells and defaults the field to zero before the halo exchange and host copy. Update the StateTest comment to note the field is now optional and defaults to zero. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Describe the new per-field optional-read capability in the Field and IOStreams dev guides and in the fill-values design doc: a field marked with setOptionalRead is skipped rather than failing when absent from an input file, keeping its fill value for the owning module to default. Update the VertCoord and OceanState guides to note that SurfacePressure is an optional read that defaults to zero when not present in the initial-condition or restart file. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

xylar · 2026-07-02T19:02:07Z

Testing

I'm testing in combination with E3SM-Project/polaris#629, which includes a rebase of E3SM-Project/polaris#615

Polaris `omega_pr` suite

Baseline workdir: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260630/omega-pr-develop
Baseline build: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260630/omega-pr-develop/build
PR build: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260702/omega-pr-surface-pressure2/build
PR workdir: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260702/omega-pr-surface-pressure2
Machine: chrysalis
Partition: compute
Compiler: oneapi-ifx
Build type: Release
Log: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260702/omega-pr-surface-pressure2/polaris_omega_pr.o1244770
Result:
- Diffs (1 of 12):
  - ocean/column/horiz_press_grad/salinity_gradient

The salinity_gradient test shows differences in the normal velocity tendency (horizontal pressure gradient acceleration) at machine precision. The differences are likely due to a small bug fix in E3SM-Project/polaris#629.

Polaris `horiz_press_grad` tests

This is:

ocean/column/horiz_press_grad/salinity_gradient
ocean/column/horiz_press_grad/surface_pressure_gradient
ocean/column/horiz_press_grad/temperature_gradient
ocean/column/horiz_press_grad/ztilde_gradient

No baseline in this case because I didn't have one handy. Since E3SM-Project/polaris#629 defines surface_pressure_gradient, there's no baseline for that one anyway. The others would likely show machine-precision differences as salinity_gradient did above.

PR build: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260702/omega-pr-surface-pressure2/build
PR workdir: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260702/horiz-press-grad-surf-press
Machine: chrysalis
Partition: compute
Compiler: oneapi-ifx
Build type: Release
Log: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260702/horiz-press-grad-surf-press/polaris_custom.o1244777
Result: All tests passed

CTest unit tests:

Machine: chrysalis
Compiler: oneapi-ifx
Build type: Release
Result: All tests passed
Log: /lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260702/omega-pr-surface-pressure2/build/ctests.log

This is exciting! They were failing before because SurfacePressure was missing from the test files, but now it is optional and defaults to zero.

xylar · 2026-07-02T19:04:03Z

As far as I'm concerned, this is ready for review and could go in once #428 is merged and this gets rebased onto develop.

xylar marked this pull request as draft June 11, 2026 17:55

xylar mentioned this pull request Jun 11, 2026

Add surfacePressure to Omega initial conditions E3SM-Project/polaris#610

Merged

2 tasks

xylar requested review from alicebarthel, andrewdnolan and cbegeman June 11, 2026 18:05

xylar commented Jun 12, 2026

View reviewed changes

Comment thread components/omega/src/ocn/OceanState.cpp Outdated

xylar commented Jun 12, 2026

View reviewed changes

Comment thread components/omega/src/ocn/OceanState.cpp Outdated

xylar force-pushed the omega/support-surface-pressure branch 2 times, most recently from 32511f1 to d18f6bd Compare June 12, 2026 20:04

xylar force-pushed the omega/support-surface-pressure branch 2 times, most recently from 6ef1c89 to 18c5e81 Compare June 13, 2026 11:50

xylar assigned sbrus89 Jun 13, 2026

cbegeman reviewed Jun 15, 2026

View reviewed changes

cbegeman mentioned this pull request Jun 16, 2026

Require state vars in write_initial_state_dataset, including surface pressure E3SM-Project/polaris#615

Open

1 task

xylar mentioned this pull request Jun 25, 2026

Add surface_pressure_gradient test to horiz_press_grad E3SM-Project/polaris#629

Draft

6 tasks

xylar force-pushed the omega/support-surface-pressure branch from 18c5e81 to cd06a3b Compare June 30, 2026 13:18

xylar and others added 5 commits July 2, 2026 07:44

Add a design document for fixing fill values

fbaf741

Add fill value constant

985b8e6

This merge adds 5 constexpr fill value constants (FillValueI4, FillValueI8, FillValueR4, FillValueR8, FillValueReal) in namespace OMEGA, wrapping PIO_FILL_* from SCORPIO.

Fill fields with fill values during attachData()

e96a067

This merge adds private fillWithValue<T>() template helper (using if constexpr on T::value_type) and a call to it at the end of attachData<T>(). Arrays are now auto-filled with the declared fill value at attach time.

Update fill values in all fields

0a7aec5

This merge removes all local fill value declarations and replaces references with the centralized constants.

xylar and others added 20 commits July 2, 2026 07:50

Reorder 2 cases so we attach first, then fill in the array

3d2c737

Fill fluxes on edges first with zeros over full physical range

4505f49

We don't want fill values at every boundary edge becausee these are also valid. Instead, we want zeros except for edges between two inactive cells (e.g. both adjacent cells are below bathymetry)

Add fill-value CTests

cc58ac2

Update the docs

48e65e6

Fix bounds on DivHU computation to avoid inactive layers

7437516

Mask NormalVelocity in time stepper

56192c7

We need to ensure that it is zero at active boundary edges.

Apply suggestions from code review

25c9a32

Co-authored-by: Maciej Waruszewski <mwarusz@igf.fuw.edu.pl>

Apply edge layer mask to NormalVelocity in OceanState constructor

3961f3d

Create applyCellLayerMask, applyVertexLayerMask

2052fa2

Use single method in OceanRun for applying fill values to OceanState …

e145479

…variables

Fixup rebase on vmix shear changes

5e06ba3

Fixup rebase on surface stress changes

35bb0b0

Correct comment on where BottomGeomDepth gets read from

bbbbfdc

Test SurfacePressure as a VertCoord field

e7db968

Retarget the SurfacePressure size and read checks in StateTest to the VertCoord-owned array and add a host-mirror size check to VertCoordTest. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

xylar force-pushed the omega/support-surface-pressure branch from cd06a3b to 3fc9579 Compare July 2, 2026 14:39

xylar and others added 4 commits July 2, 2026 12:23

xylar marked this pull request as ready for review July 2, 2026 19:02

xylar requested a review from sbrus89 July 2, 2026 19:04

Uh oh!

Conversation

xylar commented Jun 11, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Changes

Register SurfacePressure as a VertCoord IO field

Add an optional-read capability to IO fields

Define all pressure fields as relative (gauge) pressure

Uh oh!

xylar commented Jun 11, 2026

Uh oh!

xylar commented Jun 11, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Testing

Polaris omega_pr suite

Uh oh!

xylar commented Jun 11, 2026

Uh oh!

Uh oh!

Uh oh!

xylar commented Jun 12, 2026

Uh oh!

xylar commented Jun 13, 2026

Uh oh!

xylar commented Jun 13, 2026

Updated Testing

Polaris omega_pr suite

Uh oh!

xylar commented Jun 13, 2026

Uh oh!

cbegeman left a comment

Choose a reason for hiding this comment

Uh oh!

sbrus89 commented Jul 1, 2026

Uh oh!

xylar commented Jul 2, 2026

Uh oh!

xylar commented Jul 2, 2026

Testing

Polaris omega_pr suite

Polaris horiz_press_grad tests

CTest unit tests:

Uh oh!

xylar commented Jul 2, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

xylar commented Jun 11, 2026 •

edited

Loading

Register `SurfacePressure` as a `VertCoord` IO field

xylar commented Jun 11, 2026 •

edited

Loading

Polaris `omega_pr` suite

Polaris `omega_pr` suite

Polaris `omega_pr` suite

Polaris `horiz_press_grad` tests