Barotropic Stream Functions
Description¶
This recipe shows how to calculate the barotropic stream function for the North Atlantic Ocean 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 barotropic stream function is routinely used to characterise the large-scale ocean circulation $\psi_{xy}$ and can be defined using the meridional velocity field at time $t$ as follows:
$$\Psi_{xy}(\lambda, \phi, t) = \int_{x_w}^{x} \int_{-H}^{\eta} v(\lambda, \phi, z, t) \ dz \ dx$$
where the meridional velocity $v(\lambda, \phi, z, t)$ is first integrated from the sea surface $\eta$ to the seafloor $-H$ before being accumulated zonally from the western ($x_w$) to the eastern ($x_e$) boundary of our model domain.
In [1]:
Copied!
# -- Import required packages -- #
import xarray as xr
from nemo_cookbook import NEMODataTree
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
xr.set_options(display_style="text")
# -- Import required packages -- #
import xarray as xr
from nemo_cookbook import NEMODataTree
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
xr.set_options(display_style="text")
Out[1]:
<xarray.core.options.set_options at 0x3286767b0>
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 JRA-55 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 [2]:
Copied!
# Define directory path to ancillary files:
domain_filepath = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain"
# Open eORCA1 NEMO model domain_cfg:
ds_domain = (xr.open_zarr(f"{domain_filepath}/domain_cfg", consolidated=True, chunks={})
.rename({'z': 'nav_lev'})
.squeeze(drop=True)
)
ds_subbasins = xr.open_zarr(f"{domain_filepath}/subbasins", consolidated=True, chunks={})
ds_domain
# Define directory path to ancillary files:
domain_filepath = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain"
# Open eORCA1 NEMO model domain_cfg:
ds_domain = (xr.open_zarr(f"{domain_filepath}/domain_cfg", consolidated=True, chunks={})
.rename({'z': 'nav_lev'})
.squeeze(drop=True)
)
ds_subbasins = xr.open_zarr(f"{domain_filepath}/subbasins", consolidated=True, chunks={})
ds_domain
Out[2]:
<xarray.Dataset> Size: 667MB
Dimensions: (y: 331, x: 360, nav_lev: 75)
Dimensions without coordinates: y, x, nav_lev
Data variables: (12/54)
e1t (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e2v (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
bottom_level (y, x) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e2t (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e2u (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
bathy_metry (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
... ...
nav_lev (nav_lev) float32 300B dask.array<chunksize=(75,), meta=np.ndarray>
mask_csundef (y, x) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
nav_lat (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
nav_lon (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
time_counter float64 8B dask.array<chunksize=(), meta=np.ndarray>
top_level (y, x) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
Attributes:
DOMAIN_number_total: 1
DOMAIN_number: 0
DOMAIN_dimensions_ids: [1, 2]
DOMAIN_size_global: [362, 332]
DOMAIN_size_local: [362, 332]
DOMAIN_position_first: [1, 1]
DOMAIN_position_last: [362, 332]
DOMAIN_halo_size_start: [0, 0]
DOMAIN_halo_size_end: [0, 0]
DOMAIN_type: BOX
history: Mon Jun 5 12:41:32 2023: ncks -A mask.nc ORCA1_...
NCO: 4.4.7
Next, we need to import the meridional velocity and vertical grid cell thicknesses stored at V-points in a single dataset.
Typically, NEMO model outputs defined on V-grid points are stored together in netCDF files. In this case, you can replace xr.merge() with a single call to xarray's open_dataset() function passing the file path to your _gridV.nc file(s).
In [3]:
Copied!
# Define directory path to model output files:
output_dir = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/V1y"
# Construct NEMO model grid dataset, including vertical grid cell thicknesses (m) and meridional velocities (m/s):
ds_gridV = xr.merge([xr.open_zarr(f"{output_dir}/{var}", consolidated=True, chunks={})[var] for var in ['e3v', 'vo']], compat="override")
ds_gridV
# Define directory path to model output files:
output_dir = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/V1y"
# Construct NEMO model grid dataset, including vertical grid cell thicknesses (m) and meridional velocities (m/s):
ds_gridV = xr.merge([xr.open_zarr(f"{output_dir}/{var}", consolidated=True, chunks={})[var] for var in ['e3v', 'vo']], compat="override")
ds_gridV
Out[3]:
<xarray.Dataset> Size: 4GB
Dimensions: (depthv: 75, y: 331, x: 360, time_counter: 49)
Coordinates:
* depthv (depthv) float32 300B 0.5058 1.556 ... 5.698e+03 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] 392B dask.array<chunksize=(1,), meta=np.ndarray>
* time_counter (time_counter) datetime64[ns] 392B 1976-07-02 ... 2024-07-02
Dimensions without coordinates: y, x
Data variables:
e3v (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>
Attributes:
cell_methods: time: mean (interval: 3600 s)
interval_operation: 3600 s
interval_write: 1 yr
long_name: V-cell thickness
online_operation: average
standard_name: cell_thickness
units: m
Creating a NEMODataTree¶
Next, let's create a NEMODataTree to store our domain and V-grid variables for the eORCA1 model.
In [4]:
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, "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")
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, "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")
nemo
Out[4]:
<xarray.DataTree>
Group: /
│ Dimensions: (time_counter: 49)
│ Coordinates:
│ time_centered (time_counter) datetime64[ns] 392B dask.array<chunksize=(1,), meta=np.ndarray>
│ * time_counter (time_counter) datetime64[ns] 392B 1976-07-02 ... 2024-07-02
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridT
│ 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) 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) 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>
│ 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) bool 9MB False False False ... False False False
│ tmaskutil (j, i) bool 119kB False False False ... False False False
│ 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) bool 9MB False False False ... False False False
│ umaskutil (j, i) bool 119kB False False False ... False False False
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridV
│ Dimensions: (time_counter: 49, k: 75, j: 331, i: 360)
│ Coordinates:
│ * depthv (k) float32 300B 0.5058 1.556 2.668 ... 5.698e+03 5.902e+03
│ time_centered (time_counter) datetime64[ns] 392B 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 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) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360
│ Data variables:
│ e3v (time_counter, k, j, i) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
│ vo (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) bool 9MB False False False ... False False False
│ vmaskutil (j, i) bool 119kB False False False ... False False False
│ Attributes:
│ cell_methods: time: mean (interval: 3600 s)
│ interval_operation: 3600 s
│ interval_write: 1 yr
│ long_name: V-cell thickness
│ online_operation: average
│ standard_name: cell_thickness
│ units: m
│ 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) bool 9MB False False False ... False False False
│ 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) bool 9MB False False False ... False False False
fmaskutil (j, i) bool 119kB False False False ... False False False
Attributes:
nftype: F
iperio: True
Calculating the Barotropic Stream Function¶
Now we have constructed our NEMODataTree, let's calculate the time-mean barotropic stream function.
$$\psi_{xy}(t) = \sum_{i^{*}=1}^{i} \sum_{k} (e_{1v} . e_{3v}(t) . v(t))$$
In this example, our eORCA1 model uses $z^{*}$ vertical coordinates, so we will use the integral() method to perform integration along the $i$ and $k$ dimensions of the NEMO model grid. The resulting 2-dimensional barotropic stream function represents the depth-integrated volume transport in the Atlantic Ocean.
In [5]:
Copied!
atlmask = ds_subbasins['atlmsk'].rename({"x":"i", "y":"j"})
bsf_atl = nemo.integral(grid="/gridV",
var="vo",
dims=["k", "i"],
cum_dims=["i"],
dir="+1",
mask=atlmask
)
bsf_atl
atlmask = ds_subbasins['atlmsk'].rename({"x":"i", "y":"j"})
bsf_atl = nemo.integral(grid="/gridV",
var="vo",
dims=["k", "i"],
cum_dims=["i"],
dir="+1",
mask=atlmask
)
bsf_atl
Out[5]:
<xarray.DataArray 'vo' (time_counter: 49, j: 331, i: 360)> Size: 47MB
dask.array<nancumsum, shape=(49, 331, 360), dtype=float64, chunksize=(1, 331, 360), chunktype=numpy.ndarray>
Coordinates:
* time_counter (time_counter) datetime64[ns] 392B 1976-07-02 ... 2024-07-02
time_centered (time_counter) datetime64[ns] 392B 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>
* j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5
* i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360
Notice that the resulting DataArray includes a dask array, so we haven't actually computed the barotropic stream function yet. To do this, we need to call the .compute() method:
In [6]:
Copied!
# Compute barotropic stream function in Sverdrups [1 Sv = 1E6 m3/s]:
bsf_atl = 1E-6 * bsf_atl.compute()
bsf_atl
# Compute barotropic stream function in Sverdrups [1 Sv = 1E6 m3/s]:
bsf_atl = 1E-6 * bsf_atl.compute()
bsf_atl
Out[6]:
<xarray.DataArray 'vo' (time_counter: 49, j: 331, i: 360)> Size: 47MB
array([[[ 0. , 0. , 0. , ..., 0. ,
0. , 0. ],
[ 0. , 0. , 0. , ..., 0. ,
0. , 0. ],
[ 0. , 0. , 0. , ..., 0. ,
0. , 0. ],
...,
[ 0. , 0. , 0. , ..., -0.93092196,
-0.93092196, -0.93092196],
[ 0. , 0. , 0. , ..., -0.93177347,
-0.93177347, -0.93177347],
[ 0. , 0. , 0. , ..., -0.93259213,
-0.93259213, -0.93259213]],
[[ 0. , 0. , 0. , ..., 0. ,
0. , 0. ],
[ 0. , 0. , 0. , ..., 0. ,
0. , 0. ],
[ 0. , 0. , 0. , ..., 0. ,
0. , 0. ],
...
[ 0. , 0. , 0. , ..., -0.87334578,
-0.87334578, -0.87334578],
[ 0. , 0. , 0. , ..., -0.87397618,
-0.87397618, -0.87397618],
[ 0. , 0. , 0. , ..., -0.87391507,
-0.87391507, -0.87391507]],
[[ 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. , 0. , ..., 0. ,
0. , 0. ],
[ 0. , 0. , 0. , ..., 0. ,
0. , 0. ]]], shape=(49, 331, 360))
Coordinates:
* time_counter (time_counter) datetime64[ns] 392B 1976-07-02 ... 2024-07-02
time_centered (time_counter) datetime64[ns] 392B 1976-07-02 ... 2024-07-02
gphiv (j, i) float64 953kB -84.16 -84.16 -84.16 ... 50.23 50.01
glamv (j, i) float64 953kB 73.5 74.5 75.5 76.5 ... 73.0 73.0 73.0
* j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5
* i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360
Visualising the time-mean barotropic stream function¶
Finally, let's visualise the results by plotting the time-mean barotropic stream function for the North Atlantic basin:
In [8]:
Copied!
# -- Create a figure with an orthographic (globe) projection -- #
fig = plt.figure(figsize=(7, 7))
proj = ccrs.Orthographic(central_longitude=-25, central_latitude=35)
ax = plt.axes(projection=proj)
ax.coastlines(resolution='110m', linewidth=0.8)
ax.add_feature(cfeature.LAND, facecolor='0.1', edgecolor='0.1', linewidth=0.2, zorder=4)
ax.gridlines(draw_labels=False, dms=True, x_inline=False, y_inline=False)
# Plot eORCA1 JRA55v1 time-mean barotropic stream function:
c_parent = ax.pcolormesh(bsf_atl.glamv, bsf_atl.gphiv, bsf_atl.mean(dim='time_counter'),
transform=ccrs.PlateCarree(),
cmap="RdBu_r", shading="auto",
vmin=-30, vmax=30,
zorder=1)
# Plot eORCA1 JRA55v1 time-mean barotropic stream function contours:
plt.contour(bsf_atl.glamv, bsf_atl.gphiv, bsf_atl.mean(dim='time_counter'),
levels=8, colors='k', linewidths=1,
transform=ccrs.PlateCarree(),
zorder=2)
# Add colorbar with label:
cb = plt.colorbar(c_parent, orientation="horizontal", pad=0.05, shrink=0.7)
cb.set_label("Barotropic Stream Function (Sv)", fontsize=11)
# -- Create a figure with an orthographic (globe) projection -- #
fig = plt.figure(figsize=(7, 7))
proj = ccrs.Orthographic(central_longitude=-25, central_latitude=35)
ax = plt.axes(projection=proj)
ax.coastlines(resolution='110m', linewidth=0.8)
ax.add_feature(cfeature.LAND, facecolor='0.1', edgecolor='0.1', linewidth=0.2, zorder=4)
ax.gridlines(draw_labels=False, dms=True, x_inline=False, y_inline=False)
# Plot eORCA1 JRA55v1 time-mean barotropic stream function:
c_parent = ax.pcolormesh(bsf_atl.glamv, bsf_atl.gphiv, bsf_atl.mean(dim='time_counter'),
transform=ccrs.PlateCarree(),
cmap="RdBu_r", shading="auto",
vmin=-30, vmax=30,
zorder=1)
# Plot eORCA1 JRA55v1 time-mean barotropic stream function contours:
plt.contour(bsf_atl.glamv, bsf_atl.gphiv, bsf_atl.mean(dim='time_counter'),
levels=8, colors='k', linewidths=1,
transform=ccrs.PlateCarree(),
zorder=2)
# Add colorbar with label:
cb = plt.colorbar(c_parent, orientation="horizontal", pad=0.05, shrink=0.7)
cb.set_label("Barotropic Stream Function (Sv)", fontsize=11)
