Stable coastal-trapped waves with stratification, topography and mean flow in Matlab

This set of Matlab mfiles (all having names that begin with “bigr”) can be used to calculate stable, inviscid coastal-trapped wave modal structures and dispersion curves under very general circumstances. Only a real frequency is allowed, so that instability and damping cannot be accounted for directly, but computations are more efficient than for the general case, long-wave parameters can be computed for first order wave equation calculations (see Brink, 1989), and a more general perturbation decay time (Brink, 1990) can also be obtained. Modal structures and energy diagnostics are provided. Generally speaking, the code is only useful for subinertial wave frequencies.

bigromp.m This is the file that actually calls all the others and so does the calculations. Call as: >> bigromp(datnew) bigrfinch.m This is used to modify a pre-existing input array. It will give you a menu for what you want to change. Call as: For examples of usage, see the screen displays copied at the end of this document.

The Problem:
We seek to find the free wave solutions for linearized coastal-trapped waves along a straight coast. The depth h(x) is uniform in the alongshore (y) direction. The offshore direction is x (x = 0 is the coast) and the vertical direction is z (z = -h(x) is the bottom), where the origin is at the surface at the coast. Dissipative effects are limited to the infinitesmally thin bottom boundary, and are only included for calculating coefficients for use with the first order wave equation.
In all cases, subscripted independent variables represent partial differentiation. The mean velocity v0(x, z) is in geostrophic balance with the mean pressure p0 and density ρ0 + ρ1. Turbulent stresses in the x, y directions are given by τ x , τ y , and the bottom stress ( τ0 x , τ0 y ) is related to the interior bottom velocity by τ0 x = ρr uB, τ0 y = ρ r vB, (2a, 2b) where r(x) is a bottom resistance coefficient having units of velocity, and a subscript "B" denotes a quantity evaluated at the bottom. The turbulent stresses are not used for computing the modal structures, but they are used to compute the damping coefficients after the inviscid wave mode is known. We also define N 2 = -gρ0 -1 ρ1z and (3a) The problem is solved by assuming that where the frequency ω and the alongshore wavenumber l are real. The "primes" are dropped henceforth. Using this assumption, the problem can be reduced to a single partial differential equation for pressure This equation and the notation are similar to those of Mooers (1975 a, b), except that these are for hydrostatic conditions.
As noted by Brink (1981) and Dale et al. (2001), this equation admits a spurious solution at ω' = f when v0z = v0x = 0. The presence of this mode also distorts modal structures and dispersion curves for nearby frequencies. Typically, this spurious pressure mode is vertically uniform and it decays only slowly with distance offshore.
The boundary conditions are as follows: where γ = 1 for a free surface condition and γ= 0 for a rigid lid condition. In terms of pressure, where the Ekman transport is given as and the subscript "B" denotes a quantity evaluated at the bottom (along z = h(x) ). This complete form is used only in computing the long wave frictional coefficients once the inviscid mode is known. For the computing the inviscid modal structure, the boundary condition is: In terms of pressure, There is a choice of two boundary conditions at the coast.

Closed boundary condition:
This boundary condition is used when icbc = 0: Since there is no bottom friction, this boundary condition is sensible regardless of the coastal wall height. Damping coefficient calculations do account for Ekman transport here.
Stated in terms of pressure, this inviscid condition is: "Open" boundary condition: To mimic an open boundary (icbc = 1), we use which does a reasonable job for linear coastal-trapped waves at an "open" offshore boundary, where variations should have decayed (see Brink, 1982). To help make sure this is used reasonably, use of this condition requires that the bottom be flat at x = 0.
Stated in terms of pressure, this condition is: This boundary condition can only be expected to work well at subinertial frequencies. At higher frequencies, it can lead to internal wave reflection, rather than the free exit of internal waves.

