11. Utilities

11.1. MITgcmutils

This Python package includes a number of helpful functions and scripts for dealing with MITgcm output. You can install it from the model repository (in directory utils/python/MITgcmutils) or from the Python Package Index:

pip install --user MITgcmutils

The following functions are exposed at the package level:

The package also includes a standalone script for joining tiled mnc files: gluemncbig.

For more functions, see the individual modules:

11.1.1. mds

exception MITgcmutils.mds.ParseError[source]
MITgcmutils.mds.parsemeta(metafile)[source]

parses metafile (file object or filename) into a dictionary of lists of floats, ints or strings

MITgcmutils.mds.rdmds(fnamearg, itrs=-1, machineformat='b', rec=None, fill_value=0, returnmeta=False, astype=<class 'float'>, region=None, lev=(), usememmap=False, mm=False, squeeze=True, verbose=False)[source]

Read meta-data files as written by MITgcm.

Call signatures:

a = rdmds(fname,…)

a,its,meta = rdmds(fname,…,returnmeta=True)

Parameters
  • fname (string) –

    name of file to read, without the ‘.data’ or ‘.meta’ suffix. If itrs is given, the iteration number is added to fname as well. fname may contain shell wildcards, which is useful for tile files organized into directories, e.g.,

    T = rdmds(‘prefix*/T’, 2880)

    will read prefix0000/T.0000002880.*, prefix0001/T.0000002880.*, … (and any others that match the wildcard, so be careful how you name things!)

  • itrs (int or list of ints or np.NaN or np.Inf) –

    Iteration number(s). With itrs=-1, will try to read

    fname.meta or fname.001.001.meta, …

    If itrs is a list of integers of an integer, it will read the corresponding

    fname.000000iter.meta, …

    If itrs is np.NaN, it will read all iterations for which files are found. If itrs is np.Inf, it will read the highest iteration found.

  • machineformat (int) – endianness (‘b’ or ‘l’, default ‘b’)

  • rec (list of int or None) – list of records to read (default all) useful for pickups and multi-field diagnostics files

  • fill_value (float) – fill value for missing (blank) tiles (default 0)

  • astype (data type) – data type to return (default: double precision) None: keep data type/precision of file

  • region (tuple of int) – (x0,x1,y0,y1) read only this region (default (0,nx,0,ny))

  • lev (list of int or tuple of lists of int) – list of levels to read, or, for multiple dimensions (excluding x,y), tuple(!) of lists (see examples below)

  • usememmap (bool) – if True, use a memory map for reading data (default False) recommended when using lev, or region with global files to save memory and, possibly, time

Returns

  • a (array_like) – numpy array of the data read

  • its (list of int) – list of iteration numbers read (only if returnmeta=True)

  • meta (dict) – dictionary of metadata (only if returnmeta=True)

Examples

>>> XC = rdmds('XC')
>>> XC = rdmds('res_*/XC')
>>> T = rdmds('T.0000002880')
>>> T = rdmds('T',2880)
>>> T2 = rdmds('T',[2880,5760])
>>> T,its = rdmds('T',numpy.Inf)
>>> VVEL = rdmds('pickup',2880,rec=range(50,100))
>>> a5 = rdmds('diags',2880,rec=0,lev=[5])
>>> a = rdmds('diags',2880,rec=0,lev=([0],[0,1,5,6,7]))
>>> from numpy import r_
>>> a = rdmds('diags',2880,rec=0,lev=([0],r_[:2,5:8]))  # same as previous
>>> a = rdmds('diags',2880,rec=0)[0, [0,1,5,6,7], ...]  # same, but less efficient
>>> a = rdmds('diags',2880)[0, 0, [0,1,5,6,7], ...]     # even less efficient
MITgcmutils.mds.readmeta(f)[source]

read meta file and extract tile/timestep-specific parameters

MITgcmutils.mds.scanforfiles(fname)[source]

return list of iteration numbers for which metafiles with base fname exist

MITgcmutils.mds.strip_comments(text)[source]

strips C and C++ style comments from text

