Rederivation of the centroid formulation in a second-order conservative remapping scheme on spherical coordinates

The transformation of data from one grid system to another is common in climate studies. Among the many schemes used for such transformations is second-order conservative remapping. In particular, a second-order conservative remapping scheme to work on the general grids of a sphere, either directly or indirectly, has served as an important base in a variety of studies.

In this study, the author describes a fundamental problem in the derivation of the method proposed by a pioneer study relating to the treatment of the centroid used as a reference point for the second-order terms in the longitudinal direction. In principle, use of the original formulation has a potential to cause damage to the entire remapping result. However, a preprocessing procedure on the longitude coordinate suggested in the algorithm for other objectives tends to minimize or even erase the error as a side effect in many, if not most, typical applications. In this study, an alternative formulation is proposed and tested and is shown to work in a simple application.

1 Introduction

Numerical climate models commonly couple individual component models such as models for atmosphere, ocean, and land. These component models are typically developed as stand-alone models and often adopt their own grid system for efficiency. Coupling between such components involves field transformations of data from one grid system to another, while preserving key attributes of interest, e.g., global and/or local integrals. This procedure for conservative quantities is often referred to as conservative remapping (e.g. Dukowicz and Kodis, 1987). As summarized in Mahadevan et al. (2022), there have been considerable efforts to create conservative remapping algorithms for various problems.

Remapping algorithms used in global climate studies are typically based on first- and second-order conservative mesh-based schemes (Mahadevan et al., 2022). In the first-order conservative scheme, a conservative quantity assuming a constant distribution over the source grid cell is transformed into the overlapped destination grid cells with area-weighted remapping (Bryan et al., 1996). On the other hand, in the second-order conservative scheme, a linear distribution within a source grid cell is assumed, which results in a more accurate and smoother transformation than is the case for first-order schemes. In particular, a second-order algorithm works efficiently when remapping from spatially coarse resolution to fine resolution. Because of this, it is considered the preferred choice in many remapping applications. Dukowicz and Kodis (1987) (hereafter referred to as DK87) first provided a second-order conservative remapping algorithm that works for any general grid system using Gauss's divergence theorem for simplification of area integrals converted into line integrals. According to Taylor (2024), most conservative remapping algorithms are variants of this approach (there is also a good summary of the remapping method in the appendix of Taylor, 2024).

Jones (1999), hereafter referred to as J99, extends the DK87 theory to spherical coordinates (more specifically, geographical coordinates of latitude-longitude), offering an approach that can be applied to any type of grid on a sphere. Many efforts to maximize efficiency are included in the proposed algorithm, and a number of problems essentially originating from the spherical coordinate system are solved. J99 also provides the Spherical Coordinate Remapping and Interpolation Package (SCRIP), a native software to implement the algorithm (see https://github.com/SCRIP-Project, last access: 18 June 2026) in addition to four other remapping methods. SCRIP faithfully follows the algorithm of J99, but for the line integral, it adopts only one of several implementation methods: calculating along straight lines on a latitude-longitude coordinate map, rather than along geodesics on the sphere. Thus, SCRIP is applicable to any grid in the sense of the original algorithm, but it does not necessarily handle the geometric properties of arbitrary spherical grids appropriately. Despite this, due to its ease of use and simplicity, SCRIP is one of the most widely used remapping software packages in the climate community (Ullrich et al., 2009). For example, Climate Data Operators (CDO) (Schulzweida, 2023) have once included a conservative remapping option that incorporates SCRIP with rewriting the source code from Fortran to ANSI/C. Recently, CDO implementations no longer rely on the SCRIP-based schemes not only for the second-order but also for the first-order conservative remapping. Nevertheless, it appears to be used, albeit sparingly, even in recent studies (e.g., de Vries et al., 2024 explicitly mention that they use CDO REMAPCON2 command). In addition, SCRIP has been adopted by the general coupler library OASIS3-MCT_3.0 (Craig et al., 2017), which is used by many modeling groups.

The algorithms and the software proposed in J99 have, either directly or indirectly, been an important base in a variety of studies, including both observational and model data analyses (e.g. Barnes et al., 2024), as well as numerical model development (e.g. Ding et al., 2024). Recently, some softwares use implementations for the second-order conservative remapping which significantly differ from the J99 algorithm (e.g. Kritsikis et al., 2017) and the communities have been switching to the other software. However, as far as the author surveyed, some recent studies still use the J99 scheme for the second-order conservative method, as explicitly mentioned in, e.g., Ding et al. (2024), Chtirkova et al. (2024), Ren and Zhou (2024) and Damseaux et al. (2025).

Despite this widespread acceptance, however, there appears to be one distinct and fundamental problem in the derivation of core equations in J99 (Eq. 10) that, to the author's knowledge, has not previously been recognized nor reported.

The problem is in the treatment of a reference point to evaluate the second-order term in the longitudinal direction. In J99, one of the core equations is, at the very end, transformed into an invalid formulation. If one implements the J99 algorithm following the equations, in particular Eq. (10), as presented, there is a risk that serious damage will be caused to the remapping result.

Although few, if any, studies using the second-order conservative remapping scheme in SCRIP have reported strange or erroneous behavior, this is not because the derivation is valid. Rather, there is a small preprocessing block in the algorithm suggested in J99, that adjusts some of the key variables for possibly other objectives which can mask the fundamental problem as a side effect. With this adjustment, any errors originating from the invalid derivation tend to be minimized. In fact, the errors can be fully canceled when the source grid cell is a simple one, such as a regular latitude-longitude (RLL) rectangle grid.

In Sect. 2.1, a highlight of the basics of second-order conservative remapping methods is described (more detailed derivation is described in the Supplement). In Sect. 2.2, which is the core part of the present study, the fundamental problem in the J99 derivation is identified, and the reasons why the invalid derivation has not heretofore been revealed as a problem are discussed. Furthermore, a proposal for a consistent formulation of the scheme is presented in this section. In Sect. 3, the influence of the inconsistent formulation is demonstrated in simple but practical cases. An experiment showing a sample implementation of the proposed scheme is presented.

The present study has one clear limitation. More recent works (e.g. Ullrich et al., 2009) provide a geometrically exact remapping framework between regular latitude-longitude and the other grids. The proposed scheme is a correction within the J99 and SCRIP framework, rather than an alternative to fully geometric methods. It should be emphasised that the primary motivation of the present study is to demonstrate that the results of most past research can be used without significant issues, rather than proposing a new numerical scheme.

2 Description of the second-order remapping methods

This section describes the basic idea of the second-order conservative remapping scheme of DK87 and its extension to the spherical coordinate system as formulated by J99, with supplementary explanation by Jones (2024b). The original equations and terms are transformed into the formulation shown in J99. For example, the volume integral notation in DK87 is replaced by the surface integral in accordance with J99. Additionally, some new symbols unique to the present paper are introduced for description.

It is worth mentioning that the derivation of extension of J99 and thus this study is performed entirely within the latitude-longitude coordinate chart. Modern conservative remapping frameworks (e.g. Hanke et al., 2016; Ullrich and Taylor, 2015; Ullrich et al., 2016) avoid relying on latitude-longitude formulas and instead express all geometry in 3-D Cartesian coordinates.

2.1 Derivation on a general case

The object is to compute in a conservative manner, a flux term on a destination grid from the flux term on a source grid over a surface of three-dimensional Euclidean space. For any flux terms that must satisfy a constraint to preserve conservation, the flux integral over each source grid cell must be consistent with the average value in the grid cell as follows:

$$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{‾}}{f}}_n{A}_n=\underset{A_n}{\int }{f}_n\mathrm{d}A,\end{array}$$

where n is the source cell index, f_n and ${\stackrel{\mathrm{‾}}{f}}_n$ are a flux term and its average over the area of source cell n, respectively. Equation (1) corresponds to Eq. (19) in DK87. Also, Eq. (1) is identical to Eq. (4) in J99. Hereafter, Eq. (e) in J99 are referred to as Eq. (J99.e) in order to avoid confusion.

