Generation of orthogonal curvilinear grids on the sphere surface based on Laplace-Beltrami equations

In this paper, Laplace-Beltrami equations are used to generate orthogonal curvilinear grids on the sphere for ocean models. In addition to overcoming the pole-problem, the grid configuration has quasi-uniform cell-size on the whole sphere. Some quantities such as the grid length along two directions, the angle deviation from orthogonality, the area of the cell to evaluate the quality of the grid, which demonstrate the grid produced is fit to be a model grid on which the finite difference method or finite volume method can be implemented for numerical simulating of global atmosphere and ocean dynamics on large scale.


Introduction
The global climate is getting warmer, causing an increasing threat from extreme weather. It's necessary to forecast the severe change of global climate in order to make possible preparations to reduce the destruction. All this depends on the well development of the global climate model, in which the most important two components are the global atmosphere and ocean model. So far, the spherical grids used can be classified into two types [1] : one type is the structured grid, which is characterized by a regular layout, and the number of first-order neighbors of each element is constant, such as the latitude-longitude spherical grid, cubed-sphere grid, the Yin-Yang overset grid. The other is unstructured grid, which is composed by triangles or irregular polygons, the neighbors of each element are not so clear, such as the sphere icosahedral grid or geodesic grid.
The orthogonal curvilinear grid belongs to a structured grid, it has two major merits: at first, the code of mesh and neighbor-searching is obvious, and the dynamic model equations have less terms related to mesh metrics under the orthogonal curvilinear coordinate. Secondly, the finite difference method or the finite volume method can be implemented with much easy.
At present, many operational weather and climate models are based on a latitude-longitude spherical grid, however the convergence of the meridians leads to resolution clustering at the poles, which causes the well-known pole-problem [2] . Meanwhile, some methods such as zonal filtering or space-smoothing, semi-implicit or Semi-Lagrangian method have been developed to overcome this problem in the numerical simulation of global atmosphere dynamics, however it is difficult to use these methods in the ocean models [3] , therefore it's necessary to consider a new kind of orthogonal curvilinear spherical grid, which is fit for numerical modelling of large-scale atmosphere and ocean dynamics and expression of geographical space information.

Gnomonic projection
Consider the cube with sides of length 2a inscribed into a sphere of radius R such that the eight vertices of the cube exactly touch the sphere and The cube is oriented in such a way that the 3D absolute Cartesian coordinate axes ( , , ) X Y Z are normal to the faces, and the centroid of the cube is at the origin, as it can be seen in Figure 2  There are two kinds of methods to generate spherical quadrilateral grids [4]  Equidistant gnomonic mapping: The relation between the 3D Cartesian coordinates on the sphere and the local 2D Cartesian coordinate on cube surface X a = , see figure 2.2, which shows the relation between a point ( , ) x y in local 2D Cartesian coordinates and a point ( , , ) X Y Z on the sphere with radius R: where r is the Euclidean distance between ( , , ) a x y and the origin. r a x y = + +  Equiangular gnomonic mapping: Rancic′  et al [5] showed that the equiangular projection produces a more uniformly spaced grid compared to the equidistant projection and is more suitable for finite difference method. D.Nair [4] indicated that on the sphere, with the same number of grid points, an equiangular projection is more accurate than an equidistant projection. But these two grids are not orthogonal and dynamic equations would become more complex under these non-orthogonal curvilinear coordinates.

Elliptic grid generation on surface
Elliptic grid generation is the iterative relaxation of a first-guess grid (via successive over relaxation) to satisfy a quasi-linear elliptic system of partial differential equations (PDEs) while imposing smoothness and orthogonality constraints [6] . In the case of the grid on sphere surface, each face of the physical grid in spherical coordinates is mapped to a 2-dimensional rectangular u-v parametric space. The elliptic systems are solved in parametric space while imposing the desired physical constraints; orthogonality and uniform grid spacing are the most important constraints in most of practical application. There are two options for imposing orthogonality constraints: the Neumann constraint that allows the grid points to slide along the boundary of each face, and the Dirichlet constraint that holds the edge points fixed while adjusting the interior points and maintaining a quasi-uniform grid cell spacing. The Neumann approach computes a shift for each grid point along the boundary of the grid by forcing the dot product of the physical variables on two orthogonal coordinates (ξ, η) partial derivatives to be equal to zero; this approach will produce a convergence of grid cells near the corners. Orthogonality may also be imposed through the use of orthogonal control functions, i.e. the Dirichlet approach. These control functions are evaluated at the boundary to impose orthogonality, while maintaining the original grid spacing along the boundaries. The computation of two-order normal derivatives in control functions depends on a layer of ghost boundary points, getting the positions of these ghost points is the key. Then transfinite interpolation can be applied to interpolate the control functions from the boundary to the interior of the grid. In addition, a blending function can also be applied to force the interior points to remain close to the original grid spacing. In classical differential geometry, a surface S is viewed as a mapping from R 2 to R 3 .Consequently, parametric surfaces or physical variables are defined in terms of parametric variables. In grid generation, parametric variables are defined in terms of computational variables, that is: In many cases, the conformal mapping can be used to generate the Body-Fitted mesh on a simple-connected domain under 2D cartesian coordinates. Similarly, a conforming mapping of a smooth bounded surface onto a rectangular region can be constructed by establishing a mapping from a square region of the computational plane into the surface which is orthogonal and has a constant aspect ratio [7,8] . The conditions can be expressed by the system of equations: where M is the grid aspect ratio. These two equations can be rewritten as: Using the chain rule for differentiation, the above equations can be expressed in the form:

Mu av bu
Mv bv cu Assuming the relation between the computational variables , ξ η and the parametric variables , u v is: Using these quantities in equation (5) so that the parametric variables become the independent variables, the system can be expressed in the form solving the above systems for note that 2 1 b ac = − is used. This first-order system is analogous to Beltrami's system of The Beltramians 2 u ∆ and 2 v ∆ are given by: The above system is the basis for elliptic methods generating surface grids. The source terms (or control functions) Φ and Ψ are added to allow control over the distribution of the grid points on the surface. Typically, the points in the computational space are given and the points in the parametric space must be computed. Therefore, it is convenient to interchange variables so that the computational variables ξ and η are the independent variables. The transformation given by equations (10) and (11) is reduced to the following system of equations for which the parametric variables u and v are the solutions of the following quasi-linear elliptic system [8,9] In addition, the following equations could be computed beforehand to get the Beltramians 2 u ∆ and 2 v ∆ we take the finite difference method to solve the equations (13) and (14), and the terms 11 12 22 , , g g g can be computed according to the expression of the parametric surface beforehand [10] . In terms of the grid points in the interior of the region, the central difference is used to approximate the first-order and second-order derivatives: ( ) 1 1, 1, , the above systems can be solved by the iterative method known as successive over-relaxation(SOR).
where the acceleration parameter , i j w should satisfy:

Neumann orthogonal boundary condition
We require the condition of orthogonality in physical space:

The metrics of grid quality
We choose four major quantities used to analyze the grid quality: the grid length alongξ direction and η direction respectively ( h ξ and h η ), the cell-area of grid(area) and the angle deviation from orthogonality (DO), the ratio of h ξ to h η (aspect ratio ).

Numerical results
On the cube surface, the top panel Z a = was chosen as a basis, on which a local 2D Cartesian coordinates ( , ) x y can be built, then stretch each points on the panel to its circum-sphere, i.e.
where R is the radius of the sphere, Here we take 1 a = , initial grid points distribution was chosen as equi-distance gnomonic cubed sphere grid, solving the elliptic equations (13) and (14) with Neumann orthogonal boundary condition equation (18), the convergence could be achieved after about 500 times iterations, the grid in parametric space with cell-number 44*44 was given in Figure 2.3, where the left sub-figure shows the initial grid points distribution in parametric space ( , ) u v ,  the right one shows the final grid points distribution in parametric space after numerical iteration. At present, the cubed sphere grids could be also generated by other two methods, including equi-angular or equi-distance gnomonic projection and conforming mapping [5] , moreover sphere grids produced by equi-angular or equi-distance gnomonic projection are not orthogonal, conforming mapping and elliptic cubed sphere grids are orthogonal. Comparison results about the grid metrics on cubed sphere grids between the three methods were given below.
We get the results under the condition that the grid cell-number respectively. Table 2.1 gives the maximum and minimum grid length, the maximum and average grid cell angle deviation from orthogonality excluding the four corners of the panel, the grid aspect ratio with different space resolutions on equi-angular gnomonic grids; Table  2.2, 2.3 give the same information but on grids produced by conformal mapping and elliptic-smoothing methods respectively.   Once the sphere grid on the top panel was produced, we get the whole sphere grids by proper solid-body rotation. Figure 2.4 shows the grid cell-angle deviation from orthogonality and the cell area variations on equi-angular sphere grids with cell number 44*44*6; while Figure 2.5, Figure 2.6 demonstrate the same information on the sphere grids generated by conformal mapping and elliptic equations iterations method respectively. From Figure 2.5, Figure 2.6, conformal mapping and elliptic-smoothing method produce the much like results, it can be seen that most of grid cell-angle deviation from orthogonality are not more than 10  with the exception of these in the neighborhood of the eight vertices. Even the conformal mapping method, the grids adjacent to the eight corners are also not exactly orthogonal. Under the same condition, the equi-angular gnomonic projection produced grids of the most uniform size, but are non-orthogonal, which was shown in Figure 2.4; while the conformal and elliptic-smoothing method generate quasi-orthogonal curvilinear grids on the cubed sphere at the cost of gridline weak-convergence on eight corner points.

Conclusion
Basic properties of conformal mapping are utilized to improve methods of generating grids rom the solution of elliptic partial differential systems. Using the elliptic systems, A new kind of orthogonal curvilinear grid on the sphere was generated and the grid quality was evaluated according to four major quantities, including the grid length, cell-area, angle deviation from orthogonality and grid aspect ratio, moreover comparisons were made with the cubed sphere grids generated by equi-angular gnomonic projection and conformal mapping, the equi-angular gnomonic projection produced grids of the most uniform size, but are non-orthogonal; while the conformal and elliptic-smoothing method generate quasi -orthogonal curvilinear grids on the cubed sphere at the cost of gridline weak-convergence on eight corner points. This small defect would be eliminated if implicit method on time integration, which indicates that the orthogonal curvilinear spherical grid is satisfying wholly except for the neighborhood of eight vertices, so it can be chosen as a model grid on which the finite difference method or finite volume method can be implemented for numerical simulating of global atmosphere and ocean dynamics on large scale.