Support non-zero SurfacePressure

components/omega/doc/design/FillValues.md

@@ -221,6 +221,13 @@ Ordering invariants preserved by this design:
 - 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

components/omega/doc/devGuide/Field.md

@@ -192,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++

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

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` |
 |------------|-------------|

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
 

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.

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,6 +48,23 @@ 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
 

components/omega/src/infra/Field.cpp

@@ -125,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;
 
@@ -333,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.

components/omega/src/infra/Field.h

@@ -79,6 +79,11 @@ 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.
@@ -226,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
    //---------------------------------------------------------------------------

components/omega/src/infra/IOStream.cpp

@@ -1725,11 +1725,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

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,

components/omega/src/ocn/OceanInit.cpp

@@ -157,6 +157,7 @@ int ocnInit(MPI_Comm Comm ///< [in] ocean MPI communicator
    DefState->applyLayerMasks(CurTimeLevel);
 
    DefState->copyToHost(CurTimeLevel);
+   VertCoord::getDefault()->initSurfacePressure(Halo::getDefault());
 
    Forcing *DefForcing = Forcing::getDefault();
    DefForcing->exchangeHalo();

components/omega/src/ocn/OceanState.cpp

@@ -18,6 +18,7 @@
 #include "OmegaKokkos.h"
 #include "TimeStepper.h"
 #include "Tracers.h"
+#include "VertCoord.h"
 
 namespace OMEGA {
 
@@ -227,12 +228,18 @@ void OceanState::defineFields() {
                      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"))