Add MOC stream function analysis capability by brian-oneill · Pull Request #481 · E3SM-Project/Omega

brian-oneill · 2026-07-28T03:38:11Z

Overview

This PR introduces a complete MOC (Meridional Overturning Circulation) analysis capability to Omega, enabling computation of the MOC streamfunction using two complementary methods:

Latitude-binned regional MOC: Computes MOC as a function of latitude and depth for specified ocean regions
Transect-based MOC: Computes MOC across specific transects as a function of depth

The implementation adds 8 new analysis operators, enhanced analysis infrastructure with regional mask support, a new MOC analysis group, and comprehensive configuration templates.

Key Features

New Analysis Operators (8 total):

BinaryMultiplyOp: Element-wise field multiplication with vertical expansion support
BinnedAccumulatorOp: Accumulates field values into spatial bins (core MOC operator)
CoordinateBinningOp: Assigns mesh entities to bins based on coordinate values
ExtractRegionOp: Applies regional masks to fields
PrefixSumOp: Cumulative summation (integration) along specified dimension
PseudoToGeometricOp: Converts pseudo-height quantities to geometric coordinates
ScalarMultiplyOp: Multiplies field by scalar constant for unit conversion
TransectAccumulatorOp: Accumulates transport across transect edges

Infrastructure Enhancements:

Field Class Regional Mask Support: Fields can carry spatial mask information through operator chains with automatic propagation
IOStream IOName Metadata: Provides user-friendly variable names in netCDF output while maintaining internal operator chain naming
Enhanced AnalysisGroup Base Class: New setOutputIOName method and operator-specific configuration support

MOC Analysis Group:

Bundled analysis group for computing MOC streamfunction
Configurable latitude binning (number of bins, lat range)
Regional MOC computation for named ocean regions
Transect-based MOC computation across specified transects
Configurable temporal output (reduction periods and snapshots)
IOStream integration for netCDF output

MOC Computation Pipeline

Latitude-binned Regional MOC chain:

Assign cells to latitude bins (static, initialization)
Convert vertical pseudo-velocity to geometric coordinates
Compute vertical flux (velocity × area)
Apply regional mask (optional)
Accumulate into latitude bins
Horizontal integration (south→north)
Convert to Sverdrups

Transect-based MOC chain:

Convert pseudo-thickness to geometric layer thickness
Compute transport (thickness × velocity × edge width)
Accumulate across transect edges
Vertical integration (bottom→top)
Convert to Sverdrups

Technical Implementation

Design Features:

Field-level regional mask storage and automatic propagation
IOName metadata system separating internal from user-facing names
Kokkos hierarchical parallelism for 2D/3D operators
SFINAE compile-time optimization to prevent invalid template instantiations
Operator parameters passed via Config objects through chain parsing

Output Format:

MOC streamfunction in Sverdrups (1 Sv = 10⁶ m³/s)
Dimensions: latitude bins × depth levels (regional), or depth only (transect)
NetCDF files with user-friendly variable names (e.g., "MOC_streamfunction_Global")
Configurable temporal averaging and snapshot output

Limitations

Regional masks not yet implemented (operator infrastructure ready, placeholder in config)
Transect masks not yet implemented (operator infrastructure ready, placeholder in config)

Checklist

xylar · 2026-07-28T21:01:17Z

Sorry, @brian-oneill, I didn't get to this today. I'll try again tomorrow.

Also, let me know what you need from me regarding both the dynamic streams here in Omega and the Polaris support.

cbegeman · 2026-07-28T23:10:15Z

+  ReductionPeriod: [1Month]
+  SnapshotPeriod: [1Day]


Can you help me understand what sets how often the MOC is computed? Is it the SnapshotPeriod? With the options above, ReductionPeriod would then be averaging daily instantaneous MOC values over 1 month?

My understanding is that they are independent of one another. In practice, we would have:

ReductionPeriod: [1Month] SnapshotPeriod: []

since we want monthly averages and don't need snapshots.

It's hard for me to imagine very much analysis where we want both time averages and snapshots at the same time, in practice.

Yes, ReductionPeriod and SnapshotPeriod produce outputs independently. The ReductionPeriod outputs currently accumulate every timestep so the MOC is computed each timestep for time averages, but this can be extended to allow for a courser sampling frequency pretty easily.

cbegeman · 2026-07-28T23:16:29Z

+      # 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]


Do we need to test that when, e.g., LayerThickness_BinaryMultiply(NormalVelocity), is present here that the MOC chain uses the available field or is this kind of thing covered by existing CTests?

Well, I guess not present here because GlobalStats reduces spatially.

During the initial parsing, the parser checks if a Field that would be output by an operator has already been registered, to prevent building a duplicate operator. There is a unit test that checks this behavior, but it could be more robust.

