On the use of Schwarz–Christoffel conformal mappings to the grid generation for global ocean models

1Ministry of Education Key Laboratory for Earth System Modeling, Center for Earth System Science (CESS), Tsinghua University, Beijing, China 2State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics (LASG), Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing, China 3Department of Atmospheric and Environmental Sciences, University at Albany, State University of New York, Albany, NY 12222, USA


Introduction
The generation of the model grid preludes the simulation with ocean general circulation models (OGCMs) and sea ice models.The majority of OGCMs use orthogonal curvilinear dipolar or tripolar grids with the North Pole relocated to continental areas, including models in the Coupled Model Intercomparison Project, fifth phase (CMIP5; CMI, 2014), and those used for high-resolution oceanic forecast (e.g., Metzger et al., 2014;Storkey et al., 2010).The grids are usually generated with an analytical formulation, such as a Moebius transformation, or orthogonal curves, as in Murray (1996) and Murray and Reason (2001).The only land-sea distribution information that is exploited in these grids is to which continental positions the North Pole is relocated.
Contrary to OGCMs, regional ocean models (Shchepetkin and McWilliams, 2005) usually utilize orthogonal curvilinear grids that only enclose the modeling region.The design choice to align grid lines with coastlines is a well-established practice.Examples include the study of small-scale phenomena such as river plumes (Gan et al., 2009), as well as basinscale modeling such as Xu and Oey (2011).The motivations are that (1) it is beneficiary to align grid lines with coastlines or isobaths, for better simulation of the river discharge and more realistic topographic forcing on the oceanic flow; and (2) the removal of lands in the grid's domain results in lower computational overhead, since lands no longer occupy the logically rectangular index space of the grid.Grid generation tools, such as SeaGrid (Sea, 2014), utilize the numerical solution of Laplacian equations and conformal mappings to Published by Copernicus Publications on behalf of the European Geosciences Union.
ensure the orthogonality of the resulting grid.With SeaGrid, modelers need to manually specify the key vertices of the grid for the modeling region.
As an important trend for global ocean modeling, highresolution simulation has been applied in oceanic forecast (Metzger et al., 2014) and climate studies (Dennis et al., 2012).Up to 0.1 • or higher nominal resolution has been successfully utilized by using large-scale parallel computers.With the high resolution, small-scale phenomena can be resolved explicitly, such as narrow but important water channels and mesoscale eddies, but the simulation is usually very computationally costly.To alleviate this problem, ocean modelers adopt load balancing algorithms to improve computational efficiency, and multi-scale modeling with spatial refinements.The load balancing algorithms exploit the fact that many grid points are inactive during the simulation as a result of the land-sea distribution (Kerbyson and Jones, 2005).On the other hand, multi-scale simulation could exploit the spatial and temporal multi-scale characteristics of the ocean's dynamics, or the practical requirements of the spatial resolution (Mandli and Dawson, 2014;Ringler et al., 2013;Chelton et al., 1998).
In this paper, we propose two new grid generation algorithms that improve existing grids in various aspects.They are based on Schwarz-Christoffel conformal mappings.The first algorithm improves dipolar grids, with an enlarged latitudinal-longitudinal (lat-lon) portion of the grid and smooth transition of grid sizes and scaling factors.The second algorithm aims at supporting high-resolution and multiscale modeling, including (1) the removal of a major continental area from the grid, (2) enhanced spatial resolution in coastal regions, and (3) the alignment of grid lines to largescale coastlines.In the following part of this section, we first review the design of orthogonal curvilinear grids for OGCMs in Sect.1.1, and then give a short introduction to Schwarz-Christoffel conformal mappings in Sect.1.2.

