Add MOC stream function analysis capability

components/omega/configs/Default.yml

@@ -207,11 +207,37 @@ Omega:
   Analysis:
     GlobalStats:
       Enable: false
-      Fields: ["NormalVelocity", "PseudoThickness", "Temperature", "Salinity"]
-      SpatialStats: ["Max", "Min", "Mean", "StdDev"]
-      ReductionPeriod: ["1Day", "1Month"]
-      SnapshotPeriod: ["6Hours"]
-      Filename: global.stats.$Y
+      # Global Statistics (GlobalStats) analysis group
+      # Computes spatial reduction statistics (Mean, Min, Max, StdDev)
+      # for a set of ocean fields. Supports temporal reduction (time-averaged
+      # output over a window) and instantaneous snapshots (discrete sampling).
+      Fields: [NormalVelocity, PseudoThickness, Temperature, Salinity]
+                                # List of field names to compute statistics for
+      SpatialStats: [Max, Min, Mean, StdDev]
+                                # Spatial statistics to compute (one per field)
+      ReductionPeriod: [1Day, 1Month]
+                                # Temporal reduction periods (time-averaged stats)
+      SnapshotPeriod: [6Hours]
+                                # Instantaneous snapshot output intervals
+      Filename: global.stats.$Y.nc
       Stream:
         FileFreq: 1
         FileFreqUnits: years
+    MOC:
+      Enable: false
+      # Meridional Overturning Circulation (MOC) streamfunction analysis group
+      # Computes MOC as a function of latitude and depth for regions,
+      # and as a function of depth for transects
+      NumBins: 180              # Number of latitude bins (default: 180, ~1 degree)
+      MinLat: -90.0             # Minimum latitude in degrees (default: -90.0)
+      MaxLat: 90.0              # Maximum latitude in degrees (default: 90.0)
+      Regions: [Global]         # List of region names for regional MOC
+                                # NOTE: Region masks not yet implemented
+      Transects: []             # List of transect names for transect-based MOC
+                                # NOTE: Transect masks not yet implemented
+      ReductionPeriod: [1Month]       # Temporal reduction periods
+      SnapshotPeriod: []  # Instantaneous output periods
+      Filename: moc.$Y.nc
+      Stream:
+        FileFreq: 1
+        FileFreqUnits: months

components/omega/doc/devGuide/Analysis.md

@@ -200,6 +200,15 @@ field(s).
 | `SpatialMean` | 1 | 1 | scalar (`Array1DReal`, dimension `Scalar`) | `_SpatialMean` | — | Global mean of the input field. |
 | `SpatialStdDev` | 2 (the field and its `_SpatialMean`) | 1 | scalar (`Array1DReal`, dimension `Scalar`) | `_SpatialStdDev` | — | Global standard deviation of the input field. Requires the field's `SpatialMean` as an upstream input, which is added to its input list automatically. |
 | `TimeMean` | 1 | 1 | same rank and dimensions as the input (`Real`) | `_TimeMean<period>` | `Period` (string, e.g. `"1Day"`) | Time average of the input field over a configurable period (e.g. `1Day`). Accumulates every time step and finalizes the mean when the period alarm rings. Output name embeds the period, e.g. `_TimeMean1Day`. |
+| `BinaryMultiply` | 2 | 1 | same as first input | `_BinaryMultiply(<field2>)` | — | Element-wise multiplication of two fields. Supports vertical expansion (1D field replicated across vertical layers when multiplied with 2D/3D field). |
+| `HorzMean` | 1 | 1 | 1D Real (`Array1DReal`) over the input's vertical dimension | `_HorzMean` | `WeightField` (string, optional, default `AreaCell`) | Area-weighted horizontal mean per vertical layer/interface. |
+| `BinnedAccumulator` | 2 (data field, bin index field) | 1 | replaces horizontal dimension with `NumBins` | `_BinnedAccumulator(<binfield>)` | `NumBins` (I4) | Accumulates field values into spatial bins. Uses MPI reduction for global totals. Automatically applies regional mask if attached to input Field. |
+| `CoordinateBinning` | 1 | 1 | 1D integer (`Array1DI4`) | `_BinIndex` | `NumBins` (I4), `MinBin` (Real, optional), `MaxBin` (Real, optional) | Assigns mesh entities to bins based on coordinate values (e.g., latitude). Computed once during initialization and cached. |
+| `ExtractRegion` | 1 | 1 | same as input | `_<RegionName>` | `MaskName` (string) | Applies regional mask to field by multiplying values with mask. Attaches mask to output Field for downstream operators. Supports mask intersection. |
+| `PrefixSum` | 1 | 1 | same as input | `_PrefixSum` | `Dimension` (I4), `Reverse` (bool, default false) | Cumulative summation along specified dimension. Forward (start→end) or reverse (end→start) scan. |
+| `PseudoToGeometric` | 1 | 1 | same as input | `_Geometric` | — | Converts pseudo-height coordinates to geometric coordinates using specific volume. Handles vertical grid staggering automatically. |
+| `ScalarMultiply` | 1 | 1 | same as input | `_ScalarMultiply(<value>)` | `Scalar` (string, parsed to Real) | Multiplies field by scalar constant. Can be used for unit conversions (e.g., ×1e-6 for Sverdrups). |
+| `TransectAccumulator` | 3 (data field, mask field, sign field) | 1 | 1D vertical profile | `_TransectAccumulator(<transect>)` | `TransectName` (string) | Accumulates transport across transect edges. Only supports 2D Real arrays for the data field. Uses MPI reduction for global totals. |
 
 ## Operator factory and type dispatch
 