xylar

@brian-oneill, I'll do some testing but here are a few comments to keep the process moving.

This looks great! A lot of the pieces are in place and just a few tweaks would be helpful, I think. Plus a few things that might be for now or might be postponed until later.

xylar · 2026-07-29T10:57:11Z

+  ReductionPeriod: [1Month]
+  SnapshotPeriod: [1Day]


My understanding is that they are independent of one another. In practice, we would have:

ReductionPeriod: [1Month] SnapshotPeriod: []

since we want monthly averages and don't need snapshots.

It's hard for me to imagine very much analysis where we want both time averages and snapshots at the same time, in practice.

xylar · 2026-07-29T11:28:18Z

Could ScalarMultiply be generalized to take a model config option or known constant as its input, not just a hard-coded number? This would seem much more useful and general.

For the use case here, this Op gets passed 1e-6 through a "config" that is programmatically defined in the parser. So 1e-6 is hard-coded into the MOC chain construction, but the Op itself is designed to take a configurable value to be compatible with the future composable framework.

Sounds good. It's RhoSw in particular that I wanted to know about.

Gotcha, that's possible but requires a little more work. Allowing users to request a defined constant from the config file would require defining a map between string labels and the variable names in GlobalConstants.h

std::map<std::string,Real> Constants = { {"RhoSw", RhoSw}, {"Gravity", Gravity}, {"Pi", Pi}, ... };

I think we should try to figure out how to automate that somehow but, yes, that's what I was anticipating. Nothing that needs to be in this PR.

Yeah, automated would be best. Here's Claude's suggestion:

Yes, the X-macro pattern is exactly right. The idea: define a single list of constants once using a macro, then expand it in two ways — once to declare the constexpr values, and once to build the lookup.

In GlobalConstants.h, replace the individual constexpr declarations for the physical constants with:

// X-macro list: X(Name, Value) #define OMEGA_PHYSICAL_CONSTANTS(X) \ X(RhoSw, pcd::seawater_density_reference) \ X(RhoFw, pcd::pure_water_density_reference) \ X(RhoAir, pcd::dry_air_density_at_standard_temperature_and_pressure) \ X(Gravity, pcd::standard_acceleration_of_gravity) \ X(RhoIce, pcd::sea_ice_density_reference) \ /* ... all others ... */ // Expand to constexpr declarations (same as before) #define DECLARE_CONST(Name, Value) constexpr Real Name = (Value); OMEGA_PHYSICAL_CONSTANTS(DECLARE_CONST) #undef DECLARE_CONST

Then the lookup function writes itself:

