"""Core metrics to evaluate physical consistency of simulations and projections."""
# Adapted from https://github.com/Etienne-Meunier/Metrics-Ocean
# Original author: Etienne Meunier
import xarray
[docs]
def check_density(density: xarray.DataArray, epsilon: float = 1e-5) -> xarray.DataArray:
"""
Check density monotonicity violations in ocean data.
Calculate the proportion of grid points at each time step where
the density profile violates the monotonic-with-depth constraint.
Parameters
----------
density : xarray.DataArray
A 4D DataArray with dimensions including 'time_counter', 'depth',
'y', and 'x'. Represents the density field over time and space.
epsilon : float, optional
A small threshold used to determine significant non-monotonicity.
Default is 1e-5.
Returns
-------
xarray.DataArray
1D DataArray (over time_counter), with values in [0,1], representing
the proportion of grid points per time step where density increases
with depth beyond the epsilon threshold.
"""
density = density.where(density != 0)
diff = density - density.shift(depth=-1)
bad_prop = (diff > epsilon).mean(dim=["depth", "y", "x"])
return bad_prop
[docs]
def temperature_500m_30NS_metric(
temperature: xarray.DataArray, file_mask: xarray.Dataset
) -> xarray.DataArray:
"""
Metric Extraction Function.
Compute the average Temperature at 500m depth between 30N and 30S. Unit: °C.
Parameters
----------
temperature : xarray.DataArray
Temperature data with dimensions (time_counter, depth, y, x) and 2D spatial
coordinates (nav_lat, nav_lon).
file_mask : xarray.Dataset
Dataset containing the grid mask variables (e1t, e2t, e3t_0, tmask).
Returns
-------
xarray.DataArray
1D DataArray (over time_counter), representing the area-weighted mean
temperature at 500m between 30N and 30S.
"""
# Taking Temperature At 500m depth and between 30N and 30S.
DEPTH = 500 # m : upper bound for deep water box
LAT_BOUND = 30 # degrees N/S : tropical boundary
t500_30NS = temperature.sel(depth=DEPTH, method="nearest").where(
abs(temperature.nav_lat) < LAT_BOUND,
drop=False,
)
# Computing Area Weights from Mask over 30N-30S latitude zone and @500m depth
e1t = file_mask.e1t.squeeze()
e2t = file_mask.e2t.squeeze()
e3t = file_mask.e3t_0.squeeze()
tmask = file_mask.tmask.squeeze()
vol_500m_30NS = (
e1t
* e2t
* e3t
* tmask.sel(depth=DEPTH, method="nearest").where(
abs(temperature.nav_lat) < LAT_BOUND,
drop=False,
)
)
# Returning Average Temperature at 500m depth as a numpy scalar
return (t500_30NS * vol_500m_30NS).sum(dim=["y", "x", "depth"]) / vol_500m_30NS.sum(
dim=["y", "x", "depth"]
)
[docs]
def temperature_BWbox_metric(
temperature: xarray.DataArray, file_mask: xarray.Dataset
) -> xarray.DataArray:
"""
Metric Extraction in a "Bottom Water" box.
Average Temperature in a U-shaped "Bottom Water" box corresponding to waters below
3000m or beyond 30 degrees of latitude North and South.
.. code-block:: text
________________________________________________ _Surface
| . . . . |__________________________| . . . . |_500m
| . . . . | | . . . . |
| . . . . | Deep Water | . . . . |
| . . . . |__________________________| . . . . |_3000m
| . . . . . . . . Bottom Water . . . . . . . . |
|______________________________________________|_Bottom
S 30S Eq. 30N N
Figure: Schematic Representation of the Bottom Water box used in this metric.
Unit: °C
Parameters
----------
temperature : xarray.DataArray
Temperature data with dimensions (time_counter, depth, y, x) and 2D spatial
coordinates (nav_lat, nav_lon).
file_mask : xarray.Dataset
Dataset containing the grid mask variables (e1t, e2t, e3t_0, tmask).
Returns
-------
xarray.DataArray
1D DataArray (over time_counter), representing the area-weighted mean
temperature in the Bottom Water box.
"""
BD_DEPTH_MIN = 3000 # m : bottom water box lower bound
LAT_BOUND = 30 # degrees N/S : tropical boundary
condition = 1 - (temperature.depth < BD_DEPTH_MIN) * (
abs(temperature.nav_lat) < LAT_BOUND
)
t_BW = temperature.where(condition)
# Computing Area Weights from Mask over Box
e1t = file_mask.e1t.squeeze()
e2t = file_mask.e2t.squeeze()
e3t = file_mask.e3t_0.squeeze()
tmask = file_mask.tmask.squeeze()
# Return Average Temperature on Box
vol_BW = e1t * e2t * e3t * tmask.where(condition)
return (t_BW * vol_BW).sum(dim=["y", "x", "depth"]) / vol_BW.sum(
dim=["y", "x", "depth"]
)
[docs]
def temperature_DWbox_metric(
temperature: xarray.DataArray, file_mask: xarray.Dataset
) -> xarray.DataArray:
"""
Metric Extraction in a "Deep Water" box.
Average Temperature in a "Deep Water" box corresponding to waters between 500m
and 3000m depth and 30°N and 30°S.
.. code-block:: text
________________________________________________ _Surface
| |__________________________| |_500m
| | . . . . . . . . . . . . .| |
| | . . . .Deep Water . . . .| |
| |__________________________| |_3000m
| Bottom Water |
|______________________________________________|_Bottom
S 30S Eq. 30N N
Figure: Schematic Representation of the Deep Water box used in this metric.
Unit: °C
Parameters
----------
temperature : xarray.DataArray
Temperature data with dimensions (time_counter, depth, y, x) and 2D spatial
coordinates (nav_lat, nav_lon).
file_mask : xarray.Dataset
Dataset containing the grid mask variables (e1t, e2t, e3t_0, tmask).
Returns
-------
xarray.DataArray
1D DataArray (over time_counter), representing the area-weighted mean
temperature in the Deep Water box.
"""
DW_DEPTH_MID = 1750 # m : deep water mid depth
DW_DEPTH_HW = 1250 # m : half-width around mid depth
LAT_BOUND = 30 # degrees N/S : tropical boundary
e1t = file_mask.e1t.squeeze()
e2t = file_mask.e2t.squeeze()
e3t = file_mask.e3t_0.squeeze()
tmask = file_mask.tmask.squeeze()
condition_HW = (abs(temperature.depth - DW_DEPTH_MID) < DW_DEPTH_HW) * (
abs(temperature.nav_lat) < LAT_BOUND
)
t_DW = temperature.where(condition_HW)
# Computing Volume Weights from Mask over box
vol_DW = e1t * e2t * e3t * tmask.where(condition_HW)
# Returning Average Temperature on box
return (t_DW * vol_DW).sum(dim=["y", "x", "depth"]) / vol_DW.sum(
dim=["y", "x", "depth"]
)
[docs]
def ACC_Drake_metric(
uo: xarray.DataArray, file_mask: xarray.Dataset
) -> xarray.DataArray:
"""
Metric Extraction in the Drake Passage.
Antarctic Circumpolar Current Transport at the DINO equivalent of the Drake Passage
at (x=0).
Version 1 of ACC metric: Computes the flux assuming rigid lid (no sea surface height
variations), thus using the original e3u_0 variable from the mask.
Unit: Sv
Parameters
----------
uo : xarray.DataArray
Zonal velocity data for each point with dimension (t, depth, y, x)
and 2D spatial coordinates (nav_lat, nav_lon).
file_mask : xarray.Dataset
Dataset containing the grid mask variables (e2u, e3u_0, umask).
Returns
-------
xarray.DataArray
1D DataArray (over time_counter), representing the total transport across Drake
Passage (in Sv).
"""
umask_Drake = file_mask.umask.isel(x=0).squeeze()
e3u = file_mask.e3u_0.squeeze()
e2u = file_mask.e2u.squeeze()
# Masking the variables onto the Drake Passage
u_masked = uo.isel(x=0) * umask_Drake
e3u_masked = e3u.isel(x=0) * umask_Drake
e2u_masked = e2u.isel(x=0) * umask_Drake
# Multiplying zonal velocity by the sectional areas (e2u*e3u)
ubar = u_masked * e3u_masked
flux = (e2u_masked * ubar).sum(dim=["y", "depth"])
# Returning Total Transport across Drake passage as a numpy scalar (unit: Sv)
return flux / 1e6
[docs]
def ACC_Drake_metric_2(
uo: xarray.DataArray, ssh: xarray.DataArray, file_mask: xarray.Dataset
) -> xarray.DataArray:
"""
Metric Extraction of the Drake Passage.
Antarctic Circumpolar Current Transport at the DINO equivalent of the Drake Passage
at (x=0).
Version 2 of ACC metric: Computes the flux assuming varying ssh, thus needing to
recompute e3u variable from e3u_0.
Unit: Sv
Parameters
----------
uo : xarray.DataArray
Zonal velocity data with dimension (t, depth, y, x) and spatial coordinates
(nav_lat, nav_lon). Zonal velocity value for each point.
ssh : xarray.DataArray
Sea Surface Height data (t, y, x) and spatial coordinates (nav_lat, nav_lon).
file_mask : xarray.Dataset
Dataset containing the grid mask variables (e2u, e3u_0, umask).
Returns
-------
xarray.DataArray
1D DataArray (over time_counter), representing the total transport across Drake
Passage (in Sv).
"""
umask_Drake = file_mask.umask.isel(x=0).squeeze()
e3u_0 = file_mask.e3u_0.squeeze()
e2u = file_mask.e2u.squeeze()
# Recomputing e3u, using ssh to refactor the original e3u_0 cell heights)
ssh_u = (ssh + ssh.roll(x=-1)) / 2
bathy_u = e3u_0.sum(dim="depth")
ssumask = umask_Drake[:, 0]
e3u = e3u_0 * (1 + ssh_u * ssumask / (bathy_u + 1 - ssumask))
# Masking the variables onto the Drake Passage
u_masked = uo.isel(x=0) * umask_Drake
e3u_masked = e3u.isel(x=0) * umask_Drake
e2u_masked = e2u.isel(x=0) * umask_Drake
# Multiplying zonal velocity by the sectional areas (e2u*e3u)
ubar = u_masked * e3u_masked
flux = (e2u_masked * ubar).sum(dim=["y", "depth"])
# Return Total Transport across Drake passage as a numpy scalar (unit: Sv)
return flux / 1e6
[docs]
def NASTG_BSF_max(
vo: xarray.DataArray, ssh: xarray.DataArray, file_mask: xarray.Dataset
) -> xarray.DataArray:
"""
Metric Extraction of the North-Atlantic SubTropical Gyre (NASTG).
Intensity of the North-Atlantic SubTropical Gyre (NASTG) containing the Gulf-Stream
current computed from the local maximum of the Barotropic Stream Function (BSF)
Unit: Sv
Parameters
----------
vo : xarray.DataArray
Meridional velocity data with dimension (t, depth, y, x) and coordinate
(nav_lat, nav_lon). Meridional velocity value for each point.
ssh : xarray.DataArray
Sea Surface Height data with dimensions (t, y, x) and spatial coordinates
(nav_lat, nav_lon).
file_mask : xarray.Dataset
Dataset containing the grid mask variables (e1v, e2v, vmask).
Returns
-------
xarray.DataArray
1D DataArray (over time_counter), representing the maximum Barotropic Stream
Function in the NASTG region (in Sv).
"""
e3v_0 = file_mask.e3v_0.squeeze()
e1v = file_mask.e1v.squeeze()
vmask = file_mask.vmask.squeeze()
# Updating e3v from e3v_0 and SSH
ssh_v = (ssh + ssh.roll(y=-1)) / 2
bathy_v = e3v_0.sum(dim="depth")
ssvmask = vmask.isel(depth=0)
e3v = e3v_0 * (1 + ssh_v * ssvmask / (bathy_v + 1 - ssvmask))
# Integrating Meridional Transport (e3v*e1v*vo) along depth and X from Western
# boundary eastward
# to get Barotropic Stream Function with the "American continent" as reference point
# (BSF=0)
V = (vo * e3v).sum(dim="depth") # == "Barotropic Velocity" * Bathymetry
BSF = (V * e1v * ssvmask).cumsum(
dim="x"
) / 1e6 # Integrating from the West, and converting from m³/s to Sv
# Selecting 0N-40N window where to search for the maximum, which will correspond to
# the center of rotation for the gyre
BSF_NASPG = BSF.where(abs(BSF.nav_lat - 20) < 20) # noqa: PLR2004
# Selecting the maximum value of the BSF in the selected window
# and return it as a numpy scalar
return BSF_NASPG.max(dim=["y", "x"])