x = L:
The offshore open boundary condition is the same as the coastal "open" position, and is called by setting iobc = 1. This seems to work well as long as the grid boundary (x = L) occurs far offshore of where the modal solutions have substantial structure. Further, problems can occur at superinertial frequencies. The condition, as used, requires that the bottom is flat at the offshore boundary. The condition is: or, in terms of pressure: In addition, using iobc = 0 allows use of a closed offshore boundary condition (the same as for the coastal boundary).
At this point, the problem is converted into stretched coordinates that map the problem into a rectangular domain (e.g. Wang and Mooers 1976). Specifically, The problem is then solved numerically via resonance iteration. An arbitrary forcing is applied through the surface boundary condition (essentially, a localized wind stress curl near the grid center), and ω is searched to find a real frequency corresponding to a given wavenumber l. This search makes use of the Matlab function "fminsearch.m".

Using the software:
Examples of what your Matlab screen will look like as you run the different mfiles are shown at the end of this document. The main routine, "bigromp", requires an input array that can be generated by "bigrsetup" and modified by "bigrfinch". When you run "bigrsetup" or "bigrfinch", you simply need to answer the questions posed. Do not forget to use square brackets when entering more than one number on a line.
When running "bigromp", you will first be shown a few simple parameter values. You also get a section plot showing topography, stratification (minus a constant background value), and the mean flow. If you are allowing pauses during the run, you are advised of pauses, and need to hit a key to continue. At this point, the alongshore wavenumber is displayed, and the program displays the frequency and the inverse resonance parameter (rrr = inverse of the integral of |p| 2 over x and z) for each ω estimate. When you are at a resonance (free mode), rrr becomes very small. You can thus watch the output converge onto a complex free modal frequency. Convergence to your given tolerance is then announced, and the free mode pressure is plotted as a section. A pause is announced at this point. Hit any key to continue, and you will then also get section plots of u, v, w, ρ, energy diagnostics, and line plots of p, v, u at the surface. If you set the program to run more than one frequency, it will then display the next wavenumber and go on from there.
If you have set the program to compute more than one frequency, wavenumber pair, the program will also draw a dispersion curve once all pairs have been obtained.