inline std::optional<Real> getConstantByName(const std::string &Name) { #define MATCH_CONST(CName, Value) if (Name == #CName) return CName; OMEGA_PHYSICAL_CONSTANTS(MATCH_CONST) #undef MATCH_CONST return std::nullopt; }

#CName stringifies the identifier automatically, so the name in the config file ("RhoSw") matches the C++ variable name without any manual duplication.

Pros: Zero maintenance burden — adding a new constant to the list automatically makes it available by name to ScalarMultiplyOp and any future operator that calls getConstantByName.

Cons: Requires refactoring the existing constexpr declarations in GlobalConstants.h into the macro list format. The math-derived ones (TwoPi, SDay, etc.) are trickier since they depend on other constants — those would need to stay as regular constexpr or be added after the macro expansion.

The practical approach: put only the "leaf" physical constants (densities, heat capacities, etc.) in the X-macro list, and keep derived/compound ones (TwoPi, TkFrzSw, SDay) as regular constexpr declarations below. You could optionally add those to a second X-macro list if you want them accessible by name too.

Not bad, not bad!

Allow for an optional boundary condition and add a 1D compute method

Remove horizontal dimension assertion and some refactoring

xylar

I successfully ran a 5-day test with the MOC on in:
/lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260804/ec30to60-global-moc/ocean/spherical/realistic_global/EC30to60E2r2/analysis_members_test

I used:

    MOC:
      Enable: true
      # 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: [1day]       # Temporal reduction periods
      SnapshotPeriod: []  # Instantaneous output periods
      Filename: moc.$Y
      Stream:
        FileFreq: 1
        FileFreqUnits: days

So daily averaging rather than monthly for efficiency.

Here's an example plot:

However, the latitude bins and depth are missing from the output file:

$ ncdump -h moc_1dayTimeStats.0001 
netcdf moc_1dayTimeStats {
dimensions:
	MaxCellsOnEdge = 2 ;
	MaxEdges = 7 ;
	MaxEdges2 = 14 ;
	NCells = 236853 ;
	NEdges = 719506 ;
	NTracers = 2 ;
	NVertLayers = 60 ;
	NVertLayersP1 = 61 ;
	NVertices = 482371 ;
	NumBinsLatCell_BinIndex = 180 ;
	Scalar = 1 ;
	VertexDegree = 3 ;
	time = UNLIMITED ; // (5 currently)
variables:
	double MOC_streamfunction_Global_TimeMean1day(time, NumBinsLatCell_BinIndex, NVertLayersP1) ;
		MOC_streamfunction_Global_TimeMean1day:Description = "Time average of VerticalPseudoVelocity_PseudoToGeometric_BinaryMultiply(AreaCell)_BinnedAccumulator(LatCell_BinIndex)_PrefixSum_ScalarMultiply(1.0e-6)" ;
		MOC_streamfunction_Global_TimeMean1day:Name = "VerticalPseudoVelocity_PseudoToGeometric_BinaryMultiply(AreaCell)_BinnedAccumulator(LatCell_BinIndex)_PrefixSum_ScalarMultiply(1.0e-6)_TimeMean1day" ;
		MOC_streamfunction_Global_TimeMean1day:StdName = "" ;
		MOC_streamfunction_Global_TimeMean1day:Units = "" ;
		MOC_streamfunction_Global_TimeMean1day:ValidMax = 1.79769313486232e+308 ;
		MOC_streamfunction_Global_TimeMean1day:ValidMin = -1.79769313486232e+308 ;
		MOC_streamfunction_Global_TimeMean1day:_FillValue = 9.96920996838687e+36 ;
		MOC_streamfunction_Global_TimeMean1day:long_name = "Time average of VerticalPseudoVelocity_PseudoToGeometric_BinaryMultiply(AreaCell)_BinnedAccumulator(LatCell_BinIndex)_PrefixSum_ScalarMultiply(1.0e-6)" ;
		MOC_streamfunction_Global_TimeMean1day:name = "VerticalPseudoVelocity_PseudoToGeometric_BinaryMultiply(AreaCell)_BinnedAccumulator(LatCell_BinIndex)_PrefixSum_ScalarMultiply(1.0e-6)_TimeMean1day" ;
		MOC_streamfunction_Global_TimeMean1day:standard_name = "" ;
		MOC_streamfunction_Global_TimeMean1day:units = "" ;
		MOC_streamfunction_Global_TimeMean1day:valid_max = 1.79769313486232e+308 ;
		MOC_streamfunction_Global_TimeMean1day:valid_min = -1.79769313486232e+308 ;
	double time(time) ;
		time:Description = "time" ;
		time:Name = "time" ;
		time:StdName = "time" ;
		time:Units = "seconds since 0001-01-01 00:00:00" ;
		time:ValidMax = 1.e+20 ;
		time:ValidMin = 0. ;
		time:_FillValue = 9.96920996838687e+36 ;
		time:calendar = "noleap" ;
		time:long_name = "time" ;
		time:name = "time" ;
		time:standard_name = "time" ;
		time:units = "seconds since 0001-01-01 00:00:00" ;
		time:valid_max = 1.e+20 ;
		time:valid_min = 0. ;

// global attributes:
		:SimulationTime = "0001-01-06_00:00:00" ;
		:SimulationTime0 = "0001-01-02_00:00:00" ;
		:SimulationTime1 = "0001-01-03_00:00:00" ;
		:SimulationTime2 = "0001-01-04_00:00:00" ;
		:SimulationTime3 = "0001-01-05_00:00:00" ;
		:SimulationTime4 = "0001-01-06_00:00:00" ;
}

Also the file is missing a .nc extension.

Happy to rerun once this is fixed.

brian-oneill · 2026-08-04T17:36:12Z

Is refBottomDepth the appropriate depth to go with the output? Because that's not currently read in to Omega...

I've noticed the history and hifreq files that get created when running the omega ctests also no longer get the .nc suffix added. Looks like 1f92657 removed the lines that explicitly added the suffix to filenames with time templates. A comment says PIO should handle this, but that doesn't appear to be the case. Adding .nc to the end of the Filename option in the config does work.

xylar · 2026-08-04T17:39:24Z

refBottomDepth is an MPAS-Ocean concept and is not going to be available in general for Omega. So, no, it is not appropriate.

The problem with the current layer-wise approach is it implicitly assumes pure z-level layers. The best way to get depths given this approach is to get the area-weighted average of zInterface and use that as the vertical coordinate. For now, that's fine. With ice-shelf cavities, sigma coordinates, etc. in the future, we'll need to do vertical binning or interpolation to a z-level or density-level grid.

xylar · 2026-08-04T17:40:53Z

I'm fine with whatever fix to the .nc extension but we do need them to be put back somehow. I can understand the desire to support arbitrary format and not explicitly require NetCDF, though.

…streamfunction' into omega/analysis-moc-streamfunction

cbegeman · 2026-08-04T22:22:11Z

+                                # List of field names to compute statistics for
+      SpatialStats: [Max, Min, Mean, StdDev]
+                                # Spatial statistics to compute (one per field)
+      ReductionPeriod: [1Day, 1Month]


Suggested change

ReductionPeriod: [1Day, 1Month]

ReductionPeriod: []

Should we remove reductions altogether from this analysis member to save compute time?

I think we want 1Month reduction for the climatology plot, don't we?

I guess it depends on what configuration Default.yml is aimed at -- typical standalone testing or a longer production run.

cbegeman · 2026-08-04T22:23:13Z

+                                # Spatial statistics to compute (one per field)
+      ReductionPeriod: [1Day, 1Month]
+                                # Temporal reduction periods (time-averaged stats)
+      SnapshotPeriod: [6Hours]


MPAS-O is equivalent to 1Day but I like increasing the freq to 6 hours

cbegeman · 2026-08-10T15:00:47Z

However, the latitude bins and depth are missing from the output file:

@brian-oneill Is this ready for testing with the fix to this comment?

brian-oneill · 2026-08-11T17:11:44Z

However, the latitude bins and depth are missing from the output file:

@brian-oneill Is this ready for testing with the fix to this comment?

@cbegeman I added the latitude bins. The depths will require a bit more thought and work, and I won’t have time to complete that piece until I’m back at the end of the week.

brian-oneill added 21 commits July 25, 2026 16:01

Add regional mask support to Field

10daf5a

Add CoordinateBinningOp

c9e7fa9

Add PseudoToGeometricOp

f02df5d

Add BinaryMultiplyOp

8d187ab

Add BinnedAccumulatorOp

11602aa

Add PrefixSumOp

1aae717

Add ScalarMultiplyOp

54e14f1

Add ExtractRegionOp.h

e748f7a

Add register1D/2D/2DReal variant methods to AnalysisOpFactory

377432c

Register new ops

0b342ad

Add TransectAccumulatorOp

b86d187

Parsing fixes mostly

65756d1

Create MOC AnalysisGroup

61b207e

Edits to AnalysisGroup methods to pass configs for Op chains

969cf07

Update parser to handle new Ops

e1b4e37

AnalysisOperatorTest helper updates

5e2c5fb

Unit tests for new Ops

dc9f375

Add MOC node to Default.yml config file

5e8db1c

IOName capability to provide legible Analysis Field names in IOStreams

893f042

Documentation updates

693cf53

BinnedAccumulatorOp linting

bda972e

brian-oneill requested review from cbegeman and xylar July 28, 2026 04:01

cbegeman reviewed Jul 28, 2026

View reviewed changes

xylar reviewed Jul 29, 2026

View reviewed changes

brian-oneill added 3 commits August 3, 2026 21:38

Edit parser to handle optional BC in PrefixSumOp

d9311ed

Add defaults for FileFreq and FileFreqUnits to StreamParams

7a49875

Update MOC AnalysisGroup

c982478

brian-oneill added 9 commits August 3, 2026 22:02

Add OutputName optional override method for AnalysisOperator

128fed3

Fix ExtractRegionOp OutputName

4b8c155

PrefixSumOp updates

14cde2f

Allow for an optional boundary condition and add a 1D compute method

PseudoToGeometricOp Updates

bfc9c08

Remove horizontal dimension assertion and some refactoring

Fix TransectAccumulatorOp OutputName

3440bba

Updates to AnalysisOperatorTest

b17bc41

Add MOC chain test to AnalysisSystemTest

7298b61

Make some Analysis class methods public

d631f95

Update default Config

0e2fb94

xylar reviewed Aug 4, 2026

View reviewed changes

Add bin boundaries to MOC output

9a7be85

brian-oneill added 2 commits August 4, 2026 10:57

Functions with parallel dispatch cannot be private

1d570f3

Merge remote-tracking branch 'refs/remotes/origin/omega/analysis-moc-…

39fde8e

…streamfunction' into omega/analysis-moc-streamfunction

cbegeman reviewed Aug 4, 2026

View reviewed changes

Uh oh!

Conversation

brian-oneill commented Jul 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Overview

Key Features

MOC Computation Pipeline

Technical Implementation

Limitations

Uh oh!

xylar commented Jul 28, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

xylar left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

xylar left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

brian-oneill commented Aug 4, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

xylar commented Aug 4, 2026

Uh oh!

xylar commented Aug 4, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cbegeman commented Aug 10, 2026

Uh oh!

brian-oneill commented Aug 11, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

brian-oneill commented Jul 28, 2026 •

edited

Loading

xylar left a comment •

edited

Loading

brian-oneill commented Aug 4, 2026 •

edited

Loading

xylar commented Aug 4, 2026 •

edited

Loading