QUAGMIRE v1.3: a quasi-geostrophic model for investigating rotating fluids experiments

QUAGMIRE is a quasi-geostrophic numerical model for performing fast, high-resolution simulations of multi-layer rotating annulus laboratory experiments on a desktop personal computer. The model uses a hybrid finite-difference/spectral approach to numerically integrate the coupled nonlinear partial differential equations of motion in cylindrical geometry in each layer. Version 1.3 implements the special case of two fluid layers of equal resting depths. The flow is forced either by a differentially rotating lid, or by relaxation to specified streamfunction or potential vorticity fields, or both. Dissipation is achieved through Ekman layer pumping and suction at the horizontal boundaries, including the internal interface. The effects of weak interfacial tension are included, as well as the linear topographic beta-effect and the quadratic centripetal beta-effect. Stochastic forcing may optionally be activated, to represent approximately the effects of random unresolved features. A leapfrog time stepping scheme is used, with a Robert filter. Flows simulated by the model agree well with those observed in the corresponding laboratory experiments.


Introduction
For over a century, geoscientists have invoked the principles of dynamical similarity (e.g. Douglas and Gasiorek, 2000) and geometrical similarity in order to study planetary atmospheres and oceans indirectly in the laboratory. For example, the mid-latitude atmospheric flow on a rotating planet closely resembles the flow in a rotating laboratory annulus, as suggested by Fig. 1. This statement holds despite typical length and time scales for corresponding atmospheric and labo-Correspondence to: P. D. Williams (p.d.williams@reading.ac.uk) ratory phenomena differing by many orders of magnitude. What matters is equality of the relevant non-dimensional parameters, such as the Rossby number (for dynamical similarity) and the horizontal-to-vertical aspect ratio (for geometrical similarity).
Once a particular fluid flow problem has been solved by making observations in the laboratory, an infinite number of other fluid flow problems have also effectively been solved, all of which are dynamically and geometrically similar. The main benefits of laboratory experiments are that they are under the complete control of the operator; that global highresolution measurements may be taken; and that controlled experiments may be repeated as many times as required. None of these statements holds when the atmosphere and oceans are studied directly.
A review of the role of laboratory experiments in geophysical fluid dynamics has been given by Hide (1977). Laboratory investigations of non-rotating fluids began in the nineteenth century, and include the classic experiments of Reynolds (1883). At around the same time, Vettin (1884) was probably the first to exploit dynamical similarity by carrying out rotating laboratory experiments as analogues of geophysical flows. He studied the surface flow in a rotating dishpan of fluid with a lump of ice near the centre, representing a polar ice cap, and he drew meteorological conclusions from his results (to the scorn of his contemporaries).
For the most direct resemblance between annulus and planet, heating and cooling should be applied at the outer and inner sidewalls, respectively, in order to mimic the equator-to-pole temperature gradient. However, it follows from thermal (and gradient) wind balance that a radial temperature gradient will be accompanied by a vertical shear in the zonal velocity. Therefore, analogous flows may be obtained in an isothermal annulus by directly imposing a shear using a differentially-rotating lid. The continuouslystratified thermally-forced rotating annulus and the two-layer Published by Copernicus Publications on behalf of the European Geosciences Union. From Read et al. (1998).
The development of computer models for simulating the general circulation of planetary atmospheres and oceans has been accompanied by the development of computer models for simulating rotating laboratory experiments. Because of their relative simplicity, laboratory flows are generally easier to model than atmospheric and oceanic flows. Therefore, the comparison between laboratory and model flows remains an important testbed for investigating many fundamental dynamical phenomena. In this paper, we describe the development of a numerical model for simulating fluid flows in the multi-layer rotating annulus. The model is named QUAG-MIRE, the QUAsi-Geostrophic Model for Investigating Rotating fluids Experiments.
Section 2 reviews a variety of possible modelling approaches, each making different dynamical and geometrical assumptions. A balanced model with a full representation of the cylindrical geometry is the preferred approach, for a number of important reasons which are discussed. The twolayer quasi-geostrophic coupled partial differential equations in cylindrical coordinates are derived in Section 3, and the diagnostic relations between streamfunction and potential vorticity are decomposed into vertical and azimuthal normal mode form in order to simplify their solution. Suitable sidewall boundary conditions are derived by considering integral properties of the governing equations. Then, the continuous equations are carefully discretized in Section 4, in such a way as to preserve discrete analogues of the integral properties. Suitable initial conditions and numerical parameter values are given. In Section 5, the calculations are partitioned into model subroutines, the technical details of how to run the model are described, and the code is tested. The paper concludes with a summary in Section 6.

Models of the rotating annulus
In this section, the relative merits of different possible dynamical (Section 2.1) and geometrical (Section 2.2) choices will be summarised.

Possible dynamical choices
For numerically modelling the laboratory annulus, one possible approach would be a direct numerical simulation (DNS) of the Navier-Stokes equations or shallow-water equations, which are both referred to as primitive equations because both vortical and divergent eigenmodes are retained. DNS codes have been developed for the continuously-stratified thermally-forced rotating annulus (e.g. Hignett et al., 1985;White, 1986) but they are computationally expensive and generally can be used to examine a small number of casestudy flows only.
As an alternative to the existing DNS models, in this paper we develop a balanced model in which the divergent eigenmodes are filtered out by construction. The filtering is justified because the interaction between vortical and divergent eigenmodes is thought to be weak. Balanced models have fewer dynamical degrees of freedom than primitive equation models, and therefore run much more quickly, allowing larger numbers of simulations to be performed.
Three possible balanced models for multi-layer flows are those based on the quasi-geostrophic equations, the balance equations and the slow equations. These three equation sets are each derived from the shallow-water equations, which in turn are derived from the Navier-Stokes equations by assum- Fig. 1. Diagram showing the analogy between (a) the fluid in a rotating annulus experiment in the laboratory, and (b) the mid-latitude atmosphere bounded by two latitude circles on a rotating planet. From Read et al. (1998). mechanically-forced rotating annulus have both been studied extensively (e.g. Hide et al., 1977;Carrigan, 1978;King, 1979;Appleby, 1982;Lovegrove, 1997;Williams, 2003).
The development of computer models for simulating the general circulation of planetary atmospheres and oceans has been accompanied by the development of computer models for simulating rotating laboratory experiments. Because of their relative simplicity, laboratory flows are generally easier to model than atmospheric and oceanic flows. Therefore, the comparison between laboratory and model flows remains an important testbed for investigating many fundamental dynamical phenomena. In this paper, we describe the development of a numerical model for simulating fluid flows in the multi-layer rotating annulus. The model is named QUAG-MIRE, the QUAsi-Geostrophic Model for Investigating Rotating fluids Experiments.
Section 2 reviews a variety of possible modelling approaches, each making different dynamical and geometrical assumptions. A balanced model with a full representation of the cylindrical geometry is the preferred approach, for a number of important reasons that are discussed. The twolayer quasi-geostrophic coupled partial differential equations in cylindrical coordinates are derived in Sect. 3, and the diagnostic relations between streamfunction and potential vorticity are decomposed into vertical and azimuthal normal mode form in order to simplify their solution. Suitable sidewall boundary conditions are derived by considering integral properties of the governing equations. Then, the continuous equations are carefully discretized in Sect. 4, in such a way as to preserve discrete analogues of the integral properties. Suitable initial conditions and numerical parameter values are given. In Sect. 5, the calculations are partitioned into model subroutines, the technical details of how to run the model are described, and the code is tested. The paper concludes with a summary in Sect. 6.

Models of the rotating annulus
In this section, the relative merits of different possible dynamical (Sect. 2.1) and geometrical (Sect. 2.2) choices will be summarised.

Possible dynamical choices
For numerically modelling the laboratory annulus, one possible approach would be a direct numerical simulation (DNS) of the Navier-Stokes equations or shallow-water equations, which are both referred to as primitive equations because both vortical and divergent eigenmodes are retained. DNS codes have been developed for the continuously-stratified thermally-forced rotating annulus (e.g. Hignett et al., 1985;White, 1986) but they are computationally expensive and generally can be used to examine a small number of casestudy flows only.
As an alternative to the existing DNS models, in this paper we develop a balanced model in which the divergent eigenmodes are filtered out by construction. The filtering is justified because the interaction between vortical and divergent eigenmodes is thought to be weak. Balanced models have fewer dynamical degrees of freedom than primitive equation models, and therefore run much more quickly, allowing larger numbers of simulations to be performed.
Three possible balanced models for multi-layer flows are those based on the quasi-geostrophic equations, the balance equations and the slow equations. These three equation sets are each derived from the shallow-water equations, which in turn are derived from the Navier-Stokes equations by assuming hydrostatic balance and columnar flow. Discussions of these and other filtered models are given by McWilliams and Gent (1980) and McIntyre and Norton (2000).
The quasi-geostrophic equations (Charney et al., 1950) are derived by assuming that the potential vorticity is advected only by the geostrophic component of the flow, and that interface perturbations are much smaller than the mean layer depths.
The balance equations (Charney, 1955) are derived by performing a horizontal velocity decomposition into vortical and divergent components, and then truncating with respect to the divergent component. The balance described is more complicated, but also more accurate, than geostrophic balance. Efficient procedures have been developed for integrating the balance equations (Daley, 1982). However, in their most general form the balance equations have spurious non-physical wave solutions with phase speeds much larger than those of inertia-gravity waves (Moura, 1976).
The slow equations (Lynch, 1989) are derived in a similar manner to the balance equations, except that the truncation is performed more systematically (based on normal mode initialization) and the spurious solutions vanish. Numerical integrations of the slow equations show excellent agreement with initialized numerical integrations of the shallow-water equations.
Of the above three possible dynamical choices, the quasigeostrophic (Q-G) equations are used for the QUAGMIRE model developed in this paper. This is because only one scalar function of horizontal position per layer (the streamfunction) is needed to uniquely define the state of the system in a Q-G model, whereas three per layer (the streamfunction, velocity potential and geopotential) are needed in a balance equations or slow equations model. With three times fewer dependent variables, the computational advantages gained from using a Q-G model arguably outweigh the disadvantages of its slightly lower formal accuracy. Brugge et al. (1987) have developed a numerical Q-G model for simulating multi-layer flows in a rectangular channel. However, their model is not particularly suitable for simulating the flow in the laboratory annulus, for the following reasons.

Possible geometrical choices
First, the channel equations with periodic azimuthal boundary conditions are a good approximation to the annulus equations only if the width of the annular gap is much smaller than the mean radius (King, 1979). With this condition satisfied, the curvature becomes negligible and it would be possible to justify the use of a channel model to simulate the flow in an annulus. For typical laboratory annulus experiments, however, the condition is not satisfied.
Second, channel models have additional shift-reflect symmetries not present in annulus models (Cattaneo and Hart, 1990). This is because, although the annulus and periodic channel are topologically similar, the geometry of their boundaries is fundamentally different. For example, there is a reflect symmetry in the channel, in the plane that is equidistant from the sidewalls, but there is no analogous symmetry in the annulus. Kwon and Mak (1988) show that the existence of such additional symmetries in the periodic channel leads to certain vortical wave-wave interaction coefficients being identically zero. Only models in cylindrical geometry admit the complete set of wave-wave interactions.
Third, the channel and the annulus both contain background potential vorticity gradients, due to the sloping of equilibrium geopotential height surfaces in the presence of a vertical shear in horizontal velocity. In the channel, these geopotential height and potential vorticity variations are linear in the across-channel direction, but in the annulus they are quadratic because of the parabolic equilibrium interface height shape produced by centripetal effects. This is known as the quadratic β-effect and is captured only by using cylindrical geometry.
Finally, the model of Brugge et al. (1987) does not include the effects of interfacial tension, which can be significant in the laboratory. Interfacial tension is scale-selective, affecting small scales much more strongly than large scales. Therefore, although interfacial tension effects are almost always negligible in real geophysical flows, they can be important in the analogous laboratory flows.
For the above reasons, the existing multi-layer Q-G models are not particularly suitable for simulating the flow in the laboratory annulus. Therefore, the remainder of this paper describes the construction of a new model that takes into account cylindrical geometry and interfacial tension.

Continuous equations
In this section, the governing continuous equations are derived (Sect. 3.1) and decomposed into normal mode form (Sect. 3.2), and suitable boundary conditions are derived (Sect. 3.3).

Derivation from first principles
The annulus to be modelled is shown schematically in Fig. 2. We use cylindrical polar coordinates, r=(r, θ, z). The z-axis is coincident with the rotation axis. The fluid is bounded by a base of mean vertical position z=0, a lid of mean vertical position z=2H >0 and cylindrical walls at r=a and r=b>a. The base and lid linearly deviate from their mean vertical positions by d bot (r) and d top (r). We define the constants s bot =dd bot /dr and s top =dd top /dr. The two homogeneous, immiscible layers have constant densities, ρ i , kinematic viscosities, ν i and mutual interfacial tension, S. We use the oceanographic convention that i=1 refers to the upper layer and i=2 to the lower layer. The undisturbed fluid interface is at z=H and the disturbed fluid interface is at z=H +η(r, t). The acceleration due to gravity is g. The base and sidewalls rotate about the axis of symmetry with angular velocity , and the lid rotates about the axis of symmetry with angular velocity + . In the frame that rotates with the base, the four fundamental equations for the pressure, p i (r, t), and the velocity, fluid interface is at z = H + η(r, t). The acceleration due to gravity is g. The base and sidewalls rotate about the axis of symmetry with angular velocity Ω, and the lid rotates about the axis of symmetry with angular velocity Ω + ∆Ω.
In the frame which rotates with the base, the four fundamental equations for the pressure, p i (r, t), and the velocity, u i (r, t), are the Navier-Stokes equations, and the equation of volume conservation for the incompressible liquids, Defining the vorticity by ω i = ∇ × u i , we take the curl of equation (1) and use vector identities to obtain the z-component of which, in the fluid interiors (i.e. away from the boundary layers) where the flow is assumed to be columnar and inviscid, is where ξ i is the z-component of ω i , f = 2Ω is the Coriolis parameter and u i, z is the vertical velocity.
We next vertically integrate equation (4) over the fluid interiors, parameterizing vertical Ekman pumping and suction velocities (Gill, 1982) at the lid, base and interface, which are all assumed to have small slopes. Assuming that the Ekman layer depths are much smaller than the total layer depths, and making the usual quasi-geostrophic assumptions (i.e. η ≪ H, d bot ≪ H, d top ≪ H and ξ i ≪ f ), after rearrangement we obtain and where and where q i (r, θ, t)/H are the perturbation potential vorticities (PPVs) given by and To complete the derivation, we write all of the dependent variables (i.e. u i , ξ i and η) in equations (5)-(8) in terms of the streamfunctions, ψ i (r, θ, t), defined by and The streamfunctions are defined here only to within arbitrary additive constants, which will be discussed in Section 3.3.2. The vorticities are given by Assuming hydrostatic balance and nearly equal densities, the interface height perturbation is given in terms of the streamfunctions (to within an additive constant) by where g ′ = 2g(ρ 2 −ρ 1 )/(ρ 2 +ρ 1 ) is the reduced gravity. The Laplacian term involving δ m = S/[g(ρ 2 − ρ 1 )] represents the effects of interfacial tension for an interface of small curvature. δ m is the characteristic static meniscus width, as can be seen by considering solutions to equation (12) when the tank is at rest (i.e. Ω = 0) and the fluid velocities are zero (i.e. ψ i = constant). When given the ψ i , equation (12) is u i (r, t), are the Navier-Stokes equations, and the equation of volume conservation for the incompressible liquids, Defining the vorticity by ω i =∇×u i , we take the curl of Eq.
(1) and use vector identities to obtain the z-component of which, in the fluid interiors (i.e. away from the boundary layers) where the flow is assumed to be columnar and inviscid, is where ξ i is the z-component of ω i , f =2 is the Coriolis parameter and u i, z is the vertical velocity. We next vertically integrate Eq. (4) over the fluid interiors, parameterizing vertical Ekman pumping and suction velocities (Gill, 1982) at the lid, base and interface, which are all assumed to have small slopes. Assuming that the Ekman layer depths are much smaller than the total layer depths, and making the usual quasi-geostrophic assumptions (i.e. η H , d bot H , d top H and ξ i f ), after rearrangement we obtain where and where q i (r, θ, t)/H are the perturbation potential vorticities (PPVs) given by and To complete the derivation, we write all of the dependent variables (i.e. u i , ξ i and η) in Eqs. (5)-(8) in terms of the streamfunctions, ψ i (r, θ, t), defined by and u i, r = − 1 r ∂ψ i ∂θ .
The streamfunctions are defined here only to within arbitrary additive constants, which will be discussed in Sect. 3.3.2. The vorticities are given by Assuming hydrostatic balance and nearly equal densities, the interface height perturbation is given in terms of the streamfunctions (to within an additive constant) by where g =2g(ρ 2 −ρ 1 )/(ρ 2 +ρ 1 ) is the reduced gravity. The Laplacian term involving δ m = √ S/[g(ρ 2 −ρ 1 )] represents the effects of interfacial tension for an interface of small curvature. δ m is the characteristic static meniscus width, as can be seen by considering solutions to Eq. (12) when the tank is at rest (i.e. =0) and the fluid velocities are zero (i.e. ψ i =constant). When given the ψ i , Eq. (12) is a forced Helmholtz equation for η, where the boundary conditions are the slopes, ∂η/∂r, at the sidewalls, which are related to the interface contact angle. We require an explicit formula for η, and so we seek a first-order solution to the Helmholtz equation for weak interfacial tension, by estimating the ∇ 2 η term in Eq. (12) using the solution for η when δ m = 0. This approach gives where 1 and δ 2 m ∇ 2 are the first two terms in a power series solution. On simple grounds, the series would be expected to converge rapidly if δ 2 m ∇ 2 η η, which is the case if δ 2 m λ 2 for waves of wavelength λ. We expect waves to form on the scale of the internal Rossby radius, g H /|f |, and so the convergence criterion becomes δ 2 m f 2 /g H 1. This is equivalent to F I 1 where F =f 2 (b−a) 2 /g H is the Froude number and I =δ 2 m /(b−a) 2 is the interfacial tension number (Appleby, 1982).
We finally substitute Eqs. (9)-(11) into (5) and (6) to obtain the two coupled partial differential equations governing the evolution of quasi-geostrophic motions in the two-layer annulus, and D Dt 2 The total derivative operators are given by and the horizontal Laplacian operator is given by By substituting Eqs. (11) and (13) into Eqs. (7) and (8) we obtain On the right side of Eq. (14), the term in square brackets represents spin-up/down by the frictional Ekman layers at the lid (∇ 2 ψ 1 ) and interface (∇ 2 (ψ 1 −ψ 2 )). The remaining term is the (constant) forcing term, and represents generation of potential vorticity by the rotating lid, communicated to the interior by the Ekman layer. The terms on the right side of Eq. (15) have a similar interpretation, except that there is no explicit forcing term in this case.
Equations (18) and (19) are similar to the relationships between potential vorticity and streamfunction in the channel model of Brugge et al. (1987), except that our equations include interfacial tension, and their βy term has been replaced with our β * r 2 term. The β * r 2 term is the quadratic β-effect. It is equal and opposite in the upper and lower layers, corresponding to the fact that depth increases in one layer are accompanied by equal decreases in the other layer.
Upon non-dimensionalization of Eqs. (14), (15), (18) and (19), using time scale ( ) −1 and horizontal length scale (b−a), definitions of Froude number, dissipation parameter, Rossby number, Reynolds number, Ekman number and interfacial tension number appear naturally. We choose to code QUAGMIRE using dimensional units, however, and therefore do not carry out the non-dimensionalization here. The advantage of using dimensional units is that the link to the laboratory annulus is immediate and clear. The disadvantage is that the non-dimensional parameters need to be computed separately.
We now list the assumptions that were required in order to derive Eqs. (14)-(19). It is important to bear these approximations in mind, since they limit the applicability of the model. We assume: The final two assumptions are not discussed until Sect. 3.3, but are included here for completeness. P. D. Williams et al.: QUAGMIRE v1.3 There is an equilibrium solution to Eqs. (14)-(19) of the form u i, r =0 and u i, θ =r i . Substituting allows us to determine the interior solid-body rotation rates, where χ = √ ν 2 /ν 1 . The corresponding interface height (to within an additive constant) is given by Eq. (13) to be Equations (20)- (22) describe the basic equilibrium state upon which baroclinically-unstable perturbations may grow. We refer to this state as the mean flow and we label the corresponding streamfunctions and PPVs ψ i (r) and q i (r), respectively.
Governing equations for perturbations, ψ i (r, θ, t) and q i (r, θ, t), to the mean flow are obtained by substitut- and D Dt 2 where and The total derivatives in Eqs. (23) and (24) advect according to the perturbation streamfunctions, i.e.
Equations (23)-(26) are the nonlinear model equations, discretized versions of which QUAGMIRE solves numerically. The constant forcing term in Eq. (14), which represents forcing of the full flow by the lid rotation, has been replaced in Eqs. (23) and (24) with more complicated non-constant forcing terms, which represent forcing of the perturbation flow by the equilibrium state and the bottom and top topography. An analytical assessment of the stability of small perturbations could begin by linearizing Eqs. (23) and (24), i.e. discarding the advection terms, but we retain all of the nonlinear terms in QUAGMIRE.
The perturbation velocity fields are given in terms of the perturbation streamfunctions by and which are the perturbation forms of Eqs. (9) and (10). The perturbation interface height field is given (to within an additive constant) by which is the perturbation form of Eq. (13).

Normal mode decomposition of the diagnostic equations
Given ψ i and q i at any time, we may evaluate ∂q i /∂t at that time using the prognostic equations, (23) and (24). There are contributions to the PPV tendency from advection (J (ψ i , q i )), forcing (∂/∂θ) and dissipation (∇ 2 ). We may use the tendency to determine q i at a short time in the future. We may then invert the diagnostic Helmholtz equations, (25) and (26), in order to obtain ψ i at the future time. Finally, we may begin the loop again using the updated fields. This is how QUAGMIRE integrates the model equations.
The Helmholtz equations, (25) and (26), are coupled. The inversion is made easier by first re-writing them in vertical normal mode form in order to remove the coupling. We take the sum and difference of the equations to obtain, respectively, and where C itcc is an interfacial tension correction coefficient given by We know that f 2 δ 2 m /g H 1 (Sect. 3.1), and so C itcc is slightly larger than unity, and is exactly equal to unity if the interfacial tension is zero.
Defining the barotropic (bt) and baroclinic (bc) vertical normal mode variables to be then Eqs. (31) and (32) both become uncoupled Helmholtz equations of the form where m=bt, bc. The eigenvalues are λ bt =0 and λ bc =2C itcc f 2 /g H .
We now perform a second normal mode decomposition, this time in the azimuthal dimension, in order to further simplify the solution of the Helmholtz equations. At each time step, we expand and Since m (r, θ ) and Q m (r, θ) are real, the complex functionŝ n m andQ n m satisfyˆ n m =(ˆ −n m ) * andQ n m =(Q −n m ) * , where the asterisk represents complex conjugation. The n=0 term is the mean flow correction, which is a correction to the mean flow that is generated by nonlinear self interactions of the waves. It is equal to the zonal average of the perturbation quantities, as can be seen from the zonal integration of Eqs. (39) and (40). The n =0 terms represent wave (or eddy) components. Substituting Eqs. (39) and (40) into (38) gives the radial structure equation, This complex ordinary differential equation must be solved for each combination of vertical modes, m∈{bt, bc}, and azimuthal modes, n∈{0, ±1, ±2, . . .}, in order to determineˆ n m (r) when givenQ n m (r). The inversion process required in order to obtain ψ i (r, θ ) from q i (r, θ ) is therefore summarized (with the relevant equation numbers in brackets) as: We could finally perform a third normal mode decomposition, this time in the radial dimension, by projectingˆ n m (r) andQ n m (r) onto the eigenfunctions of the linear operator on the left side of Eq. (41). The baroclinic eigenfunctions are modified Bessel functions of order n in the scaled radial coordinate,r= √ λ bc r (Boas, 1983) and the barotropic eigenfunctions are of the form r ±n . However, this approach would force the streamfunction and PPV to satisfy the same boundary conditions, for which there is no justification. QUAG-MIRE therefore solves the radial structure equation directly rather than projecting onto the radial modes.

Perturbation streamfunction boundary conditions for the continuous equations
We must now choose boundary conditions to apply to the perturbation streamfunction when integrating Eq. (41). The equation was derived under the assumption of inviscid flow. Therefore, it cannot describe the thin, viscous Stewartson layers of width δ S that exist at the lateral boundaries, and applies only to the fluid interior, a+δ S <r<b−δ S . We assume that δ S a and δ S b so that we may still write the integration range as a≤r≤b, but now when we refer to r=a or r=b we mean the boundary between the fluid interior and Stewartson layer rather than the sidewall. There are a number of possible boundary conditions. To impose passive Stewartson layers that do not anywhere exchange fluid with the interior, we would apply the impermeability condition to the radial perturbation velocity, i.e. u i, r | r=a, b =0, ∀ θ, i, which in the normal mode variables corresponds to Dirichlet boundary conditions, The mean flow correction velocity (n=0) is purely zonal and so automatically satisfies impermeability. Impermeability alone is therefore not a sufficient condition to uniquely specify a solution. No-slip boundary conditions for the zonal perturbation velocity, i.e. u i, θ | r=a, b =0, ∀ θ, i, correspond in the normal mode variables to the Neumann conditions, The equilibrium solid-body rotation flow about which we perturb satisfies impermeability but does not satisfy the noslip condition.
Since we are solving a second-order differential equation, only two independent boundary conditions are required. We cannot therefore impose both impermeability and no-slip flow at both boundaries, as that would require four independent conditions. The over-constrained nature of the PPV inversion in Q-G models has been discussed by Williams (1979). A comprehensive study of the comparative effects of using no-slip boundary conditions, rather than the more traditional free-slip conditions, is described by Mundt et al. (1995).
We must use a reduced set of boundary conditions, but we must choose carefully and consistently which conditions to retain and which to abandon, in order to avoid non-physical behaviour. We are, of course, free to employ different boundary conditions for the different normal mode components specified by m and n.
The debate over suitable lateral Q-G boundary conditions has had a long and contentious history in the literature. In the classic periodic channel models of Phillips (1954Phillips ( , 1956, Eq. (42) is applied to the wave (n =0) components and Eq. (43) is applied to the mean flow correction (n=0) component. The latter condition is not imposed (but the former is retained) by Phillips (1963) and Pedlosky (1964), but McIntyre (1967) shows that relaxing the mean flow correction boundary condition leads to a spurious, unspecified energy flux through the sidewalls. The condition is included again by Pedlosky (1970), but replaced in Pedlosky (1971) and Pedlosky (1972) with an ad-hoc condition chosen for mathematical convenience. Smith (1974) points out that the resulting non-physical energy source might invalidate Pedlosky's results, and repeats Pedlosky's calculations with the proper boundary condition retained (Smith and Pedlosky, 1975;Smith, 1977). More recent studies (Appleby, 1982;Yoshida and Hart, 1986;Lewis, 1992;Stephen, 1998) have avoided the spurious energy flux by applying both conditions in full, as done by Phillips (1954Phillips ( , 1956).
An informative interpretation of Phillips' mean flow correction boundary condition has been given by Davey (1978). For non-zero zonal perturbation velocities, u i, θ | r=a, b , at the boundary between the interior and a Stewartson layer, there will be a corresponding return volume flux between the Ekman layers and the Stewartson layer due to the asymmetry of the Ekman spiral (Pedlosky, 1987), which will have a nonzero radial component proportional to u i, θ | r=a, b . We can therefore ensure that there is no net build-up of mass in the Stewartson layers by setting This condition is automatically satisfied for the wave (n =0) components, and is equivalent to Eq. (43) with n=0, which is the condition used by Phillips. With this condition, there is no net exchange of fluid due to the perturbation flow between each Ekman layer and the Stewartson layers, although local exchange is allowed. Next, we derive a consistent and plausible set of boundary conditions for the quasi-geostrophic annulus, which do not lead to non-physical behaviour, by considering integral properties of the prognostic (Sect. 3.3.1) and diagnostic (Sect. 3.3.2) model equations.

Integral properties of the prognostic equations
Consider the area-integral of the perturbation PPV tendencies over the annular domain, as given by the prognostic equations, (23) and (24). The forcing (∂/∂θ) terms integrate to give zero unconditionally. The advection (J (ψ i , q i )) terms integrate to give zero (Salmon and Talley, 1989) if The dissipation (∇ 2 ) terms integrate to give zero if The conditions (46) and (47) are equivalent to impermeability for the waves and no-slip for the mean flow correction, as originally used by Phillips (1954Phillips ( , 1956. With these conditions, the horizontal-mean PPV in each layer is conserved by the continuous equations and there is no spurious energy flux. We choose to apply these conditions in QUAGMIRE, except that the latter condition leads to an ill-posed PPV inversion for the special case of n=0 and m=bt, as we will now see.

Integral properties of the diagnostic equations
Equation (41) for the barotropic mean flow correction is Since λ bt =0 and n=0 for this case, one of the terms in the radial structure equation has vanished, making the left side an exact derivative. Equation (48) can therefore be integrated analytically between r=a and r=b to give In QUAGMIRE, we choose initial conditions for which the horizontal-mean barotropic PPV is zero, and it is then guaranteed to remain zero at all times, as shown in Sect. 3.3.1. This means that we need only explicitly set dˆ 0 bt dr r=a = 0 (50) from Eq. (49). If we explicitly impose both (50) and (51) when solving (48), we will have an underconstrained and illposed problem. We need an additional constraint to close the solution.
We have defined two streamfunctions in the model -one per layer, or equivalently, one per vertical normal modeand each of these has an integration constant associated with it (Sect. 3.1). Just because these two arbitrary constants have no physical meaning does not mean that they do not need to be defined in the numerical model. Now that we know that Eqs. (50) and (51) are not independent boundary conditions, and therefore that to explicitly impose both would lead to an underconstrained PPV inversion, we choose to explicitly impose only Eq. (50). We then take the opportunity to use the remaining degree of freedom associated with the solution of Eq. (48) to define one of the streamfunction integration constants, by arbitrarily settinĝ which completes the set of two boundary conditions for the m=bt, n=0 component, and gives a well-posed problem. Incidentally, the second streamfunction integration constant is defined by requiring the mean interface perturbation to be zero using Eq. (13), which follows from volume conservation for either layer. This requirement is imposed off-line, boundary conditions, and therefore that to explicitly impose both would lead to an underconstrained PPV inversion, we choose to explicitly impose only equation (50). We then take the opportunity to use the remaining degree of freedom associated with the solution of equation (48) to define one of the streamfunction integration constants, by arbitrarily settinĝ which completes the set of two boundary conditions for the m = bt, n = 0 component, and gives a well-posed problem. Incidentally, the second streamfunction integration constant is defined by requiring the mean interface perturbation to be zero using equation (13), which follows from volume conservation for either layer. This requirement is imposed off-line, by adding a suitably-chosen constant to one of the streamfunction fields when model diagnostics are plotted, and not as a boundary condition during the inversion.
A summary of the boundary conditions which we must explicitly set when integrating equation (41) is given in Table 1. With these conditions, the sidewall boundaries are impermeable to each component of the flow, i.e. the solid-body rotation equilibrium flow, the mean flow correction and the wave components. The boundaries are slippery to the solid-body rotation flow and the wave components, but no-slip to the mean flow correction.

Discretized equations
We derived in Section 3 a set of partial differential equations and boundary conditions which are both physically sensible and well-posed. We now discretize the equations so that they are suitable for numerical integration on a computer. We must take great care to ensure that the discretized equations and boundary conditions retain the important properties possessed by the continuous equations. In particular, it is important that they satisfy discretized analogues of the integral properties discussed in Section 3.3. by adding a suitably-chosen constant to one of the streamfunction fields when model diagnostics are plotted, and not as a boundary condition during the inversion. A summary of the boundary conditions that we must explicitly set when integrating Eq. (41) is given in Table 1. With these conditions, the sidewall boundaries are impermeable to each component of the flow, i.e. the solid-body rotation equilibrium flow, the mean flow correction and the wave components. The boundaries are slippery to the solid-body rotation flow and the wave components, but no-slip to the mean flow correction.

Discretized equations
We derived in Sect. 3 a set of partial differential equations and boundary conditions that are both physically sensible and well-posed. We now discretize the equations so that they are suitable for numerical integration on a computer. We must take great care to ensure that the discretized equations and boundary conditions retain the important properties possessed by the continuous equations. In particular, it is important that they satisfy discretized analogues of the integral properties discussed in Sect. 3.3.

The numerical grid
The regular grid on which we discretize the equations is shown in Fig. 3. The grid consists of N rad points in the radial dimension (including one point on each boundary, r=a 22 P. D. Williams et al.: QUAGMIRE v1.3 and r=b) and N azim points in the azimuthal dimension. We define and then we have and θ(j ) = j θ , j = 1, 2, . . . , N azim .
The point (i, N azim +1) is identical to the point (i, 1). We define the perturbation streamfunction ψ (i, j, k) and PPV q (i, j, k) at each of these points in each layer (k=1, 2), so that ψ and q are co-located on the grid. The area of the grid box with coordinates (i, j ) is approximately where δ * , * is the Kronecker delta function.

Prognostic equations
In the continuous case, we chose perturbation streamfunction boundary conditions such that each of the three contributions (advection, forcing and dissipation) to the areaintegrated perturbation PPV tendency was zero. We would now like to choose discretizations of these contributions, together with discretizations of the boundary conditions, for which this statement still holds exactly. If our discretization only conserves the mean PPV approximately, then there is the possibility of a non-physical and explosive increase in the PPV, even if the error is small, due to the compound effects of many time steps. Following Sect. 3.3.1, we therefore next examine the discretizations and boundary conditions necessary to ensure that for k=1, 2, where f (i, j, k) is, in turn, the discretized azimuthal derivative, Jacobian and Laplacian. The centred, second-order discretization of the azimuthal derivative, satisfies Eq. (57) unconditionally, as in the continuous case. The second-order Arakawa (1966) discretization of the Jacobian satisfies Eq. (57) if which is a discretized version of the condition, (46), for the continuous case. It is tedious but straightforward to show that the five-point discretization of the Laplacian (whose continuous definition is given in Eq. (17) for reference), with ghost point values, ψ (0, j, k) and ψ (N rad +1, j, k), given by linear extrapolation, and and which are discretized versions of the condition, (47), for the continuous case. There will be a small error in the value of the calculated discretized Laplacian at the boundaries due to the assumption of linearly-extrapolated ghost points, but there is apparently no other way to discretize the Laplacian in such a way that analogues of its integral properties are fully preserved. The assumption of linearly-extrapolated ghost points appears to be inconsequential, because we find good agreement between model and laboratory flows (see Sect. 6).

Diagnostic equations
The discretized versions of Eqs.
The summations have been truncated, compared to Eqs. (39) and (40), because there are only N azim independent Fourier components associated with the discrete Fourier transform of a series of N azim numbers. Because m (i, j ) is real, we havê for n=1, 2, . . . , N azim −1.
In terms of the normal mode variables, the discretized boundary conditions, (59), (63) and (64) We now consider the discretization of the radial structure equation, (41). Using centred three-point finite differences at the interior points, i=2, 3, . . ., N rad −1, we obtain n m (i − 1) − 2ˆ n m (i) +ˆ n m (i + 1) Re-grouping terms according to grid points gives where the dimensionless quantities α ± and γ are given by and In Cartesian geometry we would have α ± (i)=1.
The N rad −2 equations, (71), together with 2 boundary conditions, complete the set of N rad equations in the N rad unknowns,ˆ n m (i) for i=1, 2, . . ., N rad . These linear equations may be written in matrix form, where the zero elements in the tridiagonal N rad ×N rad matrix have been left blank. The two elements labelled "bdy" are boundary condition elements, dependent upon m and n.
There are two further such elements in the right-most two columns of the bottom row.

Perturbation streamfunction boundary conditions for the discretized equations
In the continuous case, we found that the boundary conditions for the barotropic mean flow correction component (m=bt, n=0) were ill-posed as originally stated, and remained so until we replaced a redundant boundary condition with an equation to define an integration constant. This happens in the discretized case, too: the square matrix in Eq. (74) is singular for the barotropic mean flow correction, when the boundary condition elements (labelled "bdy") are (−1, 1) in the top row and (1, −1) in the bottom row. The analytical proof of this, which involves showing that a certain linear combination of rows is zero, is tedious but straightforward. By analogy with the continuous case, we replace the two boundary condition elements in the bottom row with (0, 1) to define the integration constant by setting the streamfunction for this component to zero on the outer boundary, and then the matrix is no longer singular. In the continuous system, we set the n=0, m=bt normal streamfunction derivative to zero at one boundary and found that, if the mean barotropic PPV was zero, the streamfunction derivative would automatically be zero at the other boundary. Importantly, in contrast with the continuous system, this statement does not hold exactly for the discretized system. This is becauseQ n m (1) andQ n m (N rad ) do not appear in Eq. (74): we do not apply the discretized differential equation at the boundaries, because we need to use these two degrees of freedom to impose the boundary conditions. Table 2. Summary of the boundary conditions applied to the streamfunction when integrating the discretized equations. The analogous conditions for the continuous case are given in Table 1. † After the inversion,ˆ 0 bt (N rad ) is redefined byˆ 0 bt (N rad ) − 0 bt (N rad − 1)=0, as discussed in the text.
The error corresponding to this PPV leak is small, but even small errors can grow to dominate the solution after a large number of time steps. To fix this problem with the barotropic mean flow correction, we discard the outer boundary streamfunction,ˆ 0 bt (N rad ), obtained through inversion of Eq. (74) and define a new value for it by set-tingˆ 0 bt (N rad )=ˆ 0 bt (N rad −1). This ensures that the boundary conditions, (69), required for conservation of horizontalmean PPV are satisfied exactly. The consequence is that the discretized differential equation, (70), is not exactly satisfied at the point N rad −1. The imposed boundary conditions are summarized in Table 2.

Relaxation
Instead of (or in addition to) the mechanical forcing imposed by the differentially rotating lid, QUAGMIRE includes the option of relaxing the flow towards specified perturbation streamfunction or perturbation potential vorticity fields. If relaxation to a specified perturbation streamfunction is activated, then the perturbation streamfunction minus the relaxation perturbation streamfunction is used in the computation of the diffusion and hyperdiffusion terms. If relaxation to a specified perturbation potential vorticity is activated, then the perturbation potential vorticity is relaxed towards the relaxation perturbation potential vorticity at a prescribed rate.

Numerical methods
We now discuss the time-stepping scheme (Sect. 4.6.1), the need for time-lagged diffusion (Sect. 4.6.2), and the approximate representation of unresolved features using hyperdiffusion (Sect. 4.6.3) and stochastic forcing (Sect. 4.6.4).

Time stepping
For the time stepping we use a leapfrog scheme with a Robert (1966) three-level time filter applied at each time step, to suppress the computational mode splitting between even and odd steps (Mesinger and Arakawa, 1976). At each step, of size t, q t+1 is determined at each grid point using the leapfrog scheme, and then the value of q t is adjusted in such a way as to move it closer to the mean of q t−1 and q t+1 , according to The old value of q t is abandoned and the new, filtered value is used in its place. The Robert filter parameter, R>0, is chosen to be as small as possible whilst still suppressing the leapfrog decoupling.

Time-lagged diffusion
Numerical solutions of the simple diffusion equation, using the leapfrog scheme for the time discretization and a time-centred, three-point finite difference for the space discretization, are unconditionally unstable due to a computational mode (Haltiner and Williams, 1980). To avoid this in QUAGMIRE, we time-lag the diffusion terms by one time step when evaluating the right sides of the discretized analogues of Eqs. (23) and (24). This means that, when evaluating the PPV tendency, we calculate the forcing (∂/∂θ) and advection (J (ψ i , q i )) terms using the fields at the current time step, but we calculate the diffusion (∇ 2 ) terms using the fields at the previous time step.

Hyperdiffusion
To represent sub-gridscale effects we add a hyperdiffusion term to the right sides of the prognostic equations, (23) and (24), as is usual in numerical models (e.g. Lewis, 1992). At first, a fourth-order streamfunction hyperdiffusion term, ν hyper ∇ 4 ψ i , was tested, but significant gridscale features were always found to form at the lateral boundaries whenever the model was run. This is because during the PPV inversion, any gridscale features in the PPV field will give rise to corresponding gridscale features in the perturbation streamfunction field, and then the ν hyper ∇ 4 ψ contribution to the PPV tendency will tend to damp out these features in the PPV field. Unfortunately this does not happen at the boundaries in the discretized system, because boundary values of the PPV are not used when performing the inversion: as already discussed,Q n m (1) andQ n m (N rad ) are missing from Eq. (74). Values of PPV are able to feed back into the PPV tendency field only at interior points, and there is nothing to suppress gridscale features in the PPV field at the boundaries.
To avoid this, we instead use second-order hyperdiffusion applied to the PPV, by adding a term ν hyper ∇ 2 q i to the prognostic equations. This term is also time-lagged by one time step, as discussed above. The hyperdiffusion term does not exactly satisfy Eq. (57), but the error is very small: the contribution of the hyperdiffusion term to the horizontal-mean PPV tendency is typically 10 −10 s −2 . In order to keep the model solutions as close as possible to the exact solutions, we periodically reset the horizontal-mean PPV to zero in QUAGMIRE, by adding a very small constant whose value is chosen to satisfy this requirement.

Stochastic parameterization of sub-gridscale effects
As an alternative method for representing sub-gridscale effects, there is the option of a simple stochastic parameterization in QUAGMIRE. Such parameterizations are increasingly used in numerical models, with the recognition that the additional degrees of freedom they introduce may be able to compensate, at least partially, for the degrees of freedom missing from filtered models. Stochastic forcing is necessary to capture regime transitions between different baroclinic wave modes (nominally due to sub-gridscale inertiagravity waves), which are apparently not captured by the hyperdiffusion approach (e.g. Williams et al., 2004).
We choose the simplest possible form for the noise terms. At each grid point and at each time step, a random number is drawn from the uniform distribution on the interval [0, 1], and then shifted to the interval [−amp, amp] before being used as an additive contribution to the PPV tendency. At each grid point and at each time step, the added random number is equal and opposite in the upper and lower layers, nominally representing pure-baroclinic inertia-gravity waves. The constant, amp, is a given amplitude with units s −2 , which may change linearly with time in the model. The noise contains no correlations in either time or horizontal position. The horizontal-mean random number field is enforced to be zero in both layers.

Initial conditions
A feature of the leapfrog time-stepping scheme is that initial condition fields are required at two consecutive times, in order to begin the integration. We choose to specify the PPV fields as initial conditions. We use small amplitude random noise for these fields, seeding the system to permit the growth of unstable perturbations of any azimuthal and radial wavenumber. We generate random numbers from a uniform distribution, which we then shift to a chosen symmetrical interval centred on zero. We then subtract the mean PPV in each layer at both time steps, which makes the fields satisfy the zero-horizontal-mean barotropic PPV condition of Sect. 4.4.

Summary of numerical scheme
Flow charts summarizing the details of the QUAGMIRE numerical integration scheme are shown in Figs. 4 and 5. Given the PPV fields at times t−1 and t, we invert to obtain the streamfunction fields at those times, which then allows us to calculate all the contributions to the PPV tendency. We perform a leapfrog time integration to obtain the PPV field at time t+1, and then modify the PPV field at time t by applying a Robert filter. Once we have obtained q (t) and q (t+1) from q (t−1) and q (t), we discard q (t−1) and ψ (t−1), we write q (t) and ψ (t) to disk if required, then we re-label t→t−1 and begin the loop again. Note that the streamfunction and PPV must be kept in memory at three consecutive time steps.
The system state is completely determined by ψ . Note that the system state is also completely determined by q together with the boundary conditions, because Eqs. (25) and (26) are uniquely invertible. It is not necessary to write both ψ and q to disk in order to have a complete description of the system, therefore. Nevertheless, we choose to save both fields, in order to reduce the need for further calculations when plotting model diagnostics.

Suitable values for the numerical parameters
QUAGMIRE employs Fast Fourier Transform (FFT) routines, which are much faster if the only prime factors of N azim are 2, 3 and 5. A typical grid might be defined by N azim =2 5 ×3=96 and N rad =2 4 =16, as shown in Fig. 6. A suitable Robert filter parameter is usually around R=0.01. For given and , we recommend taking the amplitude of the random initial PPV perturbation to be /100, so that the growth of very small perturbations is assessed; taking the time step, t, to be such that the bulk azimuthal Courant number, 1 2 t/ θ, is 0.01; and taking the hyperdiffusion coefficient, ν hyper , to be such that the e-folding time, 1/(ν hyper k 2 Nyquist ), for damping of mid-radius gridscale waves with the Nyquist wave vector, k Nyquist =N azim /(a+b), is equal to one lid rotation period, 2π/ . By default, double numerical precision (retaining 16 significant figures) is used for the calculations and the pick-up dumps to disk, and single numerical precision (retaining 8 significant figures) is used for the regular dumps to disk. The factor by which relative errors in the perturbation streamfunction are greater than relative errors in the PPV, following solution of Eq. (74), is known as the condition number of the tridiagonal matrix in that equation. Some typical condition numbers for the matrices in Eq. (74) are shown in Table 3. The largest condition number in the system has a value of a few hundred, implying that only the last two significant figures of the inferred perturbation streamfunctions will be uncertain, and that errors due to rounding are therefore small.
In order to demonstrate insensitivity to the numerical parameters, comparative runs have been done with (separately) 26 P. D. Williams et al.: QUAGMIRE v1.3 14 P. D. Williams et al.: QUAGMIRE v1.3 MAIN LOOP

Read run parameters and initial fields (i) compute initial fields internally, or
INPUT: q ′ (t + 1) OUTPUT: ψ ′ (t + 1)  the hyperdiffusion coefficient decreased by a factor of 10, the Robert filter parameter decreased by a factor of 10, and the gridspacing doubled in both directions, but all other parameters unmodified (Williams, 2003). The equilibrated wave number was the same in each case, and the mid-radius wave amplitude and phase speed differed by at most 0.3%. We have therefore demonstrated that both rounding errors and discretization errors are small, and that the equilibrated state is insensitive to the values of the numerical parameters, im-plying that the model output gives an accurate representation of the true solutions of the continuous model equations.
In order to demonstrate insensitivity to the numerical parameters, comparative runs have been done with (separately) the hyperdiffusion coefficient decreased by a factor of 10, the Robert filter parameter decreased by a factor of 10, and the gridspacing doubled in both directions, but all other parameters unmodified (Williams, 2003). The equilibrated wave number was the same in each case, and the mid-radius wave amplitude and phase speed differed by at most 0.3%. We have therefore demonstrated that both rounding errors and discretization errors are small, and that the equilibrated state is insensitive to the values of the numerical parameters, implying that the model output gives an accurate representation

Source code
The source code is written in Fortran 95. Routines from the Numerical Algorithms Group (NAG) library are employed: nag fft for the transformations between physical and spectral space described by Eqs. (65) and (66); nag gen bnd lin sys for solving the complex band matrix equation, (74), N azim +2 times each time step; and nag math constants for the value of π.
The source code consists of 15 .f90 subroutines in the src/ directory. In total, there are 1200 lines of code in these subroutines, many of which are comments. Brief descriptions of the subroutines are now given: modules.f90 declares the global variables, categorized into five modules: precision (the numerical precisions for the calculations and dumps to disk); dyn vars (the dynamical state arrays, i.e. streamfunction and PPV, which are updated once per time step); solver vars (the permanent solver arrays, calculated once at the start of the model run); phys params (the physical parameters, including system dimensions and rotation rates); and grid params (the numerical parameters, including grid spacings and time steps).
main.f90 is the highest-level routine in the model, making one-off calls to init model.f90, init solver.f90 and init state.f90, and then calling in turn, from within the time stepping loop, jacobian.f90, forcing.f90,  of the true solutions of the continuous model equations.

Technical details
The model is available as a zip file, which contains the source code and makefile (Section 5.1), the namelist (Section 5.2), a shell script (Section 5.3) and a comprehensive Matlab diagnostics suite (Section 5.4). The zip file can be downloaded from PUBLISHERS INSERT URL HERE.

Source code
The source code is written in Fortran 95. Routines from the Numerical Algorithms Group (NAG) library are employed: nag fft for the transformations between physical and spectral space described by equations (65) and (66); nag gen bnd lin sys for solving the complex band matrix equation, (74), N azim + 2 times each time step; and nag math constants for the value of π.
The source code consists of 15 .f90 subroutines in the src/ directory. In total, there are 1200 lines of code in these subroutines, many of which are comments. Brief descriptions of the subroutines are now given: modules.f90 declares the global variables, categorized into five modules: precision (the numerical precisions for the calculations and dumps to disk); dyn vars (the dynamical state arrays, i.e. streamfunction and PPV, which are updated once per time step); solver vars (the permanent solver arrays, calculated once at the start of the model run); phys params (the physical parameters, including system dimensions and rotation rates); and grid params (the numerical parameters, including grid spacings and time steps).
main.f90 is the highest-level routine in the model, making one-off calls to init model.f90, init solver.f90 and init state.f90, and then calling in turn, from within the time stepping loop, jacobian.f90, forcing.f90, dissipation.f90, step q.f90, solver.f90 and save fields.f90.
init model.f90 initializes the model by reading the namelist (Section 5.2), and then allocating and evaluating the parameters declared in the phys params and grid params modules.
init solver.f90 initializes the solver by allocating and evaluating the arrays declared in the solver vars module.
init state.f90 initializes the model state by allocating and evaluating starting values for the arrays declared in the dyn vars module. This subroutine makes calls to read pu fields.f90 and read forcing fields.f90.
jacobian.f90 calculates the advection (J(ψ ′ i , q ′ i )) term in cylindrical polar co-ordinates, storing the result as a PPV tendency. The formula is written to minimise the number of multiplications, which are computationally more expensive than additions.
forcing.f90 calculates the forcing (∂/∂θ) terms (forcing of perturbations by the mean flow, topographic forcing, stochastic forcing and optional relaxation to a  dissipation.f90, step q.f90, solver.f90 and save fields.f90. init model.f90 initializes the model by reading the namelist (Sect. 5.2), and then allocating and evaluating the parameters declared in the phys params and grid params modules.
init solver.f90 initializes the solver by allocating and evaluating the arrays declared in the solver vars module.
init state.f90 initializes the model state by allocating and evaluating starting values for the arrays declared in the dyn vars module. This subroutine makes calls to read pu fields.f90 and read forcing fields.f90.
jacobian.f90 calculates the advection (J (ψ i , q i )) term in cylindrical polar co-ordinates, storing the result as a PPV tendency. The formula is written to minimise the number of multiplications, which are computationally more expensive than additions.
forcing.f90 calculates the forcing (∂/∂θ) terms (forcing of perturbations by the mean flow, topographic forcing, stochastic forcing and optional relaxation to a specified PPV field), adding the result to the PPV tendency. This subroutine makes calls to dtheta.f90.
dissipation.f90 calculates the dissipation (∇ 2 ) terms (Ekman layers at the upper and lower boundaries, Ekman layers at the internal interface and second-order PPV hyperdiffusion), adding the result to the PPV tendency. This subroutine makes calls to laplacian.f90.
step q.f90 uses the PPV tendency to perform the time stepping, with a Robert three-level time filter, and thereby updates the PPV field.
solver.f90 solves the Helmholtz equation in cylindrical polar co-ordinates, to update the streamfunction field given the updated PPV field.
save fields.f90 dumps the model state (streamfunction and PPV) to disk. The current state (regular dump) and/or the current and previous states (pickup dump) are saved as unformatted binary data.
read pu fields.f90 reads the initial model state from a pickup file, if required.
read forcing fields.f90 reads the forcing fields (PSI relax.bin and QGPV relax.bin) from disk, if required.
dtheta.f90 is the azimuthal derivative operator for calculating the forcing terms in the PPV equations.
laplacian.f90 is the Laplacian operator in cylindrical polar co-ordinates for calculating the dissipation terms in the PPV equations.
A Makefile is also included in the src/ directory, to build the executable from the source code files. The Makefile contains four environment variables (FC, FFLAGS, LIBS, and LM LICENSE FILE), which specify directory paths, flags and other information required by the Fortran compiler and NAG library. This information is site-specific and the default settings, in the supplied Makefile, are those appropriate to the AOPP machines at Oxford University. The default settings must be changed according to the local installations of the compiler and library, before the model is built at the user's site. Then, by typing make at a command prompt, each .f90 source file is compiled to produce a corresponding .o object file, and then the 15 object files (plus object files from the NAG library) are linked to build the executable, with filename qgam.
The azimuthal derivative (dtheta.f90), Laplacian (laplacian.f90) and advection (jacobian.f90) routines were each tested using input fields consisting of random numbers satisfying the boundary conditions. The mean PPV tendency due to each contribution was found to be zero to within numerical precision, which is a necessary (but not sufficient) condition for the code for these routines to be free from errors.
The Helmholtz solver routines (init solver.f90 and solver.f90) were tested by first using the forward formulae (25) and (26) with our discretized Laplacian (60)-(62) to calculate the PPV fields corresponding to given random perturbation streamfunction fields, and then using the routines to reconstruct the streamfunction fields from the calculated PPVs. The root-mean-square difference between the original and reconstructed streamfunction fields was around 0.1%, which is consistent with the solver code also being free from errors. The reason that the agreement is not exact, to within numerical precision, is that we assume linearly-extrapolated ghost points to evaluate the Laplacian in the forward formulae, an assumption that is not made during the inversion.

Namelist
The namelist, qgam.data, should be copied to the working directory (i.e. the directory in which the model output is to appear) and then edited to alter the physical and numerical details for the run. The entries in the namelist, which fall into five categories, are as follows. initial amplitude (REAL) = amplitude of initial PPV perturbation (s −1 ) -only used if start step is set to 1 omega (REAL) = = angular velocity of base (rad s −1 ) -must be positive lid delta omega (REAL) = = differential angular velocity of lid relative to base (rad s −1 )can be either positive (prograde) or negative (retrograde) g (REAL) = g = acceleration due to gravity (m s −2 ) new hyper (REAL) = ν hyper = second-order PPV hyperdiffusion coefficient (m 2 s −1 ) relax rate (REAL) = inverse time scale for relaxation to specified PPV field (s −1 ) -only used if relax type is set to 2 or 3 30 P. D. Williams et al.: QUAGMIRE v1.3 relax type (INTEGER) = set to 1 to relax to a specified streamfunction field (PSI relax.bin), set to 2 to relax to a specified PPV field (QGPV relax.bin), or set to 3 to do both 1 and 2 reset period (INTEGER) = number of steps between successive resetting of the mean PPV to zero -set to 1 to reset at each step internal ekman (INTEGER) = internal Ekman layers switched on? (set to 1 for 'yes') noise amp (REAL) = (starting value of) stochastic forcing amplitude (s −2 ) d dt noise amp (REAL) = rate of change of stochastic forcing amplitude (s −3 )

Shell script
A shell script, run qgam, is included with the model. To launch the model, type run qgam at the command line from a directory containing a namelist file. This deletes any data files already present, copies the current version of the executable to the local directory, creates a temporary uncommented namelist file, creates a file of parameter values in a form suitable to be read by the Matlab diagnostic script, runs the model (piping system messages to the output file, qgam.out), and finally deletes the local version of the executable and the temporary namelist.
To avoid deleting pre-existing data files (e.g. for a pick-up run), use run qgam pu.

Matlab diagnostics
A comprehensive Matlab diagnostics package (diagnostic.m, diagnostic read.m and gradient imp.m) is supplied with the model, consisting of 2100 lines of code. The package allows the model data and many other derived quantities to be plotted in cylindrical geometry. To run the package, launch Matlab from the data directory and type diagnostic. The data file created by the shell script is read, and the following options are offered. The user should select the option required and follow the on-screen instructions to display the plot on the screen. To save or print a figure, first produce it using the appropriate option, then exit the diagnostics package using option 0 and issue the appropriate Matlab print command.

Historical development
The historical development of QUAGMIRE has proceeded as follows. QUAGMIRE v1.0 did not include topographic forcing or stochastic forcing, and there were only 21 plot options in the Matlab diagnostics software. In QUAGMIRE v1.1, stochastic forcing was introduced, and 9 further options were added to the diagnostics software. In QUAGMIRE v1.2, topographic forcing was introduced, and 2 further options were added to the diagnostics software. Finally, in QUAG-MIRE v1.3, many improvements to the source code, shell script and diagnostics software have been made. QUAG-MIRE v1.3 is the first version to be released for public use. Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

Summary
The QUAGMIRE model described herein has been run extensively, and a detailed comparison between model and laboratory flows has been carried out (Williams, 2003;Williams et al., 2003Williams et al., , 2004Williams et al., , 2005. The code is very efficient: on a desktop Linux workstation with a 1.4 GHz AMD Athlon processor and 100% of the CPU usage, and with N azim =96 and N rad =16, a model integration speed of 120 time steps per second is attained. With these specifications, simulated time runs ten times faster than elapsed time. Waves in the model, which grow due to baroclinic instability if the Froude number is supercritical but otherwise decay, have phase speeds, equilibrated amplitudes and wavenumbers that agree well with those determined from the corresponding laboratory experiments. For Froude numbers that are higher still, more complicated model flows result, such as amplitude vacillations with reasonable amplitudes and periods and, ultimately, flow that is highly irregular and appears to be chaotic. The good agreement between model and laboratory provides an important validation of the model, and indicates that the numerical techniques employed are reliable.
As with any numerical model, many improvements could be made to version 1.3 of QUAGMIRE. The most obvious would be to generalize the model to apply to an unspecified number of superposed fluid layers of unspecified relative depths, rather than the implementation in version 1.3 of two layers of equal resting depths. There are plans to implement this improvement, and others, in future versions of QUAG-MIRE.