MITgcm Example

xgcm is developed in close coordination with the xmitgcm package. The metadata in datasets constructed by xmitgcm should always be compatible with xgcm’s expectations. xmitgcm is necessary for reading MITgcm’s binary MDS file format. However, for this example, the MDS files have already been converted and saved as netCDF.

Below are some example of how to make calculations on mitgcm-style datasets using xgcm.

First we import xarray and xgcm:

[1]:
import xarray as xr
import numpy as np
import xgcm
from matplotlib import pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = (10,6)
Matplotlib is building the font cache using fc-list. This may take a moment.

Now we open the example dataset, which is stored with the xgcm github repository in the datasets folder.

[2]:
# hack to make file name work with nbsphinx and binder
import os
fname = '../datasets/mitgcm_example_dataset_v2.nc'
if not os.path.exists(fname):
    fname = '../' + fname

ds = xr.open_dataset(fname)
ds
[2]:
<xarray.Dataset>
Dimensions:  (XC: 90, XG: 90, YC: 40, YG: 40, Z: 15, Zl: 15, time: 1)
Coordinates:
  * time     (time) timedelta64[ns] 11:00:00
    maskC    (Z, YC, XC) bool ...
    dyC      (YG, XC) float32 ...
    hFacC    (Z, YC, XC) float32 ...
    rA       (YC, XC) float32 ...
    hFacS    (Z, YG, XC) float32 ...
    Depth    (YC, XC) float32 ...
  * YG       (YG) float32 -80.0 -76.0 -72.0 -68.0 -64.0 ... 64.0 68.0 72.0 76.0
  * Z        (Z) float32 -25.0 -85.0 -170.0 -290.0 ... -3575.0 -4190.0 -4855.0
    PHrefC   (Z) float32 ...
    dyG      (YC, XG) float32 ...
    rAw      (YC, XG) float32 ...
    drF      (Z) float32 ...
  * YC       (YC) float32 -78.0 -74.0 -70.0 -66.0 -62.0 ... 66.0 70.0 74.0 78.0
    dxG      (YG, XC) float32 ...
  * XG       (XG) float32 0.0 4.0 8.0 12.0 16.0 ... 344.0 348.0 352.0 356.0
    iter     (time) int64 ...
    maskW    (Z, YC, XG) bool ...
  * Zl       (Zl) float32 0.0 -50.0 -120.0 -220.0 ... -3280.0 -3870.0 -4510.0
    rAs      (YG, XC) float32 ...
    rAz      (YG, XG) float32 ...
    maskS    (Z, YG, XC) bool ...
    dxC      (YC, XG) float32 ...
    hFacW    (Z, YC, XG) float32 ...
  * XC       (XC) float32 2.0 6.0 10.0 14.0 18.0 ... 346.0 350.0 354.0 358.0
Data variables:
    UVEL     (time, Z, YC, XG) float32 ...
    VVEL     (time, Z, YG, XC) float32 ...
    WVEL     (time, Zl, YC, XC) float32 ...
    SALT     (time, Z, YC, XC) float32 ...
    THETA    (time, Z, YC, XC) float32 ...
    PH       (time, Z, YC, XC) float32 ...
    Eta      (time, YC, XC) float32 ...