DK87 proposes to approximate the source flux by a combination of the average and its gradient, with assuming the flux gradient is constant across a source grid cell locally, as follows:

$$\begin{array}{}\text{(2)}& f_n={\stackrel{\mathrm{‾}}{f}}_n+{\mathrm{\nabla }}_{n}f\cdot \left(\mathit{r}-{\mathit{c}}_n\right),\end{array}$$

where r is the position vector, c_n is the position vector of a reference point (corresponding to ${\overline{\mathit{r}}}_k$ in DK87, Eq. 20) and ∇_nf is a gradient of f in source grid cell n. The reference c_n can be chosen arbitrarily in the source cell; here, it is defined such that the flux approximation Eq. (2) satisfies the condition Eq. (1). By substituting f_n into Eq. (1), given a constant gradient across source grid cell n, the following condition is obtained:

$$\begin{array}{}\text{(3)}& \underset{A_n}{\int }\mathit{r}\mathrm{d}A-\underset{A_n}{\int }{\mathit{c}}_n\mathrm{d}A=\mathrm{0}.\end{array}$$

Equation (3) is the principle condition of the reference c_n term. Conservation is preserved with second-order accuracy if the gradient is at least a first-order approximation.

At least over the three-dimensional Cartesian coordinate system, c_n term can be taken out of the integral. In this case, the reference c_n can be inverted as

$$\begin{array}{}\text{(4)}& {\mathit{c}}_n=\underset{A_n}{\int }\mathit{r}\mathrm{d}A/\underset{A_n}{\int }\mathrm{d}A={\displaystyle \frac{\mathrm{1}}{A_n}}\underset{A_n}{\int }\mathit{r}\mathrm{d}A,\end{array}$$

which is identical to the formulation of ${\overline{\mathit{r}}}_k$ in DK87, and to r_n in Eq. (J99.6). The position computed in Eq. (4) corresponds to the geometric center, often referred to as the centroid, of the source grid cell n under the geometry of the target Euclidean space.

The position provided by Eq. (4) is by definition the geometric mass centroid, that lies inside the sphere when the cell area A_n is on the surface. Chen et al. (2026) summarize a variation of “face centerpoints” and provide formulas for the mass centroid, which help to clarify the distinction between the coordinate-invariant mass centroid of a spherical face and the metric-weighted reference coordinate introduced in this study. This point will be revisited later.

Extension to the spherical coordinates requires to replace the gradient and displacement terms in Eq. (2) in the coordinate system. J99 assumes that the gradient term is fixed with the formulation as follows:

$$\begin{array}{}\text{(5)}& {\mathrm{\nabla }}_{n}f={\left({\displaystyle \frac{\partial f}{\partial \mathit{\theta }}}\right)}_n\widehat{\mathit{\theta }}+{\left({\displaystyle \frac{\mathrm{1}}{\cos \mathit{\theta }}}{\displaystyle \frac{\partial f}{\partial \mathit{\varphi }}}\right)}_n\widehat{\mathit{\varphi }},\end{array}$$

where symbols θ and ϕ are adopted for the latitude and longitude coordinates, respectively (Jones, 2024b). The position vector $\mathit{r}=[x,y,z{]}^{\mathrm{T}}$ of the Cartesian coordinate on the unit sphere is expressed using the spherical coordinate components $(\widehat{\mathit{\theta }},\widehat{\mathit{\varphi }})$ which depend on θ and ϕ. The inner product on the spherical coordinate is not simply a component-wise product as in Cartesian coordinates because the direction of the unit vectors depends on the position. Modern implementations (e.g. Ullrich et al., 2009) formulate Eq. (5) in terms of three-dimensional Euclidean space that makes all inner products coordinate-invariant. More recently, Chen et al. (2026) provide formulas for the mass centroid with boundary-integral expressions for great-circle polygons and correction terms for constant-latitude edges. In contrast, J99 instead maintains the formulation on the spherical coordinate. J99 approximates that the unit vectors are aligned over the source cell such that local orthogonality holds true for a simple approach, where the effect of a radial component ($\widehat{\mathit{\rho }}\mathrm{d}\mathit{\rho }$) is assumed to be small enough to be ignored. Formally, this assumption can be interpreted as the interior centroid being projected onto the surface of the sphere, with the subsequent remapping being derived along the surface.

The formulation of J99 is mainly derived for the area-averaged flux over the destination grid cell (after remapping), and is essentially the same as that over the source grid cell. The flux over the destination grid is formulated as follows:

$$\begin{array}{}\text{(6)}& {\stackrel{\mathrm{‾}}{F}}_k={\displaystyle \frac{\mathrm{1}}{A_k}}\sum _{n=\mathrm{1}}^N\underset{A_{nk}}{\int }{f}_n\mathrm{d}A,\end{array}$$

where ${\stackrel{\mathrm{‾}}{F}}_k$ is the average flux over the destination grid cell k, and A_nk is the area of the source grid cell n covered by the destination grid cell k. The summation is performed for all overlapped cells of N. The average flux term at the destination grid cell can be approximated with using the flux approximation of Eq. (5), as follows:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{F}}_k=& \sum _{n=\mathrm{1}}^N\left[{\stackrel{\mathrm{‾}}{f}}_n{w}_{\mathrm{1}nk}+{\left({\displaystyle \frac{\partial f}{\partial \mathit{\theta }}}\right)}_n{w}_{\mathrm{2}nk}\right.\\ & \left.+{\left({\displaystyle \frac{\mathrm{1}}{\cos \mathit{\theta }}}{\displaystyle \frac{\partial f}{\partial \mathit{\varphi }}}\right)}_n{w}_{\mathrm{3}nk}\right],\end{array}\end{array}$$

which corresponds to Eq. (J99.7). The three coefficients, w_1 nk, w_2 nk, w_3 nk, are called the remapping weights and are derived according to J99 as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }{w}_{\mathrm{1}nk}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\mathrm{d}A,\text{(9)}& {\displaystyle }{w}_{\mathrm{2}nk}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }(\mathit{\theta }-{\mathit{\theta }}_c)\mathrm{d}A,\text{(10)}& {\displaystyle }{w}_{\mathrm{3}nk}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\cos \mathit{\theta }(\mathit{\varphi }-{\mathit{\varphi }}_c)\mathrm{d}A,\end{array}$$

where the coordinates ${\mathit{c}}_n=[{\mathit{\rho }}_c,{\mathit{\theta }}_c,{\mathit{\varphi }}_c{]}^{\mathrm{T}}$, or formally, ${\mathit{c}}_n=[{\mathit{\theta }}_c,{\mathit{\varphi }}_c{]}^{\mathrm{T}}$. The reference point (θ_c,ϕ_c), with implicitly including ρ_c component, is actually called the centroid in J99 (expressed as θ_n,ϕ_n). Note that Eqs. (9) and (10) are presented as intermediate formulations (Eqs. J99.9 and J99.10) during the derivation.

It is reasonable to conclude that Eq. (10) holds for any longitudinal origin; otherwise, the remapping weight w_3 nk would change its value according to the coordinate. Thus, in the computation of the weights for each source cell n, it would be safe to rotate around the pole by ϕ_ofs, which would simply correspond to replacing the longitudinal variable with a relative one. Put formally, Eq. (10) is reformulated into

