Add sea-surface height (#755)

# Add sea-surface height
Fixes #538 

### Task List
- [x] Defined the tests that specify a complete and functioning change
(*It may help to create a [design specification & test
specification](../../../wiki/Specification-Template)*)
- [x] Implemented the source code change that satisfies the tests
- [x] Documented the feature by providing worked example
- [ ] Updated the README or other documentation
- [x] Completed the pre-Request checklist below

---
# Change Description

This PR adds the contribution of the sea-surface-height gradient to both
the mEVP and BBM codes (through VPCGDynamicsKernel and
BrittleDynamicsKernel). SSH is read in from IOceanBoundary and passed to
the dynamical core. It can also be output through ConfigOutput as "ssh".

This PR is a repeat of PR #662 with some fixes to make sure the
integration test runs.

---
# Test Description

No unit test was implemented for this feature, and no quantitative test
exists. Qualitative testing shows that setting a fixed
sea-surface-height gradient gives drift speed close to
back-of-the-envelope estimates (13 cm/s modelled, 20 cm/s from
back-of-the-envelope).

---
# Documentation Impact

None

---
# Other Details

None

---
### Pre-Request Checklist

- [x] The requirements of this pull request are fully captured in an
issue or design specification and are linked and summarised in the
description of this PR
- [x] No new warnings are generated
- [x] The documentation has been updated (or an issue has been created
to track the corresponding change)
- [x] Methods and Tests are commented such that they can be understood
without having to obtain additional context
- [x] This PR/Issue is labelled as a bug/feature/enhancement/breaking
change
- [x] File dates have been updated to reflect modification date
- [x] This change conforms to the conventions described in the README

core/src/ParaGridIO.cpp

@@ -247,7 +247,12 @@ ModelState ParaGridIO::readForcingTimeStatic(
             extentArray.push_back(ModelArray::definedDimensions.at(*riter).localLength);
         }
 
+        auto availableForcings = dataGroup.getVars();
         for (const std::string& varName : forcings) {
+            // Don't try to read non-existent data
+            if (!availableForcings.count(varName)) {
+                continue;
+            }
             netCDF::NcVar var = dataGroup.getVar(varName);
             state.data[varName] = ModelArray(ModelArray::Type::H);
             ModelArray& data = state.data.at(varName);

core/src/include/ModelComponent.hpp

@@ -50,6 +50,7 @@ namespace Protected {
         = "EXT_SSS"; // sea surface salinity from coupling or forcing, PSU
     inline constexpr TextTag EXT_SST
         = "EXT_SST"; // sea surface temperature from coupling or forcing, ˚C
+    inline constexpr TextTag SSH = "SSH"; // sea surface height, m
     inline constexpr TextTag EVAP_MINUS_PRECIP
         = "E-P"; // E-P atmospheric freshwater flux, kg s⁻¹ m⁻²
     // Derived fields, calculated once per timestep

core/src/include/constants.hpp

@@ -1,6 +1,6 @@
 /*!
  * @file constants.hpp
- * @date Sep 14, 2021
+ * @date 05 Dec 2024
  * @author Tim Spain
  */
 
@@ -31,6 +31,12 @@ const double Tt = 273.16;
 
 //! Ratio of circumference to radius
 const double tau = 6.28318530717958647652;
+
+// degrees to radians as a hex float
+const double deg2rad = 0x1.1df46a2529d39p-6;
+
+// radians to degrees
+const double rad2deg = 0x1.ca5dc1a63c1f8p+5;
 }
 
 //! Properties of water ice around 0˚C and 101.3 kPa
@@ -131,10 +137,10 @@ inline double kelvin(double celsius) { return celsius + Water::Tf; }
 inline double celsius(double kelvin) { return kelvin - Water::Tf; }
 
 //! Convert an angle from radians to degrees
-inline double degrees(double radians) { return radians * 360 / PhysicalConstants::tau; }
+inline double degrees(double radians) { return radians * PhysicalConstants::rad2deg; }
 
 //! Convert an angle from degrees to radians
-inline double radians(double degrees) { return degrees * PhysicalConstants::tau / 360; }
+inline double radians(double degrees) { return degrees * PhysicalConstants::deg2rad; }
 
 //! Convert a pressure from Pa to mbar
 inline double mbar(double pascals) { return pascals / 100; }

core/src/include/gridNames.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file gridNames.hpp
  *
- * @date Oct 24, 2022
+ * @date Aug 23, 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -28,6 +28,9 @@ static const std::string uWindName = "uwind";
 static const std::string vWindName = "vwind";
 static const std::string uOceanName = "uocean";
 static const std::string vOceanName = "vocean";
+static const std::string sshName = "ssh";
+// Mixed layer depth
+static const std::string mldName = "mld";
 
 static const std::string uIOStressName = "uiostress";
 static const std::string vIOStressName = "viostress";

core/src/modules/DynamicsModule/BBMDynamics.cpp

@@ -1,13 +1,13 @@
 /*!
  * @file BBMDynamics.cpp
  *
- * @date 20 Nov 2024
+ * @date 05 Dec 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
 
 #include "include/BBMDynamics.hpp"
-
+#include "include/constants.hpp"
 #include "include/gridNames.hpp"
 
 namespace Nextsim {
@@ -79,7 +79,7 @@ void BBMDynamics::setData(const ModelState::DataMap& ms)
 
     ModelArray coords = ms.at(coordsName);
     if (isSpherical) {
-        coords *= radians;
+        coords *= PhysicalConstants::deg2rad;
     }
     // TODO: Some encoding of the periodic edge boundary conditions
     kernel.initialise(coords, isSpherical, ms.at(maskName));
@@ -127,6 +127,7 @@ void BBMDynamics::update(const TimestepTime& tst)
     kernel.setData(vWindName, vwind.data());
     kernel.setData(uOceanName, uocean.data());
     kernel.setData(vOceanName, vocean.data());
+    kernel.setData(sshName, ssh.data());
 
     /*
      * Ice velocity components are stored in the dynamics, and not changed by the model outside the

core/src/modules/DynamicsModule/MEVPDynamics.cpp

@@ -1,24 +1,21 @@
 /*!
  * @file MEVPDynamics.cpp
  *
- * @date 20 Nov 2024
+ * @date 05 Dec 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Piotr Minakowski <piotr.minakowski@ovgu.de>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
 
 #include "include/MEVPDynamics.hpp"
-
+#include "include/constants.hpp"
 #include "include/gridNames.hpp"
 
 #include <string>
 #include <vector>
 
 namespace Nextsim {
 
-// Degrees to radians as a hex float
-static const double radians = 0x1.1df46a2529d39p-6;
-
 // TODO: We should use getName() here, but it isn't static.
 static const std::string prefix = "MEVPDynamics"; // MEVPDynamics::getName();
 static const std::map<int, std::string> keyMap = {
@@ -70,7 +67,7 @@ void MEVPDynamics::setData(const ModelState::DataMap& ms)
 
     ModelArray coords = ms.at(coordsName);
     if (isSpherical) {
-        coords *= radians;
+        coords *= PhysicalConstants::deg2rad;
     }
     // TODO: Some encoding of the periodic edge boundary conditions
     kernel.initialise(coords, isSpherical, ms.at(maskName));
@@ -97,6 +94,7 @@ void MEVPDynamics::update(const TimestepTime& tst)
     kernel.setData(vWindName, vwind.data());
     kernel.setData(uOceanName, uocean.data());
     kernel.setData(vOceanName, vocean.data());
+    kernel.setData(sshName, ssh.data());
 
     // kernel.setData(uName, uice);
     // kernel.setData(vName, vice);

core/src/modules/include/IDynamics.hpp

@@ -33,6 +33,7 @@ class IDynamics : public ModelComponent {
         , vwind(getStore())
         , uocean(getStore())
         , vocean(getStore())
+        , ssh(getStore())
         , m_usesDamage(usesDamageIn)
     {
         getStore().registerArray(Shared::DAMAGE, &damage, RW);
@@ -94,6 +95,7 @@ class IDynamics : public ModelComponent {
     ModelArrayRef<Protected::WIND_V> vwind;
     ModelArrayRef<Protected::OCEAN_U> uocean;
     ModelArrayRef<Protected::OCEAN_V> vocean;
+    ModelArrayRef<Protected::SSH> ssh;
 
     // Does this implementation of the dynamics use damage?
     bool m_usesDamage;

core/src/modules/include/ProtectedArrayNames.ipp

@@ -1,7 +1,7 @@
 /*!
  * @file ProtectedArrayNames.ipp
  *
- * @date 1 Jul 2024
+ * @date 23 Aug 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
@@ -40,5 +40,6 @@
     { "sss_slab", "SLAB_SSS" }, // Slab ocean surface salinity PSU
     { "qdw", "SLAB_QDW" }, // Slab ocean temperature nudging heat flux, W m⁻²
     { "fdw", "SLAB_FDW" }, // Slab ocean salinity nudging water flux, kg s⁻¹ m⁻²
+    { "ssh", "SSH" }, // Slab ocean salinity nudging water flux, kg s⁻¹ m⁻²
     { "taux", "IO_STRESS_X" }, // Ice-ocean stress x(east) direction, Pa
     { "tauy", "IO_STRESS_Y" }, // Ice-ocean stress y(north) direction, Pa

dynamics/src/CGDynamicsKernel.cpp

@@ -1,7 +1,7 @@
 /*!
  * @file CGDynamicsKernel.cpp
  *
- * @date Jan 26, 2024
+ * @date 06 Dec 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -36,6 +36,9 @@ void CGDynamicsKernel<DGadvection>::initialise(
     cgH.resize_by_mesh(*smesh);
     cgA.resize_by_mesh(*smesh);
 
+    xGradSeaSurfaceHeight.resize_by_mesh(*smesh);
+    yGradSeaSurfaceHeight.resize_by_mesh(*smesh);
+
     dStressX.resize_by_mesh(*smesh);
     dStressY.resize_by_mesh(*smesh);
 
@@ -118,6 +121,141 @@ template <int DGadvection> void CGDynamicsKernel<DGadvection>::prepareAdvection(
     dgtransport->prepareAdvection(u, v);
 }
 
+template <int DGadvection>
+void CGDynamicsKernel<DGadvection>::ComputeGradientOfSeaSurfaceHeight(
+    const DGVector<1>& seaSurfaceHeight)
+{
+    // First transform to CG1 Vector and do all computations in CG1
+    CGVector<1> cgSeasurfaceHeight(*smesh);
+    Interpolations::DG2CG(*smesh, cgSeasurfaceHeight, seaSurfaceHeight);
+
+    CGVector<1> uGrad(*smesh);
+    CGVector<1> vGrad(*smesh);
+    uGrad.setZero();
+    vGrad.setZero();
+
+    // parallelization in stripes
+    for (size_t p = 0; p < 2; ++p) {
+#pragma omp parallel for schedule(static)
+        for (size_t cy = 0; cy < smesh->ny; ++cy) {
+            //!< loop over all cells of the mesh
+            if (cy % 2 == p) {
+                size_t eid = smesh->nx * cy;
+                size_t cg1id = cy * (smesh->nx + 1);
+                for (size_t cx = 0; cx < smesh->nx; ++cx, ++eid, ++cg1id) {
+                    // get local CG nodes
+                    Eigen::Vector<Nextsim::FloatType, 4> loc_cgSSH = { cgSeasurfaceHeight(cg1id),
+                        cgSeasurfaceHeight(cg1id + 1), cgSeasurfaceHeight(cg1id + smesh->nx + 1),
+                        cgSeasurfaceHeight(cg1id + smesh->nx + 1 + 1) };
+
+                    // compute grad
+                    Eigen::Vector<Nextsim::FloatType, 4> tx = pmap->dX_SSH[eid] * loc_cgSSH;
+                    Eigen::Vector<Nextsim::FloatType, 4> ty = pmap->dY_SSH[eid] * loc_cgSSH;
+
+                    // add global vector
+                    uGrad(cg1id) -= tx(0);
+                    uGrad(cg1id + 1) -= tx(1);
+                    uGrad(cg1id + smesh->nx + 1) -= tx(2);
+                    uGrad(cg1id + smesh->nx + 1 + 1) -= tx(3);
+                    vGrad(cg1id) -= ty(0);
+                    vGrad(cg1id + 1) -= ty(1);
+                    vGrad(cg1id + smesh->nx + 1) -= ty(2);
+                    vGrad(cg1id + smesh->nx + 1 + 1) -= ty(3);
+                }
+            }
+        }
+    }
+
+    // scale with mass
+#pragma omp parallel for
+    for (size_t i = 0; i < uGrad.rows(); ++i) {
+        uGrad(i) /= pmap->lumpedcg1mass(i);
+        vGrad(i) /= pmap->lumpedcg1mass(i);
+    }
+
+    ///// correct boundary (just extend in last elements)
+    size_t cg1row = smesh->nx + 1;
+    size_t topleft = smesh->ny * cg1row;
+    for (size_t i = 1; i < smesh->nx; ++i) // bottom / top
+    {
+        uGrad(i) = uGrad(i + cg1row);
+        vGrad(i) = vGrad(i + cg1row);
+        uGrad(topleft + i) = uGrad(topleft + i - cg1row);
+        vGrad(topleft + i) = vGrad(topleft + i - cg1row);
+    }
+    for (size_t i = 1; i < smesh->ny; ++i) // left / right
+    {
+        uGrad(i * cg1row) = uGrad(i * cg1row + 1);
+        vGrad(i * cg1row) = vGrad(i * cg1row + 1);
+
+        uGrad(i * cg1row + cg1row - 1) = uGrad(i * cg1row + cg1row - 1 - 1);
+        vGrad(i * cg1row + cg1row - 1) = vGrad(i * cg1row + cg1row - 1 - 1);
+    }
+    // corners
+    uGrad(0) = uGrad(cg1row + 1);
+    vGrad(0) = vGrad(cg1row + 1);
+    uGrad(smesh->nx) = uGrad(smesh->nx + cg1row - 1);
+    vGrad(smesh->nx) = vGrad(smesh->nx + cg1row - 1);
+    uGrad(smesh->ny * cg1row) = uGrad((smesh->ny - 1) * cg1row + 1);
+    vGrad(smesh->ny * cg1row) = vGrad((smesh->ny - 1) * cg1row + 1);
+    uGrad((smesh->ny + 1) * cg1row - 1) = uGrad((smesh->ny) * cg1row - 1 - 1);
+    vGrad((smesh->ny + 1) * cg1row - 1) = vGrad((smesh->ny) * cg1row - 1 - 1);
+
+    // Interpolate to CG2 (maybe own function in interpolation?)
+    if (CGDEGREE == 1) {
+        xGradSeaSurfaceHeight = uGrad;
+        yGradSeaSurfaceHeight = vGrad;
+    } else {
+        // outer nodes
+        size_t icg1 = 0;
+#pragma omp parallel for
+        for (size_t iy = 0; iy <= smesh->ny; ++iy) {
+            size_t icg2 = (2 * smesh->nx + 1) * 2 * iy;
+            for (size_t ix = 0; ix <= smesh->nx; ++ix, ++icg1, icg2 += 2) {
+                xGradSeaSurfaceHeight(icg2) = uGrad(icg1);
+                yGradSeaSurfaceHeight(icg2) = vGrad(icg1);
+            }
+        }
+        // along lines
+#pragma omp parallel for
+        for (size_t iy = 0; iy <= smesh->ny; ++iy) // horizontal
+        {
+            size_t icg1 = (smesh->nx + 1) * iy;
+            size_t icg2 = (2 * smesh->nx + 1) * 2 * iy + 1;
+            for (size_t ix = 0; ix < smesh->nx; ++ix, ++icg1, icg2 += 2) {
+                xGradSeaSurfaceHeight(icg2) = 0.5 * (uGrad(icg1) + uGrad(icg1 + 1));
+                yGradSeaSurfaceHeight(icg2) = 0.5 * (vGrad(icg1) + vGrad(icg1 + 1));
+            }
+        }
+#pragma omp parallel for
+        for (size_t iy = 0; iy < smesh->ny; ++iy) // vertical
+        {
+            size_t icg1 = (smesh->nx + 1) * iy;
+            size_t icg2 = (2 * smesh->nx + 1) * (2 * iy + 1);
+            for (size_t ix = 0; ix <= smesh->nx; ++ix, ++icg1, icg2 += 2) {
+                xGradSeaSurfaceHeight(icg2) = 0.5 * (uGrad(icg1) + uGrad(icg1 + cg1row));
+                yGradSeaSurfaceHeight(icg2) = 0.5 * (vGrad(icg1) + vGrad(icg1 + cg1row));
+            }
+        }
+
+        // midpoints
+#pragma omp parallel for
+        for (size_t iy = 0; iy < smesh->ny; ++iy) // vertical
+        {
+            size_t icg1 = (smesh->nx + 1) * iy;
+            size_t icg2 = (2 * smesh->nx + 1) * (2 * iy + 1) + 1;
+            for (size_t ix = 0; ix < smesh->nx; ++ix, ++icg1, icg2 += 2) {
+                xGradSeaSurfaceHeight(icg2) = 0.25
+                    * (uGrad(icg1) + uGrad(icg1 + 1) + uGrad(icg1 + cg1row)
+                        + uGrad(icg1 + cg1row + 1));
+                yGradSeaSurfaceHeight(icg2) = 0.25
+                    * (vGrad(icg1) + vGrad(icg1 + 1) + vGrad(icg1 + cg1row)
+                        + vGrad(icg1 + cg1row + 1));
+            }
+        }
+    }
+}
+
 template <int DGadvection> void CGDynamicsKernel<DGadvection>::prepareIteration(const DataMap& data)
 {
     // interpolate ice height and concentration to local cg Variables
@@ -126,6 +264,10 @@ template <int DGadvection> void CGDynamicsKernel<DGadvection>::prepareIteration(
     Interpolations::DG2CG(*smesh, cgA, data.at(ciceName));
     VectorManipulations::CGAveragePeriodic(*smesh, cgA);
 
+    // Reinit the gradient of the sea surface height. Not done by
+    // DataMap as seaSurfaceHeight is always dG(0)
+    ComputeGradientOfSeaSurfaceHeight(DynamicsKernel<DGadvection, DGstressComp>::seaSurfaceHeight);
+
     // limit A to [0,1] and H to [0, ...)
     cgA = cgA.cwiseMin(1.0);
     cgA = cgA.cwiseMax(1.e-4);