Define fluxes needed for the ocean coupling

core/src/include/ModelComponent.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file ModelComponent.hpp
  *
- * @date 06 Dec 2024
+ * @date 10 Feb 2025
  * @author Tim Spain <[email protected]>
  * @author Einar Ólason <[email protected]>
  */
@@ -91,6 +91,8 @@ 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"; // Shortwave flux in open water W m⁻²
+    inline constexpr TextTag Q_SW_BASE = "Q_SW_BASE"; // 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
@@ -101,6 +103,31 @@ namespace Shared {
     inline constexpr TextTag NEW_ICE = "NEW_ICE"; // Volume of new ice formed [m]
 
 }
+
+namespace CouplingFields {
+    /* Ice-ocean coupling */
+    // Receive
+    inline constexpr TextTag MLD = "MLD"; // Mixed layer or slab ocean depth m
+    inline constexpr TextTag OCEAN_U = "U"; // x(east)-ward ocean current m s⁻¹
+    inline constexpr TextTag OCEAN_V = "V"; // y(north)-ward ocean current m s⁻¹
+    inline constexpr TextTag SSH = "SSH"; // sea surface height, m
+    inline constexpr TextTag SSS = "SSS"; // sea surface salinity PSU
+    inline constexpr TextTag SST = "SST"; // sea surface temperature ˚C
+    // Send
+    inline constexpr TextTag FWFLUX = "FWFLUX"; // Fresh-water flux at the ocean surface, kg m⁻²
+    inline constexpr TextTag Q_SS_NO_SW = "Q_SS_NO_SW"; // Net non-solar flux at ocean surface W m⁻²
+    inline constexpr TextTag Q_SS_SW = "Q_SS_SW"; // Net shortwave flux at ocean surface W m⁻²
+    inline constexpr TextTag SFLUX = "SFLUX"; // Salt flux at the ocean surface, kg m⁻²
+
+    /* Ice-atmosphere coupling */
+    inline constexpr TextTag EVAP = "EVAP"; // evaporation mass flux kg s⁻¹ m⁻²
+    inline constexpr TextTag RAIN = "RAIN"; // rainfall mass flux kg s⁻¹ m⁻²
+    inline constexpr TextTag SNOW = "SNOW"; // snowfall mass flux kg s⁻¹ m⁻²
+    inline constexpr TextTag SUBL = "SUBL"; // sublimation mass flux kg s⁻¹ m⁻²
+    inline constexpr TextTag WIND_U = "WIND_U"; // x-aligned wind component m s⁻¹
+    inline constexpr TextTag WIND_V = "WIND_V"; // y-aligned wind component m s⁻¹
+}
+
 /*!
  * A class encapsulating a component of the model. It also provide a method of
  * communicating data between ModelComponents using enums, static arrays of

physics/src/SlabOcean.cpp

@@ -1,7 +1,7 @@
 /*!
  * @file SlabOcean.cpp
  *
- * @date 20 Nov 2024
+ * @date 10 Feb 2025
  * @author Tim Spain <[email protected]>
  */
 
