fesomp.diag

FESOM2 Diagnostics Module.

This module provides oceanographic diagnostic calculations for FESOM2 model output.

Sea Ice Diagnostics

• ice_area : Compute sea ice area (concentration × cell area)

• ice_volume : Compute sea ice volume (thickness × cell area)

• ice_extent : Compute sea ice extent (area above threshold)

Vertical Diagnostics

• hovmoller : Area-weighted mean profile over time

• volume_mean : Volume-weighted mean over depth range

• heat_content : Ocean heat content over depth range

• total_volume : Total ocean volume over depth range

Mixed Layer

• mixed_layer_depth : Mixed layer depth using threshold criterion

• mixed_layer_depth_interpolated : MLD with linear interpolation

Meridional Overturning

• moc : Meridional overturning circulation from vertical velocity

• amoc_index : AMOC index at specified latitude

• get_basin_mask : Boolean mask for named ocean basins

• list_moc_basins : List available basin names

Utilities

• hemisphere_mask : Create mask for Northern/Southern hemisphere

• get_region_mask : Create mask for named regions

• list_available_regions : List available region names

• compute_layer_thickness : Compute layer thickness from level depths

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>>
>>> # Sea ice extent
>>> sic = xr.open_dataset("a_ice.nc")["a_ice"]
>>> extent = fesomp.diag.ice_extent(
...     sic,
...     mesh.geometry.node_area[0],
...     mesh.lat,
...     hemisphere="N"
... )
>>>
>>> # Ocean heat content in top 700m
>>> temp = xr.open_dataset("temp.nc")["temp"]
>>> ohc = fesomp.diag.heat_content(
...     temp,
...     mesh.geometry.node_area[0],
...     mesh.depth_levels,
...     depth_range=(0, 700)
... )
>>>
>>> # Atlantic MOC
>>> w = xr.open_dataset("w.nc")["w"]
>>> atl_mask = fesomp.diag.get_basin_mask(
...     mesh.lon_elem, mesh.lat_elem, "Atlantic_MOC"
... )
>>> amoc, lats = fesomp.diag.moc(
...     w,
...     mesh.geometry.elem_area,
...     mesh.lat_elem,
...     mesh.depth_levels,
...     mask=atl_mask
... )

fesomp.diag.ice_area(sic, node_area, lat, hemisphere='N', mask=None)

Compute sea ice area (ice concentration × cell area).

Sea ice area is the sum of (ice concentration × cell area) over all cells. This represents the actual area covered by ice.

Parameters:

• sic (array-like) – Sea ice concentration [0-1], shape (…, n2d). Can have any leading dimensions (time, ensemble, etc.).

• node_area (np.ndarray) – Area of each node in m², shape (n2d,) or (nlev, n2d). If 2D, uses the surface level (index 0).

• lat (np.ndarray) – Latitude of nodes in degrees, shape (n2d,).

• hemisphere ({"N", "S", "both"}) – “N” for Northern Hemisphere, “S” for Southern Hemisphere, “both” for global. Default is “N”.

• mask (np.ndarray, optional) – Additional boolean mask, shape (n2d,). True = include.

Returns:

Total ice area in m², shape (…). Same type as input (numpy or xarray).

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> sic = xr.open_dataset("a_ice.nc")["a_ice"]  # shape: (time, n2d)
>>> # Compute Northern Hemisphere ice area time series
>>> nh_area = fesomp.diag.ice_area(
...     sic,
...     mesh.geometry.node_area[0],
...     mesh.lat,
...     hemisphere="N"
... )

fesomp.diag.ice_volume(siv, node_area, lat, hemisphere='N', mask=None)

Compute sea ice volume (ice thickness × cell area).

Sea ice volume is the sum of (effective ice thickness × cell area). The input can be either: - Effective ice thickness (m): already area-weighted ice thickness - Ice volume per unit area (m): same as effective thickness

Parameters:

• siv (array-like) – Sea ice effective thickness or volume per unit area in meters, shape (…, n2d). Can have any leading dimensions.

• node_area (np.ndarray) – Area of each node in m², shape (n2d,) or (nlev, n2d). If 2D, uses the surface level (index 0).

• lat (np.ndarray) – Latitude of nodes in degrees, shape (n2d,).

• hemisphere ({"N", "S", "both"}) – “N” for Northern Hemisphere, “S” for Southern Hemisphere, “both” for global. Default is “N”.

• mask (np.ndarray, optional) – Additional boolean mask, shape (n2d,). True = include.

Returns:

Total ice volume in m³, shape (…). Same type as input (numpy or xarray).

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> siv = xr.open_dataset("m_ice.nc")["m_ice"]  # shape: (time, n2d)
>>> # Compute Arctic ice volume time series
>>> arctic_vol = fesomp.diag.ice_volume(
...     siv,
...     mesh.geometry.node_area[0],
...     mesh.lat,
...     hemisphere="N"
... )

fesomp.diag.ice_extent(sic, node_area, lat, hemisphere='N', threshold=0.15, mask=None)

Compute sea ice extent (area where concentration exceeds threshold).

Sea ice extent is the total area of cells where ice concentration exceeds the threshold (typically 15%). Unlike ice area, this doesn’t weight by concentration - a cell is either fully counted or not.

Parameters:

• sic (array-like) – Sea ice concentration [0-1], shape (…, n2d). Can have any leading dimensions (time, ensemble, etc.).

• node_area (np.ndarray) – Area of each node in m², shape (n2d,) or (nlev, n2d). If 2D, uses the surface level (index 0).

• lat (np.ndarray) – Latitude of nodes in degrees, shape (n2d,).

• hemisphere ({"N", "S", "both"}) – “N” for Northern Hemisphere, “S” for Southern Hemisphere, “both” for global. Default is “N”.

• threshold (float) – Concentration threshold (default 0.15 = 15%). Cells with concentration > threshold are counted.

• mask (np.ndarray, optional) – Additional boolean mask, shape (n2d,). True = include.

Returns:

Total ice extent in m², shape (…). Same type as input (numpy or xarray).

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> sic = xr.open_dataset("a_ice.nc")["a_ice"]
>>> # Standard ice extent with 15% threshold
>>> extent = fesomp.diag.ice_extent(
...     sic,
...     mesh.geometry.node_area[0],
...     mesh.lat,
...     hemisphere="N",
...     threshold=0.15
... )

fesomp.diag.hovmoller(data, node_area, depth=None, mask=None)

Compute area-weighted mean profile over time (Hovmöller diagram data).

Computes the horizontal mean at each depth level for each time step, weighted by node areas.

Parameters:

• data (array-like) – 3D data, shape (time, nlev, n2d) for 3D fields. For xarray, dimensions should be named appropriately.

• node_area (np.ndarray) – Area at each node and level, shape (nlev, n2d) or (n2d,). If 1D, same area is used for all levels.

• depth (np.ndarray, optional) – Depth values for the vertical coordinate. If None and data is xarray, tries to infer from coordinates.

• mask (np.ndarray, optional) – Boolean mask, shape (n2d,). True = include, False = exclude.

Returns:

Area-weighted mean, shape (time, nlev). If xarray input with time and depth coords, returns xarray with proper coordinates.

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> temp = xr.open_dataset("temp.nc")["temp"]  # (time, lev, n2d)
>>> # Compute global mean temperature profile over time
>>> hovm = fesomp.diag.hovmoller(
...     temp,
...     mesh.geometry.node_area,
...     depth=mesh.depth_layers
... )

fesomp.diag.volume_mean(data, node_area, depth_levels, layer_thickness=None, depth_range=None, mask=None)

Compute volume-weighted mean over a depth range.

Parameters:

• data (array-like) – 3D data on layers, shape (…, nlev-1, n2d). Can have any leading dimensions (time, ensemble, etc.).

• node_area (np.ndarray) – Area at each node, shape (n2d,) or (nlev, n2d). If 2D and it has one extra level compared to data layers, the extra level is dropped with a warning (levels vs layers).

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• layer_thickness (np.ndarray, optional) – Thickness of each layer in meters, shape (nlev-1,). If None, computed from depth_levels.

• depth_range (tuple[float, float], optional) – (min_depth, max_depth) in meters for averaging. If None, uses full depth range.

• mask (np.ndarray, optional) – Boolean mask for nodes, shape (n2d,). True = include.

Returns:

Volume-weighted mean, shape (…).

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> temp = xr.open_dataset("temp.nc")["temp"]  # (time, lev, n2d)
>>> # Mean temperature in top 1000m
>>> mean_temp = fesomp.diag.volume_mean(
...     temp,
...     mesh.geometry.node_area[0],
...     mesh.depth_levels,
...     depth_range=(0, 1000)
... )

fesomp.diag.heat_content(temperature, node_area, depth_levels, layer_thickness=None, depth_range=None, mask=None, rho=1025.0, cp=3985.0, reference_temp=0.0)

Compute ocean heat content over a depth range.

Heat content is computed as:

OHC = sum(rho * cp * (T - T_ref) * volume)

Parameters:

• temperature (array-like) – Temperature in °C, shape (…, nlev-1, n2d).

• node_area (np.ndarray) – Area at each node in m², shape (n2d,) or (nlev, n2d). If 2D and it has one extra level compared to layer_thickness, the extra level is dropped with a warning (levels vs layers). If 2D and it has one extra level compared to data layers, the extra level is dropped with a warning (levels vs layers).

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• layer_thickness (np.ndarray, optional) – Thickness of each layer in meters, shape (nlev-1,).

• depth_range (tuple[float, float], optional) – (min_depth, max_depth) in meters. If None, uses full depth range.

• mask (np.ndarray, optional) – Boolean mask for nodes, shape (n2d,). True = include.

• rho (float) – Reference density in kg/m³. Default: 1025.0

• cp (float) – Specific heat capacity in J/(kg·K). Default: 3985.0

• reference_temp (float) – Reference temperature in °C. Default: 0.0

Returns:

Ocean heat content in Joules, shape (…).

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> temp = xr.open_dataset("temp.nc")["temp"]
>>> # Heat content in top 700m (standard metric)
>>> ohc_700 = fesomp.diag.heat_content(
...     temp,
...     mesh.geometry.node_area[0],
...     mesh.depth_levels,
...     depth_range=(0, 700)
... )
>>> # Heat content in top 2000m
>>> ohc_2000 = fesomp.diag.heat_content(
...     temp,
...     mesh.geometry.node_area[0],
...     mesh.depth_levels,
...     depth_range=(0, 2000)
... )

fesomp.diag.total_volume(node_area, depth_levels, layer_thickness=None, depth_range=None, mask=None)

Compute total ocean volume over a depth range.

Parameters:

• node_area (np.ndarray) – Area at each node in m², shape (n2d,) or (nlev, n2d).

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• layer_thickness (np.ndarray, optional) – Thickness of each layer in meters, shape (nlev-1,).

• depth_range (tuple[float, float], optional) – (min_depth, max_depth) in meters.

• mask (np.ndarray, optional) – Boolean mask for nodes, shape (n2d,). True = include.

Returns:

Total volume in m³.

Return type:

float

fesomp.diag.mixed_layer_depth(data, depth_levels, threshold, criterion='threshold', reference_depth=10.0)

Compute mixed layer depth using threshold criterion.

The mixed layer depth is defined as the depth where the property (density or temperature) differs from the surface (or reference depth) value by more than the threshold.

Parameters:

• data (array-like) – Vertical profile data, shape (…, nlev, n2d). For density criterion: potential density in kg/m³. For temperature criterion: temperature in °C.

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• threshold (float) – Threshold value for criterion: - density: typical values 0.03 or 0.125 kg/m³ - temperature: typical values 0.2 or 0.5 °C

• criterion ({"density", "temperature", "threshold"}) – Type of criterion: - “density”: MLD where density > surface + threshold - “temperature”: MLD where temperature < surface - threshold - “threshold”: MLD where |data - surface| > threshold

• reference_depth (float) – Reference depth in meters for surface value. Default: 10.0 Uses the level closest to this depth.

Returns:

Mixed layer depth in meters, shape (…, n2d).

Return type:

array-like

Notes

Common threshold values: - de Boyer Montégut (2004): 0.03 kg/m³ density, 0.2°C temperature - Kara et al. (2000): 0.8°C temperature - Holte & Talley (2009): 0.03 kg/m³ density

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> # Using temperature criterion with 0.2°C threshold
>>> temp = xr.open_dataset("temp.nc")["temp"]  # (time, lev, n2d)
>>> mld = fesomp.diag.mixed_layer_depth(
...     temp,
...     mesh.depth_levels,
...     threshold=0.2,
...     criterion="temperature"
... )
>>> # Using density criterion with 0.03 kg/m³ threshold
>>> rho = xr.open_dataset("rho.nc")["rho"]  # (time, lev, n2d)
>>> mld = fesomp.diag.mixed_layer_depth(
...     rho,
...     mesh.depth_levels,
...     threshold=0.03,
...     criterion="density"
... )

fesomp.diag.mixed_layer_depth_interpolated(data, depth_levels, threshold, criterion='threshold', reference_depth=10.0)

Compute mixed layer depth with linear interpolation.

Similar to mixed_layer_depth but interpolates between levels to find the exact depth where threshold is crossed.

Parameters:

• data (array-like) – Vertical profile data, shape (…, nlev, n2d).

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• threshold (float) – Threshold value for criterion.

• criterion ({"density", "temperature", "threshold"}) – Type of criterion.

• reference_depth (float) – Reference depth in meters for surface value.

Returns:

Mixed layer depth in meters, shape (…, n2d). Interpolated to exact crossing depth.

Return type:

array-like

Examples

>>> mld = fesomp.diag.mixed_layer_depth_interpolated(
...     temp,
...     mesh.depth_levels,
...     threshold=0.2,
...     criterion="temperature"
... )

fesomp.diag.moc(w, node_area, lat, depth_levels, lat_range=(-90, 90), nlats=181, mask=None, cumsum_direction='bottom_up')

Compute Meridional Overturning Circulation from vertical velocity.

In FESOM2, vertical velocity (w) is defined on nodes at level interfaces. The MOC streamfunction is computed by: 1. Binning vertical velocity by latitude 2. Weighting by node area 3. Computing cumulative sum in depth

Parameters:

• w (array-like) – Vertical velocity in m/s, shape (…, nz, n2d). Positive = upward. Data is on level interfaces (nz levels).

• node_area (np.ndarray) – Node areas in m², shape (n2d,) or (nz, n2d). If 2D, uses level-dependent areas; if 1D, uses surface area for all levels.

• lat (np.ndarray) – Latitude of nodes in degrees, shape (n2d,).

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nz,).

• lat_range (tuple[float, float]) – (min_lat, max_lat) for output. Default: (-90, 90).

• nlats (int) – Number of latitude bins. Default: 181 (1° resolution).

• mask (np.ndarray, optional) – Boolean mask for nodes, shape (n2d,). True = include. Use this to select specific basins (Atlantic, Indo-Pacific, etc.).

• cumsum_direction ({"bottom_up", "top_down"}) – Direction for cumulative sum: - “bottom_up”: integrate from bottom to surface (standard) - “top_down”: integrate from surface to bottom

Return type:

tuple[ndarray | DataArray, ndarray]

Returns:

• moc (array-like) – MOC streamfunction in Sv (10⁶ m³/s), shape (…, nz, nlats).

• lat_bins (np.ndarray) – Latitude bin centers, shape (nlats,).

Notes

The MOC streamfunction Ψ(y, z) represents the zonally-integrated volume transport. At each latitude and depth:

Ψ = ∫∫ w dA

where the integral is over all nodes at that latitude, summed cumulatively from bottom (or top).

To compute Atlantic MOC, use:

mask = fesomp.diag.get_basin_mask(mesh.lon, mesh.lat, “Atlantic_MOC”) amoc, lats = fesomp.diag.moc(w, node_area, mesh.lat, depth_levels, mask=mask)

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> w = xr.open_dataset("w.nc")["w"]  # (time, nz, n2d)
>>> # Global MOC
>>> moc_global, lats = fesomp.diag.moc(
...     w,
...     mesh.geometry.node_area[0],  # surface node areas
...     mesh.lat,
...     mesh.depth_levels
... )
>>> # Atlantic MOC
>>> atl_mask = fesomp.diag.get_basin_mask(
...     mesh.lon, mesh.lat, "Atlantic_MOC"
... )
>>> amoc, lats = fesomp.diag.moc(
...     w, mesh.geometry.node_area[0], mesh.lat,
...     mesh.depth_levels, mask=atl_mask
... )

fesomp.diag.amoc_index(moc_data, lat_bins, depth_levels, lat=26.5, depth_range=(500, 2000))

Compute AMOC index as maximum streamfunction at specified latitude.

The standard AMOC index is the maximum overturning at 26.5°N (RAPID array latitude) in the depth range 500-2000m.

Parameters:

• moc_data (array-like) – MOC streamfunction in Sv, shape (…, nz, nlats).

• lat_bins (np.ndarray) – Latitude bin centers, shape (nlats,).

• depth_levels (np.ndarray) – Depth at levels, shape (nz,).

• lat (float) – Latitude for index calculation. Default: 26.5°N.

• depth_range (tuple[float, float]) – Depth range for maximum search. Default: (500, 2000).

Returns:

AMOC index in Sv, shape (…).

Return type:

array-like

Examples

>>> # Compute Atlantic MOC
>>> atl_mask = fesomp.diag.get_basin_mask(mesh.lon, mesh.lat, "Atlantic_MOC")
>>> amoc, lats = fesomp.diag.moc(w, node_area, mesh.lat, depth_levels, mask=atl_mask)
>>> # Get AMOC index at 26.5N
>>> amoc_26 = fesomp.diag.amoc_index(amoc, lats, depth_levels)

fesomp.diag.get_basin_mask(lon, lat, basin)

Get boolean mask for a named ocean basin.

Uses the MOCBasins.geojson file included with the package.

Parameters:

• lon (np.ndarray) – Longitude of points in degrees, shape (n,).

• lat (np.ndarray) – Latitude of points in degrees, shape (n,).

• basin (str) – Name of the basin. Available basins: - “Atlantic_MOC”: Atlantic Ocean (including Med, Baltic, etc.) - “IndoPacific_MOC”: Indo-Pacific Ocean - “Pacific_MOC”: Pacific Ocean - “Indian_MOC”: Indian Ocean Use list_moc_basins() to see all available basins.

Returns:

Boolean mask, shape (n,). True = inside basin.

Return type:

np.ndarray

Raises:

• ImportError – If shapely is not installed.

• ValueError – If basin name is not found.

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> # Get Atlantic mask for nodes
>>> atl_mask = fesomp.diag.get_basin_mask(
...     mesh.lon, mesh.lat, "Atlantic_MOC"
... )

fesomp.diag.list_moc_basins()

List available MOC basin names.

Returns:

List of available basin names for use with get_basin_mask().

Return type:

list[str]

fesomp.diag.hemisphere_mask(lat, hemisphere)

Create boolean mask for hemisphere selection.

Parameters:

• lat (np.ndarray) – Latitude of points in degrees, shape (n,).

• hemisphere ({"N", "S", "both"}) – “N” for Northern Hemisphere (lat >= 0), “S” for Southern Hemisphere (lat < 0), “both” for global (all points).

Returns:

Boolean mask, shape (n,). True = include.

Return type:

np.ndarray

Raises:

ValueError – If hemisphere is not one of “N”, “S”, “both”.

fesomp.diag.get_region_mask(lon, lat, region, geojson_name='MOCBasins')

Create boolean mask for points within a named region.

Requires shapely to be installed.

Parameters:

• lon (np.ndarray) – Longitude of points in degrees, shape (n,).

• lat (np.ndarray) – Latitude of points in degrees, shape (n,).

• region (str) – Name of the region (e.g., “Atlantic_Basin”, “Pacific_Basin”). Use list_available_regions() to see available names.

• geojson_name (str) – Name of GeoJSON file: “MOCBasins”, “NinoRegions”, or “oceanBasins”.

Returns:

Boolean mask, shape (n,). True = inside region.

Return type:

np.ndarray

Raises:

• ImportError – If shapely is not installed.

• ValueError – If region name is not found.

fesomp.diag.list_available_regions(geojson_name='MOCBasins')

List available region names in a GeoJSON file.

Parameters:

geojson_name (str) – Name of GeoJSON file: “MOCBasins”, “NinoRegions”, or “oceanBasins”.

Returns:

List of region names available in the file.

Return type:

list[str]

fesomp.diag.compute_layer_thickness(depth_levels)

Compute layer thickness from level interface depths.

Parameters:

depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,). Should be positive, increasing downward.

Returns:

Layer thickness in meters, shape (nlev-1,).

Return type:

np.ndarray

fesomp.diag.select_depth_indices(depth_levels, depth_range=None)

Find level indices corresponding to a depth range.

Parameters:

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• depth_range (tuple[float, float], optional) – (min_depth, max_depth) in meters. If None, returns (0, nlev-1) for full depth.

Returns:

(start_index, end_index) for slicing. Use as depth_levels[start:end+1] or layers[start:end].

Return type:

tuple[int, int]

fesomp.diag.get_surface_area(node_area)

Extract surface area from node_area array.

Parameters:

node_area (np.ndarray) – Area at each node, shape (n2d,) or (nlev, n2d).

Returns:

Surface area, shape (n2d,).

Return type:

np.ndarray

Sea Ice

Sea ice diagnostics: area, volume, and extent.

fesomp.diag.ice.ice_area(sic, node_area, lat, hemisphere='N', mask=None)

Compute sea ice area (ice concentration × cell area).

Sea ice area is the sum of (ice concentration × cell area) over all cells. This represents the actual area covered by ice.

Parameters:

• sic (array-like) – Sea ice concentration [0-1], shape (…, n2d). Can have any leading dimensions (time, ensemble, etc.).

• node_area (np.ndarray) – Area of each node in m², shape (n2d,) or (nlev, n2d). If 2D, uses the surface level (index 0).

• lat (np.ndarray) – Latitude of nodes in degrees, shape (n2d,).

• hemisphere ({"N", "S", "both"}) – “N” for Northern Hemisphere, “S” for Southern Hemisphere, “both” for global. Default is “N”.

• mask (np.ndarray, optional) – Additional boolean mask, shape (n2d,). True = include.

Returns:

Total ice area in m², shape (…). Same type as input (numpy or xarray).

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> sic = xr.open_dataset("a_ice.nc")["a_ice"]  # shape: (time, n2d)
>>> # Compute Northern Hemisphere ice area time series
>>> nh_area = fesomp.diag.ice_area(
...     sic,
...     mesh.geometry.node_area[0],
...     mesh.lat,
...     hemisphere="N"
... )

fesomp.diag.ice.ice_volume(siv, node_area, lat, hemisphere='N', mask=None)

Compute sea ice volume (ice thickness × cell area).

Sea ice volume is the sum of (effective ice thickness × cell area). The input can be either: - Effective ice thickness (m): already area-weighted ice thickness - Ice volume per unit area (m): same as effective thickness

Parameters:

• siv (array-like) – Sea ice effective thickness or volume per unit area in meters, shape (…, n2d). Can have any leading dimensions.

• node_area (np.ndarray) – Area of each node in m², shape (n2d,) or (nlev, n2d). If 2D, uses the surface level (index 0).

• lat (np.ndarray) – Latitude of nodes in degrees, shape (n2d,).

• hemisphere ({"N", "S", "both"}) – “N” for Northern Hemisphere, “S” for Southern Hemisphere, “both” for global. Default is “N”.

• mask (np.ndarray, optional) – Additional boolean mask, shape (n2d,). True = include.

Returns:

Total ice volume in m³, shape (…). Same type as input (numpy or xarray).

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> siv = xr.open_dataset("m_ice.nc")["m_ice"]  # shape: (time, n2d)
>>> # Compute Arctic ice volume time series
>>> arctic_vol = fesomp.diag.ice_volume(
...     siv,
...     mesh.geometry.node_area[0],
...     mesh.lat,
...     hemisphere="N"
... )

fesomp.diag.ice.ice_extent(sic, node_area, lat, hemisphere='N', threshold=0.15, mask=None)

Compute sea ice extent (area where concentration exceeds threshold).

Sea ice extent is the total area of cells where ice concentration exceeds the threshold (typically 15%). Unlike ice area, this doesn’t weight by concentration - a cell is either fully counted or not.

Parameters:

• sic (array-like) – Sea ice concentration [0-1], shape (…, n2d). Can have any leading dimensions (time, ensemble, etc.).

• node_area (np.ndarray) – Area of each node in m², shape (n2d,) or (nlev, n2d). If 2D, uses the surface level (index 0).

• lat (np.ndarray) – Latitude of nodes in degrees, shape (n2d,).

• hemisphere ({"N", "S", "both"}) – “N” for Northern Hemisphere, “S” for Southern Hemisphere, “both” for global. Default is “N”.

• threshold (float) – Concentration threshold (default 0.15 = 15%). Cells with concentration > threshold are counted.

• mask (np.ndarray, optional) – Additional boolean mask, shape (n2d,). True = include.

Returns:

Total ice extent in m², shape (…). Same type as input (numpy or xarray).

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> sic = xr.open_dataset("a_ice.nc")["a_ice"]
>>> # Standard ice extent with 15% threshold
>>> extent = fesomp.diag.ice_extent(
...     sic,
...     mesh.geometry.node_area[0],
...     mesh.lat,
...     hemisphere="N",
...     threshold=0.15
... )

Vertical Diagnostics

Vertical diagnostics: Hovmöller diagrams, volume means, heat content.

fesomp.diag.vertical.hovmoller(data, node_area, depth=None, mask=None)

Compute area-weighted mean profile over time (Hovmöller diagram data).

Computes the horizontal mean at each depth level for each time step, weighted by node areas.

Parameters:

• data (array-like) – 3D data, shape (time, nlev, n2d) for 3D fields. For xarray, dimensions should be named appropriately.

• node_area (np.ndarray) – Area at each node and level, shape (nlev, n2d) or (n2d,). If 1D, same area is used for all levels.

• depth (np.ndarray, optional) – Depth values for the vertical coordinate. If None and data is xarray, tries to infer from coordinates.

• mask (np.ndarray, optional) – Boolean mask, shape (n2d,). True = include, False = exclude.

Returns:

Area-weighted mean, shape (time, nlev). If xarray input with time and depth coords, returns xarray with proper coordinates.

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> temp = xr.open_dataset("temp.nc")["temp"]  # (time, lev, n2d)
>>> # Compute global mean temperature profile over time
>>> hovm = fesomp.diag.hovmoller(
...     temp,
...     mesh.geometry.node_area,
...     depth=mesh.depth_layers
... )

