6. Filip VANA : "The dynamical and physical control of kinetic energy spectra in a NWP spectral semi-Lagrangian model"

The PhD thesis will be defended in Prague in September. Here is the corresponding extended abstract.

Semi-Lagrangian advection scheme with controlled damping - an alternative way to nonlinear horizontal diffusion in a numerical weather prediction model

Introduction

Horizontal diffusion schemes have been common features of numerical models since the beginning of NWP. The reason for keeping them in models is to maintain a balance of kinetic energy in the simulated atmosphere between its generation through conversion of available potential energy and its dissipative transformation into thermal energy, in order to be in agreement with so-called "turbulence theories".

The typical representation of the horizontal diffusion contribution is through the KÑ r operator. In case the aim of the horizontal diffusion is to maintain suitable kinetic energy and enstrophy spectra, it is sufficient for the diffusion coefficient K to be a constant (Koshyk and Boer, 1995). The same condition for K was found when the horizontal diffusion is seen as a self-corrector of the model (Jakimow et al., 1992). In contrary to both above interpretations, when the horizontal diffusion scheme is seen as a kind of physical parameterisation for horizontal turbulence and molecular dissipation, the diffusion coefficient K should be made flow-dependent (Sadourny and Maynard, 1997).

The spectral representation of the linear horizontal diffusion operator is very efficient and allows infinite possibilities for its tuning (Sardeshmukh and Hoskins, 1984). Similarly, in grid-point models, the linear diffusion can also be treated in a very efficient way (Li et al., 1994, McDonald, 1994). Once a non-linear operator is required to represent the horizontal damping processes, its spectral or grid-point representation becomes relatively expensive. The problem is further complicated by the fact that such a diffusion scheme is generally only conditionally stable.

The aim of this work is to propose, introduce and validate a stable non-linear scheme for horizontal diffusion in any model with semi-Lagrangian advection. The method of this alternative technique is based upon the idea of controlling the degree of interpolation needed for this advection technique (thus the scheme will be hereafter referred as SLHD - semi-Lagrangian horizontal diffusion scheme). Most of the tuning and validation work has been done within the framework of the spectral limited area model ALADIN.

Formulation of the SLHD scheme

The general form of an adiabatic model equation for a prognostic variable Y discretised for example by a three-time-level semi-Lagrangian scheme and following the space averaging of the right hand side R proposed by Tanguay et al. (1992) can be written as :

./F_Vana_Eq1.gif

The terms in square bracket marked as I have to be evaluated at point x -2a . This point, known as the origin point, denotes the point from which a particle starting at time t -Dt arrives at time t +Dt to the grid point x. Since the origin point is generally off the model grid, some interpolation method has to be used for its evaluation. The interpolator used to be chosen in order not to be too expensive and still accurate enough with respect to the other model simplifications. The compromise is typically a cubic polynomial type interpolator (Staniforth and Côté, 1991).

The semi-Lagrangian diffusion scheme is then simply defined by evaluating the interpolator in equation (1) as a combination of two interpolators :

I =(1-k ) IA +k I D

=IA + k (ID - IA) (2)

Here IA is an interpolator of sufficient accuracy for use in the semi-Lagrangian scheme. The interpolator I D is defined in such a way as to provide a larger degree of damping than IA . Hence the difference I D - IA can be interpreted as an additional damping to the original semi-Lagrangian scheme, introduced by the interpolator ID . Following the same logic, the adimensional parameter k can then be seen as the coefficient controlling the degree of this damping.

To be able to use the additional damping of the semi-Lagrangian interpolator I defined in equation (2) as horizontal diffusion scheme, the parameter k has to be defined according the theory of the fluid dynamics. Sadourny and Maynard (1997) express the sub-grid scale contribution of the horizontal turbulence and molecular dissipation as function of the divergence and of the deformation of the flow. The latter is then fairly more important for the fluids being simulated by atmospheric models. Similarly to this, the parameter k was defined as function of the scalar quantity representing the total deformation d, defined as :

./F_Vana_Eq3.gif

