Give analysis output fields the units, standard name and cell methods of the field they reduce

components/omega/doc/devGuide/Analysis.md

@@ -144,7 +144,9 @@ A concrete derived operator is responsible for three things:
    record the input Field names, construct the output Field name(s), allocate
    the output data array, and register the output [Field](#omega-dev-field)
    with the appropriate metadata. Output Field names follow the convention
-   `<input>_<OperatorSuffix>` (for example, `Temperature_SpatialMax`).
+   `<input>_<OperatorSuffix>` (for example, `Temperature_SpatialMax`). The
+   output's CF metadata is derived from the input's, as described under
+   [Output metadata](#omega-dev-analysis-output-metadata).
 2. **Initialization**: `initialize()` is called by the orchestrator after all
    Fields exist. The base implementation stores the mesh, vertical
    coordinate, and communicator; operators may override to perform additional
@@ -245,6 +247,69 @@ normalization. The sums accumulate each owned entity's terms in a fixed
 order and reduce them with `globalSum`, so they are as reproducible as
 Omega's other global sums.
 
+(omega-dev-analysis-output-metadata)=
+### Output metadata
+
+An output Field inherits its CF metadata from the input it is derived from.
+The base class provides the helpers that do this:
+
+- `inheritMetadata(InputName, CellMethod)` returns an `InheritedMetadata`
+  struct holding the input's `units`, `standard_name` and `cell_methods`,
+  each empty if the input has none, with `CellMethod` appended to the cell
+  methods. CF requires `cell_methods` to list the reductions in the order
+  they were applied, so a time mean of a spatial mean carries
+  `area: depth: mean time: mean`.
+- `spatialCellMethod(InputName, Method)` returns the cell method for a
+  reduction over all owned mesh entities and layers of the input:
+  `area: <Method>` for a horizontal (1D) field and `area: depth: <Method>`
+  for a field with a vertical dimension. `Method` is a CF cell method name
+  (`mean`, `minimum`, `maximum`, `standard_deviation`, `sum`).
+- `createOutputField(...)` calls `Field::create()` with the inherited units
+  and standard name and adds the `cell_methods` attribute when it is not
+  empty.
+- `getFieldAttribute(FieldName, AttName)` returns any string attribute of a
+  Field, or an empty string if the Field has none.
+
+The built-in reductions (mean, minimum, maximum, standard deviation, time
+mean) all have the units and standard name of the field they reduce, so
+they pass `inheritMetadata()` straight to `createOutputField()`:
+```c++
+auto CellMethod  = spatialCellMethod(InputNames[0], "mean");
+auto Meta        = inheritMetadata(InputNames[0], CellMethod);
+auto OutputField = createOutputField(OutputNames[0],
+                                     "Spatial mean of " + InputNames[0],
+                                     Meta, ValidMin, ValidMax, NDims, DimNames);
+```
+An operator whose result has different units derives them with
+[CFUnits](#omega-dev-cfunits) before creating the output, and clears the
+standard name if the result is no longer the same physical quantity. For
+example, an area-weighted sum multiplies the input's units by those of the
+cell area:
+```c++
+auto Meta  = inheritMetadata(InputNames[0],
+                             spatialCellMethod(InputNames[0], "sum"));
+Meta.Units = CFUnits::multiply(Meta.Units, "m2");
+Meta.StdName.clear();
+```
+The other operators derive their units the same way:
+
+- `BinaryMultiply` multiplies the units of its two inputs (a velocity in
+  `m s-1` times an area in `m2` is `m3 s-1`) and has no standard name or cell
+  methods, since the product is a new quantity.
+- `ScalarMultiply` scales the units by the reciprocal of its scalar with
+  `CFUnits::scale()`, so a flux in `m3 s-1` multiplied by `1e-6` is in
+  `1e+06 m3 s-1` (Sverdrups). A negative scalar also clears the standard
+  name, since it reverses the quantity.
+- `HorzMean` inherits the units and standard name and appends `area: mean`.
+- `BinnedAccumulator`, `TransectAccumulator` and `PrefixSum` sum their input,
+  and `ExtractRegion` and `PseudoToGeometric` multiply it by a mask or a
+  dimensionless factor, so all of them inherit its units. `PrefixSum` aborts
+  if its boundary condition field has different units from its input.
+- `CoordinateBinning` outputs a bin index, which is dimensionless (`1`).
+
+A field without units yields outputs without units, never with invented ones,
+and a malformed units string aborts the run when the operator is constructed.
+
 ## Operator factory and type dispatch
 
 Fields carry runtime metadata describing their scalar type, rank, and memory
@@ -493,9 +558,11 @@ To add a new operator:
    operator type name to the base constructor, set `InputNames`, build the
    output name as `<input>_<Suffix>`, set `OutputNames` and `InstanceName`,
    allocate the output data array, and register the output Field with
-   `Field::create()` plus `attachData()`. If the operator depends on another
-   operator's output (as `SpatialStdDev` depends on `SpatialMean`), append
-   that dependency to `InputNames`.
+   `createOutputField()` plus `attachData()`, deriving its units, standard
+   name and cell methods from the input as described under
+   [Output metadata](#omega-dev-analysis-output-metadata). If the operator
+   depends on another operator's output (as `SpatialStdDev` depends on
+   `SpatialMean`), append that dependency to `InputNames`.
 2. **Implement `compute()`** to retrieve the input array(s) from the Field
    registry, perform the transformation, write the output array(s), and update
    `LastComputed`/`FieldComputed`.

components/omega/doc/devGuide/CFUnits.md

@@ -0,0 +1,58 @@
+(omega-dev-cfunits)=
+
+## CFUnits
+
+Omega fields carry their units as a string attribute in the plain form used
+by the CF standard name table: unit symbols separated by single spaces, each
+with an optional signed integer exponent (`m s-1`, `kg m-3`, `W m-2`), and `1`
+for a dimensionless quantity. When a field is derived from others, its units
+must be derived too. The `CFUnits` class parses such a string, supports the
+product, quotient and integer power of units expressions, and formats the
+result back in the plain form.
+
+The string interfaces cover the common cases:
+```c++
+#include "CFUnits.h"
+
+std::string Flux   = CFUnits::multiply("m s-1", "m2");   // "m3 s-1"
+std::string Speed  = CFUnits::divide("m3 s-1", "m2");    // "m s-1"
+std::string Energy = CFUnits::power("m s-1", 2);         // "m2 s-2"
+```
+Exponents of the same symbol are merged (`m s-1` times `m` is `m2 s-1`) and
+a symbol whose exponent reaches zero is dropped, so `kg m-3` times `m3` is
+`kg`. In the result, symbols with positive exponents come first, then those
+with negative exponents, each in the order they first appeared.
+
+A units string may also start with a positive number, the scale factor of a
+scaled unit as udunits writes it: `1e+06 m3 s-1` is a million cubic meters
+per second. Multiplying a field's values by a constant, for example to report
+a transport in units of 10^6 m3 s-1, scales its units by the reciprocal:
+```c++
+std::string Transport = CFUnits::scale("m3 s-1", 1.0 / 1.0e-6); // "1e+06 m3 s-1"
+```
+Scale factors multiply through products, quotients and powers, and are
+written with the fewest digits that reproduce them, so the rounding error in
+`1.0 / 1.0e-6` does not show. A number is only accepted before the unit
+symbols, and must be positive and finite.
+
+A field with no units attribute has *unknown* units, represented by the empty
+string. Unknown units propagate through every operation, so a field derived
+from one without units has no units either, rather than wrong ones.
+
+The parser also accepts the udunits caret and slash spellings (`m^2`, `m/s`),
+which some older field definitions still use, but always emits the plain
+form. Any other string is an error: the string interfaces abort, and the
+value interface reports it so a caller can decide:
+```c++
+CFUnits Units;
+Error Err = CFUnits::parse("m s-1", Units);  // Fail for a malformed string
+CFUnits Area;
+CFUnits::parse("m2", Area);
+std::string Flux = (Units * Area).str();     // "m3 s-1"
+```
+`CFUnits` also provides `operator/`, `pow(int)`, `scaled(double)`,
+`isUnknown()`, `isDimensionless()` and an order-insensitive `operator==`.
+
+A malformed units string is a programming error in a field definition, and
+aborting when the derived field is created stops the run at initialization
+rather than writing output with a wrong or missing attribute.

components/omega/doc/index.md

@@ -84,6 +84,7 @@ devGuide/Decomp
 devGuide/Dimension
 devGuide/Error
 devGuide/Field
+devGuide/CFUnits
 devGuide/IO
 devGuide/IOStreams
 devGuide/Halo

components/omega/doc/userGuide/Analysis.md

@@ -101,7 +101,11 @@ user-specified fields. The mean and standard deviation are weighted so that
 they do not depend on the mesh resolution: by the area of each cell, edge or
 vertex for a horizontal field, and by mass (area times pseudo-thickness) for
 a field with a vertical dimension, with an interface field weighted by half
-the mass of each adjacent layer. Only active layers take part.
+the mass of each adjacent layer. Only active layers take part. Each statistic
+is written with the `units` and `standard_name` of the field it reduces and a
+CF `cell_methods` attribute describing the reduction, for example
+`area: depth: mean` for the mean of a layered field and
+`area: depth: mean time: mean` for its time mean.
 
 **Example:**
 
@@ -157,7 +161,9 @@ Computes the Meridional Overturning Circulation (MOC) streamfunction using two
 methods: latitude-binned regional MOC and transect-based MOC. The MOC
 represents zonally integrated meridional mass transport as a function of
 latitude and depth (regional MOC) or depth alone (transect MOC). Output is in
-Sverdrups (Sv), where 1 Sv = 10⁶ m³/s.
+Sverdrups (Sv), where 1 Sv = 10⁶ m³/s, so its `units` attribute is
+`1e+06 m3 s-1`, the udunits spelling of a Sverdrup (udunits reads `Sv` as the
+sievert).
 
 **Example:**
 

components/omega/src/analysis/AnalysisOperator.cpp

@@ -87,6 +87,86 @@ const std::vector<std::string> AnalysisOperator::getOutputFieldNames() {
    return OutputNames;
 } // end getOutputFieldNames
 
+//------------------------------------------------------------------------------
+// Returns a string attribute of a Field, or an empty string if the Field has
+// no such attribute. Analysis output is derived from fields that may lack
+// units or a standard name, so absence is the empty string, not an error.
+std::string AnalysisOperator::getFieldAttribute(const std::string &FieldName,
+                                                const std::string &AttName) {
+
+   auto FieldPtr = Field::get(FieldName);
+   OMEGA_REQUIRE(FieldPtr, "AnalysisOperator: field {} not found", FieldName);
+
+   std::string Value;
+   if (FieldPtr->hasMetadata(AttName)) {
+      Error Err = FieldPtr->getMetadata(AttName, Value);
+      CHECK_ERROR_ABORT(Err,
+                        "AnalysisOperator: attribute {} of field {} is "
+                        "not a string",
+                        AttName, FieldName);
+   }
+   return Value;
+} // end getFieldAttribute
+
+//------------------------------------------------------------------------------
+// Collects the units, standard name and cell methods an output inherits from
+// its input. A new cell method is appended to the input's cell_methods so the
+// attribute lists the reductions in the order they were applied, as CF
+// requires (e.g. a time mean of a spatial mean is "area: mean time: mean").
+AnalysisOperator::InheritedMetadata
+AnalysisOperator::inheritMetadata(const std::string &InputName,
+                                  const std::string &CellMethod) {
+
+   InheritedMetadata Meta;
+   Meta.Units       = getFieldAttribute(InputName, "units");
+   Meta.StdName     = getFieldAttribute(InputName, "standard_name");
+   Meta.CellMethods = getFieldAttribute(InputName, "cell_methods");
+
+   if (!CellMethod.empty()) {
+      if (!Meta.CellMethods.empty())
+         Meta.CellMethods += " ";
+      Meta.CellMethods += CellMethod;
+   }
+   return Meta;
+} // end inheritMetadata
+
+//------------------------------------------------------------------------------
+// Returns the cell method for a reduction over all owned mesh entities and
+// layers of an input Field. The spatial operators treat a 1D field as
+// horizontal only and the last dimension of any higher-rank field as
+// vertical, so the method covers "area" alone or "area" and "depth". Both
+// are names CF allows in cell_methods without a matching coordinate.
+std::string AnalysisOperator::spatialCellMethod(const std::string &InputName,
+                                                const std::string &Method) {
+
+   auto FieldPtr = Field::get(InputName);
+   OMEGA_REQUIRE(FieldPtr, "AnalysisOperator: field {} not found", InputName);
+
+   if (FieldPtr->getNumDims() >= 2)
+      return "area: depth: " + Method;
+   return "area: " + Method;
+} // end spatialCellMethod
+
+//------------------------------------------------------------------------------
+// Creates an output Field carrying the inherited CF metadata. The cell
+// methods are stored as a std::string so the metadata can be retrieved and
+// written as a string attribute.
+std::shared_ptr<Field> AnalysisOperator::createOutputField(
+    const std::string &OutputName, const std::string &Description,
+    const InheritedMetadata &Meta, const std::any ValidMin,
+    const std::any ValidMax, const int NumDims,
+    const std::vector<std::string> &DimNames) {
+
+   auto OutputField =
+       Field::create(OutputName, Description, Meta.Units, Meta.StdName,
+                     ValidMin, ValidMax, NumDims, DimNames);
+
+   if (!Meta.CellMethods.empty())
+      OutputField->addMetadata("cell_methods", Meta.CellMethods);
+
+   return OutputField;
+} // end createOutputField
+
 //------------------------------------------------------------------------------
 // Checks whether the operator's output is valid for the given timestamp.
 // Returns true if the operator has been computed and the LastComputed

components/omega/src/analysis/AnalysisOperator.h

@@ -20,13 +20,21 @@
 /// data from the Field registry and writes results to operator-owned output
 /// arrays.
 ///
+/// An output Field inherits its CF metadata from the input it is derived
+/// from: the base class provides helpers that collect an input's units,
+/// standard name and cell methods, append the cell method describing the new
+/// reduction, and create the output Field with them. Operators whose result
+/// has different units (an area-weighted sum, a product) derive them with
+/// CFUnits before creating the output Field.
+///
 /// This file also provides helper functions (opParam, makeOpConfig) for
 /// constructing Config objects inline when instantiating operators,
 /// providing a uniform parameter-passing mechanism whether operators are
 /// created from user config or programmatically by bundled AnalysisGroups.
 ///
 //===----------------------------------------------------------------------===//
 
+#include "CFUnits.h"
 #include "Config.h"
 #include "DataTypes.h"
 #include "Dimension.h"
@@ -165,6 +173,58 @@ class AnalysisOperator {
    );
 
  protected:
+   /// CF metadata an output Field inherits from the input it is derived
+   /// from. Units and standard name are copied and are empty when the input
+   /// has none. CellMethods records every reduction applied so far, in
+   /// order, as the CF cell_methods attribute (e.g. "area: mean time: mean").
+   struct InheritedMetadata {
+      std::string Units;       ///< units of the input, empty if unknown
+      std::string StdName;     ///< CF standard name, empty if none
+      std::string CellMethods; ///< CF cell_methods, empty if none
+   };
+
+   /// Collects the CF metadata an output Field inherits from an input Field,
+   /// with the cell method of the new reduction (if any) appended to those
+   /// the input already records. The units and standard name are copied
+   /// unchanged; an operator that changes the units or invalidates the
+   /// standard name edits the result before creating its output Field.
+   static InheritedMetadata
+   inheritMetadata(const std::string &InputName,      ///< [in] input field name
+                   const std::string &CellMethod = "" ///< [in] new reduction
+   );
+
+   /// Returns the CF cell method for a spatial reduction of an input Field
+   /// over all of its owned mesh entities and layers: "area: <Method>" for a
+   /// horizontal field and "area: depth: <Method>" for a field with a
+   /// vertical dimension. Method is a CF cell method name (mean, minimum,
+   /// maximum, standard_deviation, sum).
+   static std::string
+   spatialCellMethod(const std::string &InputName, ///< [in] input field name
+                     const std::string &Method     ///< [in] CF cell method
+   );
+
+   /// Returns the value of a string attribute (units, standard_name,
+   /// cell_methods, long_name) of a Field, or an empty string if the Field
+   /// has no such attribute
+   static std::string
+   getFieldAttribute(const std::string &FieldName, ///< [in] field name
+                     const std::string &AttName    ///< [in] attribute name
+   );
+
+   /// Creates an output Field with the given description, valid range and
+   /// dimensions, carrying the inherited units, standard name and cell
+   /// methods. Empty units or standard name are passed through to
+   /// Field::create; an empty cell_methods is not added.
+   static std::shared_ptr<Field> createOutputField(
+       const std::string &OutputName,  ///< [in] output field name
+       const std::string &Description, ///< [in] long name of the output
+       const InheritedMetadata &Meta,  ///< [in] inherited CF metadata
+       const std::any ValidMin,        ///< [in] min valid value
+       const std::any ValidMax,        ///< [in] max valid value
+       const int NumDims,              ///< [in] number of dimensions
+       const std::vector<std::string> &DimNames ///< [in] dimension names
+   );
+
    const HorzMesh *Mesh;    ///< Horizontal mesh for spatial operations
    const VertCoord *VCoord; ///< Vertical coordinate for vertical ops
    MPI_Comm Comm;           ///< MPI communicator for parallel reductions

components/omega/src/analysis/analysisGroups/MOC.cpp

@@ -121,9 +121,9 @@ MOC::MOC(const std::string &GroupName, Config &AnalysisGroupOptions,
    buildTemporalChains(ChainStems, AnalysisGroupOptions, AnalysisManager,
                        ChainConfigs);
 
-   // Stamp Sverdrup units and a readable long_name on every streamfunction
-   // output field. Both the CF-compliant "units" and the legacy capitalized
-   // "Units" are set.
+   // Stamp a readable long_name on every streamfunction output field. The
+   // operators derive the units, 1e+06 m3 s-1 (Sverdrups) after the final
+   // ScalarMultiply(1e-6); "Sv" would be read by udunits as the sievert.
    for (size_t i = 0; i < ChainStems.size(); ++i) {
       const auto &StemStr    = ChainStems[i];
       const auto &RegionName = RegionList[i];
@@ -137,8 +137,6 @@ MOC::MOC(const std::string &GroupName, Config &AnalysisGroupOptions,
       for (const auto &SFName : SFFieldNames) {
          if (Field::exists(SFName)) {
             auto SFField = Field::get(SFName);
-            SFField->updateMetadata("units", std::string("Sv"));
-            SFField->updateMetadata("Units", std::string("Sv"));
             SFField->updateMetadata("long_name", LongName);
          }
       }
@@ -460,15 +458,14 @@ std::string MOC::buildDepthCoordChain(const std::string &RegionName,
    // intermediate fields shared with other chains).
    AnalysisManager->parseChainAndBuildOps(ChainStr, Cfg);
 
-   // Attach descriptive depth metadata for output (once per unique field).
+   // Attach a descriptive long_name for output (once per unique field). The
+   // HorzMean operator derives the units from GeomZInterface.
    if (Field::exists(DepthFieldName)) {
       auto DepthField = Field::get(DepthFieldName);
       DepthField->updateMetadata(
           "long_name",
           std::string("Area-weighted representative interface depth (") +
               RegionName + ")");
-      DepthField->updateMetadata("units", std::string("m"));
-      DepthField->updateMetadata("Units", std::string("m"));
    }
 
    return DepthFieldName;