$$\begin{array}{}\text{(11)}& {\displaystyle }{w}_{\mathrm{3}nk}{\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\cos \mathit{\theta }[\stackrel{\mathrm{̃}}{\mathit{\varphi }}-{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_c]\mathrm{d}\stackrel{\mathrm{̃}}{A},\end{array}$$

where $\stackrel{\mathrm{̃}}{\mathit{\varphi }}=\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{ofs}}$. This is an identity for any ϕ_ofs.

J99 suggests to adopt the source grid cell center as ϕ_ofs for each source cell instead of the globally-fixed longitude origin. The numerical library SCRIP does include this method. This suggestion is raised from the spherical coordinate system nature, where the longitude is multiple valued on one line on the sphere. Such problems can be easily avoided using this simple method.

Actually, the definition of central longitude is ambiguous for general shapes of the grid cell, which must be supplied by the user according to the source grid cell configuration. Since only the difference between the two relative longitudes adjusted by the offset longitude is used in the computation, the central value is of no particular significance. It is even possible to have the offset longitude fall outside the cell boundaries as far as it is enough to avoid the multiple-value longitude issues. This topic will be discussed later.

The final formulations of w_2 nk and w_3 nk conducted in the algorithm are obtained by expanding the reference point (θ_c,ϕ_c). Here, the reference point that corresponds to those defined in J99 is represented as (θ_n,ϕ_n). According to Jones (2024b), the position vectors in Eq. (4) are transformed into the corresponding spherical coordinates with including the metric scale factor, as follows:

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_n={\displaystyle \frac{\mathrm{1}}{A_n}}\underset{A_n}{\int }\mathit{\theta }\mathrm{d}A,\text{(13)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_n\cos \mathit{\theta }={\displaystyle \frac{\mathrm{1}}{A_n}}\underset{A_n}{\int }\mathit{\varphi }\cos \mathit{\theta }\mathrm{d}A.\end{array}$$

Introducing these formulations of the reference coordinate into Eq. (10), the final formulation of w_3 nk is as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }\underset{\mathrm{ORG}}{w_{\mathrm{3}nk}}{\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\mathit{\varphi }\cos \mathit{\theta }\mathrm{d}A-{\displaystyle \frac{w_{\mathrm{1}nk}}{A_n}}\underset{A_n}{\int }\mathit{\varphi }\cos \mathit{\theta }\mathrm{d}A,\end{array}$$

which correspond to Eqs. (J99.9) and (J99.10), respectively. As explained above, Eq. (14) is computed using $\stackrel{\mathrm{̃}}{\mathit{\varphi }}=\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{ofs}}$ term, the longitude relative to a reference longitude ϕ_ofs:

$$\begin{array}{}\text{(15)}& \underset{\mathrm{ORG}}{{\stackrel{\mathrm{̃}}{w}}_{\mathrm{3}nk}}={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\stackrel{\mathrm{̃}}{\mathit{\varphi }}\cos \mathit{\theta }\mathrm{d}\stackrel{\mathrm{̃}}{A}-{\displaystyle \frac{w_{\mathrm{1}nk}}{A_n}}\underset{A_n}{\int }\stackrel{\mathrm{̃}}{\mathit{\varphi }}\cos \mathit{\theta }\mathrm{d}\stackrel{\mathrm{̃}}{A}.\end{array}$$

The step-by-step expansion, e.g., from Eq. (10) to Eq. (14), and the related line integrals via Gauss's divergence theorem are fully detailed in the Supplement.

2.2 Inconsistency in the original formulation

However, the remapping weight for the longitudinal direction determined in Eq. (15) is invalid, which lacks important characteristics of Eqs. (10) and (11). As explained, equations of the remapping weights should hold for any longitudinal origin, which means that the remapping weights $\underset{\mathrm{ORG}}{w_{\mathrm{3}nk}}$ (Eq. 14) and $\underset{\mathrm{ORG}}{{\stackrel{\mathrm{̃}}{w}}_{\mathrm{3}nk}}$ (Eq. 15) must be identical. It is demonstrated below that this characteristic is not guaranteed on these formulations.

Substituting $\stackrel{\mathrm{̃}}{\mathit{\varphi }}=\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{ofs}}$ into Eq. (15), ${\stackrel{\mathrm{̃}}{w}}_{\mathrm{3}nk}$ can be expanded as follows:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}\underset{\mathrm{ORG}}{{\stackrel{\mathrm{̃}}{w}}_{\mathrm{3}nk}}& ={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{ofs}})\cos \mathit{\theta }\mathrm{d}A\\ & -{\displaystyle \frac{w_{\mathrm{1}nk}}{A_n}}\underset{A_n}{\int }(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{ofs}})\cos \mathit{\theta }\mathrm{d}A\\ & =\underset{\mathrm{ORG}}{w_{\mathrm{3}nk}}-{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{ofs}}}{A_k}}\left[\underset{A_{nk}}{\int }\cos \mathit{\theta }\mathrm{d}A-{\displaystyle \frac{A_{nk}}{A_n}}\underset{A_n}{\int }\cos \mathit{\theta }\mathrm{d}A\right].\end{array}\end{array}$$

Therefore, the bracket terms in Eq. (16) must be zero in order to satisfy the condition ${\stackrel{\mathrm{̃}}{w}}_{\mathrm{3}nk}\equiv w_{\mathrm{3}nk}$ for arbitrarily chosen ϕ_ofs. As long as the area element dA is a function of θ such as on the spherical coordinate system, the term cos θ cannot be taken out of the integral, and thus the terms in the bracket are not cancelled. Equation (16) is satisfied only when ϕ_ofs=0, which means that it is definitely inconsistent with its former derivation as Eqs. (10) and (11).

This inconsistency originates from invalid derivation from Eqs. (10) to (14), to substitute the reference longitude following Eq. (13). While transformation of the position vector into the spherical coordinates is conducted on Eq. (4) in J99, the same procedure should be conducted on Eq. (3) instead, because the position vector cannot be extracted from the integral on the spherical coordinate.

A new symbol is introduced to designate such a reference point as ${\mathit{p}}_n=[{\mathit{\theta }}_p,{\mathit{\varphi }}_p{]}^{\mathrm{T}}$, hereafter symbolically referred to as pivot, in order to distinguish it from the centroid described above. Then the condition of the reference coordinate (θ_p,ϕ_p), including the metric scale factor, are formulated instead of Eqs. (12) and (13) as follows:

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }\underset{A_n}{\int }{\mathit{\theta }}_p\mathrm{d}A=\underset{A_n}{\int }\mathit{\theta }\mathrm{d}A,\text{(18)}& {\displaystyle }& {\displaystyle }\underset{A_n}{\int }{\mathit{\varphi }}_p\cos \mathit{\theta }\mathrm{d}A=\underset{A_n}{\int }\mathit{\varphi }\cos \mathit{\theta }\mathrm{d}A.\end{array}$$

The elements θ_p and ϕ_p can be extracted from the integral, and the reference coordinates are formulated as:

$$\begin{array}{}\text{(19)}& {\displaystyle }{\mathit{\theta }}_p& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_n}}\underset{A_n}{\int }\mathit{\theta }\mathrm{d}A,\text{(20)}& {\displaystyle }{\mathit{\varphi }}_p& {\displaystyle }=\left(\underset{A_n}{\int }\mathit{\varphi }\cos \mathit{\theta }\mathrm{d}A\right)/\left(\underset{A_n}{\int }\cos \mathit{\theta }\mathrm{d}A\right).\end{array}$$

While the formulation in the latitudinal direction is identical (Eqs. 12 and 19), that in the longitudinal direction is different regarding the treatment of cos θ term in the denominator (Eqs. 13 and 20).

The remapping weight for the longitudinal direction are reformulated as:

$$\begin{array}{ll}{\displaystyle }{w_{\mathrm{3}nk}}_{\mathrm{P}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }(\mathit{\varphi }-{\mathit{\varphi }}_p)\cos \mathit{\theta }\mathrm{d}A\\ \text{(21)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\mathit{\varphi }\cos \mathit{\theta }\mathrm{d}A-{\displaystyle \frac{\mathrm{1}}{A_k}}{\displaystyle \frac{{\mathrm{\Omega }}_{\mathrm{3}nk}}{{\mathrm{\Omega }}_{\mathrm{3}n}}}\underset{A_n}{\int }\mathit{\varphi }\cos \mathit{\theta }\mathrm{d}A,\text{(22)}& {\displaystyle }{\mathrm{\Omega }}_{\mathrm{3}nk}& {\displaystyle }=\underset{A_{nk}}{\int }\cos \mathit{\theta }\mathrm{d}A.\end{array}$$

Equation (21) is confirmed to hold even with substituting $\stackrel{\mathrm{̃}}{\mathit{\varphi }}=\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{ofs}}$, as follows:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{{\stackrel{\mathrm{̃}}{w}}_{\mathrm{3}nk}}_{\mathrm{P}}& ={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{ofs}})\cos \mathit{\theta }\mathrm{d}A\\ & -{\displaystyle \frac{\mathrm{1}}{A_k}}{\displaystyle \frac{{\mathrm{\Omega }}_{\mathrm{3}nk}}{{\mathrm{\Omega }}_{\mathrm{3}n}}}\underset{A_n}{\int }(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{ofs}})\cos \mathit{\theta }\mathrm{d}A,\\ & ={w_{\mathrm{3}nk}}_{\mathrm{P}}-{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{ofs}}}{A_k}}\left[\underset{A_{nk}}{\int }\cos \mathit{\theta }\mathrm{d}A-{\displaystyle \frac{{\mathrm{\Omega }}_{\mathrm{3}nk}}{{\mathrm{\Omega }}_{\mathrm{3}n}}}\underset{A_n}{\int }\cos \mathit{\theta }\mathrm{d}A\right],\end{array}\end{array}$$

where the bracket terms are cancelled with definition of Eq. (22).

2.3 Influence of the longitude adjustments

How the remapping weights are influenced by the invalid formulation can be demonstrated by using a simple configuration in which both source and destination grids are set as RLL grids on a unit sphere, and the cells are equally spaced along the longitude and latitude. The latitudes and longitudes of the grid lines (cell corner coordinates) are expressed as $\mathit{\theta }=\mathrm{180}\mathrm{°}\left(j-N_{\mathit{\theta }}/\mathrm{2}\right)/N_{\mathit{\theta }},j=\mathrm{0},\mathrm{\cdots },N_{\mathit{\theta }}$; and $\mathit{\varphi }=\mathrm{360}\mathrm{°}\left(i-\mathrm{0.5}-N_{\mathit{\varphi }}/\mathrm{2}\right)/N_{\mathit{\varphi }},i=\mathrm{0},\mathrm{\cdots },N_{\mathit{\varphi }}$, respectively. The source and destination grids adopt $N_{\mathit{\theta }},N_{\mathit{\varphi }}=\mathrm{64},\mathrm{128}$ and $N_{\mathit{\theta }},N_{\mathit{\varphi }}=\mathrm{128},\mathrm{256}$, respectively, where a source cell contains $\mathrm{2}\times \mathrm{2}=\mathrm{4}$ destination cells, and a destination cell does not extend over multiple source cells. Figure 1 shows the distribution of the remapping weight w_3 nk over the example source/destination configuration. One source cell has four remapping weights for each overlapped destination cell; those for the north-west designation cells are plotted in the figure. (It is for this reason that the figure is not symmetric about the equator.) Since the relative orientation of a source cell and its overlapped destination cells is equivalent along the longitudinal direction, the remapping weight must be axisymmetric. Figure 1a displays the results for weights computed with the original formulation Eq. (14), clearly showing the breaking of symmetry. In contrast, in Fig 1b, the remapping weights were computed using the formulation satisfying the pivot condition (Eq. 21) which produces the axisymmetric results shown in the figure.

In Fig. 1, computation of remapping weights is conducted not with the formulations using longitude adjustment for each source cell (Eqs. 15 and 23), but with those using the globally-fixed longitude origin. As described above, J99 suggests to adopt the source grid cell center as reference longitude for each source cell, therefore the breaking of symmetry away from ϕ=0 shown in Fig. 1a is practically not applicable. In fact, this suggestion of adjustment in longitude in the original algorithm minimizes or even erases all the problem as a side effect. Formally, it is possible to substitute ϕ_ofs=ϕ_p in Eq. (16):

$$\begin{array}{}\text{(24)}& \begin{array}{rl}\underset{\mathrm{ORG}}{{\stackrel{\mathrm{̃}}{w}}_{\mathrm{3}nk}}& ={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }(\mathit{\varphi }-{\mathit{\varphi }}_p)\cos \mathit{\theta }\mathrm{d}A\\ & -{\displaystyle \frac{w_{\mathrm{1}nk}}{A_n}}\underset{A_n}{\int }(\mathit{\varphi }-{\mathit{\varphi }}_p)\cos \mathit{\theta }\mathrm{d}A.\end{array}\end{array}$$

Introducing the condition of ϕ_p in Eq. (20), the integral part in the second term of Eq. (24) unconditionally becomes zero. Although the second term is inconsistent overall, it is shown that only the coefficient makes the term inconsistent by comparison between Eqs. (15) and (21). Thus the cancellation of the second term in Eq. (24) provides identical solutions with Eq. (21), even though the coefficients are invalid, since these coefficients are essentially erased by the zero-valued integration term. The problem is, of course, that it does not make sense to expect a valid ϕ_p based on an inconsistent computation of the remapping weights.

However, there is, indeed, an explanation.

Although not forced, it is quite natural to set the offset longitude ϕ_ofs as the center of the longitude range of the source grid cells. One reason for this is that a sample program included in SCRIP makes the computation this way; another is that the center longitude is often used for other situations, e.g., visualization, and thus they can be easily prepared. For some special cases, such as benchmark tests, the central longitude is used to evaluate the flux gradient, which is not generally possible for practical applications. For the RLL rectangle grid cells in spherical coordinates, the center longitude is identical to the pivot longitude, and therefore the rotation helps to cancel the contribution of the pivot term. Moreover, if a cell is symmetric along a meridian, then, naturally, the pivot coordinate coincides with the center longitude. In most cases using various shapes of grid cells, the center longitude defined by the user for particular target grid cells may not be far from the pivot longitude, and the problem of the incorrect contribution of the pivot term can be rendered insignificant, as shown in the previous section. In principle, the offset longitude is left to the user's discretion, and these side effects are generally unexpected.

Figure 1 Demonstration of the remapping weight computation. (a) w_3 nk using the J99 original algorithm (Eq. 14) (b) w_3 nk by Eq. (21).

The recent version of CDO does not include the second-order conservative remapping of the SCRIP equivalent. In addition, according to an old version of CDO reference manual, the second-order conservative remapping command (REMAPCON2) is not available for unstructured source grids. Consequently, CDO users cannot encounter this inconsistency issue when working with unstructured meshes. For structured RLL grids, thanks to the mid-longitude offsetting procedure, there is virtually few risk of users suffering from the inconsistent formulation of the remapping weights. However, the author is not fully convinced, and such a conclusion should be confirmed by an expert in the area.

3 Experiment and discussion

In order to demonstrate the argument of the present paper described in the previous section, a series of sensitivity experiments are performed. The main focus of the present paper is to show how J99 and SCRIP are influenced by inconsistent reference longitudes. The evaluation of them among the other remapping packages is far beyond the scope.

All the remapping experiments are performed using SCRIP version 1.5, with minimum necessary modification relating to the remapping weight computation. The version with proposed modification is hereafter referred to as SCRIP-p to distinguish it from the official SCRIP.

The offset longitude is specified in the external input file in original SCRIP application, which is also followed by SCRIP-p. The longitude adjustment with the externally prescribed offset longitude is left as is, since it is, in any event, necessary in order to deal with the periodic boundary condition in a simple way.

It is worth mentioning that the sensitivity to the offset longitude can be examined only by replacing the values in the input data (variable src_grid_centroid_lon in the input file) Either the source code of the program to compute remapping weights or that to perform remapping can be used without any modification, even for the original SCRIP implementation. Although the test program included in the official SCRIP uses the offset longitude for computing the input source field to remapping, the input field is also prescribed by external files in the present study thus the offset longitude is not used anywhere except for the remapping weight computation.