MITgcmutils.mds.wrmds(fbase, arr, itr=None, dataprec='float32', ndims=None, nrecords=None, times=None, fields=None, simulation=None, machineformat='b', deltat=None, dimlist=None)[source]

Write an array to an mds meta/data file set.

If itr is given, the files will be named fbase.0000000itr.data and fbase.0000000itr.meta, otherwise just fbase.data and fbase.meta.

Parameters
  • fbase (string) – Name of file to write, without the ‘.data’ or ‘.meta’ suffixes, and without the iteration number if itr is give

  • arr (array_like) – Numpy array to write

  • itr (int or None) – If given, this iteration number will be appended to the file name

  • dataprec (string) – precision of resulting file (‘float32’ or ‘float64’)

  • ndims (int) – number of non-record dimensions; extra (leading) dimensions will be folded into 1 record dimension

  • nrecords (int) – number of records; will fold as many leading dimensions as necessary (has to match shape!)

  • times (float or list of floats) – times to write into meta file. Either a single float or a list of two for a time interval

  • fields (list of strings) – list of fields

  • simulation (string) – string describing the simulation

  • machineformat (string) – ‘b’ or ‘l’ for big or little endian

  • deltat (float) – time step; provide in place of either times or itr to have one computed from the other

  • dimlist (tuple) – dimensions as will be stored in file (only useful when passing meta data from an existing file to wrmds as keyword args)

11.1.2. mnc

class MITgcmutils.mnc.MNC(fpatt, layout=None, multitime=False)[source]

A file object for MNC (tiled NetCDF) data.

Should behave mostly like scipy.io.netcdf.netcdf_file in ‘r’ mode.

Parameters
  • fpatt (string) – glob pattern for tile files

  • layout (string) –

    which global layout to use:

    ’model’

    use layout implied by Nx, Ny

    ’exch2’

    use exch2 global layout

    ’faces’

    variables are lists of exch2 faces

    default is to use exch2 layout if present, model otherwise

Example

>>> nc = mnc_files('mnc_*/state.0000000000.t*.nc')
>>> temp = nc.variables['Temp'][:]
>>> salt = nv.variables['S'][:]
>>> nc.close()
temp and salt are now assembled (global) arrays of shape (Nt, Nr, Ny, Nx)
where Nt is the number iterations found in the file (in this case probably 1).

Notes

The multitime option is not implemented, i.e., MNC cannot read files split in time.

close()[source]

Close tile files

MITgcmutils.mnc.mnc_files(fpatt, layout=None)[source]

A file object for MNC (tiled NetCDF) data.

Should behave mostly like scipy.io.netcdf.netcdf_file in ‘r’ mode.

Parameters
  • fpatt (string) – glob pattern for tile files

  • layout (string) –

    which global layout to use:

    ’model’

    use layout implied by Nx, Ny

    ’exch2’

    use exch2 global layout

    ’faces’

    variables are lists of exch2 faces

    default is to use exch2 layout if present, model otherwise

Example

>>> nc = mnc_files('mnc_*/state.0000000000.t*.nc')
>>> temp = nc.variables['Temp'][:]
>>> salt = nv.variables['S'][:]
>>> nc.close()
temp and salt are now assembled (global) arrays of shape (Nt, Nr, Ny, Nx)
where Nt is the number iterations found in the file (in this case probably 1).

Notes

The multitime option is not implemented, i.e., MNC cannot read files split in time.

MITgcmutils.mnc.rdmnc(fpatt, varnames=None, iters=None, slices=Ellipsis, layout=None)[source]

Read one or more variables from an mnc file set.

Parameters
  • fpatt (string) – glob pattern for netcdf files comprising the set

  • varnames (list of strings, optional) – list of variables to read (default all)

  • iters (list of int, optional) – list of iterations (not time) to read

  • slices (tuple of slice objects) – tuple of slices to read from each variable (typically given as numpy.s_[…])

Returns

dictionary of variable arrays

Return type

dict of numpy arrays

Example

>>> S = rdmnc("mnc_*/state.0000000000.*', ['U', 'V'], slices=numpy.s_[..., 10:-10, 10:-10])
>>> u = S['U']
>>> v = S['V']

Notes

Can currently read only one file set (i.e., 1 file per tile), not several files split in time.

Consider using mnc_files for more control (and similar convenience). The same restriction about multiple files applies, however.

11.1.3. diagnostics

MITgcmutils.diagnostics.readstats(fname)[source]

locals,totals,itrs = readstats(fname)

Read a diagstats text file into record arrays (or dictionaries).

Parameters

fname (string) – name of diagstats file to read

Returns

  • locals (record array or dict of arrays) – local statistics, shape (len(itrs), Nr, 5)

  • totals (record array or dict of arrays) – column integrals, shape (len(itrs), 5)

  • itrs (list of int) – iteration numbers found in the file

Notes

  • The 5 columns of the resulting arrays are average, std.dev, min, max and total volume.

  • There is a record (or dictionary key) for each field found in the file.

11.1.4. ptracers

MITgcmutils.ptracers.iolabel(i)[source]

Map tracer number (1..3843) to 2-character I/O label:

1..99 => 01..99
100..619 => 0a..0Z,1a..9Z
620..3843 => aa..ZZ
Parameters

i (int) – ptracer number (1..3843)

Returns

2-character I/O label

Return type

string

MITgcmutils.ptracers.iolabel2num(s)[source]

Map 2-character IO label to tracer number, the inverse of iolabel()

11.1.5. jmd95

Density of Sea Water using the Jackett and McDougall 1995 (JAOT 12) polynomial

MITgcmutils.jmd95.bulkmodjmd95(s, theta, p)[source]

Compute bulk modulus

MITgcmutils.jmd95.dens(s, theta, p)

Computes in-situ density of sea water

Density of Sea Water using Jackett and McDougall 1995 (JAOT 12) polynomial (modified UNESCO polynomial).

Parameters
  • s (array_like) – salinity [psu (PSS-78)]

  • theta (array_like) – potential temperature [degree C (IPTS-68)]; same shape as s

  • p (array_like) – pressure [dbar]; broadcastable to shape of s

Returns

dens – density [kg/m^3]

Return type

array

Example

>>> densjmd95(35.5, 3., 3000.)
1041.83267

Notes

AUTHOR: Martin Losch 2002-08-09 (mlosch@mit.edu)

Jackett and McDougall, 1995, JAOT 12(4), pp. 381-388

MITgcmutils.jmd95.densjmd95(s, theta, p)[source]

Computes in-situ density of sea water

Density of Sea Water using Jackett and McDougall 1995 (JAOT 12) polynomial (modified UNESCO polynomial).

Parameters
  • s (array_like) – salinity [psu (PSS-78)]

  • theta (array_like) – potential temperature [degree C (IPTS-68)]; same shape as s

  • p (array_like) – pressure [dbar]; broadcastable to shape of s

Returns

dens – density [kg/m^3]

Return type

array

Example

>>> densjmd95(35.5, 3., 3000.)
1041.83267

Notes

AUTHOR: Martin Losch 2002-08-09 (mlosch@mit.edu)

Jackett and McDougall, 1995, JAOT 12(4), pp. 381-388

11.1.6. mdjwf

Density of Sea Water using McDougall et al. 2003 (JAOT 20) polynomial

MITgcmutils.mdjwf.dens(s, theta, p)

Computes in-situ density of sea water

Density of Sea Water using McDougall et al. 2003 (JAOT 20) polynomial (Gibbs Potential).

Parameters
  • s (array_like) – salinity [psu (PSS-78)]

  • theta (array_like) – potential temperature [degree C (IPTS-68)]; same shape as s

  • p (array_like) – pressure [dbar]; broadcastable to shape of s

Returns

dens – density [kg/m^3]

Return type

array

Example

>>> densmdjwf(35., 25., 2000.)
1031.654229

Notes

AUTHOR: Martin Losch 2002-08-09 (Martin.Losch@awi.de)

McDougall et al., 2003, JAOT 20(5), pp. 730-741

MITgcmutils.mdjwf.densmdjwf(s, theta, p)[source]

Computes in-situ density of sea water

Density of Sea Water using McDougall et al. 2003 (JAOT 20) polynomial (Gibbs Potential).

Parameters
  • s (array_like) – salinity [psu (PSS-78)]

  • theta (array_like) – potential temperature [degree C (IPTS-68)]; same shape as s

  • p (array_like) – pressure [dbar]; broadcastable to shape of s

Returns

dens – density [kg/m^3]

Return type

array

Example

>>> densmdjwf(35., 25., 2000.)
1031.654229

Notes

AUTHOR: Martin Losch 2002-08-09 (Martin.Losch@awi.de)

McDougall et al., 2003, JAOT 20(5), pp. 730-741

11.1.7. cs

MITgcmutils.cs.pcol(x, y, data, projection=None, vmin=None, vmax=None, **kwargs)[source]

Plots 2D scalar fields on the MITgcm cubed sphere grid with pcolormesh.

Parameters
  • x (array_like) – ‘xg’, that is, x coordinate of the points one half grid cell to the left and bottom, that is vorticity points for tracers, etc.

  • y (array_like) – ‘yg’, that is, y coordinate of same points

  • data (array_like) – scalar field at tracer points

  • projection (Basemap instance, optional) – used to transform if present. Unfortunatly, cylindrical and conic maps are limited to the [-180 180] range. projection = ‘sphere’ results in a 3D visualization on the sphere without any specific projection. Good for debugging.

Example

>>> from mpl_toolkits.basemap import Basemap
>>> import MITgcmutils as mit
>>> import matplotlib.pyplot as plt
>>> from sq import sq
>>>
>>> x=mit.rdmds('XG'); y=mit.rdmds('YG'); e=mit.rdmds('Eta',np.Inf)
>>> fig = plt.figure();
>>> mp = Basemap(projection='moll',lon_0 = 0.,
>>>              resolution = 'l', area_thresh = 1000.)
>>> plt.clf()
>>> h = mit.cs.pcol(x,y,sq(e), projection = mp)
>>> mp.fillcontinents(color = 'grey')
>>> mp.drawmapboundary()
>>> mp.drawmeridians(np.arange(0, 360, 30))
>>> mp.drawparallels(np.arange(-90, 90, 30))
>>> plt.show()

11.1.8. llc

MITgcmutils.llc.contour(*arguments, **kwargs)[source]

Create a contour plot of a 2-D llc array (with tricontour).

Call signatures:

contour(X, Y, C, N, **kwargs)

contour(X, Y, C, V, **kwargs)
Parameters
  • X (array-like) – x coordinates of the grid points

  • Y (array-like) – y coordinates of the grid points

  • C (array-like) – array of color values.

  • N (int) – number of levels

  • V (list of float) – list of levels

  • kwargs – passed to tricontour.

MITgcmutils.llc.contourf(*arguments, **kwargs)[source]

Create a contourf plot of a 2-D llc array (with tricontour).

Call signatures:

contourf(X, Y, C, N, **kwargs)

contourf(X, Y, C, V, **kwargs)
Parameters
  • X (array-like) – x coordinates of the grid points

  • Y (array-like) – y coordinates of the grid points

  • C (array-like) – array of color values.

  • N (int) – number of levels

  • V (list of float) – list of levels

  • kwargs – passed to tricontour.

MITgcmutils.llc.div(u, v, dxg=None, dyg=None, rac=None, hfw=None, hfs=None)[source]

Compute divergence of vector field (U,V) on llc grid

Call signatures:

divergence = div(U, V, DXG, DYG, RAC, HFW, HFS)
divergence = div(U, V)
divergence = div(U, V, DXG, DYG)
divergence = div(U, V, DXG, DYG, RAC)
divergence = div(U, V, DXG, DYG, hfw=HFW, hfs=HFS)
Parameters
  • u (array-like (timelevel,depthlevel,jpoint,ipoint)) – x-component of vector field at u-point

  • v (array-like (timelevel,depthlevel,jpoint,ipoint)) – y-component of vector field at v-point

  • dxg (array-like (jpoint,ipoint), optional) – grid spacing in x across v-point, defaults to one

  • dyg (array-like (jpoint,ipoint), optional) – grid spacing in y across u-point, defaults to one

  • rac (array-like (jpoint,ipoint), optional) – grid cell area, defaults to dxg*dyg

  • hfw (array-like (depthlevel,jpoint,ipoint), optional) – hFac at u-point, defaults to one

  • hfs (array-like (depthlevel,jpoint,ipoint), optional) – hFac at v-point, defaults to one

MITgcmutils.llc.faces(fld)[source]

convert mds multidimensional data into a list with 6 faces

MITgcmutils.llc.faces2mds(ff)[source]

convert 6 faces to mds 2D data, inverse opertation of llc.faces

MITgcmutils.llc.flat(fld, **kwargs)[source]

convert mds data into global 2D field only fields with 2 to 5 dimensions are allowed

MITgcmutils.llc.grad(X, dxc=None, dyc=None, hfw=None, hfs=None)[source]

Compute horizontal gradient of scalar field X on llc grid

Call signatures:

dXdx, dXdy = div(X, DXC, DYC, HFW, HFS)
dXdx, dXdy = div(X)
dXdx, dXdy = div(X, DXC, DYC)
dXdx, dXdy = div(X, hfw=HFW, hfs=HFS)
Parameters
  • X (array-like (timelevel,depthlevel,jpoint,ipoint)) – scalar field at c-point

  • dxc (array-like (jpoint,ipoint), optional) – grid spacing in x across u-point, defaults to one

  • dyc (array-like (jpoint,ipoint), optional) – grid spacing in y across v-point, defaults to one

  • hfw (array-like (depthlevel,jpoint,ipoint), optional) – hFac at u-point, defaults to one

  • hfs (array-like (depthlevel,jpoint,ipoint), optional) – hFac at v-point, defaults to one

MITgcmutils.llc.mds(fld, center='Atlantic')[source]

convert global ‘flat’ field into mds data; only fields with 2 to 5 dimensions are allowed

MITgcmutils.llc.pcol(*arguments, **kwargs)[source]

Create a pseudo-color plot of a 2-D llc array (with plt.pcolormesh).

Call signatures:

pcol(X, Y, C, **kwargs)

pcol(X, Y, C, m, **kwargs)
Parameters
  • X (array-like) – x coordinates of the grid point corners (G-points)

  • Y (array-like) – y coordinates of the grid point corners (G-points)

  • C (array-like) – array of color values.

  • m (Basemap instance, optional) – map projection to use. NOTE: currently not all projections work

  • kwargs – passed to plt.pcolormesh.

MITgcmutils.llc.uv2c(u, v)[source]

Average vector component (u,v) to center points on llc grid

Call signatures:

uc,vc = uv2c(U,V)
Parameters
  • U (array-like (timelevel,depthlevel,jpoint,ipoint)) – x-component of vector field at u-point

  • V (array-like (timelevel,depthlevel,jpoint,ipoint)) – y-component of vector field at v-point

11.1.9. gluemncbig

This command line script is part of MITgcmutils and provides a convenient method for stitching together NetCDF files into a single file covering the model domain. Be careful though - the resulting files can get very large.

Usage: gluemncbig [-2] [-q] [--verbose] [--help] [--many] [-v <vars>] -o <outfile> <files>

 -v <vars>  comma-separated list of variable names or glob patterns
 -2         write a NetCDF version 2 (64-Bit Offset) file allowing for large records
 --many     many tiles: assemble only along x in memory; less efficient
            on some filesystems, but opens fewer files simultaneously and
            uses less memory
 -q         suppress progress messages
 --verbose  report variables
 --help     show this help text

All files must have the same variables.
Each variable (or 1 record of it) must fit in memory.
With --many, only a row of tiles along x must fit in memory.

Examples:

gluemncbig -o ptr.nc mnc_*/ptr_tave.*.nc
gluemncbig -o BIO.nc -v 'BIO_*' mnc_*/ptr_tave.*.nc