Fix fill values in all Omega fields

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

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.

components/omega/doc/devGuide/VertCoord.md

@@ -128,3 +128,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.

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

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.

components/omega/doc/userGuide/VertCoord.md

@@ -46,3 +46,12 @@ 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 |
+
+### 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.

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

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
 
 //===----------------------------------------------------------------------===//

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;

components/omega/src/infra/Field.h

@@ -19,13 +19,15 @@
 #include "DataTypes.h"
 #include "Dimension.h"
 #include "Error.h"
+#include "FillValues.h"
 #include "Logging.h"
 #include "OmegaKokkos.h"
 #include "TimeMgr.h"
 #include <any>
 #include <map>
 #include <memory>
 #include <set>
+#include <type_traits>
 
 namespace OMEGA {
 
@@ -82,6 +84,20 @@ class Field {
    /// 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 +128,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
@@ -303,7 +318,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 +330,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 +348,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 +359,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",