There may be other problems in J99 and original SCRIP: it is reported that the treatment of parametric form for cell sides in the algorithm results in inaccuracies at intersection computation for general grid systems. In their implementation, intersection of two cell sides is computed using linear parameterization of longitudes and latitudes, which is a source of numerical errors for different edge types (Chen et al., 2026). All experiments in the present paper adopt highly simplified RLL grid systems to avoid such issues. (J99, Jones, 2024b). As noted, the main focus of the present study is not to improve the algorithm, but rather to report how the inconsistencies influence the performance in the past application. Therefore all the program source codes, except for those related to the inconsistent reference longitude, are left as they were.

3.1 Summary of original and corrected formulations

A remedy is introduced in the previous section, in order to preserve consistency during the derivation. The original formulation as well as the alternate formulation are summarised here for reference. The first-order remapping weight, w_1 nk, and the second-order remapping weight in the latitudinal direction, w_2 nk, are identical to those originally presented in J99 but are listed here for completeness.

The proposed method presented below is called “Scheme”, however, it is no more than a correction to the original method. It is not a new algorithm for the second-order conservative remapping, but rather a minor variation of the original algorithm to share most of the equations for remapping weights except for the final formulation.

3.1.1 Scheme N – original (native) method

Scheme N is the implementation of the original J99 formulation, which adopts (θ_n,ϕ_n) formulation (Eqs. 12, 13) as the reference point, with introducing the relative longitude to an offset longitude ϕ_ofs. The flux approximation is formulated as

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{F}}_k& =\sum _{n=\mathrm{1}}^N\left[{\stackrel{\mathrm{‾}}{f}}_n{w_{\mathrm{1}nk}}_{\mathrm{N}}+{\left({\displaystyle \frac{\partial f}{\partial \mathit{\theta }}}\right)}_n{w_{\mathrm{2}nk}}_{\mathrm{N}}\right.\\ & \left.+{\left({\displaystyle \frac{\mathrm{1}}{\cos \mathit{\theta }}}{\displaystyle \frac{\partial f}{\partial \mathit{\varphi }}}\right)}_n{w_{\mathrm{3}nk}}_{\mathrm{N}}\right],\end{array}\end{array}$$

and the corresponding remapping weights are formulated as

$$\begin{array}{}\text{(26)}& \left\{\begin{array}{rl}{\displaystyle }{w_{\mathrm{1}nk}}_{\mathrm{N}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\mathrm{d}A,\\ {\displaystyle }{w_{\mathrm{2}nk}}_{\mathrm{N}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }(\mathit{\theta }-{\mathit{\theta }}_n)\mathrm{d}A\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\mathit{\theta }\mathrm{d}A-{\displaystyle \frac{{w_{\mathrm{1}nk}}_{\mathrm{N}}}{A_n}}\underset{A_n}{\int }\mathit{\theta }\mathrm{d}A,\\ {\displaystyle }{w_{\mathrm{3}nk}}_{\mathrm{N}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }(\stackrel{\mathrm{̃}}{\mathit{\varphi }}-{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_n)\cos \mathit{\theta }\mathrm{d}\stackrel{\mathrm{̃}}{A}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\stackrel{\mathrm{̃}}{\mathit{\varphi }}\cos \mathit{\theta }\mathrm{d}\stackrel{\mathrm{̃}}{A}-{\displaystyle \frac{{w_{\mathrm{1}nk}}_{\mathrm{N}}}{A_k}}\underset{A_n}{\int }\stackrel{\mathrm{̃}}{\mathit{\varphi }}\cos \mathit{\theta }\mathrm{d}\stackrel{\mathrm{̃}}{A},\\ {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\varphi }}& {\displaystyle }=\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{ofs}}.\end{array}\right.\end{array}$$

As described in the previous section, the formulation of Scheme N is valid only when the offset longitude ϕ_ofs equals to the pivot longitude ϕ_p of the source cell n, defined in Eq. (20). In this particular case, the pivot term contributes virtually nothing to the remapping weights.

3.1.2 Scheme P – pivot method

Scheme P is a mostly straightforward implementation of the original J99 formulation, where only the invalid computation of the remapping weight w_3 nk is replaced according to the pivot condition. It adopts (θ_p,ϕ_p) formulation (Eqs. 19, 20) as the reference point. Formally, the centroid definition must be excluded from the beginning of the implementation as it is incompatible with this formulation. The flux approximation is formulated as

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{F}}_k& =\sum _{n=\mathrm{1}}^N\left[{\stackrel{\mathrm{‾}}{f}}_n{w_{\mathrm{1}nk}}_{\mathrm{P}}+{\left({\displaystyle \frac{\partial f}{\partial \mathit{\theta }}}\right)}_n{w_{\mathrm{2}nk}}_{\mathrm{P}}\right.\\ & \left.+{\left({\displaystyle \frac{\mathrm{1}}{\cos \mathit{\theta }}}{\displaystyle \frac{\partial f}{\partial \mathit{\varphi }}}\right)}_n{w_{\mathrm{3}nk}}_{\mathrm{P}}\right],\end{array}\end{array}$$

and the corresponding remapping weights are formulated as

$$\begin{array}{}\text{(28)}& \left\{\begin{array}{rl}{\displaystyle }{w_{\mathrm{1}nk}}_{\mathrm{P}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\mathrm{d}A,\\ {\displaystyle }{w_{\mathrm{2}nk}}_{\mathrm{P}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }(\mathit{\theta }-{\mathit{\theta }}_p)\mathrm{d}A\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\mathit{\theta }\mathrm{d}A-{\displaystyle \frac{{w_{\mathrm{1}nk}}_{\mathrm{P}}}{A_n}}\underset{A_n}{\int }\mathit{\theta }\mathrm{d}A,\\ {\displaystyle }{w_{\mathrm{3}nk}}_{\mathrm{P}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }(\mathit{\varphi }-{\mathit{\varphi }}_p)\cos \mathit{\theta }\mathrm{d}A\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A_k}}\underset{A_{nk}}{\int }\mathit{\varphi }\cos \mathit{\theta }\mathrm{d}A-{\displaystyle \frac{\mathrm{1}}{A_k}}{\displaystyle \frac{{\mathrm{\Omega }}_{\mathrm{3}nk}}{{\mathrm{\Omega }}_{\mathrm{3}n}}}\underset{A_n}{\int }\mathit{\varphi }\cos \mathit{\theta }\mathrm{d}A,\\ {\displaystyle }{\mathrm{\Omega }}_{\mathrm{3}nk}& {\displaystyle }=\underset{A_{nk}}{\int }\cos \mathit{\theta }\mathrm{d}A.\end{array}\right.\end{array}$$

The new term Ω_3 nk to be applied in the evaluation of w_3 nk is introduced. This term is not a remapping weight but is computed with the same procedure as the other three remapping weights. The integral part of Ω_3 nk is computed by transforming it into a line integral using Gauss's divergent theorem following the J99 method for the other integrals, and is formulated as

$$\begin{array}{}\text{(29)}& \underset{A_{nk}}{\int }\cos \mathit{\theta }\mathrm{d}A=\underset{C_{nk}}{\oint }-{\displaystyle \frac{\sin \mathit{\theta }\cos \mathit{\theta }+\mathit{\theta }}{\mathrm{2}}}\mathrm{d}\mathit{\varphi }.\end{array}$$

Although the replacement involves only the computation of a single variable, the source code modification would be the most substantial since the treatment of additional variable Ω_3 nk must be introduced concurrently with the three standard weights.

The formulation of scheme P is consistent for any longitude origin, thus it is not necessary to introduce the offset longitude, while it is still valid with any offset longitudes. Also, Schemes N and P produces the identical solution when the offset longitude ϕ_ofs equals to ϕ_p in Eq. (26).

3.2 Configuration of experiments

