diff --git a/components/omega/doc/design/FillValues.md b/components/omega/doc/design/FillValues.md
new file mode 100644
index 000000000000..0602e6e32bb1
--- /dev/null
+++ b/components/omega/doc/design/FillValues.md
@@ -0,0 +1,501 @@
+(omega-design-fill-values)=
+# Fill Values
+
+**Table of Contents**
+1. [Overview](#1-overview)
+2. [Requirements](#2-requirements)
+3. [Algorithmic Formulation](#3-algorithmic-formulation)
+4. [Design](#4-design)
+5. [Verification and Testing](#5-verification-and-testing)
+
+## 1 Overview
+
+Omega supports output of ocean fields to NetCDF files via I/O streams. Each field
+carries a declared fill value stored as metadata under the `_FillValue` attribute,
+which marks entries that are undefined or outside the active domain. Currently, fill
+value constants are declared locally in each module with inconsistent values
+(typically `-9.99e30` for real-valued fields, `-999` for integers). These values do
+not match the NetCDF-C standard fill values that analysis tools and visualization
+software commonly expect. Additionally, the Kokkos arrays backing these fields are
+never explicitly initialized with fill values: inactive ocean layers (model layers
+below `MaxLayerCell`) hold uninitialized memory rather than a well-defined sentinel,
+and the distinction between "no data" and "zero" is lost.
+
+This design standardizes fill values across all Omega fields to match the NetCDF-C
+standard, ensures that field arrays are automatically initialized with the correct
+fill value when data is attached to a `Field`, and verifies through CTests that the
+resulting behavior is correct for inactive layers and for edge/vertex fields at
+domain boundaries.
+
+## 2 Requirements
+
+### 2.1 Requirement: Centralized fill value constants
+
+Fill value constants must be defined in a single location in Omega that is accessible
+to all modules. Scatter of locally-declared constants has produced inconsistencies
+(e.g., `-9.99e30` vs `-9.99E+30`) and makes future maintenance error-prone.
+
+### 2.2 Requirement: Standard fill values matching NetCDF-C
+
+Fill value constants must exactly match the NetCDF-C standard values (`NC_FILL_*`)
+as exposed by the SCORPIO library (`PIO_FILL_*` in `pio.h`). Using standard values
+ensures compatibility with NetCDF-aware tools (ncview, Xarray, NCO, etc.) that
+recognize and handle fill values automatically. The required values are:
+
+| Type | Omega name       | Value                         |
+|------|-----------------|-------------------------------|
+| I4   | `FillValueI4`   | -2147483647 (`PIO_FILL_INT`)  |
+| I8   | `FillValueI8`   | `PIO_FILL_INT64`              |
+| R4   | `FillValueR4`   | 9.9692099683868690e+36f (`PIO_FILL_FLOAT`) |
+| R8   | `FillValueR8`   | 9.9692099683868690e+36 (`PIO_FILL_DOUBLE`) |
+| Real | `FillValueReal` | Alias for `FillValueR4` or `FillValueR8` depending on build |
+
+### 2.3 Requirement: Automatic array initialization at attach time
+
+When `Field::attachData()` is called, the attached Kokkos array must be
+automatically initialized with the standard fill value for the array's element
+type before any compute routine runs. The fill value is deduced from
+`typename T::value_type` at `attachData()` time using the constants in
+`FillValues.h` (`FillValueI4`, `FillValueI8`, `FillValueR4`, or `FillValueR8`),
+so the caller does not need to specify it. This guarantees that inactive ocean
+layers (below `MaxLayerCell`) and any other undefined entries contain a
+well-defined sentinel in output, without requiring each module to manually
+perform the fill.
+
+### 2.4 Requirement: No fill value argument in Field::create()
+
+The `FillValue` argument is removed from `Field::create()`. Because the fill
+value is always one of the standard constants determined by the array element
+type, specifying it at field creation time is redundant. The fill value and the
+`_FillValue` metadata entry are set automatically when `attachData()` is called.
+
+### 2.5 Requirement: Inactive layers contain fill values in output
+
+After compute routines execute, model layers with index `k > MaxLayerCell(ICell)`
+must contain the declared fill value in output. With automatic initialization at
+attach time (Requirement 2.3) and compute routines that only write to active layers,
+this is satisfied without changes to individual compute functions.
+
+### 2.6 Desired: Valid values at domain boundary edges and vertices
+
+Edge and vertex fields at valid boundaries adjacent to land or to bathymetry must
+contain computed valid values, not fill values. Switching from zero-initialization
+to fill-value initialization must not silently break any boundary computation that
+previously assumed zero defaults. This mirrors a problem encountered in MPAS-Ocean
+when initializing `NormalVelocity` with a fill value rather than zero: code that
+implicitly assumed zero at land boundaries produced incorrect results. During
+implementation, all compute routines for edge/vertex fields must be audited to
+confirm that every valid boundary entry receives an explicit computed value. As an
+example of what such an audit can reveal: `TracerHorzAdvOnCell` accumulated its
+internal working array `HighOrderFlxHorz` over a cell's full active layer range
+without clamping to each contributing edge's active range, reading fill values from
+edges shallower than the cell — see Section 4.5 for the resolution.
+
+### 2.7 Requirement: CTests
+
+CTests must verify Requirements 2.1–2.3, 2.5, and Desired Requirement 2.6. The
+three-zone boundary behavior (2.6) is exercised for `NormalVelocity` via the edge
+mask; the two-zone cell and vertex masks are verified with synthetic cell and
+vertex fields.
+
+## 3 Algorithmic Formulation
+
+No new numerical algorithm is introduced. The key behavioral change is in
+`Field::attachData<T>()`, which gains a fill step equivalent to:
+
+```c++
+Kokkos::deep_copy(InDataArray, fill_scalar);
+```
+
+where `fill_scalar` is the typed fill value extracted from the field's `FieldMeta`
+map (stored as `std::any` under key `"_FillValue"`), and the element type
+`typename T::value_type` of the Kokkos view determines the `std::any_cast` type.
+`Kokkos::deep_copy` with a scalar source broadcasts that value to every element of
+the view and dispatches correctly for both host and device memory spaces.
+
+## 4 Design
+
+### 4.1 Data types and parameters
+
+#### 4.1.1 Constants: FillValue definitions in FillValues.h
+
+The fill value constants are defined in a new lightweight header
+`components/omega/src/base/FillValues.h`. The values match the NetCDF-C `NC_FILL_*`
+constants (which are identical to the SCORPIO `PIO_FILL_*` constants), but are
+written as numeric literals so that the header can be included from any context
+without pulling in `pio.h` or `mpi.h`, which carry strict include-order requirements.
+
+A type-indexed variable template `FillValue<T>` is the primary definition. Its
+primary template is intentionally left undefined so that instantiation with an
+unsupported type produces a link error. Explicit specializations are provided for
+each supported scalar type:
+
+```c++
+template <typename T>
+constexpr T FillValue;  // primary template intentionally undefined
+
+template <> inline constexpr I4 FillValue<I4> = -2147483647;             // NC_FILL_INT
+template <> inline constexpr I8 FillValue<I8> = -9223372036854775806LL;  // NC_FILL_INT64
+template <> inline constexpr R4 FillValue<R4> = 9.9692099683868690e+36f; // NC_FILL_FLOAT
+template <> inline constexpr R8 FillValue<R8> = 9.9692099683868690e+36;  // NC_FILL_DOUBLE
+```
+
+Named aliases are provided for readability in comparison and test code:
+
+```c++
+constexpr I4   FillValueI4   = FillValue<I4>;
+constexpr I8   FillValueI8   = FillValue<I8>;
+constexpr R4   FillValueR4   = FillValue<R4>;
+constexpr R8   FillValueR8   = FillValue<R8>;
+#if defined(SINGLE_PRECISION)
+constexpr Real FillValueReal = FillValue<R4>;
+#else
+constexpr Real FillValueReal = FillValue<R8>;
+#endif
+```
+
+`FillValues.h` is included from `Field.h`. Because every module that creates or
+uses fields already includes `Field.h` (directly or transitively), no additional
+include directives are needed in module source files.
+
+#### 4.1.2 Class/struct changes
+
+No new classes or structs are introduced. The changes are confined to:
+- new header `FillValues.h` with constants (Section 4.1.1)
+- `FillValues.h` included from `Field.h`
+- the `attachData` template methods in `Field.h` / `Field.cpp` (Section 4.2)
+- removal of local fill value declarations in module files (Section 4.3)
+
+### 4.2 Methods
+
+#### 4.2.1 Field::attachData() auto-fill
+
+The `attachData<T>()` instance method and the `attachFieldData<T>()` static helper
+in `components/omega/src/infra/Field.h` gain a private fill step controlled by an
+optional `FillOnAttach` parameter (default `true`):
+
+```c++
+template <typename T>
+void attachData(const T &InDataArray, bool FillOnAttach = true) {
+    OMEGA_ASSERT(isKokkosArray<T>,
+                 "Field::attachData requires Kokkos array as input");
+    DataArray = std::make_shared<T>(InDataArray);
+    DataType  = checkArrayType<T>();
+    MemLoc    = findArrayMemLoc<T>();
+    if (FillOnAttach)
+        fillWithValue<T>(InDataArray); // auto-fill with type-deduced FillValue
+}
+
+template <typename T>
+static void attachFieldData(const std::string &FieldName,
+                            const T &InDataArray, bool FillOnAttach = true);
+```
+
+The default `FillOnAttach = true` preserves the fill behavior for the normal
+initialization workflow (fresh arrays attached before any data is written). Pass
+`FillOnAttach = false` when re-attaching an existing data-filled array where the
+fill must not overwrite computed values — see Section 4.2.3 for the canonical use case.
+
+The private helper `fillWithValue<T>()` performs the following steps:
+
+1. Derives the element type as `typename T::non_const_value_type` (using the
+   non-const form so that const view types resolve to the same underlying type).
+2. Looks up the fill value via `FillValue<ValType>`, a constexpr variable template
+   specialization from `FillValues.h`. Supporting a new type requires only adding
+   one specialization there; no changes to `fillWithValue` are needed.
+3. Stores the typed fill value in `FieldMeta["_FillValue"]` (the NetCDF/CF
+   standard attribute name), making it available to IO and metadata queries.
+4. Calls `Kokkos::deep_copy(InDataArray, scalar)` to broadcast the scalar to every
+   element of the view. This works for host views (`HostArray*`), device views
+   (`Array*`), and all dimensionalities (1D–5D) without additional dispatch.
+5. Only the view being attached is filled; no assumption is made about whether the
+   other memory space (host or device) has been attached.
+
+Ordering invariants preserved by this design:
+- `Field::create()` is always called before `attachData()`, so the field exists
+  in `AllFields` before fill time.
+- The `"_FillValue"` metadata entry is set at `attachData()` time, which precedes
+  any IO write operation in the initialization workflow.
+- Restart reads happen after `attachData()` and overwrite fill entries with data
+  from the restart file.
+- Compute routines run after initialization and overwrite active entries, leaving
+  inactive entries (below `MaxLayerCell`) at the fill value.
+
+The fill value also serves as the "not read" sentinel for optional reads. A field
+marked with `Field::setOptionalRead(true)` that is absent from an input file is
+skipped by the IOStream read and keeps the fill value from `attachData()`; the
+owning module detects the fill value afterward and substitutes a default. For
+example, `VertCoord::SurfacePressure` defaults to zero when it is not present in
+the initial-condition or restart file.
+
+#### 4.2.2 No per-module manual initialization required
+
+Because `attachData()` handles initialization automatically, no module needs to
+call an explicit fill routine. Modules may still call `Kokkos::deep_copy` or
+similar to set specific values, but this is independent of fill-value management.
+
+#### 4.2.3 Ordering invariant and anti-patterns
+
+**Invariant**: `attachData()` with `FillOnAttach = true` (the default) must always
+be called *before* any real values are written to the array. Any data written before
+such a call is silently overwritten by the auto-fill step and is dead code.
+
+Three anti-patterns discovered and corrected during implementation:
+
+**Class A — Wasted pre-fill**: A `deepCopy` (or equivalent write) is applied to an
+array before `defineFields()` calls `attachData()` on it. Because `attachData()`
+overwrites the entire array, the earlier write has no effect.
+
+```c++
+// WRONG: initial value is overwritten by attachData()
+deepCopy(SpecVol, 1.0_Real / RhoSw);
+defineFields(); // calls attachData(), fills SpecVol with FillValueReal
+```
+
+**Class B — Compute-time re-initialization**: A `deepCopy(arr, 0)` (or any scalar)
+at the top of a recurring compute method resets the entire array on every timestep,
+overwriting the correct fill values in inactive layers with zero.
+
+```c++
+// WRONG: destroys fill values in inactive layers every timestep
+void Foo::compute(...) {
+    deepCopy(LocArr, 0);
+    parallelFor({NCellsAll}, KOKKOS_LAMBDA(I4 ICell) {
+        // writes only active layers KMin..KMax
+    });
+}
+```
+
+The correct pattern for Classes A and B is to do nothing: `attachData()` sets the
+fill value once at initialization, and compute kernels that iterate only over active
+layers naturally leave inactive layers at the fill value.
+
+**Class C — Time-level pointer update**: `attachData()` is called each timestep to
+re-point the IO field to the newly-current time-level array (e.g., in
+`OceanState::updateTimeLevels()` and `Tracers::updateTimeLevels()`). The array
+already contains the just-computed state. Calling `attachData()` with the default
+`FillOnAttach = true` here would silently destroy the computed values before IO
+writes them. This is the most dangerous anti-pattern because the model runs to
+completion without aborting, but all output fields contain fill values.
+
+```c++
+// WRONG: destroys computed state before IO can write it
+void OceanState::updateTimeLevels() {
+    CurTimeIndex = (CurTimeIndex + 1) % NTimeLevels;
+    Field::attachFieldData<Array2DReal>(NormalVelocityFldName,
+                                        NormalVelocity[CurTimeIndex]); // fills!
+}
+
+// CORRECT: pointer update only, computed values preserved
+void OceanState::updateTimeLevels() {
+    CurTimeIndex = (CurTimeIndex + 1) % NTimeLevels;
+    Field::attachFieldData<Array2DReal>(NormalVelocityFldName,
+                                        NormalVelocity[CurTimeIndex], false);
+}
+```
+
+A related instance occurs in the IO write path when the time-coordinate array is
+populated before `attachFieldData("time", OutTime)` is called. The correct pattern
+is to attach first (which fills the single-element array with the fill value), then
+assign the elapsed time to `OutTime(0)`. Because Kokkos host views share backing
+memory, the assignment is visible through the stored pointer before the field is
+written to the file.
+
+### 4.3 Module-level changes: removing FillValue from Field::create() calls
+
+The `FillValue` parameter is removed from `Field::create()`. All ~50 call sites
+across Omega drop the fill value argument. The pattern is:
+
+```c++
+// BEFORE — explicit fill value argument
+Field::create(FieldName, Description, Units, StdName,
+              ValidMin, ValidMax, FillValueReal, NDims, DimNames);
+
+// AFTER — fill value is deduced automatically from the array type at attachData()
+Field::create(FieldName, Description, Units, StdName,
+              ValidMin, ValidMax, NDims, DimNames);
+```
+
+The `Tracers::define()` helper in `src/ocn/Tracers.h` / `Tracers.cpp` and
+`TracerDefs.inc` also have their `FillValue` parameter and arguments removed,
+since `Tracers::define()` calls `Field::create()` internally.
+
+### 4.4 Boundary edge and vertex handling
+
+During implementation, all compute routines that populate fields defined on mesh
+edges or vertices must be audited to confirm that every valid boundary entry —
+including edges and vertices adjacent to land cells and edges at the base of the
+water column — receives an explicit computed value rather than relying on
+zero-initialization. Where such reliance is found, an explicit assignment must be
+added. The CTest in Section 5.5 detects regressions after the change.
+
+### 4.5 Edge flux field masking
+
+Edge fields that represent fluxes normal to the edge face — such as
+`NormalVelocityTend`, `VelocityDel2Aux.Del2Edge`, and `NormalVelocity` itself —
+must carry zeros (not fill values) in the boundary layer sub-ranges
+`[MinLayerEdgeTop, MinLayerEdgeBot)` and `(MaxLayerEdgeTop, MaxLayerEdgeBot]`,
+where one neighboring cell is active but the other is not. Two helper methods
+are added to `VertCoord` (`components/omega/src/ocn/VertCoord.h/.cpp`):
+
+**`zeroEdgeField(Array2DReal &Arr, I4 NEdgesAll)`** — zeros all layers in
+`[MinLayerEdgeTop, MaxLayerEdgeBot]`. Used before recomputing a flux-type edge
+field each time step so that boundary edges show 0, not `FillValueReal`. Layers
+outside the valid range retain their fill value from `attachData()`. Called for:
+- `NormalVelocityTend` at the start of `Tendencies::computeVelocityTendenciesOnly()`
+- `VelocityDel2Aux.Del2Edge` before the `edgeAuxState2` kernel in
+  `AuxiliaryState::computeMomAux()`, because the Del4 hyperdiffusion kernel
+  accumulates `Del2Edge` over all edges sharing a cell or vertex; fill values in
+  boundary layers of a neighbouring edge would corrupt `Del2DivCell` and
+  `Del2RelVortVertex`.
+
+**`applyEdgeLayerMask(Array2DReal &Arr, I4 NEdgesAll)`** — enforces the full
+three-zone mask on an edge field after IC or restart read: layers outside
+`[MinLayerEdgeTop, MaxLayerEdgeBot]` are set to `FillValueReal`; layers inside
+`[MinLayerEdgeTop, MaxLayerEdgeBot]` but outside the active range
+`[MinLayerEdgeBot, MaxLayerEdgeTop]` are set to 0; active layers are left
+unchanged.
+
+**`applyCellLayerMask(Array2DReal &Arr, I4 NCellsAll)`** — enforces the
+two-zone mask on a cell field after IC or restart read: layers outside the
+active range `[MinLayerCell, MaxLayerCell]` are set to `FillValueReal`; active
+layers are left unchanged. Cell fields have no boundary (zero) zone because a
+cell column is either active or inactive at a given layer — there is no
+neighbor-dependent partial-activity range as there is for edges and vertices.
+
+**`applyVertexLayerMask(Array2DReal &Arr, I4 NVerticesAll)`** — enforces the
+two-zone mask on a vertex field: layers outside `[MinLayerVertexTop,
+MaxLayerVertexBot]` are set to `FillValueReal`; active layers are left
+unchanged. Unlike an edge, a vertex does **not** have a zeroed boundary zone: a
+boundary vertex with one or more active surrounding cells holds valid, generally
+non-zero data (for example, relative vorticity at such a vertex is computed from
+its active surrounding edges — see `VorticityAuxVars::computeVarsOnVertex`, which
+loops over the full `[MinLayerVertexTop, MaxLayerVertexBot]` range). Zeroing that
+boundary band would discard real signal, so the whole valid range is kept. There
+is currently no vertex-based state field read from IC/restart, so this method has
+no production caller yet; it is provided for completeness and is exercised by the
+CTest (Section 5.6).
+
+All three `apply*LayerMask` methods use the inclusive `[Min, Max]` layer
+convention consistent with `computeGeomZHeight`/`computePressure`. They are
+driven through `OceanState::applyLayerMasks(TimeLevel)`
+(`src/ocn/OceanState.cpp`), which is called from `ocnInit` (`OceanInit.cpp`)
+after `exchangeHalo` and before `copyToHost`. `applyLayerMasks` applies the
+edge mask to `NormalVelocity` and the cell mask to `PseudoThickness`,
+`Temperature`, and `Salinity` — the state and tracer fields read from
+IC/restart.
+
+Fields that are *not* fluxes — `MeanPseudoThickEdge`, `FluxPseudoThickEdge` —
+are genuinely undefined at boundary edges (interpolating thickness when one
+neighboring cell is land has no physical meaning) and therefore correctly retain
+`FillValueReal` in those layers; no zeroing is applied to them.
+`VorticityAux` on edges is already computed over the full valid range
+`[MinLayerEdgeTop, MaxLayerEdgeBot]`, so it requires no additional masking.
+
+**`HighOrderFlxHorz` in `TracerHorzAdvOnCell`** — the 3-D working array
+`(NTracers, NEdgesSize, NVertLayers)` that accumulates horizontal tracer flux
+contributions per edge requires separate treatment. The edge-pass overload of
+`TracerHorzAdvOnCell::operator()` writes `HighOrderFlxHorz(L, IEdge, K)` only
+for layers `K ∈ [0, MaxLayerEdgeTop(IEdge)]`; layers beyond that index retain
+their fill values from `attachData()`. The cell-pass overload then reads
+`HighOrderFlxHorz` over the cell's own active range `[MinLayerCell, MaxLayerCell]`,
+which for cells adjacent to shallower edges extends past `MaxLayerEdgeTop` of
+those edges — reading fill values and corrupting active-layer tracer tendencies.
+The fix is to add `MinLayerEdgeBot` and `MaxLayerEdgeTop` as members of
+`TracerHorzAdvOnCell` and clamp the cell-pass inner loop to
+`[max(KStart, MinLayerEdgeBot(IEdge)), min(KEnd, MaxLayerEdgeTop(IEdge)+1))`.
+This follows the identical pattern already used by `PseudoThicknessFluxDivOnCell`,
+`TracerDiffOnCell`, `TracerAuxVars::computeVarsOnCells`, and every other
+cell-level flux-accumulation kernel, and avoids the need to zero
+`HighOrderFlxHorz` before the edge pass each timestep.
+
+## 5 Verification and Testing
+
+A test `test/ocn/FillValueTest.cpp` is added and registered as `FILL_VALUE_TEST`
+with 8 MPI tasks (building `testFillValue.exe`) using `add_omega_test()` in the
+test `CMakeLists.txt`. The test performs a full ocean initialization (matching
+the `StateTest` pattern: `MachEnv`, logging, Config, `TimeStepper`, IO, `Field`,
+`IOStream`, `Decomp`, `Halo`, `HorzMesh`, `VertCoord`, `Tracers`,
+`AuxiliaryState`, `PressureGrad`, `Tendencies`, `OceanState`), reads the initial
+state, exchanges halos, and applies `applyEdgeLayerMask` to `NormalVelocity`
+before running its assertions. Each test counts mismatches and, on any failure,
+accumulates an `Error(ErrorCode::Fail, ...)` into a top-level `Error` whose
+final state sets the process return code.
+
+### 5.1 Test: fill constant values
+
+Verify that each Omega fill value constant exactly equals its NetCDF-C
+counterpart (`FillValueI4 == NC_FILL_INT`, `FillValueI8 == NC_FILL_INT64`,
+`FillValueR4 == NC_FILL_FLOAT`, `FillValueR8 == NC_FILL_DOUBLE`), comparing
+directly against the `<netcdf.h>` `NC_FILL_*` macros.
+
+Tests requirements: 2.1, 2.2, 2.7
+
+### 5.2 Test: attachData auto-fill
+
+Create a `Field`, allocate a Kokkos host array whose elements are all set to a
+distinct sentinel (`0`), then call `attachData()`. Verify that after the call
+every element of the array equals `FillValueR8`.
+
+Tests requirements: 2.3, 2.7
+
+### 5.3 Test: inactive layers contain fill values
+
+After `VertCoord` initialization, verify that the cell field `GeomZMid` carries
+`FillValueReal` in every inactive layer (`k >= MaxLayerCell(ICell)`) for all
+owned cells. This confirms that the auto-fill applied at `attachData()` time
+survives initialization for inactive layers.
+
+Tests requirements: 2.5, 2.7
+
+### 5.4 Test: NormalVelocity three-zone fill/zero/real pattern
+
+After full state initialization, IC read, halo exchange, and
+`VertCoord::applyEdgeLayerMask()`, verify that `NormalVelocity` satisfies the
+three-zone invariant for every owned edge:
+
+- **Fully inactive layers** (outside `[MinLayerEdgeTop, MaxLayerEdgeBot]`):
+  equal `FillValueReal`.
+- **Boundary layers** (inside `[MinLayerEdgeTop, MaxLayerEdgeBot]` but outside
+  `[MinLayerEdgeBot, MaxLayerEdgeTop]`): equal `0`.
+- **Active layers** (`[MinLayerEdgeBot, MaxLayerEdgeTop]`): not equal to
+  `FillValueReal` (for a quiescent IC the value is `0`, which is valid).
+
+This guards against two regressions: (a) the original zero-initialization
+assumption (boundary layers silently filled with zero before `attachData()`
+was changed) and (b) IC files that store zeros in inactive layers, which
+`applyEdgeLayerMask` must replace with `FillValueReal`.
+
+Tests requirements: 2.5, 2.6, 2.7
+
+### 5.5 Test: applyCellLayerMask two-zone pattern
+
+Allocate a synthetic cell field, fill it with a sentinel value (≠ 0 and
+≠ `FillValueReal`), apply `VertCoord::applyCellLayerMask()`, and verify for every
+owned cell that:
+
+- **Inactive layers** (outside `[MinLayerCell, MaxLayerCell]`): equal
+  `FillValueReal`.
+- **Active layers** (`[MinLayerCell, MaxLayerCell]`): retain the sentinel value.
+
+Using a sentinel rather than IC data gives an exact check that the mask
+overwrites only the inactive layers.
+
+Tests requirements: 2.5, 2.7
+
+### 5.6 Test: applyVertexLayerMask two-zone pattern
+
+Allocate a synthetic vertex field, fill it with the same sentinel value, apply
+`VertCoord::applyVertexLayerMask()`, and verify for every owned vertex that:
+
+- **Inactive layers** (outside `[MinLayerVertexTop, MaxLayerVertexBot]`): equal
+  `FillValueReal`.
+- **Active layers** (`[MinLayerVertexTop, MaxLayerVertexBot]`): retain the
+  sentinel value.
+
+Vertices have no zeroed boundary zone (a boundary vertex holds valid data), so
+the whole valid range is checked for the sentinel. Because no vertex-based state
+field is read from IC/restart, this synthetic test is the primary coverage for
+`applyVertexLayerMask`.
+
+Tests requirements: 2.5, 2.7
diff --git a/components/omega/doc/devGuide/Field.md b/components/omega/doc/devGuide/Field.md
index 85244430b133..b295ff9609d7 100644
--- a/components/omega/doc/devGuide/Field.md
+++ b/components/omega/doc/devGuide/Field.md
@@ -37,7 +37,6 @@ Fields are created with standard metadata using
                  StdName,     ///< [in] CF standard Name (string)
                  ValidMin,    ///< [in] min valid field value (same type as
                  ValidMax,    ///< [in] max valid field value  field data)
-                 FillValue,   ///< [in] scalar used for undefined entries
                  NumDims,     ///< [in] number of dimensions (int)
                  Dimensions,  ///< [in] dim names (vector of strings)
                  TimeDependent   ///< [in] (opt, default true) if time varying
@@ -49,9 +48,13 @@ not exist, an empty string can be provided. This is uncommon for most fields
 since the CF conventions maintain a fairly complete list, but can be the case
 for some intermediate calculations or unique analyses. If there is no
 restriction on valid range, an appropriately large range should be provided for
-the data type. Similarly, if a FillValue is not being used, a very unique
-number should be supplied to prevent accidentally treating valid data as a
-FillValue. The optional TimeDependent argument can be omitted and is assumed
+the data type. The fill value is not specified when creating a field; it is
+automatically deduced from the element type of the array passed to
+`attachData()` and set to the corresponding standard constant from `FillValues.h`
+(`FillValueI4`, `FillValueI8`, `FillValueR4`, or `FillValueR8`). These match
+the NetCDF-C `NC_FILL_*` standard values recognized by analysis tools such as
+ncview, Xarray, and NCO. The optional TimeDependent argument can be omitted
+and is assumed
 to be true by default. Fields with this attribute will be output with the
 unlimited time dimension added. Time should not be added explicitly in the
 dimension list since it will be added during I/O. Fields that do not change
@@ -113,6 +116,13 @@ is provided using the field name:
                     InDataArray  // [in] Array with data to attach
                     );
 ```
+When `attachData` is called, every element of the attached array is automatically
+initialized to the field's declared fill value. This ensures that inactive entries
+(e.g., ocean layers below `MaxLayerCell`) hold a well-defined NetCDF-standard
+sentinel in output without any explicit initialization in module code. Subsequent
+compute calls overwrite active entries; the fill value persists only where no valid
+data is written.
+
 Note that the data is assumed to reside in only one location so if a mirror
 array exists (eg if replicated on host and device), a separate Field may be
 needed. However, it is is better to define only one location and allow the
@@ -182,6 +192,18 @@ The first determines whether the unlimited time dimension should be added
 during IO operations. The second determines whether any of the dimensions
 are distributed across MPI tasks so that parallel IO is required.
 
+A field can also be marked as an optional read, which controls what happens
+when the field's variable is missing from an input file:
+```c++
+   MyField->setOptionalRead(true);
+   bool IsOptionalRead = MyField->isOptionalRead();
+```
+Fields are required by default (`isOptionalRead()` returns `false`). When a
+field is marked optional and its variable is absent from a read stream's file,
+the [IOStream](#omega-dev-iostreams) read skips it, leaving the attached array
+at its [fill value](#omega-design-fill-values) instead of failing. The owning
+module is responsible for detecting that fill value and substituting a default.
+
 The data and metadata stored in a field can be retrieved using several
 functions.  To retrieve a pointer to the full Field, use:
 ```c++
diff --git a/components/omega/doc/devGuide/IOStreams.md b/components/omega/doc/devGuide/IOStreams.md
index 19884a450ef3..f1ea8af3e23b 100644
--- a/components/omega/doc/devGuide/IOStreams.md
+++ b/components/omega/doc/devGuide/IOStreams.md
@@ -72,6 +72,20 @@ input stream is read using:
 The returned error code typically means that a field in the stream could
 not be found in the input file - most other errors abort immediately. The
 calling routine is then responsible for deciding what action to take.
+
+By default every field listed in a stream's contents must be present in the
+input file or the read returns an error. A field can instead be marked as an
+optional read using
+```c++
+   FieldPtr->setOptionalRead(true);
+```
+When such a field's variable is missing from the file, the read logs an
+informational message, leaves the field's attached array at the fill value set
+by `attachData` (see [Fill Values](#omega-design-fill-values)), and continues
+without contributing an error. The owning module is then responsible for
+detecting that fill value after the read and substituting a default value. This
+is how `VertCoord` allows `SurfacePressure` to be absent from the
+initial-condition or restart file and default to zero.
 The ReqMetadata argument is a variable of type Metadata (defined in Field but
 essentially a ``std::map<std::string, std::any>`` for the name/value pair).
 This variable should include the names of global metadata that are desired
diff --git a/components/omega/doc/devGuide/OceanState.md b/components/omega/doc/devGuide/OceanState.md
index f74bf2eef9c8..d6cb4db31d45 100644
--- a/components/omega/doc/devGuide/OceanState.md
+++ b/components/omega/doc/devGuide/OceanState.md
@@ -21,7 +21,12 @@ OceanState::create(const std::string &Name, ///< [in] Name for mesh
 );
 ```
 allocates the `NormalVelocity` and `PseudoThickness` arrays for a given number of time levels.
-The current time level is then registered with the IO infrastructure.
+The current time level of `NormalVelocity` and `PseudoThickness` are registered with the IO
+infrastructure and added to the `State` and `Restart` field groups.
+The `SurfacePressure` array is owned by [`VertCoord`](omega-dev-vert-coord) but is also added to
+these two groups so it is written to restart files and read from the initial-condition and restart
+files when present. Its read is optional, so it defaults to zero when absent (see the
+[`VertCoord`](omega-dev-vert-coord) guide).
 
 After initialization, the default state object can be retrieved via:
 ```
@@ -71,6 +76,9 @@ HostArray2DReal NormVelH = State->getNormalVelocityH(TimeLevel);
 for the host arrays. These functions return the arrays directly and will abort
 via `OMEGA_REQUIRE` if an invalid `TimeLevel` is provided.
 
+The `SurfacePressure` array used as the top boundary condition for layer pressures is owned by
+[`VertCoord`](omega-dev-vert-coord); see that guide for how it is accessed and initialized.
+
 The time level convention is:
 | time level | `TimeLevel` |
 |------------|-------------|
diff --git a/components/omega/doc/devGuide/Tracers.md b/components/omega/doc/devGuide/Tracers.md
index 5bb55c1e158a..f455df08e880 100644
--- a/components/omega/doc/devGuide/Tracers.md
+++ b/components/omega/doc/devGuide/Tracers.md
@@ -24,17 +24,23 @@ static void defineAllTracers() {
           "sea_water_potential_temperature", ///< [in] CF standard Name
           -273.15,                           ///< [in] min valid field value
           100.0,                             ///< [in] max valid field value
-          1.e33,                             ///< [in] value for undef entries
+          FillValueReal,                     ///< [in] value for undef entries
           IndxTemp);                         ///< [out] (optional) static index
 
-   define("Salt", "Salinity", "psu", "sea_water_salinity", 0.0, 50.0, 1.e33,
-          IndxSalt);
-   define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, 1.e33);
-   define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, 1.e33);
-   define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, 1.e33);
+   define("Salt", "Salinity", "psu", "sea_water_salinity", 0.0, 50.0,
+          FillValueReal, IndxSalt);
+   define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, FillValueReal);
+   define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, FillValueReal);
+   define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, FillValueReal);
 }
 ```
 
+The `FillValueReal` constant is defined in `FillValues.h` and matches the
+NetCDF-C standard fill value (`NC_FILL_DOUBLE` or `NC_FILL_FLOAT` depending on
+the build precision). It is automatically in scope in `TracerDefs.inc` via
+`Field.h`. New tracers should always use `FillValueReal` (or `FillValueI4` /
+`FillValueI8` for integer fields) rather than hardcoded literals.
+
 To add a new tracer, simply call the `define` function with the appropriate
 arguments. Index argument is optional one that allows to access the tracer
 data using the given tracer index variable.
diff --git a/components/omega/doc/devGuide/VertCoord.md b/components/omega/doc/devGuide/VertCoord.md
index 90bb3e017c84..fd2eaccceee1 100644
--- a/components/omega/doc/devGuide/VertCoord.md
+++ b/components/omega/doc/devGuide/VertCoord.md
@@ -69,6 +69,33 @@ A list of member variables along with their types and dimension sizes is below:
 | VertCoordMovementWeights | Real | NCellsSize, NVertLayers |
 | RefPseudoThickness | Real | NCellsSize, NVertLayers |
 | BottomGeomDepth | Real | NCellsSize |
+| SurfacePressure | Real | NCellsSize |
+
+### Surface pressure
+
+`SurfacePressure` is the relative pressure (gauge pressure) at the top of the ocean column and
+is the top boundary condition used by `computePressure`. It is owned by `VertCoord` and its host
+mirror is `SurfacePressureH`.
+
+Unlike the other `VertCoord` fields, whose data come from the mesh file (`InitVertCoord` group)
+during construction, `SurfacePressure` is a prognostic quantity read from the initial-condition
+or restart file. `defineFields()` therefore registers it and adds it to the `State` and `Restart`
+field groups (creating those groups if `VertCoord` initializes before `OceanState`). Because
+those streams are read later in `ocnInit`, the data are written directly into the attached device
+array.
+
+`SurfacePressure` is registered as an [optional read](omega-dev-iostreams)
+(`Field::setOptionalRead(true)`), so it is not required to be present in the initial-condition or
+restart file. When the variable is absent, the read is skipped and the array retains the fill
+value written by `attachData`. After the read, `initSurfacePressure()` must be called: it detects
+that leftover fill value on the owned cells and, if found, defaults the whole array to zero, then
+exchanges the halo and copies the device array to the host mirror:
+```c++
+VertCoord::getDefault()->initSurfacePressure(Halo::getDefault());
+```
+Eventually `SurfacePressure` will be updated each timestep via the coupler as a weighted sum of
+atmosphere, sea-ice, and land-ice pressure; that forcing update will live in a separate forcing
+class that writes into this array.
 
 ### Removal
 
@@ -128,3 +155,25 @@ This `parallel_for` iterates over vertical chunks to facilitate vectorization on
 The vector length on GPUs is set to 1 to maximize parallelism.
 The `computeGeopotential` method uses hierarchical parallelism in a very similar way to `computeTargetThickness`, except that it doesn't require a column sum.
 It has an outer `parallel_for` that splits horizontal cells into teams and an inner `parallel_for` that does vertical computations in chunks.
+
+### Layer masking
+
+Field arrays are auto-filled with `FillValueReal` when they are attached to a `Field` (see the [Fill Values design](../design/FillValues.md)).
+Compute kernels only write to active layers, so inactive layers naturally retain the fill value.
+After an initial-condition or restart read, however, state and tracer arrays come from the input file and may hold zeros (or arbitrary values) in inactive and boundary layers.
+`VertCoord` provides a family of helper methods that enforce the correct layer pattern on such arrays.
+All of them use the inclusive `[Min, Max]` active-layer convention consistent with `computeGeomZHeight`/`computePressure`, and use hierarchical parallelism (an outer `parallelForOuter` over the horizontal dimension and an inner `parallelForInner` over the vertical layers).
+
+| Method | Element | Zones applied |
+|--------|---------|---------------|
+| `zeroEdgeField(Arr, NEdgesAll)` | edge | sets `[MinLayerEdgeTop, MaxLayerEdgeBot]` to 0; layers outside retain their fill value |
+| `applyEdgeLayerMask(Arr, NEdgesAll)` | edge | three-zone: outside `[MinLayerEdgeTop, MaxLayerEdgeBot]` → `FillValueReal`; inside that but outside `[MinLayerEdgeBot, MaxLayerEdgeTop]` → 0; active layers unchanged |
+| `applyCellLayerMask(Arr, NCellsAll)` | cell | two-zone: outside `[MinLayerCell, MaxLayerCell]` → `FillValueReal`; active layers unchanged |
+| `applyVertexLayerMask(Arr, NVerticesAll)` | vertex | two-zone: outside `[MinLayerVertexTop, MaxLayerVertexBot]` → `FillValueReal`; active layers unchanged |
+
+Cell fields have no boundary (zero) zone because a cell column is either active or inactive at a given layer, with no neighbor-dependent partial-activity range.
+Vertex fields also have no zeroed boundary zone, but for a different reason: a boundary vertex with one or more active surrounding cells holds valid, generally non-zero data (for example, relative vorticity computed from its active surrounding edges over the full `[MinLayerVertexTop, MaxLayerVertexBot]` range), so the entire valid range is kept rather than zeroed.
+This is unlike a flux-type edge field, where the normal flux through a face bordering land genuinely vanishes.
+The flux-type edge fields handled by `zeroEdgeField` (for example `NormalVelocityTend` and `VelocityDel2Aux.Del2Edge`) are zeroed before being recomputed each time step so that boundary edges carry 0 rather than `FillValueReal`.
+The `apply*LayerMask` methods are driven through `OceanState::applyLayerMasks`, which is called once from `ocnInit` after the IC/restart read: it applies the edge mask to `NormalVelocity` and the cell mask to `PseudoThickness`, `Temperature`, and `Salinity`.
+`applyVertexLayerMask` has no production caller yet because no vertex-based state field is read from the IC/restart; it is provided for completeness and verified by the fill-value CTest.
diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md
index 8d1b38501df8..eb8504db8047 100644
--- a/components/omega/doc/index.md
+++ b/components/omega/doc/index.md
@@ -115,6 +115,7 @@ design/Decomp
 design/Driver
 design/EOS
 design/Error
+design/FillValues
 design/Halo
 design/HorzMeshClass
 design/Logging
diff --git a/components/omega/doc/userGuide/OceanState.md b/components/omega/doc/userGuide/OceanState.md
index b13404587134..53bc70e6d404 100644
--- a/components/omega/doc/userGuide/OceanState.md
+++ b/components/omega/doc/userGuide/OceanState.md
@@ -2,7 +2,16 @@
 
 ## Ocean State
 
-The `OceanState` class provides a container for the non-tracer prognostic variables in Omega, namely `NormalVelocity` and `PseudoThickness`.
-Upon creation of a `OceanState` instance, these variables are allocated and registered with the IO infrastructure.
+The `OceanState` class provides a container for the non-tracer prognostic variables in Omega:
+`NormalVelocity` and `PseudoThickness`.
+Upon creation of an `OceanState` instance, these variables are allocated and registered with the
+IO infrastructure as part of the `State` and `Restart` field groups, so they are read from the
+initial-condition file and written to restart files.
 The class contains a method to update the time levels for the state variables between timesteps.
 This involves a halo update, time level index update, and updating the `IOFields` data references.
+
+The relative pressure at the top of the ocean column, `SurfacePressure`, is owned by the
+[vertical coordinate](omega-user-vert-coord) rather than the `OceanState`, since it is the top
+boundary condition of the pressure field. It is still registered in the `State` and `Restart`
+field groups, so it continues to be written to restart files and read from the initial-condition
+and restart files when present; if it is absent, it defaults to zero.
diff --git a/components/omega/doc/userGuide/Tracers.md b/components/omega/doc/userGuide/Tracers.md
index d3fefb14af1e..18ba1c93a430 100644
--- a/components/omega/doc/userGuide/Tracers.md
+++ b/components/omega/doc/userGuide/Tracers.md
@@ -26,17 +26,23 @@ static void defineAllTracers() {
           "sea_water_potential_temperature", ///< [in] CF standard Name
           -273.15,                           ///< [in] min valid field value
           100.0,                             ///< [in] max valid field value
-          1.e33,                             ///< [in] value for undef entries
+          FillValueReal,                     ///< [in] value for undef entries
           IndxTemp);                         ///< [out] (optional) static index
 
-   define("Salt", "Salinity", "psu", "sea_water_salinity", 0.0, 50.0, 1.e33,
-          IndxSalt);
-   define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, 1.e33);
-   define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, 1.e33);
-   define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, 1.e33);
+   define("Salt", "Salinity", "psu", "sea_water_salinity", 0.0, 50.0,
+          FillValueReal, IndxSalt);
+   define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, FillValueReal);
+   define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, FillValueReal);
+   define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, FillValueReal);
 }
 ```
 
