ExoCubed: A Riemann-solver-based Cubed-sphere Dynamic Core for Planetary Atmospheres

The computational fluid dynamics on a sphere is relevant to global simulations of geophysical fluid dynamics. Using the conventional spherical–polar (or lat–lon) grid results in a singularity at the poles, with orders-of-magnitude-smaller cell sizes at the poles in comparison to the equator. To address this problem, we developed a general circulation model (dynamic core) with a gnomonic equiangular cubed-sphere configuration. This model is developed based on the Simulating Nonhydrostatic Atmospheres on Planets model, using a finite-volume numerical scheme with a Riemann-solver-based dynamic core and the vertical implicit correction scheme. This change of the horizontal configuration gives a 20-time acceleration of global simulations compared to the lat–lon grid with a similar number of cells at medium resolution. We presented standard tests ranging from 2D shallow-water models to 3D general circulation tests, including Earth-like planets and shallow hot Jupiters, to validate the accuracy of the model. The method described in this article is generic to transform any existing finite-volume hydrodynamic model in the Cartesian geometry to the spherical geometry.


INTRODUCTION
Circulation alters the thermal structure of the atmosphere and the distribution of chemical tracers within.However, the complex nature of solving partial differential equations (PDE) on a sphere poses a significant challenge to renovating the design of an existing GCM and staying abreast of the progress in numerical techniques.
The diverse climates on exoplanets present both challenges and opportunities for simulating planetary atmospheres.In particular, hot Jupiters exhibit enormous contrasts in temperature between their day and night sides, a phenomenon not observed within our solar system (Showman & Guillot 2002;Bell & Cowan 2018;Komacek & Showman 2016).The large pressure difference may drive supersonic jets and the resulting dissipation could be a significant energy source for the planetary interior (Fromang et al. 2016).However, Earth-centric GCMs, as well as planetary GCMs that are inherited from Earth GCMs, do not handle the formation and the instability of supersonic jets well: the maximum wind speed on Earth is much less than the sound speed (Nasr et al. 2022).Moreover, the James Webb Space Telescope (JWST) will reveal the atmospheric composition of mini-Neptunes, which are potentially hydrogen and water-rich planets.Since the mass mixing ratio of water vapor in the atmosphere may be comparable to the hydrogen envelope, the interplay of thermodynamics and atmospheric circulation could be extremely complex.A carefully constructed "moist" GCM is needed to decipher the future transit observations of mini-Neptunes (Constantinou & Madhusudhan 2022), suggested by CAMEMBERT, a GCM intercomparison project for mini-neptunes (Christie et al. 2022).
Despite the numerous GCMs developed for studying Earth's atmosphere, with at least 167 different models identified in a recent study of GCM code development (Kuma et al. 2023), they all come down to about 12 major families, which stem from weather forecasting models funded by meteorological agencies (Kuma et al. 2023).Due to the sharing of development practices, subroutines, numerical treatments, and parameterization schemes, the biases of the simulations (e.g.climate projection) are statistically correlated (Kuma et al. 2023).As a result, it is crucially important to have a variety of GCMs available to the community, each constructed with different numerical methods and physical approximations, to allow for a diversity of approaches to solving the unknown.
However, developing a GCM from scratch is a time-consuming task that requires knowledge well beyond a basic understanding of numerical methods and physics.The major difficulty is to deal with the spherical geometry, which is almost unique to geophysical fluid dynamics.The spectral method, based on the spherical harmonics, is a popular choice for solving partial differential equations on a sphere that overcomes many of the difficulties associated with the spherical geometry.Yet, the scalability of the spectral method is poor as global communication between computer processors is required.A grid-based method can achieve high scalability on parallel execution but gridding the sphere introduces many complexities that are not present in the Cartesian geometry.
One way of gridding the sphere is to use unstructured grids, which are both flexible and scalable.Unstructured grids are widely used in computational fluid dynamics (CFD), especially for simulating flows around complex geometries.Yet, the performance of unstructured grids is generally inferior to structured grids, especially for problems with simple geometry, due to the overhead of computing the connectivity between the grid points.Constructing a GCM on unstructured grids has been attempted since the 1960s (Sadourny et al. 1968).However, it did not gain popularity until the environmental model OMEGA was developed in 2000 (Bacon et al. 2000), and then THOR GCM (Mendonça et al. 2016;Deitrick et al. 2020) was developed for planetary atmospheres in 2016.The THOR GCM uses the finite-volume method on icosahedral grids and leverages the GPU computing power to achieve high performance.
Structured grids can achieve high performance on simple geometry but there is a trilemma for a GCM constructed on a structured grid using explicit time-stepping.That is, we have to sacrifice one of the following three properties: (1) devoid of numerical singularity, (2) quasi-uniform grid spacing, and (3) orthogonal coordinate lines.This is illustrated with Figure 1.Here, we regard the numerical singularity as a point where a numerical scheme ceases to be well-behaved.Special treatment is needed to deal with the point at the singularity.The traditional lat-lon grid has uniform grid spacing and orthogonal coordinate lines, but it has a singularity at each pole.Due to the converging meridians at the poles, polar filters are needed to stabilize the numerical solution.For example, a polar filter is applied to latitudes higher than 45 degrees in the EPIC GCM for planetary atmospheres (Dowling et al. 1998) and the Community Atmosphere Model (CAM) for Earth's atmosphere (Neale 2010).One can avoid the numerical singularity by employing successively coarser grid resolution toward the pole (grid-stretching) at the expense of sacrificing the uniform grid spacing (e.g.Bindle et al. (2021)).
Mapping a spherical grid to a cube is another way of removing numerical singularity.The term "cubed-sphere" refers to a method of tessellating a sphere into six faces, akin to an expanded cube.The MITgcm (Adcroft et al. 2004a) uses a conformal mapping strategy to transform the grids on a cube to the sphere while preserving the orthogonality of the coordinate lines.The conformally mapped cubed-sphere grids have eight singular points, which are the corners of the cube, although they are not as severe as the poles in the lat-lon grid.
Cubed-sphere grids can also be constructed using gnomonic projection, which is non-conformal, but the projected grids are more uniform (Chen 2021;Ullrich et al. 2010).Moreover, the gnomonic cubed-sphere grids do not have numerical singularities, i.e. no special treatment is needed at the corners between the six faces.Yet, it is a non-trivial task to work with the PDE solver on a non-orthogonal geometry.The Flexible Modeling System (FMS) developed at the Geophysical Fluid Dynamics Laboratory (GFDL) uses the gnomonic cubed-sphere grids for its atmospheric GCM (Putman & Lin 2007;Mouallem et al. 2023).Both FMS-based GCMs and the MITgcm have been used to study exoplanet atmospheres (Kaspi & Showman 2015;Komacek et al. 2021).
In this article, we discuss the development of the dynamic core of a GCM that is designed for general planetary atmospheres that may have shocks, multiple condensing species, or deep atmospheres that are as thick as the planet itself.We target the modelers who have access to a hydrodynamic model in Cartesian geometry and we present a general method to convert such existing models to spherical geometry with minimum effort.In this way, a variety of hydrodynamic codes can be converted to a GCM, enriching the diversity of GCMs available to the community.We choose to use the finite volume method on gnomonic cubed-sphere grids as the changes to the existing Cartesiangeometry-solver are minimum among all other choices.
In section 2, we begin with a brief discussion of the finite-volume method and some of its variants.This section is important for understanding the subsequent choices made in non-orthogonal geometry.In section 3, we describe the construction of the gnomonic cubed-sphere grids and the computational methods of solving the Euler equations on a cubed-sphere.In section 4, we present the results of standard test problems and the performance of the cubed-sphere GCM.In section 5, we summarize and compare our new GCM with other existing ones.

Formulation of the equations
The hydrodynamic equations can be cast into a variety of mathematically equivalent but numerically different forms.We use the momentum equation to elucidate the intricacies of the numerical methods associated with the cubed-sphere geometry, or more generally, with solving the hydrodynamic equations on a sphere.The Eulerian form of the momentum equation is: in which, the ∇v is a second-order tensor that is complex to evaluate on a non-orthogonal geometry.It can also be written in the conservative form: In the flux-form equations, the conservative variables are used, as opposed to primitive variables.Conservative variables are the quantities that are conserved in a control volume in the absence of sources and sinks, and examples include energy and momentum.Primitive variables are quantities that describe the intensive states of fluid directly, such as temperature and velocities.
One can transform the momentum equation into the vector-invariant form using the following vector identity: This formulation eliminates the need for complex metric terms, which typically arise when taking the divergence of a second-order tensor.Moreover, in describing the large-scale circulation of the atmosphere, employing the hydrostatic assumption leads to "layered atmospheres" that can be solved layer-by-layer using the "Vertically Lagrangian" approach that treats each layer as a shallow water system (Lin 2004).Coupled with equations for tracer advection (e.g. using potential temperature for the energy tracer) on a cubed-sphere (Putman & Lin 2007) and a pressure gradient solver (Lin 1997), an extended equation set based on the shallow water equations, the vertically-Lagrangian advection and conservative remapping scheme became the backbone of the Geophysical Fluid Dynamics Laboratory (GFDL)'s FV3 dynamic core (Harris et al. 2020;Mouallem et al. 2023).A similar strategy has been employed by the MITgcm (Adcroft et al. 2004a) where the vector-invariant form is preferred over the conservative momentum equations.Though the MITgcm's cubed-sphere core is restricted to a conformal mapping, and therefore, orthogonal geometry (except at the corners), while the GFDL's cubed-sphere core supports a variety of non-orthogonal mapping choices.
However, the simplicity of avoiding the metric terms comes at the cost of the complex evaluation of the vorticity in a curvilinear geometry.For instance, both the FV3 and MITgcm use the line integral method to assess scalar vorticity.In the case of FV3, velocities are computed on Arakawa-D grids and then interpolated to Arakawa-C grids.Meanwhile, MITgcm computes velocities on Arakawa-C grids and subsequently interpolates them to Arakawa-D grids.The method of calculating the vorticity becomes more intricate when dealing with higher-order methods and three-dimensional non-hydrostatic systems, in which vorticity is a vector rather than a scalar.High-order numerical evaluation of v ×(∇×v) in a three-dimensional system becomes as complex as computing the metric terms themselves.Thus, it is not obvious that the vector-invariant form of the momentum equation is superior to the Eulerian-from or the conservative flux-form of the momentum equation.
Solving the conservative flux-form of the shallow water equations has been pioneered by Ullrich et al. (2010).Let √ g = det(g µν ) 1/2 be the volume element (area element in the two-dimensional sense) of a metric tensor defined by g µν .The flux-form equation reads: where F k is the flux tensor; ϕ is the geopotential; u 1 , u 2 are contra-variant velocities; Ψ is the forcing term including the pressure gradient, the Coriolis force, and the metric terms.A different choice was made by the seminal work by Sadourny (1972), in which the covariant velocities (u 1 , u 2 ) were chosen as the prognostic variable.A similar distinction has been made by the FV3 and the MITgcm models, where FV3 uses the covariant velocities as prognostic variables and MITgcm uses the contra-variant velocities as the prognostic variables (Harris et al. 2020;Adcroft et al. 2004b).This deliberate choice hasn't been extensively detailed in official publications.So, we refrain from speculating on the authors' underlying rationale and instead present our logic here.
In a curvilinear coordinate system, contra-variant velocities are tangential to their respective coordinate lines, while covariant velocities are perpendicular to the other coordinate lines.In an orthogonal coordinate system, such as lat-lon or conformal mapping coordinates, both contra-variant and covariant velocities align in the same direction.However, in non-orthogonal systems like the gnomonic equiangle coordinate system, they are different in direction.For any field (scalar or vector), the transport velocity is dictated by the contra-variant velocity, but the off-diagonal terms of the stress tensor will only vanish if it is of the mixed type (having both upper and lower indices).Let T be the stress tensor in the Euler equations: T µν = ρv µ v ν + g µν p. (6) For a non-orthogonal coordinate system, the metric tensor g µν has off-diagonal terms, meaning that the pressure gradient in one direction forces the contra-variant velocity in another direction as well.Yet, the symmetric stress tensor reads as: where δ ν µ is the Kronecker delta, which ensures this term appears only in the diagonal entries of the 2D tensor.This ensures that the pressure gradient in one direction only contributes to the covariant velocity tendency in the same direction regardless of the orthogonality of the coordinate system.In general relativistic magnetohydrodynamics, a mixed form of contra-variant and covariant indices for the stress tensor (T ν µ ) is also favored because for matrics with an ignorable coordinate X µ , the flux divergence and the metric terms will vanish for T ν µ but not for T µν (Gammie et al. 2003).As a result, we will use the covariant momentum (ρv µ ) as the prognostic variable and diagnose the contra-variant velocity (v ν ) for transport calculation.The relation between the covariant and contra-variant velocity is v µ = g µν v ν .In this article, the Einstein summation convention is used, where an index variable appearing both in the upper and lower indices implies summation over all values of the index variable.

Choice of basis vectors
The choice of basis vector is another often overlooked feature in formulating a cubed-sphere GCM.In differential geometry, the coordinate basis is a vector describing the change of position P with respect to the change of an underlying coordinate variable X, i.e., ∂P/∂X.In the gnomonic equiangle coordinate system, the coordinate variables for each panel are the angles (ξ, η) within the range (−π/4, π/4).The coordinate basis is, thus, not normal and the corresponding velocity is the angular velocity rather than the linear velocity.Thus, there is a choice to be made between the coordinate basis and the normalized non-coordinate basis.
With the coordinate basis, the connection coefficients, or the Christoffel symbols, Γ γ αβ , are symmetric in the lower indices α and β.Moreover, the metric term does not depend on the pressure as Γ γ αβ g αβ = 0, which improves the accuracy of the model.Ullrich et al. (2010) solved the shallow water equations using the coordinate basis of the gnomonic equiangle coordinate system.Using the normalized non-coordinate basis results in the asymmetric connection coefficients and the emergence of pressure in the metric terms, as is the case with solving the Euler equations in the traditional lat-lon geometry.To achieve a stable solution, the numerical evaluation of the metric terms should be carefully examined.Chen (2021) utilized normalized non-coordinate basis, but avoided calculating the metric terms by solving the vector-invariant form of the shallow water equations.
On the other hand, using the normalized non-coordinate basis leads to a simpler metric tensor and improves the geometric calculation and interpretation.For example, the two-dimensional metric tensor of a gnomonic equiangle coordinate system using the coordinate basis (e.g.Ullrich et al. (2010); Ronchi et al. (1996); Nair et al. (2005)) is: with Here we use ( ) a to denote a power law expression rather than a superscript expression in tensor notation.
If we use normalized non-coordinate basis (Chen 2021), the two-dimensional metric tensor of a gnomonic equiangle coordinate system takes a simpler form: where cos Θ = − xy CD (13) Equation ( 12) provides a clear geometric interpretation, suggesting that the coordinate lines form an angle of Θ.This is more intuitive than what's presented in equation ( 8).Additionally, equation ( 12) demands fewer float-point operations, consumes less memory, and thus enhances the model's performance.Since the basis vector is normalized, it is the linear velocity that is solved as the prognostic variable.Area and volume preserve their normal geometric meaning.As a consequence, the numerical solver is almost exactly the same as that in the Cartesian geometry except for the non-orthogonal component indicated by Θ.Since the cross-term in the pressure gradient is eliminated by using the covariant momentum, the only extra step in the non-orthogonal cubed-sphere solver compared to a Cartesian geometry solver is to project the flux to the normal direction of the interface between two finite-volume cells.This provides the possibility of adapting many existing Cartesian-geometry-based models into a spherical GCM.In summary, no single formulation stands out as the universal superior.Each formulation has its own flavor, and every choice inevitably impacts subsequent decisions, e.g. the decision to avoid the metric terms leads to the vectorinvariant form of the momentum equations; the need to accurately calculate the vorticity in the vector-invariant form leads to advancing the velocities on Arakawa-D grids and interpolating to Arakawa-C grids.In this paper, the Arakawa-A grid is used with a normalized non-coordinate basis.

Objectives of this work
The goal pursued in this work is to offer the pathway of adapting models designed for lat-lon or Cartesian grids to the cubed-sphere grid.We describe the minimal recipes needed to facilitate such transformation.This approach effectively bridges the gap between small-scale regional simulations, typically conducted on Cartesian grids, and largescale simulations that must accommodate a planet's spherical shape.Hence, all thermodynamics and microphysical modules from small-scale models can be seamlessly integrated into large-scale frameworks, eliminating the need for redundant developments.
Our base model to build the cubed-sphere dynamic core is the SNAP model (Li & Chen 2019), which extends the Athena++ general relativistic magnetohydrodynamics code (Stone et al. 2020) for atmospheric simulation.We deliberately made subroutines for the cubed-sphere geometry non-intrusive to the original SNAP code, meaning that our model passes all original Athena++ tests and the tests published in the SNAP model.Therefore, our dynamic core, named ExoCubed, serves as a standalone set of subroutines and methods that can be incorporated into any hydrodynamic code as long as the code uses the finite volume method to solve hydrodynamics.
Particularly, the ExoCubed dynamic core assumes that the base hydrodynamic model has a partial or full implementation of the following features: 1. Three-dimensional Structured Mesh.The base model incorporates a mesh in three spatial dimensions, similar to a Cartesian geometry.
2. Dimensional Configuration.The first dimension is the radial (vertical) direction, denoted as X 1 .The remaining two dimensions, X 2 and X 3 , are two horizontal dimensions.An illustration of the geometries of the dimensions can be found in Figure 2.This configuration is compatible with a hydrodynamic model using the spherical polar coordinate system but not particularly compatible with the spherical lat-lon geometry, where the third dimension is typically vertical.
3. Grid Placement.Hydrodynamic variables are placed on collocated grids, also known as the Arakawa-A grid.
4. Reconstruction.The left/right states at the cell interface are constructed via a high-order monotonic preserving scheme.
5. Riemann Solver.First introduced by Godunov & Bohachevsky (1959), Riemann Solvers have become widely used in finite volume methods for compressible flows.The Riemann Solver gives a solution to the Riemann problem, which is an initial value problem that is piecewise constant at either side of the interface face.The shocks and rarefaction waves appear as characteristics of the solution.A Riemann Solver is applied in this model to calculate the flux given the left/right hydrodynamic states.
6. Vertically Implicit.A 1D (vertical) implicit correction/integration scheme is employed to overcome the problem of large aspect ratios.
7. Dual Variables.The dynamic equations are evolved using primitive/conserved variable pairs.
8. Variable Configuration.The primitive variable uses contra-variant velocity and the conserved variable uses covariant momentum.
The more features the base model supports, the easier it will be to adapt to the cube-sphere geometry.In the context of the finite-volume numerical method, as discussed in LeVeque ( 2002), a standard high-velocity, compressible hydrodynamic code in the Cartesian geometry generally satisfies all features.Then, extending the base model to the cubed-sphere geometry requires the following extra steps: 9. Mesh on Six Panels.The spatial dimension is divided into six panels, with the computational domain divided into two panels along one axis and three panels along another.The limits of the spatial dimension are arbitrarily defined in this case.The coordinate lines and grids are uniformly distributed within a panel but are discontinuous across panel boundaries.Such design is to fit the Cubed-Sphere geometry into the Athena++ mesh and parallelization framework.
10. Topological connection.Since the spherical shell is divided into 6 panels, the topological connections between the panels are different from their relative placement in the coordinate axes.The intra-panel communication is consistent with that in the Cartesian geometry, but the inter-panel communication needs to be specially dealt with by rotation and interpolation (Figure 2).
11. Velocity rotation and interpolation.Ghost zones are layers of extra computational cells outside of block boundaries used for reconstruction but not updated by this block.The ghost zones of one panel cross the panel boundary, and the velocities defined on the neighboring panel need to be rotated and interpolated to the ghost zones of the current panel.
12. Left/right states synchronization.At the interface between two panels, the left/right states reconstructed by the neighboring panels may differ due to the velocity rotation and interpolation.We devised a synchronization step that uses the right state reconstructed by the right panel and the left state reconstructed by the left panel to be the left/right states at the panel interface to ensure that the flux calculation is consistent across panel boundaries.
13. Projection and de-projection.The dynamic core is designed to be compatible with any existing Riemann Solver.So, we project the left/right states to the local orthonormal coordinate system defined by the cell interface and use the existing Riemann Solver to solve for the fluxes across the interface.Then, we de-project the fluxes to the original coordinate system on each panel (Figure 4).
14. Metric Terms.We implement the metric terms as the source terms for the momentum equation.We make sure that no spurious flow is generated for a uniform pressure and density field on a sphere.
In the next section, we discuss in detail the theory of solving the conservative form of the compressible Euler equations on a gnomonic equiangular cubed-sphere geometry defined by normalized (units) basis vectors.We elaborate on steps 9 ∼ 14, beginning with an introduction to the gnomonic equiangular coordinate system and ending with the metric terms.

The gnomonic equiangular coordinate system
The grids of the cubed-sphere geometry are obtained by projecting the rectangular grids on the surface of a reference cube radially onto the surface of a concentric sphere (Sadourny 1972;Ronchi et al. 1996), known as the gnomonic projection.There are a few options for choosing the discrete grids on the reference cube.One of them is called the equiangular grids, where the face of the reference cubed is tessellated into N 2 number of rectangular grids defined by two angular coordinates {ξ, η} (Ronchi et al. 1996).Suppose that the reference cube is bounded by concentric with the sphere, {ξ, η} are uniform distributed angular variables to span the range (−π/4, π/4), such that: where (j, k) are both integer numbers from 1 to N to denote the index of a finite volume cell.Thus, the horizontal span of a finite volume cell is bounded by with the center at (ξ(j), η(k)).The tangent values of the angular variables {ξ, η} define the local Cartesian coordinates of the grid points on the surface of the reference cube: with the derivatives satisfying The coordinates described above are defined locally at each face of the reference cube.Therefore, there are six local Cartesian coordinate systems.We give numbers and names to those faces and define a global Cartesian coordinate system that aligns with the local coordinates of the face at the top.Let face 1 be the face at the top, face 5 be the opposite face at the bottom and faces 2,3,4,6 be the faces at four sides.Let the z-axis point toward the top, x-axis toward the front and y-axis toward the right, the local to global transformation of the Cartesian coordinates can be obtained by vector rotation and is summarized in Table 1.Note that we have used (z, x, y) ordering rather than the traditional (x, y, z) ordering of the coordinates due to the convention that the first axis is vertical (radial).Projecting the local rectangular grids (x, y) on each face of the reference cube to a concentric sphere results in a curvilinear coordinate system defined on each panel of the cubed sphere.Each point on a sphere is thus described by the radial distance to the center r, two angular coordinates {ξ, η}, and the panel (face) number n.Similar to the notation used in Ronchi et al. (1996), we define the following auxiliary pre-computed variables for calculating the geometric relations: We present the calculation of the geometric relation of panel 1, the top panel containing the north pole.The geometric relation of the other panels can be derived similarly by rotating the coordinate axis.The starting point of constructing the cubed-sphere coordinate system is the position vector of a point on a spherical panel, whose Cartesian components are: The expression of the position vector can be similarly obtained for an oblate spheroid.In problems of rapidly rotating giant planets, the projection of an oblate spheroid becomes a more proper representation.We leave the discussion of the oblate spheroid in Appendix C and continue with the discussion of the spherical case in the main article.
From the position vector, we obtain the coordinate basis vectors {z 1 , z 2 , z 3 } by differentiating the position vector with respect to each coordinate variable: As discussed in Section 2, we favor basis vectors that are of unit length to be compatible with the existing Cartesian model.We define the lengths of the coordinate basis vectors as: Thus, the normalized non-coordinate unit basis vectors are: Now we switch to the generic tensor notation to calculate the metric and the connection coefficients defined by the unit basis vectors.Let the generic coordinate variables {X 1 , X 2 , X 3 } be: The covariant metric tensor is given by the pairwise inner products of the basis vectors: and the contra-variant metric tensor g µν is the inverse of it: We also define as the determinant of the covariance metric tensor and its square root, √ g = δ/(CD) = sin Θ, the volume element.
According to the definition of the basis vectors in equations ( 26) -( 28), the connection coefficients of the unit basis (Christoffel symbols) are: with the non-zero terms being: In a perfect spherical geometry, the radial direction aligns with the vertical direction opposite to the gravitational acceleration vector.So we will use the radial direction and vertical direction interchangeably in this article.We also define the X 1 -face as the two-dimensional surface that lies on a constant value of X 1 and likewise for the X 2 -face and the X 3 -face.

Equation of motion in a curvilinear geometry
We solve the standard compressible Euler equations in the curvilinear geometry on a rotating sphere under a potential field.The conservative form of the momentum equation reads as: where ρ is the density, v the velocity vector, Ω the angular velocity of the planet rotation, p the pressure and ϕ the potential field.Expanding the momentum equation in the vector form Equation (39) to the component form defined by the covariant basis, we have the following form after the arrangement: Here ε αβγ is the Levi-Civita symbol, which is used to write the tensor notation of a cross product: U×V = ε ijk U i V j êk (here e k is the basis vector in k-direction).The divergence theorem on the second-order tensor states that (Grinfeld 2014): Thus the momentum equations in the covariant basis are written as: where is the total stress tensor and the terms on the right-hand-side of the equation above are the Coriolis force, the metric term and the gravitational force, respectively.Note that pressure only exists in the diagonal terms due to the choice of mixed lower and upper indices (See the discussion in Section 2).
Similarly, the continuity and energy equations without external forcing are: where e is the sum of specific internal energy and kinetic energy.For example, for an ideal gas, it is defined as: with γ being the adiabatic index, calculated as the ratio of the constant-pressure heat capacity to the constant-volume heat capacity of the gas.Let q = (ρ, v 1 , v 2 , v 3 , p) denotes the hydrodynamic primitive variable, Q = (ρ, ρv 1 , ρv 2 , ρv 3 , ρe) the conserved variable and F , G and H be the flux densities in three dimensions respectively.It is illustrative to write their components explicitly as: whereis h is the sum of the specific enthalpy and kinetic energy, and for an ideal gas Substituting the fluxes in equations ( 39), ( 44), ( 45) by equations ( 47) -( 49), we arrive at a compact form of the Euler equations in the curvilinear coordinate system: where Φ represents various forcing terms on the right-hand-side of equations ( 39), ( 44), (45).

Cubed-sphere finite volume discretization
In the finite volume method, we definite integrals of a finite volume (V ijk ) and three face areas There is a √ g factor in equations ( 52) and ( 53) because the X 2 -axis and X 3 -axis are at an angle Θ, while the X 1 axis and X 3 axis are orthogonal to each other.An illustration of the relation of the three directions can be found in Figure 2.
Similarly, we can define the volume-weighted average of the conserved variable as: and the area-weighted average of the normal fluxes: We would like to emphasize that the fluxes: are perpendicular to each face area of the finite volume.This allows the use of any Riemann Solver that is supplied by the base model without rewriting a new one for the non-orthogonal geometry.
Finally, we use the superscripts {n, n + 1, * } to denote discrete time steps.The finite volume discretization of equation ( 51) with explicit spatial time stepping reads as: where the normal flux densities F ⊥ , G ⊥ , H ⊥ are evaluated at some time between n and n + 1 depending on the time integration scheme.It is clear at this moment that the finite volume model in a cubed-sphere geometry bears the same numerical scheme as any Cartesian model.This feature makes the transformation of a Cartesian model to a spherical model extremely easy if one correctly calculates the normal fluxes using a Riemann Solver and the finite volume and areas have their usual geometry meaning (facilitated by unit bases vectors).
The flux density can also be evaluated implicitly using a prediction-correction method, in which a prediction flux density is calculated using any explicit method, and the implicit correction flux density is obtained by linearization.Implicit flux densities are usually employed for the vertical dimension of a GCM because the vertical grid spacing is much smaller than the horizontal grid spacing.In ExoCubed, we use the VIC scheme developed by Ge et al. (2020) to solve the implicit flux density for the vertical dimension only, denoted as F ⊥, * and use the explicit flux densities for two horizontal directions.Because the vertical direction is perpendicular to the two horizontal directions, no change is needed for the implicit scheme.