In the present study, only the domain of the RLL grid on a unit sphere of N_θ latitudes and N_ϕ longitudes, both for the source and destination grids, is examined. The latitudes and longitudes of the grid lines (cell corner coordinates) are expressed as $\mathit{\theta }=\mathrm{180}\mathrm{°}\left(j-N_{\mathit{\theta }}/\mathrm{2}\right)/N_{\mathit{\theta }},j=\mathrm{0},\mathrm{\cdots },N_{\mathit{\theta }}$; and $\mathit{\varphi }=\mathrm{360}\mathrm{°}\left(i/N_{\mathit{\varphi }}\right)/N_{\mathit{\varphi }},i=\mathrm{0},\mathrm{\cdots },N_{\mathit{\varphi }}$, respectively (this is slightly different from the domain definition used for the demonstration in Fig. 1, which does not influence the discussion). Two series of experiments are performed in this study. The first one is one-time remapping test in order to demonstrate the influences on the standard accuracy measures (described later). The size of the source grid cell is set as $(N_{\mathit{\theta }},N_{\mathit{\varphi }})=(\mathrm{64},\mathrm{128})$. Several destination grid sizes are examined, including $(N_{\mathit{\theta }},N_{\mathit{\varphi }})=(\mathrm{90},\mathrm{180})$, (180,360), (360,720), (720,1440). The second one is 1000-time iterate remapping test to demonstrate the convergence rates. The size of the source grid cell is set to $(N_{\mathit{\theta }},N_{\mathit{\varphi }})=(\mathrm{1024},\mathrm{2048})$, and the destination grid sizes are set as above, roughly following the combination in Mahadevan et al. (2022). Three idealized experiments A1, A2 and A3 were conducted following J99 and Mahadevan et al. (2022). A1 and A2 correspond to AnalyticalFun1 and AnalyticalFun2 presented in Mahadevan et al. (2022), respectively. A2 also corresponds to that of the experiments in presented in J99 whose source field is named as $Y_{\mathrm{2}}^{\mathrm{2}}$. A3 corresponds to another experiment presented in J99, named as $Y_{\mathrm{32}}^{\mathrm{16}}$. A2 and A3 also appear in Lauritzen and Nair (2008); Ullrich et al. (2009).

The source field in experiment A1 is a combination of spherical harmonics functions with frequency wave similar to order 3, given by

$$\begin{array}{}\text{(30)}& \mathit{\psi }=Y_{\mathrm{3}}^{\mathrm{2}}+Y_{\mathrm{3}}^{\mathrm{3}},\end{array}$$

where $Y_m^l$ represents the real spherical harmonic functions evaluated for degree m and polynomial order l.

In experiments A2 and A3, a relatively smooth function resembling a spherical harmonic of order 2 and azimuthal wavenumber 2 (named as $Y_{\mathrm{2}}^{\mathrm{2}}$),

$$\begin{array}{}\text{(31)}& \mathit{\psi }=\mathrm{2}+{\cos }^{\mathrm{2}}\mathit{\theta }\cos \left(\mathrm{2}\mathit{\varphi }\right),\end{array}$$

and a relatively high-frequency wave similar to a spherical harmonic of order 32 and azimuthal wavenumber 16 (named as $Y_{\mathrm{32}}^{\mathrm{16}}$),

$$\begin{array}{}\text{(32)}& \mathit{\psi }=\mathrm{2}+{\sin }^{\mathrm{16}}\left(\mathrm{2}\mathit{\theta }\right)\cos \left(\mathrm{16}\mathit{\varphi }\right),\end{array}$$

are used as input for the source grid in each experiment. The mid-longitude and mid-latitude coordinates for each cell are used as a reference point to compute ψ and its gradient to input.

All the experiments are conducted using the test program included in the official SCRIP package with minimum modification. It is worth mentioning that a special treatment for elements around the poles is implemented in the official package, which is switched off in the present study.

The performance of the conservative remapping algorithms was evaluated using Metrics for Intercomparison of Remapping Algorithms (MIRA) package (Guerra et al., 2021; Mahadevan et al., 2022), with help of TempestRemap (Ullrich and Taylor, 2015; Ullrich et al., 2016) that is conducted to prepare the input fields. Several measures are available by MIRA. Global conservation properties are evaluated using L_g, which corresponds to relative change in the global integral of the scalar field value on the source and the destination grids. The standard accuracy measures, $\Vert E{\Vert }_{L_{\mathrm{2}}}$ and $\Vert E{\Vert }_{L_{\mathrm{\infty }}}$ are presented, which correspond to those used the second-order norm l₂ and the infinity norm l_∞, respectively. A gradient preservation measure $\Vert E{\Vert }_{H_{\mathrm{1}}}$ computed by MIRA package is also presented. The explicit definition of these norms are presented in Mahadevan et al. (2022).

As shown in Sect. 2.3, it is speculated that for RLL rectangular cell cases, the offset longitude virtually works as the pivot longitude (ϕ_p in Eq. 20) in the official SCRIP, which would erase the fundamental problem. For general shapes of grid cells, the offset longitude may not be the same as the pivot longitude. To investigate the sensitivity of this deviation, a simple experiment is presented using the official SCRIP implementation.

In the official implementation, the mid-longitude for each cell is introduced for the offset:

$$\begin{array}{}\text{(33)}& {\mathit{\varphi }}_{\mathrm{ofs}}={\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{0}}+{\mathit{\varphi }}_{\mathrm{1}}}{\mathrm{2}}},\end{array}$$

where ϕ₀, ϕ₁ are the longitude boundaries of the source cell. Using this offset to keep the difference in longitudes within 360°, it is easy to avoid the multiple-value longitude issues. Conversely, the offset can be anywhere as far as it is sufficient to avoid the multiple-value issues.

In order to demonstrate the present paper's argument, three sensitivity experiments are performed: the first one is control case, to adopt mid-longitude (Eq. 33) for each cell. The second one is cell-edge case, in which the offset matches the boundary for each cell (i.e., ϕ_ofs=ϕ₁). The third one is global case, in which the offset longitude is set as ϕ_ofs≡180° for all the source grid cells.

Additional remarks apply to the global offset case. This case is an impractical, idealised configuration, intended merely to demonstrate the insensitivity of remapping to the choice of longitude offsets. Introducing a global offset does not work for general coordinates with a multiple-value longitude issue. For such a remapping configuration, longitude offsetting for each cell is naturally applied, which is incompatible with this global adjustment. In this demonstration, the offset longitude is set to the constant value ϕ_ofs≡180° that is the only way to ensure that the experiment of the present study works correctly. It should never be regarded as a realistic solution for implementing the algorithm described in J99.

Since the pivot longitude matches the mid-longitude of the RLL rectangle cell, the same results should be obtained by Scheme N and P in the control case. The second experiment practically corresponds to an extreme case. It can be naturally expected that the pivot longitude is within the cell for general shapes of the source grid cell, therefore this can be a maximum difference of the pivot and offset longitudes for usual application. The third experiment is more than an extreme case where the equations really hold true while it may be rare for typical SCRIP application.

3.3 Results

Figure 2 shows the results of one-time remapping of experiment A2 using Scheme N, with the destination grid cells as $(N_{\mathit{\theta }},N_{\mathit{\varphi }})=(\mathrm{720},\mathrm{1440})$. Figure 2a is the remapping result (expressed with Ψ) for the mid-offset case that is expected to provide the reasonable solution. Figure 2b corresponds to the result of the edge-offset case. The difference in the remapped fields between the edge-offset and mid-offset cases (d) is comparable in magnitude to the remapping error (b). However, it is smaller than the range of the source field, so the remapped fields of both cases are mostly equivalent (a and c). On the other hand, Fig. 2e shows jagged patterns as the distance from the global longitude offset (180°) increases, which is the result of the global-offset case.