Some Notes on using the software:
The x = 0 and x = xmax boundary conditions can either be closed (as at a coast), or open. Since there is no bottom friction, the coastal wall can have any depth.
If the x = 0 or x = xmax boundary is open, the bottom should be flat at the boundary. In the case of an open boundary condition, the water does not have to be shallow. Care should be taken to make sure that the onshore and offshore boundaries are located far from the main modal structure's location of maximum value.
The numerical grid extends one point outside the x, z domain at each end (nn is the total number of grid points, and mm the total number of grid points). There is also a grid point outside each boundary. Thus, for example, nn = 5 gives one interior grid point, 2 on the boundaries and 2 outside the domain.
In choosing a grid size, caution is advised. One should experiment with grid sizes in order to be confident of a good answer. I find that nn = 80-120, mm = 20-30 seems to work well in many cases. The number of operations or amount of memory needed goes something like nn*mm*mm. There is a test for hydrostatic consistency (the ratio computed is hxΔx/(hΔθ) ), and results seem best if the ratio given is less than around 1-5. Further, steep bottom topography can cause problems. It appears that a reasonable criterion is to pick nn such that The mean alongshore flow is assumed to be Gaussian (with different onshore, offshore, up and down length scales, as well as a variable center position and depth). This provides a fair bit of flexibility. Once the velocity field is computed, the density structure is computed using the thermal wind equation. The density field can be left undisturbed at either the offshore or onshore end of the domain. If you want to use some other functional form for mean velocity, you need only change lines 13-76 of bigrvzcal2.m. The derivatives of the mean flow and the density corrections are computed elsewhere numerically, so that these are the only lines that would need changing.
Results for pressure, velocity and density are plotted.
Although the long-wave limit is well defined when the mean alongshore flow is non-zero, the resulting wave modes are generally not orthogonal. This means that you cannot compute the first order wave equations coefficients (an, bn, cn) when there is a mean flow. Thus, you are not allowed to call the long-wave normalization and coefficient routine when a mean flow is present. If you make the long wave approximation. You will only be allowed to calculate one point on the dispersion curve.
All pressure fields are normalized with the long-wave norm (although this is only valid when there be no mean flow and waves are long). Thus, all other variables are scaled consistently with the long wave normalization (Brink, 1989), corrected for the possibility of a closed offshore boundary. Units are always consistent within c.g.s., i.e. cm/sec, dyne/cm 2 , gm/cm 3 .
The search can be very inefficient if you give it 0 as the starting guess for frequency. Better to give it a nonzero value. By the same token, if you give the algorithm an initial estimate for frequency that has few places of accuracy (say 1e-5), the Matlab search code will look around at fairly large increments (perhaps ± 20%), and the search can miss nearby solutions. If you provide more places of accuracy (say 1.01e-5), the search will be initially confined to more nearby locations (perhaps ± 3%), and it becomes much more likely that nearby solutions will be found.
The search can sometime miss a frequency. If it is a good solution, the inverse resonance parameter in the output (rrr) should be several orders of magnitude lower than neighboring values. If this is not true, you may either have a bad solution, or you may only be finding the real part of a complex wavenumber (which means you should use bigcomp.m to see if the flow is stable). Using a smaller accuracy tolerance (e) will tend to eliminate false solutions.
For an accuracy estimate ("e" option in bigrfinch.m), I find that 0.0001 works out fairly well. This means accuracy to 0.01% of the absolute value of the initial frequency guess.
The bottom stress is assumed to be proportional to the bottom velocity. For the errors associated with this assumption in the presence of a mean flow, see Brink (1997).
If more than one frequency is computed, there is an option to save the "dispersion curve" information (wavenumbers and frequencies) as a "mat" file.
In the absence of a mean flow, the perturbation bottom frictional decay time is computed for each frequency (Brink, 1990), allowing for the possibility of a closed offshore boundary and variable resistance coefficient. The e-folding time is only valid for weak bottom friction (small r).
Computing this, and the long wave coupling terms, are the only reasons for providing r for the program.
The format for the input files differs slightly from the general case mfiles (bigc*.m) and these mfiles for use only with real frequencies (bigr*.m). In order to make these input conversions conveniently, there are two translation functions, "c2rinpcon.m" and "r2cinpcon.m".
Given the presence of the near-inertial spurious mode, and given that the "open" boundary condition is only reasonable for subinertial frequencies, the model results should only be trusted for subinertial frequencies, |ω| <| f | (except in the long wave limit, where ω/l is what is actually important, and the idea of super-or sub-inertial does not make sense). If, in computing a dispersion curve, you start to obtain frequencies greater than the inertial, the superinertial results should be discarded.
All frequencies are in radians/sec and all wavenumbers are in radians/cm. Because radians are nondimensional, this could also be written as "1/sec" etc.
You can get spurious results if the domain is too large relative to the natural scale of the wave under consideration. For example, with a domain size xmax = 70 km, trying to find a flat-bottom internal Kelvin wave with internal Rossby radius = 3 km yields bad results. Doing the same calculation with xmax ≤ 20 km yields reasonable results. This problem arises because the arbitrary forcing used to create the resonance is applied near xmax/2. Trouble can arise when the desired wave has negligible amplitude at this distance offshore.
Once you include a sheared mean flow, it is possible to get a mean state that is gravitationally unstable (N 2 < 0). In this case, the unstable area will be displayed by blue pluses in the density contour plot, and you will not be allowed to proceed.
When the Ertel potential vorticity q = -N 2 f* changes sign relative to f, symmetric instability is possible (Hoskins, 1974). So, the routine "bigrrhocal2.m" tests to see if f*/f is negative anywhere.
If this happens, f* is contoured in red on the same plot as mean density, and a message appears on your Matlab screen. If f*/f does not become negative, it is not contoured, and no message about it is given. Experience with the program suggests that modes calculated when f* changes sign are of questionable merit. Structures are then often complex and have short scales near the region of negative f*. It seems probable that if, in nature, a flow is symmetrically unstable, it will quickly mix to a neutral state with zero or positive f*/f. So, even a "good" modal solution for a symmetrically unstable case is of questionable value.

