Running xhycom at scale

Regridding or averaging full HYCOM archives is compute-intensive: a single year of daily fields on the TP2 grid to GLORYS is ~1.6 TB, and streaming through all daily files to build monthly means takes the same wall time regardless of how many output steps you want. For multi-year or multi-variable runs, a Slurm batch job is more practical than an interactive notebook kernel.

This page covers the most common pre-compute pattern — monthly means on a fixed grid — and points to real-world Slurm script examples.

Setup

The cells below pick up where Regridding left off — the same HYCOM dataset, GLORYS mask, and cached regrid weights.

import numpy as np
import xarray as xr
import xhycom
import dask
from dask.diagnostics import ProgressBar

import os
os.environ["XHYCOM_CACHE_DIR"] = f"/cluster/projects/nn2993k/{os.environ['USER']}/.xhycom-cache-dir"

grid = "/cluster/home/nlo043/NERSC-HYCOM-CICE/TP2a0.10/topo/regional.grid"
DATA_PATH = "/nird/datalake/NS9481K/shuang/TP2_output/expt_02.8/"

ds = xhycom.open_mfdataset(DATA_PATH + "archm.2020*", grid=grid, chunks={"time": 1}, postprocess=True)

GRIDS = "/nird/datapeak/NS9481K/MERCATOR_DATA/REGULAR_GRID_COORD"
glorys = xr.open_dataset(f"{GRIDS}/GLO-MFC_001_030_mask_bathy.nc")

Computing monthly means

The full daily regrid is large (~1.6 TB for four variables) and rarely what you want — the usual target is monthly means on the GLORYS grid. For two variables (temp and salin) that is ~26 GB for one year and takes around 10 minutes to write on Olivia.

How you build the monthly mean is a methodology question in its own right — average on the native layers (thickness-weighted, the way HYCOM averages online) or at fixed depth — and the choice matters near fronts and the seasonal pycnocline. That’s covered in Time-averaging on a moving vertical grid; here we use the recommended thickness-weighted route and write it to disk.

Even the monthly-mean route streams through all the daily fields to build the average, so runtime scales with the full time range, not just the 12 output steps. A single year interactively in a notebook is fine; a multi-decade run belongs in a Slurm script — see below.

# Recommended monthly mean: thickness-weighted on the native layers (the HYCOM
# monthly mean — see time-averaging.ipynb), then a single regrid. Cheap: ~12 remaps.
h = ds["thknss"]
hbar = h.resample(time="1MS").mean()
monthly = (ds[["temp", "salin"]] * h).resample(time="1MS").mean() / hbar.where(hbar > 0)
monthly["thknss"] = hbar
g_month = xhycom.regrid(monthly, target=glorys, grid=grid, weights=True)

g_month   # 12 monthly means on the GLORYS grid

<xarray.Dataset> Size: 26GB
Dimensions:  (time: 12, lat: 624, lon: 4320, depth: 50)
Coordinates:
  * time     (time) object 96B 2020-01-01 00:00:00 ... 2020-12-01 00:00:00
  * lat      (lat) float32 2kB 38.08 38.17 38.25 38.33 ... 89.83 89.92 90.0
  * lon      (lon) float32 17kB -180.0 -179.9 -179.8 ... 179.8 179.8 179.9
  * depth    (depth) float64 400B 0.494 1.541 2.646 ... 5.275e+03 5.728e+03
Data variables:
    temp     (time, lat, lon, depth) float64 13GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
    salin    (time, lat, lon, depth) float64 13GB dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>
Attributes:
    iversn:         23
    iexpt:          28
    yrflag:         3
    archive_type:   mean
    long_name:      layer thickness
    units:          m
    comment:        converted from Pa (factor 0.000101978)
    regrid_method:  conservative

xarray.Dataset

• Dimensions:
  • time: 12
  • lat: 624
  • lon: 4320
  • depth: 50