+The `FillValueReal` constant is defined in `FillValues.h` and matches the
+NetCDF-C standard fill value (`NC_FILL_DOUBLE` or `NC_FILL_FLOAT` depending on
+the build precision). It is automatically in scope in `TracerDefs.inc` via
+`Field.h`. New tracers should always use `FillValueReal` (or `FillValueI4` /
+`FillValueI8` for integer fields) rather than hardcoded literals.
+
 To add a new tracer, simply call the `define` function with the appropriate
 arguments. Index argument is optional one that allows to access the tracer
 data using the given tracer index variable.
diff --git a/components/omega/doc/userGuide/VertCoord.md b/components/omega/doc/userGuide/VertCoord.md
index 21d197efd6c6..80fa99f0da08 100644
--- a/components/omega/doc/userGuide/VertCoord.md
+++ b/components/omega/doc/userGuide/VertCoord.md
@@ -4,7 +4,9 @@
 
 ### Overview
 
-Omega uses pseudo height, $\tilde{z} = \frac{p}{\rho_0 g}$,  as the vertical coordinate ([V1 governing equation document](omega-design-governing-eqns-omega1)).
+Omega uses pseudo height, $\tilde{z} = \frac{p}{\rho_0 g}$,  as the vertical coordinate ([V1 governing equation document](omega-design-governing-eqns-omega1)),
+where $p$ is **relative pressure** (gauge pressure, i.e., absolute pressure minus the standard atmosphere of 101325 Pa).
+With this definition, $\tilde{z} = 0$ at the sea surface under standard atmospheric forcing.
 In practice, $\tilde{z}$ is not computed directly.
 Instead, the model tracks and evolves the pseudo-thickness, $\tilde{h} = \Delta \tilde{z}$, defined as the difference between adjacent layer interfaces.
 The pseudo height is essentially a normalized pressure coordinate, with the advantage that it has units of meters.
@@ -27,8 +29,8 @@ Multiple instances of the vertical coordinate class can be created and accessed
 | ------------- | ----------- | ----- |
 | NVertLayers   | maximum number of vertical layers | - |
 | NVertLayersP1 | maximum number of vertical layers plus 1 | - |
-| PressureInterface | pressure at layer interfaces | force per unit area at layer interfaces | kg m$^{-1}$ s$^{-2}$ |
-| PressureMid | pressure at layer mid points | force per unit area at layer mid point | kg m$^{-1}$ s$^{-2}$ |
+| PressureInterface | relative pressure (gauge pressure) at layer interfaces | Pa |
+| PressureMid | relative pressure (gauge pressure) at layer midpoints | Pa |
 | GeomZInterface | geometric height of layer interfaces | m |
 | GeomZMid | geometric height of layer midpoint | m |
 | GeopotentialMid | geopotential at layer mid points | m$^2$/s$^2$|
