A Lagrangian vertical coordinate version of the ENDGame dynamical core. Part I: Formulation, remapping strategies, and robustness

Previous work provides evidence that Lagrangian conservation and related properties of a numerical model dynamical core can be improved by the use of a Lagrangian or quasi‐Lagrangian vertical coordinate (LVC). Most previous model developments based on this idea have made the hydrostatic approximation. Here the LVC is implemented in a non‐hydrostatic compressible Euler equation dynamical core using almost identical numerical methods to ENDGame, the operational dynamical core of the Met Office atmospheric Unified Model. This enables a clean comparison of LVC and height‐coordinate versions of the dynamical core using numerical methods that are as similar as possible. Since Lagrangian surfaces distort over time, model level heights are continually reset to certain “target levels” and the values of model fields are remapped onto their new locations. Different choices for these target levels are discussed, along with remapping strategies that focus on different conservation or balance properties. Sample results from a baroclinic instability test case are presented. The LVC formulation is found to be rather less robust than the height‐coordinate version; some reasons for this are discussed.


INTRODUCTION
Weather prediction systems and climate models are based around a dynamical core, the model component that solves the differential equations describing atmospheric fluid dynamics on scales resolved by the model's grid. A number of important quantities, including entropy, water content, and potential vorticity, are approximately conserved following the fluid motion, i.e. in a Lagrangian sense, throughout much of the atmosphere. There are theoretical arguments, and some supporting evidence from numerical models, that the use of numerical methods that better capture Lagrangian conservation properties in a dynamical core can lead to improved accuracy for long-range forecasts and improved physical realism for climate simulation.
One potential way to improve Lagrangian conservation properties is through the use of a Lagrangian or quasi-Lagrangian vertical coordinate, i.e. one that moves with the fluid. Most obviously, for a Lagrangian vertical coordinate the vertical advection terms are eliminated from the governing equations, so numerical errors in vertical advection terms are eliminated from a Lagrangian vertical coordinate (LVC) numerical model. Representation of fields with a long Lagrangian time-scale, such as water vapour in the lower stratosphere (e.g. Gregory and West, 2002;Hardiman et al., 2015) is influenced by such errors, and these errors can be reduced by the use of a quasi-Lagrangian vertical coordinate (Mahowald et al., 2002;Rasch et al., 2006).
Horizontal advection errors can also be reduced. The stratification of the atmosphere is determined by the vertical profile of entropy and it is mostly stable, except in regions of strong vertical convection. Cross-isentropic mixing is therefore much weaker than the quasi-horizontal mixing, which in turn leads to isosurfaces of long-lived constituents becoming almost aligned with isentropes. Consequently, such long-lived constituents appear smoother on isentropic surfaces than on other quasi-horizontal surfaces, which should help to reduce the horizontal advection errors in an LVC model where model levels are aligned with isentropes.
The idea of using an entropy-based or hybrid entropy-based vertical coordinate was explored in various models over the years (e.g. Branković, 1981;Hsu and Arakawa, 1990;Webster et al., 1999;Konor and Arakawa, 2000;Johnson et al., 2000;Randall et al., 2000;Schaack et al., 2004). Despite some of the anticipated benefits realized in such models, the perceived limitations did not encourage their operational use. Some challenges of using a pure isentropic coordinate (based on the potential temperature, ) include difficulties in defining and handling the bottom boundary, reduced vertical resolution in near-neutral stratification, and inextricable coupling of resolved dynamics with parametrized physical processes (with playing two different roles simultaneously). All of the above can lead to reduced robustness in such models. Nevertheless, using an entropy-based vertical coordinate model may help to reveal the role of entropy and unavailable energy budgets of weather and climate models. Johnson (1997) showed that entropy budget errors might be responsible for "cold pole" problems in climate models, and there is some some evidence (e.g. Webster et al., 1999) that the use of an entropy-based vertical coordinate can reduce them.
When it comes to energy, the discussion usually focuses on the need for adequate conservation of total energy in models without further distinction between its available (potential plus kinetic) and unavailable components, which are independently conserved under adiabatic conditions. The available energy budget is crucial, for example, for determining the strength of mid-latitude eddies (e.g. Barry et al., 2002). Since only a small fraction (about 0.2%) of the total potential energy is available for conversion into kinetic energy through adiabatic processes (e.g. Peixoto and Oort, 1992), any spurious conversions between unavailable and available energy could seriously impact the simulation of these eddies. The distribution of mass as a function of specific entropy, which determines both the total entropy and the unavailable energy, should be robustly conserved in adiabatic flow (Thuburn, 2008). Spurious cross-isentrope mass fluxes are nevertheless common in typical dynamical cores (Woollings and Thuburn, 2006) so a properly formulated LVC model with coordinate surfaces close to isentropes should reduce these errors. Lin (2004) describes an LVC dynamical core that is one of the optional dynamical cores of the National Center for Atmospheric Research (NCAR) Community Atmosphere Model version 3 (CAM3). Using the Lin dynamical core with slightly retuned but otherwise identical physical parametrizations in comparison with other CAM3 dynamical cores led to improved transport characteristics such as sharper tracer gradients, reduced vertical mixing and cross-tropopause transport, and better consistency between thermodynamic and advected quantities (Rasch et al., 2006). Also, Bala et al. (2008) reported improvements in dynamical fields, including surface wind stress, with resulting improvements for example in the simulation of Arctic sea ice. A version of the Lin dynamical core with a cubed sphere horizontal grid is used in the Geophysical Fluid Dynamics Laboratory (GFDL) climate model (Donner et al., 2011).
While the above results are encouraging, the Lin (2004) dynamical core makes the hydrostatic approximation. In order to be applicable to all scales from a few km to global, a dynamical core must be based on the non-hydrostatic compressible Euler equations. Hybrid-isentropic vertical coordinates have been used in limited-area non-hydrostatic dynamical cores by Zängl (2007) and by Toy and Randall (2009);Toy (2011). A non-hydrostatic version of the Lin global dynamical core has been developed and is currently being implemented in the US National Oceanic and Atmospheric Administration (NOAA) Environmental Modelling System (NEMS; Lin et al., 2016) for operational weather prediction.
In this paper we describe the implementation of a Lagrangian vertical coordinate in a global dynamical core using almost identical numerical methods to ENDGame ("Even Newer Dynamics for General atmospheric modelling of the environment"), the fully compressible non-hydrostatic dynamical core of the Met Office's Unified Model (UM) for weather and climate prediction (Wood et al., 2014). Section 2 outlines the current operational (height-based) version of ENDGame, followed by a more detailed description of the LVC implementation. The derivation of the Helmholtz problem for the semi-implicit semi-Lagrangian discretization is emphasized, along with the form chosen for the pressure gradient term and its effect on wave propagation.
An important issue with a Lagrangian vertical coordinate is that Lagrangian surfaces can bend and fold, on time-scales as short as 1-2 h in some situations (Kent, 2006), after which they become unusable as coordinate surfaces. One way to deal with this issue is to depart from a strictly Lagrangian vertical coordinate and use an Arbitrary Lagrangian Eulerian (ALE) approach (Toy and Randall, 2009;Toy, 2011). An alternative, used in the Lin dynamical core and also here, is to reset the locations of the model levels at regular intervals and remap the model fields to the new locations. This raises two further questions. First, to what locations should the model levels be reset? For brevity we call these locations "target levels". The Lin dynamical core resets its levels at every time step using target levels defined by a terrain-following pressure-based coordinate. We test a similar option but defined by a terrain-following height-based coordinate. However, it is possible that such "strong" resetting might undo some of the benefits of the Lagrangian vertical coordinate. Therefore, we also test quasi-Lagrangian target levels, defined in such a way as to follow isentropes on short horizontal scales and specified height levels on large horizontal scales and near boundaries. The aim is keep the model levels nearly Lagrangian where possible, and so reduce any errors introduced by remapping, while avoiding the problems of folding and coarse vertical resolution in weak stratification. The choice of target levels is discussed in section 3. The second question is: how should model fields be remapped to the target levels? Remapping can introduce errors in conservation or balance. We therefore formulate different remapping strategies that maintain particular conservation or balance properties (section 4).
Some sample results from a standard baroclinic instability test case (Ullrich et al., 2014) are presented in section 5, and some effects of the different choices of target levels and remapping strategy are noted. The LVC version of ENDGame was found to be rather less robust than the height-based version; section 6 discusses some of the problems found and the modifications made to the formulation to mitigate them. Conclusions are drawn in section 7. A detailed comparison of the Lagrangian conservation properties of the HB and LVC versions is presented in Kavčič and Thuburn (2018) (hereafter Part II). More details of the LVC implementation can be found in Appendices A-C.

