API reference

API reference#

`regional_mom6`#

class regional_mom6.regional_mom6.experiment(*, longitude_extent, latitude_extent, date_range, resolution, number_vertical_layers, layer_thickness_ratio, depth, mom_run_dir, mom_input_dir, toolpath_dir, grid_type='even_spacing', repeat_year_forcing=False, read_existing_grids=False)#

The main class for setting up a regional experiment.

Everything about the regional experiment.

Methods in this class generate the various input files needed for a MOM6 experiment forced with open boundary conditions (OBCs). The code is agnostic to the user’s choice of boundary forcing, bathymetry, and surface forcing; users need to prescribe what variables are all called via mapping dictionaries from MOM6 variable/coordinate name to the name in the input dataset.

The class can be used to generate the grids for a new experiment, or to read in an existing one (when read_existing_grids=True; see argument description below).

Parameters:

longitude_extent (Tuple[float]) – Extent of the region in longitude (in degrees). For example: (40.5, 50.0).
latitude_extent (Tuple[float]) – Extent of the region in latitude (in degrees). For example: (-20.0, 30.0).
date_range (Tuple[str]) – Start and end dates of the boundary forcing window. For example: ("2003-01-01", "2003-01-31").
resolution (float) – Lateral resolution of the domain (in degrees).
number_vertical_layers (int) – Number of vertical layers.
layer_thickness_ratio (float) – Ratio of largest to smallest layer thickness; used as input in hyperbolictan_thickness_profile().
depth (float) – Depth of the domain.
mom_run_dir (str) – Path of the MOM6 control directory.
mom_input_dir (str) – Path of the MOM6 input directory, to receive the forcing files.
toolpath_dir (str) – Path of GFDL’s FRE tools (NOAA-GFDL/FRE-NCtools) binaries.
grid_type (Optional[str]) – Type of horizontal grid to generate. Currently, only 'even_spacing' is supported.
repeat_year_forcing (Optional[bool]) – When True the experiment runs with repeat-year forcing. When False (default) then inter-annual forcing is used.
read_existing_grids (Optional[Bool]) – When True, instead of generating the grids, the grids and the ocean mask are being read from within the mom_input_dir and mom_run_dir directories. Useful for modifying or troubleshooting experiments. Default: False.

FRE_tools(layout=None)#: A wrapper for FRE Tools check_mask, make_solo_mosaic, and make_quick_mosaic. User provides processor layout tuple of processing units.

_make_hgrid()#

Set up a horizontal grid based on user’s specification of the domain. The default behaviour generates a grid evenly spaced both in longitude and in latitude.

The latitudinal resolution is scaled with the cosine of the central latitude of the domain, i.e., Δφ = cos(φ_central) * Δλ, where Δλ is the longitudinal spacing. This way, for a sufficiently small domain, the linear distances between grid points are nearly identical: Δx = R * cos(φ) * Δλ and Δy = R * Δφ = R * cos(φ_central) * Δλ (here R is Earth’s radius and φ, φ_central, Δλ, and Δφ are all expressed in radians). That is, if the domain is small enough that so that cos(φ_North_Side) is not much different from cos(φ_South_Side), then Δx and Δy are similar.

Note

The intention is for the horizontal grid (hgrid) generation to be flexible. For now, there is only one implemented horizontal grid included in the package, but you can customise it by simply overwriting the hgrid.nc file in the mom_run_dir directory after initialising an experiment. To preserve the metadata, it might be easiest to read the file in, then modify the fields before re-saving.

_make_vgrid()#: Generates a vertical grid based on the number_vertical_layers, the ratio of largest to smallest layer thickness (layer_thickness_ratio) and the total depth parameters. (All these parameters are specified at the class level.)

cpu_layout(layout)#: Wrapper for the check_mask function of GFDL’s FRE Tools. User provides processor layout tuple of processing units.

initial_condition(ic_path, varnames, arakawa_grid='A', vcoord_type='height')#

Reads the initial condition from files in ic_path, interpolates to the model grid, fixes up metadata, and saves back to the input directory.

Parameters:

ic_path (Union[str, Path]) – Path to initial condition file.
varnames (Dict[str, str]) – Mapping from MOM6 variable/coordinate names to the names in the input dataset. For example, {'xq': 'lonq', 'yh': 'lath', 'salt': 'so', ...}.
arakawa_grid (Optional[str]) – Arakawa grid staggering type of the initial condition. Either 'A' (default), 'B', or 'C'.
vcoord_type (Optional[str]) – The type of vertical coordinate used in the forcing files. Either 'height' or 'thickness'.

rectangular_boundary(path_to_bc, varnames, orientation, segment_number, arakawa_grid='A')#

