Abstract

How to reduce the computational error is a key issue in numerical modeling and simulation. The higher the order of the difference scheme, the less the truncation error and the more complicated the computation. For compromise, a simple, three-point accuracy progressive (AP) finite-difference scheme for numerical calculation is proposed. The major features of the AP scheme are three-point, high-order accuracy, and accuracy progressive. The lower-order scheme acts as a “source” term in the higher-order scheme. This treatment keeps three-point schemes with high accuracy. The analytical error estimation shows the sixth-order accuracy that the AP scheme can reach. The Fourier analysis of errors indicates the accuracy improvement from lower-order to higher-order AP schemes. The Princeton Ocean Model (POM) implemented for the Japan/East Sea (JES) is used to evaluate the AP scheme. Consider a horizontally homogeneous and stably stratified JES with realistic topography. Without any forcing, initially motionless ocean will keep motionless forever; that is to say, there is a known solution (V = 0). Any nonzero model velocity can be treated as an error. The stability and accuracy are compared with those of the second-order scheme in a series of calculations of unforced flow in the JES. The three-point sixth-order AP scheme is shown to have error reductions by factors of 10–20 compared to the second-order difference scheme. Due to their three-point grid structure, the AP schemes can be easily applied to current ocean and atmospheric models.

1. Introduction

Most ocean models are hydrostatically balanced. Such a balance sets up a lowest limit, (Δx, Δy)L, for the grid spacings. If the model resolution is finer than that limit, nonhydrostatic ocean models should be used. If the hydrostatic balance is kept unchanged, the grid spacing cannot be artificially small and the reduction of discretization errors within the grid spacing limitation becomes important. A logical approach is to use high-order difference schemes.

High-order schemes can be applied to the σ-coordinate ocean models for error reduction. Let (x∗, y∗, z) denote Cartesian coordinates and (x, y, σ) σ coordinates. The conventional relationships between the two coordinates are

 
formula

where H is the water depth, and z and σ increase vertically upward such that z = σ = 0 at the surface and σ = −1 and z = −H at the bottom. The horizontal pressure gradient becomes a difference between two terms using the σ coordinates,

 
formula

which leads to a so-called hydrostatic inconsistency with a large truncation error at a steep topography (Gary 1973; Haney 1991). Mellor et al. (1994) argue that if the pressure gradient force is computed as a vertical integral of discretized Jacobian of density and z coordinate, the problem of hydrostatic inconsistency can be avoided and the spurious residual pressure gradient is drastically reduced. Furthermore, Song (1998) shows that a modified Jacobian formulation on the discredited pressure gradient term with a specially designed nonlinear weighting procedure drastically reduces the residue errors.

Thus, we have two approaches to reduce errors of the σ-coordinate ocean models: 1) bringing certain symmetries of the continuous forms into the discrete level to ensure cancellations of these terms (e.g., Mellor et al. 1994; Song 1998; Song and Wright 1998), and 2) increasing numerical accuracy (e.g., McCalpin 1994; Chu and Fan 1997, 1998, 1999, 2000). Note that the use of high-order schemes does not conflict with the use of symmetric form. We may combine the two approaches (symmetric forms with high-order schemes) to further diminish the computational errors.

Usually, the higher the accuracy of the scheme, the more the grid points are needed. For example, a C-grid fourth-order ordinary scheme needs four grid points (McCalpin 1994) and a C-grid sixth-order ordinary scheme needs six grid points (Chu and Fan 1997) for computing the first derivative. However, the current ocean models generally use three-point second-order schemes. Conversion of the existing ocean models from second- to high-order (fourth or sixth order) schemes becomes easier if the high-order scheme preserves the three-point grid.

Concurrently, Hirsh (1975) and Kreiss (1975) have proposed Hermitian compact techniques using fewer grid points (three instead of five) at each grid point to solve partial differential equations (PDEs). Later on, as pointed out by Adam (1977), the truncation errors are usually four to six times smaller than the same order noncompact schemes. Since then, much work has been done in developing compact schemes for various applications, such as an implicit compact fourth-order algorithm (Navon and Riphagen 1979), a fourth-order compact difference scheme for nonuniform grids (Goedheer and Potters 1985), and fourth-order and sixth-order compact difference schemes for the stagged grid (Chang and Shirer 1985). Recently, Chu and Fan (1998, 1999, 2000) proposed a three-point sixth-order combined compact difference (CCD) scheme for uniform, nonuniform, and staggered grids. This scheme follows earlier work on use of second derivatives in compact differences (such as Rubin and Khosla 1977),

 
formula

which is referred as the Hermitian formula. Here, f is a dependent variable. The CCD schemes use implicit solvers to calculate various differences (twin-tridiagonal solver) and to solve PDEs (triple-tridiagonal solver). Such a complicated procedure may limit the application of the CCD schemes to the ocean models.

In this study, we propose a simple, three-point accuracy progressive (AP) finite-difference scheme for numerical calculation. The major features of the AP scheme are three point, high-order accuracy, and accuracy progressive. The lower-order scheme acts as a “source” term in the higher-order scheme. This treatment keeps in three-point grid with high accuracy and is explicit in time integration.

2. Progress of accuracy

Let the dependent variable f(x) be defined on the interval, 0 ≤ xL. Use a uniform grid, 0 = x1 < x2 < x3 <  ·  ·  ·  < xN < xN+1 = L with a spacing h = xi+1xi = L/N. Let the dependent variable f(x) at any grid point xi and two neighboring points xi−1 and xi+1 be given by fi, fi−1, and fi+1 and let its derivatives at the two neighboring points xi−1 and xi+1 be given by fi−1, fi−1, … , f(k)i−1, and fi+1, fi+1, … , f(k)i+1.

A new approach is to relate higher-order schemes to lower-order schemes. For example, in the standard Pade (fourth order) scheme,

 
formula

the second-order differences of f′ and f″ at the grid point xi,

 
formula

are the source terms of the standard Pade (fourth order) scheme in the right-hand sides of (1) and (2).

Following this path, we propose a series of three-point AP finite-difference schemes with f′, f″, … , f(k) computed by a lower-order scheme as source terms in higher-order schemes. They have the following features: three-point, accuracy progressive, and only one tridiagonal solver needed. The second-order AP scheme is the same as the central difference scheme, and the fourth-order AP scheme is the same as the standard Pade scheme.

3. AP schemes

a. Standard grid

Let (Δfx, Δ2fx2), (δf/δx, δ2f/δx2), and (Df/Dx, D2f/Dx2) be the second-order, fourth-order, and sixth-order differences for (f′, f″), respectively. The accuracy progress from the second-order scheme into the fourth-order AP scheme is given by

 
formula

Taylor series expansion shows that if α1 = 1/4, β1 = 3/2, α2 = 1/10, and β2 = 6/5, the schemes (5) and (6) have the fourth-order accuracy.

The accuracy progress from the fourth-order AP scheme into the sixth-order AP scheme can be defined by

 
formula

When α1 = 7/16, β1 = 15/8, γ1 = 1/16, α2 = −1/8, β2 = 3, and γ2 = −9/4, the schemes (7) and (8) have the sixth-order accuracy (see appendix).

b. Staggered C grid

The second-order AP scheme on a staggered C grid is given by

 
formula

The fourth-order scheme is given by

 
formula

Taylor series expansion shows that if α1 = 1/22, β1 = 12/11, α2 = 1/10, and β2 = 6/5, the schemes (11) and (12) have the fourth-order accuracy.

The sixth-order AP scheme on a staggered C grid can be defined by

 
formula

When α1 = −7/254, β1 = 120/127, γ1 = −17/254, α2 = −5/94, β2 = −102/47, and γ2 = 144/47, the schemes (13) and (14) have the sixth-order accuracy.

4. Fourier analysis of errors

Fourier analysis of errors is commonly used to evaluate various difference schemes, described extensively in Orzag (1971), Kreiss and Oliger (1972), Lele (1992), and Chu and Fan (1998). As pointed out by Lele (1992), Fourier analysis provides an effective way to quantify the resolution characteristics of differencing approximations.

For the purpose of Fourier analysis, the dependent variable f(x) is assumed to be periodic over the domain [0, L] of the independent variable, that is, f1 = fN+1 and h = L/N. The dependent variable may be decomposed into Fourier series,

 
formula

where i = −1⁠. It is convenient to introduce a scaled wavenumber w = 2πkh/L = 2πk/N, and a scaled coordinate s = x/h. The Fourier modes in terms of these are simply exp(iws). The exact first-order and second-order derivatives of (15) generate a function with exact Fourier coefficients:

 
formula

However, the Fourier coefficients of the derivatives obtained from the differencing scheme might not be the same as the exact Fourier coefficients, that is,

 
formula

where w′ = w′(w) and w″ = w″(w) are the modified wavenumber (both real numbers) for the first-order and second-order differencing. The smaller the difference between the exact and modified wavenumbers, the better the difference scheme.

According to Lele (1992), if the modified wavenumbers of the current generalized difference schemes as

 
formula

the modified wavenumbers will be

 
formula

respectively.

For the second-order AP scheme (central difference scheme), we have

 
formula

for the fourth-order AP scheme (standard Pade' scheme), we have

 
formula

and for the sixth-order AP schemes (7) and (8), we have

 
formula

Among various difference schemes, the modified wavenumbers of the first-order differencing w′ (Fig. 1a) and of the second-order differencing w″ (Fig. 1b) of the AP scheme are closest to the exact wavenumber w. Besides, we computed the errors of the modified wavenumbers for the first-order derivative (Fig. 1c)

 
formula

and for the second-order derivative (Fig. 1d)

 
formula

respectively. We clearly see the accuracy increase of the AP schemes from the second-order to the fourth-order, and from the fourth-order to the sixth-order. Furthermore, the performance of the sixth-order AP scheme is very close to the CCD scheme.

Fig. 1.

Fourier analysis of error for derivative approximation: (a) second-order AP scheme (second-order central scheme), (b) fourth-order AP scheme (standard Pade scheme), (c) sixth-order combined compact scheme, (d) sixth-order AP scheme, and (e) exact differentiation

Fig. 1.

Fourier analysis of error for derivative approximation: (a) second-order AP scheme (second-order central scheme), (b) fourth-order AP scheme (standard Pade scheme), (c) sixth-order combined compact scheme, (d) sixth-order AP scheme, and (e) exact differentiation

In multidimensional problems the phase error of first-order differencing scheme appear in the form of anisotropy:

 
formula

Figure 1e shows polar plots of phase speed anisotropy of various schemes for first derivative approximations. The phase speeds for wavenumbers (magnitude) w/π = 1/50, 5/50, … , 45/50, 50/50 are plotted. Here, we also see that the sixth-order AP scheme shows improvement.

5. Boundary treatment

Due to the progressive accuracy, the AP schemes need special treatment near the boundaries x1 and xN+1. For the standard AP schemes, the derivatives can be computed at (x2, x3, … , xN) using the second-order scheme, at (x3, x4, … , xN−1) using the fourth-order scheme, and at (x4, x5, … , xN−2) using the sixth-order scheme (Fig. 2a). Thus, we use the second-order scheme for the grid points x2 and xN (next to the boundaries), the fourth-order scheme for the grid points x3 and xN−1, and the sixth-order scheme for the grid points (x4, x5, … , xN−2). For the sixth-order AP scheme on the staggered C grid, we have the similar treatment (Fig. 2b).

Fig. 2.

Boundary conditions for AP schemes: (a) standard grid and (b) staggered C grid

Fig. 2.

Boundary conditions for AP schemes: (a) standard grid and (b) staggered C grid

6. Evaluation using the Princeton Ocean Model

a. Known solution

Accuracy of any numerical scheme should be evaluated by either an analytical or known solution. For a realistic ocean model, the analytical solution is hard to find. Consider a horizontally homogeneous and stably stratified ocean with realistic topography. Without forcing, initially motionless ocean will keep motionless forever; that is to say, we have a known solution (V = 0). Any nonzero model velocity can be treated as an error.

We use the Princeton Ocean Model (POM; Blumburg and Mellor 1987), implemented for the Japan/East Sea (JES), to evaluate the AP scheme. The model contains the treatment (symmetric forms for term cancellation) proposed by Mellor et al. (1994).

b. JES geography

The Japan Sea, known as the East Sea in Korea, has a steep bottom topography that makes it a unique semi-enclosed ocean basin overlaid by a pronounced monsoon surface wind. The JES covers an area of 106 km2. It has a maximum depth in excess of 3700 m, and is isolated from open oceans except for small (narrow and shallow) straits, which connect the JES with the North Pacific through the Tsushima/Korean and Tsugaru Straits and with the Okhotsk Sea through the Soya and Tatar Straits. The smoothed JES bathymetry is shown in Fig. 3.

Fig. 3.

Geography and isobaths showing the bottom topography of the JES model

Fig. 3.

Geography and isobaths showing the bottom topography of the JES model

c. Model description

The POM is a time-dependent, σ-coordinate, primitive equation circulation model on a three-dimensional grid that includes realistic topography and a free surface (Blumburg and Mellor 1987). The model contains 94 × 100 × 15 horizontally fixed grid points. The horizontal spacing of 10′ latitude and longitude (approximately 11.54–15.18 km in the zonal direction and 18.53 km in the latitudinal direction) and 15 vertical σ-coordinate levels. The model domain is from 35.0° to 51.5°N, and from 127.0° to 142.5°E. The bottom topography is obtained from the smoothed Naval Oceanographic Office Digital Bathymetry Data Base 5 min by 5 min resolution. The horizontal diffusivities are modeled using the Smagorinsky (1963) form with the coefficient chosen to be 0.2 for this application. No atmospheric forcing is applied to the model.

Closed lateral boundaries, that is, the modeled ocean bordered by land, were defined using a free-slip condition for velocity and a zero-gradient condition for temperature and salinity. No advective or diffusive heat, salt, or velocity fluxes occur through these boundaries. At open boundaries, we use the radiative boundary condition with zero volume transport.

d. Initial conditions

The model was integrated with all three components of velocity (u, υ, w) initially set to zero, and with the same potential density used by Mellor et al. (1994) for investigating the horizontal pressure gradient error,

 
formula

which can be treated as an approximation to the area-averaged, vertical density distribution. Here z is the vertical coordinate, Hρ = 1000 m, and ρ0 = 103 kg m−3.

e. Mode splitting

For computational efficiency, the mode splitting technique (Blumburg and Mellor 1987) is applied with a barotropic time step of 25 s, based on the Courant–Friederichs–Levy (CFL; 1928) computational stability condition and the external wave speed; and a baroclinic time step of 900 s, based on the CFL condition and the internal wave speed.

f. Symmetric cancellation

The POM contains the symmetric cancellation scheme (Mellor et al. 1994). Without the cancellation scheme, we compute the scalar pressure field first by vertical integration of density and then its horizontal gradient via straightforward discretization. With the cancellation scheme, we compute the pressure gradient force as the vertical integral of discretized Jacobian of density and z coordinate. With this treatment, the problem of hydrostatic inconsistency can be avoided and the spurious residual gradient is drastically reduced (Mellor et al. 1994).

7. Horizontal pressure gradient force error

When the density field (23) is used, the horizontal pressure gradient should be zero, theoretically; however, it is usually nonzero in the σ-coordinate system due to the computational error. We may use the horizontal pressure gradient force on the three-dimensional grid cell

 
A = ΔxΔy(HΔσ)(∂p/∂x∗, ∂p/∂y∗)

at day 20 to evaluate the performance of the three schemes in computing the horizontal pressure gradient vector for the JES POM: the second-order scheme without the symmetric cancellation (Fig. 4a), the second-order scheme with the symmetric cancellation (Fig. 4b), and the sixth-order AP scheme with the symmetric cancellation (Fig. 4c) at four different σ levels (1, 5, 9, and 14). Here, Δx, Δy, and Δσ are the grid spacings in the x, y, and σ axes. The σ levels 1 and 14 are the surface and bottom, respectively. We use the cubic-spline interpolation to compute ∂p/∂σ. Using the same second-order difference scheme, we find an order of magnitude reduction in the horizontal pressure gradient error with the symmetric cancellation (Fig. 4b) compared to without the symmetric cancellation (Fig. 4a). Using the same symmetric cancellation (Mellor et al. 1994), we find another order of magnitude reduction in the horizontal pressure gradient error with the sixth-order AP scheme (Fig. 4c) compared to the second-order scheme (Fig. 4b).

Fig. 4.

Horizontal pressure gradient force error |A| (unit: N) at day 20, using (a) the second-order scheme without the symmetric cancellation, (b) the second-order scheme with the symmetric cancellation, and (c) the sixth-order AP scheme with the symmetric cancellation

Fig. 4.

Horizontal pressure gradient force error |A| (unit: N) at day 20, using (a) the second-order scheme without the symmetric cancellation, (b) the second-order scheme with the symmetric cancellation, and (c) the sixth-order AP scheme with the symmetric cancellation

We compute horizontally averaged and maximum values of |A|, over a σ level, to compare the performance of these schemes (Fig. 5). The symmetric cancellation scheme (Mellor et al. 1994) has the capability to drastically increase the accuracy using the same second-order scheme. For example, the horizontally averaged (maximum) values of |A| generally reduce by factors of 5–32 (8–23) with the symmetric cancellation scheme compared to without the symmetric cancellation scheme. Moreover, using the symmetric cancellation we see a further reduction in the horizontally averaged (maximum) values of |A| by factors of 6–9 (4–9) with the sixth-order AP scheme compared to the second-order scheme.

Fig. 5.

Vertical dependence of (a) horizontally averaged |A| and (b) maximum |A| (unit: N) at day 20 using the second-order scheme without the symmetric cancellation (a: dash–dotted), the second-order scheme with the symmetric cancellation (b: solid), and the sixth-order AP scheme with the symmetric cancellation (c: dashed)

Fig. 5.

Vertical dependence of (a) horizontally averaged |A| and (b) maximum |A| (unit: N) at day 20 using the second-order scheme without the symmetric cancellation (a: dash–dotted), the second-order scheme with the symmetric cancellation (b: solid), and the sixth-order AP scheme with the symmetric cancellation (c: dashed)

8. Error estimation using the JES POM

We integrated the diagnostic JES POM with Mellor et al.'s (1994) symmetric cancellation scheme from zero velocity and horizontally integrated density field (23) for 20 days with the second-order scheme and the sixth-order AP scheme, respectively. Without computational error, the model ocean should keep motionless. With the computational error, the model generates false motion. The absolute value of the velocity is taken as the error velocity. The smaller the error velocity, the more accurate the numerical scheme would be.

a. Volume transport streamfunction

Volume transport streamfunction (Ψ) on days 5 and 20 are presented in Fig. 6 for the second-order scheme and in Fig. 7 for the sixth-order AP scheme. All the figures show a multi-eddy structure in coincidence with the steep bottom topography (Fig. 3). The difference between two Ψ isolines indicates the volume transport at that location. The smaller the difference, the less the volume transport. The contour interval of Ψ isolines is five times smaller in Fig. 7 than in Fig. 6. This indicates the error reduction by a factor of 5 using the AP sixth-order scheme compared to using the second-order scheme.

Fig. 6.

Volume transport streamfunction (m3 s−1) using the second-order scheme on (a) day 5 and (b) day 20

Fig. 6.

Volume transport streamfunction (m3 s−1) using the second-order scheme on (a) day 5 and (b) day 20

Fig. 7.

Volume transport streamfunction (m3 s−1) using the sixth-order AP scheme on (a) day 5 and (b) day 20

Fig. 7.

Volume transport streamfunction (m3 s−1) using the sixth-order AP scheme on (a) day 5 and (b) day 20

b. Temporal variations of global and peak error velocities

Two quantities are commonly used for identifying the error evolution: peak and global error velocities. The peak error velocity is the maximum absolute value the spurious velocity generated by the pressure gradient errors. Figure 8a shows the time evolution of the peak error velocity for the first 20 days of integration with the second-order and sixth-order AP schemes. The peak error velocities increase to a maximum value and then decay as the inertial oscillation superimposed into some mean values: 4.80 × 10−3 m s−1 (the second-order scheme), and 0.62 × 10−3 m s−1 (the sixth-order AP scheme). The maximum peak error velocity is 5.38 × 10−3 m s−1 (0.76 × 10−3 m s−1) at day 12 (day 7.5) for the second-order (sixth-order AP) scheme. The steady approach of the peak error velocities to these values for the three schemes indicates the stability of the computation. Furthermore, the peak error velocity for the second-order case is roughly eight times that of the sixth-order AP case.

Fig. 8.

Temporally varying (a) peak and (b) global error velocities for the second-order (solid) and the sixth-order AP (dashed) schemes

Fig. 8.

Temporally varying (a) peak and (b) global error velocities for the second-order (solid) and the sixth-order AP (dashed) schemes

At each time step, we compute the spatially averaged error velocity over the whole JES (called the global error velocity) to represent the overall accuracy of the schemes. Figure 8b shows the time evolution of the global error velocity for the first 20 days of integration with the second-order and the sixth-order AP schemes. It increases with time at a much slower pace using the sixth-order AP scheme than the second-order scheme. On day 20 the global error velocity is 2.02 × 10−4 m s−1 for the second-order scheme and 0.34 × 10−4 m s−1 for the sixth-order AP scheme. The global error velocity for the second-order scheme is roughly six times that of the sixth-order AP scheme.

c. CPU time

We integrate the POM for 20 days on the SGI Origin-200 machine to identify the cost of pressure gradient schemes. The central processing unit (CPU) times are 21.067 h for the second-order scheme and 21.167 h for the sixth-order AP scheme, respectively. There is no noticeable CPU time increase from the second-order scheme to the sixth-order AP scheme.

9. Open boundary influence

a. Coastal geometry

One of the difficult problems in shallow water modeling is the uncertainty of the open boundary condition (OBC). At open boundaries where the numerical grid ends, the fluid motion should be unrestricted. Ideal open boundaries are transparent to motions. Two approaches, local type and inverse type, are available for determining OBC (Chu et al. 1997) with low-order conditions. Before converting any ocean model from a second-order scheme to a sixth-order scheme, it is important to verify if a high-order interior plus low-order boundary conditions (such as open boundary conditions) would degrade the interior solution. Consider a horizontally homogeneous and stably stratified coastal ocean with a longitudinal and straight coastline and three open boundaries in the south (20°N), the north (33.5°N), and the east (x = xE = 320 km). Choose coordinates so that the y axis coincides with the coast, positive x pointing offshore. The coastal water is divided into three parts (Fig. 9a): shelf (0 ≤ xx0), slope (x0 < xx1 = 200 km), and deep water (x1 < xxE).

Fig. 9.

Coastal geometry with open boundaries: (a) horizontal view, (b) cross-coastal view, and (c) three-dimensional view

Fig. 9.

Coastal geometry with open boundaries: (a) horizontal view, (b) cross-coastal view, and (c) three-dimensional view

Here, the water depth (h) is only a function of offshore distance x. At the coast (x = 0), the water has a minimum depth hmin (Fig. 9). An analytical bottom topography is proposed such that shelf and slope are arcs of two circles. The shelf has a smaller radius (r), and the slope has a larger radius (R). The two arcs are connected in such a way that the tangent of the bottom topography, dh/dx, is continuous at the shelf break (x = x0). This requirement is met using the same expanding angle (θ) for both arcs (Fig. 9b). Thus, this bottom topography has 5 degrees of freedom: hmin, hmax, r, R, and θ.

Off-shore distance x and bottom topography h are nondimensionalized by

 
formula

and in turn r and R are nondimensionalized by

 
formula

The bottom topography (Fig. 9) can be represented analytically by

 
formula

As in the previous section, initially motionless coastal ocean with an exponentially stratified density (23) and without forcing will keep motionless forever, that is, V = 0. Any nonzero model velocity can be treated as an error. We use the POM implemented for such a coastal region to test the open boundary influence for the high-order AP scheme. At open boundaries, we use the radiative boundary condition with zero volume transport.

The parameters are chosen in this study by

 
formula

This implies that the width of the shelf break (x0) is 33.33 km.

In the numerical integration, the space is discretized by

 
formula

and barotropic and baroclinic time steps are chosen as 6 and 180 s in order to satisfy the CFL condition.

b. Volume transport streamfunction

Volume transport streamfunction (Ψ) on day 10 is presented in Fig. 10a for the second-order scheme and in Fig. 10b for the sixth-order AP scheme. Both figures show a multi-eddy structure in coincidence with the continental slope (Fig. 9). The difference between two Ψ isolines indicates the volume transport at that location. The smaller the difference, the less the volume transport. The contour interval of Ψ isolines is around 200 times smaller in Fig. 10b than in Fig. 10a. This indicates the error reduction by a factor of 200 using the AP sixth-order scheme compared to using the second-order scheme.

Fig. 10.

Volume transport streamfunction (m3 s−1) on day 10 using (a) the second-order scheme and (b) the sixth-order AP scheme

Fig. 10.

Volume transport streamfunction (m3 s−1) on day 10 using (a) the second-order scheme and (b) the sixth-order AP scheme

c. Horizontal pressure gradient force error

We compute horizontally averaged and maximum values of |A| at day 10 over a σ level, to compare the performance of the sixth-order AP and second-order schemes (Fig. 11). The horizontally averaged (solid) and maximum (dashed) |A| values are quite small at the upper ocean for both the second-order scheme (Fig. 11a) and the sixth-order AP scheme (Fig. 11b). Using the second-order scheme, the errors (|A| values) increase monotonically versus σ, with largest errors next to the ocean bottom (σ = 14): 2943 N for the averaged |A| and 18 582 N for the maximum |A|. Using the sixth-order AP scheme (Fig. 11b), the errors (|A| values) are much smaller than the second-order scheme. The largest errors occur at one level higher than the second-order scheme (σ = 13): 44 N for the averaged |A| and 977 N for the maximum |A|.

Fig. 11.

Vertical dependence of mean (solid) and maximum (dashed) horizontal pressure gradient force errors |A| (unit: N) on day 10 (a) using second-order scheme and (b) using sixth-order AP scheme, (c) Error ratio between second-order scheme vs sixth-order AP scheme. The symmetric cancellation is used in both cases

Fig. 11.

Vertical dependence of mean (solid) and maximum (dashed) horizontal pressure gradient force errors |A| (unit: N) on day 10 (a) using second-order scheme and (b) using sixth-order AP scheme, (c) Error ratio between second-order scheme vs sixth-order AP scheme. The symmetric cancellation is used in both cases

Error ratios between second-order and sixth-order schemes show drastic reductions by factors of 30–514 in the horizontally averaged |A| values (solid) and by factors of 43–1160 in the horizontally maximum |A| values (dashed) (Fig. 11c). The error reduction becomes very evident in the middle level (σ = 6) of the ocean.

d. CPU time

We integrate the POM for this open boundary test case for 10 days on the SGI Origin-200 machine to identify the cost of pressure gradient schemes. The CPU times are 501.07 min for the second-order scheme and 504.19 min for the sixth-order AP scheme, respectively. There is no noticeable CPU time change from the second-order scheme to the sixth-order AP scheme.

10. Conclusions

1) This study shows a process of building a series of accuracy progress difference schemes with the lower-order scheme as the “source” term in the higher-order scheme. Such schemes have the following features: three-point, accuracy progressive, and only tridiagonal solver needed. The second-order AP scheme is the same as the central difference scheme, the fourth-order AP scheme is the same as the standard Pade scheme, and the sixth-order AP scheme has the sixth-order accuracy.

2) Fourier analysis shows the progressive error decrease in modified wavenumber from the second-order AP scheme to the fourth-order scheme, and then from the fourth-order scheme to the sixth-order scheme.

3) The Princeton Ocean Model implemented for the Japan/East Sea is used to demonstrate the benefit of using the sixth-order AP scheme. Two calculations of unforced flow are performed with the second-order scheme and the sixth-order AP scheme, respectively. The results show that the sixth-order AP scheme has error reductions by factors of 5–8 compared to the second-order difference scheme.

4) Using the sixth-order AP scheme for computing the pressure gradient does not require much more CPU time. Taking POM for the Japan/East Sea as an example, the CPU time for the sixth-order AP scheme is almost the same as for the second-order scheme.

5) Both standard and staggered AP schemes were developed. The sixth-order AP scheme on a staggered C grid can be easily implemented into current ocean models.

6) The high-order AP scheme has the capability to handle the open boundary condition. There is no degradation of the interior solution using high-order AP schemes and low-order boundary conditions.

7) Future studies include applying AP schemes to nonuniform grid systems as well as designing even higher-order schemes such as an eighth-order AP scheme.

Acknowledgments

This work was funded by the Office of Naval Research, the Naval Oceanographic Office, and the Naval Postgraduate School.

REFERENCES

REFERENCES
Adam
,
Y.
,
1977
:
Highly accurate compact implicit methods and boundary conditions.
J. Comput. Phys
,
24
,
10
22
.
Blumburg
,
A.
, and
G.
Mellor
,
1987
:
A description of a three-dimensional coastal ocean circulation model.
Three-Dimensional Coastal Ocean Models, N. S. Heaps, Ed., Amer. Geophys. Union, 1–16
.
Chang
,
H-R.
, and
H. N.
Shirer
,
1985
:
Compact spatial differencing techniques in numerical modeling.
Mon. Wea. Rev
,
113
,
409
423
.
Chu
,
P. C.
, and
C.
Fan
,
1997
:
Sixth-order difference scheme for sigma coordinate ocean models.
J. Phys. Oceanogr
,
27
,
2064
2071
.
——, and ——,
.
1998
:
A three-point combined compact difference scheme.
J. Comput. Phys
,
140
,
370
399
.
——, and ——,
.
1999
:
A three-point sixth-order nonuniform combined compact difference scheme.
J. Comput. Phys
,
148
,
663
674
.
——, and ——,
.
2000
:
A staggered three-point combined compact difference scheme.
Math. Comput. Modell
,
32
,
323
340
.
——, ——, and
Ehret
,
L. L.
,
1997
:
Determination of open boundary conditions with an optimization method.
J. Atmos. Oceanic Technol.
,
14
,
723
734
.
Courant
,
R.
,
K. O.
Friedrichs
, and
H.
Levy
,
1928
:
Uber die partiellen differenzengleichungen der mathematischen physik.
Math. Annal
,
100
,
32
74
.
Gary
,
J. M.
,
1973
:
Estimate of truncation error in transformed coordinate primitive equation atmospheric models.
J. Atmos. Sci
,
30
,
223
233
.
Goedheer
,
W. J.
, and
J. H. M.
Potters
,
1985
:
A compact finite difference scheme on a non-equidistance mesh.
J. Comput. Phys
,
61
,
269
279
.
Haney
,
R. L.
,
1991
:
On the pressure gradient force over steep topography in sigma coordinate ocean models.
J. Phys. Oceanogr
,
21
,
610
619
.
Hirsh
,
R. S.
,
1975
:
Higher order accurate difference solutions of fluid mechanics problems by a compact differencing technique.
J. Comput. Phys
,
19
,
90
109
.
Kreiss
,
H. O.
,
1975
:
Methods for the approximate solution of time dependent problem.
WMO GARP Rep. 13, 107 pp. [Available from Secretariat of the World Meteorological Organization, Case Postate No. 5, CH-1211, Geneva 20, Switzerland.]
.
——, and
Oliger
,
J.
,
1972
:
Comparison of accurate methods for the integration of hyperbolic equations.
Tellus
,
27
,
199
215
.
Lele
,
S. K.
,
1992
:
Compact finite difference schemes with spectral-like resolution.
J. Comput. Phys
,
103
,
16
42
.
McCalpin
,
J. D.
,
1994
:
A comparison of second-order and fourth-order pressure gradient algorithms in a σ-coordinate ocean model.
Int. J. Numer. Methods Fluids
,
18
,
361
383
.
Mellor
,
G. L.
,
T.
Ezer
, and
L-Y.
Oey
,
1994
:
The pressure gradient conundrum of sigma coordinate ocean models.
J. Atmos. Oceanic Technol
,
11
,
1126
1134
.
Navon
,
I. M.
, and
H. A.
Riphagen
,
1979
:
An implicit compact fourth-order algorithm for solving the shallow-water equations in conservation-law form.
Mon. Wea. Rev
,
107
,
1107
1127
.
Orzag
,
S. A.
,
1971
:
Numerical simulation of incompressible flows within simple boundaries.
J. Fluid Mech
,
49
,
75
112
.
Rubin
,
S. G.
, and
P. K.
Khosla
,
1977
:
Polynomial interpolation methods for viscous flow calculations.
J. Comput. Phys
,
24
,
217
244
.
Smagorinsky
,
J.
,
1963
:
General circulation experiments with the primitive equations. I. The basic experiment.
Mon. Wea. Rev
,
91
,
99
164
.
Song
,
Y. T.
,
1998
:
A general pressure gradient formulation for the ocean models. Part I: Scheme design and diagnostic analysis.
Mon. Wea. Rev
,
126
,
3213
3230
.
——, and
Wright
,
D. G.
,
1998
:
A general pressure gradient formulation for the ocean models. Part II: Momentum and bottom torque consistency.
Mon. Wea. Rev
,
126
,
3213
3230
.

APPENDIX

Error Estimation of AP Schemes on a Standard Grid

The truncation error is usually obtained from Taylor series expansion. The first and second derivatives at the grid point xi are represented by

 
formula

which shows the second-order accuracy of Δfx, Δ2fx2. Also, f′ and f″ at the grid point xi can be represented by (standard Pade scheme)

 
formula

which indicates the fourth-order accuracy; that is, the difference between the first-order derivative and its fourth-order difference counterpart at two grid points xi+1 and xi−1 should be

 
formula

Subtraction of (A2) from (A1) leads to

 
formula

The difference between the second-order derivative and its fourth-order difference at two grid points xi+1 and xi−1 should be

 
formula

Subtraction of (A4) from (A3) leads to

 
formula

Chu and Fan (1998) provide the following error estimation:

 
formula

Elimination of (fi+1fi−1) from (A5) and (A6) leads to

 
formula

Equations (A8) and (A7) indicate sixth-order accuracy for the AP schemes (7) and (8).

Footnotes

Corresponding author address: Dr. Peter Chu, Dept. of Oceanography, Naval Postgraduate School, Monterey, CA 93943. Email: chu@nps.navy.mil