Add sea-surface height

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 <[email protected]>
  */
 
@@ -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 <[email protected]>
  * @author Einar Ólason <[email protected]>
  */
 
 #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 <[email protected]>
  * @author Piotr Minakowski <[email protected]>
  * @author Einar Ólason <[email protected]>
  */
 
 #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 <[email protected]>
  * @author Einar Ólason <[email protected]>
  */
@@ -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 <[email protected]>
  */
 
@@ -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);