Implement, document, and add a test for -ocean ...,delta_MBP

CHANGES.rst

@@ -32,6 +32,9 @@ Changes since v1.2
 - Remove `ocean.always_grounded`. Set `sea_level.constant.value` to a large negative value
   to ensure that all ice is grounded.
 - Remove `ocean.melange_back_pressure_fraction`: it is no longer needed.
+- Add a new ocean modifier: `-ocean ...,delta_MBP`. This component reads scalar
+  time-dependent melange pressure offsets (units: Pa) and uses them in the calving front
+  boundary condition for the SSA.
 
 Changes from v1.1 to v1.2
 =========================

doc/sphinx/climate_forcing/ocean.rst

@@ -311,10 +311,27 @@ It takes the following command-line options:
   which is similar, but applies anomalies at the atmosphere or surface level,
   respectively.
 
+.. _sec-ocean-delta-mbp:
+
+Scalar melange back pressure offsets
+++++++++++++++++++++++++++++++++++++
+
+:|options|: :opt:`-ocean ...,delta_MBP`
+:|variables|: :var:`delta_MBP` [Pascal]
+:|implementation|: ``pism::ocean::Delta_MBP``
+
+The scalar time-dependent variable :var:`delta_MBP` (units: Pascal) has the meaning of
+`\bar p_{\text{melange}}` in :ref:`sec-model-melange-pressure` (equation :eq:`eq-cfbc-3`).
+
+This modifier uses the following configuration parameters:
+
+.. pism-parameters::
+   :prefix: ocean.delta_MBP.
+
 .. _sec-ocean-frac-mbp:
 
-Scalar melange back pressure fraction
-+++++++++++++++++++++++++++++++++++++
+Melange back pressure as a fraction of pressure difference
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 :|options|: :opt:`-ocean ...,frac_MBP`
 :|variables|: :var:`frac_MBP`

src/coupler/CMakeLists.txt

@@ -37,6 +37,7 @@ add_library (boundary OBJECT
   ./ocean/Anomaly.cc
   ./ocean/Delta_T.cc
   ./ocean/Delta_SMB.cc
+  ./ocean/Delta_MBP.cc
   ./ocean/Frac_MBP.cc
   ./ocean/Frac_SMB.cc
   ./ocean/Runoff_SMB.cc

src/coupler/ocean/Delta_MBP.cc

@@ -0,0 +1,67 @@
+/* Copyright (C) 2021 PISM Authors
+ *
+ * This file is part of PISM.
+ *
+ * PISM is free software; you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation; either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PISM; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+#include "Delta_MBP.hh"
+#include "pism/coupler/util/ScalarForcing.hh"
+
+namespace pism {
+namespace ocean {
+
+Delta_MBP::Delta_MBP(IceGrid::ConstPtr g, std::shared_ptr<OceanModel> in)
+  : OceanModel(g, in) {
+
+  m_forcing.reset(new ScalarForcing(g->ctx(),
+                                    "ocean.delta_MBP",
+                                    "delta_MBP",
+                                    "Pa", "Pa",
+                                    "melange back pressure"));
+
+  m_water_column_pressure = allocate_water_column_pressure(g);
+}
+
+Delta_MBP::~Delta_MBP() {
+  // empty
+}
+
+void Delta_MBP::init_impl(const Geometry &geometry) {
+
+  m_input_model->init(geometry);
+
+  m_log->message(2, "* Initializing melange back pressure forcing using scalar offsets...\n");
+
+  m_forcing->init();
+}
+
+void Delta_MBP::update_impl(const Geometry &geometry, double t, double dt) {
+  m_input_model->update(geometry, t, dt);
+
+  m_forcing->update(t, dt);
+
+  m_water_column_pressure->copy_from(m_input_model->average_water_column_pressure());
+
+  m_water_column_pressure->shift(m_forcing->value());
+}
+
+const IceModelVec2S& Delta_MBP::average_water_column_pressure_impl() const {
+  return *m_water_column_pressure;
+}
+
+
+} // end of namespace ocean
+} // end of namespace pism

src/coupler/ocean/Delta_MBP.hh

@@ -0,0 +1,53 @@
+/* Copyright (C) 2013, 2014, 2015, 2016, 2017, 2018, 2021 PISM Authors
+ *
+ * This file is part of PISM.
+ *
+ * PISM is free software; you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation; either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PISM; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+#ifndef _DELTA_MBP_H_
+#define _DELTA_MBP_H_
+
+#include "pism/coupler/OceanModel.hh"
+
+namespace pism {
+
+class ScalarForcing;
+
+namespace ocean {
+
+/**
+ * Scalar melange back-pressure offsets.
+ */
+class Delta_MBP : public OceanModel
+{
+public:
+  Delta_MBP(IceGrid::ConstPtr g, std::shared_ptr<OceanModel> in);
+  virtual ~Delta_MBP();
+
+private:
+  void init_impl(const Geometry &geometry);
+
+  void update_impl(const Geometry &geometry, double t, double dt);
+
+  const IceModelVec2S& average_water_column_pressure_impl() const;
+
+  std::unique_ptr<ScalarForcing> m_forcing;
+};
+
+} // end of namespace ocean
+} // end of namespace pism
+
+#endif /* _DELTA_MBP_H_ */

