From 804fb23041e90ecbc7b30dd9e96198ae203e2ff2 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 03:35:11 -0500 Subject: [PATCH 1/8] Correct comment on where BottomGeomDepth gets read from --- components/omega/src/ocn/VertCoord.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/omega/src/ocn/VertCoord.h b/components/omega/src/ocn/VertCoord.h index 977c931813ca..dd6a1cd9185f 100644 --- a/components/omega/src/ocn/VertCoord.h +++ b/components/omega/src/ocn/VertCoord.h @@ -165,7 +165,7 @@ 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; From 9c39d677b8a6fd879572b8e5f6a62fa95fa8b699 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 2 Jul 2026 09:34:43 -0500 Subject: [PATCH 2/8] Make VertCoord the owner of SurfacePressure SurfacePressure is the relative (gauge) pressure at the top of the ocean column and the top boundary condition of the pressure field, consumed only by VertCoord::computePressure. Give VertCoord ownership of the array (device and host mirror) instead of passing it a borrowed handle from OceanState: - VertCoord allocates and registers the SurfacePressure field and adds it to the State and Restart groups (creating them if VertCoord initializes before OceanState) so it is still read from the initial-state and restart files. - VertCoord::initSurfacePressure() exchanges the halo and copies to host after the stream read; OceanInit calls it in place of the old OceanState hook. - OceanState no longer owns SurfacePressure and guards its State-group creation since VertCoord may create the group first. Also clarify in Eos that the pressure argument is relative (gauge) pressure. Co-Authored-By: Claude Opus 4.8 --- components/omega/src/ocn/Eos.h | 15 ++-- components/omega/src/ocn/OceanInit.cpp | 1 + components/omega/src/ocn/OceanState.cpp | 11 ++- components/omega/src/ocn/VertCoord.cpp | 97 ++++++++++++++++++------- components/omega/src/ocn/VertCoord.h | 19 ++++- 5 files changed, 107 insertions(+), 36 deletions(-) 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/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index d1c7d5d227e9..dce3832a9d1e 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/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(); diff --git a/components/omega/src/ocn/OceanState.cpp b/components/omega/src/ocn/OceanState.cpp index 42764e97af36..f497d6c1a693 100644 --- a/components/omega/src/ocn/OceanState.cpp +++ b/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 StateGroup; + if (FieldGroup::exists(StateGroupName)) + StateGroup = FieldGroup::get(StateGroupName); + else + StateGroup = FieldGroup::create(StateGroupName); // Add restart group if needed if (!FieldGroup::exists("Restart")) diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp index 07cab988b895..4d4f0259a78a 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,6 +246,7 @@ void VertCoord::defineFields() { GeopotFldName.append(Name); PseudoThicknessTargetFldName.append(Name); SshFldName.append(Name); + SurfacePressureFldName.append(Name); } // Create fields for VertCoord variables @@ -296,6 +299,17 @@ void VertCoord::defineFields() { 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 + ); + NDims = 2; DimNames.resize(NDims); DimNames[1] = "NVertLayers"; @@ -334,14 +348,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::max(), // max valid value - 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::max(), // max valid value + NDims, // number of dimensions + DimNames // dimension names ); auto GeomZInterfaceField = Field::create( @@ -358,14 +372,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::max(), // max valid value - 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::max(), // max valid value + NDims, // number of dimensions + DimNames // dimension names ); auto GeomZMidField = Field::create( @@ -445,6 +459,24 @@ void VertCoord::defineFields() { PseudoThicknessTargetField->attachData(PseudoThicknessTarget); SshField->attachData(SshCell); + // Add SurfacePressure to the State and Restart groups so it is read from the + // initial-state and restart files. 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(SurfacePressure); + } // end defineFields //------------------------------------------------------------------------------ @@ -471,6 +503,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 //------------------------------------------------------------------------------ @@ -663,6 +701,14 @@ void VertCoord::setStreamArrays(const bool ReadStream, Halo *MeshHalo) { VertCoordMovementWeightsH = createHostMirrorCopy(VertCoordMovementWeights); } +//------------------------------------------------------------------------------ +// Exchange halo and copy SurfacePressure to host after the initial-state or +// restart stream has been read. +void VertCoord::initSurfacePressure(Halo *MeshHalo) { + MeshHalo->exchangeFullArrayHalo(SurfacePressure, OnCell); + deepCopy(SurfacePressureH, SurfacePressure); +} + //------------------------------------------------------------------------------ // Compute min and max layer indices for edges based on MinLayerCell and // MaxLayerCell @@ -958,11 +1004,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 diff --git a/components/omega/src/ocn/VertCoord.h b/components/omega/src/ocn/VertCoord.h index dd6a1cd9185f..aa9324288860 100644 --- a/components/omega/src/ocn/VertCoord.h +++ b/components/omega/src/ocn/VertCoord.h @@ -170,7 +170,10 @@ class VertCoord { 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(); @@ -270,9 +278,12 @@ class VertCoord { /// 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, From 5929785502ee89ef80c465c0097510dd03e96695 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 2 Jul 2026 09:35:14 -0500 Subject: [PATCH 3/8] Test SurfacePressure as a VertCoord field Retarget the SurfacePressure size and read checks in StateTest to the VertCoord-owned array and add a host-mirror size check to VertCoordTest. Co-Authored-By: Claude Opus 4.8 --- components/omega/test/ocn/StateTest.cpp | 26 +++++++++++++++++++++ components/omega/test/ocn/VertCoordTest.cpp | 10 ++++++++ 2 files changed, 36 insertions(+) diff --git a/components/omega/test/ocn/StateTest.cpp b/components/omega/test/ocn/StateTest.cpp index 09caff2e6d73..94b373298c90 100644 --- a/components/omega/test/ocn/StateTest.cpp +++ b/components/omega/test/ocn/StateTest.cpp @@ -131,6 +131,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; } @@ -270,6 +274,28 @@ int main(int argc, char *argv[]) { LOG_INFO("State: Out-of-range {} Non-zero: {}", Count2, Count1); } + // SurfacePressure is now owned by VertCoord but read from the same + // InitialState stream; verify its host mirror is sized and read + // 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); From 71b3600d49c27b46dcb99259bd1f58b92a688759 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 2 Jul 2026 09:35:40 -0500 Subject: [PATCH 4/8] Document SurfacePressure as a VertCoord field Move the SurfacePressure documentation from the OceanState user and developer guides to the VertCoord guides, describing that VertCoord owns it, that it is read from the initial-state and restart files, and that initSurfacePressure() finalizes it after the read. Note it is relative (gauge) pressure. Co-Authored-By: Claude Opus 4.8 --- components/omega/doc/devGuide/OceanState.md | 8 ++++++- components/omega/doc/devGuide/VertCoord.md | 21 +++++++++++++++++ components/omega/doc/userGuide/OceanState.md | 12 ++++++++-- components/omega/doc/userGuide/VertCoord.md | 24 +++++++++++++++++--- 4 files changed, 59 insertions(+), 6 deletions(-) diff --git a/components/omega/doc/devGuide/OceanState.md b/components/omega/doc/devGuide/OceanState.md index f74bf2eef9c8..ee56e218b5fe 100644 --- a/components/omega/doc/devGuide/OceanState.md +++ b/components/omega/doc/devGuide/OceanState.md @@ -21,7 +21,10 @@ 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 read from the initial-condition and restart files. After initialization, the default state object can be retrieved via: ``` @@ -71,6 +74,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/VertCoord.md b/components/omega/doc/devGuide/VertCoord.md index f3da8dbbd0b3..a2dfd72a5237 100644 --- a/components/omega/doc/devGuide/VertCoord.md +++ b/components/omega/doc/devGuide/VertCoord.md @@ -69,6 +69,27 @@ 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; after the read, `initSurfacePressure()` must be called to exchange the halo and copy 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 diff --git a/components/omega/doc/userGuide/OceanState.md b/components/omega/doc/userGuide/OceanState.md index b13404587134..b8596b92a348 100644 --- a/components/omega/doc/userGuide/OceanState.md +++ b/components/omega/doc/userGuide/OceanState.md @@ -2,7 +2,15 @@ ## 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 read from the initial-condition and restart files. diff --git a/components/omega/doc/userGuide/VertCoord.md b/components/omega/doc/userGuide/VertCoord.md index 7f9eb1b0b59d..d8695de2f0ec 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,6 +48,22 @@ 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 read from the initial-condition +file and written to restart files. +For a free-surface ocean run with standard atmospheric forcing, the surface pressure is +zero (0 Pa) in 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 From a1ebc8de8f47d0f7dd9decfb336d947a020fe70c Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 2 Jul 2026 12:23:34 -0500 Subject: [PATCH 5/8] Add optional-read support for IO stream fields Add a per-field "optional read" flag to Field so that a field may be absent from an input file. When a read stream does not contain an optional field's variable, IOStream::readFieldData now logs and returns success, leaving the attached array at the fill value set by attachData, instead of failing the whole stream read. Fields are required by default, so existing behavior is unchanged. Co-Authored-By: Claude Opus 4.8 --- components/omega/src/infra/Field.cpp | 15 +++++++++++++++ components/omega/src/infra/Field.h | 15 +++++++++++++++ components/omega/src/infra/IOStream.cpp | 13 ++++++++++++- 3 files changed, 42 insertions(+), 1 deletion(-) diff --git a/components/omega/src/infra/Field.cpp b/components/omega/src/infra/Field.cpp index d2e7095b5f79..ccf2d9d7984b 100644 --- a/components/omega/src/infra/Field.cpp +++ b/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. diff --git a/components/omega/src/infra/Field.h b/components/omega/src/infra/Field.h index 6683ec9d1bc8..52a3b9a315a9 100644 --- a/components/omega/src/infra/Field.h +++ b/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 //--------------------------------------------------------------------------- diff --git a/components/omega/src/infra/IOStream.cpp b/components/omega/src/infra/IOStream.cpp index cd34fa1d7a81..f1c8f4d05068 100644 --- a/components/omega/src/infra/IOStream.cpp +++ b/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 From 7630a995e47fef27c60ab24e62e42023c432fe3e Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 2 Jul 2026 12:27:16 -0500 Subject: [PATCH 6/8] Test optional-read handling in IOStream Add a case to IOStreamTest that registers a field absent from the input file, marks it as an optional read, and adds it to the InitialState stream. Verify the stream read succeeds and the field array retains the fill value on all owned cells, confirming a missing optional field is skipped rather than failing the read. Co-Authored-By: Claude Opus 4.8 --- components/omega/test/infra/IOStreamTest.cpp | 32 ++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp index b3126ad21a64..a8fe2fde2568 100644 --- a/components/omega/test/infra/IOStreamTest.cpp +++ b/components/omega/test/infra/IOStreamTest.cpp @@ -18,6 +18,7 @@ #include "Eos.h" #include "Error.h" #include "Field.h" +#include "FillValues.h" #include "Forcing.h" #include "Halo.h" #include "HorzMesh.h" @@ -179,11 +180,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(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(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); From 51a748930af2db708c3a1350942be228f5b0d57a Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 2 Jul 2026 12:29:09 -0500 Subject: [PATCH 7/8] Make SurfacePressure an optional read defaulting to zero Mark the VertCoord SurfacePressure field as an optional read so it is no longer required in the initial-condition or restart file. SurfacePressure stays in the State and Restart groups, so it is still written to history and restart files and read from the initial-state and restart files when present. When it is absent, initSurfacePressure detects the leftover fill value on the owned cells and defaults the field to zero before the halo exchange and host copy. Update the StateTest comment to note the field is now optional and defaults to zero. Co-Authored-By: Claude Opus 4.8 --- components/omega/src/ocn/VertCoord.cpp | 37 +++++++++++++++++++++---- components/omega/test/ocn/StateTest.cpp | 5 ++-- 2 files changed, 35 insertions(+), 7 deletions(-) diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp index 4d4f0259a78a..dfedfa1c32cc 100644 --- a/components/omega/src/ocn/VertCoord.cpp +++ b/components/omega/src/ocn/VertCoord.cpp @@ -310,6 +310,10 @@ void VertCoord::defineFields() { 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"; @@ -459,10 +463,12 @@ void VertCoord::defineFields() { PseudoThicknessTargetField->attachData(PseudoThicknessTarget); SshField->attachData(SshCell); - // Add SurfacePressure to the State and Restart groups so it is read from the - // initial-state and restart files. VertCoord initializes before OceanState, - // which also contributes fields to these groups, so create each group only - // if it does not already exist. + // 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); @@ -702,9 +708,30 @@ void VertCoord::setStreamArrays(const bool ReadStream, Halo *MeshHalo) { } //------------------------------------------------------------------------------ -// Exchange halo and copy SurfacePressure to host after the initial-state or +// 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); } diff --git a/components/omega/test/ocn/StateTest.cpp b/components/omega/test/ocn/StateTest.cpp index 94b373298c90..643dc54a211e 100644 --- a/components/omega/test/ocn/StateTest.cpp +++ b/components/omega/test/ocn/StateTest.cpp @@ -274,8 +274,9 @@ int main(int argc, char *argv[]) { LOG_INFO("State: Out-of-range {} Non-zero: {}", Count2, Count1); } - // SurfacePressure is now owned by VertCoord but read from the same - // InitialState stream; verify its host mirror is sized and read + // 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) { From 27f01a430883e795e75a71d7c94dd6a57e70b098 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 2 Jul 2026 12:32:40 -0500 Subject: [PATCH 8/8] Document optional-read fields and SurfacePressure default Describe the new per-field optional-read capability in the Field and IOStreams dev guides and in the fill-values design doc: a field marked with setOptionalRead is skipped rather than failing when absent from an input file, keeping its fill value for the owning module to default. Update the VertCoord and OceanState guides to note that SurfacePressure is an optional read that defaults to zero when not present in the initial-condition or restart file. Co-Authored-By: Claude Opus 4.8 --- components/omega/doc/design/FillValues.md | 7 +++++++ components/omega/doc/devGuide/Field.md | 12 ++++++++++++ components/omega/doc/devGuide/IOStreams.md | 14 ++++++++++++++ components/omega/doc/devGuide/OceanState.md | 4 +++- components/omega/doc/devGuide/VertCoord.md | 10 ++++++++-- components/omega/doc/userGuide/OceanState.md | 3 ++- components/omega/doc/userGuide/VertCoord.md | 7 ++++--- 7 files changed, 50 insertions(+), 7 deletions(-) diff --git a/components/omega/doc/design/FillValues.md b/components/omega/doc/design/FillValues.md index 20a6958a1ed8..0602e6e32bb1 100644 --- a/components/omega/doc/design/FillValues.md +++ b/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 diff --git a/components/omega/doc/devGuide/Field.md b/components/omega/doc/devGuide/Field.md index 62f2ad295d73..b295ff9609d7 100644 --- a/components/omega/doc/devGuide/Field.md +++ b/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++ 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`` 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 ee56e218b5fe..d6cb4db31d45 100644 --- a/components/omega/doc/devGuide/OceanState.md +++ b/components/omega/doc/devGuide/OceanState.md @@ -24,7 +24,9 @@ allocates the `NormalVelocity` and `PseudoThickness` arrays for a given number o 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 read from the initial-condition and restart files. +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: ``` diff --git a/components/omega/doc/devGuide/VertCoord.md b/components/omega/doc/devGuide/VertCoord.md index a2dfd72a5237..fd2eaccceee1 100644 --- a/components/omega/doc/devGuide/VertCoord.md +++ b/components/omega/doc/devGuide/VertCoord.md @@ -82,8 +82,14 @@ during construction, `SurfacePressure` is a prognostic quantity read from the in 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; after the read, `initSurfacePressure()` must be called to exchange the halo and copy the -device array to the host mirror: +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()); ``` diff --git a/components/omega/doc/userGuide/OceanState.md b/components/omega/doc/userGuide/OceanState.md index b8596b92a348..53bc70e6d404 100644 --- a/components/omega/doc/userGuide/OceanState.md +++ b/components/omega/doc/userGuide/OceanState.md @@ -13,4 +13,5 @@ This involves a halo update, time level index update, and updating the `IOFields 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 read from the initial-condition and restart files. +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/VertCoord.md b/components/omega/doc/userGuide/VertCoord.md index d8695de2f0ec..80fa99f0da08 100644 --- a/components/omega/doc/userGuide/VertCoord.md +++ b/components/omega/doc/userGuide/VertCoord.md @@ -57,10 +57,11 @@ standard atmosphere of 101325 Pa) at the top of the ocean column. It is the top 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 read from the initial-condition -file and written to restart files. +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) in the initial-condition file. +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.