• Coordinates: (4)
  • time
    (time)
    object
    2020-01-01 00:00:00 ... 2020-12-...

    long_name :

    time

    axis :

    T

    array([cftime.DatetimeGregorian(2020, 1, 1, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2020, 2, 1, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2020, 3, 1, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2020, 4, 1, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2020, 5, 1, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2020, 6, 1, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2020, 7, 1, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2020, 8, 1, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2020, 9, 1, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2020, 10, 1, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2020, 11, 1, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2020, 12, 1, 0, 0, 0, 0, has_year_zero=False)],
          dtype=object)

  • lat
    (lat)
    float32
    38.08 38.17 38.25 ... 89.92 90.0

    standard_name :

    latitude

    units :

    degrees_north

    array([38.083332, 38.166668, 38.25    , ..., 89.833336, 89.916664, 90.      ],
          shape=(624,), dtype=float32)

  • lon
    (lon)
    float32
    -180.0 -179.9 ... 179.8 179.9

    standard_name :

    longitude

    units :

    degrees_east

    array([-180.     , -179.91667, -179.83333, ...,  179.75   ,  179.83333,
            179.91667], shape=(4320,), dtype=float32)

  • depth
    (depth)
    float64
    0.494 1.541 ... 5.275e+03 5.728e+03

    long_name :

    depth

    units :

    m

    positive :

    down

    axis :

    Z

    array([4.940250e-01, 1.541375e+00, 2.645669e+00, 3.819495e+00, 5.078224e+00,
           6.440614e+00, 7.929560e+00, 9.572997e+00, 1.140500e+01, 1.346714e+01,
           1.581007e+01, 1.849556e+01, 2.159882e+01, 2.521141e+01, 2.944473e+01,
           3.443415e+01, 4.034405e+01, 4.737369e+01, 5.576429e+01, 6.580727e+01,
           7.785385e+01, 9.232607e+01, 1.097293e+02, 1.306660e+02, 1.558507e+02,
           1.861256e+02, 2.224752e+02, 2.660403e+02, 3.181274e+02, 3.802130e+02,
           4.539377e+02, 5.410889e+02, 6.435668e+02, 7.633331e+02, 9.023393e+02,
           1.062440e+03, 1.245291e+03, 1.452251e+03, 1.684284e+03, 1.941893e+03,
           2.225078e+03, 2.533336e+03, 2.865703e+03, 3.220820e+03, 3.597032e+03,
           3.992484e+03, 4.405224e+03, 4.833291e+03, 5.274784e+03, 5.727917e+03])

• Data variables: (2)
  • temp
    (time, lat, lon, depth)
    float64
    dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>

    comment :

    converted from Pa (factor 0.000101978)

    long_name :

    layer thickness

    units :

    m

    Array Chunk Bytes 12.05 GiB 1.00 GiB Shape (12, 624, 4320, 50) (1, 624, 4320, 50) Dask graph 12 chunks in 1643 graph layers Data type float64 numpy.ndarray

  • salin
    (time, lat, lon, depth)
    float64
    dask.array<chunksize=(1, 624, 4320, 50), meta=np.ndarray>

    comment :

    converted from Pa (factor 0.000101978)

    long_name :

    layer thickness

    units :

    m

    Array Chunk Bytes 12.05 GiB 1.00 GiB Shape (12, 624, 4320, 50) (1, 624, 4320, 50) Dask graph 12 chunks in 1643 graph layers Data type float64 numpy.ndarray

• Attributes: (8)

  iversn :

  23

  iexpt :

  28

  yrflag :

  3

  archive_type :

  mean

  long_name :

  layer thickness

  units :

  m

  comment :

  converted from Pa (factor 0.000101978)

  regrid_method :

  conservative

# Cap workers to keep memory bounded — the ~GB regridded chunks can outrun the
# single-locked netCDF writer at high parallelism.
with dask.config.set(num_workers=4), ProgressBar():
    g_month.to_netcdf("/cluster/work/projects/nn2993k/nlo043/hycom_on_glorys_2020_monthly.nc")

[########################################] | 100% Completed | 599.47 s

Slurm script examples

The OHC-GLORYS-HYCOM repository contains worked examples of running xhycom regridding in Slurm batch jobs on Sigma2 Olivia. Start there when setting up a multi-year regrid pipeline.