The InternalModes class can be used to quickly and accurately compute the vertical modes and SQG modes from arbitrary stratification.
The complete class hierarchy contains several implementations that include spectral methods, finite differencing, WKB approximated solutions, as well as the analytical solutions for constant stratification and exponential stratification. The details are documented in Early, Lelong, and Smith (2018).
The classes contain many options, but also try to remain simple to use. Use
help InternalModes in Matlab for a complete description, or use the Quick Start below.
If you use these classes to compute the vertical modes, please cite the following paper,
- J. Early, M. P. Lelong, K. S. Smith. Fast and accurate computation of vertical modes using Chebyshev polynomials, Ocean Modelling, In prep.
Table of contents
- Quick Start
- Convenience functions
- Constant and exponential stratification
- Built-in density profiles
- Normalization and boundary conditions
- Name/value pairs
- Numerical methods
- Class hierarchy
- Internal modes class cluster
- Unit testing
The preferred method for initializing the InternalModes class is to define density as a function, e.g.,
N0 = 5.2e-3; % reference buoyancy frequency, radians/seconds g = 9.81; rho0 = 1025; % density at the surface b = 1300; rho = @(z) rho0*(1 + b*N0*N0/(2*g)*(1 - exp(2*z/b)));
and then use that function to initialize the class,
L = -5000; zOut = linspace(L,0,200)'; latitude = 33; im = InternalModes(rho,[L 0],zOut,latitude);
InternalModes class is initialized with four arguments: the density function, an array specifying the domain bounds (lower and upper boundary), the the output grid upon which all returned functions will be given, and the latitude.
Now that the
im object is initialized, you can request the internal modes at a given wavenumber,
k, where k is 2*pi/wavelength.
[F,G,h,omega] = im.ModesAtWavenumber(2*pi/1000);
[F,G,h,k] = im.ModesAtFrequency(5*im.f0);
G contain the vertical modes for u/v and w/rho, respectively. The arrays have dimensions
size(F)=[length(zOut) length(h)], meaning that each column
i is a normal mode,
F(:,i) with corresponding eigendepth
h(i). The associated eigenfrequency or eigenwavenumber are also returned for convenience.
You can also request the surface quasigeostrophy (SQG) modes,
psi_t = im.SurfaceModesAtWavenumber(2*pi/1000); psi_b = im.BottomModesAtWavenumber(2*pi/1000);
The functions can also take arrays of wavenumbers of any shape or size.
Note that you can also request variations of the density, e.g.,
N2 = im.N2; rho_zz = im.rho_zz;
which will be on the same output grid
zOut that you specified.
Although it is best to use initialize using a density function, you can also initialize using gridded data if needed. The syntax is nearly the same,
N0 = 5.2e-3; % reference buoyancy frequency, radians/seconds g = 9.81; rho0 = 1025; % density at the surface b = 1300; z = linspace(-5000,0,500)'; rho = rho0*(1 + b*N0*N0/(2*g)*(1 - exp(2*z/b))); im = InternalModes(rho,z,zOut,latitude);
where now you pass the gridded density data,
rho, and its coordinate,
InternalModes objects is initialized there are a few notable convenience functions. You can use,
to quickly visualize the four lowest modes.
Constant and exponential stratification
Both constant and exponential stratification profiles have exact analytical solutions that have been implemented directly. Initialization follows the same basic format, but rather than specify a density profile, you specify the parameters for the stratification. For example,
N0 = 5.2e-3; imConst = InternalModesConstantStratification(N0,[L 0],zOut,latitude);
will initialize a constant stratification profile with buoyancy frequency given by
N0 = 5.2e-3; b = 1300; imExp = InternalModesExponentialStratification([N0 b],[L 0],zOut,latitude);
will initialize an exponential stratification profile with surface buoyancy frequency of
N0 and e-fold scale of
Built-in density profiles
For testing purposes and convenience, there are a number of pre-defined density profiles that you can access. Simply call,
[rho,N2,zIn] = InternalModes.StratificationProfileWithName(stratification)
where the variable
stratification is a string. The returned values
zIn can be used directly as the first two arguments in initialization. The following options are available:
'constant'A constant stratification profile.
'exponential'The standard Garrett-Munk exponential profile.
'pycnocline-constant'Constant stratification with a pycnocline, taken from Cushman-Roisin & Beckers.
'pycnocline-exponential'Exponential profile with a deep pycnocline.
'latmix-site1'Attempts to recreate the full stratification profile (down to 5000 meters) at Latmix site 1, with an intense surface mixed layer and deep pycnocline.
'latmix-site1-surface'Recreates just the near surface features at Latmix Site 1.
'latmix-site1-exponential'Recreates the surface mixed layer and associated pycnocline, and then decays exponentially below 190m.
'latmix-site1-constant'Recreates the surface mixed layer and associated pycnocline, but then goes constant below 300m.
Normalization and boundary conditions
There are four choices for normalization of the vertical modes. See in Early, et al. (2018) for details.
kConstant(default) Uses the K-constant norm, which is useful for specifying energy density of an internal wave.
omegaConstantUses the omega-constant norm, the utility of which isn't exactly clear.
uMaxSets the maximum of the F modes to 1. Useful if you want to specify the horizontal velocity of a wave.
wMaxSets the maximum of the G modes to 1. Useful if you want to specify the vertical velocity of a wave.
These constants are defined as an
enumeration class, and would be called as,
im.normalization = Normalization.uMax
before you request the modes.
The upper boundary condition can be specified as a rigid lid (default) or as a free surface (there is no reduced-gravity option yet). For example,
im.upperBoundary = UpperBoundary.freeSurface
will use the free surface condition at the upper boundary when solving the EVP.
During initialization you can specify options and set properties using a name/value pair syntax. For example, to specify the upper boundary condition during initialization, you could write
im = InternalModes(rho,[L 0],zOut,latitude,'upperBoundary',UpperBoundary.freeSurface);
You can add as many name/value pairs as needed.
There are multiple different numerical methods for solving the eigenvalue problem (EVP), some of which may work better than others, depending on the density profile in question. During initialization, you can specify which numerical method you want to use to solve the EVP.
For example, to initialize using finite differencing, you can call,
imFD = InternalModes(rho,[L 0],zOut,latitude,'method','finiteDifference');
The following different methods are available
'spectral'Uses Chebyshev polynomials on a z coordinate.
'wkbSpectral'Uses Chebyshev polynomials on a WKB stretched coordinate.
'densitySpectral'Uses Chebyshev polynomials on a density coordinate.
'wkbAdaptiveSectral'(default) Uses coupled EVPs with Chebyshev polynomials on a WKB stretched coordinate, depending on the frequency requested.
'finiteDifference'Uses finite differencing techniques.
'wkb'Uses the analytical WKB solutions (only valid for frequency, and at most one turning point).
'wkb-hydrostatic'Uses a hydrostatic form of the WKB solution as seen in Levine (2002).
Some of these methods take additional options.
All of the spectral methods have two other properties that can be specified on initialization,
nGridspecifies the number of points that should be used to compute the derivatives of the density function (which is always on the z coordinate). These computations are done in O(n log n) time, so a relatively large number of points can be used here. The default is 2^14+1.
nEVPspecifies the number of points used in the eigenvalue problem on the stretched coordinate (either z, density or wkb). This is a relatively expensive operation which scales as O(n^3). The default is 513 points.
The finite difference matrices are always computed on the given input grid, and the results are interpolated onto the output grid. However, you can specify the accuracy of the finite difference matrices.
orderOfAccuracyThe lowest appropriate value is 2 (which is the usual approach taken), and the default is 4.
InternalModes class is actually a class cluster, that is used to initialize any one of a number of concrete subclasses which inherit from the
When you initialize the
InternalModes class, internally it then initializes one of the concrete subclasses and simply passes all (most) calls to the concrete subclasses.
This means that you can choose to initialize one of the concrete subclasses directly, and simply bypass the
InternalModes class cluster.
Internal modes class cluster
InternalModes class cluster contains several methods that are not in the
Quick Initialization and Test Cases
The class cluster can be initialized quickly for a test cases with the following syntax,
im = InternalModes(stratification, method, nPoints)
where the defaults (if not given) are
If you initialize with this test case method, you can then call either
to estimate the error.
UnitTests constains two functions that use the analytical solutions from constant and exponential stratification to unit test against the numerical methods.
InternalModeRelativeErrorTests.m computes the internal mode solution at three different wavenumbers and three different frequencies, returning the number of 'good' modes (modes with less than 1 percent error). The test should be run with both surface boundary conditions.
SQGModeRelativeErrorTests.m computes the surface and bottom boundary SQG mode solutions at a range wavenumbers, returning the relative error.