where u and v are the components of the horizontal flow. To avoid any dependency of such a diffusion to the integration timestep Dt, the parameter k was defined in the following way :

./F_Vana_Eq4.gif

where F(d) stands for any monotonic function of d. In order to make the semi-Lagrangian diffusion scheme also independent of the model resolution, the function F(d) was defined as :

./F_Vana_Eq5.gif

Here a, b and d0 are tunable parameters for the proposed diffusion scheme while the ( Dhref / Dh) P term ensures the quasi-independence with respect to the chosen model resolution. Dh is a general parameter representing the horizontal model resolution (which is proposed like Dº Ö (D x 2   +  D y 2 ) for grid-point models or D º Ö (k x-2   +  k y-2) for spectral models), Dhref its reference value and P a tunable parameter generally dependent on b. Equations (2) - (5) with the tunable parameters a, b , d0 and P then fully define the proposed non-linear semi-Lagrangian diffusion scheme.

Idealised frontogenetic simulation by a NWP model

To test the skills and properties of the SLHD scheme the academic 3D Eady wave frontal development was introduced into the spectral limited area numerical weather prediction model ALADIN. Such a situation is very sensitive to the level of model dissipation (MacVean, 1983) and thus tends to be an ideal tool for testing any diffusion scheme in the model. Another advantage of this special situation is its capability to study clearly the impact of model diffusion without harmful feedbacks being generated during the moist diabatic simulation of a "full model".

Otherwise the model used for the simulations was kept as close to a typical NWP operational configuration as possible (resolution of current NWP applications, boundary conditions, coupling and so on. ...). The only differences occurred in the model geometry (ideal plane, constant Coriolis parameter) and in the definition of simulated atmosphere (dry, no diabatic forcing allowed).

The simulation was performed with the model periodic in the zonal (y) direction. The size of the model domain was defined as 3600 km along the y direction times 9000 km along the x direction with the horizontal resolutions Dx = 20 km and Dy = 50 km. Once the non-linear processes start to play a sufficiently important role in the formation of the front, the simulation is repeated on a nested domain with finer uniform horizontal mesh Dx = D y = 10 km. The distribution of the 41 vertical levels is kept the same for both simulations with the intention to fulfil the theoretical ratio between the horizontal and vertical resolution in the bottom most frontogenetic part of the atmosphere in the nested run. Like that the frontogenetic process is prevented from the generation of spurious destroying gravity waves of numerical origin (Persson and Warner, 1991).

The upper boundary of the frontogenetic area was treated by introducing a slanted tropopause and activating enough horizontal diffusion (the spectral one) for the levels above it.

Even with the nested model, the frontal gradient was not strong enough to impose the model collapse. This is fully in agreement with Garner (1989) concluding that the model horizontal resolution must be at least of order of 1 km to be harmed by numerical instability generated from a simulated frontogenesis. Consequently, the simulated kinetic energy growth rate was at a certain time deviating from its asymptotical theoretical value and the further growth of frontogenetic gradients was stopped, which is also fully in agreement with observations of Gall et al. (1987).

Since the intention was to study the academic experiment with a model similar to any typical NWP application, there was no space to create the model numerical instability causing at a given time the collapse of the model. Thus, instead of measuring the time of the model collapse as a function of level of horizontal diffusion, a special statistical tool has been developed to diagnose the impact of any horizontal diffusion scheme acting in this model. The diagnostic simulation was then defined as the 12 hour simulation of the nested model before the maximum gradient was reached. This diagnostic was found sensitive enough to any horizontal diffusion scheme used in the model.

As diagnostic output, the kinetic energy spectrum of the lowest (most frontogenetic) model level has then been used.

Method for computing the parameters of an equivalent diffusion

As a measure to quantify the observed diffusivity of any diffusive feature within the academic diagnostic simulation, a statistical tool was introduced. Its aim is to fit the unknown diffusivity by a spectral linear diffusion with the closest possible properties. This simple approach allows to express any model damping scheme by the well-describing parameters of a linear diffusion scheme: the coefficient of diffusion K and the order of diffusion r . Since the non-linear effects of horizontal diffusion start to be important with high model resolution (around 1 km of horizontal mesh) it is believed that the performance of the proposed SLHD scheme can be fitted by the linear diffusion without any loss of relevance in the framework of the academic simulation running at 10 km of horizontal mesh.

The principle of the proposed method is to decompose the resulting kinetic energy spectra of the idealised simulation as sums of the kinetic energy spectra of simulations without any horizontal diffusion scheme and of a remaining part proportional to the used diffusion. It is known that the logarithm of kinetic energy spectra is close to a linear function of the logarithm of wave-number. Hence the statistical model has been built on logarithms of the kinetic energy spectra, rather than just on kinetic energy spectra themselves. Once the spectral coefficients of the kinetic energy have been recomputed from the kinetic energy values for the elliptic bands corresponding to the zonal wave-numbers of the model, the logarithm of such spectra can be expressed as :

./F_Vana_Eq.gif

Here y0 (m) represents the resulting spectrum of the model without diffusion scheme, Dy( m) the influence on log(y (m)) caused by the tested diffusion scheme and M is the total truncation of the model. The statistical model has thus been built as the linear regression of all combinations of parameters of spectral linear diffusion r and H up to third order :

./F_Vana_Eq6.gif

The parameter H is derived from the horizontal diffusion coefficient K which is, for the ALADIN model, with a given truncation and at any given model level fully determined by just r and H .

The values of the parameters a ij (m) were then set according to the spectrum of 150 experiments with spectral linear diffusion with known parameters r and H. The statistical model correlation was for all wave-numbers at least 0.9, which allows to use this model in the following with confidence in its relevance.

Knowing the values of a ij (m), the properties of an unknown diffusion can then be characterised by most closely fitting the parameters of spectral linear diffusion r* and H* by solving the following minimising formula :

./F_Vana_Eq7.gif

Behaviour of the SLHD scheme in the ALADIN model

The SLHD scheme defined by equations (2) - (5) is not fully specified unless the interpolators IA and ID are defined. When introducing this scheme into the model ALADIN, its standard semi-Lagrangian interpolator (Ritchie at al., 1995) was set as the accurate interpolator IA. The linear interpolator was then used as the damping interpolator ID.

To study the property of the SLHD scheme, 90 experiments with different sets of tunable constants a, b and d 0 were launched in the framework of the idealised 3D frontal development. The damping characteristics of this scheme then were estimated by the best fitting parameters of spectral diffusion r* and H*. The results of 66 from all, with r* and H* falling inside the values of the statistical model were then plotted in the H -r plane.

As expected the strongest diffusion obtained by the SLHD scheme corresponds to the damping characteristics of the interpolator ID. Weakening the diffusion by allowing an increasing proportional participation of the interpolator I A to the resulting semi-Lagrangian interpolator I, the characteristic points on the H -r plane were distributed along a specific curve starting at the characteristic point of ID. Further diagnostics showed that the main factor influencing the position of the characteristic points on the H-r plane is the average value of the k coefficient. This implies that there is just one degree of freedom between the characteristics r and H of the SLHD scheme. Thus, from the damping point of view, just one tunable constant would be sufficient to control the properties of the SLHD scheme.

Another important outcome from the above is the fact, that the SLHD scheme is, within usable model damping rates, always less selective than the operational 4th order diffusion.

When testing performances in the framework of real atmospheric simulations, the following two main differences were found with respect to the academic case behaviour.

First, in order to obtain the same amount of kinetic energy as with the spectral diffusion with the characteristics of the best fitting parameters r*and H* measured in the framework of the academic experiment, the SLHD scheme has to be tuned to much stronger values. The reason for this comes from the dependency of the properties of the SLHD scheme upon the average of the field of total deformation d of the simulated atmosphere. Due to the very specific case simulated by the academic test, the values of d were about one order of magnitude greater than the ones reached in the real atmospheric simulation. Consequently the SLHD diffusion tuned for the extreme values of d in the academic simulation would act in a much weaker way in the real atmosphere characterised by lower values of the d field.