Attributes:
    Conventions:  CF-1.6
    title:        netCDF wrapper of MITgcm MDS binary data
    source:       MITgcm
    history:      Created by calling `open_mdsdataset(extra_metadata=None, ll...

Creating the grid object

Next we create a Grid object from the dataset. We need to tell xgcm that the X and Y axes are periodic. (The other axes will be assumed to be non-periodic.)

[3]:
grid = xgcm.Grid(ds, periodic=['X', 'Y'])
grid
[3]:
<xgcm.Grid>
Y Axis (periodic):
  * center   YC --> left
  * left     YG --> center
Z Axis (not periodic):
  * center   Z --> left
  * left     Zl --> center
X Axis (periodic):
  * center   XC --> left
  * left     XG --> center
T Axis (not periodic):
  * center   time

We see that xgcm identified five different axes: X (longitude), Y (latitude), Z (depth), T (time), and 1RHO (the axis generated by the output of the LAYERS package).

Velocity Gradients

The gradients of the velocity field can be decomposed as divergence, vorticity, and strain. Below we use xgcm to compute the velocity gradients of the horizontal flow.

Divergence

The divergence of the horizontal flow is is expressed as

\[\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\]

In discrete form, using MITgcm notation, the formula for the C-grid is

\[( \delta_i \Delta y_g \Delta r_f h_w u + \delta_j \Delta x_g \Delta r_f h_s v ) / A_c\]

First we calculate the volume transports in each direction:

[4]:
u_transport = ds.UVEL * ds.dyG * ds.hFacW * ds.drF
v_transport = ds.VVEL * ds.dxG * ds.hFacS * ds.drF

The u_transport DataArray is on the left point of the X axis, while the v_transport DataArray is on the left point of the Y axis.

[5]:
display(u_transport.dims)
display(v_transport.dims)
('time', 'Z', 'YC', 'XG')
('time', 'Z', 'YG', 'XC')

Now comes the xgcm magic: we take the diff along both axes and divide by the cell area element to find the divergence of the horizontal flow. Note how this new variable is at the cell center point.

[6]:
div_uv = (grid.diff(u_transport, 'X') + grid.diff(v_transport, 'Y')) / ds.rA
div_uv.dims
[6]:
('time', 'Z', 'YC', 'XC')

We plot this near the surface and observe the expected patern of divergence at the equator and in the subpolar regions, and convergence in the subtropical gyres.

[7]:
div_uv[0, 0].plot()
[7]:
<matplotlib.collections.QuadMesh at 0x7fe38a8a5128>
_images/example_mitgcm_13_1.png

Vorticity

The vertical component of the vorticity is a fundamental quantity of interest in ocean circulation theory. It is defined as

\[\zeta = - \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \ .\]

On the c-grid, a finite-volume representation is given by

\[\zeta = (- \delta_j \Delta x_c u + \delta_i \Delta y_c v ) / A_\zeta \ .\]

In xgcm, we calculate this quanity as

[8]:
zeta = (-grid.diff(ds.UVEL * ds.dxC, 'Y') + grid.diff(ds.VVEL * ds.dyC, 'X'))/ds.rAz
zeta
[8]:
<xarray.DataArray (time: 1, Z: 15, YG: 40, XG: 90)>
array([[[[-1.483329e-08, ..., -2.216229e-08],
         ...,
         [ 5.329672e-08, ..., -3.874928e-08]],

        ...,

        [[ 0.000000e+00, ...,  0.000000e+00],
         ...,
         [ 0.000000e+00, ...,  0.000000e+00]]]], dtype=float32)
Coordinates:
  * time     (time) timedelta64[ns] 11:00:00
  * Z        (Z) float32 -25.0 -85.0 -170.0 -290.0 ... -3575.0 -4190.0 -4855.0
  * YG       (YG) float32 -80.0 -76.0 -72.0 -68.0 -64.0 ... 64.0 68.0 72.0 76.0
  * XG       (XG) float32 0.0 4.0 8.0 12.0 16.0 ... 344.0 348.0 352.0 356.0
    rAz      (YG, XG) float32 ...

…which we can see is located at the YG, XG horizontal position (also commonly called the vorticity point).

We plot the vertical integral of this quantity, i.e. the barotropic vorticity:

[9]:
zeta_bt = (zeta * ds.drF).sum(dim='Z')
zeta_bt.plot(vmax=2e-4)
[9]:
<matplotlib.collections.QuadMesh at 0x7fe38a820e80>
_images/example_mitgcm_17_1.png

A different way to calculate the barotropic vorticity is to take the curl of the vertically integrated velocity. This formulation also allows us to incorporate the \(h\) factors representing partial cell thickness.

[10]:
u_bt = (ds.UVEL * ds.hFacW * ds.drF).sum(dim='Z')
v_bt = (ds.VVEL * ds.hFacS * ds.drF).sum(dim='Z')
zeta_bt_alt = (-grid.diff(u_bt * ds.dxC, 'Y') + grid.diff(v_bt * ds.dyC, 'X'))/ds.rAz
zeta_bt_alt.plot(vmax=2e-4)
[10]:
<matplotlib.collections.QuadMesh at 0x7fe38a6eda90>
_images/example_mitgcm_19_1.png

Strain

Another interesting quantity is the horizontal strain, defined as

\[s = \frac{\partial u}{\partial x} - \frac{\partial v}{\partial y} \ .\]

On the c-grid, a finite-volume representation is given by

\[s = (\delta_i \Delta y_g u - \delta_j \Delta x_g v ) / A_c \ .\]
[11]:
strain = (grid.diff(ds.UVEL * ds.dyG, 'X') - grid.diff(ds.VVEL * ds.dxG, 'Y')) / ds.rA
strain[0,0].plot()
[11]:
<matplotlib.collections.QuadMesh at 0x7fe38a629710>
_images/example_mitgcm_21_1.png

Barotropic Transport Streamfunction

We can use the barotropic velocity to calcuate the barotropic transport streamfunction, defined via

\[u_{bt} = - \frac{\partial \Psi}{\partial y} \ , \ \ v_{bt} = \frac{\partial \Psi}{\partial x} \ .\]

We calculate this by integrating \(u_{bt}\) along the Y axis using the grid object’s cumsum method:

[12]:
psi = grid.cumsum(-u_bt * ds.dyG, 'Y', boundary='fill')
psi
[12]:
<xarray.DataArray (time: 1, YG: 40, XG: 90)>
array([[[ 0.000000e+00,  0.000000e+00, ...,  0.000000e+00,  0.000000e+00],
        [-0.000000e+00, -0.000000e+00, ..., -0.000000e+00, -0.000000e+00],
        ...,
        [-1.070507e+08, -1.066121e+08, ..., -1.073358e+08, -1.071023e+08],
        [-1.059337e+08, -1.059228e+08, ..., -1.054934e+08, -1.056321e+08]]],
      dtype=float32)
Coordinates:
  * time     (time) timedelta64[ns] 11:00:00
  * YG       (YG) float32 -80.0 -76.0 -72.0 -68.0 -64.0 ... 64.0 68.0 72.0 76.0
  * XG       (XG) float32 0.0 4.0 8.0 12.0 16.0 ... 344.0 348.0 352.0 356.0

We see that xgcm automatically shifted the Y-axis position from center (YC) to left (YG) during the cumsum operation.

We convert to sverdrups and plot with a contour plot.

[13]:
(psi[0] / 1e6).plot.contourf(levels=np.arange(-160, 40, 5))
[13]:
<matplotlib.contour.QuadContourSet at 0x7fe38a5e25c0>
_images/example_mitgcm_25_1.png

This doesn’t look nice because it lacks a suitable land mask. The dataset has no mask provided for the vorticity point. But we can build one with xgcm!

[14]:
maskZ = grid.interp(ds.hFacS, 'X')
(psi[0] / 1e6).where(maskZ[0]).plot.contourf(levels=np.arange(-160, 40, 5))
[14]:
<matplotlib.contour.QuadContourSet at 0x7fe38a47def0>
_images/example_mitgcm_27_1.png

Kinetic Energy

Finally, we plot the kinetic energy \(1/2 (u^2 + v^2)\) by interpoloting both quantities the cell center point.

[15]:
ke = 0.5*(grid.interp((ds.UVEL*ds.hFacW)**2, 'X') + grid.interp((ds.VVEL*ds.hFacS)**2, 'Y'))
ke[0,0].where(ds.maskC[0]).plot()
[15]:
<matplotlib.collections.QuadMesh at 0x7fe38a339160>
_images/example_mitgcm_29_1.png