Grids for ocean general circulation models
The practice of modeling the ocean involves (1) the selection of an equation set to characterize the thermodynamic and dynamic evolution of the ocean, and (2) numerical treatments of the equation set, including spatial discretization, numerical approximation and time-domain integration.The majority of OGCMs utilize a mixture of finite difference and finite volume spatial and temporal discretization, and assume local orthogonality of the underlying grid.In the traditional latitudelongitude grid, which is a special case for general orthogonal curvilinear grids, meridians converge at the North Pole.This brings several challenges to the ocean modeling: (1) the small grid step sizes in the zonal direction near the North Pole impose a strict limitation on the maximum allowed time step sizes and computational efficiency, especially for highresolution modeling; and (2) there exists large grid step size anisotropy near the North Pole, which negatively affects the effective spatial resolution.To overcome these shortcomings, current OGCMs usually utilize angle-preserving mappings that relocate the grid's pole from the North Pole to one or several continental locations.
Figure 1 shows two examples (GX1, a dipolar grid, and TX0.1, a tripolar grid) for the POP ocean model (POP, 2014) in CESM (CES, 2014).As shown, the North Pole is relocated to lands (in Greenland for GX1, and in Eurasia and North America for TX0.1).Despite the relocated North Pole, the large portion of the grid is still lat-lon, including the low latitudes and the Southern Hemisphere.There are three benefits to the grid being lat-lon.Firstly, the latitudinal and longitudinal directions are characteristic of the large-scale geophysical fluid dynamics, and the grid lines being zonal/meridional could potentially improve the simulation accuracy, especially for models with low spatial resolutions.Secondly, it allows easy analysis of the model output for the majority part of the ocean.Thirdly, the grid variables (such as latitudes and longitudes, Coriolis parameter, cell edge sizes) could be stored by one-dimensional arrays, instead of the two-dimensional arrays for the general case.By reducing the memory footprint during the simulation, the model's computational performance could be improved.In fact, each of GX1 or TX0.1 consists of two patches that complement each other: (1) the southern patch that is lat-lon (i.e., the area with the latitude lower than a certain "turning latitude" φ); and (2) the northern patch (i.e., the area north of φ), which contains a relocated North Pole and is not lat-lon (as in Murray, 1996).On the boundary between the northern and southern patches, abrupt changes in grid edge sizes should be avoided (Roberts et al., 2006).The (local) scaling factor s (defined as the proportion of adjacent cell edge sizes) in the meridional direction should be close to 1 to ensure smooth transition of edge sizes.
To summarize, we outline the design requirements for OGCM grids as follows.They are loosely sorted according to relative importance, starting from more important or classic ones to less important or more modern ones.
2. Relocation of the grid pole to continental locations.The farther the grid poles from the ocean, the better.
3. The scaling factor is close to 1 for the whole grid; i.e., there is no jump of local grid cell sizes.
4. The major part of the grid is still lat-lon.
5. Grid cell size anisotropy should be low.
6.The grid does not induce very small time steps.
7. The grid is indexable as a Cartesian grid.
8. The grid can reduce as many unused grid points (i.e., grid points on land) as possible.9.The grid can support high-resolution and multi-scale modeling.
This list is arguably more comprehensive than that in Roberts et al. (2006), with the first four items coinciding with the whole list in Roberts et al. (2006) and added items reflecting new trends in ocean modeling.Items 6 through 8 are related to computational aspects.In the temporal domain, the maximum allowed time step size is crucial to the model's computational efficiency, which is constrained by the CFL (Courant-Friedrichs-Lewy) condition.Hence, the grid generation method should not introduce grid cell sizes that are too small, as is a shortcoming of traditional lat-lon grids.To utilize parallel computers, OGCMs usually adopt twodimensional domain decomposition in the horizontal spatial domain.A logical Cartesian grid is both intuitive and easy to implement, for both domain decomposition and communication management (item 7).We denote the logical indexation of the model grid as the grid index space.With respect to item 8, since grid cells dedicated to land are inactive for both computation and communication during the simulation, reducing their total number or proportion in the grid index space could improve efficiency and potentially exempt the need for load balancing.Item 9 is related to the aforementioned trends of high-resolution and multi-scale simulation with OGCMs.It is beneficiary that the grid could achieve higher spatial resolution where needed (such as coastal regions, comparatively shallower regions and polar regions for resolving mesoscale eddies) and alignment of coastlines to grid lines (as in regional ocean models).
It is worth noting that finite-element methods-based OGCMs utilize irregular meshes (Chen et al., 2006;Wang et al., 2014;Pain et al., 2005, e.g.), and could potentially achieve multi-scale ocean modeling.However, it is impossible for the majority of existing OGCMs to adopt irregular meshes, which potentially involves the re-formulation of the numerical treatments of the equation set from scratch.In this paper, one of the proposed algorithms aims at providing a basis for the multi-scale modeling with the large body of existing OGCMs, by using Schwarz-Christoffel mappings to achieve land removal and enhanced spatial resolution in coastal regions.

