Zonally Averaged Overturning Circulation

This notebook shows a simple example of calculation the zonally averaged global meridional overturning circulation - in density space.

Requirements: The conda/analysis3-20.01 (or later) module on the VDI (or your own up-to-date cookbook installation).

Firstly, load in the requisite libraries:

import cosima_cookbook as cc
import matplotlib.pyplot as plt
import numpy as np
import cmocean as cm
from dask.distributed import Client
client = Client(n_workers=4)
client

Client

Cluster

  • Workers: 4
  • Cores: 48
  • Memory: 202.49 GB

Next, choose an experiment. This can be any resolution, and can be with or without Gent-McWilliams eddy parameterisation. In this case, we are choosing to limit ourselves to just the last 60 years of the chosen simulation.

expt = '025deg_jra55v13_iaf_gmredi6'
session = cc.database.create_session()
start_time='2198-01-01'

Load up ty_trans_rho - and sum zonally. Also, if there is a ty_trans_rho_gm variable saved, assume that GM is switched on and load that as well.

psi = cc.querying.getvar(expt,'ty_trans_rho',session,start_time = start_time)
#psi = cc.get_nc_variable(expt, 'ocean.nc', 'ty_trans_rho',chunks={'potrho': None}, n=-10)
psi = psi.sum('grid_xt_ocean')

varlist = cc.querying.get_variables(session, expt)
if varlist['name'].str.contains('ty_trans_rho_gm').any():
    GM = True
    psiGM = cc.querying.getvar(expt,'ty_trans_rho_gm',session,start_time = start_time)
    psiGM = psiGM.sum('grid_xt_ocean')
else:
    GM = False

Most ACCESS-OM2 simulations save transport with units of kg/s - convert to Sv:

rho = 1025 # mean density of sea-water in kg/m^3
psi = psi / (1e6*rho) # converts kg/s to Sv
if GM:
    psiGM = psiGM / (1e6*rho)

Now, cumulatively sum the transport in the vertical. Note that in MOM5 the ty_trans_rho_GM variable is computed differently and does not require summing in the vertical. Once the calculation has been laid out, we then load the variable to force the computation to occur.

psi_avg = psi.cumsum('potrho').mean('time') - psi.sum('potrho').mean('time')
if GM:
    psi_avg = psi_avg + psiGM.mean('time')

psi_avg.load()
Show/Hide data repr Show/Hide attributes
xarray.DataArray
  • potrho: 80
  • grid_yu_ocean: 1080
  • 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 -4.3655746e-11
    array([[ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00, ...,
            -7.9421811e-03, -9.9680889e-03, -6.5573095e-03],
           [ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00, ...,
            -8.3506294e-03, -1.0240269e-02, -6.8268804e-03],
           [ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00, ...,
            -8.7391092e-03, -1.0530400e-02, -7.0132408e-03],
           ...,
           [ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00, ...,
             0.0000000e+00,  0.0000000e+00, -4.3655746e-11],
           [ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00, ...,
             0.0000000e+00,  0.0000000e+00, -4.3655746e-11],
           [ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00, ...,
             0.0000000e+00,  0.0000000e+00, -4.3655746e-11]], dtype=float32)
    • grid_yu_ocean
      (grid_yu_ocean)
      float64
      -81.02 -80.92 -80.81 ... 89.89 90.0
      long_name :
      ucell latitude
      units :
      degrees_N
      cartesian_axis :
      Y
      array([-81.024202, -80.918603, -80.813004, ...,  89.788884,  89.894483,
              90.      ])
    • potrho
      (potrho)
      float64
      1.028e+03 1.028e+03 ... 1.038e+03
      long_name :
      potential density
      units :
      kg/m^3
      cartesian_axis :
      Z
      positive :
      down
      edges :
      potrho_edges
      array([1028.0625, 1028.1875, 1028.3125, 1028.4375, 1028.5625, 1028.6875,
             1028.8125, 1028.9375, 1029.0625, 1029.1875, 1029.3125, 1029.4375,
             1029.5625, 1029.6875, 1029.8125, 1029.9375, 1030.0625, 1030.1875,
             1030.3125, 1030.4375, 1030.5625, 1030.6875, 1030.8125, 1030.9375,
             1031.0625, 1031.1875, 1031.3125, 1031.4375, 1031.5625, 1031.6875,
             1031.8125, 1031.9375, 1032.0625, 1032.1875, 1032.3125, 1032.4375,
             1032.5625, 1032.6875, 1032.8125, 1032.9375, 1033.0625, 1033.1875,
             1033.3125, 1033.4375, 1033.5625, 1033.6875, 1033.8125, 1033.9375,
             1034.0625, 1034.1875, 1034.3125, 1034.4375, 1034.5625, 1034.6875,
             1034.8125, 1034.9375, 1035.0625, 1035.1875, 1035.3125, 1035.4375,
             1035.5625, 1035.6875, 1035.8125, 1035.9375, 1036.0625, 1036.1875,
             1036.3125, 1036.4375, 1036.5625, 1036.6875, 1036.8125, 1036.9375,
             1037.0625, 1037.1875, 1037.3125, 1037.4375, 1037.5625, 1037.6875,
             1037.8125, 1037.9375])

