An Exception was encountered at 'In [6]'.

Ocean.stats diagnostics (ocean.stats.nc)¶

MOM6 writes a per-timestep integration-diagnostics file, ocean.stats (plain text) and ocean.stats.nc (its netCDF equivalent) -- energy, max CFL, mean sea level, total mass, mean salinity/temperature, and their conservation errors over the course of the run. This notebook finds ocean.stats.nc via the intake-esm datastore and plots each diagnostic.

See discussion in access-om3-paper-1#23 and access-nri-intake-catalog#682.

The first draft of thie notebook was written with Claude (prompted by @chrisb13)

In [1]:

parameters

# These first two cells must be in all notebooks!
# It allows us to run all the notebooks at once, this cell has a tag "parameters" which allows us to pass in
# arguments externally using papermill (see mkfigs.sh for details)

# Set esm_file to the datastore for the main experiment of interest
esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test4-d28e0359/datastore.json"

# papermill settings. *No need to modify these if running interactively.*
papermill = False                      # `cwd` and `nbname` will be populated by papermill.
cwd = None                             # current working directory
nbname = None                          # notebook name

# These first two cells must be in all notebooks! # It allows us to run all the notebooks at once, this cell has a tag "parameters" which allows us to pass in # arguments externally using papermill (see mkfigs.sh for details) # Set esm_file to the datastore for the main experiment of interest esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test4-d28e0359/datastore.json" # papermill settings. *No need to modify these if running interactively.* papermill = False # `cwd` and `nbname` will be populated by papermill. cwd = None # current working directory nbname = None # notebook name

In [2]:

injected-parameters

# Parameters
esm_file = "/jobfs/181564821.gadi-pbs/cm3-run-27-07-2026-PI-control/datastore.json"
papermill = True
cwd = "/g/data/tm70/cyb561/access-om3-live/access-om3-paper-1/notebooks/mkfigs_output_cm3-run-27-07-2026-PI-control/"
nbname = "ocean_stats.ipynb"

# Parameters esm_file = "/jobfs/181564821.gadi-pbs/cm3-run-27-07-2026-PI-control/datastore.json" papermill = True cwd = "/g/data/tm70/cyb561/access-om3-live/access-om3-paper-1/notebooks/mkfigs_output_cm3-run-27-07-2026-PI-control/" nbname = "ocean_stats.ipynb"

In [3]:

if not papermill:
    import nci_ipynb, os  # requires conda/analysis3-26.03 or later
    cwd = nci_ipynb.dir()
    nbname = nci_ipynb.name()
    os.chdir(cwd)
import mkfigs_bootstrap  # noqa: adds external/access-model-mkfigs/src to sys.path (stop-gap)
from mkfigs import MkmdWriter
mkmd = MkmdWriter(esm_file, nbname, str(cwd), pm=papermill)

if not papermill: import nci_ipynb, os # requires conda/analysis3-26.03 or later cwd = nci_ipynb.dir() nbname = nci_ipynb.name() os.chdir(cwd) import mkfigs_bootstrap # noqa: adds external/access-model-mkfigs/src to sys.path (stop-gap) from mkfigs import MkmdWriter mkmd = MkmdWriter(esm_file, nbname, str(cwd), pm=papermill)

In [4]:

import intake
import matplotlib.pyplot as plt
import xarray as xr
import numpy as np

import intake import matplotlib.pyplot as plt import xarray as xr import numpy as np

Open the intake-esm datastore and find ocean.stats.nc¶

In [5]:

COLUMNS_WITH_ITERABLES = [
    "variable",
    "variable_long_name",
    "variable_standard_name",
    "variable_cell_methods",
    "variable_units",
]

datastore = intake.open_esm_datastore(esm_file, columns_with_iterables=COLUMNS_WITH_ITERABLES)
datastore.search(filename="ocean.stats.nc").df

COLUMNS_WITH_ITERABLES = [ "variable", "variable_long_name", "variable_standard_name", "variable_cell_methods", "variable_units", ] datastore = intake.open_esm_datastore(esm_file, columns_with_iterables=COLUMNS_WITH_ITERABLES) datastore.search(filename="ocean.stats.nc").df

Out[5]:

	filename	path	file_id	frequency	start_date	end_date	variable	variable_long_name	variable_standard_name	variable_cell_methods	variable_units	realm	temporal_label

Load ocean.stats.nc¶

Loading ocean.stats.nc currently needs two workarounds:

1. The datastore's aggregation only recognises a lowercase time dimension when deciding how to join assets, but ocean.stats.nc uses Time. This means the 66 per-restart ocean.stats.nc files aren't concatenated -- they get grouped by file_id instead, based on each file's number of timesteps, collapsing everything into just two datasets. aggregate_on_keys below removes file_id from the aggregation so all matching files are joined into one dataset (see access-nri-intake-catalog#682 and this PR discussion for the full investigation).
2. Time's units attribute looks timedelta-like to CF conventions, so xarray tries to decode it as timedelta64[ns] -- which overflows for a multi-century run. decode_times=False avoids that.

Execution using papermill encountered an exception here and stopped:

In [6]:

def aggregate_on_keys(ds, keys):
    """Drop `keys` from the datastore's aggregation groupby_attrs so matching
    assets are joined into one dataset instead of being split by those keys
    (workaround for the ocean.stats.nc `Time` vs `time` capitalisation bug --
    see access-nri-intake-catalog#682)."""
    parent_keys = ds.esmcat.aggregation_control.groupby_attrs
    ds.esmcat.aggregation_control.groupby_attrs = [
        key for key in parent_keys if key not in keys
    ]
    return ds


