Meridional Overturning - Tracer Space
Description¶
This recipe shows how to calculate the Atlantic Meridional Overturning Circulation (AMOC) stream function in potential density-coordinates using annual-mean outputs from the National Oceanography Centre Near-Present-Day global eORCA1 configuration of NEMO forced using JRA55-do from 1976-2024.
For more details on this model configuration and the available outputs, users can explore the Near-Present-Day documentation here.
Background¶
The diapycnal overturning stream function is routinely used to characterise the strength and structure of the AMOC in density-space (e.g., $\sigma_{0}$ or $\sigma{2}$) as a function of latitude $\phi$ and can be defined at time $t$ as follows:
$$\Psi_{\sigma_{0}}(\phi, \sigma_{0}, t) = \int_{x_w}^{x_e} \int_{z(\lambda, \phi, \sigma_{0})}^{\eta} v(\lambda, \phi, z', t) \ dz' \ dx$$
where the meridional velocity $v(\lambda, \phi, z, t)$ is first accumulated vertically from the sea surface $\eta$ to a specified isopycnal depth $z(\lambda, \phi, \sigma_{0})$ (decreasing downward) before being integrated zonally between the western $x_w$ and eastern $x_e$ boundaries of the basin.
In [1]:
Copied!
# -- Import required packages -- #
import gsw
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from nemo_cookbook import NEMODataTree
xr.set_options(display_style="text")
# -- Import required packages -- #
import gsw
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from nemo_cookbook import NEMODataTree
xr.set_options(display_style="text")
Out[1]:
<xarray.core.options.set_options at 0x1751b8050>
Using Dask¶
Optional: Connect Client to Dask Local Cluster to run analysis in parallel.
Note that, although using Dask is not strictly necessary for this simple example using eORCA1, if we wanted to generalise this recipe to eORCA025 or eORCA12 outputs, using Dask would be essential to avoid unnecessary slow calculations using only a single process.
In [ ]:
Copied!
# -- Initialise Dask Local Cluster -- #
import os
import dask
from dask.distributed import Client, LocalCluster
# Update temporary directory for Dask workers:
dask.config.set({'temporary_directory': f"{os.getcwd()}/dask_tmp",
'local_directory': f"{os.getcwd()}/dask_tmp"
})
# Create Local Cluster:
cluster = LocalCluster(n_workers=4, threads_per_worker=3, memory_limit='5GB')
client = Client(cluster)
client
# -- Initialise Dask Local Cluster -- #
import os
import dask
from dask.distributed import Client, LocalCluster
# Update temporary directory for Dask workers:
dask.config.set({'temporary_directory': f"{os.getcwd()}/dask_tmp",
'local_directory': f"{os.getcwd()}/dask_tmp"
})
# Create Local Cluster:
cluster = LocalCluster(n_workers=4, threads_per_worker=3, memory_limit='5GB')
client = Client(cluster)
client
Accessing NEMO Model Data¶
Let's begin by loading the grid variables for our eORCA1 NEMO model from the JASMIN Object Store.
Alternatively, you can replace the domain_filepath below with a file path to your domain_cfg.nc file and read this with xarray's open_dataset() function.
In [3]:
Copied!
# Define directory path to ancillary files:
domain_url = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain_cfg"
# Open eORCA1 NEMO model domain_cfg:
ds_domain = xr.open_zarr(domain_url, consolidated=True, chunks={})
ds_domain
# Define directory path to ancillary files:
domain_url = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain_cfg"
# Open eORCA1 NEMO model domain_cfg:
ds_domain = xr.open_zarr(domain_url, consolidated=True, chunks={})
ds_domain
Out[3]:
<xarray.Dataset> Size: 709MB
Dimensions: (y: 331, x: 360, nav_lev: 75)
Coordinates:
* nav_lev (nav_lev) int64 600B 0 1 2 3 4 5 6 7 ... 68 69 70 71 72 73 74
* x (x) int64 3kB 0 1 2 3 4 5 6 7 ... 353 354 355 356 357 358 359
* y (y) int64 3kB 0 1 2 3 4 5 6 7 ... 324 325 326 327 328 329 330
Data variables: (12/49)
atlmsk (y, x) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
bathy_metry (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
bottom_level (y, x) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e1f (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e1t (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e1u (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
... ...
top_level (y, x) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
umask (nav_lev, y, x) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
umaskutil (y, x) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
vmask (nav_lev, y, x) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
vmaskutil (y, x) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
wmask (nav_lev, y, x) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
Attributes:
CfgName: UNKNOWN
CfgIndex: -999
Iperio: 1
Jperio: 0
NFold: 1
NFtype: F
VertCoord: zps
IsfCav: 0
file_name: mesh_mask.nc
TimeStamp: 01/03/2025 22:19:49 -0000
Next, we need to import the conservative temperature and absolute salinity stored at T-points and calculate the potential density referenced to the sea surface.
In [4]:
Copied!
# Define directory path to model output files:
gridT_url = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/T1y"
# Construct NEMO model T-grid dataset, including conservative temperature (degC) and absolute salinity (g/kg):
ds_gridT = xr.open_zarr(gridT_url, consolidated=True, chunks={}).sel(time_counter=slice("1976-01", "2023-12"))
# Calculate potential density anomaly referenced to the sea surface (kg/m3):
ds_gridT['sigma0'] = gsw.density.sigma0(CT=ds_gridT['thetao_con'], SA=ds_gridT['so_abs'])
ds_gridT['sigma0'].name = 'sigma0'
ds_gridT
# Define directory path to model output files:
gridT_url = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/T1y"
# Construct NEMO model T-grid dataset, including conservative temperature (degC) and absolute salinity (g/kg):
ds_gridT = xr.open_zarr(gridT_url, consolidated=True, chunks={}).sel(time_counter=slice("1976-01", "2023-12"))
# Calculate potential density anomaly referenced to the sea surface (kg/m3):
ds_gridT['sigma0'] = gsw.density.sigma0(CT=ds_gridT['thetao_con'], SA=ds_gridT['so_abs'])
ds_gridT['sigma0'].name = 'sigma0'
ds_gridT
Out[4]:
<xarray.Dataset> Size: 42GB
Dimensions: (time_counter: 48, y: 331, x: 360, deptht: 75,
axis_nbounds: 2)
Coordinates:
* deptht (deptht) float32 300B 0.5058 1.556 ... 5.902e+03
nav_lat (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
nav_lon (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
time_centered (time_counter) datetime64[ns] 384B dask.array<chunksize=(1,), meta=np.ndarray>
* time_counter (time_counter) datetime64[ns] 384B 1976-07-02 ... ...
Dimensions without coordinates: y, x, axis_nbounds
Data variables: (12/75)
berg_latent_heat_flux (time_counter, y, x) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
deptht_bounds (deptht, axis_nbounds) float32 600B dask.array<chunksize=(25, 2), meta=np.ndarray>
e3t (time_counter, deptht, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
empmr (time_counter, y, x) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
evs (time_counter, y, x) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
ficeberg (time_counter, y, x) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
... ...
vohfcisf (time_counter, deptht, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
vohflisf (time_counter, deptht, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
vowflisf (time_counter, deptht, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
zos (time_counter, y, x) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
zossq (time_counter, y, x) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
sigma0 (time_counter, deptht, y, x) float64 3GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
Attributes:
name: OUTPUT/eORCA1_1y_grid_T
description: ocean T grid variables
title: ocean T grid variables
Conventions: CF-1.6
timeStamp: 2024-Dec-09 10:28:01 GMT
uuid: 89824b34-5d44-4531-8344-fe1ce39fbb7b
Next, we need to import the meridional velocity and vertical grid cell thicknesses stored at V-points in a single dataset.
In [5]:
Copied!
# Define directory path to model output files:
gridV_url = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/V1y"
# Construct NEMO model V-grid dataset, including vertical grid cell thicknesses (m) and meridional velocities (m/s):
ds_gridV = xr.open_zarr(gridV_url, consolidated=True, chunks={}).sel(time_counter=slice("1976-01", "2023-12"))
ds_gridV
# Define directory path to model output files:
gridV_url = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/V1y"
# Construct NEMO model V-grid dataset, including vertical grid cell thicknesses (m) and meridional velocities (m/s):
ds_gridV = xr.open_zarr(gridV_url, consolidated=True, chunks={}).sel(time_counter=slice("1976-01", "2023-12"))
ds_gridV
Out[5]:
<xarray.Dataset> Size: 9GB
Dimensions: (depthv: 75, axis_nbounds: 2, time_counter: 48,
y: 331, x: 360)
Coordinates:
* depthv (depthv) float32 300B 0.5058 1.556 ... 5.902e+03
nav_lat (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
nav_lon (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
time_centered (time_counter) datetime64[ns] 384B dask.array<chunksize=(1,), meta=np.ndarray>
* time_counter (time_counter) datetime64[ns] 384B 1976-07-02 ... 2...
Dimensions without coordinates: axis_nbounds, y, x
Data variables: (12/13)
depthv_bounds (depthv, axis_nbounds) float32 600B dask.array<chunksize=(25, 2), meta=np.ndarray>
e3v (time_counter, depthv, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
hfy (time_counter, y, x) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
hfy_adv (time_counter, y, x) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
hfy_diff (time_counter, y, x) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
somesatr (time_counter, y, x) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
... ...
time_centered_bounds (time_counter, axis_nbounds) datetime64[ns] 768B dask.array<chunksize=(1, 2), meta=np.ndarray>
time_counter_bounds (time_counter, axis_nbounds) datetime64[ns] 768B dask.array<chunksize=(1, 2), meta=np.ndarray>
v2o (time_counter, depthv, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
vmo (time_counter, depthv, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
vo (time_counter, depthv, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
vo_eiv (time_counter, depthv, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
Attributes:
name: OUTPUT/eORCA1_1y_grid_V
description: ocean V grid variables
title: ocean V grid variables
Conventions: CF-1.6
timeStamp: 2024-Dec-09 10:18:04 GMT
uuid: b77cb501-a564-4451-8908-d66fe201e54a
Creating a NEMODataTree¶
Next, let's create a NEMODataTree to store our domain and T- & V-grid variables for the eORCA1 model.
In [6]:
Copied!
# Define dictionary of grid datasets defining eORCA1 parent model domain with no child/grand-child nests:
# Note: domain_cfg z-dimension is expected to be named 'nav_lev'.
datasets = {"parent": {"domain": ds_domain, "gridT": ds_gridT, "gridV": ds_gridV}}
# Initialise a new NEMODataTree whose parent domain is zonally periodic & north-folding on F-points:
nemo = NEMODataTree.from_datasets(datasets=datasets, iperio=True, nftype="F", read_mask=True)
nemo
# Define dictionary of grid datasets defining eORCA1 parent model domain with no child/grand-child nests:
# Note: domain_cfg z-dimension is expected to be named 'nav_lev'.
datasets = {"parent": {"domain": ds_domain, "gridT": ds_gridT, "gridV": ds_gridV}}
# Initialise a new NEMODataTree whose parent domain is zonally periodic & north-folding on F-points:
nemo = NEMODataTree.from_datasets(datasets=datasets, iperio=True, nftype="F", read_mask=True)
nemo
Out[6]:
<xarray.DataTree 'NEMO model'>
Group: /
│ Dimensions: (time_counter: 48, axis_nbounds: 2)
│ Coordinates:
│ time_centered (time_counter) datetime64[ns] 384B dask.array<chunksize=(1,), meta=np.ndarray>
│ * time_counter (time_counter) datetime64[ns] 384B 1976-07-02 ... 2...
│ Dimensions without coordinates: axis_nbounds
│ Data variables:
│ time_centered_bounds (time_counter, axis_nbounds) datetime64[ns] 768B dask.array<chunksize=(1, 2), meta=np.ndarray>
│ time_counter_bounds (time_counter, axis_nbounds) datetime64[ns] 768B dask.array<chunksize=(1, 2), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridT
│ Dimensions: (time_counter: 48, j: 331, i: 360, k: 75,
│ axis_nbounds: 2)
│ Coordinates:
│ * deptht (k) float32 300B 0.5058 1.556 ... 5.698e+03 5.902e+03
│ time_centered (time_counter) datetime64[ns] 384B dask.array<chunksize=(1,), meta=np.ndarray>
│ gphit (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamt (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) int64 600B 1 2 3 4 5 6 7 ... 69 70 71 72 73 74 75
│ * j (j) int64 3kB 1 2 3 4 5 6 ... 326 327 328 329 330 331
│ * i (i) int64 3kB 1 2 3 4 5 6 ... 355 356 357 358 359 360
│ Dimensions without coordinates: axis_nbounds
│ Data variables: (12/81)
│ berg_latent_heat_flux (time_counter, j, i) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ deptht_bounds (k, axis_nbounds) float32 600B dask.array<chunksize=(25, 2), meta=np.ndarray>
│ e3t (time_counter, k, j, i) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
│ empmr (time_counter, j, i) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ evs (time_counter, j, i) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ ficeberg (time_counter, j, i) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ ... ...
│ e1t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ top_level (j, i) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ bottom_level (j, i) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ tmask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
│ tmaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridU
│ Dimensions: (j: 331, i: 360, k: 75)
│ Coordinates:
│ gphiu (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamu (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
│ * j (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331
│ * i (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5
│ Data variables:
│ e1u (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2u (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ umask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
│ umaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridV
│ Dimensions: (k: 75, axis_nbounds: 2, time_counter: 48, j: 331,
│ i: 360)
│ Coordinates:
│ * depthv (k) float32 300B 0.5058 1.556 ... 5.698e+03 5.902e+03
│ time_centered (time_counter) datetime64[ns] 384B dask.array<chunksize=(1,), meta=np.ndarray>
│ gphiv (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamv (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) int64 600B 1 2 3 4 5 6 7 ... 69 70 71 72 73 74 75
│ * j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 329.5 330.5 331.5
│ * i (i) int64 3kB 1 2 3 4 5 6 ... 355 356 357 358 359 360
│ Dimensions without coordinates: axis_nbounds
│ Data variables: (12/17)
│ depthv_bounds (k, axis_nbounds) float32 600B dask.array<chunksize=(25, 2), meta=np.ndarray>
│ e3v (time_counter, k, j, i) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
│ hfy (time_counter, j, i) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ hfy_adv (time_counter, j, i) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ hfy_diff (time_counter, j, i) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ somesatr (time_counter, j, i) float32 23MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ ... ...
│ vo (time_counter, k, j, i) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
│ vo_eiv (time_counter, k, j, i) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
│ e1v (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2v (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ vmask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
│ vmaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ Attributes:
│ name: OUTPUT/eORCA1_1y_grid_V
│ description: ocean V grid variables
│ title: ocean V grid variables
│ Conventions: CF-1.6
│ timeStamp: 2024-Dec-09 10:18:04 GMT
│ uuid: b77cb501-a564-4451-8908-d66fe201e54a
│ nftype: F
│ iperio: True
├── Group: /gridW
│ Dimensions: (j: 331, i: 360, k: 75)
│ Coordinates:
│ gphit (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamt (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) float64 600B 0.5 1.5 2.5 3.5 4.5 ... 71.5 72.5 73.5 74.5
│ * j (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331
│ * i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360
│ Data variables:
│ e1t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ wmask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
│ wmaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
└── Group: /gridF
Dimensions: (j: 331, i: 360, k: 75)
Coordinates:
gphif (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
glamf (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
* k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
* j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5
* i (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5
Data variables:
e1f (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e2f (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
fmask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
fmaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
Attributes:
nftype: F
iperio: True
Calculating the AMOC diapycnal overturning stream function¶
Now we have constructed our NEMODataTree, let's calculate the diapycnal overturning stream function.
In [ ]:
Copied!
# Define Atlantic Ocean basin mask:
atlmask = ds_domain['atlmsk'].rename({"x":"i", "y":"j"}).astype(bool)
# Assign (i,j) coordinates of V-points:
atlmask['i'] = atlmask['i'] + 1
atlmask['j'] = atlmask['j'] + 1.5
# Add meridional volume transport [m3/s] to NEMO V-grid - meridional velocity [m/s] * area of meridional grid cell face [m2]:
nemo['gridV']['volume_transport'] = (nemo['gridV']['vo'] * nemo.cell_area(grid='gridV', dim='j'))
# Linearly interpolate potential density from NEMO T-grid to V-grid:
nemo['gridV/sigma0'] = nemo['gridT/sigma0'].interp_to(to='V')
# Define Atlantic Ocean basin mask:
atlmask = ds_domain['atlmsk'].rename({"x":"i", "y":"j"}).astype(bool)
# Assign (i,j) coordinates of V-points:
atlmask['i'] = atlmask['i'] + 1
atlmask['j'] = atlmask['j'] + 1.5
# Add meridional volume transport [m3/s] to NEMO V-grid - meridional velocity [m/s] * area of meridional grid cell face [m2]:
nemo['gridV']['volume_transport'] = (nemo['gridV']['vo'] * nemo.cell_area(grid='gridV', dim='j'))
# Linearly interpolate potential density from NEMO T-grid to V-grid:
nemo['gridV/sigma0'] = nemo['gridT/sigma0'].interp_to(to='V')
In [8]:
Copied!
# Define potential density bins [kg /m3]:
sigma0_bins = np.arange(22, 29, 0.01)
# Compute meridional volume transport in latitude-potential density coords:
vt_sigma0_atl = nemo.binned_statistic(grid="gridV",
vars=["sigma0"],
values="volume_transport",
keep_dims=["time_counter", "j"],
bins=[sigma0_bins],
statistic="nansum",
mask=atlmask
)
vt_sigma0_atl
# Define potential density bins [kg /m3]:
sigma0_bins = np.arange(22, 29, 0.01)
# Compute meridional volume transport in latitude-potential density coords:
vt_sigma0_atl = nemo.binned_statistic(grid="gridV",
vars=["sigma0"],
values="volume_transport",
keep_dims=["time_counter", "j"],
bins=[sigma0_bins],
statistic="nansum",
mask=atlmask
)
vt_sigma0_atl
Out[8]:
<xarray.DataArray 'volume_transport' (time_counter: 48, j: 331, sigma0_bins: 699)> Size: 89MB dask.array<reshape, shape=(48, 331, 699), dtype=float64, chunksize=(48, 331, 699), chunktype=numpy.ndarray> Coordinates: * time_counter (time_counter) datetime64[ns] 384B 1976-07-02 ... 2023-07-0... * j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5 * sigma0_bins (sigma0_bins) float64 6kB 22.01 22.02 22.03 ... 28.98 28.99
Notice that the resulting DataArray includes a dask array, so we haven't actually computed the diapycnal overturning yet. To do this, we need to call the .compute() method:
In [9]:
Copied!
# Compute diapycnal overturning stream function in Sverdrups [1 Sv = 1E6 m3/s]:
# Here, we accumulate diapycnal volume transports from the lightest to the densest
# isopycnal surface.
moc_sigma0_atl = 1E-6 * vt_sigma0_atl.cumsum(dim='sigma0_bins').compute()
moc_sigma0_atl.name = 'moc_sigma0_atl'
moc_sigma0_atl
# Compute diapycnal overturning stream function in Sverdrups [1 Sv = 1E6 m3/s]:
# Here, we accumulate diapycnal volume transports from the lightest to the densest
# isopycnal surface.
moc_sigma0_atl = 1E-6 * vt_sigma0_atl.cumsum(dim='sigma0_bins').compute()
moc_sigma0_atl.name = 'moc_sigma0_atl'
moc_sigma0_atl
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
Out[9]:
<xarray.DataArray 'moc_sigma0_atl' (time_counter: 48, j: 331, sigma0_bins: 699)> Size: 89MB
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],
...,
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
-9.30921961e-01, -9.30921961e-01, -9.30921961e-01],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
-9.45353998e-01, -9.45353998e-01, -9.45353998e-01],
[ 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],
...
-9.51691358e-01, -9.51691358e-01, -9.51691358e-01],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
-9.57064748e-01, -9.57064748e-01, -9.57064748e-01],
[ 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],
...,
[-1.24457492e-03, -1.24457492e-03, -1.24457492e-03, ...,
-9.02882356e-01, -9.02882356e-01, -9.02882356e-01],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
-9.11390452e-01, -9.11390452e-01, -9.11390452e-01],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00]]],
shape=(48, 331, 699))
Coordinates:
* time_counter (time_counter) datetime64[ns] 384B 1976-07-02 ... 2023-07-0...
* j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5
* sigma0_bins (sigma0_bins) float64 6kB 22.01 22.02 22.03 ... 28.98 28.99
Visualising the time-mean AMOC diapycnal overturning stream function¶
Finally, let's visualise the results by plotting the time-mean Atlantic Meridional Overturning stream function in potential density-coordinates:
In [10]:
Copied!
plt.figure(figsize=(8, 4))
# Plot time-mean diapycnal overturning stream function:
(moc_sigma0_atl
.mean(dim='time_counter')
.plot(y='sigma0_bins',
yincrease=False,
cbar_kwargs={"label": "$\\psi_{\\sigma_{0}}(\phi, \\sigma_{0})$ [Sv]"}
)
)
# Axes labels:
plt.title('eORCA1-JRA55v1: AMOC$_{\\sigma_{0}}$', fontdict={'size':12, 'weight':'bold'})
plt.xlabel('NEMO Model $j$ Coordinate', fontdict={'size':11, 'weight':'bold'})
plt.ylabel('Potential Density $\\sigma_{0}$ (kg m$^{-3}$)', fontdict={'size':11, 'weight':'bold'})
plt.figure(figsize=(8, 4))
# Plot time-mean diapycnal overturning stream function:
(moc_sigma0_atl
.mean(dim='time_counter')
.plot(y='sigma0_bins',
yincrease=False,
cbar_kwargs={"label": "$\\psi_{\\sigma_{0}}(\phi, \\sigma_{0})$ [Sv]"}
)
)
# Axes labels:
plt.title('eORCA1-JRA55v1: AMOC$_{\\sigma_{0}}$', fontdict={'size':12, 'weight':'bold'})
plt.xlabel('NEMO Model $j$ Coordinate', fontdict={'size':11, 'weight':'bold'})
plt.ylabel('Potential Density $\\sigma_{0}$ (kg m$^{-3}$)', fontdict={'size':11, 'weight':'bold'})
<>:8: SyntaxWarning: invalid escape sequence '\p'
<>:8: SyntaxWarning: invalid escape sequence '\p'
/var/folders/z2/j_dr250s42x34hk63_rp4bm80000gq/T/ipykernel_62322/944999173.py:8: SyntaxWarning: invalid escape sequence '\p'
cbar_kwargs={"label": "$\\psi_{\\sigma_{0}}(\phi, \\sigma_{0})$ [Sv]"}
Out[10]:
Text(0, 0.5, 'Potential Density $\\sigma_{0}$ (kg m$^{-3}$)')