Set up a boundary forcing file for a given orientation. Here the term ‘rectangular’ means boundaries along lines of constant latitude or longitude.

Parameters:

path_to_bc (str) – Path to boundary forcing file. Ideally this should be a pre cut-out netCDF file containing only the boundary region and 3 extra boundary points on either side. Users can also provide a large dataset containing their entire domain but this will be slower.
varnames (Dict[str, str]) – Mapping from MOM6 variable/coordinate names to the name in the input dataset.
orientation (str) – Orientation of boundary forcing file, i.e., 'east', 'west', 'north', or 'south'.
segment_number (int) – Number the segments according to how they’ll be specified in the MOM_input.
arakawa_grid (Optional[str]) – Arakawa grid staggering type of the boundary forcing. Either 'A' (default), 'B', or 'C'.

setup_bathymetry(*, bathymetry_path, longitude_coordinate_name='lon', latitude_coordinate_name='lat', vertical_coordinate_name='elevation', fill_channels=False, minimum_layers=3, positive_down=False, chunks='auto')#

Cut out and interpolate the chosen bathymetry and then fill inland lakes.

It’s also possible to optionally fill narrow channels (see fill_channels below), although narrow channels are less of an issue for models that are discretized on an Arakawa C grid, like MOM6.

Output is saved in the input directory of the experiment.

Parameters:

bathymetry_path (str) – Path to the netCDF file with the bathymetry.
longitude_coordinate_name (Optional[str]) – The name of the longitude coordinate in the bathymetry dataset at bathymetry_path. For example, for GEBCO bathymetry: 'lon' (default).
latitude_coordinate_name (Optional[str]) – The name of the latitude coordinate in the bathymetry dataset at bathymetry_path. For example, for GEBCO bathymetry: 'lat' (default).
vertical_coordinate_name (Optional[str]) – The name of the vertical coordinate in the bathymetry dataset at bathymetry_path. For example, for GEBCO bathymetry: 'elevation' (default).
fill_channels (Optional[bool]) – Whether or not to fill in diagonal channels. This removes more narrow inlets, but can also connect extra islands to land. Default: False.
minimum_layers (Optional[int]) – The minimum depth allowed as an integer number of layers. Anything shallower than the minimum_layers (as specified by the vertical coordinate file vcoord.nc) is deemed land. Default: 3.
positive_down (Optional[bool]) – If True, it assumes that bathymetry vertical coordinate is positive down. Default: False.
chunks (Optional Dict[str, str]) – Horizontal chunking scheme for the bathymetry, e.g., {"longitude": 100, "latitude": 100}. Use 'longitude' and 'latitude' rather than the actual coordinate names in the input file.

setup_era5(era5_path)#

Setup the ERA5 forcing files for the experiment. This assumes that all of the ERA5 data in the prescribed date range are downloaded. We need the following fields: “2t”, “10u”, “10v”, “sp”, “2d”, “msdwswrf”, “msdwlwrf”, “lsrr”, and “crr”.

Parameters:: era5_path (str) – Path to the ERA5 forcing files. Specifically, the single-level reanalysis product. For example, 'SOMEPATH/era5/single-levels/reanalysis'

setup_run_directory(surface_forcing=None, using_payu=False, overwrite=False)#

Set up the run directory for MOM6. Either copy a pre-made set of files, or modify existing files in the ‘rundir’ directory for the experiment.

Parameters:

surface_forcing (Optional[str]) – Specify the choice of surface forcing, one of: 'jra' or 'era5'. If not prescribed then constant fluxes are used.
using_payu (Optional[bool]) – Whether or not to use payu (payu-org/payu) to run the model. If True, a payu configuration file will be created. Default: False.
overwrite (Optional[bool]) – Whether or not to overwrite existing files in the run directory. If False (default), will only modify the MOM_layout file and will not re-copy across the rest of the default files.

tidy_bathymetry(fill_channels=False, minimum_layers=3, positive_down=True)#

An auxiliary function for bathymetry used to fix up the metadata and remove inland lakes after regridding the bathymetry. Having tidy_bathymetry as a separate method from setup_bathymetry() allows for the regridding to be done separately, since regridding can be really expensive for large domains.

If the bathymetry is already regridded and what is left to be done is fixing the metadata or fill in some channels, then call this function directly to read in the existing bathymetry_unfinished.nc file that should be in the input directory.

Parameters:

fill_channels (Optional[bool]) – Whether to fill in diagonal channels. This removes more narrow inlets, but can also connect extra islands to land. Default: False.
minimum_layers (Optional[int]) – The minimum depth allowed as an integer number of layers. The default value of 3 layers means that anything shallower than the 3rd layer (as specified by the vcoord) is deemed land.
positive_down (Optional[bool]) – If True (default), assume that bathymetry vertical coordinate is positive down.

