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.
nextsimhub · Sep 2, 2024 · 0da304c · 0da304c
1 parent 6eedb7e
commit 0da304c
Show file tree

Hide file tree

Showing 11 changed files with 142 additions and 71 deletions.
diff --git a/core/src/include/ModelComponent.hpp b/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

diff --git a/physics/src/SlabOcean.cpp b/physics/src/SlabOcean.cpp
@@ -1,7 +1,7 @@
 /*!
  * @file SlabOcean.cpp
  *
- * @date 7 Sep 2023
+ * @date 30 Aug 2024
  * @author Tim Spain <[email protected]>
  */
 
@@ -75,33 +75,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 */
diff --git a/physics/src/include/SlabOcean.hpp b/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;

diff --git a/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp b/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp
@@ -1,7 +1,7 @@
 /*!
  * @file FiniteElementFluxes.cpp
  *
- * @date Apr 29, 2022
+ * @date 02 Sep 2024
  * @author Tim Spain <[email protected]>
  */
 
@@ -56,7 +56,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();
@@ -137,22 +136,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;
 }

diff --git a/physics/src/modules/FluxCalculationModule/include/FiniteElementFluxes.hpp b/physics/src/modules/FluxCalculationModule/include/FiniteElementFluxes.hpp
@@ -1,18 +1,18 @@
 /*!
  * @file FiniteElementFluxes.hpp
  *
- * @date Apr 29, 2022
+ * @date 02 Sep 2024
  * @author Tim Spain <[email protected]>
  */
 
 #ifndef FINITEELEMENTFLUXES_HPP
 #define FINITEELEMENTFLUXES_HPP
 
-#include "include/ModelArrayRef.hpp"
 #include "include/Configured.hpp"
 #include "include/IFluxCalculation.hpp"
 #include "include/IIceAlbedo.hpp"
 #include "include/IIceOceanHeatFlux.hpp"
+#include "include/ModelArrayRef.hpp"
 #include "include/ModelComponent.hpp"
 
 namespace Nextsim {
@@ -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⁻²]
@@ -111,8 +109,7 @@ class FiniteElementFluxes : public IFluxCalculation, public Configured<FiniteEle
     ModelArrayRef<Protected::P_AIR> p_air;
     ModelArrayRef<Protected::WIND_SPEED> v_air;
     ModelArrayRef<Protected::H_SNOW> h_snow; // cell-averaged value
-    ModelArrayRef<Protected::HTRUE_SNOW>
-        h_snow_true; // cell-averaged value
+    ModelArrayRef<Protected::HTRUE_SNOW> h_snow_true; // cell-averaged value
     ModelArrayRef<Protected::C_ICE> cice;
     ModelArrayRef<Protected::T_ICE> tice;
     ModelArrayRef<Protected::SW_IN> sw_in;

diff --git a/physics/src/modules/OceanBoundaryModule/ConfiguredOcean.cpp b/physics/src/modules/OceanBoundaryModule/ConfiguredOcean.cpp
@@ -1,7 +1,7 @@
 /*!
  * @file ConfiguredOcean.cpp
  *
- * @date 7 Sep 2023
+ * @date 30 Aug 2024
  * @author Tim Spain <[email protected]>
  */
 
@@ -110,9 +110,11 @@ 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();
-
 }
 } /* namespace Nextsim */
diff --git a/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp b/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp
@@ -1,14 +1,14 @@
 /*!
  * @file TOPAZOcean.cpp
  *
- * @date 7 Sep 2023
+ * @date 30 Aug 2024
  * @author Tim Spain <[email protected]>
  */
 
 #include "include/TOPAZOcean.hpp"
 
-#include "include/IIceOceanHeatFlux.hpp"
 #include "include/IFreezingPoint.hpp"
+#include "include/IIceOceanHeatFlux.hpp"
 #include "include/Module.hpp"
 #include "include/ParaGridIO.hpp"
 #include "include/constants.hpp"
@@ -49,7 +49,6 @@ void TOPAZOcean::configure()
 
     getStore().registerArray(Protected::EXT_SST, &sstExt, RO);
     getStore().registerArray(Protected::EXT_SSS, &sssExt, RO);
-
 }
 
 void TOPAZOcean::updateBefore(const TimestepTime& tst)
@@ -65,22 +64,23 @@ void TOPAZOcean::updateBefore(const TimestepTime& tst)
     v = state.data.at("v");
 
     cpml = Water::rho * Water::cp * mld;
-    overElements(std::bind(&TOPAZOcean::updateTf, this, std::placeholders::_1,
-                     std::placeholders::_2),
+    overElements(
+        std::bind(&TOPAZOcean::updateTf, this, std::placeholders::_1, std::placeholders::_2),
         TimestepTime());
 
     Module::getImplementation<IIceOceanHeatFlux>().update(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();
 }
 
-
 void TOPAZOcean::setFilePath(const std::string& filePathIn) { filePath = filePathIn; }
 
 void TOPAZOcean::setData(const ModelState::DataMap& ms)

diff --git a/physics/src/modules/include/IAtmosphereBoundary.hpp b/physics/src/modules/include/IAtmosphereBoundary.hpp
@@ -1,7 +1,7 @@
 /*!
  * @file IAtmosphereBoundary.hpp
  *
- * @date Sep 22, 2022
+ * @date 30 Aug 2024
  * @author Tim Spain <[email protected]>
  */
 
@@ -15,12 +15,12 @@
 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⁻¹
+    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.
@@ -29,15 +29,17 @@ class IAtmosphereBoundary : public ModelComponent {
     IAtmosphereBoundary()
         : qia(ModelArray::Type::H)
         , dqia_dt(ModelArray::Type::H)
-    , qow(ModelArray::Type::H)
-    , subl(ModelArray::Type::H)
-    , snow(ModelArray::Type::H)
-    , rain(ModelArray::Type::H)
-    , evap(ModelArray::Type::H)
-    , emp(ModelArray::Type::H)
-    , uwind(ModelArray::Type::U)
-    , vwind(ModelArray::Type::V)
-    , penSW(ModelArray::Type::H)
+        , qow(ModelArray::Type::H)
+        , subl(ModelArray::Type::H)
+        , snow(ModelArray::Type::H)
+        , rain(ModelArray::Type::H)
+        , evap(ModelArray::Type::H)
+        , emp(ModelArray::Type::H)
+        , 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,11 +79,12 @@ class IAtmosphereBoundary : public ModelComponent {
         uwind.resize();
         vwind.resize();
         penSW.resize();
+        qswow.resize();
+        qswBase.resize();
     }
     virtual void update(const TimestepTime& tst) { }
 
 protected:
-
     ModelArrayReferenceStore& couplingArrays() { return m_couplingArrays; }
 
     HField qia;
@@ -93,6 +98,8 @@ class IAtmosphereBoundary : public ModelComponent {
     UField uwind;
     VField vwind;
     HField penSW;
+    HField qswow;
+    HField qswBase;
 
     ModelArrayReferenceStore m_couplingArrays;
 };