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.