Abstract

The Robert–Asselin time filter is widely used in numerical models of weather and climate. It successfully suppresses the spurious computational mode associated with the leapfrog time-stepping scheme. Unfortunately, it also weakly suppresses the physical mode and severely degrades the numerical accuracy. These two concomitant problems are shown to occur because the filter does not conserve the mean state, averaged over the three time slices on which it operates. The author proposes a simple modification to the Robert–Asselin filter, which does conserve the three-time-level mean state. When used in conjunction with the leapfrog scheme, the modification vastly reduces the impacts on the physical mode and increases the numerical accuracy for amplitude errors by two orders, yielding third-order accuracy. The modified filter could easily be incorporated into existing general circulation models of the atmosphere and ocean. In principle, it should deliver more faithful simulations at almost no additional computational expense. Alternatively, it may permit the use of longer time steps with no loss of accuracy, reducing the computational expense of a given simulation.

1. Introduction

From a functional perspective, the task of predicting future weather and climate may be reduced to the following iterative procedure. First, given the state of the atmosphere, ocean, and other Earth-system components at any time (the input), use the governing equations to compute the state at a slightly later time (the output). Then, repeat the loop as many times as required, always using the previous output as the next input.

The above prediction framework presents three main challenges, each of which potentially degrades the reliability of the forecast. First, Earth observations, which always contain measurement errors, are required to serve as the initial state. Second, the vast array of active physical processes and interactions is incompletely known and imperfectly represented in the spatially truncated governing equations. Third, the discrete stepping from one time level to the next is merely an approximation to the exact time-continuous evolution. This paper presents a possible avenue for progress with the third of these three challenges, which has received scant attention compared to the extensive research efforts devoted to the first two.

Pfeffer et al. (1992) have assessed the sensitivity of an atmospheric general circulation model to the time stepping, for fixed spatial discretization and physical parameterizations. They find that two different methods of time discretization—the leapfrog and Matsuno schemes—result in significant quantitative differences in the simulated climate. For example, the leapfrog scheme gives much more precipitation over the western tropical Pacific Ocean and less precipitation over the western North Atlantic Ocean. Therefore, time stepping appears to be an important contributor to model error.

Many different time-stepping methods have been proposed, including the leapfrog scheme (or centered difference scheme), the Matsuno scheme (e.g., Pfeffer et al. 1992), the Adams–Bashforth family of schemes (e.g., Durran 1991) and the Runge–Kutta family of schemes (e.g., Kar 2006). The leapfrog scheme has emerged as the method of choice in weather and climate models, despite related disciplines choosing differently (e.g., Runge–Kutta schemes are widely used in computational fluid dynamics but hardly ever used in numerical weather prediction) and despite evidence that other methods may be superior (e.g., the third-order Adams–Bashforth scheme is more accurate; Durran 1991). The leapfrog scheme is used so widely in weather and climate models probably because it is easy to implement, computationally inexpensive, and has low run-time storage requirements. Indeed, today’s widespread use of the leapfrog scheme in general circulation models is perhaps merely a legacy of computer memory having been such a severe constraint when the models were first developed.

A major problem with the leapfrog scheme is that it admits spurious computational modes (e.g., Mesinger and Arakawa 1976; Haltiner and Williams 1980; Durran 1999). In general, a differential equation that is first order in time has one degree of freedom, but an n-time-level numerical approximation to it constitutes an (n − 1)th-order difference equation with n − 1 degrees of freedom. Of these n − 1 modes, one is the physical mode and the remaining n − 2 are computational modes. The leapfrog is a three-time-level scheme, so one computational mode arises in it, in addition to the physical mode, because a second-order difference equation is used to approximate a first-order differential equation. The computational mode (or parasitic mode) is manifest as a spurious oscillation between even and odd time steps, which is referred to as time splitting.

One possible solution to time splitting is to periodically reinitialize the leapfrog scheme by applying a single step of a two-time-level scheme, which does not admit any computational modes. For example, Pfeffer et al. (1992) apply a single Matsuno step after every 11 leapfrog steps. This approach does not remove the computational mode, but merely resets its amplitude to zero periodically so that it never becomes large enough to be problematic.

The far more widely used solution to time splitting is to apply a time filter during the time-stepping procedure. Robert (1966) designed such a filter for the leapfrog scheme and Asselin (1972) showed that it selectively suppresses the computational mode but leaves the physical mode relatively undamped at low frequencies. The filter is now referred to as the Robert filter, the Asselin filter, or the Robert–Asselin filter. The behavior of the filter has been investigated not only for simple equation sets, with no space dependence, but also for the shallow-water equations (Schlesinger et al. 1983) and the hydrostatic primitive equations (Cordero and Staniforth 2004).

As testament to the filter’s success, Asselin (1972),1 has been cited over 450 times according to one citation database, mostly in journals of meteorology and the atmospheric sciences (around 300 citations) but also in journals of oceanography (around 100 citations) and fluid mechanics (around 50 citations). Examples include the use of the filter in models of regional climate (e.g., Caya and Laprise 1999), palaeoclimate (e.g., Fraedrich et al. 2005), ocean circulation (e.g., Griffies et al. 2001), geophysical fluid dynamics (e.g., Ford 1994; Bartello 2002), rotating laboratory fluids (e.g., Williams et al. 2009), and the atmosphere of Mars (e.g., Hartogh et al. 2005). André Robert’s contributions to numerical modeling, including the time filter, have been reviewed by Staniforth (1997) following a memorial symposium held at the University of Québec in 1994.

Currently, the Robert–Asselin filter is used in

  • operational numerical weather prediction models, including the Mesoscale Model (MSM) of the Japan Meteorological Agency (JMA), the global model of the Australian Bureau of Meteorology (BOM) and Research Centre (BMRC), the global model (GME) and regional model (COSMO-EU) of Deutscher Wetterdienst (DWD), and the Royal Netherlands Meteorological Institute (KNMI) model;

  • atmospheric general circulation models for climate simulation, including the ECHAM5 model of the Max-Planck-Institut für Meteorologie (MPI-M) and the Community Atmosphere Model (CAM) of the National Center for Atmospheric Research (NCAR);

  • ocean general circulation models, including Océan Parallélisé (OPA), the Nucleus for European Modeling of the Ocean (NEMO), the oceans of Met Office Hadley Centre climate models [i.e., the Hadley Centre Coupled Climate Model version 3 (HadCM3), the Hadley Centre Global Environmental Model (HadGEM), the High Resolution Global Environmental Model (HiGEM), and the Fast Met Office U.K. Universities Simulator (FAMOUS)], the Hybrid Coordinate Ocean Model (HYCOM), and (as an option) the Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model (MOM); and

  • models of the fluids in rotating annulus laboratory experiments, including the Quasi-Geostrophic Model for Investigating Rotating fluids Experiments (QUAGMIRE) and the Met Office/Oxford Rotating Annulus Laboratory Simulation (MORALS).

Despite its unquestioned success, the Robert–Asselin-filtered leapfrog scheme suffers from two related problems. First, in addition to suppressing the computational mode, the scheme also weakly suppresses the physical mode. Therefore, physical quantities (e.g., energy) that are conserved by the time-continuous equations are not necessarily conserved by the time-discretized equations when the filter is activated. The damping and nonconservation may be benign for sufficiently short integrations, but possibly not for longer ones.

Second, the Robert–Asselin filter severely degrades the leapfrog scheme’s numerical accuracy, measured as the rate at which the error tends to zero as the time step is progressively refined. Specifically, a numerical scheme is defined to be nth-order accurate if, after a given time interval, the difference between the numerical solution of the time-discretized equations with time step Δt, and the exact solution of the time-continuous equations, scales as (Δt)n as Δt → 0. Higher-order schemes are generally preferred to lower-order schemes, because they may permit the use of longer time steps with no loss of accuracy, reducing the computational expense of a given simulation. The Robert–Asselin-filtered leapfrog scheme is only first-order accurate for amplitude errors (although higher-order contributions may dominate for very small values of the filter parameter).

Because the Robert–Asselin filter is used so widely, simple-to-implement modifications that deliver more faithful simulations are very attractive. The author proposes such a modification in this paper. When used in conjunction with the leapfrog scheme, the modification vastly reduces the impacts on the physical mode and increases the numerical accuracy for amplitude errors by two orders, yielding third-order accuracy. Section 2 motivates the modified filter from a geometrical perspective. Section 3 derives analytically the amplification factor and numerical accuracy for the modified filtered leapfrog scheme, and compares them with the corresponding results for the standard Robert–Asselin-filtered leapfrog scheme. Section 4 concludes the paper with a summary and discussion.

2. The Robert–Asselin filter and proposed modification

The standard Robert–Asselin filter, and the modified filter proposed in this paper, are illustrated graphically in Fig. 1. Suppose that the values of a dependent variable, x, are given at three successive and equally spaced times, tn−1, tn, and tn+1. Then, from a geometrical perspective, the standard filter (Fig. 1a) operates by moving the inner point, with coordinates [tn, xn], a fraction ν toward the midpoint, [tn, (xn−1 + xn+1)/2], of the two outer points. Therefore, the displacement of the inner point under the influence of the filter is

 
formula

The filter parameter, ν, is usually chosen to be O(0.01–0.2).

Fig. 1.

Graphical comparison of the operation of (a) the standard Robert–Asselin filter and (b) the modified family of filters proposed in this paper. Points at three consecutive time levels are shown (marked with times signs) and a straight line is drawn between the two outer points (dashed). The standard filter moves the inner point through a displacement d, defined by (1). The modified filter moves the inner and right outer points through displacements αd and (α − 1)d, respectively, where 0 ≤ α ≤ 1. For the configuration of three points shown, d > 0.

Fig. 1.

Graphical comparison of the operation of (a) the standard Robert–Asselin filter and (b) the modified family of filters proposed in this paper. Points at three consecutive time levels are shown (marked with times signs) and a straight line is drawn between the two outer points (dashed). The standard filter moves the inner point through a displacement d, defined by (1). The modified filter moves the inner and right outer points through displacements αd and (α − 1)d, respectively, where 0 ≤ α ≤ 1. For the configuration of three points shown, d > 0.

Two relevant properties of the three points are their mean,

 
formula

and curvature,

 
formula

By displacing xn through the amount d, the standard Robert–Asselin filter reduces the magnitude of the curvature of the three points, |Cn|. When used in conjunction with the leapfrog scheme, this feature of the filter strongly suppresses the computational mode, as desired. But, crucially, the application of the filter does not conserve the three-time-level mean, Mn. The theoretical analysis of section 3 will show that, when used in conjunction with the leapfrog scheme, it is this feature of the filter that severely degrades the numerical accuracy.

In an attempt to include the possibility of conserving the three-time-level mean, the modified filter proposed in this paper (Fig. 1b) acts on the right outer point as well as the inner point. Specifically, for any α satisfying 0 ≤ α ≤ 1, the modified filter displaces xn through the amount αd and xn+1 through the amount (α − 1)d, where d is given by (1). All members of this family of modified filters reduce the magnitude of the curvature of the three points, |Cn|, with α controlling the relative contributions to the reduction from the displacements of the inner and right outer points. The special case α = 1 yields the standard Robert–Asselin filter discussed above, which displaces the inner point only. The special case α = 0 displaces the right outer point only. The special case α = 1/2 will be of particular interest in this paper, because it displaces the inner and right outer points equally and oppositely, conserving the three-time-level mean, Mn.

Before embarking upon the theoretical analysis, we briefly demonstrate the improvement that may be achieved by the proposed modification, when used in conjunction with the leapfrog scheme. We numerically integrate the equations of simple harmonic motion,

 
formula
 
formula

by alternately applying a leapfrog step and the modified filter with either (in two separate integrations) α = 1 (i.e., the standard Robert–Asselin filter) or α = 1/2. The numerical solutions so obtained are compared with each other, and with the exact solution, in Fig. 2, at the parameter values given in the caption. Amplitude errors are clearly much smaller with the modified filter than with the standard filter. One consequence is that X2 + Y2, which is conserved by the continuous equations and corresponds to the energy of the oscillation, decreases by 89% using the standard filter, but is approximately conserved using the modified filter, between the beginning and end of the integration shown in the figure.

Fig. 2.

Two numerical solutions to (4) and (5) with ω = 1 rad s−1, both obtained using the leapfrog scheme with Δt = 0.2 s. The computational mode is controlled using either the standard Robert–Asselin filter (with α = 1 and ν = 0.2) or the modified filter proposed in this paper (with α = 1/2 and ν = 0.2). A single two-time-level forward step is used to initiate the leapfrog scheme. The initial condition is X = 1, Y = 0, for which the exact solution (also plotted) is X = cosωt and Y = sinωt.

Fig. 2.

Two numerical solutions to (4) and (5) with ω = 1 rad s−1, both obtained using the leapfrog scheme with Δt = 0.2 s. The computational mode is controlled using either the standard Robert–Asselin filter (with α = 1 and ν = 0.2) or the modified filter proposed in this paper (with α = 1/2 and ν = 0.2). A single two-time-level forward step is used to initiate the leapfrog scheme. The initial condition is X = 1, Y = 0, for which the exact solution (also plotted) is X = cosωt and Y = sinωt.

3. Theoretical analysis

Amplitude and phase errors of time-stepping schemes are traditionally examined by analyzing solutions to the oscillation equation (e.g., Durran 1999), which, for the complex variable F(t), is

 
formula

where i = −1 and ω is a given (real) angular frequency. Equation (6) is related to (4) and (5) by F = X + iY. Using the modified filter proposed in section 2 to control the computational mode, the leapfrog scheme for (6), with time step Δt, is

 
formula
 
formula

and

 
formula

In this three-stage method, (7) implements the basic leapfrog scheme and (8) and (9) implement the modified filter, with 0 ≤ α ≤ 1. Here F denotes a provisional value, obtained by applying (7) during the current time step; F denotes another (singly filtered) provisional value, obtained by applying (9) during the current time step; and denotes the definitive (doubly filtered) value, obtained by applying (8) during the next time step. The occurrence of filtered values on the right sides of (7)(9) makes the scheme recursive: F is overwritten with F as soon as it is calculated, and so is F with .

It follows from (8) and (9) that the unfiltered, singly filtered, and doubly filtered values share a common complex amplification factor A, defined by

 
formula

Rewriting (7)(9) with function evaluations at time t only, using (10), yields three equations in the three unknowns: A, F(t)/F(t), and . Solving for A gives

 
formula

from which the numerical amplification factor is found to be

 
formula

Assuming 1 − ν/2 > 0, and taking the output of the square root operator to be the branch with nonnegative real part, then the positive sign (A+) corresponds to the physical mode and the negative sign (A) to the computational mode. For the special case α = 1, (12) reduces to the amplification factor derived by Asselin (1972) for the standard Robert–Asselin-filtered leapfrog scheme, as expected. The exact solution to (6) is F(t) = F(0) exp(iωt), from which the exact amplification factor is found to be

 
formula

for comparison with (12).

Figure 3 compares the amplification factor for the numerical solution, in the three cases α = 0, α = ½, and α = 1, with the amplification factor for the exact solution. The exact amplification factor (Fig. 3a) lies on the unit circle in the first quadrant and rotates anticlockwise as ωΔt increases from 0 to 1. The numerical amplification factors for the physical mode (Figs. 3b–d) also rotate anticlockwise in the first quadrant, but depart slightly from the unit circle as ωΔt increases from 0. The growing radii for α = 0 (Fig. 3b) and α = ½ (Fig. 3c), and the shrinking radius for α = 1 (Fig. 3d), correspond respectively to an artificial amplification and suppression of the physical mode. The numerical amplification factors for the computational mode (Figs. 3b–d) rotate clockwise in the second quadrant as ωΔt increases from 0 to 1. They each remain inside the unit circle, corresponding to a suppression of the computational mode, as desired.

Fig. 3.

Trajectories through the complex plane traced out by various amplification factors for the oscillation equation, (6), as ωΔt increases from 0 to 1. The plots compare (a) the exact amplification factor, given by (13), with (b)–(d) the numerical amplification factors for the modified filtered leapfrog scheme, given by (12). The filter parameters are (b) α = 0 and ν = 0.2, (c) α = 1/2 and ν = 0.2, and (d) α = 1 and ν = 0.2. The case α = 1, shown in (d), corresponds to the standard Robert–Asselin filter. In (b)–(d), solid lines denote the physical mode (A+) and dashed lines denote the computational mode (A). The unit circle, centered at the origin, is drawn in gray for reference.

Fig. 3.

Trajectories through the complex plane traced out by various amplification factors for the oscillation equation, (6), as ωΔt increases from 0 to 1. The plots compare (a) the exact amplification factor, given by (13), with (b)–(d) the numerical amplification factors for the modified filtered leapfrog scheme, given by (12). The filter parameters are (b) α = 0 and ν = 0.2, (c) α = 1/2 and ν = 0.2, and (d) α = 1 and ν = 0.2. The case α = 1, shown in (d), corresponds to the standard Robert–Asselin filter. In (b)–(d), solid lines denote the physical mode (A+) and dashed lines denote the computational mode (A). The unit circle, centered at the origin, is drawn in gray for reference.

As suggested by Fig. 3, the standard Robert–Asselin filter (α = 1) behaves qualitatively differently from all other filters in the modified family (α ≠ 1). Only the standard filter (Fig. 3d), exhibits a point in the complex plane at which the amplification factors for the physical and computational modes meet. The singularity occurs because α = 1 is the only case for which the imaginary term within the square root of (12) vanishes, allowing A+ = A at ωΔt = 1 − ν/2. For all other values of α, the presence of the imaginary term ensures that there is no value of ωΔt for which A+ = A, and the singular behavior of the standard Robert–Asselin filter is avoided.

Figure 4 shows in more detail how the magnitudes of the amplification factors depend upon ωΔt. The qualitatively different behavior between the cases α ≠ 1 and α = 1 is clearly visible. The singularity for the case α = 1 (Fig. 4c), which renders the curves nondifferentiable at ωΔt = 1 − ν/2, is replaced for the cases α = 0 (Fig. 4a) and α = ½ (Fig. 4b) with a smooth transition from the small ωΔt regime to the large ωΔt regime. The consequence for the physical mode (A+) is that, for the case α = 1 only, |A+| − 1 changes sign as ωΔt increases from 0 to 1, corresponding to a transition from artificial suppression to artificial amplification. The consequence for the computational mode (A) is that, for the case α = 1 compared to the cases α ≠ 1, the suppression is much less uniform as ωΔt increases from 0 to 1.

Fig. 4.

Magnitudes of various amplification factors for the oscillation equation, (6), plotted as functions of ωΔt. The plots compare the magnitude of the exact amplification factor, given by (13) to be unity, with the magnitudes of the numerical amplification factors for the modified filtered leapfrog scheme, given by (12). The filter parameters are (a) α = 0, (b) α = ½, and (c) α = 1, with various values of ν in each case. The case α = 1, shown in (c), corresponds to the standard Robert–Asselin filter. Solid lines denote the physical mode (A+) and dashed lines denote the computational mode (A).

Fig. 4.

Magnitudes of various amplification factors for the oscillation equation, (6), plotted as functions of ωΔt. The plots compare the magnitude of the exact amplification factor, given by (13) to be unity, with the magnitudes of the numerical amplification factors for the modified filtered leapfrog scheme, given by (12). The filter parameters are (a) α = 0, (b) α = ½, and (c) α = 1, with various values of ν in each case. The case α = 1, shown in (c), corresponds to the standard Robert–Asselin filter. Solid lines denote the physical mode (A+) and dashed lines denote the computational mode (A).

Figure 5 shows an enlarged view of how, as ωΔt → 0, the magnitudes of the numerical amplification factors for the physical mode approach the magnitude of the exact amplification factor. The limiting value, unity, appears to be approached much more rapidly for the case α = ½ than for the cases α = 0 and α = 1, suggesting a higher numerical accuracy. To confirm this suggestion, we Taylor expand the square root in (12) to obtain, for the physical mode,

 
formula

The coefficients, cn(α, ν), of the terms in the series expansion are real. The first few are given by

 
formula
 
formula
 
formula
 
formula
 
formula

It follows from (13) and (14) that

 
formula
Fig. 5.

Magnitudes of various amplification factors for the oscillation equation, (6), plotted as functions of ωΔt. The curves show how, as ωΔt → 0, the magnitudes of the numerical amplification factors for the physical mode of the modified filtered leapfrog scheme, given by (12), approach the magnitude of the exact amplification factor, given by (13) to be unity. The filter parameters are α = 0, α = ½, α = 0.53, and α = 1, with ν = 0.2 in each case. The case α = 1 corresponds to the standard Robert–Asselin filter.

Fig. 5.

Magnitudes of various amplification factors for the oscillation equation, (6), plotted as functions of ωΔt. The curves show how, as ωΔt → 0, the magnitudes of the numerical amplification factors for the physical mode of the modified filtered leapfrog scheme, given by (12), approach the magnitude of the exact amplification factor, given by (13) to be unity. The filter parameters are α = 0, α = ½, α = 0.53, and α = 1, with ν = 0.2 in each case. The case α = 1 corresponds to the standard Robert–Asselin filter.

For the special case α = 1, which corresponds to the standard Robert–Asselin filter, the coefficient of the quadratic term in (20) is negative, yielding |A+| − 1 ∼ −(ωΔt)2 as ωΔt → 0, in agreement with the leading-order behavior shown in Fig. 5. Hence, the amplitude error per time step varies as (Δt)2 and the amplitude error per unit time varies as Δt. Therefore, the numerical scheme with α = 1 is first-order accurate for amplitude errors and unconditionally stable. Alternatively, for the special case α = 0, the coefficient of the quadratic term in (20) is positive, yielding |A+| − 1 ∼ +(ωΔt)2 as ωΔt → 0, also in agreement with the leading-order behavior shown in Fig. 5. Therefore, the numerical scheme with α = 0 is first-order accurate for amplitude errors and unconditionally unstable.

In contrast, for the special case α = ½, the coefficient of the quadratic term in (20) vanishes and the coefficient of the quartic term is positive, yielding |A+| − 1 ∼ +(ωΔt)4 as ωΔt → 0, in agreement with the leading-order behavior shown in Fig. 5. Hence, the amplitude error per time step varies as (Δt)4 and the amplitude error per unit time varies as (Δt)3. Therefore, the numerical scheme with α = ½ is third-order accurate for amplitude errors and unconditionally unstable. In summary, the filter that conserves the three-time-level mean gives a numerical scheme that is two orders more accurate for amplitude errors than the standard Robert–Asselin filter. The increased accuracy may be exploited, by using it either to decrease the error for a given time step, or to increase the time step without increasing the error.

The unconditional instability of the case α = ½ may be avoided by instead choosing α ≳ ½, which almost conserves the three-time-level mean. The resulting filter is effectively a weighted blend of the third-order filter with α = ½ and the first-order filter with α = 1, the weighting of the latter being comparatively tiny. For this case, the coefficient of the quadratic term in (20) is negative, but very small, and the coefficient of the quartic term is positive. The negative quadratic term dominates for small ωΔt and the positive quartic term dominates for larger ωΔt, in agreement with the behavior shown in Fig. 5 for the case α = 0.53. Therefore, the numerical scheme with α ≳ ½ is conditionally stable. The finite stable range, for which |A+| ≤ 1, may be estimated from (20) by approximating the quartic term (but not the quadratic term) by its value when α = ½, to give

 
formula

or, for ν ≪ 1,

 
formula

These approximate formulas for the finite stable range work well. For example, for the case α = 0.53 and ν = 0.2, (21) gives 0 ≤ ωΔt ≤ 0.46, which is in good agreement with Fig. 5.

The numerical scheme with α ≳ ½ is strictly only first-order accurate for amplitude errors. However, it is clear from Fig. 5 that the error for the case α = 0.53 is much smaller than the error of the first-order scheme with α = 1. Indeed, the error for the case α − ½ ≪ 1 is comparable in magnitude to the error of the third-order scheme with α = ½ across much of the finite stable range. Therefore, for practical purposes, the numerical scheme with α ≳ ½ is as good as third-order accurate for amplitude errors.

Turning finally to consider phase errors, it follows from (13) and (14) that

 
formula

For all cases of α and ν in the practical range, the coefficient of the cubic term in (23) is positive, yielding arg(A+) − ωΔt ∼ +(ωΔt)3 as ωΔt → 0. Hence, the phase error per time step varies as (Δt)3 and the phase error per unit time varies as (Δt)2. Therefore, all numerical schemes in the modified family are second-order accurate for phase errors.

Table 1 summarizes the conservation, stability, and accuracy properties of the modified filter, when used in conjunction with the leapfrog scheme, for various values of α.

Table 1.

Summary of the conservation, stability, and accuracy properties of the modified filter proposed in this paper, when used in conjunction with the leapfrog scheme. The case α = 1 corresponds to the standard Robert–Asselin filter.

Summary of the conservation, stability, and accuracy properties of the modified filter proposed in this paper, when used in conjunction with the leapfrog scheme. The case α = 1 corresponds to the standard Robert–Asselin filter.
Summary of the conservation, stability, and accuracy properties of the modified filter proposed in this paper, when used in conjunction with the leapfrog scheme. The case α = 1 corresponds to the standard Robert–Asselin filter.

4. Summary and discussion

In the decades that have elapsed since the first general circulation models were developed, there have been major advances in the Earth observation systems from which initial conditions are derived, and in techniques for parameterizing unresolved physical processes. These advances have helped to improve the fidelity of weather and climate simulations. However, many general circulation models still use the same Robert–Asselin-filtered leapfrog time-stepping scheme as when they were first developed, despite evidence that time stepping appears to be an important contributor to model error.

This paper proposes a simple modification to the Robert–Asselin filter. The modified filter displaces the state at the future time slice as well as the current time slice. The modification yields a generalized family of filters, with a parameter controlling the relative sizes of the two displacements. The standard Robert–Asselin filter is a special case.

The behavior of the family of modified filters, when used in conjunction with the leapfrog scheme, is analyzed. The standard Robert–Asselin filter is shown to behave qualitatively differently from all other filters in the family. Each filter reduces the magnitude of the curvature at the three time slices operated upon, suppressing the computational mode. For the physical mode, each filter yields second-order accuracy for phase errors. But only the filter that conserves the three-time-level mean yields third-order accuracy for amplitude errors, whereas all other filters in the family (including the standard Robert–Asselin filter) yield only first-order accuracy. The filter that conserves the three-time-level mean yields an unconditionally unstable scheme, but conditional stability is recovered by adding a tiny amount of the standard Robert–Asselin filter, to yield a scheme that is as good as third-order accurate.

The modified filter proposed in this paper could improve weather and climate models. For example, it may permit the use of longer time steps with no loss of accuracy, reducing the computational expense of a given simulation. Alternatively, if the time step cannot be lengthened because it is constrained more strongly by other conditions (e.g., the CFL criterion) than by accuracy requirements, then the modified filter may permit an increase in accuracy at almost no additional computational expense. The modified filter would be extremely easy to implement in an existing computer model: the Robert–Asselin-filtered leapfrog routine could be upgraded by changing only a few lines of code. There may be a slight increase in the computational expense—the standard scheme is a two-stage method and the modified scheme is a three-stage method—but no extra function evaluations are required.

There are alternative methods for controlling the computational mode of the leapfrog scheme, which do not involve the application of a time filter. Kurihara (1965) proposed the leapfrog–trapezoidal method, which consists of obtaining a provisional value by applying a leapfrog predictor and then improving it by recursively applying a corrector. Being a predictor–corrector method, however, this scheme is iterative and potentially computationally demanding. Magazenkov (1980) proposed the alternate application, from one time step to the next, of a leapfrog step and a second-order Adams–Bashforth step. The need to execute a different algorithm at even- and odd-numbered time steps is cumbersome, however. In contrast, the modified filtered leapfrog scheme proposed in this paper is noniterative and nonalternating, yet still suppresses the computational mode and achieves third-order numerical accuracy for amplitude errors.

Déqué and Cariolle (1986) have shown that, despite the demonstrated ability of the standard Robert–Asselin filter to stabilize numerical solutions to atmospheric motion equations for certain combinations of temporal differencing and physical processes, in some other cases even a very weak filter may lead to an instability that can only be suppressed by a severe reduction of the time step. It remains to be seen whether the modified filter proposed in this paper also exhibits this unexpected behavior. Finally, it is possible that the modified filter could also improve the Robert–Asselin-filtered Adams–Bashforth schemes (e.g., Tandon 1987), but the exploration of this possibility is left for future work.

Acknowledgments

The author is funded through a Fellowship from the U.K. Natural Environment Research Council (NE/D009138/1). This study was supported in part by the National Science Foundation (PHY05-51164) and by a travel grant from the Royal Astronomical Society. The suggestions of two anonymous reviewers are gratefully acknowledged. The carbon footprint of this study—due to international air travel and office power consumption—is estimated by the author to be 3100 kg of CO2.

REFERENCES

REFERENCES
Asselin
,
R.
,
1972
:
Frequency filter for time integrations.
Mon. Wea. Rev.
,
100
,
487
490
.
Bartello
,
P.
,
2002
:
A comparison of time discretization schemes for two-timescale problems in geophysical fluid dynamics.
J. Comput. Phys.
,
179
,
268
285
.
Caya
,
D.
, and
R.
Laprise
,
1999
:
A semi-implicit semi-Lagrangian regional climate model: The Canadian RCM.
Mon. Wea. Rev.
,
127
,
341
362
.
Cordero
,
E.
, and
A.
Staniforth
,
2004
:
A problem with the Robert–Asselin time filter for three-time-level semi-implicit semi-Lagrangian discretizations.
Mon. Wea. Rev.
,
132
,
600
610
.
Déqué
,
M.
, and
D.
Cariolle
,
1986
:
Some destabilizing properties of the Asselin time filter.
Mon. Wea. Rev.
,
114
,
880
884
.
Durran
,
D. R.
,
1991
:
The third-order Adams–Bashforth method: An attractive alternative to leapfrog time differencing.
Mon. Wea. Rev.
,
119
,
702
720
.
Durran
,
D. R.
,
1999
:
Numerical Methods for Wave Equations in Geophysical Fluid Dynamics.
Springer-Verlag, 482 pp
.
Ford
,
R.
,
1994
:
Gravity-wave radiation from vortex trains in rotating shallow-water.
J. Fluid Mech.
,
281
,
81
118
.
Fraedrich
,
K.
,
H.
Jansen
,
E.
Kirk
,
U.
Luksch
, and
F.
Lunkeit
,
2005
:
The Planet Simulator: Towards a user friendly model.
Meteor. Z.
,
14
,
299
304
.
Griffies
,
S. M.
,
R. C.
Pacanowski
,
M.
Schmidt
, and
V.
Balaji
,
2001
:
Tracer conservation with an explicit free surface method for z-coordinate ocean models.
Mon. Wea. Rev.
,
129
,
1081
1098
.
Haltiner
,
G. J.
, and
R. T.
Williams
,
1980
:
Numerical Prediction and Dynamic Meteorology.
2nd ed. Wiley, 496 pp
.
Hartogh
,
P.
,
A.
Medvedev
,
T.
Kuroda
,
R.
Saito
,
G.
Villanueva
,
A.
Feofilov
,
A.
Kutepov
, and
U.
Berger
,
2005
:
Description and climatology of a new general circulation model of the Martian atmosphere.
J. Geophys. Res.
,
110
,
E11008
.
doi:10.1029/2005JE002498
.
Kar
,
S. K.
,
2006
:
A semi-implicit Runge–Kutta time-difference scheme for the two-dimensional shallow-water equations.
Mon. Wea. Rev.
,
134
,
2916
2926
.
Kurihara
,
Y.
,
1965
:
On the use of implicit and iterative methods for the time integration of the wave equation.
Mon. Wea. Rev.
,
93
,
33
46
.
Magazenkov
,
L. N.
,
1980
:
Time integration schemes for fluid dynamics equations, effectively damping the high frequency components (in Russian).
Tr. Gl. Geofiz. Obs.
,
410
,
120
129
.
Mesinger
,
F.
, and
A.
Arakawa
,
1976
:
Numerical Methods Used in Atmospheric Models.
Global Atmospheric Research Program (GARP) Publications Series 17, Vol. I, GARP, 64 pp
.
Pfeffer
,
R.
,
I.
Navon
, and
X.
Zou
,
1992
:
A comparison of the impact of two time-differencing schemes on the NASA GLAS climate model.
Mon. Wea. Rev.
,
120
,
1381
1393
.
Robert
,
A. J.
,
1966
:
The integration of a low order spectral form of the primitive meteorological equations.
J. Meteor. Soc. Japan
,
44
,
237
245
.
Schlesinger
,
R.
,
L.
Uccellini
, and
D.
Johnson
,
1983
:
The effects of the Asselin time filter on numerical solutions to the linearized shallow-water wave equations.
Mon. Wea. Rev.
,
111
,
455
467
.
Staniforth
,
A.
,
1997
:
André Robert (1929–1993): His pioneering contributions to numerical modeling.
Numerical Methods in Atmospheric and Oceanic Modelling: The André J. Robert Memorial Volume, C. A. Lin, R. Laprise, and H. Ritchie, Eds., NRC Research Press, 25–54
.
Tandon
,
M.
,
1987
:
Robert’s recursive frequency filter: A reexamination.
Meteor. Atmos. Phys.
,
37
,
48
59
.
Williams
,
P. D.
,
T. W. N.
Haine
,
P. L.
Read
,
S. R.
Lewis
, and
Y. H.
Yamazaki
,
2009
:
QUAGMIRE v1.3: A quasi-geostrophic model for investigating rotating fluids experiments.
Geosci. Model Dev.
,
2
,
13
32
.

Footnotes

Corresponding author address: Paul D. Williams, Department of Meteorology, University of Reading, P.O. Box 243, Earley Gate, Reading RG6 6BB, United Kingdom. Email: p.d.williams@reading.ac.uk

1

Robert (1966) predates the standard citation databases.

* Kavli Institute for Theoretical Physics Report 08-95.