Merge branch 'master' of github.com:fabm-model/fabm into scalar_stand…

…ard_default

src/c/fabm_c.F90

@@ -476,7 +476,7 @@ function get_variable(pmodel, category, index) bind(c) result(pvariable)
       case (HORIZONTAL_DIAGNOSTIC_VARIABLE)
          variable => model%p%horizontal_diagnostic_variables(index)%original
       case (CONSERVED_QUANTITY)
-         variable => model%p%conserved_quantities(index)%target
+         variable => model%p%conserved_quantities(index)%target_hz
       case (INTERIOR_DEPENDENCY, HORIZONTAL_DEPENDENCY, SCALAR_DEPENDENCY)
          select case (category)
          case (INTERIOR_DEPENDENCY);   domain = domain_interior

src/fabm_expressions.F90

@@ -184,8 +184,6 @@ function interior_temporal_mean(input, period, resolution, missing_value) result
       real(rk), optional,         intent(in)            :: missing_value
       type (type_interior_temporal_mean)                :: expression
 
-      character(len=attribute_length) :: prefix, postfix
-
       if (.not. associated(input%link)) call fatal_error('fabm_expressions::interior_temporal_mean', &
          'Input variable has not been registered yet.')
 
@@ -203,8 +201,6 @@ function horizontal_temporal_mean(input, period, resolution, missing_value) resu
       real(rk), optional,                    intent(in)            :: missing_value
       type (type_horizontal_temporal_mean)                         :: expression
 
-      character(len=attribute_length) :: prefix, postfix
-
       if (.not. associated(input%link)) call fatal_error('fabm_expressions::horizontal_temporal_mean', &
          'Input variable has not been registered yet.')
 
@@ -238,7 +234,6 @@ subroutine interior_temporal_mean_update(self, time, value _POSTARG_LOCATION_RAN
       real(rke) _ATTRIBUTES_GLOBAL_,       intent(in)    :: value
       _DECLARE_ARGUMENTS_LOCATION_RANGE_
 
-      integer  :: ibin
       real(rke) :: dt, w, dt_bin, scale
       _DECLARE_LOCATION_
 
@@ -344,7 +339,6 @@ subroutine horizontal_temporal_mean_update(self, time, value _POSTARG_HORIZONTAL
       real(rke) _ATTRIBUTES_GLOBAL_HORIZONTAL_, intent(in)    :: value
       _DECLARE_ARGUMENTS_HORIZONTAL_LOCATION_RANGE_
 
-      integer  :: ibin
       real(rke) :: dt, w, dt_bin, scale
       _DECLARE_HORIZONTAL_LOCATION_
 

src/pyfabm/__init__.py

@@ -592,7 +592,7 @@ def __getitem__(self, key: str) -> Union[float, int, bool]:
         raise KeyError
 
 
-class Variable(object):
+class Variable:
     def __init__(
         self,
         model: "Model",
@@ -688,6 +688,10 @@ def value(self, value: npt.ArrayLike):
             self.link(np.empty(self._shape, dtype=self.model.fabm.numpy_dtype))
         self._data[...] = value
 
+    @property
+    def shape(self) -> Tuple[int]:
+        return self._shape
+
     def link(self, data: np.ndarray):
         assert data.shape == self._shape, (
             f"{self.name}: shape of provided array {data.shape}"
@@ -855,6 +859,19 @@ def reset(self):
         self.model._update_configuration(settings)
 
 
+class ConservedQuantity(Variable):
+    def __init__(
+        self, model: "Model", name: str, units: str, variable_pointer: ctypes.c_void_p
+    ):
+        super().__init__(model, name, units, name.replace("_", " "))
+        self._pvariable = variable_pointer
+
+    @property
+    def missing_value(self) -> float:
+        """Value that indicates missing data, for instance, on land."""
+        return self.model.fabm.variable_get_missing_value(self._pvariable)
+
+
 class StandardVariable:
     def __init__(self, model: "Model", pointer: ctypes.c_void_p):
         self.model = model
@@ -1055,7 +1072,7 @@ def __init__(
         self.bottom_state_variables: NamedObjectList[StateVariable] = NamedObjectList()
         self.interior_diagnostic_variables: NamedObjectList[DiagnosticVariable] = NamedObjectList()
         self.horizontal_diagnostic_variables: NamedObjectList[DiagnosticVariable] = NamedObjectList()
-        self.conserved_quantities: NamedObjectList[Variable] = NamedObjectList()
+        self.conserved_quantities: NamedObjectList[ConservedQuantity] = NamedObjectList()
         self.parameters: NamedObjectList[Parameter] = NamedObjectList()
         self.interior_dependencies: NamedObjectList[Dependency] = NamedObjectList()
         self.horizontal_dependencies: NamedObjectList[Dependency] = NamedObjectList()
@@ -1364,13 +1381,13 @@ def _update_configuration(self, settings: Optional[Tuple] = None):
                 strlong_name,
                 strpath,
             )
+            ptr = self.fabm.get_variable(self.pmodel, CONSERVED_QUANTITY, i + 1)
             self.conserved_quantities._data.append(
-                Variable(
+                ConservedQuantity(
                     self,
                     strname.value.decode("ascii"),
                     strunits.value.decode("ascii"),
-                    strlong_name.value.decode("ascii"),
-                    strpath.value.decode("ascii"),
+                    ptr,
                 )
             )
         for i in range(nparameters.value):
@@ -1427,7 +1444,9 @@ def _update_configuration(self, settings: Optional[Tuple] = None):
 
         self.itime = -1.0
 
-    def getRates(self, t: Optional[float] = None, surface: bool = True, bottom: bool = True):
+    def getRates(
+        self, t: Optional[float] = None, surface: bool = True, bottom: bool = True
+    ):
         """Returns the local rate of change in state variables,
         given the current state and environment.
         """
@@ -1501,6 +1520,9 @@ def get_vertical_movement(self, out: Optional[np.ndarray] = None) -> np.ndarray:
         return out
 
     def get_conserved_quantities(self, out: Optional[np.ndarray] = None) -> np.ndarray:
+        """Return depth-integrated values across the horizontal domain for all
+        conserved quantities
+        """
         assert (
             self._cell_thickness is not None
         ), "You must assign model.cell_thickness to use get_conserved_quantities"

util/developers/run_all_testcases.py

@@ -389,9 +389,13 @@ def test_pyfabm(args, testcases: Mapping[str, str]):
                 if val != 0.0:
                     bad[var.name] = val
             assert False, f"Variability among 1D results: {r1d} (range: {bad})"
-        assert (
-            r1d[:, 0] == r0d
-        ).all(), f"Mismatch between 0D and 1D results: {r0d} vs {r1d[:, 0]}. Difference: {r1d[:, 0] - r0d}"
+        if (r1d[:, 0] != r0d).any():
+            diff = r1d[:, 0] - r0d
+            bad = {}
+            for var, val in zip(m1d.state_variables, diff):
+                if val != 0.0:
+                    bad[var.name] = val
+            assert False, f"Mismatch between 0D and 1D results: {r0d} vs {r1d[:, 0]}. Difference: {diff}. {bad}"
         print("SUCCESS")
     pyfabm_libs = ", ".join([f"{n}={l._name}" for n, l in pyfabm.name2lib.items()])
     print(f"pyfabm {pyfabm.__version__} loaded from {pyfabm.__file__} ({pyfabm_libs})")