Define the mass balance fluxes everywhere

docs/dev/generating-test-files.md

@@ -70,7 +70,7 @@ ismip7-generate-test-files --grid GrIS_16000m --scenario ctrl \
   - Flux (FL) variables: timestamp is Jul 1 of year N (mid-year), with `time_bounds` = [Jan 1 of N, Jan 1 of N+1].
   - `x,y,z,t` variables (e.g. `litemp`): ST snapshots at a sparse set of nominal years. For `historical` and `ocx`: first year of run, 1900 (if in range), last year of run. For projection scenarios: 2100, 2200, 2300 (each if within the simulation year range). The filename year range reflects the full simulation period, not the first/last snapshot year. The set is `CENTURY_SNAPSHOT_YEARS` in `isschecker.checker`, imported here so the generator and the checker cannot disagree; note that no snapshot is written at 2000, which the checker accepts but does not require (see {doc}`../user/time-encoding`).
 - An analytic ice sheet, so that the files agree with one another. Most variables are still drawn at random within the min/max range the data request allows them, but seven cannot be: `sftgif`, `sftgrf`, `sftflf`, `lithk`, `topg`, `base` and `orog` are computed from a radially symmetric dome on a bowl-shaped bed, because the checker compares them against each other and a random thickness does not agree with a random mask about where the ice is. The identities hold exactly rather than to a tolerance: grounded plus floating fraction is the ice fraction, thickness is positive exactly where the ice fraction is, surface elevation is base plus thickness, and the ice base rests on the bed where the ice is grounded and above it where it floats. The geometry is deliberately analytic rather than taken from BedMachine or Bedmap: a synthetic file's job is to be predictable, and depending on one dataset would invite depending on datasets for everything else.
-- Missing values follow each variable's `fill_policy`. A variable the request defines only where there is ice is written as the fill value wherever the ice fraction is zero, and likewise for grounded-only and floating-only variables; the ones defined only inside the computational domain are missing outside it; and the rest are defined everywhere.
+- Missing values follow each variable's `fill_policy`. A variable the request defines only where there is ice is written as the fill value wherever the ice fraction is zero, and likewise for grounded-only and floating-only variables; the ones defined only inside the computational domain are missing outside it; and the rest are defined everywhere. The basal mass balance fluxes, `libmassbfgr` and `libmassbffl`, are also set to 0 where their grounded or floating fraction is zero, as the checker requires.
 - Single precision (`float32`) for all variables and time.
 - `_FillValue` and `missing_value` set to NetCDF4 default `f4` fill value.
 - `time` is an unlimited (record) dimension.

docs/user/checks.md

@@ -71,8 +71,9 @@ required set of snapshots. Both are described in {doc}`time-encoding`.
 
 Gridded variables only. Each file is compared against the files beside it: a
 variable is missing exactly where its ice mask says there is no ice, the
