Test of a cubic spline interface for physical processes with a 1-D third-order spectral element model

Abstract A common way to introduce physical processes into numerical models of the atmosphere is to call the parameterization at every grid point. This can lead to considerable errors. A simple 1-D example is proposed to illustrate that when a physical process occurs at one grid point only, a considerable sampling error may occur, with the result that only a fraction of the true impact of this process is seen. The interface to the physical parameterization in numerical weather prediction model using a third-order 1-D spectral element method (SEM3) model is investigated by homogeneous advection. In SEM3, the grid points, called principal nodes, are at boundaries of computational intervals and two more collocation points in the interior of each cell. This article argues that it is sufficient to do the physical parameterization for principal nodes only that creating the interior grid-point values of physics schemes by linear interpolation. This is called the spline interface method. A simple condensation model of water is taken as an example. Compared to the standard paramaterization, which computes the physical processes at every grid point, the spline interface method is more accurate and has a potential to save computer time. It turns out that the standard method creates a noisy wave which can easily be filtered by hyperviscosity. In the spline interface to the condensation physics, the condensation is done at every third grid point only. Third-order spline methods are used to represent the condensation at other points. The method using a smaller grid to compute condensation represented the condensation process more accurately and produced less of the computational noise. This version could be run without hyperviscosity, as no significant computational noise mode was generated by condensation. By doing physical processes only at every third grid point computer time may be saved.


Introduction
To represent the advection process with reasonable accuracy fields of sufficient smoothness are necessary. In virtually all realistic models of the atmosphere, the physical processes are called at all grid points (Mesinger et al. 1988;Z€ angl 2002;Pielke et al. 1992;Schluenzen 1994). Some rather expensive calculations, such as radiation (Ritter and Geleyn 1992), are not called for all time levels and grid points. Apart from this, physical interfaces, where the routines are called on a space of smaller dimension than the dynamics computations, are not much in use.
The parameterized physical processes can generate rough fields, as in many current models, the parameters determining the physical processes may be chosen in an unsmooth way. In this way, the physical parameterizations create fields which the advection processes cannot handle. Therefore, it makes sense to compute the paramaterization on a reduced grid, for example only for every third grid point.
In order to explore this idea we use a one-dimensional (1-D) model of homogeneous advection of the density field h(x) using a third-order L-Galerkin method. A summary of the L-Galerkin approach is given by Steppeler and Klemp (2017). The spectral element method (SEM, see Giraldo 2001) and the third-degree method by Steppeler (1976) are examples of L-Galerkin methods. As a test example, we use the third-order spectral element method (SEM3, Taylor et al. 1997). The approximation space for this method is the set of piecewise third-order polynomials, which fit together continuously at cell boundaries. SEM3 can equivalently be performed in grid-point space, where the grid is x i (i ¼ 0, 1, 2, 3, … ). In 1-D space for SEM3, the cells are intervals ðx i ; x iþ3 Þði ¼ 0; 3; 6 . . .Þ with a length of 3dx ¼ x iþ3 Àx i : The points x i ði ¼ 0; 3; 6; . . .Þ at cell boundaries are called principal nodes and the points x iþ1 ; x iþ2 are called interior nodes. For SEM3, x iþ1 and x iþ2 must be chosen such that the points x i ; x iþ1 ; x iþ2 ; x iþ3 ði ¼ 0; 3; 6; . . .Þ form Legendre-Gauss-Lobatto points. Therefore, the grid x i ði ¼ 0; 1; 2; 3; 4; . . .Þ is irregular when the grid x i ði ¼ 0; 3; 6; 9; . . .Þ is regular.
The Legendre-Gauss-Lobatto grid x i is defined as: where i ¼ 0; 3; 6; . . . : The test equation is homogeneous advection: where the velocity u 0 is 1. For the approximation of Equation (2) with SEM3 we refer to Taylor et al. (1997). We describe only the part of the spline method to be used to create an interface to the physics routine in this article. The aim of this article is to explore the performance of physical processes on a reduced grid for physics, consisting of every third point of the full grid. The physical process chosen as a test for the reduced grid physics calculation is condensation. A very simple condensation scheme for the 1-D model will be defined in Section 2. For smooth fields, the condensation sometimes generates artificial maxima and minima, which are treated by flux adjustment and we will not investigate such schemes in connection with SEM3. The results will be shown in Section 3 and Section 4 concludes this study.

Description of the condensation procedure
As pointed out in the introduction, it is questionable that the physics is conveniently called at every grid point, or that is reasonable to allow the parameters determining the physical parameterization to change discontinuously from one grid point to the next, as is common practice in realistic models (Steppeler et al. 2003;Damrath et al. 2000). The condensation experiment is set up to test this idea.
When interpreting h(x) in Equation (2) as the density of moisture, we assume that there is a maximum moisture amplitude for hðx i Þ: This maximum density depends on the grid point i in question. The point where condensation occurs is called the point with condensation. We assume that the density has a maximum of 2 at points i with condensation. The condensation procedure is then defined as for i indicating the point with condensation in our case and i max ¼ 600. The action of the condensation at one time step for h ¼ 4.0 in grid-point space is shown in Fig. 1. When Equation (2) is applied for both the principal nodes and the interior point of a cell, we call it the fullgrid condensation. This kind of physical interface is used in nearly all atmospheric models and is the standard physical interface. A more refined interface to physics is called cubic spline or principal-node interface. The condensation Equation (2) is applied for principal nodes at cell boundaries (i ¼ 0; 3; 6; . . .) only. It would not make sense to apply Equation (2) just to the principal nodes and do nothing to the other points, the interior nodes. The definition of the spline method is inspired by the field representation used in SEM3 or third-order finite element methods (see Steppeler 1987 for a review). The basis function spaces under the serendipity grids and the associated concept of sparse grids are described by Ahlberg et al. (1967).
To define this interface, the field is be partitioned into a piecewise linear part h lin ðxÞ and a high-order contribution h ho ðxÞ: h lin ðxÞ is defined to be linear in x in all intervals ðx i ; x iþ3 Þ and continuous at x i ði ¼ 0; 3; 6; . . .Þ: The linear part of h must be computed on the irregular Legendre-Gauss-Lobatto grid defined in Equation (1) with dx ¼ 1: where i ¼ 0; 3; 6; . . . ; x i ; x iþ1 ; x iþ2 ; x iþ3 are on the irregular Legendre-Gauss-Lobatto grid and can be calculated by Equation (1). The high-order part of h can be computed by: with the spline condensation method h lin i and h ho are computed using Equations (5) and (6). The condensation Equation (3) is applied on h lin i for i ¼ 0; 3; 6; . . . : Again Equation (5) is used to compute h lin i for i ¼ 0; 1; 2; 3; 4; . . . ; after condensation at principal nodes. The field after condensation is computed using Equation (4). As the high-order part h ho ðxÞ in Equation (2) is a small contribution. It makes sense to apply the condensation procedure to h lin ðxÞ only. Steppeler and Klemp (2017) defined initial values being smooth enough to obtain a sufficiently accurate solution in 2-D. A 1-D version of this will be applied here as initial value for Equation (1).

Results
The test used here is constructed to investigate the performance of physical parameterization schemes for smallscale atmospheric structures combined with a physical parameterization which also acts on a small scale. The test will be performed on a grid with 600 points (200 grid intervals) with dx ¼ 1 and periodic boundary conditions for Equation (1). The initial conditions are given by Similar initial values were used by Steppeler and Klemp (2017) in a 2-D version and are smooth enough to allow an accurate solution with centered difference schemes. Here we use half the width of the structure than Steppeler and Klemp (2017), as we here use the highorder method SEM3, which is able to handle smaller structures accurately. For time integration, the fourthorder Runge-Kutta (RK4) scheme is used as an example described by Steppeler et al. (2008). For the comparison of the condensation method, the naive condensation method is used in the test where the condensation Equation (3) is applied on h i ði ¼ 0; 3; 6; . . . ; i max Þ while for the comparison of SEM3, the spatial centered difference scheme is also used where the semi-discreted advection equation is where i ¼ 0; 1; 2; 3; . . . ; i max with the periodic boundary condition h 0 ¼ h imax and x i is defined in Equation (1). To appreciate the information on time stepping, it may be useful to know that RK4 with centered difference in space implies a CFL condition of 2.8 and this is reduced by a factor 0.7 when centered differences are replaced by standard fourth-order spatial differencing (Ullrich 2014). For SEM3, dt ¼ 3 2 dx is the maximum stable time step, consistent with the smallest grid length implied by Equation (1). The homogeneous advection with a periodicity length of 600dx and u 0 ¼ 1 produces the initial value after 600 dt ¼ 400 time steps. The following figures show the initial value and the forecasted fields after 80; 160; 240; 320 and 400 steps, where the last approximates the initial value. Figure 2 gives the result of advection without condensation. A small decrease of the maximum of h with time occurs due to dispersion. Figure 3 shows the result of the advection for full-grid condensation. After 80 steps, the condensation point is not yet reached and the well-known accuracy of third-order advection is achieved. At this time the solution is identical to that shown in Fig. 2. The plots at all other times show the result after condensation has been experienced, which occurs at x 300 . The correct solution should show a maximum of 2. The solution in Fig. 3 is rather incorrect in the representation of the condensation. It has a maximum value of the field larger than 3, only slightly smaller than shown in Fig. 2. This can be explained, as condensation occurs at one point only and there is a sampling error, as the solution is moving by more than dx in one time-step. The sampling error goes down for smaller time steps, but not very fast. If dt is reduced from dt ¼ 1.5 to dt ¼ 0.5, the maximum of the field is still larger than 2.8 (plot not shown). The solution in Fig. 3 is contaminated by a component of 2dx noise. This noise can be eliminated by the introduction of fourth-order hyperviscosity. A coefficient of 1 10 of its maximum reasonable value is sufficient to remove the noise (plot not shown). The effect of the hyperviscosity on the solution for pure advection without condensation is marginal. It is well known that physical processes can create small scales.
The experiments are continued without hyperviscosity. Spline condensation as described above computes condensation only at every third grid point and uses splines to describe the effect of condensation at the other points. The result is shown in Fig. 4 with SEM3. The maximum of the density field is now correctly 2. The 2dx component for this interface is absent. It appears that principal node condensation is accurate enough and even more accurate than full-grid condensation. This may potentially save a significant amount of computer time.
The experiment (Fig. 4) was reproduced in Fig. 5 with second-order spatial centered difference of Equation (8), which are less accurate than the fourth-order SEM3 scheme. The maximum of the density field is also correctly 2. However, the result shows more dispersion compared to Fig. 4. Because of the increased dispersion as compared to the higher-order SEM3 scheme, after passing the condensation point, the maximum of the density still decreases marginally. In other respects, concerning the action of the parameterization, it is comparable to that using SEM3. This suggests that the third-order interface for physics can be used also with a lower-order modelling scheme. Therefore, one may hope that this physics interface could be tested with any model, not only models using the SEM3 scheme.
While the simple condensation in our play model is computationally not expensive, realistic models may have expensive condensation schemes, involving five kinds of water phases. The interface can also be applied to other physical processes, such as Ritter and Geleyn (1992). Many forecast models (see Benoit et al. 1997 for example) use more than half of the computational effort in the physics scheme (Damrath et al. 2000). SEM are often used in the horizontal only (Taylor et al. 1997). However, applying SEM for the vertical as well would increase the sparseness factor considerably. With the SEM3 method used here, the improved physics interface involves no overhead computations, as the partition between the linear and the high-order density is already done in the SEM3 scheme.

Conclusion
For the third-order L-Galerkin spectral element method SEM3, two interfaces for condensation physics were tested. For the physical parameterization, a simple condensation method together with homogeneous advection was used. A rather unsmooth physics scheme was used, where condensation occurs at one grid point only. The standard version, when condensation is done for every grid point, encountered a large sampling error. The spline condensation method would do the condensation at principal nodes only but extending it to the linear spline part of the field. It turned out that using the full-grid method created a 2dx noise, which could easily be filtered by fourth-order diffusion. While just having the condensation at the principal nodes creates no 2dx noise. Futhermore, when applying the spline condensation method for the spatial centered difference, it is found the scalability of spline condensation method to a lower-order modelling scheme with an expected larger dispersion.
One advantage of the spline condensation method is the reduced need for fourth-order diffusion. Also, this method saves computer time, as for our 1-D case the condensation is done on only 1 3 of the node points. For our 1-D model and the very simple condensation model used here, this method saves only a small amount of time. However, applying this method to realistic models for the whole physics package could lead to substantial savings. For example, the condensation in such models is a rather complicated process involving many phases of water. The physics package in such models often takes more than half of the computer resources (see Steppeler et al. 2003 as an example) and therefore doing the physics using the spline method described above has the potential of considerable improvement of computational efficiency.