fesomp.diag.vertical.volume_mean(data, node_area, depth_levels, layer_thickness=None, depth_range=None, mask=None)

Compute volume-weighted mean over a depth range.

Parameters:

• data (array-like) – 3D data on layers, shape (…, nlev-1, n2d). Can have any leading dimensions (time, ensemble, etc.).

• node_area (np.ndarray) – Area at each node, shape (n2d,) or (nlev, n2d). If 2D and it has one extra level compared to data layers, the extra level is dropped with a warning (levels vs layers).

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• layer_thickness (np.ndarray, optional) – Thickness of each layer in meters, shape (nlev-1,). If None, computed from depth_levels.

• depth_range (tuple[float, float], optional) – (min_depth, max_depth) in meters for averaging. If None, uses full depth range.

• mask (np.ndarray, optional) – Boolean mask for nodes, shape (n2d,). True = include.

Returns:

Volume-weighted mean, shape (…).

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> temp = xr.open_dataset("temp.nc")["temp"]  # (time, lev, n2d)
>>> # Mean temperature in top 1000m
>>> mean_temp = fesomp.diag.volume_mean(
...     temp,
...     mesh.geometry.node_area[0],
...     mesh.depth_levels,
...     depth_range=(0, 1000)
... )

fesomp.diag.vertical.heat_content(temperature, node_area, depth_levels, layer_thickness=None, depth_range=None, mask=None, rho=1025.0, cp=3985.0, reference_temp=0.0)

Compute ocean heat content over a depth range.

Heat content is computed as:

OHC = sum(rho * cp * (T - T_ref) * volume)

Parameters:

• temperature (array-like) – Temperature in °C, shape (…, nlev-1, n2d).

• node_area (np.ndarray) – Area at each node in m², shape (n2d,) or (nlev, n2d). If 2D and it has one extra level compared to layer_thickness, the extra level is dropped with a warning (levels vs layers). If 2D and it has one extra level compared to data layers, the extra level is dropped with a warning (levels vs layers).

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• layer_thickness (np.ndarray, optional) – Thickness of each layer in meters, shape (nlev-1,).

• depth_range (tuple[float, float], optional) – (min_depth, max_depth) in meters. If None, uses full depth range.

• mask (np.ndarray, optional) – Boolean mask for nodes, shape (n2d,). True = include.

• rho (float) – Reference density in kg/m³. Default: 1025.0

• cp (float) – Specific heat capacity in J/(kg·K). Default: 3985.0

• reference_temp (float) – Reference temperature in °C. Default: 0.0

Returns:

Ocean heat content in Joules, shape (…).

Return type:

array-like

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> temp = xr.open_dataset("temp.nc")["temp"]
>>> # Heat content in top 700m (standard metric)
>>> ohc_700 = fesomp.diag.heat_content(
...     temp,
...     mesh.geometry.node_area[0],
...     mesh.depth_levels,
...     depth_range=(0, 700)
... )
>>> # Heat content in top 2000m
>>> ohc_2000 = fesomp.diag.heat_content(
...     temp,
...     mesh.geometry.node_area[0],
...     mesh.depth_levels,
...     depth_range=(0, 2000)
... )

fesomp.diag.vertical.total_volume(node_area, depth_levels, layer_thickness=None, depth_range=None, mask=None)

Compute total ocean volume over a depth range.

Parameters:

• node_area (np.ndarray) – Area at each node in m², shape (n2d,) or (nlev, n2d).

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• layer_thickness (np.ndarray, optional) – Thickness of each layer in meters, shape (nlev-1,).

• depth_range (tuple[float, float], optional) – (min_depth, max_depth) in meters.

• mask (np.ndarray, optional) – Boolean mask for nodes, shape (n2d,). True = include.

Returns:

Total volume in m³.

Return type:

float

Mixed Layer

Mixed layer depth diagnostics.

fesomp.diag.mixed_layer.mixed_layer_depth(data, depth_levels, threshold, criterion='threshold', reference_depth=10.0)

Compute mixed layer depth using threshold criterion.

The mixed layer depth is defined as the depth where the property (density or temperature) differs from the surface (or reference depth) value by more than the threshold.

Parameters:

• data (array-like) – Vertical profile data, shape (…, nlev, n2d). For density criterion: potential density in kg/m³. For temperature criterion: temperature in °C.

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• threshold (float) – Threshold value for criterion: - density: typical values 0.03 or 0.125 kg/m³ - temperature: typical values 0.2 or 0.5 °C

• criterion ({"density", "temperature", "threshold"}) – Type of criterion: - “density”: MLD where density > surface + threshold - “temperature”: MLD where temperature < surface - threshold - “threshold”: MLD where |data - surface| > threshold

• reference_depth (float) – Reference depth in meters for surface value. Default: 10.0 Uses the level closest to this depth.

Returns:

Mixed layer depth in meters, shape (…, n2d).

Return type:

array-like

Notes

Common threshold values: - de Boyer Montégut (2004): 0.03 kg/m³ density, 0.2°C temperature - Kara et al. (2000): 0.8°C temperature - Holte & Talley (2009): 0.03 kg/m³ density

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> # Using temperature criterion with 0.2°C threshold
>>> temp = xr.open_dataset("temp.nc")["temp"]  # (time, lev, n2d)
>>> mld = fesomp.diag.mixed_layer_depth(
...     temp,
...     mesh.depth_levels,
...     threshold=0.2,
...     criterion="temperature"
... )
>>> # Using density criterion with 0.03 kg/m³ threshold
>>> rho = xr.open_dataset("rho.nc")["rho"]  # (time, lev, n2d)
>>> mld = fesomp.diag.mixed_layer_depth(
...     rho,
...     mesh.depth_levels,
...     threshold=0.03,
...     criterion="density"
... )

