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

: 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 conﬁrm 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.


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 semiimplicit 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. 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.

Laplace Transform Adjustment
Eulerian advection 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.

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].

Outline of the LaSI Scheme
We consider a general equation written in Lagrangian form where LX 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 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 where the right hand side may be computed from known quantities. We now convert Equation (2) to spectral form, multiplying by Y m n (λ, φ) and integrating over the sphere to get n is the vector of spherical harmonic coefficients of X, and n = n m n 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.

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: 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 − LL so that the term LX = L {LX} becomes L X + Γ(X), and write the transformed equation For computation of the commutator, see Appendix A.
We now convert to spectral form, integrating over the sphere, to get where x = x m n is the vector of spherical harmonic coefficients of X, x − D are the coefficients of X − D , and γ = γ m n and n = n m n are the coefficients of Γ and N. The solution may be written where, after inverse transformation to the time domain, the right hand side may be computed from known quantities.

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 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 where X − is the solution at time (n − 1)∆t. The SI approximation to (3) is Solving for the new value X + , we have 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(2iω∆t), while for the SI solution (5) the multiplier is 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 The solution for X follows immediately: .
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.

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 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 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 (dN/dt) n . Assuming that the nonlinear term varies slowly, we have dN/dt = 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.

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: 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).

Lagrangian Departure Points and Values
The Cartesian coordinates of a grid point at latitude φ and longitude λ are 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 (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].

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 /RT is, by definition, the value for an isothermal hydrostatic atmosphere. Thus, ∇π = −(1/RT)∇Φ S . We can then write The term F π involving the gradient of orography is computed in an Eulerian manner, and the continuity Equation (12) becomes The advected variable π is much smoother than the original variable π since the underlying orography has been extracted.

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: 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 The factor s 2 is included here to ensure that Γ is independent of s (see Appendix A). The transformed equations may be written Using the equations for T and π , these quantities are eliminated from the divergence equation to obtain an equation for a single variable, δ: We transfer all terms involving δ to the left:

Transforming to Vertical Eigenmodes
We now define the vertical structure matrix Since B is symmetric, its eigen-structure can be written We now transform to vertical eigen-modes by multiplying (21) by E T : 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 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: where the vectors A, B and C are Multiplying (22) by the inverse of the diagonal matrix s 2 I + Ω 2 , the equation for the k-th component is 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: The filter response function H(ω) was defined in Equation (7) above.

Inverse Vertical Transformation
Let us define four diagonal matrices (Λ D will be needed below). Then (25) can be written We can now calculate the divergence at the advanced time, For compactness, we define the propagation matrices: (P D will be used below). Then we can write the solution as The P-matrices can be pre-computed and stored, since they do not depend on the model variables.

Temperature and Pressure
We return to Equations (19) and (20): 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 Both (28) and (29) require computation of This term involves a convolution integral that may be approximated by the trapezoidal rule 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 We divide (24) by s to give We invert this using standard results for Laplace transforms to obtain Then using the Λ-matrices, we can write Noting (31) and using the P-matrices, we can now write Finally, using (28) and (29), the values of T + A and π + A are

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 ζ ) 0 M /s Dividing by s and applying L * at 2∆t yields 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.

Diffusion Stage
The governing equation for a spectral component of any of the variables δ, ζ or T may be written in the form dΨ dt = F − ςΨ We split the right hand side into two parts and integrate them separately. We first integrate the inviscid equation dΨ dt = 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 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 In the spectral domain, the damping coefficient becomes We note that the spectral equations are unchanged in form by the addition of the hyperdiffusion 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.

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 T85. 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.

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.

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.

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.

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.

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.

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.

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.

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).

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

EuSI
Eulerian advection and semi-implicit adjustment scheme EuLT Eulerian advection and Laplace transform adjustment scheme LaSI Lagrangian advection and semi-implicit adjustment scheme LaLT Lagrangian 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 . 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 0 P and C M starting from x 0 M . 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 .

Two Laplace Transforms
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 where the matrix L has an orthogonal eigenvector matrix E with LE = 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 We note that We invert this at time t n+1 = t n + ∆t to get X n+1 = e Lt n+1 X n + e Lt n+1 t n+1 t n e −Lτ N(τ) dτ (A10) 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 Lt n+1 X n + e Lt n+1 t n+1 t n e −Lτ dτ N n = e Lt n+1 X n + (−L) −1 [I − exp(L∆t)]N n . Assuming a small time-step, this reduces to X n+1 = e Lt n+1 X n + ∆tN 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 = (3N n − N n−1 )/2. Many other possibilities exist.