Skip to content

Fixes Omega ssh calculation #335

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
vanroekel wants to merge 3 commits into
base: develop
Choose a base branch
from
Open

Fixes Omega ssh calculation #335

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a10434a
Fixes ssh calculation
vanroekel Feb 5, 2026
0141956
Fixes kmin and linting
vanroekel Feb 5, 2026
52001a1
Linting fixes
vanroekel Feb 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 0 additions & 24 deletions components/omega/src/ocn/AuxiliaryState.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -189,30 +189,6 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel,
});
Pacer::stop("AuxState:cellAuxState2", 2);

Pacer::start("AuxState:cellAuxState3", 2);
Copy link
Copy Markdown

@brian-oneill brian-oneill Mar 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is fine, you should still remove the LayerThicknessAux.computeVarsOnCells call from computeMomAux. I added it to AuxiliaryState::computeAll in #328 since it is now needed for tracer tendencies, not velocity tendencies.

parallelForOuter(
"cellAuxState3", {Mesh->NCellsAll},
KOKKOS_LAMBDA(int ICell, const TeamMember &Team) {
const int KMin = MinLayerCell(ICell);
const int KMax = MaxLayerCell(ICell);
const int KRange = vertRangeChunked(KMin, KMax);

parallelForInner(
Team, KRange, INNER_LAMBDA(int KChunk) {
LocLayerThicknessAux.computeVarsOnCells(
ICell, KChunk, LayerThickCell, NormalVelEdge,
TimeStepSeconds);
});
});
Pacer::stop("AuxState:cellAuxState3", 2);

Pacer::start("AuxState:computeVerticalVelocity", 2);

const auto &FluxLayerThickEdge = LayerThicknessAux.FluxLayerThickEdge;
VAdv->computeVerticalVelocity(NormalVelEdge, FluxLayerThickEdge);

Pacer::stop("AuxState:computeVerticalVelocity", 2);

Pacer::stop("AuxState:computeMomAux", 1);
}

Expand Down
3 changes: 2 additions & 1 deletion components/omega/src/ocn/Tendencies.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -375,6 +375,7 @@ void Tendencies::computeVelocityTendenciesOnly(
OMEGA_SCOPE(LocBottomDrag, BottomDrag);
OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot);
OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop);
OMEGA_SCOPE(LocSshCell, VCoord->SshCell);

Pacer::start("Tend:computeVelocityTendenciesOnly", 1);

Expand Down Expand Up @@ -433,7 +434,7 @@ void Tendencies::computeVelocityTendenciesOnly(
}

// Compute sea surface height gradient
const Array2DReal &SSHCell = AuxState->LayerThicknessAux.SshCell;
const Array1DReal &SSHCell = LocSshCell;
if (LocSSHGrad.Enabled) {
Pacer::start("Tend:SSHGrad", 2);
parallelForOuter(
Expand Down
5 changes: 2 additions & 3 deletions components/omega/src/ocn/TendencyTerms.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -173,7 +173,7 @@ class SSHGradOnEdge {
/// The functor takes edge index, vertical chunk index, and array of
/// layer thickness/SSH, outputs tendency array
KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk,
const Array2DReal &SshCell) const {
const Array1DReal &SshCell) const {

const I4 KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge));
const I4 KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge));
Expand All @@ -184,8 +184,7 @@ class SSHGradOnEdge {
for (int KVec = 0; KVec < KLen; ++KVec) {
const I4 K = KStart + KVec;
Tend(IEdge, K) -= EdgeMask(IEdge, K) * Gravity *
(SshCell(ICell1, K) - SshCell(ICell0, K)) *
InvDcEdge;
(SshCell(ICell1) - SshCell(ICell0)) * InvDcEdge;
}
}

