Define fluxes needed for the ocean coupling

In particular, these are: the non-solar heat flux (Q_NO_SUN), the
short-wave flux (Q_SW_OCN), the fresh-water flux (FWFLUX), and the salt
flux (SFLUX), all at the ocean surface. I also defined the short wave
flux through the base of the ice, but it's zero for now.

This is an initial commit. It compiles, but is untested and still a bit
messy.

core/src/include/ModelComponent.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file ModelComponent.hpp
  *
- * @date 1 Jul 2024
+ * @date 30 Aug 2024
  * @author Tim Spain <[email protected]>
  * @author Einar Ólason <[email protected]>
  */
@@ -71,6 +71,8 @@ namespace Protected {
         = "SLAB_QDW"; // Slab ocean temperature nudging heat flux, W m⁻²
     inline constexpr TextTag SLAB_FDW
         = "SLAB_FDW"; // Slab ocean salinity nudging water flux, kg s⁻¹ m⁻²
+    inline constexpr TextTag FWFLUX = "FWFLUX"; // Fresh-water flux at the ocean surface, kg m⁻²
+    inline constexpr TextTag SFLUX = "SFLUX"; // Salt flux at the ocean surface, kg m⁻²
 }
 
 namespace Shared {
@@ -88,6 +90,9 @@ namespace Shared {
     inline constexpr TextTag DQIA_DT
         = "DQIA_DT"; // Derivative of Qᵢₐ w.r.t. ice surface temperature  W m⁻² K⁻¹
     inline constexpr TextTag Q_PEN_SW = "Q_PEN_SW"; // Penetrating shortwave flux W m⁻²
+    inline constexpr TextTag Q_SW_OW = "Q_SW_OW"; // Net shortwave flux at ocean surface W m⁻²
+    inline constexpr TextTag Q_SW_BASE
+        = "Q_SW_BASE"; // Net shortwave flux at the base of the ice W m⁻²
     // Mass fluxes
     inline constexpr TextTag HSNOW_MELT = "HSNOW_MELT"; // Thickness of snow that melted, m
     // Atmospheric conditions

physics/src/SlabOcean.cpp

@@ -74,33 +74,29 @@ std::unordered_set<std::string> SlabOcean::hFields() const { return { sstSlabNam
 void SlabOcean::update(const TimestepTime& tst)
 {
     double dt = tst.step.seconds();
+
     // Slab SST update
     qdw = (sstExt - sst) * cpml / relaxationTimeT;
     HField qioMean = qio * cice; // cice at start of TS, not updated
     HField qowMean = qow * (1 - cice); // 1- cice = open water fraction
     sstSlab = sst - dt * (qioMean + qowMean - qdw) / cpml;
+
     // Slab SSS update
     HField arealDensity = cpml / Water::cp; // density times depth, or cpml divided by cp
     // This is simplified compared to the finiteelement.cpp calculation
     // Fdw = delS * mld * physical::rhow /(timeS*M_sss[i] - ddt*delS) where delS = sssSlab - sssExt
     fdw = (1 - sssExt / sss) * arealDensity / relaxationTimeS;
-    // ice volume change, both laterally and vertically
-    HField deltaIceVol = newIce + deltaHice * cice;
-    // change in snow volume due to melting (should be < 0)
-    HField meltSnowVol = deltaSmelt * cice;
-    // Mass per unit area after all the changes in water volume
-    HField denominator
-        = arealDensity - deltaIceVol * Ice::rho - meltSnowVol * Ice::rhoSnow - (emp - fdw) * dt;
-    // Clamp the denominator to be at least 1 m deep, i.e. at least Water::rho kg m⁻²
-    denominator.clampAbove(Water::rho);
-    // Effective ice salinity is always less than or equal to the SSS
-    HField effectiveIceSal = sss;
-    effectiveIceSal.clampBelow(Ice::s);
-    sssSlab = sss
-        + ((sss - effectiveIceSal) * Ice::rho * deltaIceVol // Change due to ice changes
-              + sss * meltSnowVol
-              + (emp - fdw) * dt) // snow melt, precipitation and nudging fluxes.
-            / denominator;
+
+    /* We use the "virtual salinity flux" approach:
+     *   \partial S / \partial t = F_{fw} S / H,
+     * where S is the ocean salinity, F_{fw} is the freshwater flux and H the mixed-layer depth. In
+     * the presence of sea ice F_{fw} S must be replaced by (S-S_i)I + S(P-E+M_s), with I and S
+     * the fresh-water flux due to ice and snow melt, respectively, and P and E, as precipitation
+     * and evaporation, respectively. Using the salt (F_s) and fresh-water fluxes (F_{fw})
+     * calculated by the IOceanBoundary class, we can write the change in ocean salinity as
+     *   \partial S / \partial t = ( S F_{fw} - 10^3 F_s ) / H
+     */
+    sssSlab = sss + (sss * (fwFlux + fdw) - 1e3 * sFlux) * dt / arealDensity;
 }
 
 } /* namespace Nextsim */

physics/src/include/SlabOcean.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file SlabOcean.hpp
  *
- * @date 7 Sep 2023
+ * @date 30 Aug 2024
  * @author Tim Spain <[email protected]>
  */
 
@@ -38,9 +38,8 @@ class SlabOcean : public ModelComponent, public Configured<SlabOcean> {
         , cice(getStore())
         , qio(getStore())
         , qow(getStore())
-        , newIce(getStore())
-        , deltaHice(getStore())
-        , deltaSmelt(getStore())
+        , fwFlux(getStore())
+        , sFlux(getStore())
     {
     }
 
@@ -82,9 +81,8 @@ class SlabOcean : public ModelComponent, public Configured<SlabOcean> {
     ModelArrayRef<Shared::Q_IO, RW> qio;
     ModelArrayRef<Shared::Q_OW, RW> qow;
     // TODO ModelArrayRef to assimilation flux
-    ModelArrayRef<Shared::NEW_ICE, RW> newIce;
-    ModelArrayRef<Shared::DELTA_HICE, RW> deltaHice;
-    ModelArrayRef<Shared::HSNOW_MELT, RW> deltaSmelt;
+    ModelArrayRef<Protected::FWFLUX, RO> fwFlux;
+    ModelArrayRef<Protected::SFLUX, RO> sFlux;
 
     static const std::string sstSlabName;
     static const std::string sssSlabName;

physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp

@@ -59,7 +59,6 @@ void FiniteElementFluxes::setData(const ModelState::DataMap& ms)
     evap.resize();
     Q_lh_ow.resize();
     Q_sh_ow.resize();
-    Q_sw_ow.resize();
     Q_lw_ow.resize();
     Q_lh_ia.resize();
     Q_sh_ia.resize();
@@ -142,22 +141,29 @@ void FiniteElementFluxes::calculateIce(size_t i, const TimestepTime& tst)
     Q_lh_ia[i] = subl[i] * latentHeatIce(tice.zIndexAndLayer(i, 0));
     double dmdot_dT = dragIce_t * rho_air[i] * v_air[i] * dshice_dT[i];
     double dQlh_dT = latentHeatIce(tice.zIndexAndLayer(i, 0)) * dmdot_dT;
+
     // Sensible heat flux
     Q_sh_ia[i]
         = dragIce_t * rho_air[i] * cp_air[i] * v_air[i] * (tice.zIndexAndLayer(i, 0) - t_air[i]);
     double dQsh_dT = dragIce_t * rho_air[i] * cp_air[i] * v_air[i];
+
     // Shortwave flux
     double albedoValue, i0;
     std::tie(albedoValue, i0)
         = iIceAlbedoImpl->surfaceShortWaveBalance(tice.zIndexAndLayer(i, 0), h_snow_true[i], m_I0);
     Q_sw_ia[i] = -sw_in[i] * (1. - albedoValue) * (1. - i0);
-    penSW[i] = sw_in[i] * (1. - albedoValue) * i0;
+    const double extinction = 0.; // TODO: Replace with de Beer's law or a module
+    penSW[i] = sw_in[i] * (1. - albedoValue) * i0 * (1. - extinction);
+    Q_sw_base[i] = sw_in[i] * (1. - albedoValue) * i0 * extinction;
+
     // Longwave flux
     Q_lw_ia[i] = stefanBoltzmannLaw(tice.zIndexAndLayer(i, 0)) - lw_in[i];
     double dQlw_dT
         = 4 / kelvin(tice.zIndexAndLayer(i, 0)) * stefanBoltzmannLaw(tice.zIndexAndLayer(i, 0));
+
     // Total flux
     qia[i] = Q_lh_ia[i] + Q_sh_ia[i] + Q_sw_ia[i] + Q_lw_ia[i];
+
     // Total temperature dependence of flux
     dqia_dt[i] = dQlh_dT + dQsh_dT + dQlw_dT;
 }

physics/src/modules/FluxCalculationModule/include/FiniteElementFluxes.hpp

@@ -25,7 +25,6 @@ class FiniteElementFluxes : public IFluxCalculation, public Configured<FiniteEle
         , evap(ModelArray::Type::H)
         , Q_lh_ow(ModelArray::Type::H)
         , Q_sh_ow(ModelArray::Type::H)
-        , Q_sw_ow(ModelArray::Type::H)
         , Q_lw_ow(ModelArray::Type::H)
         , Q_lh_ia(ModelArray::Type::H)
         , Q_sh_ia(ModelArray::Type::H)
@@ -89,7 +88,6 @@ class FiniteElementFluxes : public IFluxCalculation, public Configured<FiniteEle
     HField evap; // Open water evaporative mass flux [kg  m⁻²]
     HField Q_lh_ow; // Open water latent heat flux [W m⁻²]
     HField Q_sh_ow; // Open water sensible heat flux [W m⁻²]
-    HField Q_sw_ow; // Open water incident shortwave radiative flux [W m⁻²]
     HField Q_lw_ow; // Open water net longwave radiative flux [W m⁻²]
     HField Q_lh_ia; // Ice latent heat flux [W m⁻²]
     HField Q_sh_ia; // Ice sensible heat flux [W m⁻²]

physics/src/modules/OceanBoundaryModule/ConfiguredOcean.cpp

@@ -108,6 +108,9 @@ void ConfiguredOcean::updateBefore(const TimestepTime& tst)
 
 void ConfiguredOcean::updateAfter(const TimestepTime& tst)
 {
+    overElements(
+        std::bind(&IOceanBoundary::mergeFluxes, this, std::placeholders::_1, std::placeholders::_2),
+        tst);
     slabOcean.update(tst);
     sst = ModelArrayRef<Protected::SLAB_SST, RO>(getStore()).data();
     sss = ModelArrayRef<Protected::SLAB_SSS, RO>(getStore()).data();

physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp

@@ -10,6 +10,7 @@
 #include "include/Finalizer.hpp"
 #include "include/IIceOceanHeatFlux.hpp"
 #include "include/IFreezingPoint.hpp"
+#include "include/IIceOceanHeatFlux.hpp"
 #include "include/NextsimModule.hpp"
 #include "include/ParaGridIO.hpp"
 #include "include/constants.hpp"
@@ -76,6 +77,9 @@ void TOPAZOcean::updateBefore(const TimestepTime& tst)
 
 void TOPAZOcean::updateAfter(const TimestepTime& tst)
 {
+    overElements(
+        std::bind(&IOceanBoundary::mergeFluxes, this, std::placeholders::_1, std::placeholders::_2),
+        tst);
     slabOcean.update(tst);
     sst = ModelArrayRef<Protected::SLAB_SST, RO>(getStore()).data();
     sss = ModelArrayRef<Protected::SLAB_SSS, RO>(getStore()).data();

physics/src/modules/include/IAtmosphereBoundary.hpp

@@ -38,6 +38,8 @@ class IAtmosphereBoundary : public ModelComponent {
         , uwind(ModelArray::Type::U)
         , vwind(ModelArray::Type::V)
         , penSW(ModelArray::Type::H)
+        , qswow(ModelArray::Type::H)
+        , qswBase(ModelArray::Type::H)
     {
         m_couplingArrays.registerArray(CouplingFields::SUBL, &subl, RW);
         m_couplingArrays.registerArray(CouplingFields::SNOW, &snow, RW);
@@ -55,6 +57,8 @@ class IAtmosphereBoundary : public ModelComponent {
         getStore().registerArray(Protected::WIND_U, &uwind, RO);
         getStore().registerArray(Protected::WIND_V, &vwind, RO);
         getStore().registerArray(Shared::Q_PEN_SW, &penSW, RW);
+        getStore().registerArray(Shared::Q_SW_OW, &qswow, RW);
+        getStore().registerArray(Shared::Q_SW_BASE, &qswBase, RW);
     }
     virtual ~IAtmosphereBoundary() = default;
 
@@ -75,6 +79,8 @@ class IAtmosphereBoundary : public ModelComponent {
         uwind.resize();
         vwind.resize();
         penSW.resize();
+        qswow.resize();
+        qswBase.resize();
     }
     virtual void update(const TimestepTime& tst) { }
 
@@ -92,6 +98,8 @@ class IAtmosphereBoundary : public ModelComponent {
     UField uwind;
     VField vwind;
     HField penSW;
+    HField qswow;
+    HField qswBase;
 
     ModelArrayReferenceStore m_couplingArrays;
 };

physics/src/modules/include/IFluxCalculation.hpp

@@ -23,6 +23,8 @@ class IFluxCalculation : public ModelComponent {
         , qia(getStore())
         , penSW(getStore())
         , dqia_dt(getStore())
+        , Q_sw_ow(getStore())
+        , Q_sw_base(getStore())
     {
     }
     virtual ~IFluxCalculation() = default;
@@ -53,6 +55,8 @@ class IFluxCalculation : public ModelComponent {
     ModelArrayRef<Shared::Q_IA, RW> qia;
     ModelArrayRef<Shared::Q_PEN_SW, RW> penSW;
     ModelArrayRef<Shared::DQIA_DT, RW> dqia_dt;
+    ModelArrayRef<Shared::Q_SW_OW, RW> Q_sw_ow;
+    ModelArrayRef<Shared::Q_SW_BASE, RW> Q_sw_base;
 };
 }
 #endif /* IFLUXCALCULATION_HPP */

physics/src/modules/include/IIceThermodynamics.hpp

@@ -64,6 +64,7 @@ class IIceThermodynamics : public ModelComponent {
         , tf(getStore())
         , snowfall(getStore())
         , sss(getStore())
+        , qswBase(getStore())
     {
         registerModule();
 
@@ -76,6 +77,7 @@ class IIceThermodynamics : public ModelComponent {
     ModelArrayRef<Shared::H_SNOW, RW> hsnow; // From Ice Growth
     ModelArrayRef<Shared::Q_IC, RW>
         qic; // From IceTemperature. Conductive heat flux to the ice surface.
+    ModelArrayRef<Shared::Q_SW_BASE, RW> qswBase; // Short-wave flux through the base of the ice
     ModelArrayRef<Shared::Q_IO, RW> qio; // From FluxCalculation
     ModelArrayRef<Shared::Q_IA, RO> qia; // From FluxCalculation
     ModelArrayRef<Shared::DQIA_DT, RO> dQia_dt; // From FluxCalculation

physics/src/modules/include/IOceanBoundary.hpp

@@ -9,6 +9,7 @@
 #define IOCEANBOUNDARY_HPP
 
 #include "include/ModelComponent.hpp"
+#include "include/constants.hpp"
 
 namespace Nextsim {
 
@@ -23,6 +24,14 @@ namespace CouplingFields {
 class IOceanBoundary : public ModelComponent {
 public:
     IOceanBoundary()
+        : cice(getStore())
+        , qswow(getStore())
+        , emp(getStore())
+        , newIce(getStore())
+        , deltaHice(getStore())
+        , deltaSmelt(getStore())
+        , qow(getStore())
+        , qswBase(getStore())
     {
         m_couplingArrays.registerArray(CouplingFields::SST, &sst, RW);
         m_couplingArrays.registerArray(CouplingFields::SSS, &sss, RW);
@@ -37,6 +46,8 @@ class IOceanBoundary : public ModelComponent {
         getStore().registerArray(Protected::TF, &tf, RO);
         getStore().registerArray(Protected::OCEAN_U, &u, RO);
         getStore().registerArray(Protected::OCEAN_V, &v, RO);
+        getStore().registerArray(Protected::FWFLUX, &fwFlux, RO);
+        getStore().registerArray(Protected::SFLUX, &sFlux, RO);
     }
     virtual ~IOceanBoundary() = default;
 
@@ -54,6 +65,10 @@ class IOceanBoundary : public ModelComponent {
         tf.resize();
         u.resize();
         v.resize();
+        qNoSun.resize();
+        qswNet.resize();
+        fwFlux.resize();
+        sFlux.resize();
 
         if (ms.count("sst")) {
             sst = ms.at("sst");
@@ -77,6 +92,32 @@ class IOceanBoundary : public ModelComponent {
      */
     virtual void updateAfter(const TimestepTime& tst) = 0;
 
+    /*!
+     * Merges the ice-ocean fluxes and ocean-atmosphere fluxes into a single field to be passed to a
+     * slab-ocean implementation or an ocean model through a coupler.
+     */
+    void mergeFluxes(size_t i, const TimestepTime& tst)
+    {
+        const double dt = tst.step.seconds();
+
+        qswNet(i) = cice(i) * qswBase(i) + (1 - cice(i)) * qswow(i);
+        qNoSun(i) = cice(i) * qio(i) + (1 - cice(i)) * qow(i) - qswNet(i);
+
+        // ice volume change, both laterally and vertically
+        const double deltaIceVol = newIce(i) + deltaHice(i) * cice(i);
+        // change in snow volume due to melting (should be < 0)
+        const double meltSnowVol = deltaSmelt(i) * cice(i);
+        // Effective ice salinity is always less than or equal to the SSS, and here we use the right
+        // units too
+        const double effectiveIceSal = 1e-3 * std::min(sss(i), Ice::s);
+
+        // Positive flux is up!
+        fwFlux(i)
+            = ((1 - effectiveIceSal) * Ice::rho * deltaIceVol + Ice::rhoSnow * meltSnowVol) / dt
+            + emp(i) * (1 - cice(i));
+        sFlux(i) = effectiveIceSal * Ice::rho * deltaIceVol / dt;
+    }
+
 protected:
     HField qio; // Ice-ocean heat flux, W m⁻²
     HField sst; // Coupled or slab ocean sea surface temperature, ˚C
@@ -86,8 +127,21 @@ class IOceanBoundary : public ModelComponent {
     HField cpml; // Heat capacity of the mixed layer, J K⁻¹ m²
     UField u; // x(east)-ward ocean current, m s⁻¹
     VField v; // y(north)-ward ocean current, m s⁻¹
+    HField qNoSun; // Net surface ocean heat flux, except short wave, W m⁻²
+    HField qswNet; // Net surface ocean shortwave flux, W m⁻²
+    HField fwFlux; // Net surface ocean fresh-water flux, kg m⁻²
+    HField sFlux; // Net surface ocean salt flux, kg m⁻²
 
     ModelArrayReferenceStore m_couplingArrays;
+
+    ModelArrayRef<Protected::C_ICE, RO> cice;
+    ModelArrayRef<Shared::Q_SW_OW, RO> qswow;
+    ModelArrayRef<Protected::EVAP_MINUS_PRECIP, RO> emp;
+    ModelArrayRef<Shared::NEW_ICE, RW> newIce;
+    ModelArrayRef<Shared::DELTA_HICE, RW> deltaHice;
+    ModelArrayRef<Shared::HSNOW_MELT, RW> deltaSmelt;
+    ModelArrayRef<Shared::Q_OW, RW> qow;
+    ModelArrayRef<Shared::Q_SW_BASE, RW> qswBase;
 };
 } /* namespace Nextsim */