Projection and de-projection in Riemann Solvers
Given the left/right reconstructed primitive variables (q − , q + ), the Riemann problem is usually solved in an orthonormal basis frame, which yields the fluxes normal to the face boundary between two finite volume cells.Since the coordinate line is not orthogonal to the face boundary in the horizontal directions in the cubed-sphere geometry, we project the primitive variables to a local orthonormal basis which is defined by the direction of the face boundary and the vertical direction.The projection is a specialization of solving the general relativistic MHD problem (White et al. 2016) but generalizes a two-dimensional projection (Ullrich et al. 2010) designed for the shallow water model to three dimensions.
Denote three basis vectors of the orthonormal basis frame at the cell interface as {e 1 , e 2 , e 3 }.The first basis vector e 1 is the same as the one in the original basis, b 1 : since b 1 is already normalized, we do not need further normalization.If the face boundary is spanned by {b 1 , b 3 } (X 2 -face), then its normal direction is the same as the contra-variant basis, b 2 because: So, the second basis vector is defined as the normalized b 2 such that: The third basis vector is the same as b 3 because b 3 is normalized: The orthonormal basis frame for the X 3 -face can be constructed similarly using b 3 as the direction for the face normal.
Let T be the velocity or flux density vector.Component-wise relations between the original basis and the orthonormal basis are: The components of the transformed vector O 2 T are obtained by the inner product of the T vector and the new basis vectors: So the projection operator O 2 can be written in the matrix form as: and its inverse is The component-wise de-projection formula reads as: Similarly, the projection O 3 is and its inverse: Note that the vertical component is unaffected by the projection since the vertical axis in the orthonormal frame is the same as that in the original frame.The projection matrix O 2 is also derived by Ullrich et al. (2010) in a two-dimensional shallow water manner.Now we describe the whole process of projection, Riemann Solver, and de-projection using X 2 -face as an example.A corresponding illustration of this process can be found in Figure 4. First, we obtain the left and right components of primitive variables in the original basis at the face boundary by a reconstruction method based on cell-averaged values, denoted as: Step 1: Here j − 1/2 refers to the left face of a finite volume (See Fig. 3).Second, the projected components in the local orthonormal basis are obtained by the projection matrix (70): Step where "•" indicates a matrix multiplication.
Step 2 yields two primitive variables, O 2 q ± j−1/2 , that are defined in the local orthogonal reference frame.Third, we solve the orthonormal Riemann problem in the local orthonormal frame such that the resulting flux densities are: Step Step 3 yields the normal flux densities across the cell interface face but the components are defined in the local orthogonal reference frame.Finally, we de-project the normal flux density from the local orthonormal frame to the original reference frame: Step 4: Step 4 yields normal flux densities that are expressed in the cubed-sphere reference frame.The result of step 4 is then supplied to equation (51) to calculate the flux divergence.
We have verified the validity of the projection and de-projection algorithm using a simple circular breaking dam simulation using the affine coordinates, benchmarked by the simulation results in Cartesian coordinates (not shown).

Velocity rotation and interpolation
Crossing the panel boundaries, the ghost cells of one panel do not align with the interior cells of the neighboring panel.So, when the values in the interior cells of the neighboring panel are used to fill the values in the ghost cells of the current panel, interpolation is necessary (Figure 5).In the gnomonic equiangular cubed-sphere geometry, this problem is simplified to a one-dimensional linear interpolation problem because each coordinate line is on a great circle of the sphere.Denoting the positions as angular positions along the arcs.The angular position is calculated as Here N is the number of cells along the direction of interpolation, and n is the corresponding index of the cell.k denotes the distance of the cells from the boundaries.When k < 0 or k ≥ N , the cell is within the ghost zone.
The interpolation happens in a matching manner.To interpolate the values on ghost zones (k < 0), the locations for interior points k ′ = −k − 1 are used.Note that Equation ( 80) is universal and has the same form for all the panel boundaries.
Not only do the values in the ghost cells need to be interpolated, the velocity components also need to be transformed since the coordinate lines are not continuous crossing panel boundaries.
Here R is the rotation matrix that rotates the local coordinate system counterclockwise, and The number of such rotations needed, i.e. the value of n in each panel boundary is summarized in the following table: Summary of numbers of rotations needed for transformation on panel boundaries.
For calculation, we need to first find the coordinates of interpolated points in both panels.For the target panel (where the point is within the ghost zone), the grid is apparent, while in the source panel, we can calculate the y-coordinate values from the angular position α n : and The sign of the y-coordinate abides with the sign in the target panel.Note that we have assumed the panel boundary is along y-direction in this case, but for boundaries along any direction the interpolation is similarly treated.After this step, we use Equations ( 26) -( 28) to calculate the local basis vectors and calculate the transformation matrices.

Left/right states synchronization
Since interpolation and rotation are performed to fill the ghost zones outside the panel boundary, the restricted states at the panel interface may be slightly different for the two panels sharing one cell interface.If the left/right states are different, the result of the Riemann Solver will also be different, leading to inconsistent fluxes across the panel boundary.The GFDL's FV3 model solves this issue by averaging the fluxes computed by the left and right panels (Chen 2021;Mouallem et al. 2023).We take a different approach.We implemented a communication function that exchanges the reconstructed left/right states in both panels rather than fluxes.Specifically, the right state at the cell interface gets the right values reconstructed by the right panel, and the left state at the cell interface gets the left values reconstructed by the left panel.This approach minimizes the use of the values in the ghost cell, which are inaccurate.After the exchange is done, the two states at the left/right side of the panel interface are synchronized for both panels.This ensures that the Riemman Solvers in each panel yield the same flux, and thus improves the conservation properties.

Metric terms
Due to the change of the basis vectors, the momentum equation ( 42) contains extra metric terms acting as the source term to the local momentum.The expression of the metric terms depends on the type of coordinate system, the choice of basis vectors and the variables to solve.In the classic geophysical fluid dynamics textbook (Pedlosky 2013), the metric terms only involve the components of the Reynolds stress tensor because the differential equations are written in primitive variables.Since we solve the conservative form of the equations, our metric terms are calculated from the total stress tensor which includes the contribution from the pressure.In spherical-polar coordinates, the metric terms are: where θ is the polar angle and v 1 , v 2 , v 3 are velocities in the radial, polar and azimuthal directions respectively.In cubed-sphere geometry, the metric terms are: Substitute equations ( 36) -( 38) and ( 43) into equation ( 88), we derive the metric terms in the cubed-sphere geometry: Equations ( 89) -( 91) can be compared directly to equations ( 85) -( 87).First, we find that a kinetic energy term (v 2 v 2 +v 3 v 3 ) appear in both equations ( 85) and ( 89).Second, the factor, cot θ, represents the change of face area in the spherical polar geometry whereas, the factors x/(rD) and y/(rC) represent the change of face area in the cubed-sphere geometry.Metric terms in horizontal directions are more symmetric in the cubed-sphere geometry than those in the spherical-polar geometry due to the symmetric nature of the coordinate system.

BENCHMARK TESTS
In this section, we present tests to evaluate the performance of the dynamic core.From the benchmark simulations, we can draw conclusions on the accuracy, efficiency, and stability of the cubed-sphere model presented in this work.We also evaluate how well the different physical quantities are conserved.
Tests range from shallow water tests (2D global models) to 3D general circulation models.In 3D general circulation model tests, we applied the Vertical Implicit Correction (VIC) scheme (Ge et al. 2020) to further improve the speed of simulations.It is worth noting that, with VIC applied, the relationship between the time steps and Courant-Friedrichs-Lewy (CFL) number is calculated as ∆t V IC = CFL × min(∆x, ∆y)/c s since the implicit scheme is applied along the X 1 -axis and is unconditionally stable.