LVC FORMULATION OF ENDGAME
This section presents the LVC ENDGame model set-up from the continuous governing equations to their discretization, solution procedure and numerical dispersion properties. The LVC formulation is based closely on the height-based (HB hereafter) formulation of Wood et al. (2014). As part of his collaboration on the development of ENDGame, the second author implemented a stand-alone version of ENDGame. It is simplified in that it is restricted to spherical deep-atmosphere geometry, the code is not parallel, and it stands outside the UM infrastructure; however, the numerical methods used are almost identical to those in the version used operationally in the UM. This stand-alone version of ENDGame was the starting point for the LVC implementation discussed here, and was also used to compare HB and LVC versions in section 5 below and in Part II. For the purpose of comparison with the LVC formulation, a brief overview of the HB ENDGame is given below.
• Prognostic variables are F = (u, v, w, , ), where u = (u, v, w) are the velocity components expressed relative to an orthogonal basis, is potential temperature, and is density. The Exner pressure Π is diagnosed from and where needed. • Governing equations are expressed in terms of a general, non-orthogonal vertical coordinate . In the HB formulation, the height of -surfaces is time-independent. • The equations are discretized in a semi-implicit semi-Lagrangian (SISL) scheme with an off-centring parameter . The SISL scheme is time-centred when = 0.5, however a more implicit scheme is used in operational versions with = 0.55. A full three-dimensional calculation of departure points, based oṅ, is utilized in the model. • The spatial discretization is a second-order centred finite-difference scheme on an Arakawa C-grid in the horizontal and Charney-Phillips grid in the vertical (w and staggered vertically half a grid length relative to u, v, and Π). In spherical geometry, only v is stored at the Poles. • A nested iterative approach is used to solve the resulting coupled nonlinear system of prognostic equations, based on increments of the prognostic and diagnostic vari- Replacing the fixed HB vertical coordinate with the moving LVC results in modifications in several of the above outlined points. Where modifications are not required, the existing model formulation is kept to facilitate a clear comparison between the HB and LVC formulations of ENDGame.

Governing equations and vertical coordinate
The continuous governing equations of LVC ENDGame in vector form are Du Dt where are the velocity components in spherical polar coordinates, and are potential temperature and density as before, and is Earth's rotation vector. In contrast to the HB formulation, r is now an additional prognostic variable determining the height of the Lagrangian surfaces. Note that is an alternative formulation of the pressure gradient plus geopotential gradient (hereafter PG) term, which is expressed as c p Π − g in HB ENDGame. The equation of state is now expressed through the Montgomery potential M, defined as which plays a role analogous to the Exner pressure Π in the HB formulation. Here T is the temperature and Π ≡ (R ∕p 0 ) ∕(1− ) , where p 0 is a constant reference pressure and ≡ R∕c p , with the gas constant per unit mass R and the specific heat at constant pressure c p . The apparent gravitational vector g = (0, 0, −g) is calculated as g = Φ, where the geopotential for the deep atmosphere is where a is the mean Earth's radius. The governing equations are initially expressed in spherical polar coordinates ( , , r), then, following Staniforth and Wood (2003), are transformed to a general, non-orthogonal vertical coordinate using the transformations ( ) The general vertical velocitẏis noẇ the material derivative of model variables is given as and the three-dimensional velocity divergence is where v = (u, v, 0) is the horizontal part of the velocity vector, = r 3 ∕3a 2 , and the derivatives with respect to , and t in Equations 10-12 and below are taken at constant .
The Lagrangian vertical coordinate is defined by settinġ = 0, so that the vertical advection term disappears from the material derivative (11). Note that r on a given model level is then no longer independent of time, so, in the model code, geometrical factors involving r must be recalculated as needed.

Semi-implicit semi-Lagrangian discretization
Similarly to Wood et al. (2014), we apply a semi-implicit semi-Lagrangian (SISL) discretization to the continuous governing equations in section 2.1. A general form of the SISL discretization can be written as where F is a prognostic variable and G represents non-advective terms in the equation for F, is an off-centring parameter, and = 1 − . Subscript A denotes evaluation of the left-hand side of Equation 13 at the arrival point and time step n + 1; subscript D denotes evaluation of the right-hand side at the departure point and time step n. The arrival points are chosen to lie on the model's levels and regularly spaced longitude-latitude grid. However, the positions of departure points and related quantities need to be calculated. Since the vertical coordinate remains constant along trajectories, only the horizontal position of the departure point needs to be determined. It is then convenient to express the horizontal position in terms of the unit vertical vector k and write the horizontal part of the kinematic equation as as in Thuburn and White (2013). This leads to significant algorithmic simplifications as the solution of Equation 14 and evaluation of the right-hand side of Equation 13 now require only two-dimensional interpolation instead of three-dimensional as in HB ENDGame. A SISL discretization of Equations 1-4 is where The component form of = (Ψ u , Ψ v , Ψ w ) after the coordinate transformation (Equations 8-9) is where horizontal derivatives at constant r are given by Equation 8. In the momentum equation 15, the components of the right-hand side are first evaluated in the local coordinate system at the departure point (indicated by subscript D L below); they must then be transformed into components in the local coordinate system of the arrival point. This is done using the rotation matrix M for the deep-atmosphere spherical coordinate system defined in Appendix A.1. of Wood et al. (2014). Since the resulting equations retain the same form as in HB ENDGame, only the vertical momentum equation is given here for illustration: Here M 31 , M 32 and M 33 are the relevant rotation matrix elements, and all terms on the right-hand side are evaluated at the corresponding w departure points. As in HB ENDGame, V  (u, v, w) are velocity components, is potential temperature, is density, p is pressure, r is the height of Lagrangian surfaces and M the Montgomery potential. Thermodynamic variables and w are stored at cell-centred points (circles) and horizontal velocity components u and v on the cell edges (diamonds and triangles, respectively) in both LVC and HB ENDGame. r is stored at cell-centred points (circles) in LVC ENDGame is a non-hydrostatic switch (the fully non-hydrostatic option V = 1 is used for all the results discussed below). Also as in HB ENDGame, Δtw n+1 A is an optional Rayleigh damping term, initially added to the right-hand side of the vertical component of Equation 1 as − w and then discretized with a backward Euler time step. The specific form of the profile is defined in section 5 for the Ullrich et al. (2014) test case presented in this work. Separate off-centring parameters u and w have been introduced for the horizontal and vertical components of the momentum equation.
For the kinematic equation 14, the horizontal location of the departure point should be independent of the vertical velocity; therefore we use the quasi-shallow rotation matrix treatment from Thuburn and White (2013) with the 2 × 2 rotation matrix M defined in the Appendix A.2. of Wood et al. (2014). After calculating departure points for each dependent variable, two-dimensional cubic Lagrange interpolation is used to evaluate departure point values, with polar rows omitted for v. A monotonicity filter is applied to limit the values of n D (Equation 17).

Spatial discretization and boundary conditions
The LVC and HB ENDGame have the same horizontal Arakawa C-grid with horizontal coordinates ∈ [0, 2 ) and ∈ [− ∕2, ∕2] equidistantly spaced through the domain. The density and thermodynamic variables and Π are stored at cell-centred points (Figure 1b, circles), as are r and w. Horizontal velocity components u and v are stored on the cell edges, staggered half a grid length in the and directions respectively (Figure 1b, diamonds and triangles), with the Poles chosen to be v-points.
The vertical coordinate is defined to lie in the range ∈ [0, 1], and the initial heights of the Lagrangian surfaces r 0 are given by where D is the domain depth and r surf is the surface height. Vertical velocities w and heights of Lagrangian surfaces r are placed at bottom and top surfaces of model layers and all other variables are stored at the centres of model layers. This means that is no longer stored with w as in HB ENDGame (Figure 1a for LVC and Figure 1c for HB ENDGame). As in HB ENDGame, geometrical factors are evaluated where needed and the orography is only defined at surface w points (r surf , Equation 24) although it can be averaged horizontally to other points if required. Components of gradients and rotation matrices are also evaluated where they are needed. More details about the discrete operators can be found in Appendix A, following the methodology of Wood et al. (2014), Appendices B and C.
Spatially discretized scalar equations 16-18 stay in the same form. The right-hand sides of the componentwise discretized momentum equation 15 are with the corresponding discretizations of Ψ u , Ψ v and Ψ w as Here, superscripts and subscripts , , denote averaging or differentiation in the corresponding coordinate direction (Appendix A).
Some discussion of the influence of the vertical staggering and pressure gradient discretization on numerical normal modes is given in section 2.5.
Periodic boundary conditions are applied in the direction. The Poles are treated as described in Wood et al. (2014), i.e. calculating polar wind components from the values of u and w at the nearest u-latitude. Polar values of Ψ v , required to evaluate R n u , R n v and R n w near the Poles, are also calculated in this way. For the vertical boundary conditions, the lowest level follows the height of the terrain (r surf ) and there is no flux across either bottom or top boundary. Boundary values of w are calculated from the horizontal gradient of r surf and r top and constant extrapolation of u and v from the levels immediately above or below that surface, respectively. Constant extrapolation is also used for quantities stored at M, u and v points extrapolated to departure points that are closer to the top or bottom boundary.

