Next Article in Journal
Initial-Value vs. Model-Induced Forecast Error: A New Perspective
Previous Article in Journal
Trends and Interdependence of Solar Radiation and Air Temperature—A Case Study from Germany
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Lagrange–Laplace Integration Scheme for Weather Prediction and Climate Modelling

School of Mathematics & Statistics, University College Dublin, D04 V1W8 Dublin, Ireland
Meteorology 2022, 1(4), 355-376; https://doi.org/10.3390/meteorology1040023
Submission received: 21 August 2022 / Revised: 16 September 2022 / Accepted: 21 September 2022 / Published: 27 September 2022

Abstract

:
A time integration scheme based on semi-Lagrangian advection and Laplace transform adjustment has been implemented in a baroclinic primitive equation model. The semi-Lagrangian scheme makes it possible to use large time steps. However, errors arising from the semi-implicit scheme increase with the time step size. In contrast, the errors using the Laplace transform adjustment remain relatively small for typical time steps used with semi-Lagrangian advection. Numerical experiments confirm the superior performance of the Laplace transform scheme relative to the semi-implicit reference model. The algorithmic complexity of the scheme is comparable to the reference model, making it computationally competitive, and indicating its potential for integrating weather and climate prediction models.

1. Introduction

The accuracy and efficiency of weather and climate models has been greatly enhanced by the introduction of better numerical algorithms for the solution of the equations of motion. Two of the most notable advances were the development of the semi-implicit (SI) scheme for treating the gravity-wave adjustment process and the semi-Lagrangian scheme for the advection processes. For a general review of semi-Lagrangian schemes, see [1,2]; for a review of recent and emerging time integration schemes, see [3].
Many modern operational NWP models use a semi-implicit scheme for time integration, increasing efficiency by enabling the use of a larger time step. However, this comes at a price: while stabilisation is achieved by slowing down the high-frequency gravity waves, the meteorologically significant components of the flow are also distorted by the time averaging of the SI scheme. For models that use a semi-Lagrangian advection scheme, the problem is magnified, as further increases in the time step result in greater errors.
It was pointed out in Lynch and Clancy [4] that the Laplace transform (LT) method with analytic inversion gives an exact treatment of the linear modes. This is due to the fact that the LT scheme does not involve time-averaging of the linear terms. Harney and Lynch [5] described a Laplace transform integration scheme in an Eulerian baroclinic model and showed that it yields more accurate forecasts than SI.
In this paper, we describe the implementation and performance of an integration scheme based on semi-Lagrangian advection and Laplace transform adjustment. The new scheme is incorporated in a spectral baroclinic atmospheric model Peak, described in Ehrendorfer’s book [6]. It is validated by comparison with a semi-Lagrangian semi-implicit scheme. Both schemes have been comprehensively verified against the Eulerian semi-implicit scheme originally implemented in the model and fully documented in [6].
Table 1 shows four integration schemes. EuSI is the original Eulerian semi-implicit scheme described in [6]. EuLT is the Eulerian Laplace transform scheme described in Harney and Lynch, 2019 [5]. The two Lagrangian schemes are denoted LaSI and LaLT.
Numerical tests were carried out to confirm the correctness of the model implementations. A comparison of their performance in simulating a growing baroclinically unstable disturbance is described in Section 4.1. The four simulations are quite similar, confirming the integrity of the codes.
The goal of this study is to demonstrate that Laplace transform adjustment outperforms the semi-implicit scheme for time steps typically used for semi-Lagrangian advection, and provides a means of producing more accurate weather forecasts and climate simulations. In Section 2, we review the semi-implicit semi-Lagrangian schemes and outline the Laplace transform adjustment scheme. We perform a simple analysis to compare the accuracy of the two adjustment schemes. In Section 3, we apply the Laplace transform scheme to the Peak model equations. Treatment of orography and of diffusion are given particular attention. In Section 4, a series of tests comparing simulations using the LaLT scheme and the semi-implicit scheme are described. Conclusions are presented in Section 5. Two appendices are included, one on the calculation of the commutator, and one relating the Laplace transform scheme to exponential integration schemes.

2. Outline of the Semi-Implicit and Laplace Transform Schemes

Before considering the details of the Peak model, we describe the general approach to integration using the LaSI and LaLT schemes. The original EuSI scheme is comprehensively documented by Ehrendorfer [6] and the EuLT scheme by Harney and Lynch [5].

2.1. Outline of the LaSI Scheme

We consider a general equation written in Lagrangian form
d X d t + L X + N ( X ) = 0
where L X are the linear terms and N ( X ) are the nonlinear terms. Before we convert this grid-point equation to the spectral domain, we discretise it in time using a three time-level Lagrangian scheme. For a trajectory from departure point D at time ( n 1 ) Δ t through arrival point A at ( n + 1 ) Δ t , the equations are
X A + X D 2 Δ t + L X A + + X D 2 + N M 0 = 0
where nonlinear terms are evaluated at the midpoint M and superscripts −, 0 and + denote values at times ( n 1 ) Δ t , n Δ t and ( n + 1 ) Δ t . We can write the solution at the advanced time as
X A + = [ I + Δ t L ] 1 ( I Δ t L ) X D 2 Δ t N M 0 ,
where the right hand side may be computed from known quantities. We now convert Equation (2) to spectral form, multiplying by Y n m ( λ , ϕ ) and integrating over the sphere to get
x + = ( I + Δ t L ) 1 ( I Δ t L ) x D 2 Δ t n M 0 ,
where x = x n m is the vector of spherical harmonic coefficients of X , and n = n n m is similarly related to N . The right hand side may be computed from known quantities. We note that the linear operator L commutes with the spectral transform.

2.2. Outline of the LaLT Scheme

We consider again the general equation in Lagrangian form (1). Before converting this grid-point equation to the spectral domain, we take the Laplace transform L along the trajectory starting from the departure point D at time ( n 1 ) Δ t :
s X ^ X D + L X ^ + N M 0 s = 0 ,
where s is the complex variable conjugate to t and X ^ : = L X . As for the LaSI scheme, we assume that the nonlinear term varies slowly, and approximate it by its value at the mid-point of the trajectory from ( D , ( n 1 ) Δ t ) to ( A , ( n + 1 ) Δ t ) .
To get an equation for X ^ we have to change the order of the Laplace transform and the linear operator L ; in general, these do not commute. We define the commutator
Γ : = [ L , L ] = L L L L
so that the term L X ^ = L { L X } becomes L X ^ + Γ ( X ) , and write the transformed equation
s X ^ X D + L X ^ + Γ M 0 + N M 0 s = 0 .
For computation of the commutator, see Appendix A.
We now convert to spectral form, integrating over the sphere, to get
s x ^ x D + L x ^ + γ M 0 + n M 0 s = 0 ,
where x ^ = x ^ n m is the vector of spherical harmonic coefficients of X ^ , x D are the coefficients of X D , and γ = γ n m and n = n n m are the coefficients of Γ and N . The solution may be written
x ^ = ( s I + L ) 1 x D γ M 0 n M 0 s
where, after inverse transformation to the time domain, the right hand side may be computed from known quantities.