ds_filtered = datastore.search(filename="ocean.stats.nc")
ds_filtered = aggregate_on_keys(ds_filtered, ["file_id"])

ocean_stats = ds_filtered.to_dask(xarray_open_kwargs=dict(decode_times=False))
ocean_stats

def aggregate_on_keys(ds, keys): """Drop `keys` from the datastore's aggregation groupby_attrs so matching assets are joined into one dataset instead of being split by those keys (workaround for the ocean.stats.nc `Time` vs `time` capitalisation bug -- see access-nri-intake-catalog#682).""" parent_keys = ds.esmcat.aggregation_control.groupby_attrs ds.esmcat.aggregation_control.groupby_attrs = [ key for key in parent_keys if key not in keys ] return ds ds_filtered = datastore.search(filename="ocean.stats.nc") ds_filtered = aggregate_on_keys(ds_filtered, ["file_id"]) ocean_stats = ds_filtered.to_dask(xarray_open_kwargs=dict(decode_times=False)) ocean_stats

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[6], line 16
     12 
     13 ds_filtered = datastore.search(filename="ocean.stats.nc")
     14 ds_filtered = aggregate_on_keys(ds_filtered, ["file_id"])
     15 
---> 16 ocean_stats = ds_filtered.to_dask(xarray_open_kwargs=dict(decode_times=False))
     17 ocean_stats

File /g/data/xp65/public/apps/med_conda/envs/analysis3-26.08/lib/python3.12/site-packages/intake_esm/core.py:926, in esm_datastore.to_dask(self, **kwargs)
    924 if len(self) != 1:  # quick check to fail more quickly if there are many results
    925     lens = self._aggr_lengths
--> 926     raise ValueError(
    927         f'Expected exactly one dataset. Received {len(self)} datasets. Please refine your search on {", ".join(lens.keys())} or use `.to_dataset_dict()`.'
    928     )
    929 res = self.to_dataset_dict(**{**kwargs, 'progressbar': False})
    930 if len(res) != 1:  # extra check in case kwargs did modify something

ValueError: Expected exactly one dataset. Received 0 datasets. Please refine your search on  or use `.to_dataset_dict()`.

Plot each diagnostic over the run¶

In [ ]:

# Time is raw elapsed days (decode_times=False avoids the calendar overflow).
# Convert to years for a more readable axis -- check ocean_stats["Time"].attrs
# for the exact calendar; 365.25 is a reasonable default, use 365 if it's a
# 'noleap' calendar (common for JRA-forced runs).
days_per_year = 365.25
ocean_stats = ocean_stats.assign_coords(Time_years=ocean_stats["Time"] / days_per_year)

data_vars = [v for v in ocean_stats.data_vars if "Time" in ocean_stats[v].dims]

ncols = 3
nrows = -(-len(data_vars) // ncols)  # ceil division
fig, axes = plt.subplots(nrows, ncols, figsize=(5 * ncols, 3 * nrows), squeeze=False)

for ax, var in zip(axes.flat, data_vars):
    ocean_stats[var].plot(ax=ax, x="Time_years")
    ax.set_title(var)
    ax.set_xlabel("Time (years)")

for ax in axes.flat[len(data_vars):]:
    ax.axis("off")

fig.suptitle("ACCESS-OM3 ocean.stats.nc diagnostics", y=1.02)
fig.tight_layout()

# Time is raw elapsed days (decode_times=False avoids the calendar overflow). # Convert to years for a more readable axis -- check ocean_stats["Time"].attrs # for the exact calendar; 365.25 is a reasonable default, use 365 if it's a # 'noleap' calendar (common for JRA-forced runs). days_per_year = 365.25 ocean_stats = ocean_stats.assign_coords(Time_years=ocean_stats["Time"] / days_per_year) data_vars = [v for v in ocean_stats.data_vars if "Time" in ocean_stats[v].dims] ncols = 3 nrows = -(-len(data_vars) // ncols) # ceil division fig, axes = plt.subplots(nrows, ncols, figsize=(5 * ncols, 3 * nrows), squeeze=False) for ax, var in zip(axes.flat, data_vars): ocean_stats[var].plot(ax=ax, x="Time_years") ax.set_title(var) ax.set_xlabel("Time (years)") for ax in axes.flat[len(data_vars):]: ax.axis("off") fig.suptitle("ACCESS-OM3 ocean.stats.nc diagnostics", y=1.02) fig.tight_layout()

In [ ]:

# mkmd.savefig(fig, title, caption, dpi=100)
# Save figure and append to markdown summary.
mkmd.savefig(
    fig,
    "Ocean.stats diagnostics",
    "MOM6 integration diagnostics from ocean.stats.nc (energy, max CFL, mean sea "
    "level, total mass, mean salinity/temperature, and their conservation errors) "
    "over the course of the run.",
    dpi=150,
)

# mkmd.savefig(fig, title, caption, dpi=100) # Save figure and append to markdown summary. mkmd.savefig( fig, "Ocean.stats diagnostics", "MOM6 integration diagnostics from ocean.stats.nc (energy, max CFL, mean sea " "level, total mass, mean salinity/temperature, and their conservation errors) " "over the course of the run.", dpi=150, )

In [ ]: