Meridional Overturning Stream Function - Tracer

Calculating the AMOC in latitude-tracer coordinates¶

Description¶

Recipe showing how to calculate the Atlantic Meridional Overturning 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.

In [1]:

# -- Import required Python packages -- #
import gsw
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

# -- Import nemo_cookbook tools -- #
from nemo_cookbook import compute_moc_tracer

# -- Import required Python packages -- # import gsw import numpy as np import xarray as xr import matplotlib.pyplot as plt # -- Import nemo_cookbook tools -- # from nemo_cookbook import compute_moc_tracer

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 [ ]:

# -- Initialise Dask Local Cluster -- #
import dask
from dask.distributed import Client, LocalCluster

# Update temporarty directory for Dask workers:
dask.config.set({'temporary_directory': '/home/otooth/work/Diagnostics/proj_NPD_diag/nemo_cookbook/recipes',
                 'local_directory': '/home/otooth/work/Diagnostics/proj_NPD_diag/nemo_cookbook/recipes'
                 })

# Create Local Cluster:
cluster = LocalCluster(n_workers=4, threads_per_worker=5, memory_limit='5GB')
client = Client(cluster)
client

# -- Initialise Dask Local Cluster -- # import dask from dask.distributed import Client, LocalCluster # Update temporarty directory for Dask workers: dask.config.set({'temporary_directory': '/home/otooth/work/Diagnostics/proj_NPD_diag/nemo_cookbook/recipes', 'local_directory': '/home/otooth/work/Diagnostics/proj_NPD_diag/nemo_cookbook/recipes' }) # Create Local Cluster: cluster = LocalCluster(n_workers=4, threads_per_worker=5, memory_limit='5GB') client = Client(cluster) client

Preparing 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 [ ]:

# Define directory path to ancillary files:
domain_filepath = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain"

# Open eORCA1 model grid data:
ds_domain = xr.open_zarr(domain_filepath, consolidated=True, chunks={})

# Extract zonal grid cell widths (m):
e1v = ds_domain['e1v'].squeeze()
# Extract Atlantic Ocean mask:
atl_mask = ds_domain['atlmsk'].squeeze()

# Define directory path to ancillary files: domain_filepath = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain" # Open eORCA1 model grid data: ds_domain = xr.open_zarr(domain_filepath, consolidated=True, chunks={}) # Extract zonal grid cell widths (m): e1v = ds_domain['e1v'].squeeze() # Extract Atlantic Ocean mask: atl_mask = ds_domain['atlmsk'].squeeze()

Next, we need to import the conservative temperature and absolute salinity stored at T-points

In [ ]:

# Define directory path to model output files:
output_dir = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1"

# Extract conservative temperature (C):
thetao_con = xr.open_zarr(f"{output_dir}/T1y/thetao_con", consolidated=True, chunks={})['thetao_con']
# Extract absolute salinity (g/kg):
so_abs = xr.open_zarr(f"{output_dir}/T1y/so_abs", consolidated=True, chunks={})['so_abs']

# Calculate potential density anomaly referenced to the sea surface (kg/m3):
sigma0 = gsw.density.sigma0(CT=thetao_con, SA=so_abs)
sigma0.name = 'sigma0'

# Define directory path to model output files: output_dir = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1" # Extract conservative temperature (C): thetao_con = xr.open_zarr(f"{output_dir}/T1y/thetao_con", consolidated=True, chunks={})['thetao_con'] # Extract absolute salinity (g/kg): so_abs = xr.open_zarr(f"{output_dir}/T1y/so_abs", consolidated=True, chunks={})['so_abs'] # Calculate potential density anomaly referenced to the sea surface (kg/m3): sigma0 = gsw.density.sigma0(CT=thetao_con, SA=so_abs) sigma0.name = 'sigma0'

Next, we need to import the meridional velocity field and vertical grid cell thickness stored at V-points

In [ ]:

# Extract vertical grid cell thicknesses (m):
e3v = xr.open_zarr(f"{output_dir}/V1y/e3v", consolidated=True, chunks={})['e3v']
# Extract meridional velocities (m/s):
vo = xr.open_zarr(f"{output_dir}/V1y/vo", consolidated=True, chunks={})['vo']

# Extract vertical grid cell thicknesses (m): e3v = xr.open_zarr(f"{output_dir}/V1y/e3v", consolidated=True, chunks={})['e3v'] # Extract meridional velocities (m/s): vo = xr.open_zarr(f"{output_dir}/V1y/vo", consolidated=True, chunks={})['vo']

Calculating Diapycnal Overturning Stream Function¶

Now all our input variables are ready, let's calculate the Atlantic Meridional Overturning Stream Function in density-coordinates using the compute_moc_tracer() function:

In [ ]:

# Apply the Atlantic Ocean sector mask and accumulate from the lightest to the densest isopycnal surface:
moc_sigma0_atl = compute_moc_tracer(vo=vo,
                                    e1v=e1v,
                                    e3v=e3v,
                                    tracer=sigma0,
                                    tracer_bins=np.arange(21, 29, 0.01),
                                    dir = '+1',
                                    mask=atl_mask,
                                    )

moc_sigma0_atl

# Apply the Atlantic Ocean sector mask and accumulate from the lightest to the densest isopycnal surface: moc_sigma0_atl = compute_moc_tracer(vo=vo, e1v=e1v, e3v=e3v, tracer=sigma0, tracer_bins=np.arange(21, 29, 0.01), dir = '+1', mask=atl_mask, ) moc_sigma0_atl

Out[ ]:

<xarray.DataArray 'moc_sigma0' (time_counter: 48, sigma0_bins: 799, y: 331)> Size: 102MB
dask.array<nancumsum, shape=(48, 799, 331), dtype=float64, chunksize=(48, 799, 331), chunktype=numpy.ndarray>
Coordinates:
  * time_counter  (time_counter) datetime64[ns] 384B 1976-07-02 ... 2023-07-0...
  * sigma0_bins   (sigma0_bins) object 6kB (21.0, 21.01] ... (28.980000000001...
  * y             (y) int64 3kB 0 1 2 3 4 5 6 7 ... 324 325 326 327 328 329 330
Attributes:
    units:          Sv
    long_name:      meridional overturning stream function in sigma0 coordinates
    standard_name:  moc_sigma0

xarray.DataArray

'moc_sigma0'

• time_counter: 48
• sigma0_bins: 799
• y: 331

• dask.array<chunksize=(48, 799, 331), meta=np.ndarray>

  Array Chunk Bytes 96.85 MiB 96.85 MiB Shape (48, 799, 331) (48, 799, 331) Dask graph 1 chunks in 38 graph layers Data type float64 numpy.ndarray

• Coordinates: (3)
  • time_counter
    (time_counter)
    datetime64[ns]
    1976-07-02 ... 2023-07-02T12:00:00

    axis :

    T

    bounds :

    time_counter_bounds

    long_name :

    Time axis

    standard_name :

    time

    time_origin :

    1900-01-01 00:00:00

    array(['1976-07-02T00:00:00.000000000', '1977-07-02T12:00:00.000000000',
           '1978-07-02T12:00:00.000000000', '1979-07-02T12:00:00.000000000',
           '1980-07-02T00:00:00.000000000', '1981-07-02T12:00:00.000000000',
           '1982-07-02T12:00:00.000000000', '1983-07-02T12:00:00.000000000',
           '1984-07-02T00:00:00.000000000', '1985-07-02T12:00:00.000000000',
           '1986-07-02T12:00:00.000000000', '1987-07-02T12:00:00.000000000',
           '1988-07-02T00:00:00.000000000', '1989-07-02T12:00:00.000000000',
           '1990-07-02T12:00:00.000000000', '1991-07-02T12:00:00.000000000',
           '1992-07-02T00:00:00.000000000', '1993-07-02T12:00:00.000000000',
           '1994-07-02T12:00:00.000000000', '1995-07-02T12:00:00.000000000',
           '1996-07-02T00:00:00.000000000', '1997-07-02T12:00:00.000000000',
           '1998-07-02T12:00:00.000000000', '1999-07-02T12:00:00.000000000',
           '2000-07-02T00:00:00.000000000', '2001-07-02T12:00:00.000000000',
           '2002-07-02T12:00:00.000000000', '2003-07-02T12:00:00.000000000',
           '2004-07-02T00:00:00.000000000', '2005-07-02T12:00:00.000000000',
           '2006-07-02T12:00:00.000000000', '2007-07-02T12:00:00.000000000',
           '2008-07-02T00:00:00.000000000', '2009-07-02T12:00:00.000000000',
           '2010-07-02T12:00:00.000000000', '2011-07-02T12:00:00.000000000',
           '2012-07-02T00:00:00.000000000', '2013-07-02T12:00:00.000000000',
           '2014-07-02T12:00:00.000000000', '2015-07-02T12:00:00.000000000',
           '2016-07-02T00:00:00.000000000', '2017-07-02T12:00:00.000000000',
           '2018-07-02T12:00:00.000000000', '2019-07-02T12:00:00.000000000',
           '2020-07-02T00:00:00.000000000', '2021-07-02T12:00:00.000000000',
           '2022-07-02T12:00:00.000000000', '2023-07-02T12:00:00.000000000'],
          dtype='datetime64[ns]')

  • sigma0_bins
    (sigma0_bins)
    object
    (21.0, 21.01] ... (28.9800000000...

    cell_methods :

    time: mean

    interval_operation :

    1 yr

    interval_write :

    1 yr

    long_name :

    sigma0 bins

    online_operation :

    average

    standard_name :

    sigma0

    units :

    g/kg

    array([Interval(21.0, 21.01, closed='right'),
           Interval(21.01, 21.020000000000003, closed='right'),
           Interval(21.020000000000003, 21.030000000000005, closed='right'), ...,
           Interval(28.960000000001244, 28.970000000001246, closed='right'),
           Interval(28.970000000001246, 28.980000000001247, closed='right'),
           Interval(28.980000000001247, 28.99000000000125, closed='right')],
          dtype=object)

  • y
    (y)
    int64
    0 1 2 3 4 5 ... 326 327 328 329 330

    array([  0,   1,   2, ..., 328, 329, 330])

• Indexes: (3)
  • time_counter
    PandasIndex

    PandasIndex(DatetimeIndex(['1976-07-02 00:00:00', '1977-07-02 12:00:00',
                   '1978-07-02 12:00:00', '1979-07-02 12:00:00',
                   '1980-07-02 00:00:00', '1981-07-02 12:00:00',
                   '1982-07-02 12:00:00', '1983-07-02 12:00:00',
                   '1984-07-02 00:00:00', '1985-07-02 12:00:00',
                   '1986-07-02 12:00:00', '1987-07-02 12:00:00',
                   '1988-07-02 00:00:00', '1989-07-02 12:00:00',
                   '1990-07-02 12:00:00', '1991-07-02 12:00:00',
                   '1992-07-02 00:00:00', '1993-07-02 12:00:00',
                   '1994-07-02 12:00:00', '1995-07-02 12:00:00',
                   '1996-07-02 00:00:00', '1997-07-02 12:00:00',
                   '1998-07-02 12:00:00', '1999-07-02 12:00:00',
                   '2000-07-02 00:00:00', '2001-07-02 12:00:00',
                   '2002-07-02 12:00:00', '2003-07-02 12:00:00',
                   '2004-07-02 00:00:00', '2005-07-02 12:00:00',
                   '2006-07-02 12:00:00', '2007-07-02 12:00:00',
                   '2008-07-02 00:00:00', '2009-07-02 12:00:00',
                   '2010-07-02 12:00:00', '2011-07-02 12:00:00',
                   '2012-07-02 00:00:00', '2013-07-02 12:00:00',
                   '2014-07-02 12:00:00', '2015-07-02 12:00:00',
                   '2016-07-02 00:00:00', '2017-07-02 12:00:00',
                   '2018-07-02 12:00:00', '2019-07-02 12:00:00',
                   '2020-07-02 00:00:00', '2021-07-02 12:00:00',
                   '2022-07-02 12:00:00', '2023-07-02 12:00:00'],
                  dtype='datetime64[ns]', name='time_counter', freq=None))

  • sigma0_bins
    PandasIndex

    PandasIndex(Index([                           (21.0, 21.01],
                        (21.01, 21.020000000000003],
           (21.020000000000003, 21.030000000000005],
           (21.030000000000005, 21.040000000000006],
           (21.040000000000006, 21.050000000000008],
            (21.050000000000008, 21.06000000000001],
             (21.06000000000001, 21.07000000000001],
            (21.07000000000001, 21.080000000000013],
           (21.080000000000013, 21.090000000000014],
           (21.090000000000014, 21.100000000000016],
           ...
           (28.890000000001233, 28.900000000001235],
           (28.900000000001235, 28.910000000001236],
           (28.910000000001236, 28.920000000001238],
            (28.920000000001238, 28.93000000000124],
             (28.93000000000124, 28.94000000000124],
            (28.94000000000124, 28.950000000001243],
           (28.950000000001243, 28.960000000001244],
           (28.960000000001244, 28.970000000001246],
           (28.970000000001246, 28.980000000001247],
            (28.980000000001247, 28.99000000000125]],
          dtype='object', name='sigma0_bins', length=799))

  • y
    PandasIndex

    PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
           ...
           321, 322, 323, 324, 325, 326, 327, 328, 329, 330],
          dtype='int64', name='y', length=331))

• Attributes: (3)

  units :

  Sv

  long_name :

  meridional overturning stream function in sigma0 coordinates

  standard_name :

  moc_sigma0

Notice that the output above is a dask array, so we haven't actually computed the diapycnal overturning yet. To do this, we need to call the .compute() method:

In [ ]:

moc_sigma0_atl = moc_sigma0_atl.compute()

moc_sigma0_atl = moc_sigma0_atl.compute()

Visualising Diapycnal Overturning Stream Function¶

Finally, let's visualise the results by plotting the time-mean Atlantic Meridional Overturning Stream Function in density-coordinates:

In [ ]:

moc_sigma0_atl.mean(dim='time_counter').plot(yincrease=False)

moc_sigma0_atl.mean(dim='time_counter').plot(yincrease=False)

Out[ ]:

<matplotlib.collections.QuadMesh at 0x14df5d3b2fc0>