Beside this first difference, which seems to be just a problem of tuning, the second factor responsible for a different behaviour of the SLHD scheme in real atmospheric conditions from that in the academic environment seems to be of a more systematic nature. As seen from the equation (1) the SLHD diffusion modifying the semi-Lagrangian interpolators affects just one part of the value for the prognostic field Y (x, t +Dt ). The remaining part unaffected by this diffusion is represented by the last term D t ×R( x, t ) on the right hand side. Analysing this residual for each model equation, one can find that the prognostic equation for the components of the horizontal wind contains the derivatives of the model orography. This field, obviously constant during the model integration, has been found as a source of noise that couldn't be treated by the SLHD scheme. To remove such a weakness of the proposed diffusion scheme two solutions have been proposed, from which none is ideal. First the model orography could be filtered in such a way to be in equivalent spectral slope to the expected kinetic energy spectral tail. The second solution is to keep in the model another (weak and selective) diffusion scheme controlling the tail of the kinetic energy spectra together with the SLHD one.

Since there are several numerical reasons to keep a spectral diffusion in the ALADIN model even if the physical contribution of the horizontal dissipation would be represented by another device on top of SLHD, the second solution to control the small-scale kinetic energy waves seems to be clearly preferable.

Aspects of the implementation of the SLHD scheme into the ALADIN model

On top of the described parameters controlling the performance of the SLHD diffusion scheme, some other technical parameters have been coded during the implementation of the scheme into the ALADIN NWP model. One of them works as a limitation for the SLHD "coefficient of diffusion" k prescribing a maximum value for this parameter. The reason for its implementation into the source code is the opportunity to guarantee a minimum amount of interpolation by the accurate interpolator IA. Another technical parameter specifies the additional diffusion that can be used for the boundary points for which the semi-Lagrangian trajectory originates outside the model domain.

Since the SLHD scheme depends on the used interpolators IA and ID , the spectral sensitivity of the different interpolators available in the ALADIN code - the cubic Lagrangian, the four-point spline and the linear interpolator - has been investigated within a special 1D experiment. From there the specific sensitivity of the four-point natural cubic spline has been detected. This interpolator was found to behave similarly to the linear interpolator for the long waves, with one order higher accuracy. There the cubic Lagrangian interpolator sensitivity to the wave number is different and the accuracy is about two orders higher than for the four-point spline. Once the waves from the middle of the spherical expansion with the quadratic truncation (with wavelength about twice and less the minimal wavelength represented by the last term of the expansion) are interpolated, the four-points spline starts to perform interpolation with higher accuracy as compared to the other two tested interpolators. This trend is kept till the smallest scales represented by the model spectrum, but the differences in accuracy are decreasing there for all methods. Such a result seems to be promising for giving the possibility to impact the waves from the middle of the spectral expansion if introducing the four-point spline interpolator as the accurate interpolator IA . Unfortunately there is not much space to impact the small-scale end of the model spectra by any of these interpolators.

Finally the coexistence of the spectral horizontal diffusion and of the SLHD scheme in the model is described. As already mentioned, when the SLHD scheme, acting like the model physical diffusion, is activated, there is still need to keep the spectral diffusion to remove the numerical noise. For such a case the spectral diffusion is kept untouched in the high atmosphere where it plays an important role as the damping mechanism for gravity waves of numerical origin. In the low atmosphere where the numerical reason for keeping it in the model is just its role to remove the small scale noise caused by orography in the fields of components of wind, this linear diffusion is weakened and made very selective, if not suppressed completely in the case of temperature and moisture fields for example.

Validation of the SLHD scheme

Once implemented in the ALADIN model, the new horizontal diffusion has been validated. Two kinds of validations have been proposed: the parallel run with the spectral horizontal diffusion for a random period and the verification of the scheme's ability to act in specific case studies sensitive to the level of horizontal diffusion in a model.