Figure 3 shows the results of one-time remapping of experiment A2 using Scheme P. The results clearly show that the remapped fields remain unchanged regardless of the offset chosen. Qualitatively similar results are obtained for the other destination resolution and for experiments A1 and A3 (not shown). Thus, the formulation of Scheme P is expected as the valid correction of the second-order conservative remapping scheme of J99. It is worth mentioning again that the global offset configuration is impractical, purely to demonstrate the insensitivity. The remapped fields using Scheme N are expected to be sufficiently reasonable given the choice of practical offsets in the longitude.

Figure 2 Results of sensitivity experiments A2 using Schemes N for the offset longitudes where mid, edge, global correspond to the ${\mathit{\varphi }}_{\mathrm{ofs}}=({\mathit{\varphi }}_{\mathrm{0}}+{\mathit{\varphi }}_{\mathrm{1}})/\mathrm{2}$, ϕ_ofs=ϕ₁, and ϕ_ofs=180° cases, respectively. (a) One-time remapped field for mid-offset case (Ψ_mid). (b) Error in the remapped field of the mid-offset case (Ψ_mid−ψ). (c) One-time remapped field for edge-offset case (Ψ_edge). (d) Difference between edge- and mid-offset cases (Ψ_edge−Ψ_mid). (e) Remapped field field for global-offset case (Ψ_global). (f) Difference between global- and mid-offset cases (Ψ_global−Ψ_mid). Resolution of the destination grid cells is $(N_{\mathit{\theta }},N_{\mathit{\varphi }})=(\mathrm{720},\mathrm{1440})$.

Figure 3 The same as Fig. 2 but with the Scheme P is shown. Note that figures (d) and (f) do not represent plotting errors, because the remapped fields are almost identical.

Since the algorithm discussed in the present paper is conservative remapping, it is important to check the errors in the global conservation for all the experiments. Table 1 is the summary of the metric L_g, obtained by sensitivity experiments A2. The results of other two experiments A1 and A3 are in the supplement. All norms are computed using the result of one-time forward remapping from the source grid to four variations of the destination grids. The metrics in the first row in the tables (Scheme N and mid-offset) correspond to those obtained by the official SCRIP. These are reference values of the present study, and evaluation of metrics are examined relative to these values. With floating-point arithmetic of binary64 (specified in IEEE 754-2008 standard, usually referred to as double-precision), we have around 15 significant digits.

Except for the Scheme N, global-offset cases, the global conservation properties after one-time remapping mostly agreed to one part in the first 15 digits, thus the remapping are conservative to machine accuracy, which is the same conclusion as presented in J99. For the Scheme N, global-offset cases, errors in the global conservation are prominent among the experiments. As far as the multiple-value problem is avoided, the remapping results should be insensitive to the choice of offset longitude. Thus it is confirmed that the formulation of the reference (centroid) term in the original algorithm is invalid and has a potential to damage the important properties. However, the global-offset configuration is practically more than extreme which may never happen in the typical application. Instead, Scheme N, cell-edge cases are regarded as an extreme case. Global conservation obtained by Scheme N, cell edge cases are comparable to mid-edge cases, thus no significant damages on the conservation are expected with the original algorithm.

Mathematically, Scheme P should give identical results by replacing the offset longitude, but it is not presented in the experiments. This is due to the finite-precision arithmetic. For example, difference in longitude is not computed in degree units, but in radian units after degree-to-radian conversion in the original implementation of SCRIP, which may result in slightly non-uniform values. However, Scheme P shows comparable errors in the global conservation even for global-edge cases, which confirms the expected insensitivity on the offset longitude.

Tables 2, 3, 4 are the summary the metric $\Vert E{\Vert }_{L_{\mathrm{2}}}$, $\Vert E{\Vert }_{L_{\mathrm{\infty }}}$, $\Vert E{\Vert }_{H_{\mathrm{1}}}$, respectively, obtained by sensitivity experiments A2, one-time forward remapping. The results of other two experiments A1 and A3 are summarised in the Supplement. In the tables, the first seven digits are shown for comparison. In general, all three experiments show qualitatively similar results.

As shown in the tables, the metrics of cell-edge cases are slightly deviated from those corresponding mid-cell cases. Differences relative to the mid-cell case increase according to the increase in resolution of the destination grid, however, the metric mostly maintains its order of magnitude. Since the cell-edge case is regarded as an extreme for practical application, the field after remapping may remain similar without a significant impact as far as the offset longitude is within a source cell.

It is clearly shown that the metrics of global-offset cases change their order of magnitudes. Differences relative to the mid-cell case increase according to the increase in resolution of the destination grid, and they reach around 1000-times larger value for the highest resolution in the present study. Although the magnitudes of these metrics may be still reasonably small, the significant sensitivity on the metrics to the offset longitudes confirms the statement of the present study, which the formulation of the reference (centroid) term is invalid.

The metrics obtained by Scheme P, mid-cell case show identical results with those obtained by Scheme N as expected, because the pivot longitudes match the mid-cell longitudes. Due to influence from cancellation and rounding-off during the floating point computation, the cell-edge and global-offset cases show different metrics from the mid-cell case. Nevertheless, all of them are preserved significantly better than those of Scheme N. Even in the global-offset cases, magnitude of all the metrics are unaffected. Therefore, the formulation of Scheme P is expected as the valid correction of the second-order conservative remapping scheme of J99.

Table 1 Summary of results of the sensitivity experiments A2 for the offset longitude using Schemes N and P. The second column indicates the offset longitude, where mid, edge, global correspond to the ${\mathit{\varphi }}_{\mathrm{ofs}}=({\mathit{\varphi }}_{\mathrm{0}}+{\mathit{\varphi }}_{\mathrm{1}})/\mathrm{2}$, ϕ_ofs=ϕ₁, and ϕ_ofs=180° cases, respectively. Destination grid sizes are $(N_{\mathit{\theta }},N_{\mathit{\varphi }})=(\mathrm{90},\mathrm{180})$, (180,360), (360,720), (720,1440). The absolute value of relative error in global conservation ($\left|L_{\mathrm{g}}\right|$, see Mahadevan et al., 2022 for its definition) is shown in table.

Table 2 The same as Table 1 but $\Vert E{\Vert }_{L_{\mathrm{2}}}$ of Experiment A2 is shown. The digits which are different from corresponding mid-cell experiments are written in bold. When the first different digit comes after 7 from the decimal point, it is written as “(n d)” where n denotes the largest order to show difference.

Table 3 The same as Table 2 but $\Vert E{\Vert }_{L_{\mathrm{\infty }}}$ of Experiment A2 is shown.

Table 4 The same as Table 2 but $\Vert E{\Vert }_{H_{\mathrm{1}}}$ of Experiment A2 is shown.

In order to quantify the stability of remapping, iterative two-way remapping (i.e., sequence of forward and backward remapping) is conducted. In the case of Scheme N with global-offset, the iterate remapping is extremely unstable, where the metrics explode within the first twenty steps. On the other hand, the Scheme P shows stable behaviour where the difference in metric among three variations of offsets are hardly visible on the plots. This result also supports the argument of the present study that the algorithm should hold for any longitudinal origin. The results of three experiments after 1000-time iterate remapping are summarised in the Supplement.

Figure 4 show the evolution of the metric $\Vert E{\Vert }_{L_{\mathrm{2}}}$, $\Vert E{\Vert }_{L_{\mathrm{\infty }}}$, $\Vert E{\Vert }_{H_{\mathrm{1}}}$, respectively, along iterative remapping obtained by sensitivity experiments A2. The results of other two experiments A1 and A3 are summarised in the Supplement. In general, all three experiments show qualitatively similar results also for the iterative remapping. The evolution of metrics of cell-edge case mostly overlaps those of mid-cell case for Scheme N, except for the high resolution destination grid. In experiment A2, the results of destination grid size as (720,1440) clearly deviates along the iteration steps, due to accumulation of small differences in the metrics. Practically, even for general shapes of source grid cells, the mid-cell longitude would be well close to the pivot longitude, thus no critical impact on the remapping may be expected. The difference in metrics evolution between Scheme N, mid-cell case and Scheme P is minor, whose relative difference is below around 10⁻⁷. It can be concluded that for simple application as RLL where the mid-cell longitude matches the pivot longitude for each source cell, a sufficiently reasonable remapping can be obtained.

Figure 4 Results of the sensitivity experiments A2 for the offset longitude using Schemes N and P. The metrics (a) $\Vert E{\Vert }_{L_{\mathrm{2}}}$ (b) $\Vert E{\Vert }_{L_{\mathrm{\infty }}}$ (c) $\Vert E{\Vert }_{H_{\mathrm{1}}}$ as functions of iteration number of two-way remapping are shown. Destination grid sizes are $(N_{\mathit{\theta }},N_{\mathit{\varphi }})=(\mathrm{90},\mathrm{180})$, (180,360), (360,720), (720,1440). Solid and dashed lines indicate that ${\mathit{\varphi }}_{\mathrm{ofs}}=({\mathit{\varphi }}_{\mathrm{0}}+{\mathit{\varphi }}_{\mathrm{1}})/\mathrm{2}$ and ϕ_ofs=ϕ₁ cases, respectively, and mostly they are overlapped except for some particular cases. Global offset cases (ϕ_ofs=180°) are excluded from the plot, because the metrics explode at early stage for Scheme N, and because the metrics overlaps with other cases for Scheme P.

Finally, the convergence rates for two Schemes of all the experiments are shown. These convergence rates are calculated using the script provided by MIRA Dataset (Mahadevan et al., 2021), with some minor adjustment for the present study. It is a slope obtained by the linear regression of metrics logarithms as a function of the characteristic spatial length of the destination mesh. Following Mahadevan et al. (2022), the spatial length is defined as the inverse of the square root of the number of destination grids in the present study. In Mahadevan et al. (2022) the convergence rates are computed using uniform refinements in both source and destination grids, while in the present study, only the destination grids are changed with keeping the source grid. Also, in the present study both source and destination grids are RLL, which are not discussed in Mahadevan et al. (2022). Thus the convergence rates in the present study may not be directly comparable to those in Mahadevan et al. (2022). Since the focus of the present study is to evaluate the influence of inconsistent reference longitudes on J99 and SCRIP remapping; it is considered to be sufficient just by presenting relative performance of the corrected schemes.

Table 5 summarizes the convergence rates for two Schemes of all the experiments. As noted, only relative comparison is presented here. As presented in the tables and figures above, Scheme N mid-cell and Scheme P both cases show the same convergence rates for all the experiments (Scheme P global offset cases are confirmed to be also the same, which are not shown in the table). Scheme N cell-edge cases mostly show smaller convergence rates than corresponding mid-cell cases, except for that of the metric $\Vert E{\Vert }_{L_{\mathrm{\infty }}}$ for experiment A1, which show a slightly larger value. Since the results of remapping should be independent, in principle, on even a small deviation of offset longitudes, it is not important whether the convergence rate is larger or not than the mid-cell case. Although the influence on the convergence rates may be reasonably small even in the practically extreme condition of offset longitudes, a visible sensitivity on the metrics again confirms the statement of the present study, which the formulation of the reference (centroid) term is invalid.

Table 5 Convergence rates of three metrics using Schemes N and P for experiment (a) A1, (b) A2, (c) A3.

3.4 Additional remarks – applications in past studies

Two reports relating to the second-order conservative remapping scheme based on J99 are worth noting.

Ullrich et al. (2009) present a remapping scheme called Geometrically Exact Conservative Remapping (GECoRe) and show its performance in idealized cases, comparing it with other schemes, including SCRIP. They find that the error measures in GECoRe and SCRIP deviate significantly for the second-order methods, where the former produces results one or two orders of magnitude better than the latter.

Valcke et al. (2022) present yet another intercomparison study using four remapping algorithms, including SCRIP. A few results obtained by SCRIP are analyzed in the paper, which shows no significant deviation from the other three algorithms, at least for the second-order conservative remapping. The result plots show that the misfit by SCRIP is the largest among the four algorithms for one benchmark, while it is far less for the other benchmark.

As discussed and demonstrated above, the inconsistent formulation of the reference longitude (ϕ_n in Eq. 13) has little impact on the remapping as far as the offset longitude is not far from the midpoint of source cells. It is actually introduced for the other objective, i.e., to avoid multiple-valued longitude at computing differences in longitudes, however it really works as a side effect to minimize the inconsistency.

It is possible that the inconsistencies relating to offset longitude formulation has some impact on the past studies, however, the author doubts that they explain the behaviors in the above two studies, because the impact on the results is insignificant even when the extreme offset longitude is specified. Rather, different behaviors compared to the other remapping algorithms, if any, should originate from the other part in the J99 algorithm (e.g., computation of intersection). Since, in principle, the effect of discrepancies in the offset and pivot longitudes is unexpected because the former is not under control of SCRIP, a more detailed exploration of the source code and data is required in order to determine precisely what is happening.

Chen et al. (2026) presented the issues in remapping libraries such as the geometric treatment of edge types and floating-point robustness in edge-edge intersections. All experiments in the present paper adopt highly simplified RLL grid systems to avoid such issues, thereby successfully isolating the sensitivity to the inconsistent reference longitudes. However, it should be noted that when such overlap-construction errors are present in general applications, their impacts on remapping can influence the results more significantly than the issue discussed in this study.

It may be desirable to present the closed-form expression of a set of typical line segment types for the schemes in the present study, following Ullrich et al. (2009) that present a geometrically exact regridding solution between RLL and cubed-sphere grids. It may be also beneficial to include a brief comparison between the present manuscript and the modern approaches. As mentioned, however, the main motivation is to positively support the past studies that adopted J99 algorithm, in which their remappings were largely the same as expected, and Ullrich et al. (2009) already present a detailed comparison with SCRIP. The development of the proposed schemes in the present study including such detailed discussion is left to future studies.

4 Summary and conclusion

In this paper, the second-order conservative remapping method on spherical coordinates proposed by J99 is reformulated in an effort to remove the inconsistencies discovered in the original formulation. A proposal is presented for the valid formulation of the source flux approximation and centroid (or pivot) constraints used to compute the remapping weights. The resulting weights were confirmed to be insensitive to the choice of longitude origin. Until now, the native implementation package of the original algorithm SCRIP has served to mask the inconsistency in the original formulation, as an adjustment to the relative longitude in the SCRIP code has tended to minimize or even erase the problem.

The proposed corrections apply only within the coordinate framework of J99 and do not influence the formulation of the conservative remapping algorithm over other coordinate frameworks. Also, the correction merely made the formulation valid; therefore, the proposed scheme in this study remains applicable only to the same class of structured grids as the original algorithm. Area-weighted averages of the geographical coordinate terms do not correspond to geometric centers for general spherical polygons. Their interpretation becomes ambiguous when cells are skewed, large, or irregular.

Given the adjustment in SCRIP, the author believes that in most practical cases, those using the second-order remapping algorithm in J99 will experience no significant negative impact from the inconsistency problem, especially for cases involving RLL rectangular grid cells. However, it may be prudent for those conducting studies that involve irregularly shaped grid cells or non-modest variable fields and require a high degree of accuracy to review relevant prior studies.

The present study is by no means meant to denigrate past research. To the contrary, the author truly appreciates the contributions of past studies and the accompanying programming packages, which have played an invaluable role in the efforts of the entire climate modeling community. This paper is not intended to discourage but rather to support the validity of past applications. If this were not the case, the author would have sought only to develop a new programming package without suggesting revisions to the native SCRIP package. SCRIP-p, a fork of SCRIP, can serve as a drop-in replacement for the original version, acting as a bridge until an official package revision. It should be recognized, however, that SCRIP-p was examined on a somewhat limited basis and for only a few cases. Although it may not fully resolve the fundamental problem for general cases, it is hoped that it will work well as a first trial.