@@ -220,6 +229,19 @@ are 1 through 3 (the 4D and 5D entries are present but commented out because
 Omega currently utilizes no arrays with rank $>$ 3, and each registered
 operator variant increases compile time and binary size).
 
+For operators that only support specific array ranks, the factory provides
+specialized registration methods using SFINAE (Substitution Failure Is Not An
+Error) to prevent invalid template instantiations at compile time:
+
+- `register1DVariants` — registers only 1D array variants (horizontal fields without vertical structure)
+- `register2DVariants` — registers only 2D array variants (horizontal × vertical)
+- `register2DRealVariants` — registers only 2D Real array variants (excludes integer types)
+
+These methods use compile-time type checking to filter out incompatible array
+types, reducing compile time and binary size while catching type errors early.
+The SFINAE helpers (`register1DVariantHelper`, etc.) expand over all array types
+but only register the matching variants, with no-op branches for mismatched types.
+
 At creation time, the orchestrator calls:
 ```c++
 auto Op = AnalysisOpFactory::createOp(OpType, UpstreamNames, Options);
@@ -250,6 +272,22 @@ spatial reduction; a `Time` prefix denotes a temporal reduction, with the
 period parsed from the first digit in the token) and the operator is created
 via the factory and wrapped in an `OperatorNode`.
 
+**Parenthesized arguments:** The parser supports operators with parenthesized
+arguments for passing field names, scalar values, or other parameters. The
+parser respects parenthesis depth when splitting on underscores, so
+`Field_BinaryMultiply(OtherField)_ScalarMultiply(1.0e-6)` correctly parses as
+three tokens: `Field`, `BinaryMultiply(OtherField)`, and `ScalarMultiply(1.0e-6)`.
+Different operators interpret parenthesized arguments differently:
+- `BinaryMultiply(FieldName)` — second input field name
+- `ScalarMultiply(Value)` — scalar value (parsed from string to Real)
+- `ExtractRegion(MaskName)` — mask field name passed via Config
+- `BinnedAccumulator(BinIndexField)` — bin index field name
+- `TransectAccumulator(TransectName)` — transect name for mask lookup
+
+An optional `Config` parameter can be passed to `parseChainAndBuildOps()` to
+provide operator-specific configuration that applies to all operators in the
+chain (e.g., binning parameters, integration direction).
+
 After all groups have registered their chains, the constructor resolves the
 graph edges with `buildOperatorDependencies()`, which matches each node's
 input Field names against the output Field names of every other node and
@@ -350,6 +388,61 @@ the discrete-sampling chain for instantaneous output. It then calls
 `createAnalysisGroupStreams()` to create the corresponding output streams. At
 least one of `ReductionPeriod` or `SnapshotPeriod` must be present.
 
+### MOC
+
+`MOC` computes the Meridional Overturning Circulation (MOC) streamfunction
+using two methods: latitude-binned regional MOC and transect-based MOC.
+Configuration specifies binning parameters, region/transect names, and output
+frequencies:
+
+```yaml
+MOC:
+  Enable: true
+  NumBins: 180                  # Number of latitude bins
+  MinLat: -90.0                 # Minimum latitude in degrees
+  MaxLat: 90.0                  # Maximum latitude in degrees
+  Regions: [Global, Atlantic]   # Regional MOC (latitude × depth)
+  Transects: [Drake, Atlantic26N]  # Transect MOC (depth only)
+  ReductionPeriod: [1Month]
+  SnapshotPeriod: [1Day]
+```
+
+**Latitude-binned MOC** computes the streamfunction as a function of latitude
+and depth for each specified region. The constructor builds a complex operator
+chain for each region:
+1. `LatCell_CoordinateBinning` — assigns cells to latitude bins (initialization only, shared across regions)
+2. `VerticalPseudoVelocity_PseudoToGeometric` — converts to geometric coordinates
+3. `VerticalVelocity_BinaryMultiply(AreaCell)` — computes vertical flux
+4. `[Optional] VerticalFlux_ExtractRegion(RegionMask)` — applies regional mask (for non-global regions)
+5. `VerticalFlux_BinnedAccumulator(LatCell_BinIndex)` — accumulates flux into latitude bins
+6. `BinnedFlux_PrefixSum` — horizontal integration (south to north)
+7. `MOC_ScalarMultiply(1.0e-6)` — converts to Sverdrups (Sv)
+
+**Transect-based MOC** computes the streamfunction across specific transects
+as a function of depth only. The constructor builds a chain for each transect:
+1. `PseudoThickness_PseudoToGeometric` — converts to geometric layer thickness
+2. `LayerThickness_BinaryMultiply(NormalVelocity)` — thickness × velocity
+3. `EdgeTransport_BinaryMultiply(DvEdge)` — transport × edge width
+4. `TransportField_TransectAccumulator(TransectName)` — accumulates across transect edges
+5. `TransectTransport_PrefixSum` — vertical integration (bottom to top, reverse)
+6. `TransectMOC_ScalarMultiply(1.0e-6)` — converts to Sverdrups
+
+For each chain, the constructor stores an operator configuration (`Config`)
+containing parameters like `NumBins`, `Dimension`, `Reverse`, and `TransectName`
+in a `ChainConfigs` vector. These configs are passed to `buildTemporalChains()`,
+which appends temporal operators and applies custom `IOName` metadata to
+produce user-friendly output variable names (e.g., `MOC_streamfunction_Global`
+instead of the full operator chain string). At least one of `Regions` or
+`Transects` must be specified.
+
+```{note}
+Regional and transect masks are not yet implemented in the Omega infrastructure.
+The `Regions` configuration currently only supports `[Global]`, and `Transects`
+is a placeholder. The operator infrastructure (ExtractRegionOp,
+TransectAccumulatorOp) is fully implemented and ready for use once mask fields
+become available.
+```
+
 ## Extensibility
 
 The module is designed so that new operators and new groups plug in through

components/omega/doc/devGuide/Field.md

@@ -263,6 +263,37 @@ and before exiting, all fields should be removed using:
    Field::clear();
 ```
 