@@ -73,34 +73,22 @@ std::unordered_set<std::string> SlabOcean::hFields() const { return { sstSlabNam
 
 void SlabOcean::update(const TimestepTime& tst)
 {
-    double dt = tst.step.seconds();
+    const 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;
+    sstSlab = sst - dt * (qswNet + qNoSun - qdw) / cpml;
+
     // Slab SSS update
-    HField arealDensity = cpml / Water::cp; // density times depth, or cpml divided by cp
+    const 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;
+    const HField denominator = (arealDensity - (fwFlux - fdw) * dt).clampAbove(Water::rhoOcean);
+    sssSlab = sss + (sss * fwFlux - fdw * dt) / denominator;
 }
 
 } /* namespace Nextsim */

physics/src/include/SlabOcean.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file SlabOcean.hpp
  *
- * @date 7 Sep 2023
+ * @date 10 Feb 2025
  * @author Tim Spain <[email protected]>
  */
 
@@ -23,7 +23,7 @@ namespace Nextsim {
  */
 class SlabOcean : public ModelComponent, public Configured<SlabOcean> {
 public:
-    SlabOcean()
+    SlabOcean(ModelArrayReferenceStore& coupingArrays)
         : qdw(ModelArray::Type::H)
         , fdw(ModelArray::Type::H)
         , sstSlab(ModelArray::Type::H)
@@ -32,15 +32,11 @@ class SlabOcean : public ModelComponent, public Configured<SlabOcean> {
         , sssExt(getStore())
         , sst(getStore())
         , sss(getStore())
-        , mld(getStore())
         , cpml(getStore())
-        , emp(getStore())
-        , cice(getStore())
-        , qio(getStore())
-        , qow(getStore())
-        , newIce(getStore())
-        , deltaHice(getStore())
-        , deltaSmelt(getStore())
+        , qswNet(coupingArrays)
+        , qNoSun(coupingArrays)
+        , fwFlux(coupingArrays)
+        , sFlux(coupingArrays)
     {
     }
 
@@ -75,16 +71,12 @@ class SlabOcean : public ModelComponent, public Configured<SlabOcean> {
     ModelArrayRef<Protected::EXT_SSS> sssExt;
     ModelArrayRef<Protected::SST> sst;
     ModelArrayRef<Protected::SSS> sss;
-    ModelArrayRef<Protected::MLD> mld;
     ModelArrayRef<Protected::ML_BULK_CP> cpml;
-    ModelArrayRef<Protected::EVAP_MINUS_PRECIP> emp;
-    ModelArrayRef<Protected::C_ICE> cice;
-    ModelArrayRef<Shared::Q_IO, RW> qio;
-    ModelArrayRef<Shared::Q_OW, RW> qow;
+    ModelArrayRef<CouplingFields::Q_SS_SW, RO> qswNet;
+    ModelArrayRef<CouplingFields::Q_SS_NO_SW, RO> qNoSun;
+    ModelArrayRef<CouplingFields::FWFLUX, RO> fwFlux;
+    ModelArrayRef<CouplingFields::SFLUX, RO> sFlux;
     // TODO ModelArrayRef to assimilation flux
-    ModelArrayRef<Shared::NEW_ICE, RW> newIce;
-    ModelArrayRef<Shared::DELTA_HICE, RW> deltaHice;
-    ModelArrayRef<Shared::HSNOW_MELT, RW> deltaSmelt;
 
     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

@@ -1,7 +1,7 @@
 /*!
  * @file ConfiguredOcean.cpp
  *
- * @date 06 Dec 2024
+ * @date 10 Feb 2025
  * @author Tim Spain <[email protected]>
  */
 
@@ -40,6 +40,7 @@ static const std::map<int, std::string> keyMap = {
 ConfiguredOcean::ConfiguredOcean()
     : sstExt(ModelArray::Type::H)
     , sssExt(ModelArray::Type::H)
+    , slabOcean(m_couplingArrays)
 {
 }
 
@@ -112,6 +113,7 @@ void ConfiguredOcean::updateBefore(const TimestepTime& tst)
 
 void ConfiguredOcean::updateAfter(const TimestepTime& tst)
 {
+    mergeFluxes(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

@@ -1,7 +1,7 @@
 /*!
  * @file TOPAZOcean.cpp
  *
- * @date 06 Dec 2024
+ * @date 10 Feb 2025
  * @author Tim Spain <[email protected]>
  */
 
@@ -29,6 +29,7 @@ static const std::map<int, std::string> keyMap = {
 TOPAZOcean::TOPAZOcean()
     : sstExt(ModelArray::Type::H)
     , sssExt(ModelArray::Type::H)
+    , slabOcean(m_couplingArrays)
 {
 }
 
@@ -71,7 +72,7 @@ void TOPAZOcean::updateBefore(const TimestepTime& tst)
         ssh = 0.;
     }
 
-    cpml = Water::rho * Water::cp * mld;
+    cpml = Water::rhoOcean * Water::cp * mld;
     overElements(
         std::bind(&TOPAZOcean::updateTf, this, std::placeholders::_1, std::placeholders::_2),
         TimestepTime());
@@ -81,6 +82,7 @@ void TOPAZOcean::updateBefore(const TimestepTime& tst)
 
 void TOPAZOcean::updateAfter(const TimestepTime& tst)
 {
+    mergeFluxes(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

@@ -1,7 +1,7 @@
 /*!
  * @file IAtmosphereBoundary.hpp
  *
- * @date 24 Sep 2024
+ * @date 10 Feb 2025
  * @author Tim Spain <[email protected]>
  */
 
@@ -14,15 +14,6 @@
 
 namespace Nextsim {
 
-namespace CouplingFields {
-    constexpr TextTag SUBL = "SUBL"; // sublimation mass flux kg s⁻¹ m⁻²
-    constexpr TextTag SNOW = "SNOW"; // snowfall mass flux kg s⁻¹ m⁻²
-    constexpr TextTag RAIN = "RAIN"; // rainfall mass flux kg s⁻¹ m⁻²
-    constexpr TextTag EVAP = "EVAP"; // evaporation mass flux kg s⁻¹ m⁻²
-    constexpr TextTag WIND_U = "WIND_U"; // x-aligned wind component m s⁻¹
-    constexpr TextTag WIND_V = "WIND_V"; // y-aligned wind component m s⁻¹
-
-}
 //! An interface class for the atmospheric inputs into the ice physics.
 class IAtmosphereBoundary : public ModelComponent {
 public:

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