Sample: Running bigrsetup.m
The following is a copy of what you see on the screen as you run bigrsetup.m.

>> datnew= bigrsetup;
This mfile will ask you a sequence of questions that are used to build the input file. Comment: For the depth h and bottom resistance coefficient r, if the last value you enter is for an x that is less than the domain size, the last value you enter will be used to fill out all values farther offshore. That is, the value is constant beyond the last one entered.

Comment:
-For r, if you enter 0 values to read, it will not ask for any more information on r. If you enter one or more values to read, it will prompt you for the x, r pairs. -Likewise, if you enter that the peak v0 = 0, then it will ask for no further information about the mean flow. If you enter that the peak v0 is nonzero, then you will be asked for the location (x, depth) of the maximum flow, and the horizontal (onshore, offshore) and vertical (up, down) Gaussian scales. -In setting up the dispersion curve in the input file, if you only want one point (one wavenumber: npts = 1), you will only be asked for that wavenumber. If you are setting up a curve with more than one point (npts > 1), you will also be asked for the wavenumber increment along that curve.

Mfiles in this package:
This collection of mfiles includes the following:

bigromp.m
The driver for the coastal trapped wave calculations. Carries out all input functions, calls the functions for computing mean density and velocity fields, and for searching for resonant frequencies. Also does some output and graphics functions.

bigrsetup.m
This mfile sets up an input array for "bigromp", based on user responses to questions. An example of its use is given above.

bigrfinch.m
This mfile changes an existing "bigromp" input array, without having to generate an all new array. An example of its use is given above.

c2rinpcon.m
The "bigcomp" (complex frequency) and "bigromp" (strictly real frequency) model systems call for slightly different input array structures. Use this mfile to convert a "bigcomp" input file into a "bigromp" file.

r2cinpcon.m
The "bigcomp" (complex frequency) and "bigromp" (strictly real frequency) model systems call for slightly different input array structures. Use this mfile to convert a "bigromp" input file into a "bigcomp" file.

bigrdep.m
Computes a gridded depth profile, along with its derivatives.

bigrconchk.m
Grid consistency check to see if the bottom slope is large relative to vertical spacing.

Bigrnorm.m
Normalizes pressure by the long wave norm for every calculation.

bigrvzcal2.m
Calculates the mean velocity field and its vertical derivative, from inputs.

bigrrhocal2.m
Calculates the modifications to the basic density field caused by the mean flow.

bigrdrver.m
Mainly called from fminsearch. This sets up the matrix equations that are equivalent to the field differential equation and its boundary conditions. Also solves the equation, and computes an integral of the response to evaluate convergence.

bigrsigtoz.m
Converts a field from sigma coordinates to z coordinates.

bigrengdiag.m
Does the energy diagnostics for each solution.

bigrztosig.m
Converts z coordinates to sigma coordinates.

bigrlong.m
Computes coupling coefficients and normalizations for the long-wave limit.

bigrtf.m
Computes perturbation estimate of frictional damping time.
The following are called by bigrdrver.m:

Bigrcc.m
Computes coefficients for the coastal boundary condition.

bigroff.m
Computes coefficients for the offshore boundary condition.

bigrbot.m
Computes coefficients for the bottom boundary condition.

bigrsbc.m
Computes coefficients for the surface boundary condition.

bigrmatco.m
Computes coefficients for interior points on the finite difference grid.

bigrint.m
Computes the integral of the response magnitude.
Note: Most of these files were originally written for use with the 'bigcomp" system, and then modified for use with strictly real frequency. When they were changed, the names were all changed to begin with "bigr".