Computing Meridional Overturning

Meriodional Overturning Circulation is a key diagnostic for ocean circulation. This notebook demonstrate how to compute various flavors of it using the xoverturning package.

import xarray as xr # requires >= 0.15.1
import numpy as np
import warnings
%matplotlib inline

Load the sample dataset and geolon/geolat from static file (without holes on eliminated processors). Note that this is working for both symetric and non-symetric grids but the dataset and grid file must be consistent.

dataurl = ''

ds = xr.open_dataset(f'{dataurl}/ocean_monthly_z.200301-200712.nc4',
                     chunks={'time':1, 'z_l': 1}, engine='pydap')
dsgrid = xr.open_dataset('./data/')

A first overturning computation

from xoverturning import calcmoc

This computes the meridional overturning circulation (moc) on the whole ocean:

moc = calcmoc(ds, dsgrid=dsgrid)
generating basin codes

The result is given at the vertical interfaces since the streamfunction is the result of the integration of transports over layers and on the q-point because V-velocities are located on (xh, yq) and integrated over the x-axis.

<xarray.DataArray 'psi' (time: 60, z_i: 36, yq: 576)>
dask.array<truediv, shape=(60, 36, 576), dtype=float32, chunksize=(1, 36, 576), chunktype=numpy.ndarray>
  * yq       (yq) float64 -77.82 -77.63 -77.45 -77.26 ... 89.37 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
  * time     (time) object 2003-01-16 12:00:00 ... 2007-12-16 12:00:00
moc.sel(time='2003-01').plot(vmin=-50,vmax=50, yincrease=False, cmap='bwr')
<matplotlib.collections.QuadMesh at 0x2ad48499b3a0>

MOC and Basins

xoverturning can compute the MOC over a specified basin. Default is basin='global', but pre-defined available options are basin='atl-arc' and basin='indopac'. If the grid file contain a basin DataArray produced by cmip_basins, it will use it in combination with predefined basin codes to create the basin mask. If the array does not exist, xoverturning will generate it on-the-fly using cmip_basin code. It is also possible to use an arbitrary list of basin codes as argument to basin, which could be particularly handy if you have custom build basin codes.

# atlantic
amoc = calcmoc(ds, dsgrid=dsgrid, basin='atl-arc')
generating basin codes
amoc.sel(time='2003-01').plot(vmin=-30,vmax=30, yincrease=False, cmap='bwr')
<matplotlib.collections.QuadMesh at 0x2ad484d80ac0>
# indopacific
pmoc = calcmoc(ds, dsgrid=dsgrid, basin='indopac')
generating basin codes
pmoc.sel(time='2003-01').plot(vmin=-50,vmax=50, yincrease=False, cmap='bwr')
<matplotlib.collections.QuadMesh at 0x2ad4851f9fd0>

Masking (only in geopotential coordinates)

You can also mask out the ocean floor to make it more plot-friendly with mask_output.

amoc = calcmoc(ds, dsgrid=dsgrid, basin='atl-arc', mask_output=True)
amoc.sel(time='2003-01').plot(vmin=-30,vmax=30, yincrease=False, cmap='bwr',
                              subplot_kws={'facecolor': 'grey'})
generating basin codes
<matplotlib.collections.QuadMesh at 0x2ad484dde040>

Density coordinates

Because xoverturning is written in a non-specific vertical coordinate system, it can also compute the MOC in other coordinates system, such as the rho2 output by using the vertical='rho2' option.

pproot = '/archive/Raphael.Dussin/FMS2019.01.03_devgfdl_20201120/'
run = 'CM4_piControl_c96_OM4p25_half_kdadd'
ppdir = f'{pproot}/{run}/gfdl.ncrc4-intel18-prod-openmp/pp/ocean_annual_rho2'

ds_rho2 = xr.open_mfdataset([f"{ppdir}/ts/annual/10yr/",

dsgrid_rho2 = xr.open_dataset(f"{ppdir}/")
amoc_rho2 = calcmoc(ds_rho2, dsgrid=dsgrid_rho2, vertical='rho2', basin='atl-arc')
generating basin codes
<xarray.DataArray 'psi' (time: 10, rho2_i: 36, yq: 1081)>
dask.array<truediv, shape=(10, 36, 1081), dtype=float32, chunksize=(10, 36, 1081), chunktype=numpy.ndarray>
  * rho2_i   (rho2_i) float64 999.5 1.028e+03 1.029e+03 ... 1.037e+03 1.038e+03
  * time     (time) object 0021-07-02 12:00:00 ... 0030-07-02 12:00:00
  * yq       (yq) float64 -80.43 -80.35 -80.27 -80.19 ... 89.68 89.78 89.89 90.0
amoc_rho2.sel(time='0022-07').plot(vmin=-30,vmax=30, yincrease=False, cmap='bwr')
<matplotlib.collections.QuadMesh at 0x2ad484a5d430>

Derived quantities

Since the output of xoverturning is a xarray.DataArray, computing derived quantities is a breeze thanks to xarray built-in functions:

maxmoc_atl = amoc.max(dim=['yq', 'z_i'])
[<matplotlib.lines.Line2D at 0x2ad4868e1ee0>]
maxmoc_atl = amoc.sel(yq=slice(20.0, 80.0), z_i=slice(500.0, 2500.0)).max(dim=['yq', 'z_i'])
[<matplotlib.lines.Line2D at 0x2ad484d68340>]

Misc options

xoverturning also have more specialized options such as remove_hml that is used to remove the signal ov vhml and a rotate option attempting to rotate to true north. This latest option is still in progress and will be updated.