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: degCA first look, in nominal coordinates¶
The quickest way to see the data is to hand the array to xarray and let it use the
1-D x and y coordinates the file ships with. No projection, no map. It is worth
doing first: the land mask, the equator-to-pole temperature gradient and the shape
of the array are all visible straight away.
Those 1-D coordinates are only nominal, though. CMIP6 stores them for convenience and they approximate where the cells are. Away from the Arctic they are close enough to longitude and latitude to read the plot by; near the top of the grid they are not. This is a tripolar grid, and its northernmost row is the seam joining two poles that sit over Arctic land.
The top row has a single nominal y of 89.95, as though it were a line of constant
latitude at 89.95°N. Its true 2-D lat runs from 64.11° to 89.94°, since the row
actually tracks down to Canada and Siberia and back over the pole.
For a map we therefore need the 2-D curvilinear lon/lat arrays and a
projection.
thetao.plot(size=4.5, aspect=1.6)
<matplotlib.collections.QuadMesh at 0x169f58e30>
Defining 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"}},
padding={"X": "periodic", "Y": "extend"},
autoparse_metadata=False,
)
grid
<xgcm.Grid> X Axis (periodic, padding='periodic'): * center x_c --> right * right x_f --> center Y Axis (not periodic, padding='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. The X difference
wraps around the periodic seam and the Y difference uses the 'extend' boundary
we gave the grid, so neither call needs a padding of its own.
The 'extend' is not right at the top of the domain. This is a tripolar grid, so
its northernmost row is the bipolar fold, where $\partial_y u$ reaches across the
pole rather than off the edge of the grid. For this global map we leave it alone
and do not look closely at the Arctic; the
north fold example picks the problem up and fixes it.
zeta = (grid.diff(ds.vo, "X") - grid.diff(ds.uo, "Y")).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 MOM6 Analysis Cookbook and the COSIMA cosima-recipes for more MOM6-specific recipes.