Getting started with xgcm for MOM6¶
MOM6 discretises the ocean on a staggered
Arakawa C-grid with a north-east index convention: tracers sit at cell
centres, zonal velocity uo on the east face (U-point), and meridional
velocity vo on the north face (V-point). Global MOM6 configurations use a
tripolar grid, so the true geographic longitude/latitude are 2-D curvilinear
fields — not the 1-D index coordinates.
Here we use real GFDL-CM4 output (MOM6 ocean) from the CMIP6 archive, read
anonymously from the Pangeo Google Cloud Zarr store. We take the surface level
of one month of the historical run.
Dependencies — reading the cloud store needs
zarrandgcsfs(pip install zarr gcsfs); the maps usecartopy.
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from xgcm import Grid
base = ("gs://cmip6/CMIP6/CMIP/NOAA-GFDL/GFDL-CM4/historical/"
"r1i1p1f1/Omon/{var}/gn/v20180701/")
so = {"storage_options": {"token": "anon"}}
def surface(var):
da = xr.open_dataset(base.format(var=var), engine="zarr", backend_kwargs=so)
return da[var].isel(time=0, lev=0)
thetao = surface("thetao") # sea-water potential temperature (centre)
uo = surface("uo") # zonal velocity (U-point, east face)
vo = surface("vo") # meridional velocity (V-point, north face)
# Geographic cell-centre coordinates for plotting (2-D, curvilinear).
lon = thetao["lon"].values
lat = thetao["lat"].values
thetao
<xarray.DataArray 'thetao' (y: 1080, x: 1440)> Size: 6MB
[1555200 values with dtype=float32]
Coordinates:
* y (y) float64 9kB -80.39 -80.31 -80.23 -80.15 ... 89.73 89.84 89.95
* x (x) float64 12kB -299.7 -299.5 -299.2 -299.0 ... 59.53 59.78 60.03
lat (y, x) float32 6MB -79.81 -79.81 -79.8 -79.8 ... 64.33 64.22 64.11
lon (y, x) float32 6MB -299.7 -299.5 -299.2 -299.0 ... 60.0 60.0 60.0
lev float64 8B 2.5
time object 8B 1850-01-16 12:00:00
Attributes:
cell_measures: area: areacello volume: volcello
cell_methods: area: mean where sea time: mean
long_name: Sea Water Potential Temperature
original_name: thetao
standard_name: sea_water_potential_temperature
units: degCDefining the xgcm Grid¶
The CMIP6 file stores uo, vo and thetao on the same (y, x) index, but
they live at different staggered positions. We give each axis a center and a
(half-cell-shifted) right position, place uo on the X-right (east) face and
vo on the Y-right (north) face, and treat X as periodic.
ny, nx = thetao.shape
coords = dict(x_c=np.arange(nx), x_f=np.arange(nx), y_c=np.arange(ny), y_f=np.arange(ny))
T = xr.DataArray(thetao.values, dims=["y_c", "x_c"])
u = xr.DataArray(uo.values, dims=["y_c", "x_f"]).assign_coords(x_f=coords["x_f"], y_c=coords["y_c"])
v = xr.DataArray(vo.values, dims=["y_f", "x_c"]).assign_coords(x_c=coords["x_c"], y_f=coords["y_f"])
ds = xr.Dataset(dict(thetao=T, uo=u, vo=v), coords=coords)
grid = Grid(
ds,
coords={"X": {"center": "x_c", "right": "x_f"},
"Y": {"center": "y_c", "right": "y_f"}},
boundary={"X": "periodic", "Y": "extend"},
autoparse_metadata=False,
)
grid
<xgcm.Grid> X Axis (periodic, boundary='periodic'): * center x_c --> right * right x_f --> center Y Axis (not periodic, boundary='extend'): * center y_c --> right * right y_f --> center
A staggered-grid operation: surface current speed¶
A common need is to bring the C-grid velocity components to a common point.
We interpolate uo and vo from their faces to the cell centre and form the
surface current speed $\sqrt{u^2 + v^2}$ there — a one-liner with xgcm.
uc = grid.interp(ds.uo, "X") # east face -> centre
vc = grid.interp(ds.vo, "Y") # north face -> centre
speed = np.hypot(uc, vc).load()
speed
<xarray.DataArray 'uo' (y_c: 1080, x_c: 1440)> Size: 6MB
array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]],
shape=(1080, 1440), dtype=float32)
Coordinates:
* y_c (y_c) int64 9kB 0 1 2 3 4 5 6 ... 1074 1075 1076 1077 1078 1079
* x_c (x_c) int64 12kB 0 1 2 3 4 5 6 ... 1434 1435 1436 1437 1438 1439fig = plt.figure(figsize=(11, 6))
ax = plt.axes(projection=ccrs.Robinson(central_longitude=-150))
pm = ax.pcolormesh(lon, lat, speed.values, transform=ccrs.PlateCarree(),
cmap="magma_r", vmin=0, vmax=0.6)
ax.coastlines(linewidth=0.4)
fig.colorbar(pm, ax=ax, shrink=0.6, label="surface current speed [m s$^{-1}$]")
ax.set_title("GFDL-CM4 surface current speed (xgcm interp to cell centre)")
plt.show()
Relative vorticity (curl of the velocity)¶
Differencing the velocity components lands the relative vorticity
$\zeta = \partial_x v - \partial_y u$ on the cell corner. Near the North Pole
this crosses the bipolar north fold — handled properly in
06_tripolar_fold.ipynb — so for this simple global
map we just look away from the Arctic seam.
zeta = (grid.diff(ds.vo, "X", boundary="fill")
- grid.diff(ds.uo, "Y", boundary="fill")).load()
# corner geographic coords ~ the centre coords (half-cell offset, fine for plotting)
lim = float(np.nanpercentile(np.abs(zeta), 99))
fig = plt.figure(figsize=(11, 6))
ax = plt.axes(projection=ccrs.Robinson(central_longitude=-150))
pm = ax.pcolormesh(lon, lat, zeta.values[:lat.shape[0], :lat.shape[1]],
transform=ccrs.PlateCarree(), cmap="RdBu_r", vmin=-lim, vmax=lim)
ax.coastlines(linewidth=0.4)
fig.colorbar(pm, ax=ax, shrink=0.6, label="relative vorticity [s$^{-1}$]")
ax.set_title("GFDL-CM4 surface relative vorticity")
plt.show()
See the MITgcm examples for more of what xgcm can do, the
MOM6 Analysis Cookbook
for more MOM6-specific recipes, and 06_tripolar_fold.ipynb
for fold-aware operations across the Arctic seam.