Helmholtz equation and back-substitution
The discretized equations 15-18 (with from Equations 28-30) comprise a system of nonlinear equations for the unknown fields at time step n + 1. The system is solved using a quasi-Newton method in which only the stiffest terms in the Jacobian are retained, and the resulting linear system for the Newton update is solved only approximately at each Newton iteration. This approach leads to the following linear system for the Newton increments (indicated by ′ ): Here, P ′ = M ′ ∕ * is a modified Montgomery increment, and the various H coefficients are defined in Appendix B. The right-hand sides of Equations 31-36 are the residuals in the components of Equation 15 and in Equations 16-18 at the current iteration, that is, the amount by which those equations fail to be satisfied with the latest estimates for the unknowns.
Starred quantities (e.g. * , * ) indicate reference thermodynamic profiles about which the system has been linearized to obtain the approximate Jacobian. Following the operational HB ENDGame, they are taken to be the solutions from the previous time step; thus they vary in time as well as in all three space dimensions. Also following the operational ENDGame, the horizontal derivatives on the left-hand sides of Equations 31, 32 are approximated by neglecting the bent terms from the curvilinear coordinate transformation in order to reduce the computational stencil of the Helmholtz operator defined below.
The complexity of the linear system is reduced by elimi- where the vertical derivative operators D 1 and D 2 , the various Helmholtz coefficients, and the right-hand side termR rw are defined in Appendix B. As in Wood et al. (2014), two nested loops are used to obtain the model variables at the new time level. The outer loop solves the kinematic equation 14 using the latest available estimates for the velocity components and after that interpolates the required fields to the departure points. The right-hand sides of the (35) and (34) equations are also updated here to improve convergence.
Within the inner loop the Coriolis terms are updated and the residuals in the velocity (31)-(33) and r (36) equations are computed. Finally, the linear Helmholtz equation 37 is solved for P ′ , and the increments of all prognostic variables are found by back-substitution. The number of iterations in the inner and outer loops is chosen to manage the trade-off between the expense of the outer loop calculations and the diminishing returns of taking too many inner iterations. In the operational HB ENDGame, two outer and two inner iterations are used. In the LVC ENDGame, to minimize any effects on stability, four outer and one inner iterations are used for the results shown below.
For the solution of the Helmholtz equation 37, we use a column by column vertical line relaxation combined with a horizontal geometric multigrid method similar to that described by Buckeridge and Scheichl (2010). It uses a conditional semi-coarsening strategy to make it suitable for the longitude-latitude grid.

Treatment of pressure gradient term, and normal modes
The unusual form of the combined pressure gradient and geopotential gradient terms in Equation 19 deserves some discussion. First note that, for the continuous equations, so the form used in Equation 19 is equivalent to the more familiar form. However, these expressions are no longer equivalent when discretized, and the form used can affect the accuracy of wave propagation. Following Thuburn and Woollings (2005) and Thuburn (2006), for small perturbations to a resting, hydrostatic basic state, optimal wave dispersion is obtained by placing the pressure variable at the same levels as u and v and placing the buoyancy variable at the same levels as w. The vertical pressure gradient term must also be evaluated in such a way as to avoid averaging of the buoyancy perturbation. If perturbations on model levels can be neglected, then the LVC case becomes equivalent to the isentropic vertical coordinate case. In that case M plays the role of the pressure variable and z or Φ plays the role of the buoyancy variable, and the form (M∕ ) + Φ (ln ) of the pressure gradient term was shown by Thuburn (2006) to be optimal. The form used in Equation 19 is a slight variant on this one with the same dispersion properties, but giving rise to a Helmholtz problem that is very similar to that arising in the HB case.
More generally, the LVC case is not quite the same as the isentropic coordinate case. The LVC governing equations have an additional prognostic field compared with the height coordinate or isentropic coordinate governing equations, which leads to an extra branch of normal modes in the dispersion relation. These extra modes correspond to the freedom to describe a given atmospheric state with the coordinate levels in different positions. We might therefore call them "remapping modes" to distinguish them from the physical acoustic, inertio-gravity, and Rossby modes. For perturbations to a resting hydrostatic state, the remapping modes have zero frequency, since the physical state remains at rest and in hydrostatic balance.
In the discrete case, the remapping modes can be characterized as follows. (a) The horizontal pressure gradient force should vanish. Since the basic state must be horizontally uniform, this implies where * indicates basic state fields and ′ indicates perturbations.
(b) The vertical pressure gradient force should also vanish: (40) Substituting Equation 39 in Equation 40 and using the formula for the discrete derivative of a product leads to This is of the general form (q ′ ≈ r ′ r q * for any thermodynamic variable q) expected for remapping modes. However, because of the vertical averaging of ′ , the amplitude of r ′ will be suppressed for modes with fine vertical scale. In particular, there exists a mode with a vertical two-grid-length oscillation in ′ while r ′ ≡ 0. This mode may justifiably be called a computational mode, since it involves a non-zero pattern in ′ (and ′ ∕ * = − ′ ∕ * ), with level heights unchanged, that is invisible to the dynamics.
The numerical normal modes were computed directly from the vertically discrete linearized equations, following the method of Thuburn and Woollings (2005) and Thuburn (2006). Figure 2 shows the frequencies of the westward propagating acoustic, inertio-gravity, and Rossby modes, along with the corresponding frequencies for the normal modes of the continuous equations. The figure is virtually identical to that for other discretizations with optimal wave propagation, confirming the optimal wave dispersion properties of the scheme based on Equation 19 for these physical modes. There are three other branches to the numerical dispersion relation not shown on the figure, corresponding to eastward propagating acoustic and inertio-gravity modes, and, as predicted, a set of zero-frequency remapping modes. There is some freedom in the vertical placing of in the LVC formulation as it does not influence the dispersion properties of the physical modes the way it does in the HB formulation. Placing together with the predicted mass variable has the advantage of allowing , and specific humidity to be collocated when water is predicted. This, in turn, would facilitate conservation of water substance at the same time as an accurate representation of moist thermodynamics.

TARGET LEVELS
The LVC formulation of ENDGame requires frequent resetting of coordinate surfaces to ensure that they do not fold over.
When the model levels are reset, some choice must be made for the new values of r on model levels; we refer to these as "target levels". The simplest choice is to define the target levels to equal the initial values r 0 (Equation 24), and this is one option that we test (labelled R0). However, remapping to r 0 at every time step might result in as much numerical diffusion as the use of a vertical advection scheme in a height-based coordinate, and so be just as detrimental to Lagrangian conservation properties. This motivates us to investigate also a quasi-Lagrangian (QL) choice of target levels; the idea is that the target levels will remain close to Lagrangian levels, so that remapping will introduce less numerical diffusion. Specifically, the idea is to choose the target levels to minimize a certain cost function: where indicates the horizontal gradient along levels, and the integral is a global horizontal integral over an level. The second term in the integral penalizes variations in along model levels and so will cause the target levels to tend to follow isentropes on small horizontal scales. The first term in the integral penalizes departures of the target levels from the initial level locations and will cause the target levels to follow r 0 on large horizontal scales. The coefficients 2 ( ) are a set of tunable parameters, one per model level, that determine what is meant by "small" and "large" horizontal scales here; see below. Note that this cost function does not penalize departures from any specific isentropes, only horizontal gradients of . Thus it does not tend to produce coarse vertical resolution in regions of weak stratification such as the tropical upper troposphere, a known weakness of isentropic coordinate models. The cost function J is minimized if the variation J vanishes for arbitrary small changes r in level location: Thus J will vanish for arbitrary r provided In practice we have found it helpful to include some additional terms in Equation 45 to ensure that the level locations remain smooth on small horizontal and vertical scales: where, at level k + 1∕2, D r (r) is given by where D is the domain depth.
For the results presented below, the parameters used to define the quasi-Lagrangian target levels are set as follows. A base length-scale over which the target levels should follow isentropes is set to l b = 2 × 10 6 m. Near the top and bottom boundaries, target levels should be close to r 0 to avoid problems with levels crossing the boundaries or each other, so a modified length-scale is defined by where z 0 is the initial height of the level above the lower boundary and z l = 5000 m is a transition height scale over which the target levels become more quasi-isentropic. The coefficient 2 is then given by where † r is an estimate for the vertical gradient of based on an isothermal atmosphere at 300 K. The remaining parameters are l s = 2 × 10 5 m and h s = 10 3 m.
Note that, whereas Equation 45 is a set of two-dimensional elliptic problems for r, one per model level, the inclusion of the D r term introduces some vertical coupling and makes Equation 46 a single three-dimensional problem. Also, because of the 2 term, the problem is nonlinear. In fact, the nonlinearity is not too strong, and Equation 46 can be solved using a quasi-Newton method similar to that used for the semi-implicit time integration. After l iterations, Equation 46 will not be satisfied but will have some residual given by where the superscript (l) indicates the estimate after l iterations. An improved value for r is then given by where r ′ satisfies The left-hand side of Equation 52 is an approximate linearization of the left-hand side of Equation 46. It is solved using a multigrid method similar to that used to solve the Helmholtz problem associated with the semi-implicit time integration. Because the target levels are recomputed at each time step, they are never very far from the actual model levels.
Thus it is sufficient to use a single multigrid V-cycle to solve Equation 52 and a single quasi-Newton iteration to solve Equation 46. The first guess is given by the actual model level locations.
Finally, it is worth noting that the parametrization packages of some models (though not the UM) require specific locations for the model levels. The QL option for target levels would not be suitable for use with such packages.

REMAPPING
After selecting target levels, model fields must be transferred to the new coordinate surfaces. This gives rise to a one-dimensional remapping problem in each column of the dynamical core. 1 It is possible to design this remapping to have certain desirable properties, such as conservation of mass, energy, or entropy, or preservation of hydrostatic balance; however, not all of these properties can be obtained simultaneously. We therefore compare a number of options in section 5.

Global diagnostics
This section defines some of the integral diagnostics that are conserved by certain remapping options, or are used as global diagnostics for the performance of the model and comparison with HB ENDGame. The total mass  is calculated as where dA is the area element on the unit sphere. The total energy  is where E denotes the energy per unit mass divided into internal plus potential, E P , and kinetic, E K , parts. Finally, the total entropy  is The total imbalance is a global measure of the departure from hydrostatic balance. It is calculated as 1 The term "remapping" is sometimes used to mean the conservative transfer of data from one grid to another. Here we use "conservative remapping" to explicitly identify the conservative cases.
To compute these diagnostics in the model, the global integrals are approximated as sums over grid cells.

Remapping of mass and thermodynamic variables
All four remapping strategies described below conserve mass by conservative remapping of , as described in section 4.4. The remaining freedom in the remapping of may be used to conserve energy or entropy or to preserve hydrostatic balance.