Hydrodynamic blast in one cubed-sphere panel
For the first test, the experiment used in Chen ( 2021) is adopted to validate the dispersion and dissipation properties.This test involves simulating a splash of a sinusoidal droplet on a non-rotational sphere.The initial perturbation to a uniform height field is limited to a small region at the pole.The initial geopotential height is written as: Here r is the distance from the north pole, R is 500 km, Φ 0 = 50 m × g, and Φ ′ = 1 m × g.These values give a gravity wave that reaches the south pole in 10 days.In Figure 6, the height fields at day 5, when the wave reaches the equator, are shown in the left column.
For the tests, we conducted the simulation in multiple different horizontal resolutions.In this paper, we use the letter "C" followed by a number, which is the number of grids along each direction for each of the six panels, to label the resolution of each configuration.For example, "C48" denotes the resolution of 48×48×6 on the surface of the sphere.
The quantitative comparison between grids and Chen (2021) (gnomonic equiangular coordinate with rk3 scheme results are shown here) is shown in the right column of Figure 6.As expected, the solutions are less dissipative and preserve the spherical symmetry better by having a smaller zonal wind magnitude with finer grids.Intriguingly, our cubed-sphere model is less dissipative in comparison to Chen (2021) at finer grids.Williamson et al. (1992) provided a widely-used standard test set for testing the accuracy of shallow water equations on spherical geometry.For shallow-water tests, we performed tests on test cases 2 (global steady state nonlinear zonal geostrophic flow) and 6 (Rossby-Haurwitz Wave) to evaluate the order of accuracy of our model.Note that in Williamson et al. (1992), there are four different configurations for case 2 with jets going along different directions, but here we show only one case, where the jet is tilted 0.05 radians from the perfectly zonal direction, since the results from all four configurations are essentially similar.

Rossby-Haurwitz Wave and Steady Zonal Jets
These two tests are conducted with different resolutions, from C48 to C384, and we use the results from C768 to calculate the normalized global errors l 1 and l 2 in layer thickness: Here h T is the true solution, and I[] is the operator that integrates over the surface of the sphere.The respective results are shown in Figure 7.The log-linear fits are shown in the plots, and the slope values correspond to the order of accuracy.In the tests, we used CFL=0.9, while we tested CFL=0.1 and found no significant difference in both errors and the fitted slope values.This insensitivity to a nine-fold difference in time step has important implications for the 3D simulations.The vertical implicit scheme ensures unconditional stability in the vertical direction, albeit at the cost of damping fast waves like sound waves.Yet, for long-period waves like Rossby waves, the time step is significantly smaller than the wave period, resulting in minimal distortion of wave propagation by the vertical implicit scheme.Non-hydrostatic models with vertical implicit schemes should resolve slow waves as well as any hydrostatic models.This contention is validated by the findings in Ge et al. (2020), where 3D simulations with the implicit scheme turned on and off showed no discernible differences.This supports the use of larger time steps enabled by the vertical implicit scheme to enhance computational efficiency.
We found that the model's convergence order in the two test cases is about 2, even though the 5-th order numerical scheme (WENO5) is applied in the reconstruction.Currently, the convergence order is limited by the second-order cell averaging scheme when converting between the conserved and primitive variables.Stone et al. (2020) similarly reported second order of accuracy in general, even if a Piecewise Parabolic Method (PPM) is used for the reconstruction.While the second-order cell averaging scheme poses a second-order convergence rate to the whole model, applying higherorder reconstruction schemes has the advantage of reducing both the amplitude and the phase error at a fixed model resolution.Specifically, in the hydrodynamic blast test, we noted that the lower-order reconstruction schemes are more dissipative than higher-order schemes such as WENO5.Since the viscosity of the atmospheric motion is small, using a higher-order reconstruction scheme is essential for resolving hydrodynamic instability in a weakly dissipative system.For example, in a standard rising bubble test (Robert 1993), where the temperature anomaly is only 0.5 K, a 5th-order reconstruction scheme is necessary to resolve the Kelvin-Helmholtz instability caused by the rising bubble.Yet, in the standard sinking bubble test (Straka et al. 1993), where the temperature anomaly is about 20 K, a 2nd-order reconstruction scheme suffices to resolve the instability.
In Figure 8, we show the height field generated by different levels of resolutions, at days 40, 80, and 100, respectively.The Rossby-Haurwitz wave is supposed to be a wave that travels from east to west without changing the shape, but  the error accumulates, and the wave breaks up at later dates.For finite volume models, the truncation error makes the wave more susceptible to dynamic instability, and it usually breaks up at about days 30-40 (Smith & Dritschel 2006;Thuburn & Li 2000).The spectral semi-Lagrangian models have a chance to keep stable till around day 100 (Thuburn & Li 2000).In Figure 8, even the coarsest grid C48 survives beyond day 80, and for various resolutions, the cases all become unstable at about day 90-100.Such solutions show a slower error accumulation than other finite volume methods previously reported.