+Fields can optionally have a regional mask attached to support spatial subsetting
+in analysis operations. A regional mask is a 1D integer array (`Array1DI4`) over
+the horizontal dimension (cells, edges, or vertices) where mask values indicate
+inclusion (1) or exclusion (0) of each location.
+
+Regional masks are set, queried, and retrieved using:
+```c++
+   // Set a regional mask for a field (shallow copy - shares data with source)
+   MyField->setRegionalMask(MaskArray);
+
+   // Check if a field has a regional mask attached
+   bool HasMask = MyField->hasRegionalMask();
+
+   // Retrieve the regional mask (returns Array1DI4)
+   Array1DI4 Mask = MyField->getRegionalMask();
+```
+
+Key implementation details:
+- The mask is stored as a shallow copy (Kokkos view), so the original mask array and the field's mask share the same data
+- `hasRegionalMask()` returns true if `setRegionalMask()` was called, even if the pointer is null (used for lifetime bug detection)
+- `getRegionalMask()` validates pointer integrity and aborts if mask was set but deallocated
+- Regional masks are automatically propagated through analysis operator chains:
+  - Binary operators (e.g., BinaryMultiplyOp) propagate mask from first input
+  - ExtractRegionOp attaches masks and can compute mask intersections
+  - Spatial reduction operators (e.g., BinnedAccumulatorOp) apply masks automatically
+
+Regional masks are primarily used by the Analysis framework for regional statistics
+and diagnostics (e.g., Atlantic MOC, basin-averaged temperatures). The mask mechanism
+enables efficient spatial subsetting without modifying field data or creating separate
+regional field copies.
+
 As mentioned above, Fields can be assigned to groups to provide an easy way
 to reference fields that commonly appear together, especially when listing
 contents of fields in IO files. Internally, a field group is implemented as