Mass-energy remapping (-)
In this option, the internal plus potential energy density E P is diagnosed from the model predicted variables. Then column values of mass density and energy density E P , considered as functions of the vertical coordinate , are conservatively remapped to the target levels. Finally, the new values of are inferred from the new values of , E P and r. In practice we did not include the kinetic part of the total energy, as the potential energy is always dominant, and this greatly simplifies the remapping problem.

Mass-entropy remapping (-)
In this option, the specific entropy c p ln is first diagnosed from . Then column values of mass density and entropy density c p ln , considered as functions of the vertical coordinate , are conservatively remapped to the target levels. Finally, the new values of are inferred from the new values of and c p ln .

Hydrostatic remapping (-H )
As in the first two options, is conservatively remapped. The new values of are computed so as to avoid, as far as possible, the introduction of any vertical imbalance through the remapping process. A hydrostatically balanced potential temperature profile, H , is calculated from and r for old and new model levels; to remove the ambiguity introduced by the computational mode, the hydrostatic profile is chosen to have the same internal plus potential energy as the original profile. For the original levels, the non-hydrostatic contribution to the potential temperature NH = − H is diagnosed and then interpolated to the target levels using cubic Lagrange interpolation. Finally the new value of is reconstructed using the values of NH and H on the new levels: = H + NH .

Interpolation of (-I )
Again, as in the first two options, is conservatively remapped. However, a simple cubic Lagrange interpolation is used for (linear near the boundaries) with a limiter to ensure that the interpolated value is bounded by the values above and below.

Remapping of velocity
We also explored different strategies for transfer of the velocity components to target levels. Our initial tests used cubic Lagrange interpolation (without and with a limiter) to interpolate u, v and w. However, in some situations this scheme could significantly alter the three-dimensional divergence leading to large imbalances in the mass equation.
The results presented below all use an alternative scheme designed to prevent this problem. The horizontal velocities u and v are conservatively remapped so as to preserve the column integral of horizontal divergence. The three-dimensional divergence, given by Equation 12, is also conservatively remapped, and then the divergence calculation is inverted to diagnose the new w. The same remapping methods (section 4.4) used for transfer of thermodynamic variables are also used here.

Conservative remapping methods
The conservative remapping strategies for the column values of model variables described above use methods based on one-dimensional piecewise polynomial (White and Adcroft, 2008) and spline (Zerroukat et al., 2006) reconstructions of cell averages. Piecewise polynomial methods presented in White and Adcroft (2008) use fourth-degree (Piecewise Quartic Method or PQM) and second-degree (Piecewise Parabolic Method or PPM) polynomials, with a limiter for PQM which ensures its monotonicity. Estimates of cell edge values and slopes are obtained from explicit or implicit schemes of degrees from two to five, with the implicit schemes shown to be more accurate but also computationally more expensive. Boundary values can be calculated from one-sided polynomials of the same order or gradually decreasing the order down to a constant extrapolation at the boundary. Neale et al. (2012) use PPM for vertical remapping in the Community Earth System (CESM) model (their section 3.1.6. gives more details).
The Parabolic Spline Method (PSM) from Zerroukat et al. (2006) constructs a piecewise parabolic fit to the cell average data that gives continuous derivatives at cell edges. Although PSM formally has a similar order of accuracy as PPM, it is more accurate in practice for problems with slow decay of energy spectrum such as atmospheric flows.
After testing the model with all the above listed methods, we found that they influenced the model runs much less than the specific remapping strategy within which they were used (sections 4.2 and 4.3). The results presented in section 5 all use PSM as specified above.

Frequency of remapping
All of the results presented below and in Part II use remapping at every step. We have also tested hourly remapping, but found little sensitivity to this choice.

Test case set-up
Here we present the behaviour of LVC ENDGame tested on the deep-atmosphere baroclinic instability test of Ullrich et al. (2014). The basic state comprises a zonally symmetric, steady and balanced but unstable eastward jet in each hemisphere. We use an iterative scheme to improve the discrete balance of the basic state. The baroclinic instability is triggered by a horizontally non-divergent perturbation to the zonal and meridional wind components localized to the lower atmosphere; this reduces the unbalanced motions resulting from the trigger perturbation compared to the earlier, shallow atmosphere, baroclinic instability test case of Jablonowski and Williamson (2006). Following Ullrich et al. (2014), a vertically stretched grid is used to provide finer vertical resolution near the bottom boundary: A time step Δt = 1200 s is used. The departure points from the kinematic equation 14 are evaluated by a purely centred in time scheme ( x = x = 0.5). A slight off-centring is used for the other variables ( u = w = = 0.51). As mentioned in section 2.4, the iterative solution procedure uses n out = 4 outer and n in = 1 inner iterations to minimize effects on stability. With the choice of two outer and two inner iterations (as used in the operational HB ENDGame) the LVC ENDGame runs failed in day 1 for all remapping strategies.
In contrast to Wood et al. (2014) (but similarly to the operational set-up of ENDGame), the damping coefficient in the w-momentum equation 23 is increased near the Poles as well as near the model top: where D = 1∕ D and P = 1∕ P define values at the domain top and Poles, respectively, and the values of the parameters are D = P = 7200 s and h D = 3000 m. As mentioned in section 2.2, a monotonicity filter is applied for the interpolation of n D . Also a strong zonal filtering is applied at the latitude row nearest the Pole to the zonal wind and scalar  fields, retaining only the zonal mean and zonal wavenumber 1 contributions. The model is run for N T = 1080 time steps, that is 15 days. Remapping is done column by column at the end of every time step. The effect of different choices of target levels (section 3) and remapping strategies (section 4) on model results, including global conservation and balance, is compared. In both cases the -I remapping strategy is used, that is, conservative PSM remapping of mass and cubic Lagrange interpolation of potential temperature.