For the parallel test, the three configurations of horizontal diffusion have been used, running with the ALADIN/LACE operational domain (Janousek, 1999). As reference run, the spectral linear horizontal diffusion representing the current operational scheme was taken. The SLHD diffusion scheme with the additional weak and selective spectral diffusion controlling the tail of the model spectra was used as the second test. Finally the third test used just the SLHD diffusion scheme in the troposphere with the orography field filtered by a relatively strong 5th order spectral diffusion (the last wave-number had been reduced to 40 % of its original value). As testing period, the 20 days between March 3rd and March 22nd 2000 were randomly chosen.

The evolution of scores (BIAS and RMSE) behaved for all three runs in very similar way during the testing period. The differences were very small if ever detectable. This is not surprising with the 12.18 km of horizontal resolution of the ALADIN/LACE model, which is far from the scales where the non-linear effects in horizontal dissipation are of theoretical importance. Anyway a few deviations have been detected between the SLHD and the spectral diffusion scheme. First the BIAS of the MSL pressure is systematically worsening during the model simulation for the SLHD scheme compared to the reference. This feature is a consequence of the local character of the grid-point diffusion that is not conserving the mean of the diffused fields contrary to the globally acting spectral diffusion. On the other hand the SLHD scheme is formulated in a more accurate way with respect to the orographic features than the spectral diffusion acting along the quasi-horizontal terrain-following model levels. Hence the RMSE of the MSL pressure is nearly the same for all three tests. The lower conservative properties of the grid-point SLHD are compensated by the more accurate formulation of horizontal dissipation over the orographic features. Second a slight tendency of cooling the atmosphere is detected in the areas where the SLHD scheme acts as the main damping scheme. This (positive) impact on the model scores is probably a consequence of the 3D formulation of the SLHD diffusion, which is smoothing as well along the vertical and thus completing the vertical diffusion as parameterised by the model physics. The role of controlling the temperature dissipation through the flow deformation when the SLHD scheme is acting might be important as well.

The above differences were highlighted by the test of the SLHD diffusion without spectral supplement in the low atmosphere. Surprisingly, this run was not resulting in bad scores for the upper atmosphere. The smoothing of the model orography affected just the near surface variable (mainly wind field) and the geopotential in the low troposphere (between surface and 700 hPa).

Contrary to the parallel tests meant to detect some systematic tendencies, case study simulations can show the ability of a studied scheme with respect to situations for which it is specially suited. As shown for example by MacVean (1983) and Hello et al.(1999), cyclogeneses simulated by a numerical model are sensitive to the model diffusion scheme. To test the benefits of the proposed SLHD scheme several cyclogenetic situations as well as events on very small scales have been studied. Most of them showed no impact with respect to the spectral diffusion scheme used as the reference. The reason was mainly due to the fact that cyclogenesis is a consequence of the large scale forcing, thus its simulation in the limited area model is extremely dependent upon the initial and lateral boundary conditions. However a few cases of very local character or from meso-cyclonic development and sensitive to the model horizontal diffusion scheme could be found, even when using the simulation within a relatively small domain like for the operational ALADIN models.

The first presented case study concerns the storm from December 26th 1999. This time the strong meso-cyclone was created just to the west of the French Atlantic coast. The cyclone was moving eastward with a speed of about 100 km/h. This event was very well forecast by the operational ALADIN applications covering the affected area, as well as by the coupling model of ALADIN - the French global model ARPEGE. The simulation was provided with the ALADIN/LACE operational domain. The cyclone enters this domain already developed but continuing its development. Both diffusion schemes (SLHD and spectral diffusion) produced very similar results. The maximum difference between the two simulations appeared 9 hours after the cyclone entered the model domain. The difference in the MSL pressure was less than 2 hPa. After another 12 hours of integration the minimum value of MSL pressure in the centre of the cyclone was about the same in both cases.

The case just mentioned exhibits that the two different diffusion schemes can act very similarly for an event that has been successfully forecast by the model in both cases.