components/omega/doc/devGuide/IOStreams.md

@@ -103,6 +103,28 @@ The Metadata corresponding to ForcingTime will then be read from the file
 and inserted as the Metadata value. If no metadata is to be read from the
 file, then an empty ReqMetadata variable can be passed.
 
+### IOName Metadata for User-Friendly Output Variable Names
+
+Fields written to netCDF files can optionally specify a custom output variable
+name using the `IOName` metadata field. By default, IOStream uses the field's
+internal name (which may be a long operator chain string like
+`VerticalPseudoVelocity_PseudoToGeometric_BinaryMultiply_...`) as the netCDF
+variable name. Setting `IOName` metadata allows specifying a more readable name
+for end users:
+
+```c++
+   MyField->addMetadata("IOName", "MOC_streamfunction_Global");
+```
+
+During output, IOStream checks for the presence of `IOName` metadata and uses it
+as the netCDF variable name if found, falling back to the field's internal name
+otherwise. The `IOName` metadata is automatically filtered out of CF attribute
+writing since it is an internal control key, not a CF-compliant attribute.
+
+This feature is particularly useful for analysis outputs where operator chain
+names are long and technical, but output files should have concise, descriptive
+variable names.
+
 As described in the [User Guide](#omega-user-iostreams), all streams are
 defined in the input configuration file and most other IOStream functions
 are associated either with that initialization or to support the read/write

components/omega/doc/userGuide/Analysis.md

@@ -147,6 +147,80 @@ Automatically created and named:
  - Instantaneous output: `GlobalStats_FreqInstants` (e.g.,
  `GlobalStats_6HourInstants`)
 
+### MOC
+
+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.
+
+**Example:**
+
+```yaml
+Omega:
+  Analysis:
+    MOC:
+      Enable: true
+      NumBins: 180
+      MinLat: -90.0
+      MaxLat: 90.0
+      Regions: [Global]
+      Transects: []
+      ReductionPeriod: [1Month]
+      SnapshotPeriod: [1Day]
+      Filename: moc.$Y-$M
+      Stream:
+        FileFreq: 1
+        FileFreqUnits: months
+```
+
+**Group-Specific Parameters:**
+
+- **NumBins:** Optional integer specifying the number of latitude bins for
+  regional MOC (default: 180, approximately 1-degree resolution).
+
+- **MinLat:** Optional minimum latitude in degrees for binning (default: -90.0).
+
+- **MaxLat:** Optional maximum latitude in degrees for binning (default: 90.0).
+
+- **Regions:** Optional list of region names for regional MOC computation
+  (default: `[Global]`). Regional MOC computes the streamfunction as a function
+  of latitude and depth. Currently only `Global` is supported; region masks are
+  not yet implemented.
+
+- **Transects:** Optional list of transect names for transect-based MOC
+  computation (default: `[]`). Transect MOC computes the streamfunction as a
+  function of depth only by accumulating transport across specified transect
+  edges. Transect masks are not yet implemented.
+
+- **ReductionPeriod:** Optional list of time periods for temporal reduction.
+  Each period must divide evenly into the restart interval. At least one of
+  `ReductionPeriod` or `SnapshotPeriod` must be specified.
+
+- **SnapshotPeriod:** Optional list of intervals for instantaneous output. At
+  least one of `ReductionPeriod` or `SnapshotPeriod` must be specified.
+
+**Output fields:**
+
+- Regional MOC: `MOC_streamfunction_RegionName` (2D: latitude bins × depth)
+- Transect MOC: `MOC_streamfunction_transect_TransectName` (1D: depth only)
+
+For temporal reduction, `_TimeMeanPeriod` is appended (e.g.,
+`MOC_streamfunction_Global_TimeMean1Month`).
+
+**Output streams:**
+Automatically created and named:
+ - Time reduction: `MOC_FreqTimeStats` (e.g., `MOC_1MonthTimeStats`)
+ - Instantaneous output: `MOC_FreqInstants` (e.g., `MOC_1DayInstants`)
+
+```{note}
+Regional and transect masks are not yet implemented. The `Regions` parameter
+currently only supports `[Global]`, and `Transects` is a placeholder for future
+functionality. The underlying operator infrastructure is complete and ready for
+use once the mask fields are implemented in Omega.
+```
+
 ## Usage Notes
 
 ### Temporal Reduction Period Constraint