regional_mom6.regional_mom6.hyperbolictan_thickness_profile(nlayers, ratio, total_depth)#

Generate a hyperbolic tangent thickness profile with nlayers vertical layers and total depth of total_depth whose bottom layer is (about) ratio times larger than the top layer.

The thickness profile transitions from the top-layer thickness to the bottom-layer thickness via a hyperbolic tangent proportional to tanh(2π * (k / (nlayers - 1) - 1 / 2)), where k = 0, 1, ..., nlayers - 1 is the layer index with k = 0 corresponding to the top-most layer.

The sum of all layer thicknesses is total_depth.

Positive parameter ratio prescribes (approximately) the ratio of the thickness of the bottom-most layer to the top-most layer. The final ratio of the bottom-most layer to the top-most layer ends up a bit different from the prescribed ratio. In particular, the final ratio of the bottom over the top-most layer thickness is (1 + ratio * exp(2π)) / (ratio + exp(2π)). This slight departure comes about because of the value of the hyperbolic tangent profile at the end-points tanh(π), which is approximately 0.9963 and not 1. Note that because exp(2π) is much greater than 1, the value of the actual ratio is not that different from the prescribed value ratio, e.g., for ratio values between 1/100 and 100 the final ratio of the bottom-most layer to the top-most layer only departs from the prescribed ratio by ±20%.

Parameters:

nlayers (int) – Number of vertical layers.
ratio (float) – The desired value of the ratio of bottom-most to the top-most layer thickness. Note that the final value of the ratio of bottom-most to the top-most layer thickness ends up (1 + ratio * exp(2π)) / (ratio + exp(2π)). Must be positive.
total_depth (float) – The total depth of grid, i.e., the sum of all thicknesses.

Returns:

An array containing the layer thicknesses.

Return type:

numpy.array

Examples

The spacings for a vertical grid with 20 layers, with maximum depth 1000 meters, and for which the top-most layer is about 4 times thinner than the bottom-most one.

>>> from regional_mom6 import hyperbolictan_thickness_profile
>>> nlayers, total_depth = 20, 1000
>>> ratio = 4
>>> dz = hyperbolictan_thickness_profile(nlayers, ratio, total_depth)
>>> dz
array([20.11183771, 20.2163053 , 20.41767549, 20.80399084, 21.53839043,
       22.91063751, 25.3939941 , 29.6384327 , 36.23006369, 45.08430684,
       54.91569316, 63.76993631, 70.3615673 , 74.6060059 , 77.08936249,
       78.46160957, 79.19600916, 79.58232451, 79.7836947 , 79.88816229])
>>> dz.sum()
1000.0
>>> dz[-1] / dz[0]
3.9721960481753706

If we want the top layer to be thicker then we need to prescribe ratio < 1.

>>> from regional_mom6 import hyperbolictan_thickness_profile
>>> nlayers, total_depth = 20, 1000
>>> ratio = 1/4
>>> dz = hyperbolictan_thickness_profile(nlayers, ratio, total_depth)
>>> dz
array([79.88816229, 79.7836947 , 79.58232451, 79.19600916, 78.46160957,
       77.08936249, 74.6060059 , 70.3615673 , 63.76993631, 54.91569316,
       45.08430684, 36.23006369, 29.6384327 , 25.3939941 , 22.91063751,
       21.53839043, 20.80399084, 20.41767549, 20.2163053 , 20.11183771])
>>> dz.sum()
1000.0
>>> dz[-1] / dz[0]
0.25174991059652

Now how about a grid with the same total depth as above but with equally-spaced layers.

>>> from regional_mom6 import hyperbolictan_thickness_profile
>>> nlayers, total_depth = 20, 1000
>>> ratio = 1
>>> dz = hyperbolictan_thickness_profile(nlayers, ratio, total_depth)
>>> dz
array([50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50.,
       50., 50., 50., 50., 50., 50., 50.])

regional_mom6.regional_mom6.longitude_slicer(data, longitude_extent, longitude_coords)#

Slice longitudes while handling periodicity and the ‘seams’, that is the longitude values where the data wraps around in a global domain (for example, longitudes are defined, usually, within domain [0, 360] or [-180, 180]).

The algorithm works in five steps:

Determine whether we need to add or subtract 360 to get the middle of the longitude_extent to lie within data’s longitude range (hereby old_lon).
Shift the dataset so that its midpoint matches the midpoint of longitude_extent (up to a multiple of 360). Now, the modified old_lon does not increase monotonically from West to East since the ‘seam’ has moved.
Fix old_lon to make it monotonically increasing again. This uses the information we have about the way the dataset was shifted/rolled.
Slice the data index-wise. We know that |longitude_extent[1] - longitude_extent[0]| / 360 multiplied by the number of discrete longitude points in the global input data gives the number of longitude points in our slice, and we’ve already set the midpoint to be the middle of the target domain.
Finally re-add the correct multiple of 360 so the whole domain matches the target.

Parameters:

data (xarray.Dataset) – The global data you want to slice in longitude.
longitude_extent (Tuple[float, float]) – The target longitudes (in degrees) we want to slice to. Must be in increasing order.
longitude_coords (Union[str, list[str]) – The name or list of names of the longitude coordinates(s) in data.

Returns:

The sliced data.

Return type:

xarray.Dataset

regional_mom6.regional_mom6.rectangular_hgrid(λ, φ)#

Construct a horizontal grid with all the metadata required by MOM6, based on arrays of longitudes (λ) and latitudes (φ) on the supergrid. Here, ‘supergrid’ refers to both cell edges and centres, meaning that there are twice as many points along each axis than for any individual field.

Caution

It is assumed the grid’s boundaries are lines of constant latitude and longitude. Rotated grids need to be handled differently.

It is also assumed here that the longitude array values are uniformly spaced.

Ensure both λ and φ are monotonically increasing.

Parameters:

λ (numpy.array) – All longitude points on the supergrid. Must be uniformly spaced.
φ (numpy.array) – All latitude points on the supergrid.

Returns:

An FMS-compatible horizontal grid (hgrid) that includes all required attributes.

Return type:

xarray.Dataset

class regional_mom6.regional_mom6.segment(*, hgrid, infile, outfolder, varnames, segment_name, orientation, startdate, arakawa_grid='A', time_units='days', tidal_constituents=None, repeat_year_forcing=False)#

Class to turn raw boundary segment data into MOM6 boundary segments.

Boundary segments should only contain the necessary data for that segment. No horizontal chunking is done here, so big fat segments will process slowly.

Data should be at daily temporal resolution, iterating upwards from the provided startdate. Function ignores the time metadata and puts it on Julian calendar.

Note

Only supports z-star (z*) vertical coordinate.

Parameters:

hgrid (xarray.Dataset) – The horizontal grid used for domain.
infile (Union[str, Path]) – Path to the raw, unprocessed boundary segment.
outfolder (Union[str, Path]) – Path to folder where the model inputs will be stored.
varnames (Dict[str, str]) – Mapping between the variable/dimension names and standard naming convention of this pipeline, e.g., {"xq": "longitude, "yh": "latitude", "salt": "salinity", ...}. Key “tracers” points to nested dictionary of tracers to include in boundary.
segment_name (str) – Name of the segment, e.g., 'segment_001'.
orientation (str) – Cardinal direction (lowercase) of the boundary segment, i.e., 'east', 'west', 'north', or 'south'.
startdate (str) – The starting date to use in the segment calendar.
arakawa_grid (Optional[str]) – Arakawa grid staggering type of the boundary forcing. Either 'A' (default), 'B', or 'C'.
time_units (str) – The units used by the raw forcing files, e.g., hours, days (default).
tidal_constituents (Optional[int]) – An integer determining the number of tidal constituents to be included from the list: M₂, S₂, N₂, K₂, K₁, O₂, P₁, Q₁, Mm, Mf, and M₄. For example, specifying 1 only includes M₂; specifying 2 includes M₂ and S₂, etc. Default: None.
repeat_year_forcing (Optional[bool]) – When True the experiment runs with repeat-year forcing. When False (default) then inter-annual forcing is used.

rectangular_brushcut()#: Cut out and interpolate tracers. rectangular_brushcut assumes that the boundary is a simple Northern, Southern, Eastern, or Western boundary.

`utils`#

regional_mom6.utils.angle_between(v1, v2, v3)#

Return the angle v2-v1-v3 (in radians), where v1, v2, v3 are 3-vectors. That is, the angle that is formed between vectors v2 - v1 and vector v3 - v1.

Example

>>> from regional_mom6.utils import angle_between
>>> v1 = (0, 0, 1)
>>> v2 = (1, 0, 0)
>>> v3 = (0, 1, 0)
>>> angle_between(v1, v2, v3)
1.5707963267948966
>>> from numpy import rad2deg
>>> rad2deg(angle_between(v1, v2, v3))
90.0

regional_mom6.utils.latlon_to_cartesian(lat, lon, R=1)#

Convert latitude and longitude (in degrees) to Cartesian coordinates on a sphere of radius R. By default R = 1.

Parameters:

lat (float) – Latitude (in degrees).
lon (float) – Longitude (in degrees).
R (float) – The radius of the sphere; default: 1.

Returns:

Tuple with the Cartesian coordinates x, y, z

Return type:

tuple

Examples

Find the Cartesian coordinates that correspond to point with (lat, lon) = (0, 0) on a sphere with unit radius.

>>> from regional_mom6.utils import latlon_to_cartesian
>>> latlon_to_cartesian(0, 0)
(1.0, 0.0, 0.0)

Now let’s do the same on a sphere with Earth’s radius

>>> from regional_mom6.utils import latlon_to_cartesian
>>> R = 6371e3
>>> latlon_to_cartesian(0, 0, R)
(6371000.0, 0.0, 0.0)

regional_mom6.utils.quadrilateral_area(v1, v2, v3, v4)#

Return the area of a spherical quadrilateral on the unit sphere that has vertices on the 3-vectors v1, v2, v3, v4 (counter-clockwise orientation is implied). The area is computed via the excess of the sum of the spherical angles of the quadrilateral from 2π.

Example

Calculate the area that corresponds to half the Northern hemisphere of a sphere of radius R. This should be 1/4 of the sphere’s total area, that is π R².

>>> from regional_mom6.utils import quadrilateral_area, latlon_to_cartesian
>>> R = 434.3
>>> v1 = latlon_to_cartesian(0, 0, R)
>>> v2 = latlon_to_cartesian(0, 90, R)
>>> v3 = latlon_to_cartesian(90, 0, R)
>>> v4 = latlon_to_cartesian(0, -90, R)
>>> quadrilateral_area(v1, v2, v3, v4)
592556.1793298927
>>> from numpy import pi
>>> quadrilateral_area(v1, v2, v3, v4) == pi * R**2
True

regional_mom6.utils.quadrilateral_areas(lat, lon, R=1)#

Return the area of spherical quadrilaterals on a sphere of radius R. By default, R = 1. The quadrilaterals are formed by constant latitude and longitude lines on the lat-lon grid provided.

Parameters:

lat (numpy.array) – Array of latitude points (in degrees).
lon (numpy.array) – Array of longitude points (in degrees).
R (float) – The radius of the sphere; default: 1.

Returns:

Array with the areas of the quadrilaterals defined by the lat-lon grid provided. If the provided lat, lon arrays are of dimension m \(\times\) n then returned areas array is of dimension (m-1) \(\times\) (n-1).

Return type:

numpy.array

Example

Let’s construct a lat-lon grid on the sphere with 60 degree spacing. Then we compute the areas of each grid cell and confirm that the sum of the areas gives us the total area of the sphere.

>>> from regional_mom6.utils import quadrilateral_areas
>>> import numpy as np
>>> λ = np.linspace(0, 360, 7)
>>> φ = np.linspace(-90, 90, 4)
>>> lon, lat = np.meshgrid(λ, φ)
>>> lon
array([[  0.,  60., 120., 180., 240., 300., 360.],
       [  0.,  60., 120., 180., 240., 300., 360.],
       [  0.,  60., 120., 180., 240., 300., 360.],
       [  0.,  60., 120., 180., 240., 300., 360.]])
>>> lat
array([[-90., -90., -90., -90., -90., -90., -90.],
       [-30., -30., -30., -30., -30., -30., -30.],
       [ 30.,  30.,  30.,  30.,  30.,  30.,  30.],
       [ 90.,  90.,  90.,  90.,  90.,  90.,  90.]])
>>> R = 6371e3
>>> areas = quadrilateral_areas(lat, lon, R)
>>> areas
array([[1.96911611e+13, 1.96911611e+13, 1.96911611e+13, 1.96911611e+13,
        1.96911611e+13, 1.96911611e+13],
       [4.56284230e+13, 4.56284230e+13, 4.56284230e+13, 4.56284230e+13,
        4.56284230e+13, 4.56284230e+13],
       [1.96911611e+13, 1.96911611e+13, 1.96911611e+13, 1.96911611e+13,
        1.96911611e+13, 1.96911611e+13]])
>>> np.isclose(areas.sum(), 4 * np.pi * R**2, atol=np.finfo(areas.dtype).eps)
True

regional_mom6.utils.vecdot(v1, v2)#: Return the dot product of vectors v1 and v2. v1 and v2 can be either numpy vectors or numpy.ndarrays in which case the last dimension is considered the dimension over which the dot product is taken.

API reference

Contents

API reference#

regional_mom6#

utils#

`regional_mom6`#

`utils`#