@@ -46,3 +48,29 @@ Multiple instances of the vertical coordinate class can be created and accessed
 | VertCoordMovementWeights | weights to specify how total column thickness changes are distributed across layers | - |
 | RefPseudoThickness | reference pseudo-thickness used to distribute total column thickness changes | m |
 | BottomGeomDepth | positive down distance from the reference geoid to the bottom | m |
+| SurfacePressure | relative pressure (gauge pressure) at the top of the ocean column | Pa |
+
+### Surface pressure
+
+`SurfacePressure` is the relative pressure (gauge pressure, i.e., absolute pressure minus the
+standard atmosphere of 101325 Pa) at the top of the ocean column. It is the top boundary
+condition used when integrating layer pressures downward in `computePressure`, and is
+effectively $\tilde{z}$ at the top of the water column (times a scaling factor).
+Unlike the other pressure variables, which are recomputed each timestep, `SurfacePressure` is
+registered in the `State` and `Restart` field groups, so it is written to restart files and read
+from the initial-condition and restart files when present. The read is optional: if
+`SurfacePressure` is not present in the file, it defaults to zero rather than causing an error.
+For a free-surface ocean run with standard atmospheric forcing, the surface pressure is
+zero (0 Pa), so it may simply be omitted from the initial-condition file.
+For a coupled ocean–atmosphere run, `SurfacePressure` holds the anomaly of the atmospheric
+surface pressure from the standard atmosphere; eventually it will be updated each timestep
+via the coupler as a weighted sum of atmosphere, sea-ice, and land-ice pressure.
+
+### Inactive and boundary layers in output
+
+Each output field declares a fill value (matching the NetCDF-C standard `_FillValue`) that marks entries outside the active domain.
+In a cell column only the active layers (`[MinLayerCell, MaxLayerCell]`) carry valid data; deeper, inactive layers contain the fill value.
+Edge and vertex fields have an intermediate boundary range where some, but not all, of the surrounding cells are active.
+For flux-type edge fields (such as `NormalVelocity`), the flux through a face bordering land vanishes, so those boundary layers carry `0`; layers fully outside the valid range carry the fill value.
+Vertex fields, by contrast, keep their computed (generally non-zero) values throughout the valid range, since a boundary vertex still has active surrounding cells contributing to quantities like vorticity; only layers outside the valid range carry the fill value.
+After an initial-condition or restart read, `VertCoord` enforces this pattern on the state and tracer fields so that analysis tools (ncview, Xarray, NCO, etc.) recognize and mask the inactive entries automatically.
diff --git a/components/omega/src/base/FillValues.h b/components/omega/src/base/FillValues.h
new file mode 100644
index 000000000000..21b6caf9ed1b
--- /dev/null
+++ b/components/omega/src/base/FillValues.h
@@ -0,0 +1,41 @@
+#ifndef OMEGA_FILL_VALUES_H
+#define OMEGA_FILL_VALUES_H
+//===-- base/FillValues.h - standard fill value constants -------*- C++ -*-===//
+///
+/// \file
+/// \brief Standard fill value constants matching NetCDF-C NC_FILL_* values.
+///
+/// These values exactly match the SCORPIO PIO_FILL_* constants (which wrap
+/// NC_FILL_* from netcdf.h). Defined here as literals so this header can be
+/// included without pulling in pio.h or mpi.h, which have strict include-order
+/// requirements that make them unsafe to include from generic headers.
+//
+//===----------------------------------------------------------------------===//
+
+#include "DataTypes.h"
+
+namespace OMEGA {
+
+/// Type-indexed fill value. Instantiating FillValue<T> for an unsupported
+/// type T produces a link error (primary template is intentionally undefined).
+template <typename T> constexpr T FillValue;
+
+template <> inline constexpr I4 FillValue<I4> = -2147483647; ///< NC_FILL_INT
+template <>
+inline constexpr I8 FillValue<I8> = -9223372036854775806LL; ///< NC_FILL_INT64
+template <>
+inline constexpr R4 FillValue<R4> = 9.9692099683868690e+36f; ///< NC_FILL_FLOAT
+template <>
+inline constexpr R8 FillValue<R8> = 9.9692099683868690e+36; ///< NC_FILL_DOUBLE
+
+/// Named aliases for readability in comparison and test code.
+constexpr I4 FillValueI4     = FillValue<I4>;
+constexpr I8 FillValueI8     = FillValue<I8>;
+constexpr R4 FillValueR4     = FillValue<R4>;
+constexpr R8 FillValueR8     = FillValue<R8>;
+constexpr Real FillValueReal = FillValue<Real>;
+
+} // end namespace OMEGA
+
+//===----------------------------------------------------------------------===//
+#endif // OMEGA_FILL_VALUES_H
diff --git a/components/omega/src/base/IO.h b/components/omega/src/base/IO.h
index a8ec08e73f4a..4aa606558500 100644
--- a/components/omega/src/base/IO.h
+++ b/components/omega/src/base/IO.h
@@ -364,6 +364,7 @@ void writeNDVar(void *Variable, ///< [in] variable to be written
 );
 
 } // end namespace IO
+
 } // end namespace OMEGA
 
 //===----------------------------------------------------------------------===//
diff --git a/components/omega/src/infra/Field.cpp b/components/omega/src/infra/Field.cpp
index 5c24a9efd01d..ccf2d9d7984b 100644
--- a/components/omega/src/infra/Field.cpp
+++ b/components/omega/src/infra/Field.cpp
@@ -54,8 +54,8 @@ void Field::init(const Clock *ModelClock // [in] default model clock
    std::string CalName      = CalendarCFName[CalKind];
    std::vector<std::string> DimNamesTmp; // empty dim names vector
    std::shared_ptr<Field> TimeField =
-       create("time", "time", UnitString, "time", 0.0, 1.e20, -9.99e30, 0,
-              DimNamesTmp, true, true);
+       create("time", "time", UnitString, "time", 0.0, 1.e20, 0, DimNamesTmp,
+              true, true);
    TimeField->addMetadata("calendar", CalName);
 }
 
