Getting started with xgcm for MOM6¶
MOM6 variables are staggered according to the Arakawa C-grid
It uses a north-east index convention
center points are labelled (xh, yh) and corner points are labelled (xq, yq)
important: variables xh/yh, xq/yq that are named “nominal” longitude/latitude are not the true geographical coordinates and are not suitable for plotting (more later)
See indexing for details.
[1]:
import xarray as xr
from xgcm import Grid
import warnings
import matplotlib.pylab as plt
from cartopy import crs as ccrs
import numpy as np
[2]:
%matplotlib inline
warnings.filterwarnings("ignore")
_ = xr.set_options(display_style='text')
For this tutorial, we use sample data for the \(\frac{1}{2}^{\circ}\) global model OM4p05 hosted on a GFDL thredds server:
[3]:
dataurl = 'http://35.188.34.63:8080/thredds/dodsC/OM4p5/'
ds = xr.open_dataset(f'{dataurl}/ocean_monthly_z.200301-200712.nc4',
chunks={'time':1, 'z_l': 1}, drop_variables=['average_DT',
'average_T1',
'average_T2'],
engine='pydap')
[4]:
ds
[4]:
<xarray.Dataset>
Dimensions: (nv: 2, time: 60, xh: 720, xq: 720, yh: 576, yq: 576, z_i: 36, z_l: 35)
Coordinates:
* nv (nv) float64 1.0 2.0
* xh (xh) float64 -299.8 -299.2 -298.8 -298.2 ... 58.75 59.25 59.75
* xq (xq) float64 -299.5 -299.0 -298.5 -298.0 ... 59.0 59.5 60.0
* yh (yh) float64 -77.91 -77.72 -77.54 -77.36 ... 89.47 89.68 89.89
* yq (yq) float64 -77.82 -77.63 -77.45 -77.26 ... 89.58 89.79 90.0
* z_i (z_i) float64 0.0 5.0 15.0 25.0 ... 5.75e+03 6.25e+03 6.75e+03
* z_l (z_l) float64 2.5 10.0 20.0 32.5 ... 5.5e+03 6e+03 6.5e+03
* time (time) object 2003-01-16 12:00:00 ... 2007-12-16 12:00:00
Data variables:
Coriolis (yq, xq) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
areacello (yh, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
areacello_bu (yq, xq) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
areacello_cu (yh, xq) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
areacello_cv (yq, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
deptho (yh, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
dxCu (yh, xq) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
dxCv (yq, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
dxt (yh, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
dyCu (yh, xq) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
dyCv (yq, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
dyt (yh, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
geolat (yh, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
geolat_c (yq, xq) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
geolat_u (yh, xq) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
geolat_v (yq, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
geolon (yh, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
geolon_c (yq, xq) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
geolon_u (yh, xq) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
geolon_v (yq, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
hfgeou (yh, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
sftof (yh, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
thkcello (z_l, yh, xh) float32 dask.array<chunksize=(1, 576, 720), meta=np.ndarray>
wet (yh, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
wet_c (yq, xq) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
wet_u (yh, xq) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
wet_v (yq, xh) float32 dask.array<chunksize=(576, 720), meta=np.ndarray>
so (time, z_l, yh, xh) float32 dask.array<chunksize=(1, 1, 576, 720), meta=np.ndarray>
time_bnds (time, nv) timedelta64[ns] dask.array<chunksize=(1, 2), meta=np.ndarray>
thetao (time, z_l, yh, xh) float32 dask.array<chunksize=(1, 1, 576, 720), meta=np.ndarray>
umo (time, z_l, yh, xq) float32 dask.array<chunksize=(1, 1, 576, 720), meta=np.ndarray>
uo (time, z_l, yh, xq) float32 dask.array<chunksize=(1, 1, 576, 720), meta=np.ndarray>
vmo (time, z_l, yq, xh) float32 dask.array<chunksize=(1, 1, 576, 720), meta=np.ndarray>
vo (time, z_l, yq, xh) float32 dask.array<chunksize=(1, 1, 576, 720), meta=np.ndarray>
volcello (time, z_l, yh, xh) float32 dask.array<chunksize=(1, 1, 576, 720), meta=np.ndarray>
zos (time, yh, xh) float32 dask.array<chunksize=(1, 576, 720), meta=np.ndarray>
Attributes:
filename: ocean_monthly.200301-200712.zos.nc
title: OM4p5_IAF_BLING_CFC_abio_csf_mle200
associated_files: areacello: 20030101.ocean_static.nc
grid_type: regular
grid_tile: N/A
external_variables: areacello
DODS_EXTRA.Unlimited_Dimension: timexgcm grid definition¶
The horizontal dimensions are a combination of (xh or xq) and (yh or yq) corresponding to the staggered point. In the vertical z_l refers to the depth of the center of the layer and z_i to the position of the interfaces, such as len(z_i) = len(z_l) +1.
the geolon/geolat family are the TRUE geographical coordinates and are the longitude/latitude you want to use to plot results. The subscript correspond to the staggered point (c: corner, u: U-point, v: V-point, no subscript: center)
the areacello family is the area of the ocean cell at various points with a slightly naming convention (bu: corner, cu: U-point, cv: V-point, no subscript: center). Warning, because of the curvilinear grid:
\[areacello \neq dxt * dyt\]the dx/dy family has the following naming convention: dx(Cu: U-point, Cv: V-point, no suffix: center)
thkcello is the layer thickness for each cell (variable). volcello is the volume of the cell, such as:
\[volcello = areacello * thkcello\]
The MOM6 output can be written in Symetric (len(Xq) = len(Xh) + 1) or Non-symetric mode (len(Xq) = len(Xh)), where X is a notation for both x and y. In Symetric mode, one would define the grid for the global as:
grid = Grid(ds, coords={'X': {'center': 'xh', 'outer': 'xq'},
'Y': {'center': 'yh', 'outer': 'yq'},
'Z': {'inner': 'z_l', 'outer': 'z_i'} }, periodic=['X'])
and in Non-symetric mode:
grid = Grid(ds, coords={'X': {'center': 'xh', 'right': 'xq'},
'Y': {'center': 'yh', 'right': 'yq'},
'Z': {'inner': 'z_l', 'outer': 'z_i'} }, periodic=['X'])
Of course, don’t forget to drop the periodic option if you’re running a regional model. Our data is written in Non-symetric mode hence:
[5]:
grid = Grid(ds, coords={'X': {'center': 'xh', 'right': 'xq'},
'Y': {'center': 'yh', 'right': 'yq'},
'Z': {'inner': 'z_l', 'outer': 'z_i'} }, periodic=['X'])
A note on geographical coordinates¶
MOM6 uses land processor elimination, which creates blank holes in the produced geolon/geolat fields. This can result in problems while plotting. It is recommended to overwrite them by the full arrays that are produced by running the model for a few steps without land processor elimination. Here we copy one of these files.
[6]:
!curl -O https://raw.githubusercontent.com/raphaeldussin/MOM6-AnalysisCookbook/master/docs/notebooks/data/ocean_grid_sym_OM4_05.nc
% Total % Received % Xferd Average Speed Time Time Time Current
Dload Upload Total Spent Left Speed
100 6512k 100 6512k 0 0 11.2M 0 --:--:-- --:--:-- --:--:-- 11.2M
[7]:
!ncdump -h ocean_grid_sym_OM4_05.nc
netcdf ocean_grid_sym_OM4_05 {
dimensions:
yh = 576 ;
xh = 720 ;
yq = 577 ;
xq = 721 ;
variables:
float geolat(yh, xh) ;
geolat:long_name = "Latitude of tracer (T) points" ;
geolat:units = "degrees_north" ;
geolat:missing_value = 1.e+20f ;
geolat:_FillValue = 1.e+20f ;
geolat:cell_methods = "time: point" ;
float geolat_c(yq, xq) ;
geolat_c:long_name = "Latitude of corner (Bu) points" ;
geolat_c:units = "degrees_north" ;
geolat_c:missing_value = 1.e+20f ;
geolat_c:_FillValue = 1.e+20f ;
geolat_c:cell_methods = "time: point" ;
geolat_c:interp_method = "none" ;
float geolon(yh, xh) ;
geolon:long_name = "Longitude of tracer (T) points" ;
geolon:units = "degrees_east" ;
geolon:missing_value = 1.e+20f ;
geolon:_FillValue = 1.e+20f ;
geolon:cell_methods = "time: point" ;
float geolon_c(yq, xq) ;
geolon_c:long_name = "Longitude of corner (Bu) points" ;
geolon_c:units = "degrees_east" ;
geolon_c:missing_value = 1.e+20f ;
geolon_c:_FillValue = 1.e+20f ;
geolon_c:cell_methods = "time: point" ;
geolon_c:interp_method = "none" ;
double xh(xh) ;
xh:long_name = "h point nominal longitude" ;
xh:units = "degrees_east" ;
xh:cartesian_axis = "X" ;
double xq(xq) ;
xq:long_name = "q point nominal longitude" ;
xq:units = "degrees_east" ;
xq:cartesian_axis = "X" ;
double yh(yh) ;
yh:long_name = "h point nominal latitude" ;
yh:units = "degrees_north" ;
yh:cartesian_axis = "Y" ;
double yq(yq) ;
yq:long_name = "q point nominal latitude" ;
yq:units = "degrees_north" ;
yq:cartesian_axis = "Y" ;
// global attributes:
:filename = "19000101.ocean_static.nc" ;
:title = "OM4_SIS2_cgrid_05" ;
:grid_type = "regular" ;
:grid_tile = "N/A" ;
:history = "Tue Mar 3 13:41:58 2020: ncks -v geolon,geolon_c,geolat,geolat_c /archive/gold/datasets/OM4_05/mosaic_ocean.v20180227.unpacked/ocean_static_sym.nc -o ocean_grid_sym_OM5_05.nc" ;
:NCO = "4.0.3" ;
}
[8]:
ocean_grid_sym = xr.open_dataset('ocean_grid_sym_OM4_05.nc')
[9]:
ocean_grid_sym
[9]:
<xarray.Dataset>
Dimensions: (xh: 720, xq: 721, yh: 576, yq: 577)
Coordinates:
* xh (xh) float64 -299.8 -299.2 -298.8 -298.2 ... 58.75 59.25 59.75
* xq (xq) float64 -300.0 -299.5 -299.0 -298.5 ... 58.5 59.0 59.5 60.0
* yh (yh) float64 -77.91 -77.72 -77.54 -77.36 ... 89.47 89.68 89.89
* yq (yq) float64 -78.0 -77.82 -77.63 -77.45 ... 89.37 89.58 89.79 90.0
Data variables:
geolat (yh, xh) float32 ...
geolat_c (yq, xq) float32 ...
geolon (yh, xh) float32 ...
geolon_c (yq, xq) float32 ...
Attributes:
filename: 19000101.ocean_static.nc
title: OM4_SIS2_cgrid_05
grid_type: regular
grid_tile: N/A
history: Tue Mar 3 13:41:58 2020: ncks -v geolon,geolon_c,geolat,geol...
NCO: 4.0.3I used here a symetric grid to highlight the differences with the non-symetric. Since MOM6 uses the north-east convention, we can obtain the non-symetric grid from the symetric by removing the first row and column in our arrays. This overwrites our “gruyere” coordinates in our Non-symetric dataset:
[10]:
ds['geolon_c'] = xr.DataArray(data=ocean_grid_sym['geolon_c'][1:,1:], dims=('yq', 'xq'))
ds['geolat_c'] = xr.DataArray(data=ocean_grid_sym['geolat_c'][1:,1:], dims=('yq', 'xq'))
ds['geolon'] = xr.DataArray(data=ocean_grid_sym['geolon'], dims=('yh', 'xh'))
ds['geolat'] = xr.DataArray(data=ocean_grid_sym['geolat'], dims=('yh', 'xh'))
Vorticity computation¶
We can borrow the expression for vorticity from the MITgcm example and adapt it for MOM6:
[11]:
vorticity = ( - grid.diff(ds.uo * ds.dxCu, 'Y', boundary='fill')
+ grid.diff(ds.vo * ds.dyCv, 'X', boundary='fill') ) / ds.areacello_bu
[12]:
vorticity
[12]:
<xarray.DataArray (time: 60, z_l: 35, yq: 576, xq: 720)> dask.array<truediv, shape=(60, 35, 576, 720), dtype=float32, chunksize=(1, 1, 575, 719), chunktype=numpy.ndarray> Coordinates: * time (time) object 2003-01-16 12:00:00 ... 2007-12-16 12:00:00 * z_l (z_l) float64 2.5 10.0 20.0 32.5 ... 5e+03 5.5e+03 6e+03 6.5e+03 * yq (yq) float64 -77.82 -77.63 -77.45 -77.26 ... 89.37 89.58 89.79 90.0 * xq (xq) float64 -299.5 -299.0 -298.5 -298.0 ... 58.5 59.0 59.5 60.0
Let’s take the surface relative vorticity at the first time record:
[13]:
data_plot = 1e5 * vorticity.isel(time=0, z_l=0)
_ = data_plot.load()
Plotting¶
Here we want to be careful and make sure we use the right set of coordinates (geolon_c/geolat_c). Since they are not present in the DataArray, we can add them easily with:
[14]:
data_plot = data_plot.assign_coords({'geolon_c': ds['geolon_c'],
'geolat_c': ds['geolat_c']})
One thing worth noting is that geolon_c is not monotonic in the uppermost row. Hence this row needs to be removed for cartopy to properly plot. Another option is to subsample x in the MOM6 supergrid, usually named ocean_hgrid.nc.
[15]:
data_plot = data_plot.isel(xq=slice(0,-1), yq=slice(0,-1))
Now let’s define a function that will produce a publication-quality plot:
[16]:
subplot_kws=dict(projection=ccrs.PlateCarree(),
facecolor='grey')
plt.figure(figsize=[12,8])
p = data_plot.plot(x='geolon_c', y='geolat_c',
vmin=-1, vmax=1,
cmap='bwr',
subplot_kws=subplot_kws,
transform=ccrs.PlateCarree(),
add_labels=False,
add_colorbar=False)
# add separate colorbar
cb = plt.colorbar(p, ticks=[-1,-0.5,0,0.5,1], shrink=0.6)
cb.ax.tick_params(labelsize=18)
# optional add grid lines
p.axes.gridlines(color='black', alpha=0.5, linestyle='--')
[16]:
<cartopy.mpl.gridliner.Gridliner at 0x7fd8892d2b80>
Please look at the MITgcm examples for more about what xgcm can do. Also for MOM6 analysis examples using xarray and its companion software, please visit the MOM6 Analysis Cookbook.
[ ]: