test up OHC restarts

src/fair/fair.py

@@ -378,8 +378,9 @@ def allocate(self):
             coords=(self.timebounds, self.scenarios, self.configs, self.species),
             dims=("timebounds", "scenario", "config", "specie"),
         )
+        # init with NaNs better than zeros, but makes code cleaner later
         self.ocean_heat_content_change = xr.DataArray(
-            np.ones((self._n_timebounds, self._n_scenarios, self._n_configs)) * np.nan,
+            np.zeros((self._n_timebounds, self._n_scenarios, self._n_configs)),
             coords=(self.timebounds, self.scenarios, self.configs),
             dims=("timebounds", "scenario", "config"),
         )
@@ -1941,8 +1942,15 @@ def run(self, progress=True, suppress_warnings=True):
         )
 
         # 18. Ocean heat content change
-        ocean_heat_content_change_array = (
-            np.cumsum(toa_imbalance_array * self.timestep, axis=TIME_AXIS)
+        # if a restart value is present, include it and subtract off the TOA
+        # imbalance from the restart as this would be double counting
+        # in non-restart runs, both OHC and N are zero in first timebound
+        ocean_heat_content_change_array = self.ocean_heat_content_change[0:1, ...] + (
+            (
+                np.cumsum(toa_imbalance_array, axis=TIME_AXIS)
+                - toa_imbalance_array[0:1, ...]
+            )
+            * self.timestep
             * earth_radius**2
             * 4
             * np.pi

tests/integration_tests/restart_test.py

@@ -1,13 +1,14 @@
 """Module for integration test of fair."""
 
 import copy
-from fair.interface import initialise
+
 import numpy as np
+
 from fair import FAIR
+from fair.interface import initialise
 from fair.io import read_properties
 
 
-
 def minimal_ghg_run(start=1750, end=2020, timestep=27, stochastic_run=False, seed=37):
     fair_obj = FAIR()
     species = ["CO2", "CH4", "N2O"]
@@ -55,6 +56,4 @@ def test_ohc_restart():
     initialise(ffut.ocean_heat_content_change, ohc_out)
     ffut.run()
     # check if first OHC timebound in restart is same as last in historical
-    np.testing.assert_allclose(
-        ffut.ocean_heat_content_change[0, ...], ohc_out
-    )
+    np.testing.assert_allclose(ffut.ocean_heat_content_change[0, ...], ohc_out)