-variables of the computational domain cover the ice, the grounded and
-floating fractions sum to the ice fraction, thickness agrees with the ice
+variables of the computational domain cover the ice, the basal mass balance
+fluxes are zero where there is none of their ice, the grounded and floating
+fractions sum to the ice fraction, thickness agrees with the ice
 mask, and surface elevation, ice base and bed agree with each other. See
 [Checks that compare files](errors-and-warnings.md#checks-that-compare-files).
 

docs/user/errors-and-warnings.md

@@ -66,12 +66,17 @@ every variable in its fill_policy column; the table is under
 
 Two points modelers ask about. **Ice thickness is zero where there is no
 ice**, not missing, including outside your computational domain, and so are
-the three masks and the calving, grounding-line and ice-front fluxes. **The
+the three masks. **All mass fluxes are defined everywhere**: acabf,
+libmassbfgr, libmassbffl and the calving, grounding-line and ice-front
+fluxes hold a number in every cell, never a fill value. Where there is no
+ice, that number is usually 0, but acabf may be nonzero there. **The
 masks are not restricted to 0 and 1**: any fraction in [0, 1] is accepted,
 because conservative interpolation from your native grid legitimately
 produces intermediate values.
 
-A `forbidden` variable holding any fill value is an error. The other
+A `forbidden` variable holding any fill value is an error. For acabf,
+libmassbfgr and libmassbffl it is a warning for now ([issue #39]); see
+[Severities that will change](#severities-that-will-change). The other
 policies say where a field sits relative to the ice masks, which takes more
 than one file to check; see
 [Checks that compare files](#checks-that-compare-files).
@@ -87,6 +92,7 @@ A variable whose fill_policy cell is blank is unconstrained, and nothing
 about its missing values is checked. No shipped row is blank today.
 
 [issue #23]: https://github.com/ismip/ISM_SimulationChecker/issues/23
+[issue #39]: https://github.com/ismip/ISM_SimulationChecker/issues/39
 
 ## Checks that compare files
 
@@ -106,6 +112,9 @@ the name but the variable, and compares the two.
   domain is. That is a warning, because a field taken from a forcing or
   reference dataset may legitimately cover more of the grid than the ice
   model does.
+- libmassbfgr is 0 where sftgrf is 0, and libmassbffl is 0 where sftflf is
+  0. Where there is some grounded or floating ice, any value is accepted,
+  including 0. For now this is a warning (see below).
 - sftgrf + sftflf equals sftgif; lithk is greater than zero exactly where
   sftgif is; orog equals base + lithk; and the ice base rests on the bed
   where sftgrf is 1 and lies above it where sftflf is 1. None of these needs
@@ -121,17 +130,20 @@ there. The grounded and floating comparisons against the bed are the
 exception: they are made only in cells that are wholly one or the other,
 since in a half-and-half cell the mean ice base sits somewhere between.
 
-### One severity that will change
+### Severities that will change
 
-Two findings turn on where a model puts the ice margin: a variable being
-*missing where there is ice*, and thickness disagreeing with the ice mask. A
+Three findings turn on where a model puts the ice margin: a variable being
+*missing where there is ice*, thickness disagreeing with the ice mask, and a
+basal mass balance flux that is not 0 where its mask is 0. A
 conservatively interpolated mask puts fractions like 1e-6 in a ring along
 the edge, and a model that writes fill from its own native-grid mask will
 disagree in every one of those cells.
 
-One round of real submissions will settle which is right, so **those two
+One round of real submissions will settle which is right, so **those three
 findings are warnings for the first round and will become errors
-afterwards**. The definitions above are fixed, and a model can be written to
+afterwards**. A fill value in acabf, libmassbfgr or libmassbffl is also a
+warning for now: these fluxes were allowed missing values in the first
+round. The definitions above are fixed, and a model can be written to
 them today. The finding reports how much of itself is margin, "300 of them
 have sftgif below 0.01", so the decision can be made on evidence.
 

isschecker/checker.py

@@ -32,9 +32,11 @@
 #      an error whatever the variable is.
 #    - A variable whose 'fill_policy' is 'forbidden' holds no fill values at
 #      all: it is defined over the whole grid, and where there is no ice the
-#      value is 0 rather than missing.  The other policies -- outside_domain,
-#      no_ice, no_grounded_ice, no_floating_ice -- say where a field sits
-#      relative to the ice masks, which takes more than one file to check.
+#      value is 0 rather than missing.  That is an error, or a warning for a
+#      variable whose 'fill_severity' says so -- see _fill_severity.  The
+#      other policies -- outside_domain, no_ice, no_grounded_ice,
+#      no_floating_ice -- say where a field sits relative to the ice masks,
+#      which takes more than one file to check.
 #    - All values lie within the allowed min/max range for the relevant region.
 #      A few cells outside it are a warning; RANGE_ERROR_FRACTION or more of
 #      the values are an error, unless the variable's 'range_severity' says
@@ -52,6 +54,8 @@
 #    - An outside_domain variable is defined wherever there is ice (error), and
 #      the outside_domain variables agree about where the domain is (warning:
 #      a forcing or reference dataset may legitimately cover more).
+#    - libmassbfgr is 0 where sftgrf is 0, and libmassbffl is 0 where sftflf
+#      is 0, at the severity 'margin_severity' names.
 #    - sftgrf + sftflf equals sftgif; lithk is greater than 0 exactly where
 #      sftgif is; orog equals base + lithk; and the ice base rests on the bed
 #      where sftgrf is 1 and lies above it where sftflf is 1.  No densities:
@@ -578,6 +582,20 @@ def _margin_severity(value) -> str:
     return "error"
 
 
+def _fill_severity(value) -> str:
+    """How to report a fill value in a variable that may not have any.
+
+    Issue #39 moved the surface and basal mass balance fluxes to `forbidden`
+    after the first round of submissions had been written to the old policies,
+    so for those variables a fill value is a warning until the later rounds.
+    As with the other severity columns, anything the column does not say means
+    error.
+    """
+    if value is not None and pd.notna(value) and str(value).strip().lower() == "warning":
+        return "warning"
+    return "error"
+
+
 def _fill_policy(value) -> str | None:
     """Where a variable is defined, from the fill_policy column of the request.
 
@@ -644,6 +662,7 @@ def _load_criteria(variable_list: str):
             "range_severity": _range_severity(row.get("range_severity")),
             "fill_policy": _fill_policy(row.get("fill_policy")),
             "margin_severity": _margin_severity(row.get("margin_severity")),
+            "fill_severity": _fill_severity(row.get("fill_severity")),
         }
         for col in df.columns:
             lc = str(col).lower()
@@ -1824,7 +1843,12 @@ def _check_missing_values(reporter, ivar, criteria, is_fill, is_nonfinite, fill)
     # whole grid, so a hole in it is a hole in the archive.
     n_fill = int(is_fill.sum())
     if n_fill and criteria.get("fill_policy") == "forbidden":
-        reporter.error(
+        report = (
+            reporter.warning
+            if criteria.get("fill_severity") == "warning"
+            else reporter.error
+        )
+        report(
             f"variable '{ivar}' holds a fill value in"
             f" {_count_phrase(n_fill, total)}. The data request does not permit"
             f" missing values in this variable: where there is no ice the value"
@@ -2133,6 +2157,13 @@ def _open_companion(reporter, source_path, file_name, variable, years,
     "no_floating_ice": "sftflf",
 }
 
+# The basal mass balance fluxes are defined everywhere, but are zero wherever
+# there is none of the ice they are the flux beneath (issue #39).
+ZERO_WITHOUT_MASK = {
+    "libmassbfgr": "sftgrf",
+    "libmassbffl": "sftflf",
+}
+
 # Below this fraction a cell is the partly glaciated margin rather than the
 # interior.  Used only to describe a finding, never to decide one: the rule is
 # exact, and this says how much of what it found is margin.
@@ -2214,6 +2245,45 @@ def _check_ice_extent(reporter, ds, ivar, criteria, companion) -> None:
         )
 
 
+def _check_zero_without_ice(reporter, ds, ivar, criteria, companion) -> None:
+    """A basal flux is 0 wherever there is none of the ice it lies beneath.
+
+    Only that direction: a cell with some of that ice may still have no melt.
+    Fill cells have been reported already and are left out here.
+    """
+    values, missing, steps = _field_and_steps(ds, ivar)
+    nonzero_without_ice = 0
+    total = 0
+
+    for step in range(steps):
+        here = _step_of(values, step)
+        fraction, unknown = _companion_slice(companion, step)
+        if here.shape != fraction.shape:
+            reporter.note(
+                f"Not checked against '{companion.variable}': its grid differs"
+                f" from this file's."
+            )
+            return
+        no_ice = np.logical_and(
+            np.logical_and(~unknown, ~_step_of(missing, step)), fraction == 0.0
+        )
+        nonzero_without_ice += int(np.logical_and(no_ice, here != 0.0).sum())
+        total += here.size
+
+    if nonzero_without_ice:
+        report = (
+            reporter.warning
+            if criteria.get("margin_severity") == "warning"
+            else reporter.error
+        )
+        report(
+            f"variable '{ivar}' is not 0 in"
+            f" {_count_phrase(nonzero_without_ice, total)} where"
+            f" '{companion.variable}' is 0. The flux must be 0 where there is"
+            f" no ice for it to be beneath."
+        )
+
+
 def _companion_slice(companion, step):
     """The companion's values for one time step, and where it says nothing.
 
@@ -2282,8 +2352,8 @@ def _check_footprint_matches(reporter, ds, ivar, companion) -> None:
     They are missing in exactly one place -- outside it -- so their missing
     masks ought to be identical.  A warning rather than an error because some
     of them come from forcing or reference datasets whose footprint is
-    legitimately wider than the ice model's: acabf over ice-free ground,
-    hfgeoubed and refgeoid over the whole grid.
+    legitimately wider than the ice model's: hfgeoubed and refgeoid over the
+    whole grid.
     """
     is_fill, is_nonfinite = _missing_masks(ds, ivar)
     missing = np.logical_or(is_fill, is_nonfinite)
@@ -2555,7 +2625,7 @@ def _check_consistency(reporter, ds, ivar, criteria, source_path, file_name,
     # defined -- and some by name, because a particular pair or triple of them
     # has to agree about the ice sheet.  A variable can be both: orog is
     # defined throughout the domain and is also base plus thickness.
-    by_name = ivar in ("sftgif", "lithk", "orog", "base")
+    by_name = ivar in ("sftgif", "lithk", "orog", "base", *ZERO_WITHOUT_MASK)
     if ivar not in ds or (
         mask_name is None and policy != "outside_domain" and not by_name
     ):
@@ -2615,6 +2685,11 @@ def against(names, check):
             ["topg", "sftgrf", "sftflf"],
             lambda b, g, f: _check_base_against_bed(reporter, ds, ivar, b, g, f),
         )
+    elif ivar in ZERO_WITHOUT_MASK:
+        against(
+            [ZERO_WITHOUT_MASK[ivar]],
+            lambda c: _check_zero_without_ice(reporter, ds, ivar, criteria, c),
+        )
 
     if reporter.total_errors + reporter.total_warnings == findings_before:
         reporter.ok("Consistent with the files beside it: OK")