Schwarz-Christoffel conformal mapping
Conformal mappings (Nehari, 1975) are angle-preserving transformations on the complex plane.As the simplest conformal mapping, the Moebius transformation has been used to generate dipolar grids for OGCMs (Murray, 1996).A Moebius transformation z = a * z+b c * z+d is constructed to (1) map the northern patch (under stereographic projection) to the unit disk, and to (2) map the grid's North Pole (a userspecified position, e.g., in Greenland) to the origin of the unit disk.Once the conformal mapping is constructed (i.e., the variables a, b, c and d are computed), a polar coordinate could be constructed on the unit disk and mapped back to form the grid on the northern patch.Since (1) the Moebius transformation is angle-preserving, (2) the polar coordinate is orthogonal, and (3) forward and backward stereographic projections preserve angles, it is guaranteed that the grid for the northern patch is orthogonal.Figure 2 shows the grid constructed by a Moebius transformation for the northern patch of a turning latitude of 20 • N and the conformal center in Greenland (45 • N, 40 • W).Note that the maximum and minimum of scaling factors of the grid are d2/d3 and d1/d3, which are 1.66 and 0.66, respectively, which are far from 1.
For a single-connected region R with an irregular boundary (e.g., not a circle for a uniform turning latitude), instead of a Moebius transformation, a Schwarz-Christoffel (SC) conformal mapping (Driscoll and Trefethen, 2002) is needed for the mapping between a unit circle and R. For a polygon boundary with n vertices, {v i |1 ≤ i ≤ n}, and internal angles of {φ i |1 ≤ i ≤ n}, the SC mapping f from the unit disk to the region enclosed by the polygon could be defined as where A and C are scalar constants, and z i s correspond to v i s under the SC mapping.The values of A, C and z i s could be numerically computed through a construction process (Driscoll and Trefethen, 2002).Figure 3 shows an example in which an SC mapping is generated for a user-defined region with a polygonal boundary, generated with SCTool-Box (SCT, 2012).Vertices v i s of the polygon and the corresponding vertices on the unit circle (z i s) are shown.The polar coordinate on the unit disk is constructed and mapped with f to form an orthogonal grid on the polygonal region (potentially a northern patch).
For an area with a connectivity of n (2 < n < ∞), e.g., a complex plane with n−1 non-intersecting areas, there exist canonical forms that this region could be mapped to (Chapter 7 of Nehari, 1975).Figure 4 shows a sample region with connectivity of 5, and its canonical forms under Schwarz-Christoffel mappings.For a complex plane with polygonal regions as in Fig. 4a, an SC mapping exists that maps this plane to other planes on which the polygonal regions are mapped to (1) circles, as in Fig. 4b, (2) parallel line segments, as in Fig. 4c, (3) radial line segments, as in Fig. 4d, or (4) circular line segments in Fig. 4d.The conformal invariant point, i.e., the point whose position remains unchanged under the SC mapping, is marked out by an asterisk on each complex plane.
Despite the fact that the canonical forms of multipleconnected regions have long been recognized, the construction theory for such a Schwarz-Christoffel mapping between a user-specified polygonal region and its canonical forms is a recent discovery since Delillo et al. (2004).In this paper, we omit the technical details for SC mapping construction as introduced in DeLillo and Kropf (2011) and DeLillo et al. (2013), and focus on their application in grid generation for OGCMs.
Hereby we denote SCSC as the Schwarz-Christoffel mapping for single-connected regions, and MCSC the Schwarz-Christoffel mapping for multiple-connected regions.In Sect.2, we apply SCSC to improve existing dipolar grids.Specifically, we improve the patching scheme by designing a North Polar cap with an irregular boundary instead of a uniform turning latitude, to enlarge the lat-lon portion of the dipolar grid.By using SCToolBox for implementation, we construct an SCSC mapping to generate the orthogonal grid for the North Polar cap.In Sect.3, we treat the continental part of the earth as a multiple-connected region (with irregular boundaries), denoted as C, and remove C by constructing an MCSC mapping to map C to a set of slits.Hence, a grid can be constructed with continents removed and grid lines aligned to continental boundaries.The open-source software of MCSC (MCS, 2013) is used for the construction of the conformal mapping.For each grid generation algorithm, we provide a sample grid for the static evaluation and carry out idealized simulation with the POP model, to validate that they can be easily adopted by OGCMs.In Sect.4, extensions to the proposed methods and the relationship with existing methods are discussed.Finally, Sect. 5 concludes the article.

