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:
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.
# 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
<xarray.Dataset> Size: 2MB
Dimensions: (time: 1, Z: 15, YC: 40, XG: 90, YG: 40, XC: 90, Zl: 15)
Coordinates: (12/25)
* time (time) int64 8B 39600
iter (time) int64 8B ...
* Z (Z) float32 60B -25.0 -85.0 -170.0 ... -4.19e+03 -4.855e+03
PHrefC (Z) float32 60B ...
drF (Z) float32 60B ...
* YC (YC) float32 160B -78.0 -74.0 -70.0 -66.0 ... 66.0 70.0 74.0 78.0
... ...
hFacS (Z, YG, XC) float32 216kB ...
Depth (YC, XC) float32 14kB ...
dxG (YG, XC) float32 14kB ...
rAs (YG, XC) float32 14kB ...
maskS (Z, YG, XC) bool 54kB ...
* Zl (Zl) float32 60B 0.0 -50.0 -120.0 ... -3.28e+03 -3.87e+03 -4.51e+03
Data variables:
UVEL (time, Z, YC, XG) float32 216kB ...
VVEL (time, Z, YG, XC) float32 216kB ...
WVEL (time, Zl, YC, XC) float32 216kB ...
SALT (time, Z, YC, XC) float32 216kB ...
THETA (time, Z, YC, XC) float32 216kB ...
PH (time, Z, YC, XC) float32 216kB ...
Eta (time, YC, XC) float32 14kB ...
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.)
grid = xgcm.Grid(ds, periodic=['X', 'Y'])
grid
<xgcm.Grid> X Axis (periodic, boundary='periodic'): * center XC --> left * left XG --> center T Axis (periodic, boundary='periodic'): * center time Y Axis (periodic, boundary='periodic'): * center YC --> left * left YG --> center Z Axis (periodic, boundary='periodic'): * center Z --> left * left Zl --> 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:
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.
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.
div_uv = (grid.diff(u_transport, 'X') + grid.diff(v_transport, 'Y')) / ds.rA
div_uv.dims
('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.
div_uv[0, 0].plot()
<matplotlib.collections.QuadMesh at 0x178043140>
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
zeta = (-grid.diff(ds.UVEL * ds.dxC, 'Y') + grid.diff(ds.VVEL * ds.dyC, 'X'))/ds.rAz
zeta
<xarray.DataArray (time: 1, Z: 15, YG: 40, XG: 90)> Size: 216kB
array([[[[-1.48332937e-08, -1.02609041e-08, -6.38644559e-09, ...,
-2.02205257e-08, -2.36162538e-08, -2.21622898e-08],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 3.81209766e-08, 3.97684836e-08, 4.13142871e-08, ...,
1.29681979e-07, 3.34656640e-08, 3.61909045e-08],
...,
[ 4.37701573e-08, 9.16013079e-08, -2.54561030e-08, ...,
1.26603963e-08, -6.90892721e-09, 5.01466104e-08],
[ 6.99142433e-08, 3.79757061e-08, 6.82349199e-09, ...,
6.06992998e-08, 6.45270646e-08, 4.38413430e-08],
[ 5.32967235e-08, 3.35858523e-08, 7.44277244e-08, ...,
-4.63408085e-08, 1.05465695e-07, -3.87492776e-08]],
[[-5.91518168e-09, -5.92747895e-09, -4.68231898e-09, ...,
4.79974616e-10, -5.15723109e-09, -8.35224956e-09],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 2.93285876e-08, 2.88299127e-08, 2.76142540e-08, ...,
4.49574458e-08, 2.74702163e-08, 2.87617432e-08],
...
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, ...,
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, ...,
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]]]],
shape=(1, 15, 40, 90), dtype=float32)
Coordinates:
* time (time) int64 8B 39600
* Z (Z) float32 60B -25.0 -85.0 -170.0 ... -4.19e+03 -4.855e+03
* YG (YG) float32 160B -80.0 -76.0 -72.0 -68.0 ... 64.0 68.0 72.0 76.0
* XG (XG) float32 360B 0.0 4.0 8.0 12.0 16.0 ... 344.0 348.0 352.0 356.0
rAz (YG, XG) float32 14kB ...
Attributes:
standard_name: cell_area_at_f_location
long_name: cell area
units: m
coordinate: YG XG...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:
zeta_bt = (zeta * ds.drF).sum(dim='Z')
zeta_bt.plot(vmax=2e-4)
<matplotlib.collections.QuadMesh at 0x178331fd0>
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.
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)
<matplotlib.collections.QuadMesh at 0x17864a450>
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 \ . $$
strain = (grid.diff(ds.UVEL * ds.dyG, 'X') - grid.diff(ds.VVEL * ds.dxG, 'Y')) / ds.rA
strain[0,0].plot()
<matplotlib.collections.QuadMesh at 0x1786d1f70>
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:
psi = grid.cumsum(-u_bt * ds.dyG, 'Y', boundary='fill')
psi
<xarray.DataArray (time: 1, YG: 40, XG: 90)> Size: 14kB
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]]],
shape=(1, 40, 90), dtype=float32)
Coordinates:
* time (time) int64 8B 39600
* YG (YG) float32 160B -80.0 -76.0 -72.0 -68.0 ... 64.0 68.0 72.0 76.0
* XG (XG) float32 360B 0.0 4.0 8.0 12.0 16.0 ... 344.0 348.0 352.0 356.0
Attributes:
mate: VVEL
units: m
coordinate: YC XGWe 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.
(psi[0] / 1e6).plot.contourf(levels=np.arange(-160, 40, 5))
<matplotlib.contour.QuadContourSet at 0x17834df10>
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!
maskZ = grid.interp(ds.hFacS, 'X')
(psi[0] / 1e6).where(maskZ[0]).plot.contourf(levels=np.arange(-160, 40, 5))
<matplotlib.contour.QuadContourSet at 0x1788c5b80>
Kinetic Energy¶
Finally, we plot the kinetic energy $1/2 (u^2 + v^2)$ by interpoloting both quantities the cell center point.
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()
<matplotlib.collections.QuadMesh at 0x178a02450>
