diff --git a/core/src/include/ModelComponent.hpp b/core/src/include/ModelComponent.hpp index 73e944c94..26f8ae0a0 100644 --- a/core/src/include/ModelComponent.hpp +++ b/core/src/include/ModelComponent.hpp @@ -1,7 +1,7 @@ /*! * @file ModelComponent.hpp * - * @date 1 Jul 2024 + * @date 06 Dec 2024 * @author Tim Spain * @author Einar Ólason */ @@ -64,6 +64,8 @@ namespace Protected { inline constexpr TextTag WIND_V = "WIND_V"; // y(north)-ward component of wind, m s⁻¹ inline constexpr TextTag ICE_U = "ICE_U"; // x(east)-ward ice velocity, m s⁻¹ inline constexpr TextTag ICE_V = "ICE_V"; // y(north)-ward ice velocity, m s⁻¹ + inline constexpr TextTag IO_STRESS_X = "IO_STRESS_X"; // x(east)-ward ice-ocean stress, Pa + inline constexpr TextTag IO_STRESS_Y = "IO_STRESS_Y"; // y(north)-ward ice-ocean stress, Pa // Slab ocean fields inline constexpr TextTag SLAB_SST = "SLAB_SST"; // Slab ocean sea surface temperature, ˚C inline constexpr TextTag SLAB_SSS = "SLAB_SSS"; // Slab ocean sea surface salinity, ˚C diff --git a/core/src/include/gridNames.hpp b/core/src/include/gridNames.hpp index 82ddcde61..6fc8a6f7a 100644 --- a/core/src/include/gridNames.hpp +++ b/core/src/include/gridNames.hpp @@ -29,6 +29,9 @@ static const std::string vWindName = "vwind"; static const std::string uOceanName = "uocean"; static const std::string vOceanName = "vocean"; +static const std::string uIOStressName = "uiostress"; +static const std::string vIOStressName = "viostress"; + static const std::string coordsName = "coords"; static const std::string latitudeName = "latitude"; static const std::string longitudeName = "longitude"; diff --git a/core/src/modules/DynamicsModule/BBMDynamics.cpp b/core/src/modules/DynamicsModule/BBMDynamics.cpp index 935d8aecc..f6ada9cd0 100644 --- a/core/src/modules/DynamicsModule/BBMDynamics.cpp +++ b/core/src/modules/DynamicsModule/BBMDynamics.cpp @@ -141,6 +141,9 @@ void BBMDynamics::update(const TimestepTime& tst) uice = kernel.getDG0Data(uName); vice = kernel.getDG0Data(vName); + + taux = kernel.getDG0Data(uIOStressName); + tauy = kernel.getDG0Data(vIOStressName); } // All data for prognostic output diff --git a/core/src/modules/DynamicsModule/MEVPDynamics.cpp b/core/src/modules/DynamicsModule/MEVPDynamics.cpp index 823fc8f16..f8cdde95f 100644 --- a/core/src/modules/DynamicsModule/MEVPDynamics.cpp +++ b/core/src/modules/DynamicsModule/MEVPDynamics.cpp @@ -108,6 +108,9 @@ void MEVPDynamics::update(const TimestepTime& tst) uice = kernel.getDG0Data(uName); vice = kernel.getDG0Data(vName); + + taux = kernel.getDG0Data(uIOStressName); + tauy = kernel.getDG0Data(vIOStressName); } ModelState MEVPDynamics::getStateRecursive(const OutputSpec& os) const diff --git a/core/src/modules/include/IDynamics.hpp b/core/src/modules/include/IDynamics.hpp index 29d9423d2..82c8614f4 100644 --- a/core/src/modules/include/IDynamics.hpp +++ b/core/src/modules/include/IDynamics.hpp @@ -1,7 +1,7 @@ /*! * @file IDynamics.hpp * - * @date 7 Sep 2023 + * @date 06 Dec 2024 * @author Tim Spain */ @@ -23,6 +23,8 @@ class IDynamics : public ModelComponent { : uice(ModelArray::Type::H) , vice(ModelArray::Type::H) , damage(ModelArray::Type::H) + , taux(ModelArray::Type::H) + , tauy(ModelArray::Type::H) , hice(getStore()) , cice(getStore()) , hsnow(getStore()) @@ -34,6 +36,8 @@ class IDynamics : public ModelComponent { , m_usesDamage(usesDamageIn) { getStore().registerArray(Shared::DAMAGE, &damage, RW); + getStore().registerArray(Protected::IO_STRESS_X, &taux, RO); + getStore().registerArray(Protected::IO_STRESS_Y, &tauy, RO); } virtual ~IDynamics() = default; @@ -76,12 +80,14 @@ class IDynamics : public ModelComponent { HField vice; // Updated damage array HField damage; + // Ice-ocean stress (for the coupler, mostly) + HField taux; + HField tauy; // References to the DG0 finite volume data arrays ModelArrayRef hice; ModelArrayRef cice; ModelArrayRef hsnow; ModelArrayRef damage0; - // ModelArrayRef damage; // References to the forcing velocity arrays ModelArrayRef uwind; diff --git a/core/src/modules/include/ProtectedArrayNames.ipp b/core/src/modules/include/ProtectedArrayNames.ipp index f273e9d9f..7781eb7ad 100644 --- a/core/src/modules/include/ProtectedArrayNames.ipp +++ b/core/src/modules/include/ProtectedArrayNames.ipp @@ -40,3 +40,5 @@ { "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⁻² + { "taux", "IO_STRESS_X" }, // Ice-ocean stress x(east) direction, Pa + { "tauy", "IO_STRESS_Y" }, // Ice-ocean stress y(north) direction, Pa diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index 64d6d974b..b72c3d79c 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -98,6 +98,16 @@ ModelArray CGDynamicsKernel::getDG0Data(const std::string& name) co DGVector vtmp(*smesh); Nextsim::Interpolations::CG2DG(*smesh, vtmp, v); return DGModelArray::dg2ma(vtmp, data); + } else if (name == uIOStressName) { + ModelArray data(ModelArray::Type::U); + DGVector utmp(*smesh); + Nextsim::Interpolations::CG2DG(*smesh, utmp, getIceOceanStress(name)); + return DGModelArray::dg2ma(utmp, data); + } else if (name == vIOStressName) { + ModelArray data(ModelArray::Type::V); + DGVector vtmp(*smesh); + Nextsim::Interpolations::CG2DG(*smesh, vtmp, getIceOceanStress(name)); + return DGModelArray::dg2ma(vtmp, data); } else { return DynamicsKernel::getDG0Data(name); } diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp index 99254598a..ffe9cd353 100644 --- a/dynamics/src/include/BrittleCGDynamicsKernel.hpp +++ b/dynamics/src/include/BrittleCGDynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file BrittleCGDynamicsKernel.hpp * - * @date 19 Nov 2024 + * @date 06 Dec 2024 * @author Tim Spain * @author Einar Ólason */ @@ -163,6 +163,22 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern } } + double getIceOceanStressElement(const std::string& name, const int i) const override + { + const double FOcean = params.COcean * params.rhoOcean; + + const double uOceanRel = uOcean(i) - avgU(i); + const double vOceanRel = vOcean(i) - avgV(i); + const double cPrime = FOcean * std::hypot(uOceanRel, vOceanRel); + + if (name == uIOStressName) + return cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle); + else if (name == vIOStressName) + return cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle); + else + return std::numeric_limits::quiet_NaN(); + } + protected: CGVector avgU; CGVector avgV; diff --git a/dynamics/src/include/CGDynamicsKernel.hpp b/dynamics/src/include/CGDynamicsKernel.hpp index 7c50f1ef6..53ece021e 100644 --- a/dynamics/src/include/CGDynamicsKernel.hpp +++ b/dynamics/src/include/CGDynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file CGDynamicsKernel.hpp * - * @date Jan 31, 2024 + * @date 06 Dec 2024 * @author Tim Spain */ @@ -47,6 +47,23 @@ class CGDynamicsKernel : public DynamicsKernel { void applyBoundaries() override; void prepareAdvection() override; + virtual inline double getIceOceanStressElement(const std::string& name, const int i) const = 0; + CGVector getIceOceanStress(const std::string& name) const + { + if (name != uIOStressName && name != vIOStressName) + throw std::logic_error(std::string(__func__) + " called with an unknown argument " + + name + ". Only " + uIOStressName + " and " + vIOStressName + " are supported\n"); + + CGVector tau; + tau.resizeLike(u); + +#pragma omp parallel for + for (int i = 0; i < tau.rows(); ++i) + tau(i) = getIceOceanStressElement(name, i); + + return tau; + } + protected: void addStressTensorCell(const size_t eid, const size_t cx, const size_t cy); void dirichletZero(CGVector&) const; diff --git a/dynamics/src/include/FreeDriftDynamicsKernel.hpp b/dynamics/src/include/FreeDriftDynamicsKernel.hpp index f007ecfb7..ccacacc08 100644 --- a/dynamics/src/include/FreeDriftDynamicsKernel.hpp +++ b/dynamics/src/include/FreeDriftDynamicsKernel.hpp @@ -4,7 +4,7 @@ * Implementation of "classic free drift", where we ignore all \rho h terms in the momentum * equation. This is equivalent to assuming that the ice is very thin. * - * @date 19 Nov 2024 + * @date 06 Dec 2024 * @author Tim Spain * @author Einar Ólason */ @@ -69,6 +69,22 @@ template class FreeDriftDynamicsKernel : public CGDynamicsKern + NansenNumber * (-uAtmos(i) * sinOceanAngle + vAtmos(i) * cosOceanAngle); } } + + double getIceOceanStressElement(const std::string& name, const int i) const override + { + const double FOcean = params.COcean * params.rhoOcean; + + const double uOceanRel = uOcean(i) - u(i); + const double vOceanRel = vOcean(i) - v(i); + const double cPrime = FOcean * std::hypot(uOceanRel, vOceanRel); + + if (name == uIOStressName) + return cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle); + else if (name == vIOStressName) + return cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle); + else + return std::numeric_limits::quiet_NaN(); + } }; } /* namespace Nextsim */ diff --git a/dynamics/src/include/VPCGDynamicsKernel.hpp b/dynamics/src/include/VPCGDynamicsKernel.hpp index badc3c33c..8bb7149a4 100644 --- a/dynamics/src/include/VPCGDynamicsKernel.hpp +++ b/dynamics/src/include/VPCGDynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file VPCGDynamicsKernel.hpp * - * @date 19 Nov 2024 + * @date 06 Dec 2024 * @author Tim Spain */ @@ -89,6 +89,22 @@ template class VPCGDynamicsKernel : public CGDynamicsKernel::update(tst); } + double getIceOceanStressElement(const std::string& name, const int i) const override + { + const double FOcean = params.COcean * params.rhoOcean; + + double uOcnRel = u(i) - uOcean(i); + double vOcnRel = v(i) - vOcean(i); + double absocn = sqrt(SQR(uOcnRel) + SQR(vOcnRel)); + + if (name == uIOStressName) + return FOcean * absocn * uOcnRel; + else if (name == vIOStressName) + return FOcean * absocn * vOcnRel; + else + return std::numeric_limits::quiet_NaN(); + } + protected: StressUpdateStep& stressStep; const VPParameters& params;