fesomp.diag.mixed_layer.mixed_layer_depth_interpolated(data, depth_levels, threshold, criterion='threshold', reference_depth=10.0)

Compute mixed layer depth with linear interpolation.

Similar to mixed_layer_depth but interpolates between levels to find the exact depth where threshold is crossed.

Parameters:

• data (array-like) – Vertical profile data, shape (…, nlev, n2d).

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• threshold (float) – Threshold value for criterion.

• criterion ({"density", "temperature", "threshold"}) – Type of criterion.

• reference_depth (float) – Reference depth in meters for surface value.

Returns:

Mixed layer depth in meters, shape (…, n2d). Interpolated to exact crossing depth.

Return type:

array-like

Examples

>>> mld = fesomp.diag.mixed_layer_depth_interpolated(
...     temp,
...     mesh.depth_levels,
...     threshold=0.2,
...     criterion="temperature"
... )

Meridional Overturning

Meridional Overturning Circulation (MOC) diagnostics.

fesomp.diag.moc.moc(w, node_area, lat, depth_levels, lat_range=(-90, 90), nlats=181, mask=None, cumsum_direction='bottom_up')

Compute Meridional Overturning Circulation from vertical velocity.

In FESOM2, vertical velocity (w) is defined on nodes at level interfaces. The MOC streamfunction is computed by: 1. Binning vertical velocity by latitude 2. Weighting by node area 3. Computing cumulative sum in depth

Parameters:

• w (array-like) – Vertical velocity in m/s, shape (…, nz, n2d). Positive = upward. Data is on level interfaces (nz levels).

• node_area (np.ndarray) – Node areas in m², shape (n2d,) or (nz, n2d). If 2D, uses level-dependent areas; if 1D, uses surface area for all levels.

• lat (np.ndarray) – Latitude of nodes in degrees, shape (n2d,).

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nz,).

• lat_range (tuple[float, float]) – (min_lat, max_lat) for output. Default: (-90, 90).

• nlats (int) – Number of latitude bins. Default: 181 (1° resolution).

• mask (np.ndarray, optional) – Boolean mask for nodes, shape (n2d,). True = include. Use this to select specific basins (Atlantic, Indo-Pacific, etc.).

• cumsum_direction ({"bottom_up", "top_down"}) – Direction for cumulative sum: - “bottom_up”: integrate from bottom to surface (standard) - “top_down”: integrate from surface to bottom

Return type:

tuple[ndarray | DataArray, ndarray]

Returns:

• moc (array-like) – MOC streamfunction in Sv (10⁶ m³/s), shape (…, nz, nlats).

• lat_bins (np.ndarray) – Latitude bin centers, shape (nlats,).

Notes

The MOC streamfunction Ψ(y, z) represents the zonally-integrated volume transport. At each latitude and depth:

Ψ = ∫∫ w dA

where the integral is over all nodes at that latitude, summed cumulatively from bottom (or top).

To compute Atlantic MOC, use:

mask = fesomp.diag.get_basin_mask(mesh.lon, mesh.lat, “Atlantic_MOC”) amoc, lats = fesomp.diag.moc(w, node_area, mesh.lat, depth_levels, mask=mask)

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> w = xr.open_dataset("w.nc")["w"]  # (time, nz, n2d)
>>> # Global MOC
>>> moc_global, lats = fesomp.diag.moc(
...     w,
...     mesh.geometry.node_area[0],  # surface node areas
...     mesh.lat,
...     mesh.depth_levels
... )
>>> # Atlantic MOC
>>> atl_mask = fesomp.diag.get_basin_mask(
...     mesh.lon, mesh.lat, "Atlantic_MOC"
... )
>>> amoc, lats = fesomp.diag.moc(
...     w, mesh.geometry.node_area[0], mesh.lat,
...     mesh.depth_levels, mask=atl_mask
... )

fesomp.diag.moc.get_basin_mask(lon, lat, basin)

Get boolean mask for a named ocean basin.

Uses the MOCBasins.geojson file included with the package.

Parameters:

• lon (np.ndarray) – Longitude of points in degrees, shape (n,).

• lat (np.ndarray) – Latitude of points in degrees, shape (n,).

• basin (str) – Name of the basin. Available basins: - “Atlantic_MOC”: Atlantic Ocean (including Med, Baltic, etc.) - “IndoPacific_MOC”: Indo-Pacific Ocean - “Pacific_MOC”: Pacific Ocean - “Indian_MOC”: Indian Ocean Use list_moc_basins() to see all available basins.

Returns:

Boolean mask, shape (n,). True = inside basin.

Return type:

np.ndarray

Raises:

• ImportError – If shapely is not installed.

• ValueError – If basin name is not found.

Examples

>>> import fesomp
>>> mesh = fesomp.load_mesh("path/to/mesh.nc")
>>> # Get Atlantic mask for nodes
>>> atl_mask = fesomp.diag.get_basin_mask(
...     mesh.lon, mesh.lat, "Atlantic_MOC"
... )

fesomp.diag.moc.list_moc_basins()

List available MOC basin names.

Returns:

List of available basin names for use with get_basin_mask().

Return type:

list[str]