src/coupler/ocean/Factory.cc

@@ -1,4 +1,4 @@
-/* Copyright (C) 2015, 2017, 2018, 2019 PISM Authors
+/* Copyright (C) 2015, 2017, 2018, 2019, 2021 PISM Authors
  *
  * This file is part of PISM.
  *
@@ -26,6 +26,7 @@
 #include "GivenClimate.hh"
 #include "Delta_T.hh"
 #include "Delta_SMB.hh"
+#include "Delta_MBP.hh"
 #include "Frac_MBP.hh"
 #include "Frac_SMB.hh"
 #include "Runoff_SMB.hh"
@@ -51,6 +52,7 @@ Factory::Factory(IceGrid::ConstPtr g)
   add_modifier<Frac_SMB>("frac_SMB");
   add_modifier<Delta_T>("delta_T");
   add_modifier<Runoff_SMB>("runoff_SMB");
+  add_modifier<Delta_MBP>("delta_MBP");
   add_modifier<Frac_MBP>("frac_MBP");
 }
 

src/coupler/ocean/OceanModel.cc

@@ -51,8 +51,8 @@ IceModelVec2S::Ptr OceanModel::allocate_water_column_pressure(IceGrid::ConstPtr
                                               "average_water_column_pressure",
                                               WITHOUT_GHOSTS));
   result->set_attrs("diagnostic",
-                    "vertically-integrated water column pressure",
-                    "Pa m", "Pa m", "", 0);
+                    "vertically-averaged water column pressure",
+                    "Pa", "Pa", "", 0);
 
   return result;
 }

src/pism_config.cdl

@@ -1606,6 +1606,20 @@ netcdf pism_config {
     pism_config:ocean.constant.melt_rate_type = "number";
     pism_config:ocean.constant.melt_rate_units = "m / year";
 
+    pism_config:ocean.delta_MBP.file = "";
+    pism_config:ocean.delta_MBP.file_doc = "Name of the file containing melange back-pressure offsets";
+    pism_config:ocean.delta_MBP.file_type = "string";
+
+    pism_config:ocean.delta_MBP.period = 0;
+    pism_config:ocean.delta_MBP.period_doc = "Length of the period of the climate forcing data (set to zero to disable)";
+    pism_config:ocean.delta_MBP.period_type = "integer";
+    pism_config:ocean.delta_MBP.period_units = "years";
+
+    pism_config:ocean.delta_MBP.reference_year = 0;
+    pism_config:ocean.delta_MBP.reference_year_doc = "Reference year to use when :config:`ocean.delta_MBP.period` is active";
+    pism_config:ocean.delta_MBP.reference_year_type = "integer";
+    pism_config:ocean.delta_MBP.reference_year_units = "years";
+
     pism_config:ocean.delta_T.file = "";
     pism_config:ocean.delta_T.file_doc = "Name of the file containing temperature offsets.";
     pism_config:ocean.delta_T.file_option = "ocean_delta_T_file";

src/pythonbindings/pism_ocean.i

@@ -92,3 +92,6 @@
 %shared_ptr(pism::ocean::sea_level::Factory)
 %rename(SeaLevelFactory) pism::ocean::sea_level::Factory;
 %include "coupler/ocean/sea_level/Factory.hh"
+
+%rename(OceanDeltaMBP) pism::ocean::Delta_MBP;
+pism_class(pism::ocean::Delta_MBP, "pism/coupler/ocean/Delta_MBP.hh")

test/regression/ocean_models.py

@@ -310,6 +310,48 @@ def test_ocean_anomaly(self):
     def tearDown(self):
         os.remove(self.filename)
 
+class DeltaMBP(TestCase):
+    def setUp(self):
+        self.filename = "ocean_delta_MBP_input.nc"
+        self.grid = shallow_grid()
+        self.geometry = PISM.Geometry(self.grid)
+        self.model = PISM.OceanConstant(self.grid)
+        self.dP = 100.0
+        self.H = 1000.0
+
+        self.geometry.ice_thickness.set(self.H)
+        self.geometry.bed_elevation.set(-2 * self.H)
+        self.geometry.ensure_consistency(0.0)
+
+        create_scalar_forcing(self.filename, "delta_MBP", "Pa", [self.dP], [0])
+
+    def test_ocean_delta_mpb(self):
+        "Modifier Delta_MBP"
+
+        modifier = PISM.OceanDeltaMBP(self.grid, self.model)
+
+        config.set_string("ocean.delta_MBP.file", self.filename)
+
+        modifier.init(self.geometry)
+        modifier.update(self.geometry, 0, 1)
+
+        model = self.model
+
+        check_difference(modifier.shelf_base_temperature(),
+                         model.shelf_base_temperature(),
+                         0.0)
+
+        check_difference(modifier.shelf_base_mass_flux(),
+                         model.shelf_base_mass_flux(),
+                         0.0)
+
+        check_difference(modifier.average_water_column_pressure(),
+                         model.average_water_column_pressure(),
+                         self.dP)
+
+    def tearDown(self):
+        os.remove(self.filename)
+
 class FracMBP(TestCase):
     def setUp(self):
         self.filename = "ocean_frac_MBP_input.nc"