Model behaviour, remapping and conservation properties
The overall agreement between HB and LVC ENDGame is good for R0 target levels at both day 8 and day 10. However, the QL target levels seem to delay the time development of the baroclinic wave. Further investigation (not shown here) revealed that this behaviour is consistent for all four remapping strategies, and that the time delay for QL target levels is about half a day (e.g. the QL runs at 8.5 days show similar structure to HB and R0 fields at day 8). For the initial state, the QL target levels differ significantly from the initial model levels in some locations. Therefore, for the QL runs, there is a significant jump in level heights at the end of the first time step when remapping is first done. A plausible hypothesis might be that this first remapping results in a damping of the trigger perturbation leading to the observed delay in the development of the instability. However, when we tested this hypothesis by initializing the model on QL levels before adding the trigger perturbation, the same delay in development was observed, implying that the hypothesis is not correct. Thus, the explanation for the delay in development in the QL version remains unknown. Other effects of this jump are noted below.
Locations of model levels and cross-sections of potential temperature for the QL -I case at day 10 are shown in Figure 6. The figure confirms that using Equation 46 to define target levels does indeed give suitable quasi-Lagrangian levels. In areas of stable stratification away from the lower boundary and on short horizontal scales, the target levels deform with the isentropes, but without noise or folding. On large horizontal scales and near the boundaries, the target levels remain close to their initial positions, and avoid producing regions of very coarse vertical resolution.
In the LVC ENDGame, the choice of target levels and remapping strategy had only a small effect on the fields shown in Figures 3-5, other than the delay in evolution resulting from the use of QL target levels. However, the stability of the model is very sensitive to these choices. We have been able to run the model for 15 days with a reasonable time step of Δt = 1200 s in most R0 remapping configurations, while for  QL target levels only the -I configuration completed the run (Table 1). The runs with the -H remapping strategy are particularly unstable for both R0 and QL target levels, ending near day 11 and day 13, respectively (Table 1, bold text). The robustness of the LVC formulation is discussed further in section 6. The effect of different target levels and remapping strategies on global conservation and balance is demonstrated by Figures 7 and 8, which present the relative changes in total mass  (Figure 7a,b), total energy  (Figure 7c,d), total entropy  (Figure 8a,b) and imbalance  (Figure 8c,d) for all combinations of target levels and remapping strategies.
As long as the model remains stable, the global conservation properties of LVC ENDGame are comparable to those of HB ENDGame. The LVC formulation with the QL remapping does nearly as well as HB ENDGame (Figures 7 and 8,  crosses and black circles, respectively), while the LVC R0 remapping formulation preserves global diagnostic quantities even better than HB ENDGame (relative changes of ∼ 0.5 × 10 −4 for LVC R0 compared to ∼ 1.5-2 × 10 −4 for HB ENDGame). As noted above, for the QL choice of target levels, there is a significant jump in level heights at the end of the first time step. For the -, -, remapping strategies, one effect of this jump is the introduction of a significant imbalance (Figure 8c,d). This imbalance is much greater than any imbalance that arises during the later evolution, though it fades away over the first 2 days of the integration. The -H remapping strategy, on the other hand, is designed to avoid the introduction of spurious imbalance by the remapping. Figure 8d shows that the imbalance in the QL -H run remains similar in magnitude to that in the the HB run and LVC R0 runs over the first 2 days, confirming that it behaves as intended. This interpretation was confirmed by running further integrations that switched from R0 target levels to QL target levels after a few hours; imbalance was generated at the time of the switch, except with -H remapping.
Another interesting feature of the runs with QL target levels is the occurrence of small-amplitude oscillations in the total mass, accompanied by synchronous oscillations in total energy and entropy (Figures 7b,d and 8b,d). Such oscillations are not present for R0 remapping strategies. Since all remapping strategies are mass conserving, the mass change must occur via the dynamical integration of the mass equation 16 rather than the remapping. However, it is not clear what aspect of the dynamics is different for QL target levels compared with R0 target levels. Note these oscillations in ,  and  are apparently not directly linked to the large initial imbalance discussed above, since the QL -H configuration does not suffer such large initial imbalance.
In Part II a variety of diagnostics are used to make a detailed comparison of the Lagrangian conservation properties of the the HB and LVC model versions.

Computational cost
Despite carrying an additional prognostic field and the need to recompute geometrical factors as r evolves, and despite the need for remapping, the LVC formulation is considerably cheaper than the height-based formulation. Table 2 shows run times on a single processor for various model configurations. The difference in cost between 4 × 1 and 2 × 2 Newton iteration schemes confirms that the departure point calculations and semi-Lagrangian interpolation are dominant costs in both LVC and HB versions. The reduction of these calculations from three dimensions to two dimensions explains why the LVC version is cheaper than the HB version. The table also shows that the additional cost of the remapping calculations is rather small. Relative change in (a, b) total entropy  and (c, d) imbalance . The initial value of total entropy for HB ENDGame is  0,HB = 3.0157 × 10 22 J K −1 kg −1 and for LVC ENDGame is  0,LVC = 3.0167 × 10 22 J K −1 kg −1 for all remapping strategies with ratio  0,LVC ∕ 0,HB = 1.0003; the respective initial value of total imbalance for both HB and LVC ENDGame is  0 = 0 m s −2 . Other details are as in Figure 7