Expand Down
58 changes: 44 additions & 14 deletions components/omega/src/ocn/VertCoord.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,20 +160,24 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord
BottomDepth = Array1DReal("BottomDepth", NCellsSize);
PressureInterface =
Array2DReal("PressureInterface", NCellsSize, NVertLayersP1);
PressureMid = Array2DReal("PressureMid", NCellsSize, NVertLayers);
ZInterface = Array2DReal("ZInterface", NCellsSize, NVertLayersP1);
ZMid = Array2DReal("ZMid", NCellsSize, NVertLayers);
PressureMid = Array2DReal("PressureMid", NCellsSize, NVertLayers);
ZInterface = Array2DReal("ZInterface", NCellsSize, NVertLayersP1);
ZMid = Array2DReal("ZMid", NCellsSize, NVertLayers);
SshCell = Array1DReal("SshCell", NCellsSize);

GeopotentialMid = Array2DReal("GeopotentialMid", NCellsSize, NVertLayers);
LayerThicknessTarget =
Array2DReal("LayerThicknessTarget", NCellsSize, NVertLayers);
RefLayerThickness =
Array2DReal("RefLayerThickness", NCellsSize, NVertLayers);

// Make host copies for device arrays not being read from file
PressureInterfaceH = createHostMirrorCopy(PressureInterface);
PressureMidH = createHostMirrorCopy(PressureMid);
ZInterfaceH = createHostMirrorCopy(ZInterface);
ZMidH = createHostMirrorCopy(ZMid);
PressureInterfaceH = createHostMirrorCopy(PressureInterface);
PressureMidH = createHostMirrorCopy(PressureMid);
ZInterfaceH = createHostMirrorCopy(ZInterface);
ZMidH = createHostMirrorCopy(ZMid);
SshCellH = createHostMirrorCopy(SshCellH);