The ensuing case was, contrary to the previous one, unrealistically exaggerated in the operational configuration of ALADIN/LACE. The small meso-cyclone appearing during July 20th 2001 in the Adriatic sea approached the Dalmatian cost with a lowest value of MSL pressure of about 1006 hPa. The ALADIN/LACE forecast with the spectral diffusion simulates this event with a lowest MSL pressure value of 990 hPa. Of course, such a low was in the model simulation accompanied by extreme weather events, which didn't occur in reality. When applying the SLHD scheme to this forecast the depth of the low was reduced to 1004 hPa. Further analysis of the model simulation showed that this kind of locally acting diffusion prevents a strong and unrealistic stratiform precipitation generation in the model behind the moving cyclone. The analysis of the potential vorticity in the atmosphere clearly exhibited that for this case the SLHD scheme with an equivalent strength to the operational spectral diffusion significantly improved all the tropospheric part of the atmosphere in the vicinity of the simulated low.

The last presented situation concerns a small-scale event again not really well simulated by the ALADIN/FRANCE operational forecast. During the August 24th 2001 12 UTC simulation, 3 hours after the model start an extreme and unrealistic downburst developed above the Atlantic ocean westwards from Gibraltar. This event was of very local character (horizontal scale of about 100 km). Even if the reason for this event can be traced to some other part of the model (mainly the stratiform precipitation), the (spectral) horizontal diffusion scheme was felt to act too weakly. When the SLHD scheme of similar impact on the model kinetic energy spectra was used to recompute this situation, the descending motions were reduced by 35 % . A similar impact could be reached with the spectral horizontal diffusion tuned to be eight times stronger than the operational one. But this re-tuning of course affected drastically the whole kinetic energy spectrum (and thus as well the quality of the model forecast), which was not the case for the locally acting SLHD diffusion.

Conclusion

A new alternative way to the horizontal diffusion treatment in numerical models using semi-Lagrangian advection has been proposed. The scheme was diagnosed in the framework of a specific academic testing tool specially developed for the purpose of this work. Then the behaviour of the scheme was studied in the real atmosphere. Finally the scheme was validated with respect to the spectral linear diffusion scheme in a parallel test as well as for several specific case studies.

The properties of the proposed diffusion scheme are of such nature that they make this scheme to be advantageous mainly for the numerical weather prediction applications. The new scheme offers the following pleasant technical aspects for horizontal diffusion in a model. First to incorporate such scheme into a semi-Lagrangian model is from the technical point of view very easy. Moreover, except the semi-Lagrangian model advection, the SLHD diffusion implies no other restriction for the model dynamics. Second the computational cost of the SLHD scheme is mainly limited by the time spent for additional interpolations of the semi-Lagrangian amounts by the diffusive interpolator ID which is typically a cheap (tri-)linear interpolator. Similarly to the spectral diffusion scheme, the computational cost of the SLHD scheme is fixed with respect to its chosen selectivity. Third the proposed diffusion scheme is unconditionally stable in the sense of damping characteristics. Its stability restriction is given by the stability criteria for the semi-Lagrangian advection. Consequently the model timestep is generally not more restricted when the SLHD scheme is activated in a model than without it.

For spectral models the SLHD scheme offers other clear advantages compared with the spectral linear diffusion schemes typically used there. First the SLHD scheme proposes a stable and efficient way for a non-linear and thus more realistic description of the effects of horizontal dissipation processes. Although it is believed that non-linear effects don't play any important role on horizontal scales like 10 km and bigger, it has been demonstrated that for some extreme situations they can be of specific importance. However with the increasing model resolution, the importance of non-linear horizontal dissipation will be constantly growing. Second the grid point formulation of the SLHD scheme allows to use this scheme for representing the horizontal damping of those variables which are not passing through spectral space (i.e. atmospheric liquid water and ice, pollutants, ...). Third the grid point character of the SLHD diffusion scheme enables efficiency and accuracy when forcing diffusion to be homogeneous over a variable mesh, by making it easily resolution dependent.

As usual, the proposed scheme brings also some disadvantages when comparing it to the separate sophisticated diffusion schemes. The main two of them are as follows. The way by which the proposed SLHD scheme makes the spectra of kinetic energy (or the horizontal components of the flow field) dependent on the small scale spectra of the model orography. To be able to control these model variables therefore requires either to significantly filter the model orography or to add another horizontal diffusion scheme removing the small-scale noise of the kinetic energy spectra. The second limitation of the proposed SLHD scheme is its limited tuning possibilities. In case the influence of a diffusive scheme is characterised by the coefficient of diffusion and the order of diffusion, the SLHD scheme makes one of them to be fully expressed as a function of the other. To change the property of a given SLHD characteristic function is indeed possible just by changing the semi-Lagrangian interpolators. This is from a practical point of view rather uncomfortable.

The proposed SLHD scheme however seems to be a viable alternative to the currently used diffusion schemes in numerical models. Its benefit arises mainly from the fact that it reverts the known disadvantage of the semi-Lagrangian advection - the additional non-excessive damping - to an advantageous control of the complex separation between the model signal keeping predictive skills and the noise.

References

R. L.  Gall, R.  T.  Williams, and T.  L. Clark : On the minimum scale of surface fronts. J. Atmos. Sci., 44 (18), 2562-2574, September 1987.

S. T.  Garner : Fully Lagrangian numerical solutions of unbalanced frontogenesis and frontal collapse. J. Atmos. Sci., 46 (6), 717-739, March 1989.

G. Hello, F.  Lalaurette, and J.-N.  Thépaut : Combined use of sensitivity information and observations to improve meteorological forecast: A feasibility study applied to the "Christmas storm" case. Quart. J. Roy. Meteor. Soc., 126, 621-647, 2000.

G. Jakimow, E. Yakimiw, and A.  Robert : An implicit formulation for horizontal diffusion in gridpoint models. Mon. Wea. Rev. , 120, 124-130, January 1992.

M. Janousek : National status report of Czech Hydrometeorological Institute. In LAM Newsletter , number 29, pages 55-56, 1999.

J. N.  Koshyk and G.  J.  Boer : Parameterization of dynamical subgrid-scale processes in a spectral GCM. J. Atmos. Sci., 52 (7), 965-976, April 1995.

Yong Li, S.  Moorthi, and J.  R.  Bates : Direct solution of the implicit formulation of fourth order horizontal diffusion for gridpoint models on the sphere. Tech. Memo 109466 Vol. 2, NASA, December 1994.

M. K.  MacVean : The effects of horizontal diffusion on baroclinic development in a spectral model. Quart. J. Roy. Meteor. Soc., 109, 771-783, 1983.

A. Mc  Donald : Using second, fourth and sixth order implicit horizontal diffusion to control noise in three dimensional semi-{L}agrangian, semi-implicit, limited area, gridpoint models of the primitive equations. Technical report, Irish Meteorological Service, 1994. Unpublished note.

P. O.  G.  Persson and T.  T. Warner : Model generation of spurious gravity waves due to inconsistency of the vertical and horizontal resolution. Mon. Wea. Rev., 119, 917-935, April 1991.

H. Ritchie, C. Temperton, A.  Simmons, M.  Hortal, T.  Davies, D. Dent, and M. Hamrud : Implementation of the semi-Lagrangian method in a high-resolution version of the ECMWF forecast model. Mon. Wea. Rev., 123, 489-514, February 1995.

R. Sadourny and K. Maynard : Formulations of lateral diffusion in geophysical fluid dynamics models. In Numerical methods in atmospheric and oceanic modelling - The André J. Robert memorial volum , pages 547-556. NRC Research Press, 1997.

P. D.  Sardeshmukh and B.  J.  Hoskins : Spatial smoothing on the sphere. Mon. Wea. Rev. , 112, 2524-2529, December 1984.

A. Staniforth and J. Côté : Semi-Lagrangian integration schemes for atmospheric models - a review. Mon. Wea. Rev., 119, 2206-2223, September 1991.

M. Tanguay, E. Yakimiw, H.  Ritchie, and A.   Robert : Advantages of spatial averaging in semi-implicit semi-Lagrangian schemes. Mon. Wea. Rev., 120, 113-123, 1992.