ROBUSTNESS
The LVC formulation of ENDGame has been found to be much more susceptible to instability than the HB formulation.
Because model levels can move, the model has additional degrees of freedom and, it might seem intuitively, more scope for things to go wrong. This idea seems to be born out in practice. This section summarizes some of the stability issues The horizontal resolution is 2 p ×2 p−1 grid points with 32 vertical levels. The code was run for 100 time steps on a single processor. n out and n in indicate the number of outer and inner Newton iterations. encountered and the adjustments made to the formulation to mitigate them. We believe it is valuable to document these issues for the benefit of those attempting similar developments in the future.
It is important to note that the Ullrich et al. (2014) test case is very challenging for a dynamical core. In the absence of physical damping processes such as a realistic boundary layer, strong fronts form, initially at the surface, and collapse to the grid scale, resulting in very strong horizontal and vertical gradients of velocity and temperature, with strong vertical velocities and overturning circulations. Nevertheless, the HB ENDGame is robust on this test case and it was hoped that a similar robustness could be obtained with the LVC ENDGame.
Note also that the time step is quite large for a non-hydrostatic dynamical core, and that very little numerical damping is used beyond that inherent in the advection and remapping schemes. We wish to avoid, as far as possible, the pitfalls of adding extra numerical damping to control an unstable scheme.
1. For the Coriolis terms we initially used the same scheme as in HB ENDGame. It is formulated in terms of mass fluxes to ensure energy conservation (Wood et al., 2014, based on the shallow-water scheme recommended by Thuburn and Staniforth, 2004). However, the LVC formulation was found to be very sensitive to the influence of layer thickness on the mass fluxes and suffered from an instability with short vertical and meridional scale in regions of strong zonal wind. In fact the Thuburn and Staniforth (2004) shallow-water formulation has recently been shown to be unstable when the mean depth is very small (Peixoto et al., 2018); it seems likely that this result is related to the problem found in LVC ENDGame. In order to avoid instabilities we therefore use a simple non-energy-conserving scheme as written in Equations 28-30. As shown in Figure 7, the model nevertheless retains good energy-conserving properties. 2. Instead of Equation 2, we initially used a more "natural" formulation of the continuity equation that does not involve w: where is the density with respect to the vertical coordinate , defined as = . (60) The modified divergence L ⋅ is defined by where H ⋅ is the contribution to the three-dimensional divergence operator that involves only horizontal derivatives of velocity. However, this version of the model was susceptible to an instability that occurred in strong vertical shear. It appears to arise through an inconsistency between the horizontal advection of r, which uses vertically averaged u and v, and the horizontal advection of , which uses u and v that have not been vertically averaged. Consequently, we changed to the form of the continuity equation. 3. The LVC ENDGame appears to be more susceptible to slow solver convergence and instability at the Poles than the HB version. One specific problem is that the calculation of the three-dimensional divergence becomes very sensitive to the level locations r when the horizontal advective Courant number is large. To mitigate these problems, some w-damping and polar filtering is used, as discussed in section 5.1. Such measures are not needed in HB ENDGame to run the baroclinic instability test case, though some w-damping, increasing at the Poles, is used operationally at high resolution. 4. When strong vertical temperature gradients form at fronts, most of the remapping strategies tested can lead to overshoots in the remapped potential temperature; the model then responds by convecting vigorously on resolved scales before failing. Only the simple cubic interpolation of (-I ) with a limiter to prevent overshooting consistently avoided this problem. 5. Our initial formulation used cubic Lagrange interpolation for remapping of u and v. When strong vertical shear formed at fronts, the remapping of u and v could lead to overshoots and instability. Including a limiter in the interpolation of u and v eliminated the overshoots, but revealed a further problem. The three-dimensional divergence field is very sensitive to small changes in the velocity components, and could change sign locally as a result of the remapping, again leading to instability. To avoid this problem, the divergence-conserving scheme described in section 4.3 is used instead.
Despite all these adjustments to the formulation, and with a time step of Δt = 1200 s, the LVC formulation is still relatively fragile, as seen in section 5. At the same resolution the HB ENDGame is robust on this test case with Δt = 1800 s and without w-damping or polar filtering.

SUMMARY
A Lagrangian vertical coordinate formulation of the ENDGame dynamical core has been presented, solving the deep-atmosphere, non-hydrostatic, compressible Euler equations on the sphere. The formulation is as similar as possible to the original height-based formulation of Wood et al. (2014), allowing a clean comparison of costs and performance. As in the height-based formulation, a semi-implicit semi-Lagrangian time integration scheme is used with a centred-difference spatial discretization; the resulting nonlinear system for the unknown fields is solved iteratively via a Helmholtz problem for increments to a pressure-like variable at each iteration. The key differences are that (a) model level height is now a predicted field (section 2.1), (b) model level heights must be periodically reset and model fields remapped to these new levels (sections 3 and 4), and (c) a modified form of the pressure gradient term is needed to ensure optimal wave propagation of physical normal modes (section 2.5).
The LVC formulation is considerably cheaper than the HB formulation because the departure point calculations and semi-Lagrangian interpolations are reduced from three dimensions to two dimensions.
Two options have been investigated for the target levels for resetting model level heights: (a) the initial model level heights (R0), and (b) quasi-Lagrangian (QL) levels that follow isentropes on small horizontal scales and initial model level heights on large horizontal scales and near boundaries. Also, several strategies have been investigated for remapping model fields to target levels, including options to conserve internal plus potential energy or entropy or to preserve hydrostatic balance.
Results from the baroclinic instability test case of Ullrich et al. (2014) showed that, while the model remained stable, the LVC and HB versions gave very similar results, though with a time lag of about half a day when the QL option for target levels is used. The results confirmed that, with the QL option, model levels do indeed behave in a quasi-Lagrangian manner. They also confirmed that the different remapping strategies display the conservation and balance properties they were designed to have. However, the most significant result from the baroclinic instability test case is that the LVC ENDGame is much less robust than the HB version, and that its stability is strongly sensitive to both the choice of target levels and the choice of remapping strategy.
Despite the poor stability of the model, and the use of a simple non-conservative Coriolis discretization, the LVC formulation shows global conservation properties as good as those of HB ENDGame for QL or even superior to it for R0 remapping strategies (section 5.2). Part II shows that there are small but measurable improvements in some diagnostics of Lagrangian conservation with the LVC-QL version, though also some degradation in other diagnostics. These results, together with the reduced computational cost, encourage us to continue efforts to understand and improve the stability of the LVC ENDGame formulation.