2.3. Accuracy Analysis

An accuracy analysis of the Laplace transform (LT) and semi-implicit (SI) schemes was presented in [5]. The main conclusions are reviewed here. Considering a simple oscillation equation
X t = i ω X + N ( X )
that represents a component of the full system used in numerical weather prediction, the analysis showed that LT is more accurate than SI for both the linear and nonlinear terms. Assuming that the nonlinear term varies slowly, we take it to be a constant N ¯ . The exact solution of (3) at time ( n + 1 ) Δ t is then
X + = exp ( 2 i ω Δ t ) X + exp ( 2 i ω Δ t ) 1 i ω N ¯
where X is the solution at time ( n 1 ) Δ t .
The SI approximation to (3) is
X + X 2 Δ t = i ω X + + X 2 + N ¯
Solving for the new value X + , we have
X + = 1 + i ω Δ t 1 i ω Δ t X + 1 1 i ω Δ t 2 Δ t N ¯
Comparing this with the exact solution (4), we find that both the linear and nonlinear components of the solution are misrepresented. For the exact solution (4), X is multiplied by exp ( 2 i ω Δ t ) , while for the SI solution (5) the multiplier is
ρ = 1 + i ω Δ t 1 i ω Δ t .
Thus, although there is no amplification, the phase error increases with ω Δ t . For the nonlinear term, the SI scheme has both modulus and phase errors; for details, see Section 2 in [5]. The errors in the SI scheme are significant when the time step is large.
We now apply the Laplace Transform L to Equation (3), taking the origin of time to be ( n 1 ) Δ t . The transform of (3) is
s X ^ X = i ω X ^ + N ¯ s .
The solution for X ^ follows immediately:
X ^ = X s i ω + N ¯ s ( s i ω ) .
The inverse Laplace transform L 1 with time set to 2 Δ t yields X + equal to the exact solution (4) at time ( n + 1 ) Δ t . Thus, to the extent that the nonlinear term can be regarded as constant, the LT scheme is free from error.

2.4. Filtering with the LT Scheme

The LT scheme filters high frequency components by using a modified inversion operator L * : this can be done numerically by distorting the Bromwich contour for the inversion integral to a closed curve excluding poles associated with the high frequencies, as in [7,8,9,10,11]. In the present study, as in those of Lynch and Clancy [4] and Harney and Lynch [5], we invert the transform analytically, explicitly eliminating components with frequency greater than a specified cut-off frequency ω c .
Filtering may be done with a sharp cut-off at ω c , or with a smooth function such as a Butterworth filter having frequency response function
H ( ω ) = 1 1 + ( ω / ω c ) L .
Formally, we can define the modified inversion operator as the composition of the filter and the inverse Laplace transform: L * = L 1 H . In the numerical tests reported in Section 4 we set L = 16 and choose the cut-off period τ c = 1 h and ω c = 2 π / τ c . These values were chosen on the basis of many tests which showed that they yielded the best results for the given model resolution. Assuming that | ω | ω c , the inverse Laplace transform of (6) at time ( n + 1 ) Δ t gives
X + = exp ( 2 i ω Δ t ) X + exp ( 2 i ω Δ t ) 1 i ω N ¯ .
This agrees with the exact analytic result (4) for both the linear and nonlinear terms.
Equation (8) is formally identical to Equation (4) in Cox and Matthews [12], for which the local truncation error is 1 2 Δ t 2 ( d N / d t ) n . Assuming that the nonlinear term varies slowly, we have d N / d t = O ( Δ t ) , and the local truncation error is third order in time. We note that (8) is also equivalent to a component of the vector Equation (2) in Clancy and Pudykiewicz [13]. For the relationship between Laplace transform integration and exponential integrators, see Appendix B.

3. Laplace Transforming the Peak Equations

We now describe the application of the Laplace transform scheme to the Peak model equations. We begin with Equations (9.38)–(9.41) in Section 9.3 of Ehrendorfer [6], but written in Lagrangian form:
d ζ d t = F ζ ς ζ
d δ d t = F δ σ 2 ( R T ¯ π + Φ S + G T ) ς δ
d T d t = F T H δ ς T
p T d Π d t = p T δ
The notation is generally conventional and equations similar to these have appeared frequently, going back to Hoskins and Simmons [14]. The dependent variables are vorticity ( ζ ), divergence ( δ ), temperature ( T ) and log surface pressure π = log ( p S / p 00 ), where p 00 = 10 5 Pa is the reference pressure. All variables are in the spectral domain but, for compactness, we suppress the spectral indices so that, for example, ζ m is written simply as ζ . Bold-face variables are vectors with values at all M model levels, which are equally spaced in σ -coordinates. The explicit expressions for the matrices G and H are given in Ehrendorfer [6] and vectors are defined as follows: p T = ( Δ σ 1 , Δ σ 2 , Δ σ M ) , Π = ( π , π , , π ) and T ¯ = ( T ¯ , T ¯ , , T ¯ ) . Linear damping with coefficient ς may be applied to all variables except the surface pressure. The surface orography is represented by Φ S .
Diffusion is generally required to control spurious oscillations. Explicit diffusion schemes may become unstable for large time steps. To avoid this problem, we may use a scheme that is fully implicit. However, when combined with the Laplace transform scheme, this leads to additional complications. To circumvent these, we employ an operator splitting method. In the first stage, the inviscid Equations (9)–(12) with ς = 0 , are advanced to time ( n + 1 ) Δ t . This is then followed by a second stage in which the damping terms are integrated analytically (see Section 3.8).

3.1. Lagrangian Departure Points and Values

The Cartesian coordinates of a grid point at latitude ϕ and longitude λ are
X A = r cos λ cos ϕ , Y A = r sin λ cos ϕ , Z A = r sin ϕ .
We compute ( X A , Y A , Z A ) for each grid-point and store them; they are independent of the model level. We use a method described by McGregor [15] to determine the departure points ( X D , Y D , Z D ) . The inverse transformation is
r = X D 2 + Y D 2 + Z D 2 , ϕ = arcsin Z D / r , λ = arctan ( Y D / X D ) .
(We use the Fortran function atan2 ( Y D , X D ) here). The inverse transform must be done every time step, as the departure points change. They also depend on the model level.
Values of the variables at the departure points are obtained by interpolation in ( λ , ϕ ) space. A bicubic interpolation scheme is used. To facilitate interpolation, the grid point arrays are extended by two rows or columns around the boundaries. In the east–west direction, the extension is cyclic. At the poles, we repeat the two boundary rows, reversing their order and shifting the longitude by 180 to allow for crossing the poles [15].