Held-Suarez test of Earth's climate
The Held-Suarez test of Earth's climate is designed as a common problem for comparisons between different dynamical cores of GCMs (Held & Suarez 1994).Focused on the global-scale atmospheric features and the long-term state of an earth-like planet, the comparisons are usually made on a final quasi-steady state, and averaged quantities are used.The final quasi-equilibrium state of this test presents a latitudinal temperature gradient and large scale winds.This is similar to the atmosphere of the earth.
In addition to the Coriolis force corresponding to the rotation rate of the earth, the setup of this simulation involves two major forcing mechanisms.A reference thermal structure is imposed, and the atmosphere is relaxed to the reference structure by a Newtonian cooling scheme.This forcing simplifies the radiative forcing processes: Following Held & Suarez (1994), we use the normalized pressure scale σ = p/p 0 here, and p 0 = 1bar is the surface pressure.In addition, c v is the isochoric specific heat of air, and k T (in the unit of day −1 ) and T e q (in the unit of K) are the temperature damping strength and the reference thermal structure is defined as k T (ϕ, σ) = 0.025 + 0.225 max 0, σ − 0.7 0.3 cos 4 ϕ (96) T eq (ϕ, σ) = max 200, 315 − 60 sin 2 ϕ − 10 log σ cos 2 ϕ σ κ (97 where κ = R/c p = 2/7.The Rayleigh drag is applied to the lower boundary to simplify the boundary effects from the interface with the ground.It appears as an additional source term in the momentum equation The parameter k v , the Rayleigh drag strength, is calculated as (in the unit of day −1 ) The simulation domain is set to have a 40×40×6 resolution along the horizontal direction.This has a similar number of cells as a spherical-polar grid with an angular resolution of 2.8 • along both longitudinal and latitudinal directions.The domain is 25 km in height, with 40 layers, giving a vertical resolution of 625 meters.Different from Ge et al. (2020); Mendonça et al. (2016), we have not applied a sponge layer at the top of the atmosphere to absorb the vertically propagating waves, as the stability of the model is not compromised by these waves.For the timestep, we use a high CFL number of 0.9 for high efficiency and also highlighting the stability of the dynamic core.Similar to Ge et al. (2020) and the original work, the simulation reached statistical equilibrium at around day 200, and the simulation continued until day 1200.The statistics of the flow fields for the 1000 days between day 200 and day 1200 are presented here.The timestep ∆t for a spherical-polar grid with 2.8 • angular resolution is about 25 seconds, while the ExoCubed dynamic core gives 425 seconds, providing an acceleration over an order of magnitude.Delays due to additional MPI communications needed for the flux synchronization steps exist but are negligible.
The averaged quantities are plotted in Figure 9. From the average meridional wind v, the Hadley cells and Ferrel cells are clearly defined.The maximum prograde jet has a magnitude of 30.41 m s −1 , and the maximum meridional wind magnitude is 2.18 m s −1 .This magnitude of wind maximum is only 1% from the values reported in Ge et al. (2020), and the average wind speeds as well as the average temperature structure agree with the results reported in previous studies (e.g.Ge et al. (2020);Lin (2004); Ullrich & Jablonowski (2012); Mendonça et al. (2016)).In addition to the averaged large-scale circulations, the eddy-related statics are also produced.The meridional and vertical eddy flux combined presents the transport related to the well established Hadley and Ferrel cells.Additionally, the eddy fluxes and the eddy kinetic energy patterns agree with the results in Ge et al. (2020).
To better demonstrate the patterns of the general circulation, the streamfunction is calculated based on the averaged velocity fields.There are several pairs of blobs marking major circulation patterns.There are four cells below σ = 0.2, which demonstrate the contours followed by Hadley cells (closer to the equator) and Ferrel cells (farther from the equator).The results agree qualitatively with the results reported by Mendonça et al. (2016).
A key to note in this model presented is that, at the top of the atmosphere, a sponge layer was not applied, and the top boundary is treated as a solid boundary.This generates an additional pair of convection cells above the Hadley cells, at σ < 0.2.Such structure is seen in both the streamfunction plot (Figure 10) and the meridional eddy flux plot (Figure 9).Bretherton (1969) analytically investigated the scale height of the Lamb wave pattern H Lamb , and if we compare it to the pressure scale height H p , the ratio is only associated with the heat capacity ratio γ: This agrees with the observation here that the additional cells extend from σ = 0.01 to σ ≈ 0.2, so the results are reasonable interpreted by the lamb waves.Where R p is the radius of the earth (6371 km).The change of total AAM and the total energy is approximately 1.1% and 0.2% respectively over the period of simulation.The change of total mass is 10 orders of magnitude smaller than the initial value.Hence, we conclude that the conservation laws are well followed in this simulation.

Shallow Hot Jupiter simulation
We present a final test with global circulation model on the hot jupiter atmospheres.Observational selection effects lead to most known Jovian-size extrasolar planets being close to their stars.These hot jupiters are likely tidally locked to their host stars, orbiting around them synchronously.This brings about an interesting radiative forcing pattern, as the dayside is always irradiated.This permanent day-night irradiation patterns bring a climate state distinct from those observed in the solar system, making it an interesting problem to simulate.We use this as a test case to validate the capability of our model to simulate an atmosphere that is in a completely different regime.It has been shown that the equator super-rotating jet could have a velocity of about 1000 m s −1 , which is in the subsonic to sonic regime, and the temperature difference between the day and night sides can exceed 1000 K (Showman et al. 2015).
The benchmark hot jupiter test is introduced by Heng et al. (2011).It follows from a shallow three-dimensional model of the hot jupiter's atmospheric circulation (Menou & Rauscher 2009).In this simulation, we focus on a planet with radius R p = 10 5 km and gravitational acceleration of 8 m s −2 .The initial bottom pressure is set to be 1 bar, and the initial atmosphere is isothermal at 1600 K. Similar to the Held-Suarez test, the radiative forcing is simulated by relaxing the atmosphere to a prescribed temperature structure with Newtonian cooling, and the reference temperature profile is given by (Menou & Rauscher 2009) Here T vert and β trop are parameters in the same form as introduced in Menou & Rauscher (2009).T E−P = 300 K is the temperature difference between the equator and the pole.The number of grids is the same as used in the Held-Suarez test, with 40×40×6 in the horizontal direction, and vertically 4500 km with 40 layers.The simulation takes longer (3000 earth days) to reach statistical equilibrium, which is apparent from Figure 13.Hence, the average flow fields shown are average values from days 3000-4000.Similar to Figure 9, the same set of flow field quantities are shown in Figure 11.The maximum speed of the prograde zonal wind is 1175.13m s −1 .The general zonal mean structures, especially the well-reported average temperature and zonal wind profiles, agree well with previous results (Ge et al. 2020;Mendonça et al. 2016;Heng et al. 2011).Subtle differences exist, and differences can be attributed to the different capabilities of resolving structures by different grids.
In Figure 12, we present the streamfunction of the averaged flow field.The streamfunction of the shallow hot Jupiter simulation was presented in Mendonça et al. (2016).The structure of the Hadley cells is similar, but the relative strength of the Ferrel cells is much stronger than the previous results.In Mendonça et al. (2016), an additional convection cell structure is present at σ > 0.9.This structure is not observed in our model.
In Figure 13, the conservation of physical quantities clearly shows the long spin-up of the model.The time taken for the model to reach equilibrium is about 3000 days, an order of magnitude longer compared to the Held-Suarez test.As expected for a finite volume method, the total mass conserves well, and no change is observed up to the numerical precision.During the simulation, the total AAM and total energy rise by 34.1% and 1.9%, respectively.

CONCLUSION
In this paper, we introduced a cutting-edge Riemann-Solver-based Finite-Volume Cubed-Sphere dynamic core, Ex-oCubed.We presented the details and validated the model with various tests ranging from shallow-water tests to general circulation models.We validated the accuracy of the model with standard tests, with or without Coriolis force on a shallow-water sphere.For the general circulation models, we simulated two very different planets, one earth-like (Held-Suarez test) and the other being a hot Jupiter.We provide a variety of Riemann solvers that adeptly manage large temperature/pressure gradients, shock waves, and refine the depiction of rapid atmospheric transitions and high-velocity wind events.The shallow hot Jupiter tests have validated that our model has the capability of dealing with flows in subsonic/sonic regime.Our model can be seamlessly applied to other planetary atmosphere simulations, varying from thin atmospheres such as Mars, Titan, Pluto, super earths, and lava planets, to extended atmospheres relevant to Jupiter, Saturn, Uranus, Neptune, extrasolar gas giants, and brown dwarfs, where projection to an oblate spheroid may be needed (Appendix C).
The ExoCubed dynamic core employs the Arakawa A-grid (Arakawa & Lamb 1977), which is an unstaggered grid where all variables are defined at the center of a finite volume.There was a misconception that Arakawa C-grid yielded better phase-propagation properties than Arakawa A-grid.The above statement is only true when the numerical scheme is second-order and when the flow is smooth (Harris et al. 2021).Xu et al. (2021) showed that the phase error of shallow-water waves in unstaggered grids (Arakawa-A) is greatly reduced or eliminated for higher-order numerics.Moreover, when the flow is discontinuous, the C-grid staggering produces far more numerical noise at a discontinuity than an unstaggered grid (Chen et al. 2018).
We use the finite volume method to solve the hydrodynamic equations.The design of finite volume methods is based on the control volume and fluxes, which has superior performance in the conservation of the physical quantities, particularly the variables directly solved in the governing equations, such as mass, momentum, and energy.For this reason, the finite volume method is widely used by many existing GCMS (e.g.GFDL FV3, MITgcm).However, unlike the prevailing GCMs, we pioneer the use of Riemann Solvers to solve for fluxes across the cell interface.Using the projection and de-projection method, the ExoCubed dynamic core is agnostic with the selection of the Riemann Solver.For example, the Roe solver (Roe 1981), the Harten-Lax-van Leer-Einfeldt (HLLE) solver (Harten et al. 1983), and its variation HLLC solver (Toro et al. 1994) are all viable choices.The Riemman Solver applies minimum numerical viscosity to stabilize the numerical solution but is computationally more expensive than traditional methods such as using hyper-viscosity.For pure hydrodynamic problems, the ExoCubed core is slower than other GCMs using hyperviscosity to dissipate grid-scale noise.Yet, in many cases, the latency is dominated by the computation of physical modules rather than the hydrodynamics itself.So in a full-fledged GCM with comprehensive physics, radiation, and chemistry modules, the ExoCubed-based GCM is expected to perform similarly in terms of speed compared to other models.
To tessellate a sphere, we choose to use the gnomonic equiangular cubed-sphere geometry.This allows a seamless integration with an existing Cartesian geometry based hydrodynamic model.Other options have been explored in the literature, such as using triangular and hexagonal grids (e.g.Mendonça et al. (2016)).The quadrilateral formulation employed here has the advantage of being simple and fast for solving problems in simple geometry.With the transformation of grids to a gnomonic equiangular cubed-sphere, each panel on a sphere is essentially rectangular, and the calculations and parallelization become almost the same as those in Cartesian geometry (see equation 61).Utilizing the quadrilateral formulation simplifies the implementation of high-order reconstruction methods.This is achieved by extending the width of the stencil, in consonance with the decision to use the Arakawa A-grid.

B. CORIOLIS FORCE
The Coriolis acceleration differs depending on the cubed-sphere panel.In the Cartesian coordinate with (z, x, y) ordering, the planetary vorticity is 2Ω = (2Ω, 0, 0).In a matrix form, the Coriolis acceleration for the polar panels is: and for the equatorial panels is: For problems related to giant planets and brown dwarfs, it may be useful to use the projection onto an oblate spheroid instead of a sphere.In a spherical geometry, the vertical (radial) direction is special because it is normal to the other two horizontal directions.Thus, the covariant and the contra-variant basis vectors, b 1 and b 1 are the same.In an oblate spheroid, it is crucial to distinguish between the "vertical" direction, which is normal to a geopotential surface, and the "radial" direction.The "vertical" direction is along the contra-variant basis vector, b 1 and the "radial" direction is along the covariant basis vector, b 1 .This definition is consistent with the definition of basis vectors in the spherical geometry (Equations 26 -28).The methods described in the main article can be easily extended to an oblate spheroid by generalizing the relation between the radial direction and two horizontal directions.We outline the generalization in the following paragraphs.
For an oblate spheroid with eccentricity e, equation ( 20) and the following derivations take a different form.Here we consider a "unit spheroid" in the Cartesian coordinate (z, x, y): With the formula of the position vector determined, following the same line as Equations ( 20)-( 28), the basis vectors can be determined.Subsequently, the metric tensors and the Christoffel symbols can be found in Equations ( 33)-(38).
It is worth noting that the direction of the effective gravitational acceleration is along the "vertical" direction rather than the "radial" direction in an oblate spheroid geometry.So, effective gravitational acceleration has only one component in the contra-variant basis vectors and has all three components in the covariant basis vectors, i.e.: Equation (C31) reveals an additional advantage of solving for covariant momentum (ρv α ) instead of contra-variant momentum (ρv α ): in a covariant formulation, all three momentum equations remain unchanged, with gravitational acceleration having no impact on the horizontal covariant momenta, even in the context of oblate spheroid geometry.Consequently, the equations of motion, Equations ( 39)-( 61), remain unchanged.
Next, in the projection step within the Riemann Solver, the basis vectors b 1 and b 1 are no longer equivalent.Hence, projections are needed for all surfaces in three directions.The projection and deprojection relationships can be inferred in the same way as Equations ( 62)-( 75), projecting to a local coordinate basis formed by the normalized contra-variant basis normal to this surface and other two covariant bases.After the projection, the Riemann Solver step described in Equations ( 76)-( 79) also remains the same.
Despite that the oblate spheroid has a more complex geometry, the convenience that only a one-dimensional interpolation is needed for ghost zones is preserved.The formula and intuition provided by Equation (80) and Figure 5 remain valid.The transformation of vectors across the panel boundary can be similarly derived, as b 1 is unchanged when a panel boundary is crossed, so the transformation is still from Finally, the metric terms can be calculated according to Equation (88), with more complex Christoffel symbols calculated for the oblate spheroid coordinate basis.
One final generalization involves the "vertical" implicit solver.In the spherical geometry, the "vertical" velocity is identical in both its covariant and contra-variant forms.However, in an oblate spheroid, the velocity differs between the "vertical" direction (contravariant component) and the "radial" direction (covariant component).The original vertical implicit solver was developed for a spherical-polar grid.New implicit solver should be developed using covariant momentum in an oblate spheroid.

Figure 1 .
Figure1.The trilemma for a GCM constructed on a structured gird using explicit time-stepping.

Figure 2 .
Figure 2. The flattened geometry (left), the folded cube (middle), and the cubed sphere geometry (right) of the six panels of the cubed-sphere.The outer panel boundaries noted with the same color are in touch.In the enlarged view, the coordinate basis of panel 1 (blue) and panel 2 (black) are shown.
Unit basis vectors {b 1 , b 2 , b 3 } constitute the non-coordinate basis at each point on the spherical panel.On each panel, b 1 points at the radial direction normal to b 2 and b 3 , but b 2 and b 3 are at an angle Θ to each other, which is:

Figure 3 .
Figure 3. Illustration of non-orthogonal curvilinear geometry.Cell centers are represented by black dots, while face centers by colored dots.Cell centers are marked with integer indices, while cell faces by half-integer indices.Curvilinear coordinate lines are presented in pink and green colors, with the non-orthogonal angle highlighted in blue.

Figure 4 .
Figure 4. Illustration of the steps of the projection operations for Riemann Solver calculations.We convert the velocity on the covariant basis to the local orthonormal basis before the Riemann Solver calculation (left) and convert the flux from the local orthonormal basis to the covariant basis.

Figure 5 .
Figure 5. Grid visualization of ghost cells at the panel boundary.Panel A's coordinate lines and cells are depicted in pink, while Panel B's are in green.The delineating boundary between Panels A and B is shown in black.Cell centers are represented by dots.Ghost cells in Panel B, designated by open green circles, will be populated using cell-averaged values from Panel A. Bold pink and green lines indicate meshblock boundaries.Within each panel, four dots bordered in black can be observed: the lower two dots correspond to the top meshblock's intra-panel ghost zones and the upper two to the bottom meshblock's.Values in the open green circles are determined through interpolation from the pink dots column by column from left to right.

Figure 6 .
Figure 6.The flow height field at day 5 (left), and the convergence of the peak flow height and maximum zonal flow speed at day 5 when we increase our grid resolution (right).

Figure 7 .
Figure7.The average l1 and l2 error relative to high-resolution solution for cases 2 and 6(Williamson et al. 1992).The log-linear fits of the error values are plotted, and the slopes are labeled.

Figure 8 .
Figure 8.The flow height field at day 40 (left column), day 80 (middle column), and day 100 (right column) for resolutions C48 to C192.

Figure 9 .
Figure 9.The mean flow fields averaged both zonally and temporally over day 200-1200 for the Held-Suarez test.Note that the y-axis is selected to be the normalized pressure scale σ = p/p0.

Figure 10 .
Figure 10.The stream function of Held-Suarez test of the averaged flow field.

Figure 12 .
Figure 12.The stream function of shallow hot Jupiter test of the averaged flow field.

Figure 13 .
Figure 13.The flow height field at day 5 (left), and the convergence of the peak flow height and maximum zonal flow speed at day 5 when we increase our grid resolution (right).
x and y are as previously defined x = tan ξ and y = tan η.Due to the loss of symmetry in an oblate spheroid, the values of δ(n) depend on the panel number n.In polar panels (top and bottom), δ(1, 5) = (x) 2 + (y)

Table 1 .
Local to global coordinate system transformation