fesomp.diag.moc.amoc_index(moc_data, lat_bins, depth_levels, lat=26.5, depth_range=(500, 2000))

Compute AMOC index as maximum streamfunction at specified latitude.

The standard AMOC index is the maximum overturning at 26.5°N (RAPID array latitude) in the depth range 500-2000m.

Parameters:

• moc_data (array-like) – MOC streamfunction in Sv, shape (…, nz, nlats).

• lat_bins (np.ndarray) – Latitude bin centers, shape (nlats,).

• depth_levels (np.ndarray) – Depth at levels, shape (nz,).

• lat (float) – Latitude for index calculation. Default: 26.5°N.

• depth_range (tuple[float, float]) – Depth range for maximum search. Default: (500, 2000).

Returns:

AMOC index in Sv, shape (…).

Return type:

array-like

Examples

>>> # Compute Atlantic MOC
>>> atl_mask = fesomp.diag.get_basin_mask(mesh.lon, mesh.lat, "Atlantic_MOC")
>>> amoc, lats = fesomp.diag.moc(w, node_area, mesh.lat, depth_levels, mask=atl_mask)
>>> # Get AMOC index at 26.5N
>>> amoc_26 = fesomp.diag.amoc_index(amoc, lats, depth_levels)

Utilities

Utility functions for diagnostics module.

fesomp.diag.utils.hemisphere_mask(lat, hemisphere)

Create boolean mask for hemisphere selection.

Parameters:

• lat (np.ndarray) – Latitude of points in degrees, shape (n,).

• hemisphere ({"N", "S", "both"}) – “N” for Northern Hemisphere (lat >= 0), “S” for Southern Hemisphere (lat < 0), “both” for global (all points).

Returns:

Boolean mask, shape (n,). True = include.

Return type:

np.ndarray

Raises:

ValueError – If hemisphere is not one of “N”, “S”, “both”.

fesomp.diag.utils.compute_layer_thickness(depth_levels)

Compute layer thickness from level interface depths.

Parameters:

depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,). Should be positive, increasing downward.

Returns:

Layer thickness in meters, shape (nlev-1,).

Return type:

np.ndarray

fesomp.diag.utils.select_depth_indices(depth_levels, depth_range=None)

Find level indices corresponding to a depth range.

Parameters:

• depth_levels (np.ndarray) – Depth at level interfaces in meters, shape (nlev,).

• depth_range (tuple[float, float], optional) – (min_depth, max_depth) in meters. If None, returns (0, nlev-1) for full depth.

Returns:

(start_index, end_index) for slicing. Use as depth_levels[start:end+1] or layers[start:end].

Return type:

tuple[int, int]

fesomp.diag.utils.get_surface_area(node_area)

Extract surface area from node_area array.

Parameters:

node_area (np.ndarray) – Area at each node, shape (n2d,) or (nlev, n2d).

Returns:

Surface area, shape (n2d,).

Return type:

np.ndarray

fesomp.diag.utils.weighted_sum(data, weights, mask=None, axis=-1)

Compute weighted sum along specified axis.

Handles both numpy arrays and xarray DataArrays (including dask).

Parameters:

• data (array-like) – Data array, arbitrary shape.

• weights (np.ndarray) – Weights, must be broadcastable to data shape along axis.

• mask (np.ndarray, optional) – Boolean mask for points to include. Same shape as weights.

• axis (int) – Axis along which to sum.

Returns:

Weighted sum, with axis removed.

Return type:

array-like

fesomp.diag.utils.weighted_mean(data, weights, mask=None, axis=-1)

Compute weighted mean along specified axis.

Handles both numpy arrays and xarray DataArrays (including dask).

Parameters:

• data (array-like) – Data array, arbitrary shape.

• weights (np.ndarray) – Weights, must be broadcastable to data shape along axis.

• mask (np.ndarray, optional) – Boolean mask for points to include. Same shape as weights.

• axis (int) – Axis along which to compute mean.

Returns:

Weighted mean, with axis removed.

Return type:

array-like

fesomp.diag.utils.load_geojson(name)

Load a GeoJSON file from the package data directory.

Parameters:

name (str) – Name of the GeoJSON file (without .geojson extension). Available: “MOCBasins”, “NinoRegions”, “oceanBasins”

Returns:

Parsed GeoJSON data.

Return type:

dict

Raises:

FileNotFoundError – If the GeoJSON file does not exist.

fesomp.diag.utils.list_available_regions(geojson_name='MOCBasins')

List available region names in a GeoJSON file.

Parameters:

geojson_name (str) – Name of GeoJSON file: “MOCBasins”, “NinoRegions”, or “oceanBasins”.

Returns:

List of region names available in the file.

Return type:

list[str]

fesomp.diag.utils.get_region_mask(lon, lat, region, geojson_name='MOCBasins')

Create boolean mask for points within a named region.

Requires shapely to be installed.

Parameters:

• lon (np.ndarray) – Longitude of points in degrees, shape (n,).

• lat (np.ndarray) – Latitude of points in degrees, shape (n,).

• region (str) – Name of the region (e.g., “Atlantic_Basin”, “Pacific_Basin”). Use list_available_regions() to see available names.

• geojson_name (str) – Name of GeoJSON file: “MOCBasins”, “NinoRegions”, or “oceanBasins”.

Returns:

Boolean mask, shape (n,). True = inside region.

Return type:

np.ndarray

Raises:

• ImportError – If shapely is not installed.

• ValueError – If region name is not found.