3.2. Orography

Preliminary tests with real data produced instabilities that were insensitive to damping and to spectral truncation. The problem first became manifest over the Andes and Himalayas. As is well known, orography can lead to problems with semi-Lagrangian integration schemes. Ritchie and Tanguay [16] proposed a modification that alleviates this problem. An orographic term is subtracted from the pressure, yielding a much smoother field that is more accurately interpolated to the departure points.
The method of Ritchie and Tanguay is used in the IFS model and described in the documentation [17]. The variable π = log p S / p 00 in (12) is separated into a constant part involving the orography and a variable part independent of orography:
π = π ¯ + π
where π ¯ = Φ S / R T ¯ is, by definition, the value for an isothermal hydrostatic atmosphere. Thus, π ¯ = ( 1 / R T ¯ ) Φ S . We can then write
d π d t k = d π d t k + ( F π ) k , where ( F π ) k = 1 R T ¯ V k · Φ S .
The term F π involving the gradient of orography is computed in an Eulerian manner, and the continuity Equation (12) becomes
d π d t = p T ( F π δ ) .
The advected variable π is much smoother than the original variable π since the underlying orography has been extracted.

3.3. Inviscid Stage

We integrate the equations using a split scheme, where the diffusion is omitted in the first stage and integrated analytically in the second. The nonlinear terms are evaluated at the central time and averaged between the departure and arrival points:
F M 0 = 1 2 ( F D 0 + F A 0 ) .
This averaging is similar to the approach adopted in the ECMWF IFS model [17].
We now take the Laplace transform of (10)–(12), after applying (15) to change the pressure variable. Since the Laplace operator L does not commute with the spacial Laplacian σ 2 , we introduce the commutator
Γ : = s 2 [ L , σ 2 ] ( R T ¯ π + G T )
The factor s 2 is included here to ensure that Γ is independent of s (see Appendix A). The transformed equations may be written
s δ ^ = ( F δ ) M 0 / s Γ M 0 / s 2 + δ D σ 2 ( R T ¯ π ^ + G T ^ )
s T ^ = ( F T ) M 0 / s H δ ^ + T D
s π ^ = p T [ ( F π ) A 0 δ ^ ] + p T Π D
Using the equations for T ^ and π ^ , these quantities are eliminated from the divergence equation to obtain an equation for a single variable, δ ^ :
s 2 δ ^ = ( F δ ) M 0 σ 2 ( R T ¯ s π ^ + G s T ^ ) Γ M 0 / s + s δ D = ( F δ ) M 0 σ 2 [ R T ¯ p T δ ^ + p T Π D ) ] σ 2 [ G ( ( F T ) M 0 / s H δ ^ + T D ) ] Γ M 0 / s + s δ D
We transfer all terms involving δ ^ to the left:
[ s 2 I σ 2 ( R T ¯ p T + G H ) ] δ ^ = ( F δ ) M 0 Γ M 0 / s + s δ D σ 2 [ R T ¯ p T Π D + G ( F T ) M 0 / s + G T D ]

3.4. Transforming to Vertical Eigenmodes

We now define the vertical structure matrix
B : = R T ¯ p T + G H .
Since B is symmetric, its eigen-structure can be written
B E = E Λ o r E T B E = Λ o r B = E Λ E T
We now transform to vertical eigen-modes by multiplying (21) by E T :
s 2 I σ 2 ( E T B E ) ( E T δ ^ ) = E T ( F δ ) M 0 Γ M 0 / s + s δ D σ 2 [ R T ¯ p T Π D + G ( F T ) M 0 / s + G T D ] .
We note that, for a specific spectral component with total wavenumber , the Laplacian has a simple form σ 2 = ( + 1 ) / a 2 . Then, defining Ω 2 : = σ 2 Λ , the matrix on the left hand side of (22) is
[ s 2 I σ 2 ( E T B E ) ] = [ s 2 I σ 2 Λ ] = [ s 2 I + Ω 2 ]
Since this is a diagonal matrix, the equation falls into M separate scalar equations, one for each vertical mode. We group the right hand terms of (22) according to powers of s:
RHS = E T A × s + B × 1 + C × ( 1 / s )
where the vectors A , B and C are
A = δ D B = ( F δ ) M 0 σ 2 ( R T ¯ p T Π D + G T D ) C = σ 2 ( G ( F T ) M 0 ) Γ M 0 + σ 2 ( R T ¯ F π ) A 0 .
Multiplying (22) by the inverse of the diagonal matrix s 2 I + Ω 2 , the equation for the k-th component is
( E T δ ^ A + ) k = s s 2 + Ω k 2 ( E T A ) k + 1 s 2 + Ω k 2 ( E T B ) k + 1 s ( s 2 + Ω k 2 ) ( E T C ) k
We apply the operator L * to (24), noting that the vertical transform and Laplace transform commute. The terms can be inverted using standard results from Laplace transform theory [18]. The value at time ( n + 1 ) Δ t is denoted by a + superscript:
( E T δ A + ) k = H ( Ω k ) cos 2 Ω k Δ t ( E T A ) k + H ( Ω k ) sin 2 Ω k Δ t Ω k ( E T B ) k + 1 H ( Ω k ) cos 2 Ω k Δ t Ω k 2 ( E T C ) k
The filter response function H ( ω ) was defined in Equation (7) above.

3.5. Inverse Vertical Transformation

Let us define four diagonal matrices
Λ A = diag H ( Ω k ) cos 2 Ω k Δ t Λ B = diag H ( Ω k ) sin 2 Ω k Δ t Ω k Λ C = diag 1 H ( Ω k ) cos 2 Ω k Δ t Ω k 2 Λ D = diag 2 Ω k Δ t H ( Ω k ) sin 2 Ω k Δ t Ω k 3
( Λ D will be needed below). Then (25) can be written
E T δ A + = Λ A E T A + Λ B E T B + Λ C E T C
We can now calculate the divergence at the advanced time,
δ A + = E ( E T δ A + ) = E Λ A E T A + E Λ B E T B + E Λ C E T C
For compactness, we define the propagation matrices:
P A = E Λ A E T P B = E Λ B E T P C = E Λ C E T P D = E Λ D E T
( P D will be used below). Then we can write the solution as
δ A + = P A A + P B B + P C C
The P -matrices can be pre-computed and stored, since they do not depend on the model variables.

3.6. Temperature and Pressure

We return to Equations (19) and (20):
s T ^ = ( F T ) M 0 / s H δ ^ + T D s π ^ = p T [ ( F π ) A 0 δ ^ ] + p T Π D
Noting that L * { 1 / s } = 1 and L * { 1 / s 2 } = t , dividing these equations by s and applying the operator L * at time 2 Δ t , we have
T A + = T D + 2 Δ t ( F T ) M 0 H L * δ ^ / s
π A + = p T Π D + 2 Δ t p T ( F π ) A 0 L * δ ^ / s .
Both (28) and (29) require computation of
δ ˜ = L * δ ^ / s .
This term involves a convolution integral that may be approximated by the trapezoidal rule
L * δ ^ / s = 0 2 Δ t δ d t 2 Δ t δ D + δ A + 2 = 2 Δ t δ ¯
where δ ¯ is the average of old and new values of δ . However, this method was found not to perform in a satisfactory manner. We therefore employed an alternative strategy: the term (30) was computed by noting that
δ ˜ = L * δ ^ / s = E L * E T δ ^ / s .
We divide (24) by s to give
E T δ ^ s k = 1 s 2 + Ω k 2 ( E T A ) k + 1 s ( s 2 + Ω k 2 ) ( E T B ) k + 1 s 2 ( s 2 + Ω k 2 ) ( E T C ) k
We invert this using standard results for Laplace transforms to obtain
L * E T δ ^ s k = H ( Ω k ) sin 2 Ω k Δ t Ω k ( E T A ) k + 1 H ( Ω k ) cos 2 Ω k Δ t Ω k 2 ( E T B ) k + 2 Ω k Δ t H ( Ω k ) sin 2 Ω k Δ t Ω k 3 ( E T C ) k
Then using the Λ -matrices, we can write
L * E T δ ^ s = Λ B E T A + Λ C E T B + Λ D E T C
Noting (31) and using the P -matrices, we can now write
δ ˜ = P B A + P C B + P D C .
Finally, using (28) and (29), the values of T A + and π A + are
T A + = T D + 2 Δ t [ ( F T ) M 0 H δ ˜ / 2 Δ t ]
π A + = p T Π D + 2 Δ t p T ( F π ) A 0 δ ˜ / 2 Δ t .

3.7. Integrating the Vorticity

To complete the inviscid stage of the time step, we take the Laplace transform of the vorticity Equation (9)
s ζ ^ = ζ D + ( F ζ ) M 0 / s
Dividing by s and applying L * at 2 Δ t yields
ζ A + = ζ D + 2 Δ t ( F ζ ) M 0 ,
which is a standard centred Lagrangian step along the trajectory. As the poles are at s = 0 , the Laplace transform has no filtering effect here.
As is usual with the leapfrog model, a Robert-Asselin filter [19] is applied to the prognostic variables to prevent separation of the solutions at odd and even time steps. The coefficient is fixed at ϵ = 0.03 in all cases.

3.8. Diffusion Stage

The governing equation for a spectral component of any of the variables δ , ζ or T may be written in the form
d Ψ d t = F ς Ψ
We split the right hand side into two parts and integrate them separately. We first integrate the inviscid equation
d Ψ d t = F
over a time interval 2 Δ t with initial condition Ψ and denote the result as Ψ * . We then integrate the equation
Ψ t = ς Ψ
analytically with the initial condition Ψ * to get
Ψ + = exp ( 2 Δ t ς ) Ψ *
which is the required solution at time ( n + 1 ) Δ t . This second stage is applied to the divergence, vorticity and temperature; the surface pressure is not damped in this way.
The parameter ς depends only upon the horizontal scale. The diffusion is assumed to be of the form
Q t = ν 2 2 Q + ν 6 6 Q
In the spectral domain, the damping coefficient becomes
ς = ν 2 ( + 1 ) a 2 + ν 6 ( + 1 ) a 2 3
We note that the spectral equations are unchanged in form by the addition of the hyper-diffusion term ( ν 6 ); only the value of the coefficient is changed. The ν 6 -term more strongly damps the smaller scales. Having applied diffusion, we have all the model variables at the advanced time, and a new time step can be taken.

4. Numerical Evaluation of the Integration Schemes

In this section, we describe a series of tests comparing simulations using the Laplace transform scheme (LaLT) and the semi-implicit scheme (LaSI). The LaSI and LaLT models were run with 20, 40 and 60 min time steps. For LaLT, a cut-off value τ c = 1 h was set for all time steps. In most cases, the reference forecast was an integration of the LaSI model with a time step Δ t = 10 min.
Eulerian models are subject to a Courant-Friedrichs-Lewy stability condition. For an advection speed u ¯ = 100 m/s, the non-dimensional stability ratio u ¯ Δ t / δ x is unity for a time step Δ t = 1500 s or 25 min. In fact, both the Eulerian models, EuSI and EuLT, were found to be unstable for a time step of 24 min. The Lagrangian models are not subject to this limitation.
The horizontal resolution of the model was at triangular truncation T 85 . The colocation grid corresponding to this has 256 × 129 grid points, with a grid interval of approximately 150 km. In all cases, there were 20 vertical levels, uniformly spaced in σ -coordinates.
The default setting of the diffusion coefficient was ν 2 = 7 × 10 5 m 2 s 1 : the damping of a component of total wavenumber is ς = ν 2 ( + 1 ) / a 2 s 1 . The default value implies an e-folding time of 2.2 h for the shortest waves represented at truncation T85. Sixth order diffusion ( ν 6 = 10 7 m 2 s 1 ) was also applied for runs with a 60 min time step. We note that LaLT consistently required less explicit diffusion than LaSI.

4.1. Initial Validation Tests

Numerous tests were carried out to confirm the correct operation of the model codes. For short time steps, the Eulerian and Lagrangian advection produced similar results, as did the semi-implicit and Laplace Transform adjustment. As an example of the performance of the four models, simulations of a growing baroclinically unstable disturbance [20] are shown in Figure 1. The initial field is a zonally symmetric flow with a small perturbation on the Greenwich meridian. The four models were integrated for 12 days, each with a time step Δ t = 5 min. The damping coefficient was ν 2 = 7 × 10 5 and all other parameter settings were equal for the four runs. The four simulations are quite similar, although we notice a slight damping for the Lagrangian runs, associated with the interpolations involved in the treatment of advection.

4.2. Kelvin Waves

Kelvin waves are eastward propagating waves that play an important role in atmospheric dynamics. Clancy and Lynch [7] showed that the LT scheme had a significantly smaller phase error than the semi-implicit scheme for the integration of these waves. Exact Kelvin Wave initial conditions can be generated using the method of Kasahara [21]. In this study, we use a simple analytical approximation described in [22]. We examine the solutions for zonal wave numbers 1 and 4. The wave amplitude is 100 m in both cases. The theoretical period for the Kelvin wave with zonal wavenumber m = 1 is about 32 h and for m = 4 is about 8.3 h ([21], Figure 9). Figure 2 shows that the root mean square errors for wavenumber 1 for LaLT (blue lines) are significantly smaller than for LaSI (red lines). Forecasts were run with time steps of 20, 40 and 60 min (solid, dashed and dotted lines, respectively).
For the LaSI forecast of wave number 4 (Figure 3) with a 20 min step, the propagation of the wave lags by a full wavelength by the end of the integration, reducing the rms error. The errors for the LaSI scheme with the larger time steps oscillate wildly as the wave moves into and out of phase with the reference solution. The errors for LaLT are much smaller and behave in a more realistic and steady manner.

4.3. The Five-Day Wave

The Five-day Wave RO(1,2) is the gravest symmetric rotational Hough mode of zonal wavenumber 1. It is closely related to the initial state chosen by Lewis Fry Richardson for their preliminary shallow-water forecast experiment ([23], Section 4.1). The initial conditions for a three-dimensional Five-day Wave were implemented in Peak. The initial pressure amplitude was 10 hPa, with mean pressure 1000 hPa. As no zonal mean flow was included, the wave has a period close to 5 days. Figure 4 shows the root mean square errors in surface pressure for the LaSI scheme (red) and the LaLT scheme (blue), with time steps of 20 min (solid lines), 40 min (dashed lines) and 60 min (dotted lines). The reference is an SI forecast with time step of 10 min. The error level for LaLT is significantly less than that of the LaSI scheme, especially for the longer time step. The scores for vorticity at 250 hPa (Figure 5) confirm the superior performance of LaLT.

4.4. Rossby-Haurwitz Wave

Rossby-Haurwitz (RH) waves are exact solutions of the nonlinear barotropic vorticity equation. While they are not eigenfunctions of the shallow water equations, they have frequently been used as test cases. Following Phillips [24], the RH(4,5) wave was chosen as Test Case 6 by Williamson et al. [25]. This test case has been extended to three dimensions; the initial vorticity field is as in the barotropic case, the divergence is zero and a vertical temperature profile and surface pressure field are defined; for details, see [26].
The LaSI and LaLT schemes, with time steps of 20, 40 and 60 min, were compared to a reference run of LaSI using a time step of 10 s. No diffusion was used for the reference or LaLT runs, but the LaSI runs were unstable. This was overcome by applying horizontal diffusion with a coefficient ν 2 = 3 × 10 6 m 2 s 1 . Figure 6 shows the root mean square errors in surface pressure for the LaSI scheme (red) and the LaLT scheme (blue). The error level for LaLT is substantially less than for the LaSI scheme. Figure 7 shows scores for vorticity near the tropopause (250 hPa). The forecasts for LaSI and LALT have comparable errors, but there is a slight advantage for the Laplace transform scheme.

4.5. Flow over a Mountain

Test Case 5 of Williamson et al. [25] treats a zonal flow over an isolated mountain. The mountain is centred at ( 90   E , 30   N ) with maximum height 2000 m. No analytic solution is known so, as usual, we take the LaSI run with Δ t = 10 min as a reference. Figure 8 shows the root mean square errors in surface pressure for the LaSI scheme (red) and LaLT scheme (blue), with time steps of 20 min (solid), 40 min (dashed) and 60 min (dotted lines). The error levels for the two schemes are very close in value.

4.6. A Baroclinically Unstable Wave

Polvani et al. [20] devised a test case for baroclinic instability. The initial conditions consist of a non-divergent zonal flow with constant surface pressure. A small perturbation is added to the temperature to trigger the development of baroclinic instability. The test case of Polvani et al. was used by Ehrendorfer [6] to validate the Peak model. With a fixed value for the diffusion coefficient, the initial conditions are ‘numerically convergent’ as shown in [20] using two different numerical models. A test case quite similar to that of Polvani et al. was constructed by Jablonowski and Williamson [27].
We use the test case of Polvani et al. [20] to show that the LaLT scheme can accurately simulate baroclinic development. Using this case, Harney and Lynch [5] showed that the EuLT scheme can accurately simulate baroclinic development. In Figure 1, above, we showed the 12-day forecasts for all four schemes all with a small time step Δ t = 5 min. There was no substantive difference between the four schemes. They are also indistinguishable from the results plotted in ([20], Figure 4). Thus, all four schemes are capable of forecasting baroclinic development with high precision.
Quantitative scores confirm that the differences in performance are small: Figure 9 shows the root mean square errors in surface pressure (hPa) for forecasts with the LaSI scheme (red) and LaLT scheme (blue), with time steps 20, 40 and 60 min. As usual, the reference forecast is LaSI with a 10 min time step. It is clear that the time truncation error grows with forecast range. It is also clear that the error is greater for larger time steps. The important point is that the errors for the two integration schemes are very similar in magnitude.

4.7. Real Data Test

The simple wave tests described above indicate superior accuracy for the LT scheme compared to the semi-implicit scheme. The ultimate conclusion on superiority of the scheme must involve comprehensive comparisons for a large range of meteorological conditions. As a first step, a single test using real atmospheric data is described here.
Data was retrieved from the European Centre for Medium-Range Weather Forecasts MARS archive. The date chosen was 00 UTC on 15 October 2017, the day before a major storm, Ophelia, reached Ireland. This data comprised temperature, divergence and vorticity fields on 25 pressure levels, surface pressure and the relevant orography field. These fields where interpolated onto the 20 sigma levels and reduced to the spectral resolution T85 used for the Peak forecasts. The process of interpolation introduced noise, which was removed by initialization, as described by Harney and Lynch [5]. Using initialized data, six forecasts were performed using the LaSI and LaLT schemes with time steps of 10, 20 and 40 min.
Figure 10 shows the root mean square error for surface pressure. The reference is a forecast using LaSI with a time step of 5 min. The red curves are for LaSI and the blue ones for LaLT. For the 10 min step, the errors are comparable for the two models, although the error during the initial day is smaller for LaLT. For the larger time steps, the Laplace transform scheme is clearly superior to the semi-implicit scheme. Scores for mid-troposphere vorticity (Figure 11) confirm the superior performance of LaLT.

5. Conclusions

An integration scheme using a Laplace transform for adjustment combined with a semi-Lagrangian advection scheme (LaLT) has been found to yield results at least as accurate as the popular semi-implicit semi-Lagrangian scheme. The numerical tests described in Section 4 clearly indicate the superior accuracy for the LaLT scheme compared to the semi-implicit scheme LaSI. The single experiment with real data reinforces this advantage.
Ultimate conclusions on the superiority of the LaLT scheme require more comprehensive comparisons for a large range of meteorological conditions. The potential for operational implementation of LaLT would depend upon more exhaustive testing with higher spatial resolution and incorporating a full package of physical processes.
There are well-known advantages of using a two time level scheme for Lagrangian advection. There appear to be no difficulties in principle combining such a scheme with Laplace transform adjustment.
The efficient formulation of the LaLT scheme, with analytical inversion of the Laplace transform, is made possible through the use of a spectral model. An active debate on the future of spectral models has been ongoing for decades. The global simulation of an entire season of the Earth’s atmosphere, with a 1 km grid and upper boundary of 80 km [28], suggests that spectral models will continue to be competitive in the future.
For practical reasons, the semi-Lagrangian method used in this study was applied only to the horizontal advection. However, there is no difficulty to include vertical advection in the scheme. The algorithmic complexity of the LaLT scheme is comparable to that of LaSI. Averaging over a range of integrations, the processing time for LaLT was just 6% more than for LaSI, so the Laplace transform scheme is computationally competitive.
In summary, the main result of this study is that the LaLT scheme is clearly superior in accuracy to the LaSI scheme for the large time steps typically used with semi-Lagrangian advection. The evidence presented gives a clear indication of the practical potential of the LaLT scheme for integrating weather and climate prediction models.

Funding

This research received no external funding.

Data Availability Statement

The FORTRAN code of the Peak model is available on the following website in precisely the same form as it is reprinted in Martin Ehrendorfer’s book [6]: https://sites.google.com/site/spectralnwpmodels/home/code (accessed on 23 September 2022).

Acknowledgments

My sincere thanks go to Martin Ehrendorfer for his comprehensive documentation of the Peak model in his book [6] and for making the model code freely available. I also thank Colm Clancy and Eoghan Harney of Met Éireann for several valuable scientific discussions on Zoom during the Covid lockdown. Finally, I thank the reviewers for their insightful comments.

Conflicts of Interest

The author declares no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
EuSIEulerian advection and semi-implicit adjustment scheme
EuLTEulerian advection and Laplace transform adjustment scheme
LaSILagrangian advection and semi-implicit adjustment scheme
LaLTLagrangian advection and Laplace transform adjustment scheme

Appendix A. Calculating the Commutator

It was assumed in earlier work on implementing the Laplace Transform Integration Scheme in a Semi-Lagrangian context that the Laplace Transform operator L along a trajectory commutes with spatial differential operators such as the gradient operator and Laplacian 2 . This is not the case, as may be shown by simple counter-examples. In this appendix we derive expressions for the commutator of the Laplace transform with the Laplacian operator 2 .

Two Laplace Transforms

For a function f ( x , t ) of space and time, we define two distinct Laplace transforms. The Euler–Laplace transform (ELT) is evaluated at a fixed point in space
L [ f ] ( x 0 , s ) = 0 C 00 e s t f ( x 0 , t ) d t .
Here, C 00 is a line in the x -t space parallel to the time axis and passing through the point ( x 0 , 0 ) .
The Lagrange–Laplace transform (LLT) is evaluated along a trajectory of the motion:
L [ f ] ( x 0 , s ) = 0 C 0 e s t f ( x ( t ) , t ) d t .
Here, C 0 is the trajectory of the motion starting at the point ( x 0 , 0 ) .
The Euler–Laplace transform commutes with spatial operators like x and 2 . This is not the case for the Lagrange–Laplace transform. We write
L = L + [ L , ]
and, to replace L by L , we require an expression for the commutator [ L , ] = ( L L ) .

Evaluating the Commutator [,2]

In Figure A1, we show the trajectory C 0 , starting at point x 0 0 , along which the transform L { f ( x 0 ) } is evaluated (superscript 0 indicates the initial time t = 0 ). Shown also are the trajectories C P starting from x P 0 and C M starting from x M 0 . The contours C + and C are not trajectories, but are parallel to C 0 , shifted to x + = x 0 + Δ x 0 and x = x 0 Δ x 0 .
Figure A1. Trajectories (solid) and pseudo-trajectories (dashed) in ( x , y , t ) -space. See text for details.
Figure A1. Trajectories (solid) and pseudo-trajectories (dashed) in ( x , y , t ) -space. See text for details.
Meteorology 01 00023 g0a1
The transform of x [ f ( x 0 ) ] is given by
L { x [ f ( x 0 ) ] } = L { f ( x + ) } L { f ( x ) } x P 0 x M 0 + H . O . T .
( H . O . T . denotes higher order terms). The derivative of the transform L { f ( x 0 ) } is given by
x [ L { f ( x 0 ) } ] = L { f ( x P ) } L { f ( x M ) } x P 0 x M 0 + H . O . T .
We expand the variables as follows
f + = f ( x + ) f + 0 + ( u 0 f x + v 0 f y ) 0 t f = f ( x ) f 0 + ( u 0 f x + v 0 f y ) 0 t f P = f ( x P ) f + 0 + ( u P f x + v P f y ) 0 t f M = f ( x M ) f 0 + ( u M f x + v M f y ) 0 t
Using these expansions in (A3) and (A4) we get
[ L , x ] f ( x 0 0 , s ) = 1 s 2 u x f x + v x f y 0 0 .
It is easy now to obtain the following commutators:
L , f = 1 s 2 [ u x f x + v x f y , u y f x + v y f y ] 0 0
L , · f = 1 s 2 [ u x f x x + v x f x y + u y f x y + v y f y y ] 0 0
We can also establish the identity
[ L , 2 ] = [ L , · ] f + · [ L , ] f
which leads to the result
[ L , 2 ] f = 1 s 2 2 u x f x x + v x f x y + u y f x y + v y f y y + u x x f x + v x x f y + u y y f x + v y y f y 0 0 .

Appendix B. LTI and Exponential Integrators

The LT method with analytic inversion gives an exact treatment of the linear modes. This is due to the fact that the LT scheme does not involve time-averaging of the linear terms. An alternative way of achieving accuracy is to use an exponential integrator (see, for example, [29,30]. In this appendix we demonstrate the relationship between Laplace transform integration and exponential integrators.
We may write the model equations in the form
X t = L X + N ( X ) .
where the matrix L has an orthogonal eigenvector matrix E with L E = E Λ . Assuming that the solution of (A8) at time t n = n Δ t is known, the Laplace transform with this initial time is
s X ^ X n = L X ^ + N ^
where L { X } = X ^ is the Laplace transform of the state vector. Solving for this, we get
X ^ = ( s I L ) 1 [ X n + N ^ ] .
We note that
( s I L ) 1 = E ( s I Λ ) 1 E T
and also note the transforms
( s I Λ ) 1 = L { exp ( Λ t ) } and ( s I L ) 1 = L { exp ( L t ) } .
We can write the nonlinear term as
( s I L ) 1 N ^ = L { exp ( L t ) } · L { N } .
The convolution theorem allows this to be written
( s I L ) 1 N ^ = L t n t exp ( L ( t τ ) N ( τ ) d τ .
The transformed Equation (A9) now becomes
X ^ = L { exp ( L t ) } X n + L t n t exp ( L ( t τ ) N ( τ ) d τ
We invert this at time t n + 1 = t n + Δ t to get
X n + 1 = e L t n + 1 X n + e L t n + 1 t n t n + 1 e L τ N ( τ ) d τ
We note that (A10) is formally identical to Equation (8) of Peixoto and Schreiber [30], which they call the variation-of-constants formula. We have thus established a close relationship between the Laplace transform scheme and exponential integrators.

Approximating the Nonlinear Term

The convolution term must be evaluated by approximate means, since it involves unknown quantities. Suppose we evaluate the nonlinear term at time t n and assume that it is constant throughout the time step ( t n , t n + 1 ) . Then the convolution integral can be evaluated, giving
X n + 1 = e L t n + 1 X n + e L t n + 1 t n t n + 1 e L τ d τ N n = e L t n + 1 X n + ( L ) 1 [ I exp ( L Δ t ) ] N n .
Assuming a small time-step, this reduces to
X n + 1 = e L t n + 1 X n + Δ t N n .
This is perhaps the simplest version of an exponential integrator. There is a wide range of more sophisticated and accurate approximations of the convolution integral. For example, we might estimate N at the centre of the time step by extrapolation N n + 1 / 2 = ( 3 N n N n 1 ) / 2 . Many other possibilities exist.
The time-averaging of the SI scheme also results in an error in the nonlinear term, even when this term is constant (see Harney and Lynch ([5], Equation (3)). In this ideal case, the LT scheme has no error in the nonlinear term ([5], Equation (4)).

References

  1. Coiffier, J. Fundamentals of Numerical Weather Prediction; Cambridge University Press: Cambridge, UK, 2011; 340p, ISBN 978-1-1070-0103-9. [Google Scholar]
  2. Durran, D.R. Numerical Methods for Wave Equations in Geophysical Fluid Dynamics; Springer: Berlin, Germany, 1999; 465p, ISBN 0-387-98376-7. [Google Scholar]
  3. Mengaldo, G.; Wyszogrodzki, A.; Diamantakis, M.; Lock, S.-J.; Giraldo, F.X.; Wedi, N.P. Current and Emerging Time-Integration Strategies in Global Numerical Weather and Climate Prediction. Arch. Comp. Methods Eng. 2019, 26, 663–684. [Google Scholar] [CrossRef]
  4. Lynch, P.; Clancy, C. Improving the Laplace transform integration method. Q. J. R. Met. Soc. 2016, 142, 1196–1200. [Google Scholar] [CrossRef]
  5. Harney, E.; Lynch, P. Laplace Transform Integration of a Baroclinic Model. Q. J. R. Met. Soc. 2019, 145, 347–355. [Google Scholar] [CrossRef]
  6. Ehrendorfer, M. Spectral Numerical Weather Prediction Models; Society for Industrial & Applied Mathematics (SIAM); Cambridge University Press: Cambridge, UK, 2012; pp. xxv+482; ISBN 978-1-61197-198-9. [Google Scholar]
  7. Clancy, C.; Lynch, P. Laplace transform integration of the shallow water equations. Part 1: Eulerian formulation and Kelvin waves. Q. J. R. Met. Soc. 2011, 137, 792–799. [Google Scholar] [CrossRef]
  8. Clancy, C.; Lynch, P. Laplace transform integration of the shallow water equations. Part 2: Lagrangian formulation and orographic resonance. Q. J. R. Met. Soc. 2011, 137, 800–809. [Google Scholar] [CrossRef]
  9. Lynch, P. Initialization using Laplace transforms. Q. J. R. Met. Soc. 1985, 111, 243–258. [Google Scholar] [CrossRef]
  10. Lynch, P. Initialization of a barotropic limited area model using the Laplace Transform technique. Mon. Weather Rev. 1985, 113, 1338–1344. [Google Scholar] [CrossRef]
  11. Van Isacker, J.; Struylaert, W. Laplace Transform Applied to a Baroclinic Model. In Short- and Medium-Range Numerical Weather Prediction, Proceedings of IUGG NWP Symposium, Tokyo, Japan, 4–8 August 1986; Matsuno, T., Ed.; Meteorological Society of Japan: Tokyo, Japan, 1986; pp. 247–253. [Google Scholar]
  12. Cox, S.M.; Matthews, P.C. Exponential time differencing for stiff systems. J. Comp. Phys. 2002, 176, 430–455. [Google Scholar] [CrossRef]
  13. Clancy, C.; Pudykiewicz, J.A. On the use of exponential time integration methods in atmospheric models. Tellus 2013, 65, 1–16. [Google Scholar] [CrossRef]
  14. Hoskins, B.J.; Simmons, A.J. A multi-layer spectral model and the semi-implicit method. Q. J. R. Met. Soc. 1975, 101, 637–655. [Google Scholar] [CrossRef]
  15. McGregor, J.L. Economical determination of departure points for semi-Lagrangian models. Mon. Weather Rev. 1993, 121, 221–230. [Google Scholar] [CrossRef]
  16. Ritchie, H.; Tanguay, M. A comparison of spatially averaged Eulerian and semi-Lagrangian treatment of mountains. Mon. Weather Rev. 1996, 124, 167–181. [Google Scholar] [CrossRef]
  17. ECMWF. IFS Documentation–Cy47r3. Part III: Dynamics and Numerical Procedures. Operational implementation 12 October 2021. 2021. Available online: https://www.ecmwf.int/en/elibrary/20197-ifs-documentation-cy47r3-part-iii-dynamics-and-numerical-procedures (accessed on 21 August 2022).
  18. Doetsch, G. Guide to the Applications of the Laplace and Z Transforms; Van Nostrand Reinhold: New York, NY, USA, 1971; 240p. [Google Scholar]
  19. Asselin, R. Frequency filter for time integrations. Mon. Weather Rev. 1972, 100, 487–490. [Google Scholar] [CrossRef]
  20. Polvani, L.M.; Scott, R.K.; Thomas, S.J. Numerically Converged Solutions of the Global Primitive Equations for Testing the Dynamical Core of Atmospheric GCMs. Mon. Weather Rev. 2004, 132, 2539–2552. [Google Scholar] [CrossRef]
  21. Kasahara, A. Normal Modes of Ultra-long Waves in the Atmosphere. Mon. Weather Rev. 1976, 104, 669–690. [Google Scholar] [CrossRef]
  22. Harney, E. Laplace Transform Integration Scheme in Barotropic and Baroclinic Atmospheric Models. Ph.D. Thesis, University College Dublin, Dublin, Ireland, 2018. [Google Scholar]
  23. Lynch, P. The Emergence of Numerical Weather Prediction: Richardson’s Dream; Cambridge University Press: Cambridge, UK, 2006. [Google Scholar]
  24. Phillips, N.A. Numerical integration of the primitive equations on the hemisphere. Mon. Weather Rev. 1959, 87, 333–345. [Google Scholar] [CrossRef]
  25. Williamson, D.L.; Drake, J.B.; Hack, J.J.; Jakob, R.; Swarztrauber, P.N. A standard test set for numerical approximation to the shallow water equations in spherical geometry. J. Comp. Phys. 1992, 102, 211–224. [Google Scholar] [CrossRef]
  26. Jablonowski, C.; Lauritzen, P.; Nair, R.; Taylor, M. 2008: Idealized Test Cases for the Dynamical Cores of Atmospheric General Circulation Models: A Proposal for the NCAR ASP 2008 Summer Colloquium. Available online: http://www-personal.umich.edu/~cjablono/NCAR_ASP_2008_idealized_testcases_29May08.pdf (accessed on 21 August 2022).
  27. Jablonowski, C.; Williamson, D. A baroclinic instability test case for atmospheric model dynamical cores. Q. J. R. Met. Soc. 2006, 132, 2943–2975. [Google Scholar] [CrossRef]
  28. Wedi, N.P.; Polichtchouk, I.; Dueben, P.; Anantharaj, V.G.; Bauer, P.; Boussetta, S.; Browne, P.; Deconinck, W.; Gaudin, W.; Hadade, I.; et al. A Baseline for Global Weather and Climate Simulations at 1 km Resolution. J. Adv. Mod. Earth Sys. 2020, 12, e2020MS002192. [Google Scholar] [CrossRef]
  29. Pudykiewicz, J.A.; Clancy, C. Convection experiments with the exponential time integration scheme. J. Comp. Phys. 2022, 449. [Google Scholar] [CrossRef]
  30. Peixoto, P.S.; Schreiber, M. Semi-Lagrangian Exponential Integration with application to the Rotating Shallow Water Equations. SIAM J. Sci. Comp. 2019, 41, B903–B928. [Google Scholar] [CrossRef]
Figure 1. Vorticity field at model level 20 for 12-day forecasts with four integration schemes. Top left: EuSI. Top right: EuLT. Bottom left: LaSI. Bottom right: LaLT. All runs had time step 5 min. Units 10 5 s 1 .
Figure 1. Vorticity field at model level 20 for 12-day forecasts with four integration schemes. Top left: EuSI. Top right: EuLT. Bottom left: LaSI. Bottom right: LaLT. All runs had time step 5 min. Units 10 5 s 1 .
Meteorology 01 00023 g001
Figure 2. RMS errors for six integrations of Kelvin waves with m = 1 . Red: LaSI scheme. Blue: LaLT scheme.
Figure 2. RMS errors for six integrations of Kelvin waves with m = 1 . Red: LaSI scheme. Blue: LaLT scheme.
Meteorology 01 00023 g002
Figure 3. RMS errors for six integrations of Kelvin waves with m = 4 . Red: LaSI scheme. Blue: LaLT scheme.
Figure 3. RMS errors for six integrations of Kelvin waves with m = 4 . Red: LaSI scheme. Blue: LaLT scheme.
Meteorology 01 00023 g003
Figure 4. Five-day wave RO(1,2). RMS errors of surface pressure (hPa) for 5 day forecasts.
Figure 4. Five-day wave RO(1,2). RMS errors of surface pressure (hPa) for 5 day forecasts.
Meteorology 01 00023 g004
Figure 5. Five-day wave RO(1,2). RMS errors of 250 hPa vorticity (units 10 6 s 1 ) for 5 day forecasts.
Figure 5. Five-day wave RO(1,2). RMS errors of 250 hPa vorticity (units 10 6 s 1 ) for 5 day forecasts.
Meteorology 01 00023 g005
Figure 6. Rossby-Haurwitz wave RH(4,5). RMS errors of surface pressure (hPa) for 6 day forecasts.
Figure 6. Rossby-Haurwitz wave RH(4,5). RMS errors of surface pressure (hPa) for 6 day forecasts.
Meteorology 01 00023 g006
Figure 7. Rossby-Haurwitz wave RH(4,5). RMS errors of 250 hPa vorticity (units 10 6 s 1 ) for 6 day forecasts.
Figure 7. Rossby-Haurwitz wave RH(4,5). RMS errors of 250 hPa vorticity (units 10 6 s 1 ) for 6 day forecasts.
Meteorology 01 00023 g007
Figure 8. RMS errors in surface pressure for Test Case 5 (Williamson et al. [25]). Red lines: LaSI. Blue lines: LaLT.
Figure 8. RMS errors in surface pressure for Test Case 5 (Williamson et al. [25]). Red lines: LaSI. Blue lines: LaLT.
Meteorology 01 00023 g008
Figure 9. RMS errors in surface pressure for baroclinic wave (Polvani et al. [20]). Red lines: LaSI. Blue lines: LaLT.
Figure 9. RMS errors in surface pressure for baroclinic wave (Polvani et al. [20]). Red lines: LaSI. Blue lines: LaLT.
Meteorology 01 00023 g009
Figure 10. Real data: rms error for surface pressure (hPa) over 3 days for LaSI forecasts (red) and LaLT forecasts (blue) with time steps 10, 20 and 40 min.
Figure 10. Real data: rms error for surface pressure (hPa) over 3 days for LaSI forecasts (red) and LaLT forecasts (blue) with time steps 10, 20 and 40 min.
Meteorology 01 00023 g010
Figure 11. Real data: rms error for 500 hPa vorticity (units 10 6 s 1 ) over 3 days for LaSI forecasts (red) and LaLT forecasts (blue) with time steps 10, 20 and 40 min.
Figure 11. Real data: rms error for 500 hPa vorticity (units 10 6 s 1 ) over 3 days for LaSI forecasts (red) and LaLT forecasts (blue) with time steps 10, 20 and 40 min.
Meteorology 01 00023 g011
Table 1. Four numerical schemes for Peak. The Eulerian schemes, EuSI and EuLT, are described in earlier work, [5,6]. The Lagrangian schemes, LaSI and LaLT, are described in this paper.
Table 1. Four numerical schemes for Peak. The Eulerian schemes, EuSI and EuLT, are described in earlier work, [5,6]. The Lagrangian schemes, LaSI and LaLT, are described in this paper.
Semi - Implicit Adjustment Laplace Transform Adjustment
Eulerian advection EuSI ( Eh 12 ) EuLT ( HL 19 )
Lagrangian advection LaSI ( here ) LaLT ( here )
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lynch, P. A Lagrange–Laplace Integration Scheme for Weather Prediction and Climate Modelling. Meteorology 2022, 1, 355-376. https://doi.org/10.3390/meteorology1040023

AMA Style

Lynch P. A Lagrange–Laplace Integration Scheme for Weather Prediction and Climate Modelling. Meteorology. 2022; 1(4):355-376. https://doi.org/10.3390/meteorology1040023

Chicago/Turabian Style

Lynch, Peter. 2022. "A Lagrange–Laplace Integration Scheme for Weather Prediction and Climate Modelling" Meteorology 1, no. 4: 355-376. https://doi.org/10.3390/meteorology1040023

APA Style

Lynch, P. (2022). A Lagrange–Laplace Integration Scheme for Weather Prediction and Climate Modelling. Meteorology, 1(4), 355-376. https://doi.org/10.3390/meteorology1040023

Article Metrics

Back to TopTop