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)

Now we open the example dataset which is stored in a zenodo archive.

[2]:

# download the data
import urllib.request
import shutil

url = 'https://zenodo.org/record/4421428/files/'
file_name = 'mitgcm_example_dataset_v2.nc'
with urllib.request.urlopen(url + file_name) as response, open(file_name, 'wb') as out_file:
    shutil.copyfileobj(response, out_file)

# open the data
ds = xr.open_dataset(file_name)
ds

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>
T Axis (not periodic, boundary=None):
  * center   time
Z Axis (not periodic, boundary=None):
  * center   Z --> left
  * left     Zl --> center
Y Axis (periodic, boundary=None):
  * center   YC --> left
  * left     YG --> center
X Axis (periodic, boundary=None):
  * center   XC --> left
  * left     XG --> center

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 0x7f93445afa00>

../../../_images/repos_xgcm_xgcm-examples_02_mitgcm_13_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 0x7f933b1c4e80>

../../../_images/repos_xgcm_xgcm-examples_02_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.00000000e+00,  0.00000000e+00,  0.00000000e+00, ...,
          0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
        [-0.00000000e+00, -0.00000000e+00, -0.00000000e+00, ...,
         -0.00000000e+00, -0.00000000e+00, -0.00000000e+00],
        [-0.00000000e+00, -0.00000000e+00, -0.00000000e+00, ...,
         -0.00000000e+00, -0.00000000e+00, -0.00000000e+00],
        ...,
        [-1.06634920e+08, -1.06666912e+08, -1.05922952e+08, ...,
         -1.07724152e+08, -1.07341128e+08, -1.06941400e+08],
        [-1.07050680e+08, -1.06612056e+08, -1.06356232e+08, ...,
         -1.07048232e+08, -1.07335792e+08, -1.07102304e+08],
        [-1.05933672e+08, -1.05922800e+08, -1.05856504e+08, ...,
         -1.05376272e+08, -1.05493376e+08, -1.05632128e+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

xarray.DataArray

time: 1
YG: 40
XG: 90

0.0 0.0 0.0 0.0 ... -105376270.0 -105493380.0 -105632130.0

array([[[ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, ...,
          0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
        [-0.00000000e+00, -0.00000000e+00, -0.00000000e+00, ...,
         -0.00000000e+00, -0.00000000e+00, -0.00000000e+00],
        [-0.00000000e+00, -0.00000000e+00, -0.00000000e+00, ...,
         -0.00000000e+00, -0.00000000e+00, -0.00000000e+00],
        ...,
        [-1.06634920e+08, -1.06666912e+08, -1.05922952e+08, ...,
         -1.07724152e+08, -1.07341128e+08, -1.06941400e+08],
        [-1.07050680e+08, -1.06612056e+08, -1.06356232e+08, ...,
         -1.07048232e+08, -1.07335792e+08, -1.07102304e+08],
        [-1.05933672e+08, -1.05922800e+08, -1.05856504e+08, ...,
         -1.05376272e+08, -1.05493376e+08, -1.05632128e+08]]],
      dtype=float32)

Coordinates: (3)

time
(time)
timedelta64[ns]
11:00:00
standard_name :
time
long_name :
Time
axis :
T
calendar :
gregorian
```
array([39600000000000], dtype='timedelta64[ns]')
```

(YG)

float32

-80.0 -76.0 -72.0 ... 72.0 76.0

standard_name :: latitude_at_f_location
long_name :: latitude
units :: degrees_north
axis :: Y
c_grid_axis_shift :: -0.5

array([-80., -76., -72., -68., -64., -60., -56., -52., -48., -44., -40., -36.,
       -32., -28., -24., -20., -16., -12.,  -8.,  -4.,   0.,   4.,   8.,  12.,
        16.,  20.,  24.,  28.,  32.,  36.,  40.,  44.,  48.,  52.,  56.,  60.,
        64.,  68.,  72.,  76.], dtype=float32)

(XG)

float32

0.0 4.0 8.0 ... 348.0 352.0 356.0

standard_name :: longitude_at_f_location
long_name :: longitude
units :: degrees_east
coordinate :: YG XG
axis :: X
c_grid_axis_shift :: -0.5

array([  0.,   4.,   8.,  12.,  16.,  20.,  24.,  28.,  32.,  36.,  40.,  44.,
        48.,  52.,  56.,  60.,  64.,  68.,  72.,  76.,  80.,  84.,  88.,  92.,
        96., 100., 104., 108., 112., 116., 120., 124., 128., 132., 136., 140.,
       144., 148., 152., 156., 160., 164., 168., 172., 176., 180., 184., 188.,
       192., 196., 200., 204., 208., 212., 216., 220., 224., 228., 232., 236.,
       240., 244., 248., 252., 256., 260., 264., 268., 272., 276., 280., 284.,
       288., 292., 296., 300., 304., 308., 312., 316., 320., 324., 328., 332.,
       336., 340., 344., 348., 352., 356.], dtype=float32)

Attributes: (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 0x7f933b0efc40>

../../../_images/repos_xgcm_xgcm-examples_02_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 0x7f933b1f72b0>

../../../_images/repos_xgcm_xgcm-examples_02_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 0x7f9344a4c5e0>

../../../_images/repos_xgcm_xgcm-examples_02_mitgcm_29_1.png

[ ]:

[ ]:

[ ]:

[ ]: