Total energy and potential enstrophy conserving schemes for the shallow water equations using Hamiltonian methods – Part 1: Derivation and properties

The shallow water equations provide a useful analogue of the fully compressible Euler equations since they have similar characteristics: conservation laws, inertiagravity and Rossby waves, and a (quasi-) balanced state. In order to obtain realistic simulation results, it is desirable that numerical models have discrete analogues of these properties. Two prototypical examples of such schemes are the 1981 Arakawa and Lamb (AL81) C-grid total energy and potential enstrophy conserving scheme, and the 2007 Salmon (S07) Z-grid total energy and potential enstrophy conserving scheme. Unfortunately, the AL81 scheme is restricted to logically square, orthogonal grids, and the S07 scheme is restricted to uniform square grids. The current work extends the AL81 scheme to arbitrary non-orthogonal polygonal grids and the S07 scheme to arbitrary orthogonal spherical polygonal grids in a manner that allows for both total energy and potential enstrophy conservation, by combining Hamiltonian methods (work done by Salmon, Gassmann, Dubos, and others) and discrete exterior calculus (Thuburn, Cotter, Dubos, Ringler, Skamarock, Klemp, and others). Detailed results of the schemes applied to standard test cases are deferred to part 2 of this series of papers.


Introduction
Consider the motion of a (multi-component) fluid on a rotating spheroid under the influence of gravity and radiation. This is the fundamental subject of inquiry for geophysical fluid dynamics, covering fields such as weather prediction, climate dynamics, and planetary atmospheres. Central to our current understanding of these subjects is the use of numerical models to solve the otherwise intractable equations (such as the fully compressible Euler equations) that result. As a first step towards developing a numerical model for simulating geophysical fluid dynamics, schemes are usually developed for the rotating shallow water equations (RSWs). The RSWs provide a useful analogue of the fully compressible Euler equations since they have similar conservation laws, many of the same types of waves, and a similar (quasi-) balanced state. It is desirable that a numerical model possesses at least some these same properties (see Fig. 1, and the discussion in Staniforth and Thuburn, 2012).
In fact, there exists some evidence (Dubinkina and Frank, 2007) that schemes without the appropriate conservation properties can fail to correctly capture long-term statistical behavior, at least for simplified models without any dissipative effects. However, questions remain as to the relative importance of various conservation properties for a full atmospheric model, especially in the presence of forcing and dissipation (Thuburn, 2008). This subject deserves further study, but a key first step is the development of a numerical scheme that possesses the relevant conserved quantities, and is capable of being run at realistic resolutions on the types of grids that are used in operational weather and climate models.
A pioneering scheme developed over 30 years ago possesses many of these properties (including both total energy and potential enstrophy conservation): the 1981 Arakawa and Lamb scheme (AL81, Arakawa and Lamb, 1981). Unfortunately, this scheme is restricted to logically square, orthog-onal grids (an orthogonal grid has a dual grid with edges orthogonal to the primal grid edges) such as the latitudelongitude (lat.-long.) or conformal cubed-sphere grid. These grids are not quasi-uniform under refinement of resolution, and this leads to clustering at typical target resolutions for next generation weather and climate models (such as 2-3 km for weather; and 10-15 km for climate). Such clustering will introduce strong CFL limits, and in the case of the lat.-long. grid requires polar filtering (which is not scalable on current computational architectures) in order to take practical time steps. For these reasons, it is desirable to be able to use quasi-uniform grids such as the icosahedral grid (orthogonal but non-square) or gnomic cubed sphere (square but non-orthogonal). In addition to the restriction to logically square, orthogonal grids, the AL81 scheme also suffers from poor wave dispersion properties when the Rossby radius is under-resolved (Randall, 1994). In fact, the unavoidable averaging required for the Coriolis term in a C-grid scheme is expected to lead to poor wave dispersion properties for an under-resolved Rossby radius regardless of the specific discretization employed.
Recently, there has been an effort to extend the AL81 scheme to more general grids, using tools from discrete exterior calculus (commonly referred to as the TRiSK scheme; Thuburn et al., 2009;Ringler et al., 2010;Thuburn and Cotter, 2012;Weller, 2014;Thuburn et al., 2014). This has lead to the development of a family of schemes on general nonorthogonal (spherical) polygonal meshes that possess all of the desirable properties of AL81 except for extra modes branches on non-quadrilateral meshes, which are unavoidable for C-grid schemes, and lack of either total energy or potential enstrophy conservation. It is possible to obtain one or the other, but not both at the same time. Along different lines, Salmon (2004) showed that AL81 and other doubly conservative schemes (such as Takano and Wurtele, 1982) are all members of another family of schemes on logically square, orthogonal meshes. This was done using tools from Hamiltonian methods, which are an area of active research in atmospheric model development.
Rather than using finite differences, it is also possible to extend the AL81 scheme to arbitrary grids by using compatible finite elements (see McRae and. Methods based on compatible finite elements are capable of achieving most of the properties listed in Fig. 1, and with the appropriate selection of spaces and careful mass lumping can have good linear properties as well (see Cotter and Shipton, 2012, Staniforth et al., 2013, and Melvin et al., 2014. However, these methods require the solution of elliptic problems at every time step, and so it is still useful to investigate a finite difference-based approach. As an alternative to the AL81 scheme, which preserves many of its valuable mimetic properties, but has good wave dispersion properties independent of Rossby radius, Randall (1994) introduced a scheme for uniform square grids based on the vorticity-divergence formulation (termed the Z grid) of the continuous equations. Subsequently, this approach was extended to arbitrary (spherical) orthogonal polygonal grids with a triangular dual in Heikes and Randall (1995a, b), which included the important case of an icosahedralhexagonal grid. Although this scheme possesses many of the desirable properties from AL81, it does not conserve total energy or potential enstrophy. However, a similar Z grid scheme based on a Helmholtz decomposition of the momentum instead of the wind, which does conserve both total energy and potential enstrophy, was developed by Salmon (2005Salmon ( , 2007 using techniques from Hamiltonian mechanics (specifically, Nambu brackets). The idea of using Hamiltonian mechanics to derive conservative models for atmospheric dynamical cores has seen a great deal of interest and progress in the past 10 years (see Gassmann and Herzog, 2008;Gassmann, 2013;Sommer and Névir, 2009;Nevir and Sommer, 2009;Dubos and Tort, 2014;Dubos et al., 2015;Tort et al., 2015;Salmon, 1988;Shepherd, 2003). With the recent development of Hamiltonian formulations for essentially all of the equation sets and vertical coordinates used in atmospheric dynamics, it seems likely that this approach will continue to be employed in the future. Unfortunately, the scheme in S07 is defined only for planar grids, and in the key case of general polygonal grids no expression for discrete Hamiltonian or Casimirs was given. This precludes its further development for implementation into an operational dynamical core.
This work combines the discrete exterior calculus approach from Thuburn and Cotter (2012) and the Hamiltonian approach from Salmon (2004) to extend AL81 to general non-orthogonal (spherical) polygonal grids in a manner that conserves both total energy and potential enstrophy, and to extend S07 to arbitrary (spherical) orthogonal polygonal grids. The extension of AL81 is done through the development of a new Q (the discretization of qk×, which is also known as the nonlinear potential vorticity flux) operator, using tools from Hamiltonian methods. S07 is extended by combining the Nambu bracket-based approach from Salmon (2007) with the discrete exterior calculus tools introduced in Thuburn and Cotter (2012). It should be noted that this work deals only with spatially conservative discretization. Conservation errors introduced due to time discretization are typically much smaller than those due to space discretization. However, the extension of this approach to fully conservative discretization would be a useful contribution.
The remainder of this paper is structured as follows: Sect. 2 introduces the rotating shallow water equations in both their familiar vector-invariant form and the less familiar Hamiltonian forms. Section 3 presents a family of Cgrid numerical schemes that possess many of the desirable properties, and discusses the specific member of this family introduced here. Section 4 introduces the new operator Q that enables the conservation of both total energy and potential enstrophy in the C-grid scheme. Section 5 presents the Z grid scheme and discusses its key mimetic and con-servation properties. Section 6 discusses an implementation of these schemes, and shows some limited results. Finally, some conclusions (Sect. 7) are drawn. The appendices discuss various ancillary topics such as the computational grid used (Appendix A), the specific discrete operators employed (Appendices B, C, and D), and the discrete variables used in the C-and Z grid schemes (Appendices E and F).

Rotating shallow water equations
The RSWs for both planar and spherical domains are presented below in several forms: the vector-invariant formulation, the vorticity-divergence formulation, the symplectic Hamiltonian formulation based on the vector-invariant form and both Poisson bracket, and Nambu bracket formulations based on the vorticity-divergence formulations. Although all of these formulations are equivalent in the continuous case, they lead to very different discretizations.

Vector-invariant formulation
The mass continuity equation for the RSWs is expressed in vector-invariant form as where h is the fluid height and u is the fluid velocity. Similarly, the momentum equation is expressed as where F = hu is the mass flux, q = η h is the potential vorticity, η = ζ + f is the absolute vorticity, ζ =k · ∇ × u is the relative vorticity, f is the Coriolis force, = gh + K + gh s is the Bernoulli function, h s is the topography height, and g is gravity and K = u·u 2 is the kinetic energy.

Poisson bracket formulation (vector invariant)
As discussed in Salmon (2004), let the Hamiltonian H be given by and let x = (h, u). Note that denotes the entire domain of interest, which is restricted here to either a doubly periodic plane or the sphere. Therefore, there are no boundary condition to consider. The Hamiltonian formulation can be used in the presence of boundaries, but it becomes more complicated and is not treated here. The time evolution of an arbitrary functional F can be written as where the Poisson bracket {F, H} (which is a bilinear, antisymmetric operator that satisfies the Jacobian identity) is It is useful to split this into two separate brackets as where encompasses the gradient and divergence terms and encompasses the nonlinear potential vorticity (PV) flux term. The functional derivatives δH δx of the Hamiltonian are given by This formulation is useful for development of a scheme that possesses discrete conservation properties, as discussed below. A functional derivative of some functional F[x] is defined as

Conserved quantities
Since the rotating shallow water equations form a (noncanonical) Hamiltonian system, we know from Noether's theorem and other considerations (such as the singular nature of the symplectic operator) that there are at least two categories of conserved quantities: Hamiltonian and Casimirs.

Energy (Hamiltonian)
The first is simply the Hamiltonian itself. In this case, the Hamiltonian is the total energy of the system. Conservation of the Hamiltonian arises due to the skew-symmetric nature of the Poisson bracket. In particular, using Eq. (4) the evolution of H is given by dH dt = {H, H} = −{H, H} = 0 (11) Figure 1. A diagram of some desirable model properties for the shallow water equations, organized thematically into groups. Similar considerations apply for the Euler, hydrostatic primitive, and other equation sets used in atmospheric models. There is vigorous discussion in the literature and between model designers about the importance of various properties for different applications (such as weather forecasting or long-term climate prediction). The schemes presented here satisfy all of these properties, with the exception of accuracy. There are additional desirable model properties, such as consistent physics-dynamics coupling, compatible and accurate tracer advection, and tractable treatment of acoustic waves that are not presented.
since {, } is skew symmetric. For the rotating shallow water equations, the Hamiltonian is the total energy of the system. The elegant derivation of energy conservation and its simplicity (relying ONLY on the skew symmetry of {, }) motivates the use of the Hamiltonian formulation for development of numerical schemes that conserve energy.

Casimirs
The second category of conserved quantities consists of Casimir invariants. Since the rotating shallow water equations are a non-canonical Hamiltonian system, the Poisson bracket {, } is singular and thus it possesses Casimir invariants C that satisfy for any functional F. Note that from above, this implies that For the rotating shallow water equations, the Casimirs take the form where F (q) is an arbitrary function of the potential vorticity and where ∇ ⊥ is the skew-gradient operator. On the plane it is k×∇, and it has a coordinate-independent definition on more general manifolds such as the sphere. Important cases for F include F = 1 (mass conservation), F = q (circulation or mass-weighted potential vorticity), and F = q 2 2 (potential enstrophy).
(2), we obtain the vorticity-divergence form of the equations where µ = ∇ · u is the divergence. The mass flux can then be split into rotational and divergent components (i.e., a Helmholtz decomposition) as where (hu) div = ∇χ and (hu) rot = ∇ ⊥ ψ. The stream function ψ and velocity potential χ can be related to the vorticity and divergence as where J (a, b) = ∇·(a∇ T b) = ∇ T ·(a∇b) is the Jacobian operator. The Helmholtz decomposition connects the vorticitydivergence formulation and the vector-invariant formulations. In the preceding, we have neglected the possibility of a harmonic component (a component A for which ∇ 2 A = 0), which works because the harmonic component on the sphere is zero. On the doubly periodic plane, it would be possible to have a constant harmonic component. Finally, Eqs. (1) and (2) can be re-written in terms of χ and ψ directly as

Poisson bracket formulation (vorticity-divergence)
As shown in Salmon (2007), the preceding equations (Eqs. 21, 22, and 23) can be also be written in terms of a Poisson bracket. Let x = (h, ζ, µ) and define the Hamiltonian Note that which gives (this is the functional derivative of the Hamiltonian with respect to x), and define a Poisson bracket (which is bilinear, anti-symmetric and satisfies the Jacobian identity) as where for arbitrary functionals A and B. As before, the time evolution of an arbitrary functional A is then given by It is easy to see that Eqs. (21), (22), and (23) are recovered when A is set equal to h, ζ or µ, respectively. Note that each of the brackets in Eqs. (29), (30), and (31) are antisymmetric, and that the Casimirs C = hF (q)d satisfy {A, C} = 0 (where F is an arbitrary function and A is an arbitrary functional) independently for each bracket. The use of the Poisson (and Nambu) bracket formulation of the shallow water equations is motivated by the intimate connection between these formulations and the conserved quantities. As is well known, the conservation of energy H rests solely on the anti-symmetry of the Poisson bracket, and a numerical scheme that retains this feature will automatically conserve energy. However, potential enstrophy is a Casimir, and therefore developing a numerical scheme using the Poisson formulation that conserves it requires that the discrete potential enstrophy lies in the null space of the resulting discrete bracket. This can be difficult, especially on arbitrary grids, and this motivates the use of a continuous formulation that does not contain a null space, which is discussed below.

Nambu bracket formulation (vorticity-divergence)
Fortunately, there is a closely related formulation of the shallow water equations in terms of Nambu brackets (see Salmon, 2007): where cyc is a cyclic permutation, Z = d h q 2 2 is the potential enstrophy, and the multi-part dot product is simply the product of the individual components, summed over each basis (for example, in two-dimensional (2-D) doubly periodic flow the first term is . The time evolution of an arbitrary functional A is now given by These brackets are useful because they are triply antisymmetric (which ensures the conservation of H and Z) and non-degenerate (they have no Casimirs). In fact, discrete conservation of both total energy and potential enstrophy requires only the triply anti-symmetric nature is retained. It is also possible to generalize these brackets to any Casimir (as shown in Salmon, 2005), but since we are interested mostly in potential enstrophy conservation this is not necessary. These brackets will form the basis of the Z grid discretization method discussed below.

C-grid scheme
Following Thuburn and Cotter (2012), the prognostic variables for the C-grid scheme are the fluid height integrated over a primal grid cell (equivalent to the mass in that grid cell) M and the wind integrated over a dual edge (equivalent to the circulation along a dual edge)Û . Letting x = (M,Û ), the vector-invariant Poisson bracket can be discretized in a manner that preserves its anti-symmetric character (which ensures total energy conservation) and a subset of the Casimir invariants (specifically mass, potential vorticity, and potential enstrophy). Combined with a choice for the discrete Hamiltonian, this constitutes a complete discretization for the nonlinear rotating shallow water equations. Ideally, one would use a Nambu bracket formulation of the vector-invariant shallow water equations rather than the Poisson bracket formulation in order to avoid the difficulties associated with developing a discretization that has the correct Casimirs, since in the Nambu bracket case only antisymmetry must be enforced. Unfortunately, the only known Nambu bracket for the vector-invariant shallow water equations possesses intractable singularities and is not suitable as the basis for developing a discretization (Salmon, 2005). In what follows, uppercase letters will denote the entire (column) vector of degrees of freedom, while lowercase letters will denote a specific degree of freedom. A hat on a variable indicates that the quantity is defined on the dual grid.
Specifically, the brackets in Eqs. (7) and (8) are discretized using the operators from Appendices C and B as where the functional derivatives are simply the partial derivatives with respect to the appropriate quantity. Note that these discrete brackets are only bilinear and anti-symmetric, they do not satisfy the Jacobian identity. In addition, they possess only a subset of the Casimirs of the continuous brackets. Therefore, they should be properly be termed quasi-Poisson brackets. The brackets given in Eqs. (37) and (38) are essentially a generalization of the brackets introduced in S04 from uniform square grids to arbitrary polygonal grids, using the operators from Thuburn and Cotter (2012). The Hamiltonian H is discretized as where g is the acceleration due to gravity, B is the topographic height integrated over a primal grid cell,Ĉ = M eÛ is the mass flux on dual edges, and M e = φIM, with φ an interpolation operator. Taking functional derivatives yields whereˆ is the Bernoulli function sampled at dual vertices and F is the mass flux integrated over primal edges. Computing actual values yieldsˆ = I(K + gM + gB) with K = φ TÛ ×HÛ 2 , where K is the kinetic energy integrated over primal grid cells and F = HĈ. A detailed description of these discrete variables and their staggering on the computational grid can be found in Appendix E, and a diagram of their staggering is in Fig. 2. The resulting discrete evolution equations are In fact, by making alternative choices for F , Q, andˆ (along with the operators discussed below) it is possible to recover a wide range of C-grid schemes present in the literature (such as Ringler et al., 2010, and Weller, 2014; see Eldred, 2015 for more details). The operators D 2 , D 1 , D 1 , D 2 , I, J, R, W, and H are defined in Appendices B and C (and can also be found in a general form in Thuburn and Cotter, 2012). The novelty of the current scheme is a new definition of Q, such that the properties of total energy conservation, potential enstrophy conservation, and steady geostrophic modes hold simultaneously. This is the subject of Sect. 4. A potential vorticity equation can be obtained from Eq. (42) by taking D 2 to yield where D 2 D 1 = 0 has been used to eliminate the gradient term, and D 2Û =ˆ is the relative vorticity integrated over dual grid cells;Υ =ˆ + f =M × Q, whereM = RM is the mass integrated over dual grid cells,Υ is the absolute vorticity integrated over dual grid cells, f is the Coriolis parameter integrated over dual grid cells, and Q is the potential vorticity sampled at primal grid vertices.

Relationship to discrete exterior calculus
As discussed in Thuburn and Cotter (2012), these operators have an interpretation in terms of discrete exterior calculus. In fact, D 2 , D 1 , D 2 , and D 1 are discrete exterior derivatives, I, J, and H are Hodge stars and the various prognostic and diagnostic quantities can be interpreted as discrete differential forms. This connection is further explored in Eldred (2015).

Linearized scheme
As is well known, the linearized version of a Hamiltonian system about a steady state can be found by evaluating the brackets at that state and using the quadratic approximation to the associated pseudo-energy as the Hamiltonian (Shepherd, 1993). Following this procedure and letting the Coriolis parameter f be a constant, B = 0, and assuming a background state of x = (H, 0), we obtain for the brackets (where W = Q q v =1 is the linearized version of Q) and for the Hamiltonian, which has associated functional derivatives of The resulting evolution equations are Figure 2. A subset of discrete variables and their staggering on the computational grid for the C-grid scheme. A subscript i indicates quantities defined at primal grid cells or dual grid vertices, a subscript e indicates quantities defined at primal or dual grid edges, and a subscript v indicates quantities defined at primal grid vertices or dual grid cells. The prognostic (red) quantities are the mass integrated over primal grid cells m i and the wind integrated along dual grid edges u e , the other quantities are diagnostic (blue). More details can be found in Appendix E.

Properties of scheme
This scheme has many important properties, including the following: 1. Mass and potential vorticity conservation: both mass M and mass-weighted potential vorticityM × Q are conserved in both a local (flux-form) and global (integral) sense.
2. No spurious vorticity production: by construction, D 2 D 1 = 0 and there is no spurious production of vorticity due to the gradient term in the wind equation.
3. Linear stability (pressure gradient force and Coriolis force conserve energy): this is due to the fact that I, J, and H are all symmetric positive definite; D T 2 = −D 1 , D 2 T = D 1 , and W = −W T . 4. Steady geostrophic modes: by construction, −RD 2 = WD 2 (noting that W is the same for all members of this family), which gives steady geostrophic modes.
5. PV Compatibility: again by construction −RD 2 = WD 2 with Q q v =c → cW, and therefore the potential vorticity equation is compatible with the diagnostic mass equation (a constant PV field remains constant). Note that this is same as the condition required for steady geostrophic modes.
6. Other conservation properties: see below for a discussion on total energy and potential enstrophy conservation. Table 1. Summary of required operator properties for obtaining the desirable mimetic properties along with total energy and potential enstrophy conservation. A example of operators that satisfy these properties can be found in Appendix C. More details can be found in Thuburn and Cotter (2012) or Eldred (2015). Note that the mapping column indicates which types of quantities the operator accepts as inputs, and what it produces as output, with "p" or "d" denoting the primal and dual grids; while the number (0,1,2) denotes the geometric entity the quantity is integrated over. For example, the R operator takes as input quantities integrated over primal grid cells, and produces as output quantities integrated over dual grid cells.  Table 1 shows a summary of the required properties in order for the resulting scheme to have all of the mimetic and conservation properties discussed above.

Total energy conservation
Following S04, total energy will be conserved for any choice of H if the discrete brackets retain their anti-symmetric character. This requires that D T 2 = −D 1 , and that Q = −Q T . The first condition is satisfied by construction of the discrete derivative operators D 2 and D 1 . The second condition is satisfied only for certain choices of Q. One example is Q = 1 2 Q e W + 1 2 WQ e (as used in Ringler et al., 2010), where Q e is any function that, given the set of q v at primal vertices, computes a unique Q e at primal edges (such as Q e = 1 2 v∈VE(e) q v ). Flexibility in the choice of Q e allows for a wide variety of stabilization methods such as CLUST or APVM (Weller, 2012;Weller et al., 2012). Unfortunately, this choice does not conserve potential enstrophy.

Potential enstrophy conservation
Following S04, potential enstrophy is a Casimir and therefore will be conserved when holds for any choice of functional A. Potential enstrophy is defined as Its functional derivatives are Using the chain rule for functional derivatives, it suffices to show that Eq. (50) holds for A = i m i and A = e u e . Therefore, Eq. (50) reduces to which must hold for any choice of Q. The first of these is again satisfied by construction for D 2 and D 1 . The second is much trickier, and is the main subject of Sect. 4. One example is Q = Q e W (as used in Ringler et al., 2010), where Q e = 1 2 v∈VE(e) q v . Unfortunately, this choice does not conserve total energy. It would be possible to explore alternative definitions of Z, but these would lead to different, less natural stencils for Q.

Arakawa and Lamb 1981
In the case of a uniform square grid, the C-scheme grid above reduces to the well-known Arakawa and Lamb 1981 total energy and potential enstrophy scheme (modified to prognose m i and u e ) if their choice of Q is used. Unfortunately, the definition of Q presented in AL81 works only for logically square, orthogonal grids. For more general, non-orthogonal polygonal grids, a new operator Q must be found. This is the subject of the next section.

Hollingsworth instability
Since this is an extension of Arakawa and Lamb 1981 scheme, it seems extremely likely that the proposed scheme will suffer from the Hollingsworth instability, especially if applied in a height coordinate framework using a Lorenz staggering in the vertical (as discussed in Bell et al., 2016 andHollingsworth et al., 1983). However, other similar schemes have been able to mitigate the Hollingsworth instability when used with an isentropic or Lagrangian vertical coordinate, or when a Charney-Phillips staggering is used in the vertical. On a uniform square grid using AL81, it is possible to rigorously modify the kinetic energy stencil to eliminate the non-cancellation error that is at the heart of the instability. Furthermore, this modification can be done in such a way as to conserve total energy, by expressing it as a modification to the Hamiltonian H itself and then deriving the associated consistent mass flux F e . A similar modification of the kinetic energy stencil for a similar C-grid scheme on nonsquare grids (Gassmann, 2013) has been shown to mitigate the Hollingsworth instability even without rigorous elimination of the non-cancellation error. Therefore, given the many possible mitigation strategies, the possible presence of the instability is not expected to prevent use of this scheme in a model solving the hydrostatic or non-hydrostatic equations.

Operator Q
The principal novelty of the new C-grid scheme is the specification of a Q operator that simultaneously conserves total energy and potential enstrophy, and also supports PV compatibility. Previous work found choices for Q that conserved either total energy or potential enstrophy, but not both. The key lies in S04, showing that the AL81 approach could be extended to more general stencils (although retaining a logically square, orthogonal grid). This work takes the Salmon 2004 approach in a different direction, keeping the same stencil as AL81 but considering a general polygonal grid.

Definition of Q
Loosely following S04, Q is defined as where the weights α e,e ,v are themselves a linear combination of the potential vorticity q v at vertices v ∈ VC(i) (blue), and i is the cell shared between edges e (green) and e (red). By choosing the weights α e,e ,v appropriately, an operator Q can be found that simultaneously conserves both total energy and potential enstrophy and supports steady geostrophic modes.
where i is the primal grid cell covered by both e and e . A diagram of this operator is shown in Fig. 3. An equivalent alternative form for Q given in terms of the Poisson bracket that closely mimics the one found in S04 can be found in Eldred (2015). It is easy to see that in the case of a logically square, orthogonal grid, this approach reduces to the same stencil considered by AL81. At this point, the coefficients α e,e ,v are undetermined.

Linear system for α
It remains to determine the coefficients α e,e ,v in a manner such that the resulting operator Q conserves both total energy and potential enstrophy, and satisfies PV consistency.

Requirements introduced by energy conservation
Following S04, in order for Q to be energy conserving then Q = −Q T . In terms of the coefficients, this implies that α e,e ,v = −α e ,e,v , or in other words, they are anti-symmetric under an interchange of e and e .

Requirements introduced by potential enstrophy conservation
From Eq. (54), in order for Q to conserve potential enstrophy, −D 1 R Q×Q 2 + QD 1 Q = 0 must hold for any choice of Q. Expanding this out yields for every e, which must hold for any choice of q v . For a given edge e, the vertices in question are v ∈ CVE(e) (shown in Fig. 4), where CVE(e) = VE(i1) ∪ VE(i2) and (i1, i2) = CE(e). Both the left-and right-hand side of these equations are a quadratic form in this set of vertices, and for this to hold for arbitrary q v the coefficients in these two quadratic forms must be equal. These coefficients are linear combinations of the α's, and therefore the equality of these quadratic forms implies a set of linear equations for the α's. Specifically, for each grid cell i with n e edges and n v vertices (note that n e = n v for a polygonal grid cell, but it is useful to keep distinct notation to ease exposition), there are n e n v (n v +1) 2 equations (coefficients in the quadratic forms) and n v n e (n e −1) 2 unknowns (the coefficients α e,e ,v ). Since n v = n e , this is therefore an overdetermined system, and the coefficient will be found through a least-squares procedure. The equations come from equating the coefficients in the two quadratic forms: there are n v (n v +1) 2 independent vertex pairs, and n e edges. The unknowns are the coefficients α e,e ,v that are associated with the grid cell: there are n e (n e −1) 2 independent unique edge pairs, and n v vertices. Note that this has already taken into account the fact that α e,e ,v = −α e ,e,v (hence the wording unique edge pair), which reduces the number of independent coefficients in half. Letting v and v loop over the vertices in the cell (they are the unique members of VC(i) × VC(i)), the equations are given by where the sum for B v,v occurs only when v ∈ VE(e); and where e loops over each edge in i and EVE(v, e, i) = EC(i)∩ EV(v) − e; and sgn(e, e ) = 1 = −sgn(e , e) (which ensures that the scheme is also energy conservative). A diagram of EVE(v, e, i) is provided in Fig. 5. Note that coefficients in one cell are coupled with adjacent cells when v ∈ VE(e) or v ∈ VE(e); that is to say, the equations involve coefficients that are associated with other grid cells. On a non-uniform mesh, this means that the entire set of coefficients must be solved at the same time.
The solution procedure outlined above gives a large matrix system where each row in A represents an equation obtained by equating coefficients in the quadratic forms, and α is the vector of unknown coefficients. This system can be solved (via a least-squares approach) to yield a set of coefficients α such that Q conserves potential enstrophy, provided the system can be solved exactly. This procedure is essentially identical to the one employed in S04. In addition, the coefficients can be computed once, and then stored for later use. Unfortunately, the system that results directly from this procedure is impractical to solve for realistic non-uniform meshes: it is too large and ill conditioned. For example, on an icosahedral-hexagonal mesh with O(1 million) grid cells, there will be O(90 million) coupled coefficients that need to be solved for.

Practical solution
Instead, following Thuburn et al. (2009), the coefficients can be uncoupled by defining when v ∈ VE(e) or v ∈ VE(e), where C = −1/6. This will produce an independent subsystem for each grid cell, with for example 90 unknowns on a hexagonal grid cell and 24 unknowns on a square grid cell. When this procedure is applied to a uniform square grid it reproduces the AL81 scheme, and on a uniform hexagonal grid it produces a total energy and potential enstrophy conserving scheme (not shown, verified numerically). In all cases, including the nonuniform meshes tested (icosahedral and cubed-sphere grids), the least-squares problem is solved exactly, in the sense that the coefficients exactly satisfy the relationships for potential enstrophy and total energy conservation. Since each subsystem is still overdetermined, this implies the existence of an associated solvability condition. It seems likely that the solvability condition is the key to writing down an explicit formula for the coefficients in terms of R i,v and n e,i . Unfortunately, we were unable to derive such a condition. However, this does not prevent the numerical solution of the leastsquare problems, which is sufficient for practical use of the scheme. We were able to solve the systems on cubed-sphere meshes with up to 884 736 grid cells and on icosahedralhexagonal meshes with up to 655 363 grid cells in a few hours using an unoptimized, serial algorithm on a laptop computer. Furthermore, the uncoupled nature of the problem (one small independent least-squares problem per grid cell) would facilitate trivial parallelism if needed for larger meshes.

PV compatibility
The astute reader will note that nothing has been said yet about enforcing PV compatibility (Q q v =c = cW. It was originally believed that PV compatibility would have to added as additional equations in the matrix-vector system. However, it was found that enforcing potential enstrophy conservation (even using the uncoupled form) was sufficient to ensure that Q was PV compatible. This corresponds with the results of S04 (Salmon, 2004), who did not explicitly add PV compatibility, yet all of his schemes had this property.
for every edge pair (e, e ); which could be easily added to the independent system of equations solved in each grid cell.
Although enforcing Q q v =c = cW ensures PV compatibility, it also requires the use of the W operator from Thuburn et al. (2009). Therefore, Q will share the same limitations, including inconsistency on general grids. The consequences of this are explored more in Eldred and Randall (2017a).

Z grid scheme
Unlike the C-grid scheme, the Z grid scheme starts with Nambu brackets rather than Poisson brackets. This greatly simplifies the derivation, since only the triply anti-symmetric nature of the brackets must be retained to ensure total energy and potential enstrophy conservation: there is no consideration of Casimirs. Start by defining a set of collocated discrete variables which are pointwise values of h, ζ , and µ at primal grid centers. We will also use the More details about the grid; discrete operators and discrete variables can be found in Appendices A, D, and F.

Functional derivatives
The functional derivative of a general functional F with respect to discrete variable x i is then defined as where A i is the area of primal grid cell i. The diagnostic variables i , χ i , ψ i , and q i are defined through the functional derivatives of the discrete Hamiltonian H and discrete Potential Enstrophy Z as At this point the discrete Hamiltonian H and discrete potential enstrophy Z are left unspecified.

Discrete Nambu brackets
Following Salmon (2007), the general discretization starts from the Nambu brackets in Eqs. (33), (34), and (35) for the shallow water equations in vorticity-divergence form. As long as these brackets retain their triply anti-symmetric structure when discretized, total energy and potential enstrophy will be automatically conserved for any definition of the total energy and potential enstrophy (with one caveat explained below). In addition, the bracket structure ensures that this conservation is local as well as global. That is, the evolution of a conserved quantity can be written in flux form for each grid cell, where cancellation of fluxes between adjacent cells leads to the global integral being invariant. This is in contrast to a method that conserves the global integral, but cannot be written in flux form for each grid cell. In what follows below, we will consider only the case where Z is the potential enstrophy, although this approach could be easily generalized to arbitrary Casimirs (see Salmon, 2005 for an example of this on a uniform square grid). In discretizing the Nambu bracket, the operators D 1 , D 2 , D 1 , and D 2 from the C-grid scheme are needed. In addition to these, the additional operators J (A, B), K, X v , and X e are also needed, and they are given in Appendix D.

Jacobian brackets
Loosely following S07, the {F, H, Z} ζ ζ ζ bracket can be discretized as Note that this bracket is triply anti-symmetric (due to the cyclic permutation), as required. The {F, H, Z} µµζ bracket can be similarly discretized as This bracket is only doubly anti-symmetric (in H and F due to the anti-symmetry of J ), but it will conserve Z as well provided that δZ δµ i = 0 (since J (A, B) = 0 when either A = 0 or B = 0). These brackets are essentially those encountered when discretizing the Arakawa Jacobian, as detailed in Salmon (2005).

Mixed brackets
The mixed bracket is trickier since it contains an apparent singularity 1 ∇q . On closer inspection, in the continuous case this singularity cancels out when combined with the functional derivative of the potential enstrophy. This is the caveat mentioned above; i.e., the discrete mixed bracket must be constructed such that the apparent singularity cancels out with the discrete functional derivative of the potential enstrophy. With this in mind, the general form of the discrete mixed bracket is chosen as where, from before, q i ≡ δZ δζ i . The quantities le and de are the edge lengths on the primal and dual grid, as defined in Appendix A. This bracket is triply anti-symmetric (again due to the cyclic permutation), and the apparent singularity will cancel if Z is chosen with care.

Conservation
Since the {F, H, Z} ζ ζ ζ and {F, H, Z} µζ h brackets are triply anti-symmetric, and the {F, H, Z} ζ µµ bracket is doubly antisymmetric, both total energy and potential enstrophy will be conserved for any choice of H and Z, provided that the caveats mentioned above are obeyed. Those are 1. δZ δµ i = 0 (ensures that the {F, H, Z} ζ µµ bracket conserves potential enstrophy) 2. Z chosen such that the apparent singularity ( D 1 (Z h ) D 1 q i term + cyc(F, H, Z) terms) in the {F, H, Z} µζ h bracket cancels out .
These are fairly minimal requirements, and many reasonable choices for Z satisfy them.

Discrete Hamiltonian and Helmholtz decomposition
The Hamiltonian H can be split into three parts: H FD , H J , and H PE , where the first two are the kinetic energy due to flux-divergence terms and Jacobian terms, and the last is the potential energy. In the continuous system we have where These can be discretized as

Helmholtz decompositions and Bernoulli function
By taking variations of H, we obtain − edges le de After a lot of algebra, these can be grouped (half of each term involving δh i goes to i and half to µ i /ζ i ) to obtain where (using the definition of functional derivative) The latter two equations (Eqs. 86 and 87) are the discrete version of the Helmholtz decomposition, and form a pair of non-singular elliptic equations. They can be combined into a single equation as where, for example, Note that (without the 1 A i factors) FDH is symmetric and JA is anti-symmetric, which means that A = −A T (i.e., A itself is skew symmetric). Also note that when h i = H is a constant (and therefore h e = H ), they reduce to where L = 1 A i D 2 le de D 1 , which is the expected linearization behavior.

Discrete potential enstrophy
A natural definition of the discrete potential enstrophy is where η i = ζ i + f i . Taking variations of this yields Then the natural definition for q i = η i h i works, and the above simplifies to By plugging these back into the {F, H, Z} µζ h bracket, it is seen that this choice of Z also ensures that the singularity cancels.

Independence between choices for H/Z and Nambu Brackets
As noted before, the mimetic and conservation properties of the discrete scheme are completely independent of the choice of discrete Hamiltonian H, provided the Hamiltonian is positive definite and produces invertible elliptic equations for the Helmholtz decomposition. If the resulting elliptic equations were singular, then the scheme would have a computational mode (as discussed in Salmon, 2007). Additionally, the discrete Helmholtz decomposition should also simplify to a pair of uncoupled Poisson problems when linearized. The mimetic and conservation properties are also independent of the specific choice of Z, provided that the singularity in the mixed bracket cancels and Z µ = 0. The given choices of H and Z were selected to have these properties, and also correspond with those in S07 for the special cases of a uniform planar square grid and an orthogonal polygonal planar grid with a triangular dual.

Discrete evolution equations
By setting F = (h i , ζ i , µ i ) in turn, the following evolution equations are obtained: where L is the Laplacian, FD is the flux divergence and J is the Jacobian. On an icosahedral hexagonal-pentagonal grid these operators are the same as those from Heikes and Randall (1995a), and will therefore share the same limitations as those operators, such as the inconsistency of the Jacobian on general grids. The consequences of this are explored more in Eldred and Randall (2017a). The differences between the schemes arise from the use of different arguments to the operators (q i instead of η i ) and the use of different definitions for χ i and ψ i (which in turns induces a different Poisson problem and different expression for i ).

Laplacian and flux-divergence operators
The Laplacian and flux-divergence operators (which come from the mixed bracket) can be written as where α e = i∈CE(e) α i 2 .

Jacobian operators
The Jacobian operators (which come from the Jacobian brackets) can be written as Note that on a polygonal grid with a purely triangular dual (including the important case of an icosahedral grid), J δ = J ζ .

Linearized version
Under the assumption of linear variations around a state of rest (h i = H , ζ i = µ i = 0, q i = f H ) on a f plane, this scheme reduces to where the Helmholtz equations given by Eqs. (89) and (90) have been used to simplify the scheme (to the point that it no longer requires solving any elliptic equations). In the case of a uniform square grid (uniform hexagonal grid), this scheme is identical to the one studied in Randall (1994) (Ničković et al., 2002), and it shares the same excellent linear wave properties found for those schemes.

Relation to Salmon schemes
For the cases of a uniform planar square grid and a general orthogonal planar polygonal grid with triangular dual, the general discretization scheme presented above reduces to the schemes given in S07. However, this discretization scheme is more general, and it also makes specific choices for the total energy H and potential enstrophy Z when using a general polygonal grid.

Properties of scheme
The discrete scheme as outlined above possesses the following (among others) key properties: 1. Linear stability (Coriolis and pressure gradient forces conserve energy): provided that L = L T (which is satisfied for the L given above, and the majority of discrete Laplacians), the scheme will conserve energy in the linear case.
2. No spurious vorticity production: by construction, the pressure gradient term does not produce spurious vorticity since the curl is taken in the continuous system, prior to discretization.
3. Conservation: by construction, this scheme conserves mass, potential vorticity, total energy, and potential enstrophy in both a local (flux-form) sense and global (integral) sense.
4. PV compatibility and consistency: by inspection, the mass-weighted potential vorticity equation is a fluxform equation that ensures both local and global conservation of mass-weighted potential vorticity. In addition, an initially uniform potential vorticity field will remain uniform. This rests on the fact that J ζ (q i , ψ i ) = 0 and FD(q i , χ i ) = cLχ i when q i = c is constant.
5. Steady geostrophic modes: since the same divergence µ i appears in both the linearized vorticity and continuity equations, the scheme possesses steady geostrophic modes.
6. Linear properties (dispersion relations, computational modes): as expected, the scheme possesses the same linear mode properties on uniform planar grids as those presented in Randall (1994) and Ničković et al. (2002); and it does not have any computational modes. More details of the linear mode properties of the scheme on both uniform planar and quasi-uniform spherical grids can be found in a forthcoming paper (Eldred and Randall, 2017b).
7. Accuracy: unfortunately, as shown in Heikes et al. (2013), the Jacobian operator as given is inconsistent on general grids. Even more unfortunately, the fix proposed in that paper breaks key properties of the Jacobian necessary to retain total energy and potential enstrophy conservation. Surprisingly, as shown in Eldred and Randall (2017a), the inconsistency of the Jacobian operator does not appear to cause issues in the test cases that were run. More details on possible fixes to the accuracy issue are discussed in Eldred and Randall (2017a).
6 Implementation and results

Implementation
To test the utility of the C-and Z grid schemes developed above, they were implemented in a combination of Python as a driver language along with Fortran kernels for the numerics. Although only tested on quasi-uniform grids that admit a structured approach, for simplicity the code uses an unstructured mesh with indirect addressing. Due to this highly unoptimized implementation, no cost comparisons were made with other codes; instead, we simply note that both the C-and Z grid schemes are structurally similar to other schemes used in existing models such as MPAS, Dynamico, and UZIM, and can be expected to share similar performance characteristics.

Results
As a short preview of the more detailed results in Eldred (2015) and Eldred and Randall (2017a), a run of the Galewksy et al. (2004) test case using the C-and Z grid scheme is presented below. Third-order Adams-Bashford time stepping was used, and no dissipation beyond the inherent damping in the time scheme was applied. The C-grid scheme was run on both a cubed-sphere grid (with 884 736 grid cells, approximately 26 km resolution, t = 15 s) and an icosahedral grid (with 655 362 grid cells, approximately 34 km resolution, t = 22.5 s), whereas the Z grid scheme was run only on the icosahedral grid. The absolute vorticity at day 6 for all three schemes is shown in Fig. 6, and the results from both the C-and Z grid schemes are broadly similar to both each other and to other results in the literature. Some differences can be found in the inactive region of the jet, especially when comparing the cubed sphere to the icosahedral grid simulations. It is believed that these differences are due to the underlying grid structure, since the C-and Z grid scheme on the icosahedral grid produce the same pattern for the inactive region (and a very similar Cgrid scheme on the cubed sphere that conserves only enstrophy produces an extremely similar pattern to the fully conservative C-grid scheme on the cubed sphere, as shown in Eldred, 2015). In contrast to the results in Weller (2014), we did not encounter any issues in using the C-grid scheme on the cubed-sphere grid. Plots of the time series of total energy and potential enstrophy are available in Eldred and Randall (2017a) and Eldred (2015), and verify that the schemes are conserving both energy and potential enstrophy in the spatial semi-discretization limit.

Conclusions
This paper presents an extension of AL81 to arbitrary nonorthogonal (spherical) polygonal grids in a manner that preserves almost all of the desirable properties of that scheme (including both total energy and potential enstrophy conservation) through a new Q operator. Unfortunately, on nonquadrilateral grids such as the icosahedral grid there will be extra branches of the dispersion relationship due to a mismatch in the number of degrees of freedom in the wind and mass fields inherent to the C-grid approach. Switching from a C-grid type staggering (to an A grid staggering, for example) is undesirable for many reasons, foremost among them being the natural association of physical variables with geometric entities in a staggered grid as suggested by exterior calculus and differential geometry (see Tonti, 2014 andBlair Perot andZusi, 2014). Fortunately, other than these extra mode branches on the icosahedral grid the proposed Cgrid scheme does not possess any additional computational modes. Furthermore, extensive testing has thus far been unable to show negative impacts from this extra mode branch, especially when running full-physics simulations with realistic topography and initial conditions (John Thuburn and Bill Skamarock, personal communication, 2016). This work has also presented an extension of the total energy and potential enstrophy conserving Z grid scheme in S07 from planar grids to arbitrary orthogonal (spherical) polygonal grids, using the same toolkit of Nambu brackets and Hamiltonian methods. The restriction to orthogonal grids rather than more general non-orthogonal grids (such as a cubed sphere) is a drawback. However, the major motivations for using a cubed-sphere grid are the ability to properly balance degrees of freedom when using a staggered Cgrid methods (and therefore avoid spurious branches of the dispersion relationship), a tensor-product grid structure for spectral or finite element type methods (which ensures a diagonal mass matrix for spectral element methods and effi-cient implementation of finite element methods) and higherorder finite volume methods (enabling easy dimension splitting), and an underlying piecewise continuous coordinate system for higher-order finite volume methods (allowing extended stencils). None of these considerations apply to a Z grid method, so the restriction to icosahedral grids is not anticipated to be a significant hurdle.
A detailed comparison of the two schemes, including an analysis of the accuracy of the operators used and results from a variety of test cases, can be found in second part of this series (Eldred and Randall, 2017a). In addition, an analysis of the linear mode properties of these two schemes on various quasi-uniform grids is undertaken in the third part of this paper series (Eldred and Randall, 2017b).

Code availability
The schemes described in this manuscript have been implemented in a Python/Fortran mixed language code, and are freely available at https://bitbucket.org/chris_eldred/phd_ thesis under a GNU Lesser General Public License Version 3.
Appendix A: Discrete grid The schemes described above are designed to work on arbitrary (spherical) polygonal grids along with an associated dual grid. In the case of the C-grid scheme, the grid can be either orthogonal or non-orthogonal, while the Z grid scheme is restricted to orthogonal grids. A description of the grid framework is given in what follows.

A1 General non-orthogonal polygonal grid
Consider a (primal) conformal grid constructed of polygons (or spherical polygons). A dual grid is constructed such that there is a unique one to one relationship between elements of the primal grid and element of the dual grid: primal grid cells are associated with dual grid vertices, primal grid edges are associated with dual grid edges, and primal grid vertices are associated with dual grid cells. This grid configuration covers the majority of grids that are used in current and upcoming atmospheric dynamical cores, including cubed sphere and icosahedral grids (both hexagonal-pentagonal and triangular variants). Once the dual grid vertices have been placed, there are several important geometric quantities that are needed in order to construct the discrete operators (shown graphically in Fig. A1). Specifically, we need the primal cell area A i , the dual cell area A v , the distance between primal grid centers "le", the distance between dual grid centers "de", and the overlap areas A iv and A ie . On a planar grid, these are easily defined using the standard Euclidean metric and formulas. On a spherical grid, distances must be calculated using geodesic arcs, and areas are calculated by subdividing into spherical triangles as needed and then applying the relevant spherical area formulas. See the discussion in Weller (2014) for more details.

Appendix B: Discrete derivative operators
Following Thuburn and Cotter (2012), a set of discrete derivative operators can be defined as where n e,i is an indicator that is 1 when e is oriented out of a primal grid cell and −1 when e is oriented into a primal grid cell, and t e,v is an indicator that is 1 when e is oriented into a dual grid cell and −1 when e is oriented out of a dual grid cell. D 2 is the divergence, D 2 is the curl, D 1 is the skew gradient, and D 1 is the gradient; by the Gauss theorem these are Figure A1. The geometric quantities on a planar grid. Primal grid edge lengths are denoted as "de", dual grid edge lengths are denoted as "le", the area associated with an edge by A e , the overlap between primal grid cell i and edge e by A ie , and the overlap between dual grid cell v and edge e by A iv . Note that the same definitions can be used on a spherical grid, provided the appropriate measures are used (such as geodesic lengths for distances, and spherical polygonal areas for areas). See Weller (2014) for more details.
exact (since they operate on integrated quantities). By construction, these satisfy D 2 D 1 = 0, D 2 D 1 = 0, D T 2 = −D 1 and D 2 T = D 1 for arbitrary polygonal grids. These identities are the discrete analogues of ∇ ·∇ ⊥ = 0, ∇ ⊥ ·∇ = 0, and adjointness between divergence and gradient, and curl and skew gradient. The operators can also be identified as the discrete exterior derivative operators from discrete exterior calculus.

Appendix C: Specific choices for various C-grid operators
In order to close the C-grid scheme presented in Sect. 3, specific choices must be made for I, J, H, R, φ, and W. The ones used here (and in Ringler et al., 2010 andThuburn et al., 2014) are where I, J, and H 0 are diagonal matrices, H O is used on orthogonal grids such as the icosahedral grid and H NO is used on non-orthogonal grids such as the cubed-sphere grid. The details of the construction of this operator, including the stencil S(e) and the weights H e,e , can be found in Thuburn et al., 2014). The weights W e,e are chosen such that W = −W T and −RD 2 = D 2 W (the details for this operator can be found in Thuburn et al., 2009). I, H, and J transform quantities between the primal and the dual grids, and are in fact what are known as discrete Hodge star operators. R is a discrete analogue of the identity operator that maps quantities integrated over primal cells to quantities integrated over dual cells, while W can in fact be identified as a discrete analogue of the interior product (or contraction) operator. On an orthogonal grid, the choices given for I, J, and H correspond to the Voronoi Hodge star from discrete exterior calculus.

Appendix D: Specific choices for various Z-grid operators
For the Z grid scheme, the following operators are needed: which is the sum of edges for a given grid cell, and J (A, B) = n e,2 A 2 B 1 + n e,1 A 1 B 2 , which is used to build a discrete Jacobian operator. Note that J (A, B) is anti-symmetric (J (A, B) = −J (B, A)) and satisfies J (A, 0) = J (B, 0) = J (A, A) = 0. In addition, two different interpolations (from cell centers to vertices and to edges, respectively) are defined: where C-is a constant given by 1 n , and n is the size of CV(v) (equal to 4 for quadrilateral dual grid cells and 3 for triangular dual grid cells). Appendix E: Discrete variables for the C-grid scheme Table E1 gives the discrete variables used in the C-grid scheme, their type (which indicates the staggering on the grid), and their diagnostic equation (where applicable). For the type, the first designator indicates the location on the grid type (primal or dual) and the second designator indicates the degree of the geometric entity the quantity is integrated over (0,1 or 2). For example, C-is a quantity on the dual grid integrated over edges.
Appendix F: Discrete variables for the Z-grid scheme Table F1 gives the discrete variables used in the Z grid scheme and their type (either prognostic or diagnostic).