@@ -85,8 +85,7 @@ Field::create(const std::string &FieldName,   // [in] Name of variable/field
               const std::string &StdName,     // [in] CF standard Name
               const std::any ValidMin,        // [in] min valid field value
               const std::any ValidMax,        // [in] max valid field value
-              const std::any FillValue, // [in] scalar for undefined entries
-              const int NumDims,        // [in] number of dimensions
+              const int NumDims,              // [in] number of dimensions
               const std::vector<std::string> &Dimensions, // [in] dim names
               bool TimeDependent,  // [in] flag for time dependent field
               bool RetainPrecision // [in] flag to retain full prec in IO
@@ -118,8 +117,6 @@ Field::create(const std::string &FieldName,   // [in] Name of variable/field
    ThisField->FieldMeta["valid_min"]     = ValidMin;
    ThisField->FieldMeta["ValidMax"]      = ValidMax;
    ThisField->FieldMeta["valid_max"]     = ValidMax;
-   ThisField->FieldMeta["FillValue"]     = FillValue;
-   ThisField->FieldMeta["_FillValue"]    = FillValue;
 
    // Set the time-dependent flag
    ThisField->TimeDependent = TimeDependent;
@@ -128,6 +125,11 @@ Field::create(const std::string &FieldName,   // [in] Name of variable/field
    // if reduced precision for the stream/file is requested
    ThisField->RetainPrecision = RetainPrecision;
 
+   // Fields are required in input files by default; owning modules may mark a
+   // field optional via setOptionalRead if it should default to its fill value
+   // when absent from a read stream.
+   ThisField->OptionalRead = false;
+
    // Number of dimensions for the field
    ThisField->NDims = NumDims;
 
@@ -336,6 +338,16 @@ bool Field::isDistributed() const { return Distributed; }
 // when reduced precision for the stream/file is requested
 bool Field::retainPrecision() const { return RetainPrecision; }
 
+//------------------------------------------------------------------------------
+// Determines whether the field may be absent from an input file. When true, a
+// read stream that does not contain the field's variable skips it and leaves
+// the attached array at its fill value rather than failing.
+bool Field::isOptionalRead() const { return OptionalRead; }
+
+//------------------------------------------------------------------------------
+// Sets whether the field may be absent from an input file (see isOptionalRead)
+void Field::setOptionalRead(bool OptRead) { OptionalRead = OptRead; }
+
 //------------------------------------------------------------------------------
 // Returns a vector of dimension names associated with each dimension
 // of an array field. Returns an error code.
diff --git a/components/omega/src/infra/Field.h b/components/omega/src/infra/Field.h
index 313883184edc..52a3b9a315a9 100644
--- a/components/omega/src/infra/Field.h
+++ b/components/omega/src/infra/Field.h
@@ -19,6 +19,7 @@
 #include "DataTypes.h"
 #include "Dimension.h"
 #include "Error.h"
+#include "FillValues.h"
 #include "Logging.h"
 #include "OmegaKokkos.h"
 #include "TimeMgr.h"
@@ -26,6 +27,7 @@
 #include <map>
 #include <memory>
 #include <set>
+#include <type_traits>
 
 namespace OMEGA {
 
@@ -77,11 +79,30 @@ class Field {
    /// operations when reduced precision is requested
    bool RetainPrecision;
 
+   /// Flag for whether this field may be absent from an input file. When true,
+   /// a read stream that does not contain this field's variable skips it
+   /// (leaving the attached array at its fill value) instead of failing.
+   bool OptionalRead;
+
    /// Data attached to this field. This will be a pointer to the Kokkos
    /// array holding the data. We use a void pointer to manage all the
    /// various types and cast to the appropriate type when needed.
    std::shared_ptr<void> DataArray;
 
+   /// Fills every element of InDataArray with the standard fill value for the
+   /// array's element type, then records that value in the field metadata.
+   /// Called automatically by attachData() so inactive entries are initialized
+   /// to a well-defined sentinel before any compute routine runs.
+   template <typename T> void fillWithValue(const T &InDataArray) {
+      using ValType             = typename T::non_const_value_type;
+      constexpr ValType FillVal = FillValue<ValType>;
+      // _FillValue is the NetCDF/CF standard attribute recognized by analysis
+      // tools (ncview, Xarray, NCO, ...); store only that, no non-standard
+      // alias
+      FieldMeta["_FillValue"] = FillVal;
+      deepCopy(InDataArray, FillVal);
+   }
+
  public:
    //---------------------------------------------------------------------------
    // Initialization
@@ -112,7 +133,6 @@ class Field {
           const std::string &StdName,     ///< [in] CF standard Name
           const std::any ValidMin,        ///< [in] min valid field value
           const std::any ValidMax,        ///< [in] max valid field value
-          const std::any FillValue,       ///< [in] scalar for undefined entries
           const int NumDims,              ///< [in] number of dimensions
           const std::vector<std::string> &Dimensions, ///< [in] dim names
           const bool TimeDependent   = true, ///< [in] opt flag for unlim time
@@ -211,6 +231,16 @@ class Field {
    /// a reduced precision stream is requested.
    bool retainPrecision() const;
 
+   /// Determine whether this field may be absent from an input file. When true,
+   /// a read stream that does not contain the field's variable skips it and
+   /// leaves the attached array at its fill value rather than failing.
+   bool isOptionalRead() const;
+
+   /// Set whether this field may be absent from an input file (see
+   /// isOptionalRead). Fields are required (not optional) by default.
+   void setOptionalRead(bool OptRead ///< [in] optional-read flag value
+   );
+
    //---------------------------------------------------------------------------
    // Metadata functions
    //---------------------------------------------------------------------------
@@ -303,7 +333,8 @@ class Field {
    /// templated based on the array data type so a template argument with
    /// the proper array data type (eg <Array2DR4>) must be supplied.
    template <typename T>
-   void attachData(const T &InDataArray ///< [in] Array with data to attach
+   void attachData(const T &InDataArray,    ///< [in] Array with data to attach
+                   bool FillOnAttach = true ///< [in] fill array with fill value
    ) {
       OMEGA_ASSERT(isKokkosArray<T>,
                    "Field::attachData requires Kokkos array as input");
@@ -314,6 +345,13 @@ class Field {
       // Determine type and location
       DataType = checkArrayType<T>();
       MemLoc   = findArrayMemLoc<T>();
+
+      // Initialize every element to the declared fill value so inactive
+      // entries (e.g. layers below MaxLayerCell) are well-defined in output.
+      // Pass FillOnAttach=false when re-attaching an existing data-filled array
+      // (e.g. time-level pointer updates) to preserve the computed values.
+      if (FillOnAttach)
+         fillWithValue<T>(InDataArray);
    };
 
    //---------------------------------------------------------------------------
@@ -325,7 +363,8 @@ class Field {
    template <typename T>
    static void
    attachFieldData(const std::string &FieldName, ///< [in] Name of Field
-                   const T &InDataArray ///< [in] Array with data to attach
+                   const T &InDataArray,    ///< [in] Array with data to attach
+                   bool FillOnAttach = true ///< [in] fill array with fill value
    ) {
       OMEGA_ASSERT(isKokkosArray<T>,
                    "Field::attachFieldData requires Kokkos array as input");
@@ -335,7 +374,7 @@ class Field {
          // Retrieve the field
          auto ThisField = AllFields[FieldName];
          // Attach the data array
-         ThisField->attachData<T>(InDataArray);
+         ThisField->attachData<T>(InDataArray, FillOnAttach);
 
       } else { // field has not yet been defined
          ABORT_ERROR("Field: error attaching data to {}. Field not defined",
diff --git a/components/omega/src/infra/IOStream.cpp b/components/omega/src/infra/IOStream.cpp
index 95e082480ae9..cd7444e853e9 100644
--- a/components/omega/src/infra/IOStream.cpp
+++ b/components/omega/src/infra/IOStream.cpp
@@ -1019,8 +1019,8 @@ void IOStream::writeFieldData(
       DataI4.resize(LocSize);
       DataPtr = DataI4.data();
       // get fill value
-      Err += FieldPtr->getMetadata("FillValue", FillValI4);
-      CHECK_ERROR_ABORT(Err, "Error retrieving FillValue for Field {}",
+      Err += FieldPtr->getMetadata("_FillValue", FillValI4);
+      CHECK_ERROR_ABORT(Err, "Error retrieving _FillValue for Field {}",
                         FieldName);
       FillValPtr = &FillValI4;
 
@@ -1161,8 +1161,8 @@ void IOStream::writeFieldData(
       DataI8.resize(LocSize);
       DataPtr = DataI8.data();
       // Get fill value
-      Err += FieldPtr->getMetadata("FillValue", FillValI8);
-      CHECK_ERROR_ABORT(Err, "Error retrieving FillValue for Field {}",
+      Err += FieldPtr->getMetadata("_FillValue", FillValI8);
+      CHECK_ERROR_ABORT(Err, "Error retrieving _FillValue for Field {}",
                         FieldName);
       FillValPtr = &FillValI8;
 
@@ -1302,8 +1302,8 @@ void IOStream::writeFieldData(
       DataR4.resize(LocSize);
       DataPtr = DataR4.data();
       // Get fill value
-      Err += FieldPtr->getMetadata("FillValue", FillValR4);
-      CHECK_ERROR_ABORT(Err, "Error retrieving FillValue for Field {}",
+      Err += FieldPtr->getMetadata("_FillValue", FillValR4);
+      CHECK_ERROR_ABORT(Err, "Error retrieving _FillValue for Field {}",
                         FieldName);
       FillValPtr = &FillValR4;
 
@@ -1441,8 +1441,8 @@ void IOStream::writeFieldData(
    case ArrayDataType::R8:
 
       // Get fill value
-      Err += FieldPtr->getMetadata("FillValue", FillValR8);
-      CHECK_ERROR_ABORT(Err, "Error retrieving FillValue for Field {}",
+      Err += FieldPtr->getMetadata("_FillValue", FillValR8);
+      CHECK_ERROR_ABORT(Err, "Error retrieving _FillValue for Field {}",
                         FieldName);
       DataR8.resize(LocSize);
       if (ReducePrecision and !RetainPrecision) {
@@ -1716,11 +1716,22 @@ Error IOStream::readFieldData(
       } else {
          Err = IO::readNDVar(DataPtr, OldFieldName, FileID, FieldID);
       }
-      if (Err.isFail()) // still cannot find field, return with error
+      if (Err.isFail()) {
+         // If the field is optional, a missing variable is not an error. Leave
+         // the attached array at the fill value set by Field::attachData and
+         // return success so the remaining fields in the stream still read.
+         if (FieldPtr->isOptionalRead()) {
+            LOG_INFO("IOStream::readFieldData: optional field {} not found in "
+                     "stream {}; leaving fill value",
+                     FieldName, Name);
+            return Error{}; // success
+         }
+         // still cannot find field, return with error
          RETURN_ERROR(
              Err, ErrorCode::Fail,
              "IOStream::readFieldData: Field {} or {} not found in stream {}",
              FieldName, OldFieldName, Name);
+      }
    }
 
    // Unpack vector into array based on type, dims and location
@@ -2482,13 +2493,15 @@ void IOStream::writeStream(
    if (Multiframe and !FileAlarm.isRinging())
       FileTime = *(FileAlarm.getRingTimePrev());
 
-   // Update the time field with elapsed time since simulation start
+   // Update the time field with elapsed time since simulation start.
+   // Attach before setting the value so the fill-on-attach does not
+   // overwrite the elapsed time (Kokkos host views share backing memory).
    TimeInterval ElapsedTime = SimTime - StartTime;
    R8 ElapsedTimeR8;
    ElapsedTime.get(ElapsedTimeR8, TimeUnits::Seconds);
    HostArray1DR8 OutTime("OutTime", 1);
-   OutTime(0) = ElapsedTimeR8;
    Field::attachFieldData("time", OutTime);
+   OutTime(0) = ElapsedTimeR8;
 
    // Reset alarms and flags
    if (OnStartup)
diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp
index c2705a3b9cc6..7c0bc3d71ca5 100644
--- a/components/omega/src/ocn/AuxiliaryState.cpp
+++ b/components/omega/src/ocn/AuxiliaryState.cpp
@@ -117,10 +117,8 @@ void AuxiliaryState::computeMomAux(const OceanState *State,
 
    OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell);
    OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell);
-   OMEGA_SCOPE(MinLayerVertexBot, VCoord->MinLayerVertexBot);
    OMEGA_SCOPE(MinLayerVertexTop, VCoord->MinLayerVertexTop);
    OMEGA_SCOPE(MaxLayerVertexBot, VCoord->MaxLayerVertexBot);
-   OMEGA_SCOPE(MaxLayerVertexTop, VCoord->MaxLayerVertexTop);
    OMEGA_SCOPE(MinLayerEdgeTop, VCoord->MinLayerEdgeTop);
    OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot);
    OMEGA_SCOPE(MaxLayerEdgeBot, VCoord->MaxLayerEdgeBot);
@@ -169,6 +167,12 @@ void AuxiliaryState::computeMomAux(const OceanState *State,
    const auto &VelocityDivCell = KineticAux.VelocityDivCell;
    const auto &RelVortVertex   = VorticityAux.RelVortVertex;
 
+   // Zero Del2Edge at boundary layers before accumulating over neighbouring
+   // edges in computeVarsOnCell/Vertex. Del4 hyperdiffusion reads Del2Edge at
+   // all edges sharing a cell or vertex; fill values in boundary layers would
+   // corrupt Del2DivCell and Del2RelVortVertex.
+   VCoord->zeroEdgeField(VelocityDel2Aux.Del2Edge, Mesh->NEdgesAll);
+
    Pacer::start("AuxState:edgeAuxState2", 2);
    parallelForOuter(
        "edgeAuxState2", {Mesh->NEdgesAll},
@@ -204,8 +208,12 @@ void AuxiliaryState::computeMomAux(const OceanState *State,
    parallelForOuter(
        "vertexAuxState2", {Mesh->NVerticesAll},
        KOKKOS_LAMBDA(int IVertex, const TeamMember &Team) {
-          const int KMin   = MinLayerVertexBot(IVertex);
-          const int KMax   = MaxLayerVertexTop(IVertex);
+          // Del2RelVortVertex is computed over the full vertex valid range
+          // [MinLayerVertexTop, MaxLayerVertexBot] so that boundary-vertex
+          // layers read by the biharmonic velocity tendency are valid rather
+          // than fill values (see VelocityDel2AuxVars::computeVarsOnVertex).
+          const int KMin   = MinLayerVertexTop(IVertex);
+          const int KMax   = MaxLayerVertexBot(IVertex);
           const int KRange = vertRangeChunked(KMin, KMax);
 
           parallelForInner(
diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp
index 8101e6bab42f..16ea4da9f5a2 100644
--- a/components/omega/src/ocn/Eos.cpp
+++ b/components/omega/src/ocn/Eos.cpp
@@ -49,9 +49,6 @@ Eos::Eos(const std::string &Name, ///< [in] Name for eos object
    BruntVaisalaFreqSq = Array2DReal("BruntVaisalaFreqSq", Mesh->NCellsSize,
                                     VCoord->NVertLayersP1);
 
-   deepCopy(SpecVol, 1.0_Real / RhoSw);
-   deepCopy(SpecVolDisplaced, 1.0_Real / RhoSw);
-
    defineFields();
 }
 
@@ -142,8 +139,6 @@ void Eos::computeSpecVol(const Array2DReal &ConservTemp,
    OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell);
    OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell);
 
-   deepCopy(LocSpecVol, 0); /// Initialize local specific volume to zero
-
    I4 KDisp = 0; /// No displacement in this case
 
    /// Dispatch to the correct EOS calculation
@@ -208,9 +203,6 @@ void Eos::computeSpecVolDisp(const Array2DReal &ConservTemp,
    OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell);
    OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell);
 
-   deepCopy(LocSpecVolDisplaced,
-            0); /// Initialize local specific volume to zero
-
    /// Dispatch to the correct EOS calculation
    /// If EosChoice is Linear, the displaced specific
    /// volume is the same as the specific volume
@@ -273,9 +265,6 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp,
    OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell);
    OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell);
 
-   deepCopy(LocBruntVaisalaFreqSq,
-            0); /// Initialize local Brunt-Vaisala frequency to zero
-
    /// Dispatch to the correct squared Brunt-Vaisala frequency calculation
    if (EosChoice == EosType::LinearEos) {
       /// If Linear EOS, use linear squared Brunt-Vaisala frequency calculation
@@ -368,8 +357,7 @@ void Eos::defineFields() {
    }
 
    /// Create fields for state variables
-   const Real FillValue = -9.99e30;
-   int NDims            = 2;
+   int NDims = 2;
    std::vector<std::string> DimNames(NDims);
    DimNames[0] = "NCells";
    DimNames[1] = "NVertLayers";
@@ -382,9 +370,8 @@ void Eos::defineFields() {
                      "sea_water_specific_volume",      // CF-ish Name
                      0.0,                              // Min valid value
                      std::numeric_limits<Real>::max(), // Max valid value
-                     FillValue, // Scalar used for undefined entries
-                     NDims,     // Number of dimensions
-                     DimNames   // Dimension names
+                     NDims,                            // Number of dimensions
+                     DimNames                          // Dimension names
        );
    /// Create and register the displaced specific volume field
    auto SpecVolDisplacedField =
@@ -395,9 +382,8 @@ void Eos::defineFields() {
                      "sea_water_specific_volume_displaced", // CF-ish Name
                      0.0,                                   // Min valid value
                      std::numeric_limits<Real>::max(),      // Max valid value
-                     FillValue, // Scalar used for undefined entried
-                     NDims,     // Number of dimensions
-                     DimNames   // Dimension names
+                     NDims,   // Number of dimensions
+                     DimNames // Dimension names
        );
 
    // Brunt-Vaisala frequency is located at interfaces
@@ -411,9 +397,8 @@ void Eos::defineFields() {
                      "sea_water_brunt_vaisala_frequency_squared", // CF-ish Name
                      std::numeric_limits<Real>::min(), // Min valid value
                      std::numeric_limits<Real>::max(), // Max valid value
-                     FillValue, // Scalar used for undefined entries
-                     NDims,     // Number of dimensions
-                     DimNames   // Dimension names
+                     NDims,                            // Number of dimensions
+                     DimNames                          // Dimension names
        );
 
    // Create a field group for the eos-specific state fields
diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h
index a530549d493d..5e4fd89cae21 100644
--- a/components/omega/src/ocn/Eos.h
+++ b/components/omega/src/ocn/Eos.h
@@ -36,8 +36,10 @@ class Teos10Eos {
 
    //   The functor takes the full arrays of specific volume (inout),
    //   the indices ICell and KChunk, and the ocean tracers (conservative)
-   //   temperature, and (absolute) salinity as inputs, and outputs the
-   //   specific volume according to the Roquet et al. 2015 75 term expansion.
+   //   temperature, (absolute) salinity, and relative pressure (gauge pressure,
+   //   i.e., absolute pressure minus the standard atmosphere) as inputs, and
+   //   outputs the specific volume according to the Roquet et al. 2015 75 term
+   //   expansion.
    KOKKOS_FUNCTION void operator()(Array2DReal SpecVol, I4 ICell, I4 KChunk,
                                    const Array2DReal &ConservTemp,
                                    const Array2DReal &AbsSalinity,
@@ -353,7 +355,9 @@ class Teos10Eos {
    }
 
    /// Calculates freezing Conservative Temperature using TEOS-10 polynomial
-   /// (polynomial error in [-5e-4, 6e-4] K, from GSW package)
+   /// (polynomial error in [-5e-4, 6e-4] K, from GSW package).
+   /// P is relative pressure (gauge pressure in Pa, i.e., absolute pressure
+   /// minus the standard atmosphere).
    KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P,
                                        const Real SaturationFract) const {
       constexpr Real Sso = 35.16504;
@@ -482,8 +486,9 @@ class Teos10BruntVaisalaFreqSq {
 
    //   The functor takes the full arrays of squared Brunt-Vaisala frequency
    //   (inout) the index ICell, and the ocean tracers (conservative)
-   //   temperature, (absolute) salinity, pressure, specific volume as inputs,
-   //   and outputs the squared Brunt-Vaisala frequency.
+   //   temperature, (absolute) salinity, relative pressure (gauge pressure,
+   //   i.e., absolute pressure minus the standard atmosphere), and specific
+   //   volume as inputs, and outputs the squared Brunt-Vaisala frequency.
    KOKKOS_FUNCTION void operator()(Array2DReal BruntVaisalaFreqSq, I4 ICell,
                                    I4 KChunk, const Array2DReal &ConservTemp,
                                    const Array2DReal &AbsSalinity,
diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp
index 9e47b81e350a..bb8b909aad42 100644
--- a/components/omega/src/ocn/HorzMesh.cpp
+++ b/components/omega/src/ocn/HorzMesh.cpp
@@ -559,10 +559,9 @@ void HorzMesh::defineMeshFields() {
                      "x",                                 // CF standard Name
                      0.0,                                 // min valid value
                      7.0E+6,                              // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,                               // num of dimensions
+                     DimNames,                            // dimension names
+                     false                                // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, XCell);
@@ -576,10 +575,9 @@ void HorzMesh::defineMeshFields() {
                      "y",                                 // CF standard Name
                      0.0,                                 // min valid value
                      7.0E+6,                              // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,                               // num of dimensions
+                     DimNames,                            // dimension names
+                     false                                // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, YCell);
@@ -593,10 +591,9 @@ void HorzMesh::defineMeshFields() {
                      "z",                                 // CF standard Name
                      0.0,                                 // min valid value
                      7.0E+6,                              // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,                               // num of dimensions
+                     DimNames,                            // dimension names
+                     false                                // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, ZCell);
@@ -610,10 +607,9 @@ void HorzMesh::defineMeshFields() {
                      "latitude",                             // CF standard Name
                      -3.1415927,                             // min valid value
                      3.1415927,                              // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,    // num of dimensions
+                     DimNames, // dimension names
+                     false     // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, LatCell);
@@ -627,7 +623,6 @@ void HorzMesh::defineMeshFields() {
                      "longitude", // CF standard Name
                      0.0,         // min valid value
                      6.28319,     // max valid value
-                     -9.99E+30,   // scalar for undefined entries
                      NDims,       // num of dimensions
                      DimNames,    // dimension names
                      false        // not time dependent
@@ -646,10 +641,9 @@ void HorzMesh::defineMeshFields() {
                      "x",                               // CF standard Name
                      0.0,                               // min valid value
                      7.0E+6,                            // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,                             // num of dimensions
+                     DimNames,                          // dimension names
+                     false                              // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, XEdge);
@@ -663,10 +657,9 @@ void HorzMesh::defineMeshFields() {
                      "y",                               // CF standard Name
                      0.0,                               // min valid value
                      7.0E+6,                            // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,                             // num of dimensions
+                     DimNames,                          // dimension names
+                     false                              // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, YEdge);
@@ -680,10 +673,9 @@ void HorzMesh::defineMeshFields() {
                      "z",                               // CF standard Name
                      0.0,                               // min valid value
                      7.0E+6,                            // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,                             // num of dimensions
+                     DimNames,                          // dimension names
+                     false                              // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, ZEdge);
@@ -697,10 +689,9 @@ void HorzMesh::defineMeshFields() {
                      "latitude",                           // CF standard Name
                      -3.1415927,                           // min valid value
                      3.1415927,                            // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,                                // num of dimensions
+                     DimNames,                             // dimension names
+                     false                                 // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, LatEdge);
@@ -714,10 +705,9 @@ void HorzMesh::defineMeshFields() {
                      "longitude",                           // CF standard Name
                      0.0,                                   // min valid value
                      6.28319,                               // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,                                 // num of dimensions
+                     DimNames,                              // dimension names
+                     false // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, LonEdge);
@@ -733,10 +723,9 @@ void HorzMesh::defineMeshFields() {
                      "x",                                  // CF standard Name
                      0.0,                                  // min valid value
                      7.0E+6,                               // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,                                // num of dimensions
+                     DimNames,                             // dimension names
+                     false                                 // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, XVertex);
@@ -750,10 +739,9 @@ void HorzMesh::defineMeshFields() {
                      "y",                                  // CF standard Name
                      0.0,                                  // min valid value
                      7.0E+6,                               // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,                                // num of dimensions
+                     DimNames,                             // dimension names
+                     false                                 // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, YVertex);
@@ -767,10 +755,9 @@ void HorzMesh::defineMeshFields() {
                      "z",                                  // CF standard Name
                      0.0,                                  // min valid value
                      7.0E+6,                               // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // num of dimensions
-                     DimNames,  // dimension names
-                     false      // not time dependent
+                     NDims,                                // num of dimensions
+                     DimNames,                             // dimension names
+                     false                                 // not time dependent
        );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, ZVertex);
@@ -784,7 +771,6 @@ void HorzMesh::defineMeshFields() {
                      "latitude", // CF standard Name
                      -3.1415927, // min valid value
                      3.1415927,  // max valid value
-                     -9.99E+30,  // scalar for undefined entries
                      NDims,      // num of dimensions
                      DimNames,   // dimension names
                      false       // not time dependent
@@ -801,7 +787,6 @@ void HorzMesh::defineMeshFields() {
                      "longitude", // CF standard Name
                      0.0,         // min valid value
                      6.28319,     // max valid value
-                     -9.99E+30,   // scalar for undefined entries
                      NDims,       // num of dimensions
                      DimNames,    // dimension names
                      false        // not time dependent
@@ -820,10 +805,9 @@ void HorzMesh::defineMeshFields() {
                                       "cell_area", // CF standard name
                                       0.0,         // min valid value
                                       9.99E+30,    // max valid value
-                                      -9.99E+30, // scalar for undefined entries
-                                      NDims,     // num of dimensions
-                                      DimNames,  // dimension names
-                                      false      // not time dependent
+                                      NDims,       // num of dimensions
+                                      DimNames,    // dimension names
+                                      false        // not time dependent
    );
    MeshGroupIn->addField(FieldName);
    Field::attachFieldData<Array1DReal>(FieldName, AreaCell);
@@ -838,7 +822,6 @@ void HorzMesh::defineMeshFields() {
                      "",              // CF standard name
                      0.0,             // min valid value
                      9.99E+30,        // max valid value
-                     -9.99E+30,       // scalar for undefined entries
                      NDims,           // num of dimensions
                      DimNames,        // dimension names
                      false            // not time dependent
@@ -856,7 +839,6 @@ void HorzMesh::defineMeshFields() {
                      "cell_area", // CF standard name
                      0.0,         // min valid value
                      9.99E+30,    // max valid value
-                     -9.99E+30,   // scalar for undefined entries
                      NDims,       // num of dimensions
                      DimNames,    // dimension names
                      false        // not time dependent
@@ -875,7 +857,6 @@ void HorzMesh::defineMeshFields() {
                      "",                    // CF standard name
                      0.0,                   // min valid value
                      9.99E+30,              // max valid value
-                     -9.99E+30,             // scalar for undefined entries
                      NDims,                 // num of dimensions
                      DimNames,              // dimension names
                      false                  // not time dependent
@@ -893,7 +874,6 @@ void HorzMesh::defineMeshFields() {
                      "",                 // CF standard name
                      0.0,                // min valid value
                      9.99E+30,           // max valid value
-                     -9.99E+30,          // scalar for undefined entries
                      NDims,              // num of dimensions
                      DimNames,           // dimension names
                      false               // not time dependent
@@ -911,7 +891,6 @@ void HorzMesh::defineMeshFields() {
                      "",           // CF standard name
                      -3.1415927,   // min valid value
                      3.1415927,    // max valid value
-                     -9.99E+30,    // scalar for undefined entries
                      NDims,        // num of dimensions
                      DimNames,     // dimension names
                      false         // not time dependent
@@ -934,7 +913,6 @@ void HorzMesh::defineMeshFields() {
                      "",                         // CF standard name
                      0.0,                        // min valid value
                      9.99E+30,                   // max valid value
-                     -9.99E+30,                  // scalar for undefined entries
                      NDims,                      // num of dimensions
                      DimNames,                   // dimension names
                      false                       // not time dependent
@@ -955,7 +933,6 @@ void HorzMesh::defineMeshFields() {
                      "",             // CF standard name
                      -1.0,           // min valid value
                      1.0,            // max valid value
-                     -9.99E+30,      // scalar for undefined entries
                      NDims,          // num of dimensions
                      DimNames,       // dimension names
                      false           // not time dependent
@@ -976,7 +953,6 @@ void HorzMesh::defineMeshFields() {
                      "coriolis_parameter", // CF standard name
                      -1.0E-3,              // min valid value
                      1.0E-3,               // max valid value
-                     -9.99E+30,            // scalar for undefined entries
                      NDims,                // num of dimensions
                      DimNames,             // dimension names
                      false                 // not time dependent
@@ -994,7 +970,6 @@ void HorzMesh::defineMeshFields() {
                      "coriolis_parameter", // CF standard name
                      -1.0E-3,              // min valid value
                      1.0E-3,               // max valid value
-                     -9.99E+30,            // scalar for undefined entries
                      NDims,                // num of dimensions
                      DimNames,             // dimension names
                      false                 // not time dependent
@@ -1012,7 +987,6 @@ void HorzMesh::defineMeshFields() {
                      "coriolis_parameter", // CF standard name
                      -1.0E-3,              // min valid value
                      1.0E-3,               // max valid value
-                     -9.99E+30,            // scalar for undefined entries
                      NDims,                // num of dimensions
                      DimNames,             // dimension names
                      false                 // not time dependent
diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp
index ceabc0ede033..880f9fd11307 100644
--- a/components/omega/src/ocn/OceanInit.cpp
+++ b/components/omega/src/ocn/OceanInit.cpp
@@ -149,7 +149,14 @@ int ocnInit(MPI_Comm Comm ///< [in] ocean MPI communicator
    OceanState *DefState = OceanState::getDefault();
    I4 CurTimeLevel      = 0;
    DefState->exchangeHalo(CurTimeLevel);
+
+   // Enforce layer masks on state and tracer variables: fully-inactive layers
+   // get FillValueReal, boundary layers get 0, active layers keep their
+   // IC/restart value.
+   DefState->applyLayerMasks(CurTimeLevel);
+
    DefState->copyToHost(CurTimeLevel);
+   VertCoord::getDefault()->initSurfacePressure(Halo::getDefault());
 
    Forcing *DefForcing = Forcing::getDefault();
    DefForcing->exchangeHalo();
diff --git a/components/omega/src/ocn/OceanState.cpp b/components/omega/src/ocn/OceanState.cpp
index e696d06a5a7e..f497d6c1a693 100644
--- a/components/omega/src/ocn/OceanState.cpp
+++ b/components/omega/src/ocn/OceanState.cpp
@@ -17,6 +17,8 @@
 #include "MachEnv.h"
 #include "OmegaKokkos.h"
 #include "TimeStepper.h"
+#include "Tracers.h"
+#include "VertCoord.h"
 
 namespace OMEGA {
 
@@ -115,6 +117,14 @@ OceanState::OceanState(
                                        NEdgesSize, NVertLayers);
    }
 
+   // Apply edge layer mask to NormalVelocity for all time levels so that
+   // fully-inactive layers are initialized with FillValueReal and boundary
+   // layers with 0.
+   VertCoord *DefVCoord = VertCoord::getDefault();
+   for (int I = 0; I < NTimeLevels; I++) {
+      DefVCoord->applyEdgeLayerMask(NormalVelocity[I], NEdgesAll);
+   }
+
    // Register fields and metadata for IO
    defineFields();
 
@@ -202,9 +212,8 @@ void OceanState::defineFields() {
                      "sea_water_velocity",                // CF standard Name
                      -9.99E+10,                           // min valid value
                      9.99E+10,                            // max valid value
-                     -9.99E+30, // scalar for undefined entries
-                     NDims,     // number of dimensions
-                     DimNames   // dimension names
+                     NDims,   // number of dimensions
+                     DimNames // dimension names
        );
 
    DimNames[0] = "NCells";
@@ -215,17 +224,22 @@ void OceanState::defineFields() {
                      "cell_thickness", // CF standard Name
                      0.0,              // min valid value
                      9.99E+30,         // max valid value
-                     -9.99E+30,        // scalar used for undefined entries
                      NDims,            // number of dimensions
                      DimNames          // dimension names
        );
 
-   // Create a field group for state fields
+   // Create a field group for state fields. VertCoord initializes before
+   // OceanState and owns SurfacePressure, which it also adds to this "State"
+   // group, so the group may already exist; reuse it if so.
    StateGroupName = "State";
    if (Name != "Default") {
       StateGroupName.append(Name);
    }
-   auto StateGroup = FieldGroup::create(StateGroupName);
+   std::shared_ptr<FieldGroup> StateGroup;
+   if (FieldGroup::exists(StateGroupName))
+      StateGroup = FieldGroup::get(StateGroupName);
+   else
+      StateGroup = FieldGroup::create(StateGroupName);
 
    // Add restart group if needed
    if (!FieldGroup::exists("Restart"))
@@ -295,6 +309,21 @@ void OceanState::copyToHost(const I4 TimeLevel) {
    deepCopy(NormalVelocityH[TimeIndex], NormalVelocity[TimeIndex]);
 } // end copyToHost
 
+//------------------------------------------------------------------------------
+// Apply layer masks to NormalVelocity, PseudoThickness, Temperature, and
+// Salinity at the given time level after reading an IC or restart file.
+void OceanState::applyLayerMasks(const I4 TimeLevel) {
+
+   VertCoord *DefVCoord = VertCoord::getDefault();
+   DefVCoord->applyEdgeLayerMask(getNormalVelocity(TimeLevel), NEdgesAll);
+   DefVCoord->applyCellLayerMask(getPseudoThickness(TimeLevel), NCellsAll);
+   DefVCoord->applyCellLayerMask(Tracers::getByName(TimeLevel, "Temperature"),
+                                 NCellsAll);
+   DefVCoord->applyCellLayerMask(Tracers::getByName(TimeLevel, "Salinity"),
+                                 NCellsAll);
+
+} // end applyLayerMasks
+
 //------------------------------------------------------------------------------
 // Perform state halo exchange
 // TimeLevel == [1:new, 0:current, -1:previous, -2:two times ago, ...]
@@ -320,11 +349,12 @@ void OceanState::updateTimeLevels() {
    // Update current time index for pseudo-thickness and normal velocity
    CurTimeIndex = (CurTimeIndex + 1) % NTimeLevels;
 
-   // Update IOField data associations
+   // Update IOField data associations. Pass false to avoid overwriting the
+   // just-computed state with fill values; this is a pointer update only.
    Field::attachFieldData<Array2DReal>(NormalVelocityFldName,
-                                       NormalVelocity[CurTimeIndex]);
+                                       NormalVelocity[CurTimeIndex], false);
    Field::attachFieldData<Array2DReal>(PseudoThicknessFldName,
-                                       PseudoThickness[CurTimeIndex]);
+                                       PseudoThickness[CurTimeIndex], false);
 
 } // end updateTimeLevels
 
diff --git a/components/omega/src/ocn/OceanState.h b/components/omega/src/ocn/OceanState.h
index 3632449ccc32..fe3891e05b7c 100644
--- a/components/omega/src/ocn/OceanState.h
+++ b/components/omega/src/ocn/OceanState.h
@@ -122,6 +122,10 @@ class OceanState {
    /// Get normal velocity host array at given time level
    HostArray2DReal getNormalVelocityH(const I4 TimeLevel) const;
 
+   /// Apply layer masks to NormalVelocity, PseudoThickness, Temperature, and
+   /// Salinity after reading an IC or restart file
+   void applyLayerMasks(const I4 TimeLevel);
+
    /// Exchange halo
    void exchangeHalo(const I4 TimeLevel);
 
diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp
index ba0df5dd692a..eca1252826b7 100644
--- a/components/omega/src/ocn/Tendencies.cpp
+++ b/components/omega/src/ocn/Tendencies.cpp
@@ -338,7 +338,7 @@ void Tendencies::defineFields() {
    auto PseudoThicknessTendField =
        Field::create(PseudoThicknessTendFieldName, "Pseudo-thickness tendency",
                      "m/s", "cell_thickness_tendency", -9.99E+10, 9.99E+10,
-                     -9.99E+30, NDims, DimNamesThickness);
+                     NDims, DimNamesThickness);
    NDims = 3;
    std::vector<std::string> DimNamesTracer(NDims);
    DimNamesTracer[0]    = "NTracers";
@@ -346,7 +346,7 @@ void Tendencies::defineFields() {
    DimNamesTracer[2]    = "NVertLayers";
    auto TracerTendField = Field::create(
        TracerTendFieldName, "Tracer tendency", "kg/m^3/s", "tracer_tendency",
-       -9.99E+10, 9.99E+10, -9.99E+30, NDims, DimNamesTracer);
+       -9.99E+10, 9.99E+10, NDims, DimNamesTracer);
    NDims = 2;
    std::vector<std::string> DimNamesVelocity(NDims);
    DimNamesVelocity[0] = "NEdges";
@@ -354,7 +354,7 @@ void Tendencies::defineFields() {
    auto NormalVelocityTendField =
        Field::create(NormalVelocityTendFieldName, "Normal velocity tendency",
                      "m/s^2", "sea_water_velocity_tendency", -9.99E+10,
-                     9.99E+10, -9.99E+30, NDims, DimNamesVelocity);
+                     9.99E+10, NDims, DimNamesVelocity);
 
    std::string TendGroupName = "Tendencies";
    if (Name != "Default") {
@@ -517,15 +517,12 @@ void Tendencies::computeVelocityTendenciesOnly(
 
    Pacer::start("Tend:computeVelocityTendenciesOnly", 1);
 
-   parallelForOuter(
-       {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) {
-          const int KMin = MinLayerEdgeBot(IEdge);
-          const int KMax = MaxLayerEdgeTop(IEdge);
-
-          parallelForInner(
-              Team, Range{KMin, KMax},
-              INNER_LAMBDA(int K) { LocNormalVelocityTend(IEdge, K) = 0; });
-       });
+   // Zero NVTend over all layers where at least one neighboring cell is active:
+   // [MinLayerEdgeTop, MaxLayerEdgeBot]. This includes boundary edges adjacent
+   // to land or bathymetry steps, which have no computed tendency but should
+   // show 0, not FillValueReal. Fully inactive layers (outside this range)
+   // keep their FillValueReal from attachData().
+   VCoord->zeroEdgeField(NormalVelocityTend, Mesh->NEdgesAll);
 
    // Compute potential vorticity horizontal advection
    const Array2DReal &FluxPseudoThickEdge =
diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp
index 4d7ce499941f..5b142888142c 100644
--- a/components/omega/src/ocn/TendencyTerms.cpp
+++ b/components/omega/src/ocn/TendencyTerms.cpp
@@ -88,9 +88,7 @@ TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh,
                        Mesh->NEdgesAll, VCoord->NVertLayers),
       NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
       CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell),
-      DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell) {
-   deepCopy(HighOrderFlxHorz, 0);
-}
+      DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell) {}
 
 TracerDiffOnCell::TracerDiffOnCell(const HorzMesh *Mesh,
                                    const VertCoord *VCoord)
diff --git a/components/omega/src/ocn/TracerDefs.inc b/components/omega/src/ocn/TracerDefs.inc
index fb2326e9c203..adafaa945f11 100644
--- a/components/omega/src/ocn/TracerDefs.inc
+++ b/components/omega/src/ocn/TracerDefs.inc
@@ -27,16 +27,12 @@ static void defineAllTracers() {
           "sea_water_conservative_temperature", ///< [in] CF standard Name
           -273.15,                              ///< [in] min valid field value
           100.0,                                ///< [in] max valid field value
-          1.e33,                                ///< [in] value of undef entries
           IndxTemp);                            ///< [out] (optional) static index
 
    define("Salinity", "Absolute Salinity", "g kg-1",
-          "sea_water_absolute_salinity", 0.0, 50.0, 1.e33, IndxSalt);
+          "sea_water_absolute_salinity", 0.0, 50.0, IndxSalt);
 
-   define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, 1.e33,
-          IndxDebug1);
-   define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, 1.e33,
-          IndxDebug2);
-   define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, 1.e33,
-          IndxDebug3);
+   define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, IndxDebug1);
+   define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, IndxDebug2);
+   define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, IndxDebug3);
 }
diff --git a/components/omega/src/ocn/Tracers.cpp b/components/omega/src/ocn/Tracers.cpp
index efa7b00b635f..c84fa5b33454 100644
--- a/components/omega/src/ocn/Tracers.cpp
+++ b/components/omega/src/ocn/Tracers.cpp
@@ -177,8 +177,7 @@ void Tracers::init() {
 //---------------------------------------------------------------------------
 I4 Tracers::define(const std::string &Name, const std::string &Description,
                    const std::string &Units, const std::string &StdName,
-                   const Real ValidMin, const Real ValidMax,
-                   const Real FillValue, I4 &Index) {
+                   const Real ValidMin, const Real ValidMax, I4 &Index) {
 
    // Do nothing if this tracer is not selected
    if (TracerIndexes.find(Name) == TracerIndexes.end()) {
@@ -201,9 +200,9 @@ I4 Tracers::define(const std::string &Name, const std::string &Description,
 
    // create a tracer field
    std::string TracerFieldName = Name;
-   auto TracerField = Field::create(TracerFieldName, Description, Units,
-                                    StdName, ValidMin, ValidMax, FillValue,
-                                    TracerDimNames.size(), TracerDimNames);
+   auto TracerField =
+       Field::create(TracerFieldName, Description, Units, StdName, ValidMin,
+                     ValidMax, TracerDimNames.size(), TracerDimNames);
    if (!TracerField) {
       LOG_ERROR("Tracers: Tracer field '{}' is not created", TracerFieldName);
       return -2;
@@ -447,7 +446,8 @@ void Tracers::updateTimeLevels() {
 
       Array2DReal TracerSubview = Kokkos::subview(
           TracerArrays[CurTimeIndex], TracerIndex, Kokkos::ALL, Kokkos::ALL);
-      TracerField->attachData<Array2DReal>(TracerSubview);
+      // Pass false to avoid overwriting just-computed tracer values with fill.
+      TracerField->attachData<Array2DReal>(TracerSubview, false);
    }
 
    return;
diff --git a/components/omega/src/ocn/Tracers.h b/components/omega/src/ocn/Tracers.h
index 6480880a6235..a1ee899056b8 100644
--- a/components/omega/src/ocn/Tracers.h
+++ b/components/omega/src/ocn/Tracers.h
@@ -81,7 +81,6 @@ class Tracers {
           const std::string &StdName,     ///< [in] CF standard Name
           const Real ValidMin,            ///< [in] min valid field value
           const Real ValidMax,            ///< [in] max valid field value
-          const Real FillValue,           ///< [in] value for undef entries
           I4 &Index = IndxInvalid         ///< [out] (optional) index value
    );
 
diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp
index 04e363a0ad3a..0a29ec7f6b17 100644
--- a/components/omega/src/ocn/VertAdv.cpp
+++ b/components/omega/src/ocn/VertAdv.cpp
@@ -128,8 +128,7 @@ void VertAdv::defineFields() {
       LowOrderVertFluxFldName.append(Name);
    }
 
-   const Real FillValueReal = -9.99e30;
-   I4 NDims                 = 2;
+   I4 NDims = 2;
    std::vector<std::string> DimNames(NDims);
    DimNames[0] = "NCells";
    DimNames[1] = "NVertLayersP1";
@@ -143,7 +142,6 @@ void VertAdv::defineFields() {
        "",                               // CF standard Name
        std::numeric_limits<Real>::min(), // min valid value
        std::numeric_limits<Real>::max(), // max valid value
-       FillValueReal,                    // scalar for undefined entries
        NDims,                            // number of dimensions
        DimNames                          // dimension names
    );
@@ -156,9 +154,8 @@ void VertAdv::defineFields() {
                      "",                      // CF standard Name
                      std::numeric_limits<Real>::min(), // min valid value
                      std::numeric_limits<Real>::max(), // max valid value
-                     FillValueReal, // scalar for undefined entries
-                     NDims,         // number of dimensions
-                     DimNames       // dimension names
+                     NDims,                            // number of dimensions
+                     DimNames                          // dimension names
        );
 
    NDims = 3;
@@ -174,7 +171,6 @@ void VertAdv::defineFields() {
        "",                               // CF standard Name
        std::numeric_limits<Real>::min(), // min valid value
        std::numeric_limits<Real>::max(), // max valid value
-       FillValueReal,                    // scalar for undefined entries
        NDims,                            // number of dimensions
        DimNames                          // dimension names
    );
@@ -187,7 +183,6 @@ void VertAdv::defineFields() {
        "",                               // CF standard Name
        std::numeric_limits<Real>::min(), // min valid value
        std::numeric_limits<Real>::max(), // max valid value
-       FillValueReal,                    // scalar for undefined entries
        NDims,                            // number of dimensions
        DimNames                          // dimension names
    );
@@ -384,6 +379,7 @@ void VertAdv::computeVerticalPseudoVelocity(
    OMEGA_SCOPE(LocEOnC, Mesh->EdgesOnCell);
    OMEGA_SCOPE(LocDvE, Mesh->DvEdge);
    OMEGA_SCOPE(LocESOnC, Mesh->EdgeSignOnCell);
+   OMEGA_SCOPE(LocMaxLayerEdgeTop, VCoord->MaxLayerEdgeTop);
 
    // Loop over all cells owned by the task
    parallelForOuter(
@@ -407,12 +403,16 @@ void VertAdv::computeVerticalPseudoVelocity(
                  const I4 KLen            = chunkLength(KChunk, KStart, KMax);
 
                  for (int J = 0; J < LocNEOnC(ICell); ++J) {
-                    const I4 JEdge = LocEOnC(ICell, J);
+                    const I4 JEdge    = LocEOnC(ICell, J);
+                    const I4 MaxKEdge = LocMaxLayerEdgeTop(JEdge);
                     for (int KVec = 0; KVec < KLen; ++KVec) {
                        const I4 K = KStart + KVec;
-                       DivHUTmp[KVec] -= LocDvE(JEdge) * LocESOnC(ICell, J) *
-                                         FluxPseudoThickEdge(JEdge, K) *
-                                         NormalVelocity(JEdge, K) * InvAreaCell;
+                       if (K <= MaxKEdge) {
+                          DivHUTmp[KVec] -= LocDvE(JEdge) * LocESOnC(ICell, J) *
+                                            FluxPseudoThickEdge(JEdge, K) *
+                                            NormalVelocity(JEdge, K) *
+                                            InvAreaCell;
+                       }
                     }
                  }
                  for (int KVec = 0; KVec < KLen; ++KVec) {
diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp
index 146ceb414efc..dfedfa1c32cc 100644
--- a/components/omega/src/ocn/VertCoord.cpp
+++ b/components/omega/src/ocn/VertCoord.cpp
@@ -143,6 +143,7 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord
    MaxLayerCell    = Array1DI4("MaxLayerCell", NCellsSize);
    MinLayerCell    = Array1DI4("MinLayerCell", NCellsSize);
    BottomGeomDepth = Array1DReal("BottomGeomDepth", NCellsSize);
+   SurfacePressure = Array1DReal("SurfacePressure", NCellsSize);
    PressureInterface =
        Array2DReal("PressureInterface", NCellsSize, NVertLayersP1);
    PressureMid     = Array2DReal("PressureMid", NCellsSize, NVertLayers);
@@ -157,10 +158,6 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord
    VertCoordMovementWeights =
        Array1DReal("VertCoordMovementWeights", NVertLayers);
 
-   // TODO: Temporary handling of SurfacePressure
-   SurfacePressure = Array1DReal("SurfacePressure", NCellsSize);
-   deepCopy(SurfacePressure, 0);
-
    // Make host copies for device arrays not being read from file
    PressureInterfaceH     = createHostMirrorCopy(PressureInterface);
    PressureMidH           = createHostMirrorCopy(PressureMid);
@@ -170,6 +167,10 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord
    GeopotentialMidH       = createHostMirrorCopy(GeopotentialMid);
    PseudoThicknessTargetH = createHostMirrorCopy(PseudoThicknessTarget);
 
+   // SurfacePressure is read from the initial-state/restart stream later, in
+   // OceanInit; its host mirror is refreshed then by initSurfacePressure().
+   SurfacePressureH = createHostMirrorCopy(SurfacePressure);
+
    // Define field metadata
    defineFields();
 
@@ -230,6 +231,7 @@ void VertCoord::defineFields() {
    SshFldName                   = "SshCell";
    GeopotFldName                = "GeopotentialMid";
    PseudoThicknessTargetFldName = "PseudoThicknessTarget";
+   SurfacePressureFldName       = "SurfacePressure";
 
    if (Name != "Default") {
       MinLayerCellFldName.append(Name);
@@ -244,12 +246,11 @@ void VertCoord::defineFields() {
       GeopotFldName.append(Name);
       PseudoThicknessTargetFldName.append(Name);
       SshFldName.append(Name);
+      SurfacePressureFldName.append(Name);
    }
 
    // Create fields for VertCoord variables
-   const I4 FillValueI4     = -999;
-   const Real FillValueReal = -9.99e30;
-   int NDims                = 1;
+   int NDims = 1;
    std::vector<std::string> DimNames(NDims);
    DimNames[0] = "NCells";
 
@@ -260,9 +261,8 @@ void VertCoord::defineFields() {
        "",                                          // CF standard Name
        0,                                           // min valid value
        std::numeric_limits<I4>::max(),              // max valid value
-       FillValueI4, // scalar for undefined entries
-       NDims,       // number of dimensions
-       DimNames     // dimension names
+       NDims,                                       // number of dimensions
+       DimNames                                     // dimension names
    );
 
    auto MaxLayerCellField = Field::create(
@@ -272,9 +272,8 @@ void VertCoord::defineFields() {
        "",                                         // CF standard Name
        -1,                                         // min valid value
        std::numeric_limits<I4>::max(),             // max valid value
-       FillValueI4, // scalar for undefined entries
-       NDims,       // number of dimensions
-       DimNames     // dimension names
+       NDims,                                      // number of dimensions
+       DimNames                                    // dimension names
    );
 
    auto BottomGeomDepthField = Field::create(
@@ -285,7 +284,6 @@ void VertCoord::defineFields() {
        "sea_floor_depth_below_geoid",    // CF standard Name
        0.0,                              // min valid value
        std::numeric_limits<Real>::max(), // max valid value
-       FillValueReal,                    // scalar for undefined entries
        NDims,                            // number of dimensions
        DimNames                          // dimension names
    );
@@ -297,11 +295,25 @@ void VertCoord::defineFields() {
        "sea_surface_height",                // CF standard Name
        std::numeric_limits<Real>::min(),    // min valid value
        std::numeric_limits<Real>::max(),    // max valid value
-       FillValueReal,                       // scalar for undefined entries
        NDims,                               // number of dimensions
        DimNames                             // dimension names
    );
 
+   auto SurfacePressureField = Field::create(
+       SurfacePressureFldName,                             // field name
+       "Relative pressure at the top of the ocean column", // long name
+       "Pa",                                               // units
+       "",                                                 // CF standard Name
+       -9.99E+10,                                          // min valid value
+       9.99E+10,                                           // max valid value
+       NDims,                                              // number of dims
+       DimNames                                            // dimension names
+   );
+
+   // SurfacePressure is optional in the initial-state/restart file; when it is
+   // absent the read is skipped and it defaults to zero in initSurfacePressure.
+   SurfacePressureField->setOptionalRead(true);
+
    NDims = 2;
    DimNames.resize(NDims);
    DimNames[1] = "NVertLayers";
@@ -314,7 +326,6 @@ void VertCoord::defineFields() {
        "",                               // CF standard Name
        0.0,                              // min valid value
        std::numeric_limits<Real>::max(), // max valid value
-       FillValueReal,                    // scalar for undefined entries
        NDims,                            // number of dimensions
        DimNames                          // dimension names
    );
@@ -331,7 +342,6 @@ void VertCoord::defineFields() {
        "",                               // CF standard Name
        0.0,                              // min valid value
        std::numeric_limits<Real>::max(), // max valid value
-       FillValueReal,                    // scalar for undefined entries
        NDims,                            // number of dimensions
        DimNames                          // dimension names
    );
@@ -342,15 +352,14 @@ void VertCoord::defineFields() {
    DimNames[1] = "NVertLayersP1";
 
    auto PressureInterfaceField = Field::create(
-       PressInterfFldName,                      // field name
-       "Pressure at vertical layer interfaces", // long name or description
-       "Pa",                                    // units
-       "sea_water_pressure",                    // CF standard Name
-       0.0,                                     // min valid value
-       std::numeric_limits<Real>::max(),        // max valid value
-       FillValueReal,                           // scalar for undefined entries
-       NDims,                                   // number of dimensions
-       DimNames                                 // dimension names
+       PressInterfFldName,                                       // field name
+       "Relative pressure (gauge pressure) at layer interfaces", // long name
+       "Pa",                                                     // units
+       "",                               // CF standard Name
+       -AtmRefP,                         // min valid value
+       std::numeric_limits<Real>::max(), // max valid value
+       NDims,                            // number of dimensions
+       DimNames                          // dimension names
    );
 
    auto GeomZInterfaceField = Field::create(
@@ -360,7 +369,6 @@ void VertCoord::defineFields() {
        "height",                               // CF standard Name
        std::numeric_limits<Real>::min(),       // min valid value
        std::numeric_limits<Real>::max(),       // max valid value
-       FillValueReal,                          // scalar for undefined entries
        NDims,                                  // number of dimensions
        DimNames                                // dimension names
    );
@@ -368,15 +376,14 @@ void VertCoord::defineFields() {
    DimNames[1] = "NVertLayers";
 
    auto PressureMidField = Field::create(
-       PressMidFldName,                        // field name
-       "Pressure at vertical layer midpoints", // long name or description
-       "Pa",                                   // units
-       "sea_water_pressure",                   // CF standard Name
-       0.0,                                    // min valid value
-       std::numeric_limits<Real>::max(),       // max valid value
-       FillValueReal,                          // scalar for undefined entries
-       NDims,                                  // number of dimensions
-       DimNames                                // dimension names
+       PressMidFldName,                                         // field name
+       "Relative pressure (gauge pressure) at layer midpoints", // long name
+       "Pa",                                                    // units
+       "",                               // CF standard Name
+       -AtmRefP,                         // min valid value
+       std::numeric_limits<Real>::max(), // max valid value
+       NDims,                            // number of dimensions
+       DimNames                          // dimension names
    );
 
    auto GeomZMidField = Field::create(
@@ -386,7 +393,6 @@ void VertCoord::defineFields() {
        "height",                              // CF standard Name
        std::numeric_limits<Real>::min(),      // min valid value
        std::numeric_limits<Real>::max(),      // max valid value
-       FillValueReal,                         // scalar for undefined entries
        NDims,                                 // number of dimensions
        DimNames                               // dimension names
    );
@@ -398,7 +404,6 @@ void VertCoord::defineFields() {
        "geopotential",                    // CF standard Name
        std::numeric_limits<Real>::min(),  // min valid value
        std::numeric_limits<Real>::max(),  // max valid value
-       FillValueReal,                     // scalar for undefined entries
        NDims,                             // number of dimensions
        DimNames                           // dimension names
    );
@@ -411,7 +416,6 @@ void VertCoord::defineFields() {
        "",                                    // CF standard Name
        0.0,                                   // min valid value
        std::numeric_limits<Real>::max(),      // max valid value
-       FillValueReal,                         // scalar for undefined entries
        NDims,                                 // number of dimensions
        DimNames                               // dimension names
    );
@@ -459,6 +463,26 @@ void VertCoord::defineFields() {
    PseudoThicknessTargetField->attachData<Array2DReal>(PseudoThicknessTarget);
    SshField->attachData<Array1DReal>(SshCell);
 
+   // Add SurfacePressure to the State and Restart groups so it is written to
+   // the history/restart files and read from the initial-state and restart
+   // files when present. The read is optional (see setOptionalRead above): if
+   // the variable is absent, SurfacePressure defaults to zero. VertCoord
+   // initializes before OceanState, which also contributes fields to these
+   // groups, so create each group only if it does not already exist.
+   std::string StateGrpName = "State";
+   if (Name != "Default") {
+      StateGrpName.append(Name);
+   }
+   if (!FieldGroup::exists(StateGrpName))
+      FieldGroup::create(StateGrpName);
+   FieldGroup::addFieldToGroup(SurfacePressureFldName, StateGrpName);
+
+   if (!FieldGroup::exists("Restart"))
+      FieldGroup::create("Restart");
+   FieldGroup::addFieldToGroup(SurfacePressureFldName, "Restart");
+
+   SurfacePressureField->attachData<Array1DReal>(SurfacePressure);
+
 } // end defineFields
 
 //------------------------------------------------------------------------------
@@ -485,6 +509,12 @@ VertCoord::~VertCoord() {
       FieldGroup::destroy(GroupName);
    }
 
+   // SurfacePressure belongs to the shared "State"/"Restart" groups (owned by
+   // OceanState), so destroy only the field here, not those groups.
+   if (Field::exists(SurfacePressureFldName)) {
+      Field::destroy(SurfacePressureFldName);
+   }
+
 } // end destructor
 
 //------------------------------------------------------------------------------
@@ -524,15 +554,6 @@ void VertCoord::setStreamArrays(const bool ReadStream, Halo *MeshHalo) {
    // uninitialized and will need to be initialized explicitly if needed.
    if (ReadStream) {
 
-      I4 FillValueI4     = -1;
-      Real FillValueReal = -999._Real;
-
-      deepCopy(MinLayerCell, FillValueI4);
-      deepCopy(MaxLayerCell, FillValueI4);
-      deepCopy(BottomGeomDepth, FillValueReal);
-      deepCopy(RefPseudoThickness, FillValueReal);
-      deepCopy(VertCoordMovementWeights, FillValueReal);
-
       // Fetch input stream and validate
       std::string StreamName = "InitialVertCoord";
       if (Name != "Default") {
@@ -558,76 +579,93 @@ void VertCoord::setStreamArrays(const bool ReadStream, Halo *MeshHalo) {
          Err = IOStream::read(StreamName);
          if (Err.isFail()) {
             LOG_INFO("VertCoord: Error while reading {} stream", StreamName);
-            I4 Sum1 = 0;
+            I4 NFill1 = 0;
             parallelReduce(
                 {MinLayerCell.extent_int(0)},
-                KOKKOS_LAMBDA(int I, int &Accum) {
-                   Accum += LocMinLayerCell(I);
+                KOKKOS_LAMBDA(int I, I4 &Count) {
+                   if (LocMinLayerCell(I) == FillValueI4)
+                      ++Count;
                 },
-                Sum1);
-            if (Sum1 < 0) {
+                NFill1);
+            if (NFill1 > 0) {
                LOG_INFO("VertCoord: Error reading MinLayerCell from {}, "
                         "using MinLayerCell = 0",
                         StreamName);
                deepCopy(MinLayerCell, 1);
             }
-            I4 Sum2 = 0;
+            I4 NFill2 = 0;
             parallelReduce(
                 {MaxLayerCell.extent_int(0)},
-                KOKKOS_LAMBDA(int I, int &Accum) {
-                   Accum += LocMaxLayerCell(I);
+                KOKKOS_LAMBDA(int I, I4 &Count) {
+                   if (LocMaxLayerCell(I) == FillValueI4)
+                      ++Count;
                 },
-                Sum2);
-            if (Sum2 < 0) {
+                NFill2);
+            if (NFill2 > 0) {
                LOG_INFO("VertCoord: Error reading MaxLayerCell from {}, "
                         "using MaxLayerCell = NVertLayers - 1",
                         StreamName);
                deepCopy(MaxLayerCell, NVertLayers);
             }
-            Real Sum3 = 0.;
+            I4 NFill3 = 0;
             parallelReduce(
                 {BottomGeomDepth.extent_int(0)},
-                KOKKOS_LAMBDA(int I, Real &Accum) {
-                   Accum += LocBottomGeomDepth(I);
+                KOKKOS_LAMBDA(int I, I4 &Count) {
+                   if (LocBottomGeomDepth(I) == FillValueReal)
+                      ++Count;
                 },
-                Sum3);
-            if (Sum3 < 0.) {
+                NFill3);
+            if (NFill3 > 0) {
                ABORT_ERROR("VertCoord: Error reading BottomGeomDepth from {}",
                            StreamName);
             }
-            Real Sum4 = 0.;
+            I4 NFill4 = 0;
             parallelReduce(
                 {RefPseudoThickness.extent_int(0),
                  RefPseudoThickness.extent_int(1)},
-                KOKKOS_LAMBDA(int I, int J, Real &Accum) {
-                   Accum += LocPseudoThick(I, J);
+                KOKKOS_LAMBDA(int I, int J, I4 &Count) {
+                   if (LocPseudoThick(I, J) == FillValueReal)
+                      ++Count;
                 },
-                Sum4);
-            if (Sum4 < 0.) {
+                NFill4);
+            if (NFill4 > 0) {
                ABORT_ERROR("VertCoord: Error reading RefPseudoThickness "
                            "from {}",
                            StreamName);
             }
-            Real Sum5     = 0.;
-            I4 NumNonZero = 0;
+            I4 NFill5 = 0;
             parallelReduce(
                 {VertCoordMovementWeights.extent_int(0)},
-                KOKKOS_LAMBDA(int I, Real &Accum, I4 &NonZeroCount) {
-                   Real W = LocVCoordMvmtWgts(I);
-                   Accum += W;
-                   if (W != 0) {
-                      NonZeroCount += 1;
-                   }
+                KOKKOS_LAMBDA(int I, I4 &Count) {
+                   if (LocVCoordMvmtWgts(I) == FillValueReal)
+                      ++Count;
                 },
-                Sum5, NumNonZero);
-            if (Sum5 < 0.) {
-               ABORT_ERROR("VertCoord: Error reading VertCoordMovementWeights "
-                           "from {}",
-                           StreamName);
-            }
-            if (NumNonZero == 0) {
-               // TODO: ABORT_ERROR when all weights equal 0
+                NFill5);
+            if (NFill5 > 0) {
+               // Stream did not populate weights — use default
                deepCopy(VertCoordMovementWeights, 1._Real);
+            } else {
+               Real Sum5     = 0.;
+               I4 NumNonZero = 0;
+               parallelReduce(
+                   {VertCoordMovementWeights.extent_int(0)},
+                   KOKKOS_LAMBDA(int I, Real &Accum, I4 &NonZeroCount) {
+                      Real W = LocVCoordMvmtWgts(I);
+                      Accum += W;
+                      if (W != 0) {
+                         NonZeroCount += 1;
+                      }
+                   },
+                   Sum5, NumNonZero);
+               if (Sum5 < 0.) {
+                  ABORT_ERROR(
+                      "VertCoord: Error reading VertCoordMovementWeights "
+                      "from {}",
+                      StreamName);
+               } else if (NumNonZero == 0) {
+                  // TODO: ABORT_ERROR when all weights equal 0
+                  deepCopy(VertCoordMovementWeights, 1._Real);
+               }
             }
          }
       } else {
@@ -669,6 +707,35 @@ void VertCoord::setStreamArrays(const bool ReadStream, Halo *MeshHalo) {
    VertCoordMovementWeightsH = createHostMirrorCopy(VertCoordMovementWeights);
 }
 
+//------------------------------------------------------------------------------
+// Apply the zero default when SurfacePressure was absent from the input, then
+// exchange halo and copy SurfacePressure to host after the initial-state or
+// restart stream has been read.
+void VertCoord::initSurfacePressure(Halo *MeshHalo) {
+   // SurfacePressure is an optional read. If it was not present in the
+   // initial-state or restart file, its owned cells still hold the fill value
+   // set by attachData; in that case default the whole array to zero. Only the
+   // owned cells are checked because halo cells are not populated until the
+   // exchange below.
+   OMEGA_SCOPE(LocSurfacePressure, SurfacePressure);
+   I4 NFill = 0;
+   parallelReduce(
+       {NCellsOwned},
+       KOKKOS_LAMBDA(int ICell, I4 &Count) {
+          if (LocSurfacePressure(ICell) == FillValueReal)
+             ++Count;
+       },
+       NFill);
+   if (NFill > 0) {
+      LOG_INFO("VertCoord: SurfacePressure not found in input, "
+               "using SurfacePressure = 0");
+      deepCopy(SurfacePressure, 0._Real);
+   }
+
+   MeshHalo->exchangeFullArrayHalo(SurfacePressure, OnCell);
+   deepCopy(SurfacePressureH, SurfacePressure);
+}
+
 //------------------------------------------------------------------------------
 // Compute min and max layer indices for edges based on MinLayerCell and
 // MaxLayerCell
@@ -964,11 +1031,12 @@ void VertCoord::setMasks() {
 } // end setMasks()
 
 //------------------------------------------------------------------------------
-// Compute the pressure at each layer interface and midpoint given the
-// PseudoThickness and SurfacePressure. Hierarchical parallelism is used with a
-// parallel_for loop over all cells and a parallel_scan performing a prefix sum
-// in each column to compute pressure from the top-most active layer to the
-// bottom-most active layer.
+// Compute the relative pressure (gauge pressure, i.e., absolute pressure minus
+// the standard atmosphere of 101325 Pa) at each layer interface and midpoint
+// given PseudoThickness and SurfacePressure (also relative). Hierarchical
+// parallelism is used with a parallel_for loop over all cells and a
+// parallel_scan performing a prefix sum in each column from the top-most active
+// layer to the bottom-most active layer.
 void VertCoord::computePressure(
     const Array2DReal &PseudoThickness, // [in] pseudo-thickness
     const Array1DReal &SurfacePressure  // [in] surface pressure
@@ -1186,6 +1254,87 @@ VertCoord *VertCoord::get(const std::string Name ///< [in] Name of VertCoord
 
 } // end get VertCoord
 
+//------------------------------------------------------------------------------
+// Zero all layers in [MinLayerEdgeTop, MaxLayerEdgeBot] of an edge field.
+void VertCoord::zeroEdgeField(Array2DReal Arr, I4 NEdgesAll) const {
+   OMEGA_SCOPE(LocMinLayerEdgeTop, MinLayerEdgeTop);
+   OMEGA_SCOPE(LocMaxLayerEdgeBot, MaxLayerEdgeBot);
+   parallelForOuter(
+       {NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) {
+          const int KTop = LocMinLayerEdgeTop(IEdge);
+          const int KBot = Kokkos::max(0, LocMaxLayerEdgeBot(IEdge));
+          parallelForInner(
+              Team, Range{KTop, KBot},
+              INNER_LAMBDA(int K) { Arr(IEdge, K) = 0._Real; });
+       });
+}
+
+//------------------------------------------------------------------------------
+// Enforce masking on a cell field after IC or restart read.
+void VertCoord::applyCellLayerMask(Array2DReal Arr, I4 NCellsAll) const {
+   OMEGA_SCOPE(LocMinLayerCell, MinLayerCell);
+   OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell);
+   I4 LocNVertLayers = NVertLayers;
+   parallelForOuter(
+       {NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) {
+          const int KMin = LocMinLayerCell(ICell);
+          const int KMax = LocMaxLayerCell(ICell);
+          parallelForInner(
+              Team, Range{0, LocNVertLayers - 1}, INNER_LAMBDA(int K) {
+                 if (K < KMin || K > KMax)
+                    Arr(ICell, K) = FillValueReal;
+                 // else: active layer - keep IC/restart value
+              });
+       });
+}
+
+//------------------------------------------------------------------------------
+// Enforce masking on a vertex field after IC or restart read.
+void VertCoord::applyVertexLayerMask(Array2DReal Arr, I4 NVerticesAll) const {
+   OMEGA_SCOPE(LocMinLayerVertexTop, MinLayerVertexTop);
+   OMEGA_SCOPE(LocMaxLayerVertexBot, MaxLayerVertexBot);
+   I4 LocNVertLayers = NVertLayers;
+   parallelForOuter(
+       {NVerticesAll}, KOKKOS_LAMBDA(int IVertex, const TeamMember &Team) {
+          const int KMin = LocMinLayerVertexTop(IVertex);
+          const int KMax = LocMaxLayerVertexBot(IVertex);
+          parallelForInner(
+              Team, Range{0, LocNVertLayers - 1}, INNER_LAMBDA(int K) {
+                 if (K < KMin || K > KMax)
+                    Arr(IVertex, K) = FillValueReal;
+                 // else: active layer - keep IC/restart value
+                 // Unlike edges, boundary vertices (those with one or more
+                 // active surrounding cells) hold valid, generally non-zero
+                 // data and are not zeroed.
+              });
+       });
+}
+
+//------------------------------------------------------------------------------
+// Enforce 3-zone masking on an edge field after IC or restart read.
+void VertCoord::applyEdgeLayerMask(Array2DReal Arr, I4 NEdgesAll) const {
+   OMEGA_SCOPE(LocMinLayerEdgeTop, MinLayerEdgeTop);
+   OMEGA_SCOPE(LocMaxLayerEdgeBot, MaxLayerEdgeBot);
+   OMEGA_SCOPE(LocMinLayerEdgeBot, MinLayerEdgeBot);
+   OMEGA_SCOPE(LocMaxLayerEdgeTop, MaxLayerEdgeTop);
+   I4 LocNVertLayers = NVertLayers;
+   parallelForOuter(
+       {NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) {
+          const int KTop = LocMinLayerEdgeTop(IEdge);
+          const int KBot = Kokkos::max(0, LocMaxLayerEdgeBot(IEdge));
+          const int KMin = LocMinLayerEdgeBot(IEdge);
+          const int KMax = LocMaxLayerEdgeTop(IEdge);
+          parallelForInner(
+              Team, Range{0, LocNVertLayers - 1}, INNER_LAMBDA(int K) {
+                 if (K < KTop || K > KBot)
+                    Arr(IEdge, K) = FillValueReal;
+                 else if (K < KMin || K > KMax)
+                    Arr(IEdge, K) = 0._Real;
+                 // else: active layer - keep IC/restart value
+              });
+       });
+}
+
 } // end namespace OMEGA
 
 //===----------------------------------------------------------------------===//
diff --git a/components/omega/src/ocn/VertCoord.h b/components/omega/src/ocn/VertCoord.h
index 1e478a04c373..aa9324288860 100644
--- a/components/omega/src/ocn/VertCoord.h
+++ b/components/omega/src/ocn/VertCoord.h
@@ -165,12 +165,15 @@ class VertCoord {
 
    HostArray2DReal BoundaryVertexH; ///< Mask to determine boundary vertex
 
-   // BottomGeomDepth read from mesh file
+   // BottomGeomDepth read from initial vertical coordinate file
    Array1DReal BottomGeomDepth;
 
    HostArray1DReal BottomGeomDepthH;
 
-   // TODO: Temporary handling of SurfacePressure
+   // Relative pressure (gauge pressure) at the top of the ocean column. Read
+   // from the initial-state/restart file via the "State"/"Restart" groups and
+   // used as the top boundary condition of the pressure field in
+   // computePressure().
    Array1DReal SurfacePressure;
 
    HostArray1DReal SurfacePressureH;
@@ -194,8 +197,9 @@ class VertCoord {
    std::string GeomZMidFldName; ///< Field name for midpoint geometric height
    std::string GeopotFldName;   ///< Field name for geopotential
    std::string
-       PseudoThicknessTargetFldName; ///< Field name for target thickness
-   std::string SshFldName;           ///< Field name for sea surface height
+       PseudoThicknessTargetFldName;   ///< Field name for target thickness
+   std::string SshFldName;             ///< Field name for sea surface height
+   std::string SurfacePressureFldName; ///< Field name for SurfacePressure
 
    // methods
 
@@ -216,6 +220,10 @@ class VertCoord {
    /// Copy member arrays from device to host
    void copyToHost();
 
+   /// Exchange halo and copy SurfacePressure to host after the initial-state
+   /// or restart stream has been read
+   void initSurfacePressure(Halo *MeshHalo);
+
    /// Copy member arrays from host to device
    void copyToDevice();
 
@@ -242,11 +250,40 @@ class VertCoord {
    /// Initialize computational masks
    void setMasks();
 
+   /// Zero all layers in [MinLayerEdgeTop, MaxLayerEdgeBot] of an edge field.
+   /// Call before computing a flux-type edge field so boundary edges show 0,
+   /// not FillValueReal. Layers outside the valid range retain their fill
+   /// value.
+   void zeroEdgeField(Array2DReal Arr, I4 NEdgesAll) const;
+
+   /// Enforce 3-zone masking on an edge field after IC or restart read:
+   ///   K < MinLayerEdgeTop or K > MaxLayerEdgeBot -> FillValueReal
+   ///   [MinLayerEdgeTop, MaxLayerEdgeBot] but outside
+   ///   [MinLayerEdgeBot, MaxLayerEdgeTop] -> 0
+   ///   [MinLayerEdgeBot, MaxLayerEdgeTop] -> unchanged (IC/restart values)
+   void applyEdgeLayerMask(Array2DReal Arr, I4 NEdgesAll) const;
+
+   /// Enforce masking on a cell field after IC or restart read:
+   ///   K < MinLayerCell(ICell) or K > MaxLayerCell(ICell) -> FillValueReal
+   ///   [MinLayerCell(ICell), MaxLayerCell(ICell)] -> unchanged (IC/restart
+   ///   values)
+   void applyCellLayerMask(Array2DReal Arr, I4 NCellsAll) const;
+
+   /// Enforce masking on a vertex field after IC or restart read. Unlike
+   /// edges, a boundary vertex with one or more active surrounding cells holds
+   /// valid (generally non-zero) data, so there is no zeroed boundary zone:
+   ///   K < MinLayerVertexTop or K > MaxLayerVertexBot -> FillValueReal
+   ///   [MinLayerVertexTop, MaxLayerVertexBot] -> unchanged (IC/restart values)
+   void applyVertexLayerMask(Array2DReal Arr, I4 NVerticesAll) const;
+
    /// Sums the mass thickness times g from the top layer down, starting with
    /// the surface pressure
+   /// Computes the relative pressure (gauge pressure, i.e., absolute pressure
+   /// minus the standard atmosphere) at layer interfaces and midpoints by
+   /// summing mass thickness times g from the top layer down.
    void computePressure(
        const Array2DReal &PseudoThickness, ///< [in] pseudo-thickness
-       const Array1DReal &SurfacePressure  ///< [in] surface pressure
+       const Array1DReal &SurfacePressure  ///< [in] relative surface pressure
    );
 
    /// Sum the mass thickness times specific volume from the bottom layer up,
diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp
index 815dab5e6814..2b988e8ddb19 100644
--- a/components/omega/src/ocn/VertMix.cpp
+++ b/components/omega/src/ocn/VertMix.cpp
@@ -332,8 +332,7 @@ void VertMix::defineFields() {
    }
 
    /// Create fields for state variables
-   const Real FillValue = -9.99e30;
-   int NDims            = 2;
+   int NDims = 2;
    std::vector<std::string> DimNames(NDims);
    DimNames[0] = "NCells";
    DimNames[1] = "NVertLayersP1";
@@ -347,9 +346,8 @@ void VertMix::defineFields() {
                      "vertical_diffusivity",           // CF-ish Name
                      0.0,                              // Min valid value
                      std::numeric_limits<Real>::max(), // Max valid value
-                     FillValue, // Scalar used for undefined entries
-                     NDims,     // Number of dimensions
-                     DimNames   // Dimension names
+                     NDims,                            // Number of dimensions
+                     DimNames                          // Dimension names
        );
    /// Create and register the VertVisc field
    auto VertViscField =
@@ -360,9 +358,8 @@ void VertMix::defineFields() {
                      "vertical_viscosity",             // CF-ish Name
                      0.0,                              // Min valid value
                      std::numeric_limits<Real>::max(), // Max valid value
-                     FillValue, // Scalar used for undefined entried
-                     NDims,     // Number of dimensions
-                     DimNames   // Dimension names
+                     NDims,                            // Number of dimensions
+                     DimNames                          // Dimension names
        );
    /// Create and register the GradRichNum field
    auto GradRichNumField =
@@ -372,9 +369,8 @@ void VertMix::defineFields() {
                      "sea_water_gradient_richardson_number", // CF-ish Name
                      std::numeric_limits<Real>::min(),       // Min valid value
                      std::numeric_limits<Real>::max(),       // Max valid value
-                     FillValue, // Scalar used for undefined entries
-                     NDims,     // Number of dimensions
-                     DimNames   // Dimension names
+                     NDims,   // Number of dimensions
+                     DimNames // Dimension names
        );
    /// Create and register the GradRichNumSmoothed field
    auto GradRichNumSmoothedField = Field::create(
@@ -384,9 +380,8 @@ void VertMix::defineFields() {
        "sea_water_gradient_richardson_number_smoothed", // CF-ish Name
        std::numeric_limits<Real>::min(),                // Min valid value
        std::numeric_limits<Real>::max(),                // Max valid value
-       FillValue, // Scalar used for undefined entries
-       NDims,     // Number of dimensions
-       DimNames   // Dimension names
+       NDims,                                           // Number of dimensions
+       DimNames                                         // Dimension names
    );
 
    // Create a field group for the vertmix-specific state fields
diff --git a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp
index 5e5bb08d4714..5347699acdcc 100644
--- a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp
+++ b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp
@@ -23,8 +23,7 @@ void KineticAuxVars::registerFields(
 ) const {
 
    // Create fields
-   const Real FillValue = -9.99e30;
-   int NDims            = 2;
+   int NDims = 2;
    std::vector<std::string> DimNames(NDims);
    std::string DimSuffix;
    if (MeshName == "Default") {
@@ -43,9 +42,8 @@ void KineticAuxVars::registerFields(
        "specific_kinetic_energy_of_sea_water",           // CF standard Name
        0,                                                // min valid value
        std::numeric_limits<Real>::max(),                 // max valid value
-       FillValue, // scalar for undefined entries
-       2,         // number of dimensions
-       DimNames   // dim names
+       2,                                                // number of dimensions
+       DimNames                                          // dim names
    );
 
    // Velocity divergence on cells
@@ -56,7 +54,6 @@ void KineticAuxVars::registerFields(
        "",                                  // CF standard Name
        std::numeric_limits<Real>::min(),    // min valid value
        std::numeric_limits<Real>::max(),    // max valid value
-       FillValue,                           // scalar used for undefined entries
        NDims,                               // number of dimensions
        DimNames                             // dimension names
    );
diff --git a/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp
index f22e80968447..2a608c176dc1 100644
--- a/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp
+++ b/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp
@@ -25,8 +25,7 @@ void PseudoThicknessAuxVars::registerFields(const std::string &AuxGroupName,
                                             const std::string &MeshName) const {
 
    // Create/define fields
-   const Real FillValue = -9.99e30;
-   int NDims            = 2;
+   int NDims = 2;
    std::vector<std::string> DimNames(NDims);
    std::string DimSuffix;
    if (MeshName == "Default") {
@@ -47,7 +46,6 @@ void PseudoThicknessAuxVars::registerFields(const std::string &AuxGroupName,
        "",                                       // CF standard Name
        0,                                        // min valid value
        std::numeric_limits<Real>::max(),         // max valid value
-       FillValue,                                // scalar for undefined entries
        NDims,                                    // number of dimensions
        DimNames                                  // dimension names
    );
@@ -61,7 +59,6 @@ void PseudoThicknessAuxVars::registerFields(const std::string &AuxGroupName,
        "",                               // CF standard Name
        0,                                // min valid value
        std::numeric_limits<Real>::max(), // max valid value
-       FillValue,                        // scalar used for undefined entries
        NDims,                            // number of dimensions
        DimNames                          // dimension names
    );
@@ -74,7 +71,6 @@ void PseudoThicknessAuxVars::registerFields(const std::string &AuxGroupName,
        "",                                       // CF standard Name
        0,                                        // min valid value
        std::numeric_limits<Real>::max(),         // max valid value
-       FillValue,                                // scalar for undefined entries
        NDims,                                    // number of dimensions
        DimNames                                  // dimension names
    );
diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp
index 191c22e05043..cc075b8ab634 100644
--- a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp
+++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp
@@ -19,8 +19,7 @@ void SurfTracerRestAuxVars::registerFields(
 ) const {
 
    // Create fields
-   const Real FillValue = -9.99e30;
-   int NDims            = 2;
+   int NDims = 2;
    std::vector<std::string> DimNames(NDims);
    std::string DimSuffix;
    if (MeshName == "Default") {
@@ -39,7 +38,6 @@ void SurfTracerRestAuxVars::registerFields(
                      "",                                   // CF standard Name
                      std::numeric_limits<Real>::min(),     // min valid value
                      std::numeric_limits<Real>::max(),     // max valid value
-                     FillValue, // scalar for undefined entries
                      NDims,     // number of dimensions
                      DimNames); // dim names
 
diff --git a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp
index 9c2645441877..6cbaffdf4dce 100644
--- a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp
+++ b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp
@@ -21,8 +21,7 @@ void TracerAuxVars::registerFields(const std::string &AuxGroupName,
                                    const std::string &MeshName) const {
 
    // Create fields
-   const Real FillValue = -9.99e30;
-   int NDims            = 3;
+   int NDims = 3;
    std::vector<std::string> DimNames(NDims);
    std::string DimSuffix;
    if (MeshName == "Default") {
@@ -43,7 +42,6 @@ void TracerAuxVars::registerFields(const std::string &AuxGroupName,
        "",                               // CF standard name
        std::numeric_limits<Real>::min(), // min valid value
        std::numeric_limits<Real>::max(), // max valid value
-       FillValue,                        // scalar for undefined entries
        NDims,                            // number of dimensions
        DimNames                          // dimension names
    );
diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp
index ffb8046ad5af..e67da008c01b 100644
--- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp
+++ b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp
@@ -24,16 +24,17 @@ VelocityDel2AuxVars::VelocityDel2AuxVars(const std::string &AuxStateSuffix,
       AreaTriangle(Mesh->AreaTriangle), VertexDegree(Mesh->VertexDegree),
       MinLayerEdgeBot(VCoord->MinLayerEdgeBot),
       MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop),
-      MinLayerVertexBot(VCoord->MinLayerVertexBot),
-      MaxLayerVertexTop(VCoord->MaxLayerVertexTop),
+      MinLayerEdgeTop(VCoord->MinLayerEdgeTop),
+      MaxLayerEdgeBot(VCoord->MaxLayerEdgeBot),
+      MinLayerVertexTop(VCoord->MinLayerVertexTop),
+      MaxLayerVertexBot(VCoord->MaxLayerVertexBot),
       MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {}
 
 void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName,
                                          const std::string &MeshName) const {
 
    // Create/define fields
-   const Real FillValue = -9.99e30;
-   int NDims            = 2;
+   int NDims = 2;
    std::vector<std::string> DimNames(NDims);
    DimNames[1] = "NVertLayers";
    std::string DimSuffix;
@@ -52,7 +53,6 @@ void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName,
        "",                                          // CF standard Name
        std::numeric_limits<Real>::min(),            // min valid value
        std::numeric_limits<Real>::max(),            // max valid value
-       FillValue,                                   // scalar for undef entries
        NDims,                                       // number of dimensions
        DimNames                                     // dimension names
    );
@@ -67,9 +67,8 @@ void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName,
                      "",          // CF standard Name
                      std::numeric_limits<Real>::min(), // min valid value
                      std::numeric_limits<Real>::max(), // max valid value
-                     FillValue, // scalar used for undefined entries
-                     NDims,     // number of dimensions
-                     DimNames   // dimension names
+                     NDims,                            // number of dimensions
+                     DimNames                          // dimension names
        );
 
    // Del2 of relative vorticity on vertices
@@ -81,9 +80,8 @@ void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName,
        "",                                            // CF standard Name
        std::numeric_limits<Real>::min(),              // min valid value
        std::numeric_limits<Real>::max(),              // max valid value
-       FillValue, // scalar used for undefined entries
-       NDims,     // number of dimensions
-       DimNames   // dimension names
+       NDims,                                         // number of dimensions
+       DimNames                                       // dimension names
    );
 
    // Add fields to Aux Field group
diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h
index b6f05f9c0d8b..3e9e57782536 100644
--- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h
+++ b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h
@@ -75,8 +75,17 @@ class VelocityDel2AuxVars {
    }
 
    KOKKOS_FUNCTION void computeVarsOnVertex(int IVertex, int KChunk) const {
-      const int KStart = chunkStart(KChunk, MinLayerVertexBot(IVertex));
-      const int KLen = chunkLength(KChunk, KStart, MaxLayerVertexTop(IVertex));
+      // Compute over the full vertex valid range [MinLayerVertexTop,
+      // MaxLayerVertexBot] so that boundary-vertex layers (where only some
+      // surrounding cells are active) receive a valid value. Each edge's
+      // contribution is clamped to that edge's valid range [MinLayerEdgeTop,
+      // MaxLayerEdgeBot], where Del2Edge has been computed or zeroed; this
+      // matches VorticityAuxVars::computeVarsOnVertex and avoids reading
+      // uninitialized (fill-value) layers of Del2Edge for deeper edges.
+      const int KStartVertex = chunkStart(KChunk, MinLayerVertexTop(IVertex));
+      const int KLenVertex =
+          chunkLength(KChunk, KStartVertex, MaxLayerVertexBot(IVertex));
+      const int KEndVertex = KStartVertex + KLenVertex - 1;
 
       const Real InvAreaTriangle = 1._Real / AreaTriangle(IVertex);
 
@@ -84,16 +93,19 @@ class VelocityDel2AuxVars {
 
       for (int J = 0; J < VertexDegree; ++J) {
          const int JEdge = EdgesOnVertex(IVertex, J);
-         for (int KVec = 0; KVec < KLen; ++KVec) {
-            const int K = KStart + KVec;
+         const int KStartEdge =
+             Kokkos::max(KStartVertex, MinLayerEdgeTop(JEdge));
+         const int KEndEdge = Kokkos::min(KEndVertex, MaxLayerEdgeBot(JEdge));
+         for (int K = KStartEdge; K <= KEndEdge; ++K) {
+            const int KVec = K - KStartVertex;
             Del2RelVortVertexTmp[KVec] += InvAreaTriangle * DcEdge(JEdge) *
                                           EdgeSignOnVertex(IVertex, J) *
                                           Del2Edge(JEdge, K);
          }
       }
 
-      for (int KVec = 0; KVec < KLen; ++KVec) {
-         const int K                   = KStart + KVec;
+      for (int KVec = 0; KVec < KLenVertex; ++KVec) {
+         const int K                   = KStartVertex + KVec;
          Del2RelVortVertex(IVertex, K) = Del2RelVortVertexTmp[KVec];
       }
    }
@@ -118,8 +130,10 @@ class VelocityDel2AuxVars {
    I4 VertexDegree;
    Array1DI4 MinLayerEdgeBot;
    Array1DI4 MaxLayerEdgeTop;
-   Array1DI4 MinLayerVertexBot;
-   Array1DI4 MaxLayerVertexTop;
+   Array1DI4 MinLayerEdgeTop;
+   Array1DI4 MaxLayerEdgeBot;
+   Array1DI4 MinLayerVertexTop;
+   Array1DI4 MaxLayerVertexBot;
    Array1DI4 MinLayerCell;
    Array1DI4 MaxLayerCell;
 };
diff --git a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp
index fabb2d53e0a1..062aeac80cf6 100644
--- a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp
+++ b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp
@@ -35,8 +35,7 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName,
                                       const std::string &MeshName) const {
 
    // Create fields with metadata
-   const Real FillValue = -9.99e30;
-   int NDims            = 2;
+   int NDims = 2;
    std::vector<std::string> DimNames(NDims);
    std::string DimSuffix;
    if (MeshName == "Default") {
@@ -56,9 +55,8 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName,
        "ocean_relative_vorticity",                         // CF standard Name
        std::numeric_limits<Real>::min(),                   // min valid value
        std::numeric_limits<Real>::max(),                   // max valid value
-       FillValue, // scalar for undefined entries
-       NDims,     // number of dimensions
-       DimNames   // dimension names
+       NDims,   // number of dimensions
+       DimNames // dimension names
    );
 
    // Normalized relative vorticity on vertices
@@ -69,7 +67,6 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName,
        "",                               // CF standard Name
        std::numeric_limits<Real>::min(), // min valid value
        std::numeric_limits<Real>::max(), // max valid value
-       FillValue,                        // scalar used for undefined entries
        NDims,                            // number of dimensions
        DimNames                          // dimension names
    );
@@ -83,7 +80,6 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName,
        "",                               // CF standard Name
        std::numeric_limits<Real>::min(), // min valid value
        std::numeric_limits<Real>::max(), // max valid value
-       FillValue,                        // scalar used for undefined entries
        NDims,                            // number of dimensions
        DimNames                          // dimension names
    );
@@ -98,7 +94,6 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName,
        "",                               // CF standard Name
        std::numeric_limits<Real>::min(), // min valid value
        std::numeric_limits<Real>::max(), // max valid value
-       FillValue,                        // scalar used for undefined entries
        NDims,                            // number of dimensions
        DimNames                          // dimension names
    );
@@ -112,9 +107,8 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName,
        "",                                           // CF standard Name
        std::numeric_limits<Real>::min(),             // min valid value
        std::numeric_limits<Real>::max(),             // max valid value
-       FillValue, // scalar used for undefined entries
-       NDims,     // number of dimensions
-       DimNames   // dimension names
+       NDims,                                        // number of dimensions
+       DimNames                                      // dimension names
    );
 
    // Add fields to Aux field group
diff --git a/components/omega/src/ocn/forcingVars/SfcStressForcingVars.cpp b/components/omega/src/ocn/forcingVars/SfcStressForcingVars.cpp
index 4375f79f63e3..f61cee3ccc26 100644
--- a/components/omega/src/ocn/forcingVars/SfcStressForcingVars.cpp
+++ b/components/omega/src/ocn/forcingVars/SfcStressForcingVars.cpp
@@ -36,7 +36,6 @@ void SfcStressForcingVars::registerFields(
                      "",                               // CF standard Name
                      std::numeric_limits<Real>::min(), // min valid value
                      std::numeric_limits<Real>::max(), // max valid value
-                     FillValue,                        // scalar for undefined
                      NDims,                            // number of dimensions
                      DimNames                          // dim names
        );
@@ -48,7 +47,6 @@ void SfcStressForcingVars::registerFields(
                      "",                          // CF standard Name
                      std::numeric_limits<Real>::min(), // min valid value
                      std::numeric_limits<Real>::max(), // max valid value
-                     FillValue,                        // scalar used undefined
                      NDims,                            // number of dimensions
                      DimNames                          // dimension names
        );
diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp
index 22a6c71a843b..334e3bb3a0fe 100644
--- a/components/omega/src/timeStepping/TimeStepper.cpp
+++ b/components/omega/src/timeStepping/TimeStepper.cpp
@@ -501,6 +501,10 @@ void TimeStepper::updateVelocityByTend(OceanState *State1, int TimeLevel1,
                                         CoeffSeconds * NormalVelTend(IEdge, K);
               });
        });
+   // Zero boundary layers; updateVelocityByTend only writes [KMin,KMax], so
+   // layers outside that range retain fill values that would corrupt
+   // tendencies.
+   VCoord->applyEdgeLayerMask(NormalVel1, Mesh->NEdgesAll);
 }
 
 //------------------------------------------------------------------------------
diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt
index b30760fea72b..adb8d8fda1d3 100644
--- a/components/omega/test/CMakeLists.txt
+++ b/components/omega/test/CMakeLists.txt
@@ -523,6 +523,17 @@ add_omega_test(
     "-n;8"
 )
 
+##################
+# Fill Value test
+##################
+
+add_omega_test(
+    FILL_VALUE_TEST
+    testFillValue.exe
+    ocn/FillValueTest.cpp
+    "-n;8"
+)
+
 ##########################
 # Vert Mix test
 ##########################
diff --git a/components/omega/test/infra/FieldTest.cpp b/components/omega/test/infra/FieldTest.cpp
index f3032cea8a77..9be90e8038ec 100644
--- a/components/omega/test/infra/FieldTest.cpp
+++ b/components/omega/test/infra/FieldTest.cpp
@@ -158,44 +158,44 @@ void initFieldTest() {
    DimNames[0] = "NCells";
    auto Test1DI4H =
        Field::create("Test1DI4H", "Test 1DI4 field on host", "Units1DI4H",
-                     "var_name_1DI4", 0, 100000, 999, 1, DimNames);
+                     "var_name_1DI4", 0, 100000, 1, DimNames);
 
    DimNames[0] = "NEdges";
    auto Test1DI8H =
        Field::create("Test1DI8H", "Test 1DI8 field on host", "Units1DI8H",
-                     "var_name_1DI8", 0, 100000, 999, 1, DimNames);
+                     "var_name_1DI8", 0, 100000, 1, DimNames);
 
    DimNames[0] = "NVertices";
    auto Test1DR4H =
        Field::create("Test1DR4H", "Test 1DR4 field on host", "Units1DR4H",
-                     "var_name_1DR4", 0.0, 100000.0, 999, 1, DimNames);
+                     "var_name_1DR4", 0.0, 100000.0, 1, DimNames);
 
    DimNames[0] = "NVertLayers";
    auto Test1DR8H =
        Field::create("Test1DR8H", "Test 1DR8 field on host", "Units1DR8H",
-                     "var_name_1DR8", 0.0, 100000.0, 999, 1, DimNames);
+                     "var_name_1DR8", 0.0, 100000.0, 1, DimNames);
 
    // 1D Fields on device
 
    DimNames[0] = "NCells";
    auto Test1DI4 =
        Field::create("Test1DI4", "Test 1DI4 field on device", "Units1DI4",
-                     "var_name_1DI4", 0, 100000, 999, 1, DimNames);
+                     "var_name_1DI4", 0, 100000, 1, DimNames);
 
    DimNames[0] = "NEdges";
    auto Test1DI8 =
        Field::create("Test1DI8", "Test 1DI8 field on device", "Units1DI8",
-                     "var_name_1DI8", 0, 100000, 999, 1, DimNames);
+                     "var_name_1DI8", 0, 100000, 1, DimNames);
 
    DimNames[0] = "NVertices";
    auto Test1DR4 =
        Field::create("Test1DR4", "Test 1DR4 field on device", "Units1DR4",
-                     "var_name_1DR4", 0.0, 100000.0, 999, 1, DimNames);
+                     "var_name_1DR4", 0.0, 100000.0, 1, DimNames);
 
    DimNames[0] = "NVertLayers";
    auto Test1DR8 =
        Field::create("Test1DR8", "Test 1DR8 field on device", "Units1DR8",
-                     "var_name_1DR8", 0.0, 100000.0, 999, 1, DimNames);
+                     "var_name_1DR8", 0.0, 100000.0, 1, DimNames);
 
    // 2D Fields on host
    DimNames.resize(2);
@@ -204,25 +204,25 @@ void initFieldTest() {
    DimNames[1] = "NVertLayers";
    auto Test2DI4H =
        Field::create("Test2DI4H", "Test 2DI4 field on host", "Units2DI4H",
-                     "var_name_2DI4", 0, 100000, 999, 2, DimNames);
+                     "var_name_2DI4", 0, 100000, 2, DimNames);
 
    DimNames[0] = "NEdges";
    DimNames[1] = "NVertLayers";
    auto Test2DI8H =
        Field::create("Test2DI8H", "Test 2DI8 field on host", "Units2DI8H",
-                     "var_name_2DI8", 0, 100000, 999, 2, DimNames);
+                     "var_name_2DI8", 0, 100000, 2, DimNames);
 
    DimNames[0] = "NVertices";
    DimNames[1] = "NVertLayers";
    auto Test2DR4H =
        Field::create("Test2DR4H", "Test 2DR4 field on host", "Units2DR4H",
-                     "var_name_2DR4", 0.0, 100000.0, 999, 2, DimNames);
+                     "var_name_2DR4", 0.0, 100000.0, 2, DimNames);
 
    DimNames[0] = "NCells";
    DimNames[1] = "NVertLayers";
    auto Test2DR8H =
        Field::create("Test2DR8H", "Test 2DR8 field on host", "Units2DR8H",
-                     "var_name_2DR8", 0.0, 100000.0, 999, 2, DimNames);
+                     "var_name_2DR8", 0.0, 100000.0, 2, DimNames);
 
    // 2D Fields on device
 
@@ -230,25 +230,25 @@ void initFieldTest() {
    DimNames[1] = "NVertLayers";
    auto Test2DI4 =
        Field::create("Test2DI4", "Test 2DI4 field on device", "Units2DI4",
-                     "var_name_2DI4", 0, 100000, 999, 2, DimNames);
+                     "var_name_2DI4", 0, 100000, 2, DimNames);
 
    DimNames[0] = "NEdges";
    DimNames[1] = "NVertLayers";
    auto Test2DI8 =
        Field::create("Test2DI8", "Test 2DI8 field on device", "Units2DI8",
-                     "var_name_2DI8", 0, 100000, 999, 2, DimNames);
+                     "var_name_2DI8", 0, 100000, 2, DimNames);
 
    DimNames[0] = "NVertices";
    DimNames[1] = "NVertLayers";
    auto Test2DR4 =
        Field::create("Test2DR4", "Test 2DR4 field on device", "Units2DR4",
-                     "var_name_2DR4", 0.0, 100000.0, 999, 2, DimNames);
+                     "var_name_2DR4", 0.0, 100000.0, 2, DimNames);
 
    DimNames[0] = "NEdges";
    DimNames[1] = "NVertLayers";
    auto Test2DR8 =
        Field::create("Test2DR8", "Test 2DR8 field on device", "Units2DR8",
-                     "var_name_2DR8", 0.0, 100000.0, 999, 2, DimNames);
+                     "var_name_2DR8", 0.0, 100000.0, 2, DimNames);
 
    // Higher dimension fields on device
 
@@ -258,7 +258,7 @@ void initFieldTest() {
    DimNames[2] = "NTracers";
    auto Test3DI4 =
        Field::create("Test3DI4", "Test 3DI4 field on device", "Units3DI4",
-                     "var_name_3DI4", 0, 100000, 999, 3, DimNames);
+                     "var_name_3DI4", 0, 100000, 3, DimNames);
 
    DimNames.resize(4);
    DimNames[0] = "NCells";
@@ -267,7 +267,7 @@ void initFieldTest() {
    DimNames[3] = "NTime";
    auto Test4DI8 =
        Field::create("Test4DI8", "Test 4DI8 field on device", "Units4DI8",
-                     "var_name_4DI8", 0, 100000, 999, 4, DimNames);
+                     "var_name_4DI8", 0, 100000, 4, DimNames);
 
    DimNames.resize(5);
    DimNames[0] = "NCells";
@@ -277,7 +277,7 @@ void initFieldTest() {
    DimNames[4] = "NStuff";
    auto Test5DR4 =
        Field::create("Test5DR4", "Test 5DR4 field on device", "Units5DR4",
-                     "var_name_5DR4", 0, 100000, 999, 5, DimNames);
+                     "var_name_5DR4", 0, 100000, 5, DimNames);
 
    /// Create two field groups for 1D and 2D fields
    auto FieldGroup1D = FieldGroup::create("FieldGroup1D");
@@ -322,6 +322,33 @@ void initFieldTest() {
    Array5DR4 Data5DR4("Test5DR4", NCellsSize, NVertLayers, NTracers, NTime,
                       NStuff);
 
+   // Attach arrays before populating so the fill-on-attach does not overwrite
+   // the reference values written below. Kokkos views share backing memory, so
+   // assignments through the local handle are visible via the stored pointer.
+   Test1DI4H->attachData<HostArray1DI4>(Data1DI4H);
+   Test1DI8H->attachData<HostArray1DI8>(Data1DI8H);
+   Test1DR4H->attachData<HostArray1DR4>(Data1DR4H);
+   Test1DR8H->attachData<HostArray1DR8>(Data1DR8H);
+
+   Field::attachFieldData<HostArray2DI4>("Test2DI4H", Data2DI4H);
+   Field::attachFieldData<HostArray2DI8>("Test2DI8H", Data2DI8H);
+   Field::attachFieldData<HostArray2DR4>("Test2DR4H", Data2DR4H);
+   Field::attachFieldData<HostArray2DR8>("Test2DR8H", Data2DR8H);
+
+   Field::attachFieldData("Test1DI4", Data1DI4);
+   Field::attachFieldData("Test1DI8", Data1DI8);
+   Field::attachFieldData("Test1DR4", Data1DR4);
+   Field::attachFieldData("Test1DR8", Data1DR8);
+
+   Test2DI4->attachData(Data2DI4);
+   Test2DI8->attachData(Data2DI8);
+   Test2DR4->attachData(Data2DR4);
+   Test2DR8->attachData(Data2DR8);
+
+   Test3DI4->attachData(Data3DI4);
+   Test4DI8->attachData(Data4DI8);
+   Test5DR4->attachData(Data5DR4);
+
    // Host arrays vertical vector
    for (int K = 0; K < NVertLayers; ++K) {
       Data1DR8H(K) = RefR8 + K;
@@ -388,33 +415,6 @@ void initFieldTest() {
           }
        });
 
-   // Attach data arrays
-   // Use member function for some and name interface for others
-
-   Test1DI4H->attachData<HostArray1DI4>(Data1DI4H);
-   Test1DI8H->attachData<HostArray1DI8>(Data1DI8H);
-   Test1DR4H->attachData<HostArray1DR4>(Data1DR4H);
-   Test1DR8H->attachData<HostArray1DR8>(Data1DR8H);
-
-   Field::attachFieldData<HostArray2DI4>("Test2DI4H", Data2DI4H);
-   Field::attachFieldData<HostArray2DI8>("Test2DI8H", Data2DI8H);
-   Field::attachFieldData<HostArray2DR4>("Test2DR4H", Data2DR4H);
-   Field::attachFieldData<HostArray2DR8>("Test2DR8H", Data2DR8H);
-
-   Field::attachFieldData("Test1DI4", Data1DI4);
-   Field::attachFieldData("Test1DI8", Data1DI8);
-   Field::attachFieldData("Test1DR4", Data1DR4);
-   Field::attachFieldData("Test1DR8", Data1DR8);
-
-   Test2DI4->attachData(Data2DI4);
-   Test2DI8->attachData(Data2DI8);
-   Test2DR4->attachData(Data2DR4);
-   Test2DR8->attachData(Data2DR8);
-
-   Test3DI4->attachData(Data3DI4);
-   Test4DI8->attachData(Data4DI8);
-   Test5DR4->attachData(Data5DR4);
-
    // End of init
 
 } // End initialization Fields
diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp
index 1678a08470f6..7052fa94e689 100644
--- a/components/omega/test/infra/IOStreamTest.cpp
+++ b/components/omega/test/infra/IOStreamTest.cpp
@@ -17,6 +17,7 @@
 #include "Dimension.h"
 #include "Error.h"
 #include "Field.h"
+#include "FillValues.h"
 #include "Forcing.h"
 #include "Halo.h"
 #include "HorzMesh.h"
@@ -171,11 +172,42 @@ int main(int argc, char **argv) {
       I4 NCellsOwned    = DefDecomp->NCellsOwned;
       Array1DI4 CellID  = DefDecomp->CellID;
 
+      // Register a field that is absent from the input file and mark it as an
+      // optional read, then add it to the InitialState stream. Reading the
+      // stream should succeed and leave the array at its fill value rather than
+      // failing because the variable is missing from the file.
+      Array1DReal OptTestData("OptionalTestField", NCellsSize);
+      auto OptTestField = Field::create("OptionalTestField", // field name
+                                        "optional read test field", // long name
+                                        "none",                     // units
+                                        "",         // CF standard name
+                                        -9.99E+10,  // min valid value
+                                        9.99E+10,   // max valid value
+                                        1,          // number of dimensions
+                                        {"NCells"}, // dimension names
+                                        false);     // not time dependent
+      OptTestField->setOptionalRead(true);
+      OptTestField->attachData<Array1DReal>(OptTestData);
+      IOStream::get("InitialState")->addField("OptionalTestField");
+
       // Read restart file for initial temperature and salinity data
       Metadata ReqMetadata; // leave empty for now - no required metadata
       Err = IOStream::read("InitialState", ModelClock, ReqMetadata);
       CHECK_ERROR_ABORT(Err, "IOStreamTest: Error reading initial state");
 
+      // The optional field was not in the file, so it should have been skipped
+      // and retain the fill value set by attachData on all owned cells.
+      I4 NOptFill     = 0;
+      auto OptReducer = Kokkos::Sum<I4>(NOptFill);
+      parallelReduce(
+          {NCellsOwned},
+          KOKKOS_LAMBDA(int Cell, I4 &Acc) {
+             if (OptTestData(Cell) == FillValueReal)
+                ++Acc;
+          },
+          OptReducer);
+      TestEval("Optional read fill retained", NOptFill, NCellsOwned, Err);
+
       // Overwrite salinity array with values associated with global cell
       // ID to test proper indexing of IO
       Array2DReal Test("Test", NCellsSize, NVertLayers);
diff --git a/components/omega/test/ocn/FillValueTest.cpp b/components/omega/test/ocn/FillValueTest.cpp
new file mode 100644
index 000000000000..d0b268a9e721
--- /dev/null
+++ b/components/omega/test/ocn/FillValueTest.cpp
@@ -0,0 +1,457 @@
+//===-- Test driver for Omega fill value standardization --------*- C++ -*-===//
+//
+/// \file
+/// \brief Test driver for Omega fill value constants and auto-fill behavior
+///
+/// Tests:
+///   1. FillValueI4/I8/R4/R8 constants match the NetCDF-C NC_FILL_* values
+///   2. Field::attachData() auto-fills the array with the declared fill value
+///   3. Inactive layers (k >= MaxLayerCell) contain FillValueReal after
+///      VertCoord initialization and compute
+///   4. NormalVelocity has the correct 3-zone pattern after IC read and
+///   masking:
+///      fully-inactive layers = FillValueReal, boundary layers = 0, active
+///      layers
+///      != FillValueReal (quiescent IC gives 0, which is valid)
+///   5. applyCellLayerMask sets inactive cell layers to FillValueReal and
+///      leaves active layers ([MinLayerCell, MaxLayerCell]) unchanged
+///   6. applyVertexLayerMask sets inactive vertex layers to FillValueReal and
+///      leaves active layers ([MinLayerVertexTop, MaxLayerVertexBot])
+///      unchanged (no zeroed boundary zone, unlike edges)
+//
+//===----------------------------------------------------------------------===//
+
+#include "AuxiliaryState.h"
+#include "DataTypes.h"
+#include "Decomp.h"
+#include "Dimension.h"
+#include "Error.h"
+#include "Field.h"
+#include "Halo.h"
+#include "HorzMesh.h"
+#include "IO.h"
+#include "IOStream.h"
+#include "Logging.h"
+#include "MachEnv.h"
+#include "OceanState.h"
+#include "OmegaKokkos.h"
+#include "PGrad.h"
+#include "Pacer.h"
+#include "Tendencies.h"
+#include "TimeMgr.h"
+#include "TimeStepper.h"
+#include "Tracers.h"
+#include "VertCoord.h"
+#include "mpi.h"
+
+#include <netcdf.h>
+
+using namespace OMEGA;
+
+//------------------------------------------------------------------------------
+void logNormalVelocityFillDiagnostics(const char *Stage) {
+
+   auto *DefEnv       = MachEnv::getDefault();
+   auto *DefMesh      = HorzMesh::getDefault();
+   auto *DefVertCoord = VertCoord::getDefault();
+   auto *DefState     = OceanState::getDefault();
+
+   I4 NEdgesOwned = DefMesh->NEdgesOwned;
+   I4 NVertLayers = DefVertCoord->NVertLayers;
+
+   auto NVH = createHostMirrorCopy(DefState->NormalVelocity[0]);
+   const auto &MinLayerEdgeTopH = DefVertCoord->MinLayerEdgeTopH;
+   const auto &MaxLayerEdgeBotH = DefVertCoord->MaxLayerEdgeBotH;
+   const auto &MinLayerEdgeBotH = DefVertCoord->MinLayerEdgeBotH;
+   const auto &MaxLayerEdgeTopH = DefVertCoord->MaxLayerEdgeTopH;
+
+   int ErrActiveLocal = 0;
+   int NLogged        = 0;
+   for (int IEdge = 0; IEdge < NEdgesOwned; ++IEdge) {
+      int KTop = MinLayerEdgeTopH(IEdge);
+      int KBot = std::max(0, (int)MaxLayerEdgeBotH(IEdge));
+      int KMin = MinLayerEdgeBotH(IEdge);
+      int KMax = MaxLayerEdgeTopH(IEdge);
+
+      for (int K = 0; K < NVertLayers; ++K) {
+         Real Val    = NVH(IEdge, K);
+         bool Active = (KMax >= 0 && K >= KMin && K <= KMax);
+
+         if (Active && Val == FillValueReal) {
+            ++ErrActiveLocal;
+            if (NLogged < 5) {
+               LOG_INFO("FillValueTest: {} active fill sample: IEdge={} K={} "
+                        "KTop={} KBot={} KMin={} KMax={} Val={}",
+                        Stage, IEdge, K, KTop, KBot, KMin, KMax, Val);
+               ++NLogged;
+            }
+         }
+      }
+   }
+
+   int ErrActiveGlobal = 0;
+   MPI_Allreduce(&ErrActiveLocal, &ErrActiveGlobal, 1, MPI_INT, MPI_SUM,
+                 DefEnv->getComm());
+   LOG_INFO("FillValueTest: {} active NormalVelocity fill count: local={} "
+            "global={}",
+            Stage, ErrActiveLocal, ErrActiveGlobal);
+}
+
+//------------------------------------------------------------------------------
+// Full initialization needed for OceanState (matches StateTest.cpp pattern)
+void initFillValueTest() {
+
+   MachEnv::init(MPI_COMM_WORLD);
+   MachEnv *DefEnv  = MachEnv::getDefault();
+   MPI_Comm DefComm = DefEnv->getComm();
+
+   initLogging(DefEnv);
+
+   Config("Omega");
+   Config::readAll("omega.yml");
+
+   TimeStepper::init1();
+   TimeStepper *DefStepper = TimeStepper::getDefault();
+   Clock *ModelClock       = DefStepper->getClock();
+
+   IO::init(DefComm);
+
+   Field::init(ModelClock);
+   IOStream::init(ModelClock);
+
+   Decomp::init();
+
+   int Err = Halo::init();
+   if (Err != 0)
+      ABORT_ERROR("FillValueTest: error initializing default halo");
+
+   HorzMesh::init(ModelClock);
+   VertCoord::init();
+   Tracers::init();
+   AuxiliaryState::init();
+   PressureGrad::init();
+   Tendencies::init();
+   TimeStepper::init2();
+
+   Err = OceanState::init();
+   if (Err != 0)
+      ABORT_ERROR("FillValueTest: error initializing OceanState");
+
+   bool StreamsValid = IOStream::validateAll();
+   if (!StreamsValid)
+      ABORT_ERROR("FillValueTest: error validating IO streams");
+}
+
+//------------------------------------------------------------------------------
+int main(int argc, char *argv[]) {
+
+   Error ErrAll;
+
+   MPI_Init(&argc, &argv);
+   Kokkos::initialize();
+   Pacer::initialize(MPI_COMM_WORLD);
+   Pacer::setPrefix("Omega:");
+   {
+      initFillValueTest();
+
+      // Read initial state so NormalVelocity has real values (quiescent IC:
+      // NV=0 for all active layers). Mirrors the pattern in ocnInit /
+      // StateTest.
+      {
+         Clock *ModelClock = TimeStepper::getDefault()->getClock();
+         Metadata ReqMeta;
+         Error Err1 = IOStream::read("InitialState", ModelClock, ReqMeta);
+         CHECK_ERROR_ABORT(Err1, "FillValueTest: error reading initial state");
+         logNormalVelocityFillDiagnostics("after InitialState read");
+
+         OceanState::getDefault()->exchangeHalo(0);
+         logNormalVelocityFillDiagnostics("after NormalVelocity halo exchange");
+
+         // Apply 3-zone mask after IC read: boundary layers -> 0,
+         // inactive layers -> FillValueReal, active layers unchanged.
+         VertCoord::getDefault()->applyEdgeLayerMask(
+             OceanState::getDefault()->NormalVelocity[0],
+             HorzMesh::getDefault()->NEdgesAll);
+         logNormalVelocityFillDiagnostics("after edge layer mask");
+
+         OceanState::getDefault()->copyToHost(0);
+      }
+
+      auto *DefMesh      = HorzMesh::getDefault();
+      auto *DefVertCoord = VertCoord::getDefault();
+      auto *DefState     = OceanState::getDefault();
+
+      // ------------------------------------------------------------------
+      // Test 1: FillValue constants match NetCDF-C NC_FILL_* values
+      // ------------------------------------------------------------------
+      {
+         int Err = 0;
+         if (FillValueI4 != static_cast<I4>(NC_FILL_INT))
+            ++Err;
+         if (FillValueI8 != static_cast<I8>(NC_FILL_INT64))
+            ++Err;
+         if (static_cast<double>(FillValueR4) !=
+             static_cast<double>(NC_FILL_FLOAT))
+            ++Err;
+         if (FillValueR8 != static_cast<R8>(NC_FILL_DOUBLE))
+            ++Err;
+
+         if (Err == 0) {
+            LOG_INFO("FillValueTest: fill constant values PASS");
+         } else {
+            ErrAll += Error(ErrorCode::Fail,
+                            "FillValueTest: fill constant values FAIL");
+         }
+      }
+
+      // ------------------------------------------------------------------
+      // Test 2: attachData() auto-fills the array with the declared fill value
+      // ------------------------------------------------------------------
+      {
+         // Create a small non-distributed test dimension and field
+         I4 NTest = 10;
+         Dimension::create("FVTestDim", NTest);
+         auto TestField =
+             Field::create("FVTestField", "fill value test", "1", "", -1.0e40,
+                           1.0e40, 1, {"FVTestDim"}, false);
+
+         // Allocate a host array set to 0.0 (not the fill value)
+         HostArray1DR8 TestArr("FVTestArr", NTest);
+         Kokkos::deep_copy(TestArr, 0.0);
+
+         // attachData() should automatically fill with FillValueR8
+         TestField->attachData<HostArray1DR8>(TestArr);
+
+         int Err = 0;
+         for (int I = 0; I < NTest; ++I) {
+            if (TestArr(I) != FillValueR8)
+               ++Err;
+         }
+
+         if (Err == 0) {
+            LOG_INFO("FillValueTest: attachData auto-fill PASS");
+         } else {
+            ErrAll += Error(ErrorCode::Fail,
+                            "FillValueTest: attachData auto-fill FAIL");
+         }
+      }
+
+      // ------------------------------------------------------------------
+      // Test 3: Inactive layers (k >= MaxLayerCell) contain FillValueReal
+      //         in a VertCoord field (GeomZMid) after initialization
+      // ------------------------------------------------------------------
+      {
+         I4 NCellsOwned = DefMesh->NCellsOwned;
+         I4 NVertLayers = DefVertCoord->NVertLayers;
+
+         auto GeomZMidH = createHostMirrorCopy(DefVertCoord->GeomZMid);
+
+         int Err = 0;
+         for (int ICell = 0; ICell < NCellsOwned; ++ICell) {
+            I4 KMax = DefVertCoord->MaxLayerCellH(ICell);
+            for (int K = KMax; K < NVertLayers; ++K) {
+               if (GeomZMidH(ICell, K) != FillValueReal)
+                  ++Err;
+            }
+         }
+
+         if (Err == 0) {
+            LOG_INFO("FillValueTest: inactive layers contain fill value PASS");
+         } else {
+            ErrAll +=
+                Error(ErrorCode::Fail,
+                      "FillValueTest: inactive layers contain fill value FAIL");
+         }
+      }
+
+      // ------------------------------------------------------------------
+      // Test 4: NormalVelocity has the correct 3-zone pattern after IC read
+      //         and applyEdgeLayerMask():
+      //   - Fully inactive layers (outside [MinLayerEdgeTop, MaxLayerEdgeBot]):
+      //     FillValueReal
+      //   - Boundary layers ([MinLayerEdgeTop, MaxLayerEdgeBot] but outside
+      //     [MinLayerEdgeBot, MaxLayerEdgeTop]): 0
+      //   - Active layers ([MinLayerEdgeBot, MaxLayerEdgeTop]):
+      //     != FillValueReal (quiescent IC gives 0, which is valid)
+      // ------------------------------------------------------------------
+      {
+         I4 NEdgesOwned = DefMesh->NEdgesOwned;
+         I4 NVertLayers = DefVertCoord->NVertLayers;
+
+         auto NVH = createHostMirrorCopy(DefState->NormalVelocity[0]);
+         const auto &MinLayerEdgeTopH = DefVertCoord->MinLayerEdgeTopH;
+         const auto &MaxLayerEdgeBotH = DefVertCoord->MaxLayerEdgeBotH;
+         const auto &MinLayerEdgeBotH = DefVertCoord->MinLayerEdgeBotH;
+         const auto &MaxLayerEdgeTopH = DefVertCoord->MaxLayerEdgeTopH;
+
+         int ErrInactive = 0; // non-fill in fully-inactive zone (bad)
+         int ErrBoundary = 0; // non-zero in boundary zone (bad)
+         int ErrActive   = 0; // fill value in active zone (bad)
+
+         for (int IEdge = 0; IEdge < NEdgesOwned; ++IEdge) {
+            int KTop = MinLayerEdgeTopH(IEdge);
+            int KBot = std::max(0, (int)MaxLayerEdgeBotH(IEdge));
+            int KMin = MinLayerEdgeBotH(IEdge);
+            int KMax = MaxLayerEdgeTopH(IEdge); // -1 for land-adjacent edges
+
+            for (int K = 0; K < NVertLayers; ++K) {
+               Real Val      = NVH(IEdge, K);
+               bool valid    = (K >= KTop && K <= KBot);
+               bool active   = (KMax >= 0 && K >= KMin && K <= KMax);
+               bool boundary = valid && !active;
+
+               if (!valid && Val != FillValueReal)
+                  ++ErrInactive;
+               if (boundary && Val != 0._Real)
+                  ++ErrBoundary;
+               if (active && Val == FillValueReal)
+                  ++ErrActive;
+            }
+         }
+
+         if (ErrInactive == 0 && ErrBoundary == 0 && ErrActive == 0) {
+            LOG_INFO(
+                "FillValueTest: NormalVelocity fill/zero/real pattern PASS");
+         } else {
+            ErrAll += Error(
+                ErrorCode::Fail,
+                "FillValueTest: NormalVelocity fill/zero/real pattern FAIL"
+                " (ErrInactive={}, ErrBoundary={}, ErrActive={})",
+                ErrInactive, ErrBoundary, ErrActive);
+         }
+      }
+
+      // ------------------------------------------------------------------
+      // Test 5: applyCellLayerMask sets inactive cell layers to FillValueReal
+      //         and leaves active layers ([MinLayerCell, MaxLayerCell])
+      //         unchanged.
+      //   A synthetic cell field is filled with a sentinel value (!= 0 and
+      //   != FillValueReal); after masking:
+      //   - Inactive layers (K < MinLayerCell or K > MaxLayerCell):
+      //     FillValueReal
+      //   - Active layers ([MinLayerCell, MaxLayerCell]): sentinel (unchanged)
+      // ------------------------------------------------------------------
+      {
+         I4 NCellsOwned = DefMesh->NCellsOwned;
+         I4 NVertLayers = DefVertCoord->NVertLayers;
+
+         const Real Sentinel = 42.0_Real;
+
+         Array2DReal SynthCell("SynthCell", DefMesh->NCellsSize, NVertLayers);
+         deepCopy(SynthCell, Sentinel);
+         DefVertCoord->applyCellLayerMask(SynthCell, DefMesh->NCellsAll);
+
+         auto SynthCellH           = createHostMirrorCopy(SynthCell);
+         const auto &MinLayerCellH = DefVertCoord->MinLayerCellH;
+         const auto &MaxLayerCellH = DefVertCoord->MaxLayerCellH;
+
+         int ErrInactive = 0; // non-fill in inactive zone (bad)
+         int ErrActive   = 0; // not sentinel in active zone (bad)
+
+         for (int ICell = 0; ICell < NCellsOwned; ++ICell) {
+            int KMin = MinLayerCellH(ICell);
+            int KMax = MaxLayerCellH(ICell); // -1 for fully-land cells
+
+            for (int K = 0; K < NVertLayers; ++K) {
+               Real Val    = SynthCellH(ICell, K);
+               bool active = (K >= KMin && K <= KMax);
+
+               if (active && Val != Sentinel)
+                  ++ErrActive;
+               if (!active && Val != FillValueReal)
+                  ++ErrInactive;
+            }
+         }
+
+         if (ErrInactive == 0 && ErrActive == 0) {
+            LOG_INFO("FillValueTest: applyCellLayerMask pattern PASS");
+         } else {
+            ErrAll += Error(ErrorCode::Fail,
+                            "FillValueTest: applyCellLayerMask pattern FAIL"
+                            " (ErrInactive={}, ErrActive={})",
+                            ErrInactive, ErrActive);
+         }
+      }
+
+      // ------------------------------------------------------------------
+      // Test 6: applyVertexLayerMask sets inactive vertex layers to
+      //         FillValueReal and leaves active layers
+      //         ([MinLayerVertexTop, MaxLayerVertexBot]) unchanged.
+      //   Unlike edges, a boundary vertex with one or more active surrounding
+      //   cells holds valid (generally non-zero) data, so there is no zeroed
+      //   boundary zone.
+      //   A synthetic vertex field is filled with a sentinel value (!= 0 and
+      //   != FillValueReal); after masking:
+      //   - Inactive layers (outside [MinLayerVertexTop, MaxLayerVertexBot]):
+      //     FillValueReal
+      //   - Active layers ([MinLayerVertexTop, MaxLayerVertexBot]): sentinel
+      //     (unchanged)
+      // ------------------------------------------------------------------
+      {
+         I4 NVerticesOwned = DefMesh->NVerticesOwned;
+         I4 NVertLayers    = DefVertCoord->NVertLayers;
+
+         const Real Sentinel = 42.0_Real;
+
+         Array2DReal SynthVtx("SynthVtx", DefMesh->NVerticesSize, NVertLayers);
+         deepCopy(SynthVtx, Sentinel);
+         DefVertCoord->applyVertexLayerMask(SynthVtx, DefMesh->NVerticesAll);
+
+         auto SynthVtxH                 = createHostMirrorCopy(SynthVtx);
+         const auto &MinLayerVertexTopH = DefVertCoord->MinLayerVertexTopH;
+         const auto &MaxLayerVertexBotH = DefVertCoord->MaxLayerVertexBotH;
+
+         int ErrInactive = 0; // non-fill in inactive zone (bad)
+         int ErrActive   = 0; // not sentinel in active zone (bad)
+
+         for (int IVertex = 0; IVertex < NVerticesOwned; ++IVertex) {
+            int KMin = MinLayerVertexTopH(IVertex);
+            int KMax = MaxLayerVertexBotH(IVertex); // -1 for invalid vertices
+
+            for (int K = 0; K < NVertLayers; ++K) {
+               Real Val    = SynthVtxH(IVertex, K);
+               bool active = (K >= KMin && K <= KMax);
+
+               if (active && Val != Sentinel)
+                  ++ErrActive;
+               if (!active && Val != FillValueReal)
+                  ++ErrInactive;
+            }
+         }
+
+         if (ErrInactive == 0 && ErrActive == 0) {
+            LOG_INFO("FillValueTest: applyVertexLayerMask pattern PASS");
+         } else {
+            ErrAll += Error(ErrorCode::Fail,
+                            "FillValueTest: applyVertexLayerMask pattern FAIL"
+                            " (ErrInactive={}, ErrActive={})",
+                            ErrInactive, ErrActive);
+         }
+      }
+
+      // ------------------------------------------------------------------
+      // Finalize
+      // ------------------------------------------------------------------
+      OceanState::clear();
+      Tendencies::clear();
+      PressureGrad::clear();
+      AuxiliaryState::clear();
+      Tracers::clear();
+      IOStream::finalize();
+      TimeStepper::clear();
+      HorzMesh::clear();
+      VertCoord::clear();
+      Field::clear();
+      Dimension::clear();
+      Halo::clear();
+      Decomp::clear();
+      MachEnv::removeAll();
+   }
+   Pacer::finalize();
+   Kokkos::finalize();
+   int RetVal = ErrAll.isFail() ? 1 : 0;
+   MPI_Finalize();
+
+   return RetVal;
+}
+//===----------------------------------------------------------------------===//
diff --git a/components/omega/test/ocn/StateTest.cpp b/components/omega/test/ocn/StateTest.cpp
index 252056ef2074..8ca90c95bba4 100644
--- a/components/omega/test/ocn/StateTest.cpp
+++ b/components/omega/test/ocn/StateTest.cpp
@@ -123,6 +123,10 @@ void initStateTest() {
    DefState->exchangeHalo(0);
    DefState->copyToHost(0);
 
+   // SurfacePressure is owned by VertCoord but read from the same InitialState
+   // stream; finish its initialization (halo exchange + copy to host).
+   VertCoord::getDefault()->initSurfacePressure(Halo::getDefault());
+
    return;
 }
 
@@ -262,6 +266,29 @@ int main(int argc, char *argv[]) {
          LOG_INFO("State: Out-of-range {} Non-zero: {}", Count2, Count1);
       }
 
+      // SurfacePressure is owned by VertCoord and an optional read from the
+      // InitialState stream. It is absent from the test input file, so it
+      // defaults to zero; verify its host mirror is sized and defaulted
+      // correctly.
+      VertCoord *DefVCoord = VertCoord::getDefault();
+      if (DefVCoord->SurfacePressureH.extent_int(0) == DefState->NCellsSize) {
+         LOG_INFO("State: SurfacePressureH size PASS");
+      } else {
+         RetVal += 1;
+         LOG_INFO("State: SurfacePressureH size FAIL");
+      }
+      int Count3 = 0;
+      for (int Cell = 0; Cell < NCellsAll; Cell++) {
+         if (DefVCoord->SurfacePressureH(Cell) != 0.0)
+            Count3++;
+      }
+      if (Count3 == 0) {
+         LOG_INFO("State: SurfacePressure read PASS");
+      } else {
+         RetVal += 1;
+         LOG_INFO("State: SurfacePressure read FAIL");
+      }
+
       // Test time swapping with 2 and higher numbers of time levels
       for (int NTimeLevels = 2; NTimeLevels < 5; NTimeLevels++) {
 
diff --git a/components/omega/test/ocn/VertCoordTest.cpp b/components/omega/test/ocn/VertCoordTest.cpp
index d9f6b0f94bb9..69a627617307 100644
--- a/components/omega/test/ocn/VertCoordTest.cpp
+++ b/components/omega/test/ocn/VertCoordTest.cpp
@@ -128,6 +128,16 @@ int main(int argc, char *argv[]) {
              Error(ErrorCode::Fail, "VertCoordTest: Bathy min/max test FAIL");
       }
 
+      // SurfacePressure is owned by VertCoord and its host mirror is sized on
+      // cells. Its value is read from the InitialState stream (exercised in
+      // StateTest); here we just verify the host mirror allocation.
+      if (DefVertCoord->SurfacePressureH.extent_int(0) == NCellsSize) {
+         LOG_INFO("VertCoordTest: SurfacePressureH size PASS");
+      } else {
+         ErrAll += Error(ErrorCode::Fail,
+                         "VertCoordTest: SurfacePressureH size FAIL");
+      }
+
       // Tests for computePressure
 
       Array2DReal PseudoThickness("PseudoThickness", NCellsSize, NVertLayers);