Pole relocation with SCSC mappings
In this section we apply the Schwarz-Christoffel mapping for single-connected regions (SCSC) to the generation of dipolar grids for OGCMs.We mainly address the following grid design objectives: (1) the enlargement of the lat-lon portion of the grid; (2) the mitigation of scale changes across patch boundaries; and (3) reduced grid cell size anisotropy in polar regions.Traditional dipolar grids (Murray, 1996;Roberts et al., 2006) usually contain (1) a northern patch and a southern patch, divided by a turning latitude, (2) a regular lat-lon grid for the southern patch, and (3) a non-lat-lon grid in the northern patch.To reduce the non-lat-lon portion of the grid (50 % for GX1 and 28 % for TX0.1 as shown in Table 1), we design a new patching scheme, which shrinks it to as low as 6.2 % while maintaining the overall smooth transition of  scaling factors.The patching scheme is shown in Fig. 5, in which the earth's surface is divided into five patches.Central to the grid generation algorithm is a northern patch with an irregular, non-longitudinal boundary.We denote it the North Polar cap patch (NP).It includes a relocated pole in Greenland and a smooth but irregular non-longitudinal boundary with four segments, as listed below.
-Longitudinal segment across the Atlantic Ocean at a certain latitude (at 47 • N), which extends into the Eurasia and North America continents on both ends; -Longitudinal segment across the Pacific Ocean at a certain latitude (across the Bering Strait), which extends into the Eurasian and North American continents on both ends; -Smooth linkage for the first two segments in Eurasia; -Smooth linkage for the first two segments in North America.
The last two segments were constructed to ensure that the overall boundary is smooth.Due to the irregular boundary of NP, an SCSC mapping is constructed to map (1) a unit disk to the stereographic projection of NP and (2) the origin of the unit disk to the prescribed continental position in Greenland (i.e., the location for the grid's pole).
The remaining patches are (1) the southern patch (SP), to cover the middle and high latitudes in the Southern Hemisphere; (2) the equatorial band patch (EBP), to cover lowlatitude areas (with meridional refinement); (3) the North Pacific patch (NPP), covering the area between EBP and NP on the Pacific Ocean and Indian Ocean; and (4) the North Atlantic patch (NAP), covering the leftover area between EBP and NP (which mainly corresponds to the North Atlantic Ocean).The boundaries between NPP and NAP are in the meridional direction and in continental areas, i.e., Eurasia and North America.
Furthermore, we construct the oceanic part of NPP and NAP to be lat-lon (and orthogonal).The non-oceanic areas in NPP and NAP are filled with non-orthogonal grid cells.Since the grid points in these areas are not active in the simulation, it is guaranteed that the loss of orthogonality of these points will not affect simulation results.During the construction of each patch, meridional refinements are introduced where the grid size anisotropy is large.
The grid generation algorithm is outlined as follows.
1. Generation of the NP boundary.
2. Grid generation for NP, by constructing an SCSC mapping from a unit disk to the stereographic projection of NP.An orthogonal polar coordinate is generated for the unit disk and mapped back by the SCSC mapping and backward stereographic projection.The meridional grid cell sizes along the boundary of NP are computed.
3. Generation of the Pacific and North Atlantic basin patches, i.e., NAP and NPP.Linkage between (1) NAP and NPP on the eastern and western boundaries; (2) NAP (or NPP) and NP are also constructed.4. Generation of the equatorial and southern grid patches, i.e., EBP and SP, according to grid cell anisotropy requirements.
5. Assembly of the patches into a global grid.
6. Generation of land and depth masks.
As a consequence of the irregular boundary of NP, the grid cell sizes in the meridional direction along its boundary are not uniform.Because the Atlantic Ocean and the Pacific Ocean are not directly connected between 40 and 70 • N, in order to ensure the same meridional step count, we utilize different meridional grid edge sizes in NAP and NPP to mitigate the difference in latitude range in NAP and NPP.
Before introducing the detailed design, we show a sample grid with nominal 1 • resolution (360 × 306) in Fig. 7, with a detailed view of the North Polar region.One in every five grid points is shown in both directions.As is shown, a large portion of the Pacific and Atlantic oceans is lat-lon.The boundaries between NPP and NAP are not smooth, since the orthogonality in these continental regions is not ensured.

Settings for the NP
The choice for the boundary of NP is a trade-off between several factors: the reduction of the length of the boundary on the Pacific and Atlantic oceans, the smoothness of the scaling factor in both meridional and zonal directions, etc.In the sample grid, we choose (1) on the Pacific side, the longitudinal segment crossing the Bering Strait, i.e, at about 66 • N, from 170 • N to 160 • W; (2) on the Atlantic side, the longitudinal segment at about 48 • N, from 70 to 0 • N; (3) the smooth linkage between the two longitudinal segments on the North American and Eurasian continents, respectively; and (4) the relocation of the grid's pole to Greenland (78 • N, 42 • N).The smoothness is ensured by constructing a cosine function-shaped curve in the latitude-longitude space.Suppose that we need to construct the linkage between a point at (lat 1 , lon 1 ) and another point at (lat 2 , lon 2 ), with smooth linkage to longitudinal lines on both ends.The latitude lat on a specific point on the link at a certain longitude lon could be written as a function of the longitude: The scheme is shown in Fig. 6a.Approaching either end of the link, i.e., (lat 1 , lon 1 ) and (lat 2 , lon 2 ), the line is gradually parallel to and joined with the corresponding longitudinal line.Hence, the overall smoothness is attained.We use a discretized boundary (nominal resolution in the zonal direction) as the polygonal boundary (shown in Fig. 5 by blue lines) to construct the SCSC mapping.The grid lines are orthogonal within NP and perpendicular to the boundary of NP.

Control of grid anisotropy
For dipolar grids, grid cells in the polar regions tend to feature very large anisotropy in cell sizes.Meanwhile, equatorial regions are often modeled with higher meridional resolution for purposes such as higher accuracy in the simulation of tropical waves and ENSO (Griffies et al., 2000).For GX1 (shown in Fig. 1a), the grid cell anisotropy (meridional edge size divided by zonal edge size) in the tropics is about 4.24.
In the proposed grid generation method, we introduce a bespoken threshold value to limit the maximum anisotropy in polar regions.This value is also used for the meridional refinement in EBP.Hence, the maximum anisotropy of the whole grid is kept below this value.For SP, due to the fact that it is purely lat-lon, the latitudes and longitudes of the grid points could be computed as one-dimensional arrays.We start from the lowest latitude (the southern boundary of EBP) and numerically integrate to higher latitudes by gradually decreasing the sizes of latitudinal steps.The gradual decrease in meridional step sizes is designed to ensure that the maximum anisotropy does not increase beyond the predefined threshold.For NP, a similar strategy to gradually reduce the latitudinal steps is used, except that due to the uneven edge sizes for any circle in the zonal direction, the anisotropy of a certain zonal circle of the grid is computed as the average meridional edge size divided by the average zonal edge size on the circle.
One important property of the anisotropy control scheme is that although it increases the number of unknowns, it has a limited effect on the largest allowed time step size T max for the simulation.The value of T max is the global minimum of the local largest allowed time step sizes (T max (i, j ), where i and j are the logical indices of the model grid).As constrained by the local stability condition of gravity wave dispersion, T max (i, j ) is proportional to the value of , where x and y are the local grid cell sizes in the zonal and meridional directions.Therefore, in the polar regions, the value of T max (i, j ) is mainly dominated by x s (i.e., the zonal direction).Since the grid anisotropy in these regions is above 1, the control of the anisotropy (with a threshold of 3 in the sample grid) will induce a limited decrease in T max (i, j ).

Mitigation of scaling factors
In order to maintain the accuracy of finite difference operators, we should keep the local scaling factors close to 1.
To ensure that there are no abrupt scale changes across the boundary of NP, after the grid generation for NP, we compute the cell edge sizes along the two oceanic segments of the NP boundary as the average meridional cell edge size of all the cells in each of the segments.Then, NAP and NPP are treated separately to ensure a smooth transition of meridional edge sizes, starting from the uniform meridional edge size in the south (i.e., their boundaries with EBP), and gradually changing to the meridional edge sizes on its northern boundary with NP.We use a cosine function to construct the meridional edge sizes: (1) on both the southern and northern ends the transition of cell sizes is smooth, and (2) the numerical integration (i.e., the sum) of all the meridional edges equals the latitude difference between EBP and NP.This scheme is shown in Fig. 6b.The count of meridional points for NAP and NPP should be the same, to ensure that the grid is still addressable as a Cartesian grid.Within EBP, since the meridional refinement is adopted, a similar scheme for smooth transition is also used, as shown by the example grid in Fig. 7.
The meridional and zonal grid edge sizes of the sample grid are shown in Fig. 8a and b in grid index space, and the anisotropy (local meridional edge size divided by local zonal edge size) is shown in Fig. 8c.The value of the anisotropy threshold for constructing the sample grid is 3.As is shown, the overall anisotropy is kept under 3.As a result, more grid points are dedicated to the Arctic Ocean, with enhanced resolution to resolve the important passages such as Nares Strait and Lancaster Sound.The effect of scaling factor mitigation is shown in Fig. 8d (with areas of the scaling factor larger than 1.1 marked out by black contours).The North Atlantic Ocean is shown, and the largest scaling factor of the whole grid is present (on the boundary between NAP and NP).As is shown, the scaling factor is kept lower than 1.1 for the largest part of the area, and only exceeds 1.1 at the eastern end of the patch boundary.
In the sample grid, 93.8 % of the global oceanic area is still lat-lon.Compared with traditional dipolar grids (as in Roberts et al., 2006), in which about 75.3, 79.0 and 82.2 % of the earth's oceanic part are covered by a regular latitudelongitude part for the turning angles of 20, 25 and 30 • N, respectively, the proposed dipolar grid achieves a much higher portion of pure lat-lon regions.Actually this proportion is similar to that achieved by a tripolar grid (96 %) with highlatitude transition (HLT, as in Murray, 1996) with a polar patch starting at 66 • N. It is worth noting that there exist large scaling factors in the HLT grid, despite the fact that it features a high lat-lon portion.The advantages of the sample grid over existing dipolar grids are achieved through (1) the irregular boundary that allows more of the oceanic part to be included by the lat-lon portion of the grid, and (2) nonorthogonal grid cells with abrupt scale changes that reside on the continents and do not affect simulation.The proposed grid generation method exploits more land-sea distribution information for the grid construction.

Evaluation with OGCM simulation
Since the sample grid is an orthogonal curvilinear grid, it can be utilized by the majority of OGCMs.We evaluate the sample grid from the following aspects: (1) the static evaluation in terms of the maximum allowed time steps and comparison with existing dipolar grids, and (2) the application of the grid in POP.
With the assumption of an explicit (split) formulation, we show the global map of the maximum allowed count of time steps per simulation hour for the external gravity wave (i.e., the barotropic mode) in Fig. 10a.The local maximum allowed time step count is defined as 3600 s gravity acceleration, h the ocean depth at the cell, and dist defined in Sect.2.2.Similar to traditional dipolar grids (e.g., GX1), the critical area with respect to barotropic time step size is in the Arctic, near Greenland, as the result of (1) small zonal edge sizes and (2) the large external gravity wave speed.The sample grid is comparable to GX1 in terms of the maximum allowed time step count per simulation hour (87.7 vs. 58.5 for GX1, which has a lower zonal resolution of 1.125 • and a shallower Arctic basin).
We further implemented the grid in POP.We generate the grid's depth mask field by interpolation and discretization on 46 depth levels.The depth for each vertical layer starts from 3 m for the surface layer to 250 m in the deep ocean.The global maximum depth is 6000 m.Simulation with idealized forcing is carried out to demonstrate that the sample grid can be easily utilized by OGCMs.The configuration of the simulation is as follows: (1) the potential temperature and salinity are initiated to a climatological profile (POP, 2014); (2) the model is forced by analytical, latitude-dependent wind stress and surface heat forcing (with SST restoring) that is constant in time, shown in Fig. 9a and b, respectively; and (3) the spin-up period is 50 years.In Fig. 10b and c, we show the mean surface temperature (SST) field and sea surface height (SSH) field of the first month after the spin-up, respectively.It is shown that under the analytical forcings, the model could reproduce reasonable SST and SSH distributions.This provides validation that the grid could be integrated as a swapin option in existing OGCMs.The further application of the grid in climate studies serves as future work, including (1) refined boundary information for the grid construction, (2) the long-term spin-up with realistic atmospheric forcings, and (3) the simulation with coupled models.

Grid generation with MCSC mappings
In this section we propose the second grid generation method for OGCMs.It targets the new trends of high-resolution and multi-scale ocean modeling (items 8 and 9 in the list in Sect.1).By utilizing Schwarz-Christoffel conformal mappings for multiple-connected regions (MCSC), we remove major continental masses from the grid by mapping them to slits with no area.Similar to the sample grid in Sect.2, the resulting grid is orthogonal curvilinear, and can be utilized by the majority of OGCMs.It can potentially achieve (1) the removal of major continental masses from the grid index space, (2) higher spatial resolution in coastal regions, and (3) the alignment of the large-scale coastlines to the grid lines.
The outline of the grid generation method is as follows.
1.The manual selection of the polygonal boundary for each continental mass that is to be mapped to slits.The area enclosed by the polygon should be strictly in-land.
2. The construction of an MCSC mapping f between the following two complex planes: (a) the complex plane p e , which is the stereographic projection of the earth, and (b) a complex plane p s on which the aforementioned polygonal continental areas are mapped to slit regions.
3. According to the grid resolution requirements, a polar grid coordinate system is generated on p s and mapped back to p e with f , and further onto the globe through backward stereographic projection to form the model grid.4. The generation of land and depth masks.
The mapping scheme between p e and p s is shown in Fig. 11a.
Because of (1) the orthogonality of the polar coordinates on p s , (2) the angle-preserving property of f , and (3) the properties of the stereographic projection, the orthogonality of the grid is guaranteed.The resulting grid features several unique properties.On p s , the grid points close to the slits are mapped to physical positions close to continental boundaries on p e .On p s , the circular slits are parallel to circular grid lines of the polar coordinate, and the radial slits to radial grid lines.The parallelism between grid lines and slits is kept on p e by f .The closer the grid points are to a slit on p s , the better the alignment of grid lines (in either direction) to the corresponding continental boundary that is attained, as is the result of the conformal mapping.Since slits have zero area on p s , the grid lines that cross any slit on p s will be long, potentially curved lines on p e .Because only the inner boundaries of the continental masses are chosen to be mapped to slits, the grid lines crossing slits and their corresponding grid cells only reside on land, and are inactive during the simulation of the ocean.Hence, the removal of a major continental mass in the grid index space (i.e,I -J space) is achieved.Note that as the lateral boundary of the oceans (hence the lateral boundary condition for the model), the coastlines are still present in the grid.
For the numerical implementation, we use an adapted version of the MCSC open-source software (MCS, 2013) for the generation of f , with code changes and augmentations to accommodate the mixed type of slit maps (i.e., supporting radial and circular slits simultaneously).We evaluate the grid generation method by constructing a sample grid with a set of basic, manually chosen continental areas.In the following part of this section, we cover the details of the grid design, the basic evaluation of the sample grid, and the verification with OGCM simulation.

Continental boundary and slit choices
For the sample grid, we limit the choice of the polygonal boundaries for continents to those enclosed by manually picked points.The number of points per continental mass Radial 64.42 % (i.e., the vertex count for the corresponding polygon) is kept small, so that it is feasible for a manual choosing process.It also ensures that the region enclosed by the polygon is strictly in-land, to guarantee that no oceanic region is mapped to slits.The scheme presented here is basic, and only used to demonstrate the grid generation methodology.More advanced schemes are also possible but beyond the scope of this paper, including (1) a spline-based smooth boundary generated from manually picked in-land points, or (2) automatically retrieved continental boundaries.Each polygonal region, representing a continental mass, is mapped to either a circular or radial slit on p s , as introduced in Sect. 1.The slit type of each area can be specified by the user by subjective judgments based on the shape of the area.For example, in Fig. 11, Eurasia, Africa, North America and South America are mapped to a set of two circular slits (for Eurasia and North America) and two radial slits (for Africa and South America).
We construct the sample grid with the four major continental masses as listed in Table 2.The conformal center for p s and p c is in Greenland (77.5 • N, 41 • N), which serves as the geographic location of the grid's North Pole.The corresponding polygonal areas are shown as blue or red colored patches in Fig. 12a, with blue ones mapped to circular slits and red ones to radial slits on p s .For the Eurasia continent, we omit the Black Sea and the Caspian Sea from the ocean modeling, so they are included in the polygonal area for the Eurasia continent.For North America, Greenland is omitted since it contains the conformal center.Africa is divided into two parts as the result of its special shape: the northern part mapped to a circular slit and the southern part mapped to a radial slit.Along the south-eastern coastline of China, several extra points are added, to demonstrate that the grid lines can be forced to follow coastlines more closely if needed.

Basic evaluation of the sample grid
We focus on two aspects for the basic evaluation of the sample grid: (1) the effect of the continental area removal, and (2) the alignment of grid lines to large-scale coastlines.As shown in Table 2, the proportion of the area mapped to slits on North America is 56.55 % and above 64 % for other continents.It is worth noting that this is the result of the specific choice of polygons for the sample grid.If a high-resolution grid is to be constructed, the choice of in-land points could be refined to fit the coastlines more closely.For the sample grid, the in-land points are chosen to be at least 100 km from any oceanic area.On the earth, these areas (at least 100 km away from any seas) account for 80 % of the total land area.Hence, the removed area is over 80 % of the total removable area with respect to the 100 km threshold.In addition, we do not consider peninsulas and archipelagos, which also lowers the proportion of the continental area removal.Figure 12 shows the global sample grid and details in specific regions.The sample grid is nominal 0.25 • in resolution, and 1 in every 10, 10, 2, 2 and 4 points is plotted in the subfigures, respectively.As shown in Fig. 12b, the grid lines in the North Polar region are well constrained by the two major continental boundaries, in particular Hudson Bay.Since regions such as the Scandinavian Peninsula, the Canadian Arctic Archipelago, Indochina and the Iberian Peninsula are not included, grid lines are not forced to align with the continental boundaries in these regions.It is worth noting that (1) the land-sea distributions in these regions are still present in the grid, providing the necessary lateral boundary condition for the OGCM, and that (2) the resolutions in these regions, especially among archipelagos, are enhanced, since the grid points that were dedicated to lands are now redistributed to enhance the spatial resolution of coastal regions.The alignment of grid lines to the continental boundaries is a collateral result of the alignment of grid lines to the boundaries of the removed polygonal regions.Due to (1) the hand-picked poly-gon boundary only reflecting the large-scale land-sea distribution information and (2) the physical coastlines featuring multi-scale and fractal characteristics, the alignment of grid lines to coastlines of the sample grid is only at the large scale.Figure 12c and d show that coastlines can be forced to align with grid lines by adding more points.
The actual land-sea distribution and grid cell edge sizes in the grid index space are shown in Fig. 13.The slits on p s correspond to one-dimensional arrays (constant-I or constant-J ) in the grid index space.They are marked out by corresponding colors as in Fig. 12.It is shown that major continental masses are removed from the grid index space.Although Fig. 13 is not directly recognizable as global maps, it is still clear which part the oceans correspond to.There are some large cell sizes present in the oceanic areas, especially at the slit tips of each continental area.This could be improved through modifications to the positions of hand-picked in-land points.
Through the static evaluation of the grid, we show that the sample grid satisfies the three major aforementioned features: (1) the removal of major continental masses from the grid index space, (2) the alignment of large-scale coastlines to grid lines, and (3) enhanced spatial resolution in coastal areas.The features are achieved at the same time, as a result of the conformal mapping and its harmonics behavior.

Evaluation with OGCM simulation
We evaluate the sample grid under the same protocol as in Sect.2, i.e., by examination of the maximum allowed time step count per simulation hour, and the simulation result with POP under idealized forcings.The sample grid used for simulation is nominal 0.5 • (720 by 330).
Figure 14a shows that the most critical area in restricting the time step sizes is also the Arctic region close to Greenland, with respect to an explicit (split) model formulation.Figure 14b and c show the mean SST field and SSH field after 50 years of the model spin-up.Similar to the sample grid in Sect.2, the sample grid generated with MCSC mapping yields reasonable SST and SSH distribution.The simulation result provides verification that the grid could be integrated with existing OGCMs as a swap-in option.The further work with the proposed grid generation method includes the grid generation with a set of refined continental boundaries for the MCSC mapping, and long-term ocean simulation with realistic atmospheric forcings.

Related work and discussion
In this article, advanced conformal mapping techniques, namely Schwarz-Christoffel mappings, are used to generate orthogonal grids for ocean general circulation models.These curvilinear orthogonal grids are indexable in a regular Cartesian manner, and can be easily integrated with existing OGCMs that already support orthogonal curvilinear grids, such as POP or NEMO (NEM, 2014).
Spatial discretization with irregular and non-orthogonal grids is adopted by algorithms such as finite-element methods, which are widely used in computational fluid dynamics and structure design.In recent years, these non-structured grids have also been adopted in ocean modeling, such as FVCOM (Chen et al., 2006), MPAS-Ocean (Ringler et al., 2013), FESOM (Wang et al., 2014), ICOM (Pain et al., 2005), etc.Despite certain limitations as listed in Griffies et al. (2010), these models have shown advantages in the context of multi-scale modeling (Chen et al., 2009;Scholz et al., 2013, e.g.,).On the contrary, traditional OGCMs based on structured meshes usually provide multi-scale simulation capabilities through nested grids with one-way forcing or two-way coupling.The proposed grid generation method in Sect. 3 provides a basis for multi-scale simulation by using existing OGCMs.Grid points that were dedicated to lands are redistributed to oceanic areas, with coastal areas with potentially higher density of grid points.By combining with other dynamic or static spatial refinement strategies, modelers can carry out multi-scale ocean simulation with respect to the re- quired spatial resolution, such as Rossby deformation radii for resolving mesoscale eddies.
The sample grid in Sect. 3 does not have a latitudelongitude grid along the Equator.It is worth noting that this is not due to the limitation of the methodology.Extra polygons on p e could be added for the equatorial lines on both the Atlantic and Pacific oceans (corresponding to circular slits on p e ), and the conformal map could be constructed to map these polygons to two circular slits on p s .Hence, the alignment of grid lines to the Equator is achieved.
In Murray (1996), a grid generation method was proposed to utilize the conformal equivalence of the orthogonal grid to an electrostatic field.This approach, denoted the multi-polar algorithm (MP), is able to remove some land area from the grid by placing two points with positive and negative charges on the same continent (e.g., Eurasia, Antarctica).Although different in motivation, this approach and the one proposed in Sect. 3 bear certain similarity in the capability to remove continental areas.With MP, the modeler does not have direct control of the continental area to be removed.For effective removal of continental masses, the modelers need to make extra efforts, including (1) a variable number of charged points, with variable values of charge, and (2) the numerical optimization of these values in order to approximate the boundary of removed areas with the real continental boundaries.Instead, with the proposed method based on MCSC mapping, modelers have direct control over the areas to be removed, by specifying the polygonal areas and the specific slit type for each area.
The invariance of the solution of a Laplacian equation with first-or second-type boundary conditions under conformal mappings could be further utilized for the grid generation for OGCMs.Instead of the construction of the conformal mapping, a numerical solution of a boundary value problem is needed.The first-and second-type boundary conditions are equivalent to the circular and radial slits as in the approach in Sect.3, assuming that the grid's poles in Greenland and the South Pole have the first-type boundary condition.
The proposed method in Sect. 3 could also be applied to the generation of grids for regional ocean models such as ROMS (Shchepetkin and McWilliams, 2005), especially for regions with complex land-sea distribution such as the Canadian Arctic Archipelago or Southeast Asia.For traditional grid generation methods, the multiple connectivity of land areas in these regions poses a special challenge, as well as narrow but hydrologically important water channels.The proposed algorithm in Sect. 3 potentially serves as a framework to overcome both problems.
In the supplements, we provide the grid input files for POP in binary format for the two sample grids, with nominal 1 and 0.5 • resolution, respectively.Integration with other OGCMs requires the conversion of grid input formats and the interpolation of initial conditions for the ocean, etc.The two grid generation methods in this paper are implemented in MATLAB, which utilize open-source softwares of SC-ToolBox and MCSC.We plan to provide the two grid generation methods in the format of open-source software, with which users could specify various grid settings such as different spatial resolutions, etc.The generation of the SC mappings, as compared with analytical forms-based grid generation methods, is much more computationally intensive.One notable shortcoming of MATLAB is the limited computational performance as compared with the C language or FORTRAN-based implementations.Since the grid generation task is only carried out once for a whole set of simulations (the spin-up run or the coupled model simulations, etc.), currently we only focus on the functionality aspects of the methods.The optimization in terms of the computational performance is postponed to future works.

Conclusions
In this paper, we propose two new grid generation methods for global ocean generation circulation models.Contrary to conventional dipolar or tripolar grids based on analytical formulations, these new methods are based on Schwarz-Christoffel (SC) conformal mappings with userdefined boundary information.The first method improves the conventional dipolar grids.With SCSC mappings, we construct an orthogonal North Pole patch with a smooth but irregular southern boundary.By utilizing the disconnectedness of major ocean basins in the mid-latitudes in the Northern Hemisphere, the scaling factors across patch boundaries are kept low.In the sample grid, 93.8 % of the oceanic area is still covered by a regular latitudinal-longitudinal grid, which is higher than conventional dipolar grids and comparable to commonly used tripolar grids.
The second grid generation method aims at more modern topics in ocean modeling: (1) the removal of major continental masses from the global grid, (2) better resolution at coastal regions with the alignment of large-scale coastlines to grid lines, and (3) the support for multi-scale ocean modeling.This is achieved by constructing an MCSC mapping, which maps the continental areas to slit regions.In the grid index space, these areas correspond to one-dimensional grid cells.The conformal mapping ensures the alignment of grid lines with the boundaries of these areas, hence achieving the approximate alignment with coastlines.Oceanic regions near these boundaries feature a higher density of grid points, which corresponds to better spatial resolution in coastal regions.Compared with conventional dipolar or tripolar grids, this method exploits more information on land-sea distribution in the form of user-prescribed continental boundaries.
Through static evaluation and simulation with the POP ocean model, we show that the sample grids can serve as swap-in replacements for existing grids for the majority of OGCMs that already support orthogonal curvilinear grids.The MCSC-based grid generation method could also be used in combination with other dynamic or static spatial refinement methods to achieve multi-scale ocean modeling.For the first aspect of the future work, we plan to further apply the proposed methods with refined continental boundary information for grid construction and long-term spin-up runs with realistic atmospheric forcings.As the second aspect, we plan to formulate the proposed methods into complete, opensource grid generation software for both global and regional ocean modeling.
The Supplement related to this article is available online at doi:10.5194/gmd-8-3471-2015-supplement.

Figure 2 .
Figure 2. Northern patch constructed from a Moebius transformation.

Figure 8 .
Figure 8. Scales of the 1 • global grid with North Pole relocation.

Figure 9 .
Figure 9. Analytical wind stress and surface heat forcing.

Figure 10 .
Figure 10.Model evaluation of sample grids with SCSC mapping.

Figure 11 .
Figure 11.Conformal mapping between p e and p s .
The minimum time step count per simulation hour is 116.1.As compared with the sample grid generated with SCSC mapping (1 • nominal zonal resolution with 87.7 steps per hour) and GX1 (1.125 • nominal zonal resolution with the global maximum depth of 5500 m and 58.5 steps per hour), the sample grid is comparable in terms of the maximum allowed time step size.

Figure 13 .
Figure 13.Meridional and zonal grid step sizes in grid space.

Figure 14 .
Figure 14.Model evaluation of sample grids with MCSC mapping.

Table 1 .
Statistics of the sample orthogonal curvilinear grids.

Table 2 .
Continental boundaries and slit information.