GeopotentialMidH = createHostMirrorCopy(GeopotentialMid);
LayerThicknessTargetH = createHostMirrorCopy(LayerThicknessTarget);
RefLayerThicknessH = createHostMirrorCopy(RefLayerThickness);
Expand Down Expand Up @@ -229,13 +233,15 @@ VertCoord *VertCoord::create(
void VertCoord::defineFields() {

// Set field names (append Name if not default)
MinLayerCellFldName = "MinLayerCell";
MaxLayerCellFldName = "MaxLayerCell";
BottomDepthFldName = "BottomDepth";
PressInterfFldName = "PressureInterface";
PressMidFldName = "PressureMid";
ZInterfFldName = "ZInterface";
ZMidFldName = "ZMid";
MinLayerCellFldName = "MinLayerCell";
MaxLayerCellFldName = "MaxLayerCell";
BottomDepthFldName = "BottomDepth";
PressInterfFldName = "PressureInterface";
PressMidFldName = "PressureMid";
ZInterfFldName = "ZInterface";
ZMidFldName = "ZMid";
SshFldName = "SshCell";

GeopotFldName = "GeopotentialMid";
LyrThickTargetFldName = "LayerThicknessTarget";

Expand All @@ -249,6 +255,7 @@ void VertCoord::defineFields() {
ZMidFldName.append(Name);
GeopotFldName.append(Name);
LyrThickTargetFldName.append(Name);
SshFldName.append(Name);
}

// Create fields for VertCoord variables
Expand Down Expand Up @@ -295,6 +302,18 @@ void VertCoord::defineFields() {
DimNames // dimension names
);

auto SshField = Field::create(
SshFldName, // field name
"sea surface height at cell center", // long Name or description
"m", // units
"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
);

NDims = 2;
DimNames.resize(NDims);
DimNames[1] = "NVertLayersP1";
Expand Down Expand Up @@ -402,6 +421,7 @@ void VertCoord::defineFields() {
VCoordGroup->addField(ZMidFldName);
VCoordGroup->addField(GeopotFldName);
VCoordGroup->addField(LyrThickTargetFldName);
VCoordGroup->addField(SshFldName);

// Associate Field with data
PressureInterfaceField->attachData<Array2DReal>(PressureInterface);
Expand All @@ -410,6 +430,7 @@ void VertCoord::defineFields() {
ZMidField->attachData<Array2DReal>(ZMid);
GeopotentialMidField->attachData<Array2DReal>(GeopotentialMid);
LayerThicknessTargetField->attachData<Array2DReal>(LayerThicknessTarget);
SshField->attachData<Array1DReal>(SshCell);

} // end defineFields

Expand All @@ -431,6 +452,7 @@ VertCoord::~VertCoord() {
Field::destroy(ZMidFldName);
Field::destroy(GeopotFldName);
Field::destroy(LyrThickTargetFldName);
Field::destroy(SshFldName);
FieldGroup::destroy(GroupName);
}

Expand Down Expand Up @@ -872,6 +894,7 @@ void VertCoord::computeZHeight(
OMEGA_SCOPE(LocZInterf, ZInterface);
OMEGA_SCOPE(LocZMid, ZMid);
OMEGA_SCOPE(LocBotDepth, BottomDepth);
OMEGA_SCOPE(LocSshCell, SshCell);

parallelForOuter(
"computeZHeight", {NCellsAll},
Expand All @@ -891,6 +914,9 @@ void VertCoord::computeZHeight(
LocZInterf(ICell, KLyr) = -LocBotDepth(ICell) + Accum;
LocZMid(ICell, KLyr) =
-LocBotDepth(ICell) + Accum - 0.5 * DZ;
if (KLyr == KMin) {
LocSshCell(ICell) = LocZInterf(ICell, KLyr);
}
}
});
});
Expand Down Expand Up @@ -995,6 +1021,8 @@ void VertCoord::copyToHost() {
deepCopy(PressureMidH, PressureMid);
deepCopy(ZInterfaceH, ZInterface);
deepCopy(ZMidH, ZMid);
deepCopy(SshCellH, SshCell);

deepCopy(GeopotentialMidH, GeopotentialMid);
deepCopy(LayerThicknessTargetH, LayerThicknessTarget);
deepCopy(RefLayerThicknessH, RefLayerThickness);
Expand All @@ -1008,6 +1036,8 @@ void VertCoord::copyToDevice() {
deepCopy(PressureMid, PressureMidH);
deepCopy(ZInterface, ZInterfaceH);
deepCopy(ZMid, ZMidH);
deepCopy(SshCell, SshCellH);

deepCopy(GeopotentialMid, GeopotentialMidH);
deepCopy(LayerThicknessTarget, LayerThicknessTargetH);
deepCopy(RefLayerThickness, RefLayerThicknessH);
Expand Down
3 changes: 3 additions & 0 deletions components/omega/src/ocn/VertCoord.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,13 +107,15 @@ class VertCoord {
Array2DReal ZMid;
Array2DReal GeopotentialMid;
Array2DReal LayerThicknessTarget;
Array1DReal SshCell;

HostArray2DReal PressureInterfaceH;
HostArray2DReal PressureMidH;
HostArray2DReal ZInterfaceH;
HostArray2DReal ZMidH;
HostArray2DReal GeopotentialMidH;
HostArray2DReal LayerThicknessTargetH;
HostArray1DReal SshCellH;

// Vertical loop bounds
Array1DI4 MinLayerCell;
Expand Down Expand Up @@ -181,6 +183,7 @@ class VertCoord {
std::string ZMidFldName; ///< Field name for midpoint Z height
std::string GeopotFldName; ///< Field name for geopotential
std::string LyrThickTargetFldName; ///< Field name for target thickness
std::string SshFldName; ///< Field name for sea surface height

// methods

Expand Down
20 changes: 0 additions & 20 deletions components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,11 @@ LayerThicknessAuxVars::LayerThicknessAuxVars(const std::string &AuxStateSuffix,
Mesh->NEdgesSize, VCoord->NVertLayers),
MeanLayerThickEdge("MeanLayerThickEdge" + AuxStateSuffix,
Mesh->NEdgesSize, VCoord->NVertLayers),
SshCell("SshCell" + AuxStateSuffix, Mesh->NCellsSize,
VCoord->NVertLayers),
ProvThickness("ProvThickness" + AuxStateSuffix, Mesh->NCellsSize,
VCoord->NVertLayers),
AreaCell(Mesh->AreaCell), DvEdge(Mesh->DvEdge),
NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
EdgeSignOnCell(Mesh->EdgeSignOnCell), CellsOnEdge(Mesh->CellsOnEdge),
BottomDepth(VCoord->BottomDepth),
MinLayerEdgeBot(VCoord->MinLayerEdgeBot),
MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop),
MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {}
Expand Down Expand Up @@ -69,20 +66,6 @@ void LayerThicknessAuxVars::registerFields(const std::string &AuxGroupName,
DimNames // dimension names
);

// Sea surface height
DimNames[0] = "NCells" + DimSuffix;
auto SshCellField = Field::create(
SshCell.label(), // field name
"sea surface height at cell center", // long Name or description
"m", // units
"sea_surface_height", // 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
);

// Provisional Thickness
auto ProvThicknessField = Field::create(
ProvThickness.label(), // field name
Expand All @@ -99,20 +82,17 @@ void LayerThicknessAuxVars::registerFields(const std::string &AuxGroupName,
// Add fields to Aux field group
FieldGroup::addFieldToGroup(FluxLayerThickEdge.label(), AuxGroupName);
FieldGroup::addFieldToGroup(MeanLayerThickEdge.label(), AuxGroupName);
FieldGroup::addFieldToGroup(SshCell.label(), AuxGroupName);
FieldGroup::addFieldToGroup(ProvThickness.label(), AuxGroupName);

// Attach field data
FluxLayerThickEdgeField->attachData<Array2DReal>(FluxLayerThickEdge);
MeanLayerThickEdgeField->attachData<Array2DReal>(MeanLayerThickEdge);
SshCellField->attachData<Array2DReal>(SshCell);
ProvThicknessField->attachData<Array2DReal>(ProvThickness);
}

void LayerThicknessAuxVars::unregisterFields() const {
Field::destroy(FluxLayerThickEdge.label());
Field::destroy(MeanLayerThickEdge.label());
Field::destroy(SshCell.label());
Field::destroy(ProvThickness.label());
}

Expand Down
7 changes: 0 additions & 7 deletions components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@ class LayerThicknessAuxVars {
public:
Array2DReal FluxLayerThickEdge;
Array2DReal MeanLayerThickEdge;
Array2DReal SshCell;
Array2DReal ProvThickness;

FluxThickEdgeOption FluxThickEdgeChoice;
Expand Down Expand Up @@ -73,11 +72,6 @@ class LayerThicknessAuxVars {
const int KStart = chunkStart(KChunk, MinLayerCell(ICell));
const int KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell));

for (int KVec = 0; KVec < KLen; ++KVec) {
const int K = KStart + KVec;
SshCell(ICell, K) = LayerThickCell(ICell, K) - BottomDepth(ICell);
}

Real TmpProv[VecLength] = {0.};

Real DtInvAreaCell = Dt / AreaCell(ICell);
Expand Down Expand Up @@ -109,7 +103,6 @@ class LayerThicknessAuxVars {
Array2DI4 EdgesOnCell;
Array2DReal EdgeSignOnCell;
Array2DI4 CellsOnEdge;
Array1DReal BottomDepth;
Array1DI4 MinLayerEdgeBot;
Array1DI4 MaxLayerEdgeTop;
Array1DI4 MinLayerCell;
Expand Down
3 changes: 1 addition & 2 deletions components/omega/test/ocn/TendencyTermsTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -533,8 +533,7 @@ int testSSHGrad(int NVertLayers, Real RTol) {
},
ExactSSHGrad, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No);

// Set input array
Array2DReal SSHCell("SSHCell", Mesh->NCellsSize, NVertLayers);
Array1DReal SSHCell("SSHCell", Mesh->NCellsSize);

Err += setScalar(
KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalar(X, Y); }, SSHCell,
Expand Down
Loading