Now we are ready to plot. We usually plot the streamfunction over a reduced range of density levels to highlight the deep ocean contribution…

plt.figure(figsize=(10, 5))
clev = np.arange(-25,25,2)
plt.contourf(psi_avg.grid_yu_ocean,psi_avg.potrho, psi_avg, cmap=cm.cm.curl, levels=clev, extend='both')
cb=plt.colorbar(orientation='vertical', shrink = 0.7)

cb.ax.set_xlabel('Sv')
plt.contour(psi_avg.grid_yu_ocean, psi_avg.potrho, psi_avg, levels=clev, colors='k', linewidths=0.25)
plt.contour(psi_avg.grid_yu_ocean, psi_avg.potrho, psi_avg, levels=[0.0,], colors='k', linewidths=0.5)
plt.gca().invert_yaxis()

plt.ylim((1037.5,1034))
plt.ylabel('Potential Density (kg m$^{-3}$)')
plt.xlabel('Latitude ($^\circ$N)')
plt.xlim([-75,80])
plt.title('Overturning in %s' % expt);
../_images/Zonally_Averaged_Global_Meridional_Overturning_Circulation_0.png

Alternatively, you may want to stretch your axes to minimise the visual impact of the surface circulation, while showing the full-depth ocean.

fig,ax = plt.subplots(1,1,figsize=(10, 5))
clev = np.arange(-25,25,2)
yticks = np.array([1030, 1032, 1033, 1034, 1035, 1036,1036.5, 1037])
scfac = 4  ## A power to set teh stretching
p1=ax.contourf(psi_avg.grid_yu_ocean,(psi_avg.potrho-1028)**scfac, psi_avg, cmap=cm.cm.curl, levels=clev, extend='both')
cb=plt.colorbar(p1,orientation='vertical', shrink = 0.7)

cb.ax.set_xlabel('Sv')
ax.contour(psi_avg.grid_yu_ocean,(psi_avg.potrho-1028)**scfac, psi_avg, levels=clev, colors='k', linewidths=0.25)
ax.contour(psi_avg.grid_yu_ocean,(psi_avg.potrho-1028)**scfac, psi_avg, levels=[0.0,], colors='k', linewidths=0.5)

ax.set_yticks((yticks-1028)**scfac)
ax.set_yticklabels(yticks)
ax.set_ylim([0.5**scfac, 9.2**scfac])
ax.invert_yaxis()
ax.set_ylabel('Potential Density (kg m$^{-3}$)')
ax.set_xlabel('Latitude ($^\circ$N)')
ax.set_xlim([-75,80])
ax.set_title('Overturning in %s' % expt);
../_images/Zonally_Averaged_Global_Meridional_Overturning_Circulation_1.png

Notes: * We have not included the submesoscale contribution to the meridional transport in these calculations, as it tends to be relatively unimportant for the deep circulation, which we where we are primarily interested. * These metrics do not use mathematically correct zonal averaging in the tripole region, north of 65°N.