An online trajectory module ( version 1 . 0 ) for the nonhydrostatic numerical weather prediction model

A module to calculate online trajectories has been implemented into the nonhydrostatic limited-area weather prediction and climate model COSMO. Whereas offline trajectories are calculated with wind fields from model output, which is typically available every one to six hours, online trajectories use the simulated resolved wind field at every model time step (typically less than a minute) to solve the trajectory equation. As a consequence, online trajectories much better capture the short-term temporal fluctuations of the wind field, which is particularly important for mesoscale flows near topography and convective clouds, and they do not suffer from temporal interpolation errors between model output times. The numerical implementation of online trajectories in the COSMO-model is based upon an established offline trajectory tool and takes full account of the horizontal domain decomposition that is used for parallelization of the COSMO-model. Although a perfect workload balance cannot be achieved for the trajectory module (due to the fact that trajectory positions are not necessarily equally distributed over the model domain), the additional computational costs are found to be fairly small for the high-resolution simulations described in this paper. The computational costs may, however, vary strongly depending on the number of trajectories and trace variables. Various options have been implemented to initialize online trajectories at different locations and times during the model simulation. As a first application of the new COSMO-model module, an Alpine north foehn event in summer 1987 has been simulated with horizontal resolutions of 2.2, 7 and 14 km. It is shown that lowtropospheric trajectories calculated offline with oneto sixhourly wind fields can significantly deviate from trajectories calculated online. Deviations increase with decreasing model grid spacing and are particularly large in regions of deep convection and strong orographic flow distortion. On average, for this particular case study, horizontal and vertical positions between online and offline trajectories differed by 50–190 km and 150–750 m, respectively, after 24 h. This first application illustrates the potential for Lagrangian studies of mesoscale flows in high-resolution convection-resolving simulations using online trajectories.


Introduction
The Lagrangian depiction of atmospheric processes has a long tradition in atmospheric sciences: the first dynamic studies date back to the beginning of the 20th century, when Shaw et al. (1903) and Shaw and Lempfert (1906) used trajectories to describe the motion of air parcels in cyclones.These first trajectories were so-called surface trajectories, which took into account only the wind fields close to the Earth's surface.Later trajectories were calculated on isobaric surfaces and, when the importance of vertical motion of air parcels became evident, also on isentropic surfaces (Danielsen, 1961).As trajectories calculated on isentropic surfaces can only represent adiabatic motions, three-dimensional trajectories using all three components of the wind field were introduced (e.g., Reap, 1972).These three-dimensional trajectories can either be calculated taking into account subgrid-scale velocities, for instance turbulent motions most relevant in the planetary boundary layer, or only use the resolved-scale wind.Accordingly the models are termed "Lagrangian particle dispersion model" (e.g., FLEXPART, Stohl et al., 2005) and "Lagrangian parcel model" (e.g., LAGRANTO, Wernli and Davies, 1997).Both types of three-dimensional kinematic trajectories have become especially popular since the late 1980s, when good-quality and well-resolved gridded wind fields became available.They are nowadays considered to be the most accurate type of trajectories in the troposphere (Stohl and Seibert, 1998).
From the early studies on, the Lagrangian perspective had a large impact on the advance of the understanding of atmospheric processes, as it for instance allowed the identification of coherent air streams within extratropical cyclones, most prominently the warm conveyor belt (e.g., Whitaker et al., 1988;Wernli and Davies, 1997).Lagrangian studies were also important to connect ozone-rich episodes in the troposphere to stratospheric intrusions (e.g., Buzzi et al., 1984) or to illustrate the flow blocking at mountain ranges during the passage of fronts (e.g., Buzzi and Tibaldi, 1978;Steinacker, 1984;Kljun et al., 2001).Since trajectory calculations are computationally cheap, the Lagrangian approach is also very valuable to assess the climatological frequency and geographical distribution of atmospheric flow features like warm conveyor belts (e.g., Eckhardt et al., 2004) or stratospheretroposphere exchange (e.g., Sprenger and Wernli, 2003).But trajectory analysis has also made its way into other fields of atmospheric sciences: for instance, it has been successfully used for analyzing the detailed microphysical evolution of clouds along the flow (e.g., Haag and Kärcher, 2004;Hoyle et al., 2005;Brabec et al., 2012), identifying evaporative water sources for precipitation (e.g., Bertò et al., 2004;Sodemann et al., 2008), interpreting measurements of stable water isotopes (Pfahl and Wernli, 2008) and in studies on atmospheric chemistry (e.g., Coude-Gaussen et al., 1987;Miller, 1987;Forrer et al., 2000;Methven et al., 2006).Lagrangian particle dispersion models have been successfully used for modeling the dispersion of atmospheric tracers and pollutants (e.g., Gross et al., 1987;Bellasio et al., 2012;Weil et al., 2012).
Besides these large-scale applications of trajectory analysis, Lagrangian parcel models have also been applied to study single events in high-resolution numerical models, e.g., for studying hurricane eye dynamics (Stern and Zhang, 2013), the origin of air parcels feeding convective cells (Wang and Xue, 2012) or nocturnal equatorial oceanic squall lines (Fierro et al., 2009).The Lagrangian approach has also been successfully used for analyzing stratocumulus clouds and entrainment rates in large eddy simulations (LES) (e.g., Stevens et al., 1996;Kogan, 2006;Yamaguchi and Randall, 2012;Yeo and Romps, 2013).
While frequently employed in all fields of meteorology, the accuracy of trajectories derived from measured or modeled wind data is a long-standing point of discussion (e.g., Danielsen, 1961;Kahl, 1993;Stohl, 1998;Brioude et al., 2012).The comparison of computed trajectories to the actual path of an air parcel in the atmosphere is difficult, but several attempts have been made, for instance, with the aid of tracer experiments (e.g., Draxler, 1987;Haagenson et al., 1990;van Dop et al., 1998) and balloon or tetroon flights (e.g., Djuric, 1961;Reisinger and Mueller, 1983;Stohl and Koffi, 1998).
In addition to these experiments different trajectory models have been compared (e.g., Stohl et al., 2001), the sensitivity to the input data frequency has been tested (e.g., Rolph and Draxler, 1990;Doty and Perkey, 1993;Stohl et al., 1995) and the errors associated with the numerical scheme have been investigated (Seibert, 1993).In the case of three-dimensional kinematic trajectories, which are usually based on wind data from a reanalysis data set or a weather forecast, several error sources can be identified (Stohl, 1998): 1. truncation errors -due to neglecting higher order terms in the Taylor expansion of the trajectory equation; 2. interpolation errors -due to the interpolation of the wind field data from the model grid and model output times to the actual trajectory position in space and time 3. wind field errors -due to the nonrepresentativity of the model wind field, prediction errors and errors in the initial conditions.
The third error source depends largely on the quality of the entire forecasting system and the precision of the initial conditions, which are independent of the trajectory model.The truncation error depends on the numerical scheme used for solving the trajectory equation.All of today's frequently used trajectory models -FLEXPART (Stohl et al., 2005), HYS-PLIT (Draxler and Hess, 1998), TRAJKS (Scheele et al., 1996), and LAGRANTO (Wernli and Davies, 1997) -employ the Petterssen scheme (Petterssen, 1940), which is a second-order second-order scheme with a truncation error proportional to t 2 .For further discussion of the numerical properties the reader is referred to Seibert (1993).Seibert (1993) argued that the truncation error should clearly be smaller than the overall trajectory uncertainty, and that the time step should be small enough to fulfill the Courant-Friedrich-Levy criterion, which is a prerequisite for convergent solutions.With a strongly decreasing grid spacing in recent numerical weather prediction models this constrains the time step to some minutes to seconds.All of the abovementioned models compute the trajectories based on winds from weather forecasts or reanalysis data sets.Typically, output from global models (analyses and forecasts) is available every six hours, and from regional models every hour.
As integration time steps of this order would induce unacceptably large errors, temporal interpolation between these output times is required for calculating trajectories.The errors introduced by temporal and spatial interpolation depend strongly on the spatial and temporal resolution of the wind data.While the spatial resolution of numerical weather prediction models has strongly increased over the past years, the output frequency has increased only slowly and therefore often constitutes the limiting factor for a reduction of the interpolation error.
Both the truncation and the temporal interpolation error can be reduced to a minimum if the trajectory equation is solved "online", i.e., during the integration of the Eulerian numerical weather prediction model.In this case the wind fields are available at every time step of the Eulerian model, which is some tens of seconds in state-of-the-art regional weather prediction models.Consequently, with this approach a large increase in data resolution is obtained compared to the standard approach of calculating "offline" trajectories.However, while offline trajectories can be computed forward or backward in time, with the online approach the trajectory equation can only be solved forward in time.
The first implementation of online trajectories we know of has been accomplished by Rössler et al. (1992).More recently the online computation of air parcel trajectories has been utilized in chemistry models within a Lagrangian advection scheme for trace gases (e.g., Becker and Keuler, 2001;Reithmeier and Sausen, 2002).In some studies air parcel trajectories are constructed based on passive tracer fields, which also makes use of the resolved and subgrid-scale wind fields at each time step of the Eulerian model (e.g.Gheusi and Stein, 2002;Duffourg and Ducrocq, 2011).Rössler et al. (1992) showed in a case study that the increased data input frequency indeed significantly alters the pathways of air parcels, particularly over strongly structured topography.In addition to the expected improvement of the trajectory position found by Rössler et al. (1992), the online computation of trajectories may better capture the small-scale variability of the vertical wind field as it is represented on the grid of the Eulerian numerical model.As shown in recent studies by Grell et al. (2004) and Brioude et al. (2012), this is especially relevant for regional weather prediction models with a high spatial resolution.A good representation of the smallscale structure of the wind field is important for investigations of orographic flow and deep convection and for studies using the Lagrangian approach to investigate the evolution of clouds or chemical substances it is essential to capture the vertical wind fluctuations: for instance, the number of homogeneously nucleating ice crystals is very sensitive to the cooling rates along the trajectory (e.g., Kärcher and Lohmann, 2002;Hoyle et al., 2005;Spichtinger and Krämer, 2013).
In this paper we describe a new implementation of the online computation of trajectories based on grid-scale wind velocities in the regional weather prediction model COSMO (Sect.2).To illustrate the capabilities of the online trajectory approach, the module is used in a simulation of an Alpine foehn event in July 1987, which is described in Sect.3. We provide a comparison of the online trajectories to offline trajectories based on COSMO-model output at different output frequencies between one and six hours, with a particular emphasis on the representation of the foehn flow.In addition we investigate the dependence of the detected differences between online and offline trajectories on the horizontal grid spacing of the COSMO-model (14, 7 and 2.2 km).Finally, in Sect.4, potentials of and challenges for the online computation of trajectories are discussed.

The online trajectory module
The goal of our work has been to construct a new module for the numerical weather prediction model COSMO, which allows calculating forward online trajectories and tracing user-specified variables along these trajectories.For the implementation of the online trajectory module, two existing model codes are combined: on the one hand we select the numerical weather prediction model COSMO as the Eulerian model into which the online trajectory calculation should be embedded.The COSMO-model is a limited-area, nonhydrostatic model (Baldauf et al., 2011) that is used for highresolution operational weather forecasting by several mainly European weather services (e.g., the German Meteorological Service (DWD) and MeteoSwiss).On the other hand we base the trajectory calculation procedure (time stepping as well as interpolation) on the trajectory tool LAGRANTO (Wernli and Davies, 1997), which has been employed in numerous studies to compute offline trajectories (e.g., Wernli and Davies, 1997;Stohl et al., 2001;Lefohn et al., 2011;Cirisan et al., 2013;Grams et al., 2013).The LAGRANTO source code has been modified to be consistent with the COSMO-model, the temporal interpolation has been omitted and it has been parallelized to account for the spatial domain decomposition of the COSMO-model.
The major technical criteria for the design of the online trajectory module are (i) to make as few changes to the existing COSMO-model source code as possible, (ii) to write a module that does not rely heavily on existing COSMO-model source code to make its adaption to other numerical weather prediction models straightforward and (iii) to obtain a reasonable computational performance of the COSMO-model with the online trajectory module.In the present version the trajectory module inherits from the COSMO-model only the spatial grid decomposition and partly relies on the IO structure of the COSMO-model for the output of trajectory data.The implementation of the interpolation procedure is based on a staggered Arakawa C grid and terrain-following, rotated spherical coordinates as used in the COSMO-model.Therefore it should be rather easy to port the online trajectory module to other numerical weather prediction models, which use a similar coordinate system and employ a spatial domain decomposition.With the setup described below, only a few additional lines have to be added to the existing code.In the namelist used for starting the COSMO-model, an additional switch is introduced that allows for switching on or off the trajectory calculation, and an additional namelist block allows the user to specify essential parameters of the module like the starting region, the time step for the output and the traced variables.A detailed description of the namelist is provided in the user guide, which is published as a supplement to this article.The complete source code of the described module is available from the authors upon request.The general workflow inside the module and its embedding in the existing model is illustrated by the flow chart in Fig. 1: after initialization of the COSMO-model and the trajectory module (Sect.2.3), the model loops through the time steps of the integration.At the end of each time step the trajectory module is called to calculate the new trajectory positions (Sect.2.1).Then the necessary interprocessor communication takes place (Sect.2.2), and finally, if the time step is an output time step, the trace variables are interpolated to the trajectory positions and are together with the trajectory position written to the output files.The module is written in Fortran 90.For the inter-processor communication the MPICH implementation (www.mpich.org) of the MPI library is used.

Trajectory integration scheme
The physical core has to solve the trajectory equation Dx Dt = u(x, t).For the solution of this equation we employ the Petterssen scheme (Petterssen, 1940), which computes the new trajectory position at time t 1 = t 0 + t from the externally specified velocity u(x, t) using a second-order semi-implicit discretization in space and time: This implicit equation for the new position x(t i+1 ) can be rewritten in fix-point iteration form: For the required starting values, x 0 (t i+1 ) and u(x 0 (t i+1 ), t i+1 ), x(t i ) and u(x(t i ), t i ) are used in the Petterssen scheme.Therefore, the Petterssen scheme can also be viewed as a predictor-corrector method that approximates the velocity for the forward integration by the mean wind between two successive trajectory locations.Expanding the iteration and dropping the explicit time dependency of x yields (3) . . . (5) The new trajectory position is calculated with this scheme at every main model time step of the COSMO-model, which is usually 20 s for simulations with a grid spacing of 2.2 km and 40 s with a grid spacing of 7 km.We choose the Petterssen scheme because it is used by all of the frequently used trajectory models, because it is accurate to the second order and because several studies have shown that higher order schemes do not perform essentially better (e.g., Seibert, 1993).It should be noted that the Petterssen scheme is a socalled "constant acceleration" solution; that is, it neglects the change in acceleration of the air parcel during an integration time step.However, we think that this assumption is justified in the online computation even more than in the offline calculation because of the very small integration time step.
The convergence of the Petterssen scheme depends on the properties of the flow field and the starting values for the iteration.Seibert (1993) showed that fulfilling the Courant-Friedrich-Levy (CFL) criterion is important for obtaining a convergent solution of the trajectory equation.Assuming a maximum wind velocity of 50 m s −1 , the horizontal CFL criterion requires the time step to be smaller than 280 s for a grid spacing of 14 km and below 40 s for a grid spacing of 2 km.The main model time step, at which the trajectory integration is performed, is smaller than these values in the standard COSMO-model setup.In the vertical the CFL criterion is also almost always fulfilled in our test simulations (Fig. 2, left panel).
The number of iterations required for convergence depends on the flow situation and the time step.We analyzed the convergence behavior of the online trajectories in our test simulations (see Sect. 3) by interrupting the iteration either if the trajectory position changes less than a tenth grid spacing in the horizontal and less than 1 m in the vertical between single iterations or if 50 iterations have been performed.The results are summarized in Fig. 2 (right panel) for simulations with three different horizontal resolutions: in almost all cases the solution converges in the first few iteration steps.Only in about 0.1 ‰ of all time steps more than seven iterations are required, but in these instances the solution is probably truly nonconvergent as the iteration does not stop during the first 50 cycles.Based on these results we choose a default value of three iterations.However, it is possible to adapt this number via the namelist.
The trajectory position is written to the output at a userspecified multiple of the major COSMO-model time step together with the variables traced along the trajectories.The values of the trace variables are obtained by interpolation, as described in the next section, from the Eulerian grid to the trajectory position.The trace variables can be defined by the user via the name list.

Interpolation, lower boundary condition and terrain intersection problem
The Petterssen scheme requires the velocity u(x, t) at the parcel location at each integration time step, and therefore an interpolation of the wind data from the Eulerian grid to the parcel position is necessary.We decided to use a threedimensional linear interpolation as it is, for instance, also used in LAGRANTO.The interpolation is performed between the eight neighboring grid points of the trajectory position along the coordinate axes of the COSMO-model; that is, the horizontal interpolation is done along terrainfollowing surfaces.Higher order interpolation is used by some other trajectory models (e.g., FLEXPART, optional in LAGRANTO), but may not always give better results.Moreover the errors introduced by the linear interpolation are much easier to interpret.In contrast to existing offline models, no temporal interpolation of the wind field is required as it is available at each model time step, which is for online trajectories identical to the trajectory integration time step.
Since the COSMO-model uses a staggered Arakawa C grid, the question arises as to how to derive the wind field as well as other parameters close to the surface, i.e. below the lowest model level but above the surface.This question is also strongly linked to the formulation of the lower boundary condition.In the current version we decided to use a linear extrapolation of the horizontal wind velocity components from the two lowest model levels as it is for instance used, as one option, in the fast wave solver of the COSMO-model.The vertical velocity is calculated by terrain-following linear interpolation.
During the calculation of trajectories it can happen that trajectories intersect the topography.One would expect that the number of terrain intersecting trajectories decreases with the time step used for the solution of the trajectory equation.However, we observe the opposite effect for the online trajectories in our case study presented in Sect.3. In the COSMO simulation with a horizontal resolution of 7 km (see Sect. 3) about 5 % more trajectories hit the ground than in the offline trajectory data set based on one-hourly model output.There are three potential reasons for this: one candidate is the formulation of the lower boundary condition, because the above-described linear extrapolation lower boundary condition might be unrealistic very close to the ground.Another potential reason could be the starting values used in the fixpoint iteration of the Petterssen scheme.Finally, the horizontal interpolation tends to smooth the wind field, which can lead to an artificial motion of the trajectory towards the surface, particularly at isolated terrain peaks.In this situation, the high temporal resolution increases the frequency at which such terrain intersecting trajectories are detected.The effect of horizontal interpolation can be illustrated by a simple model assuming a "zig-zag" topography and a terrain For the calculation a zig-zag-shaped topography is assumed (black line).The wind is assumed to be terrain following and linearly increasing with elevation above ground: the horizontal velocity is 3 m s −1 at the first level, which is 10 m above the surface and 6 m s −1 40 m above the surface.The vertical velocity is calculated such that the wind is terrain following on all levels.The horizontal velocity at the surface is obtained by linear vertical extrapolation from the two levels above, and the vertical velocity at the surface is calculated from this extrapolated velocity.If the air parcel is advected below the topography, the wind components at the surface are used.The blue line shows the calculated online trajectory (time step of 1 s) and the cyan line an offline trajectory (time step of 5.3 min).
following wind field, which increases with height (Fig. 3).The velocity at the lowest model level is assumed to be nonzero but terrain following.As is well visible in Fig. 3, the horizontal interpolation smooths the wind field, and therefore the "online" (blue) as well as the "offline" (cyan) trajectory do not follow the terrain at a constant distance.In the example of the "offline" trajectory, this does not lead to a detected terrain intersection as the trajectory points up-and downwind of the peak are above topography.The online trajectory has, due to the shorter time step, much more points and therefore the intersection even with this narrow peak is detected.
To investigate our third point from above in more detail, we have constructed a composite of the terrain, the trajectory elevation and the height of the lowest model level for all trajectories that hit the topography in the case study described in Sect. 3 (Fig. 4).In the last 200 time steps before terrain intersection, the surface elevation as well as the lowest model level are traced along these trajectories.In addition the elevation of the surface, which would be beneath the trajectory if it continued in the same direction and with the same speed as during the last time steps before the terrain intersection, was computed for the next 200 time steps.The composite was constructed by normalizing the data with the maximum elevation along each trajectory segment and averaging.The situation revealed by the composite analysis (Fig. 4) shows a trajectory passing over a narrow peak, as the 400 time steps correspond to an average travel distance of about 3 to 6 grid points.Hence the composite closely resembles the picture constructed above with the idealized "zig-zag" topography, which indicates that horizontal interpolation is one factor contributing to terrain intersections.The strong reduction of the time step from offline to online calculation of trajectories increases the frequency of detection of such events and therefore explains the higher percentage of terrain hitting trajectories.
The composite analysis indicates that the smoothing of the wind field by horizontal interpolation may play a role, but other possible explanations such as the formulation of the lower boundary condition and the choice of the starting values used for the fix-point iteration in the Petterssen scheme cannot be ruled out.Another option for the lower boundary condition would be to use u = v = w = 0 m s −1 at z = 0 m.This no-slip lower boundary condition may reduce the number of trajectories hitting the topography, but it also may decrease the velocities of the trajectories very close to the surface to almost zero.This would stop the trajectories and therefore result in an undesired virtual loss of air parcels similar to terrain intersections.In order to test the effect of an altered formulation of the lower boundary condition, we implemented the no-slip lower boundary condition in the trajectory module.In our case study the number of trajectories that are lost is about the same as with the original extrapolation lower boundary condition.In addition, increasing the number of iterations in the Petterssen scheme and using the surface wind for the iteration, in case the trajectory is advected below the surface in an early iteration step, has little effect on the number of lost trajectories.Further investigations of this issue could consider splitting the time step close to the surface or different choices of the starting value for the fix-point iteration in the Petterssen scheme.As it is rather unsatisfactory to "lose" trajectories during the computation, for the moment being we have adopted the approach of other offline trajectory tools to artificially place every parcel hitting the topography at a point 10 m above the surface.

Parallelization and communication
The COSMO-model employs a spatial grid decomposition for the computation on multiprocessor machines, which constitutes a major difficulty for introducing the trajectory calculation: in contrast to the Eulerian model, for which the grid points have a fixed spatial position -that is, they remain associated with a certain processor during the entire integration -trajectories have no fixed spatial position and hence may pass from the spatial domain associated with a certain processor to another domain.This problem increases the interprocessor communication significantly, which may deteriorate the model performance on multiprocessor machines.In addition the trajectories may not be equally distributed in space at each time instant, which makes it impossible to obtain a workload balance without strongly modifying the existing COSMO-model structure.However, the effects of the workload imbalance may not strongly affect the model performance because usually the number of operations needed on a certain processor to compute the new trajectory positions is much smaller than the number of operations required at the Eulerian grid points.To ensure a perfect workload balance, it would be necessary to either perform the trajectory integration on processors that are not associated with the domain in which the trajectory is located or to reassign Eulerian grid points to different processors at each time step.In both cases parts of the Eulerian variable fields have to be exchanged between processors during each time step.This additional communication is computationally expensive and would offset the effect of a perfect workload balance in almost all applications.We therefore decided to keep a fixed spatial domain decomposition and a corresponding association of trajectories.
The trajectory locations are stored in an array that is a priori known to all processors.A certain processor works only on the entries corresponding to trajectories inside its domain.However, in order to minimize the communication, the trajectory array is only updated at the processor that performs the integration, and in case a trajectory passes to the domain of another processor, the information is transferred to this other processor.As illustrated in Fig. 1 all communications are performed when all processors have finished the forward integration.Therefore a minimum number of communication operations is obtained, and the communication overhead is kept small.Additional communication is required if output has to be written at a certain time step because then the entire trajectory array on the processor responsible for the IO has to be updated.In some cases with very few trajectories passing between processors, a complete update of the trajectory array on all processors after each time step may be faster due to a smaller relative communicational overhead.Therefore this option is also implemented.For all communications the MPI library is used as in the entire COSMO-model.

Selection of trajectory starting points
An essential choice made by the user of the online trajectory module is the specification of the starting points of the trajectories.This is not a trivial task as no backward computation is possible and hence an a priori knowledge of the interesting starting regions and times is required.Starting trajectories at all grid points and time steps is not feasible for highresolution models due to storage limitations and because it will very strongly increase the runtime of the model.In the present version several options to specify the starting region are available: -start trajectories once at each grid point inside a rectangular box (via namelist); -start trajectories once at user-specified coordinates (via external file); -start trajectories repeatedly at fixed locations at userspecified times (via namelist); -start trajectories repeatedly at fixed locations at a regular time interval (via namelist); and -start trajectories at different locations at different times (via external file).
A detailed description is provided in the Supplement.

Results from an Alpine case study
For a first application of the new online trajectory module we simulate the meteorological evolution over central Europe from 25 to 29 July 1987, which has already been investigated by Rössler et al. (1992), Buzzi and Alberoni (1992) and Paccagnella et al. (1992).

Meteorological situation
Between the 25 July and 27 July 1987 an upper-level trough was located over central Europe and the Mediterranean, which propagated slowly eastward, and mostly northwesterly flow prevailed in this region (Fig. 5).On 26 July 1987 the associated surface cold front reached the Alps and was strongly deformed due to the influence of the Alpine orography: along the Rhone Valley a strong mistral was observed, some portions of the cold air spilled over the Alpine ridge and induced north foehn in Ticino and northern Italy, and finally, along the eastern edge of the Alps a low-level jet formed (Buzzi and Alberoni, 1992).As the cold air propagating around the eastern edge of the Alps met the warmer air over the eastern Po Valley, deep convection developed along the convergence line in the afternoon of 26 July (Buzzi and Alberoni, 1992).In addition a moderate Alpine lee cyclone developed over the Adriatic Sea, which was influenced by the retardation and deformation of the cold front by the Alpine ridge (Buzzi and Alberoni, 1992).For the trajectory analysis we focus on the time period during which the cold front passes the Alps, i.e., the afternoon and evening of 26 July 1987.

Modeling framework
For the numerical simulations of the case study, the COSMOmodel version 4.17 (Baldauf et al., 2011) is used with three different spatial resolutions: two with 40 vertical levels and a horizontal grid spacing of 14 km (0.125 To assess the computational performance, each simulation is performed twice, once with and once without the online trajectory module.
At 02:00 UTC on 25 July 1987 in total 24 615 trajectories are started over the British Isles at each grid point between 50 • and 54 • N and −5 • and 2 • E and each model level from the surface up to 5 km.According to the distribution of vertical levels in the Gal-Chen hybrid coordinate system this gives about twice as many trajectories starting below 2 km than above.
A comparison of the sea-level pressure, temperature and precipitation evolution simulated by the COSMOmodel with the analysis by Buzzi and Alberoni (1992) and Paccagnella et al. (1992) reveals a reasonable performance for all simulations (not shown).The developing lee cyclone is slightly shallower than in the observations and the convection over the eastern Po Valley starts about three hours later.Nevertheless, the essential mesoscale phenomena like the foehn flow, a strong mistral and the strong low-level jet around the eastern edge of the Alpine ridge are well captured in the COSMO simulations.
For the evaluation of the trajectory module we also compute offline trajectories with LAGRANTO for each model simulation.The offline trajectories are started at the same points and time as the online trajectories.For the integration of the offline trajectories the wind fields from the COSMO simulations are used at output intervals of one, three and six hours and the integration time step is set to one-twelfth of this time interval.Similar to the online trajectory calculation, air parcels are placed 10 m above the surface if they are advected below ground.

Computational performance of the online trajectory module
To assess the computational performance of the COSMOmodel with the online trajectory module, we use the model simulations performed for the Alpine case study described in the previous section.The number of processors varies with the spatial resolution of the simulation to obtain reasonable runtimes: 16 for the 14 and 7 km simulations and 128 for the 2.2 km simulation.In addition to the position, 10 additional variables are traced along the online trajectories, and all variables are written to the output files every model time step.Twenty-two three-dimensional and 14 two-dimensional Eulerian variables are written to output files every model hour.
The results from this performance test are summarized in Table 1.For this specific setup the runtime increase due to Table 1.Relative runtime increase with respect to the named reference simulation for COSMO-model simulations with the online trajectory module for 24615 trajectories and 10 trace variables.The simulations were run using 16 processors (COSMO14 and COSMO7) or 128 processors (COSMO2.2).
x = 14 km x = 7 km x = 2.2 km without trajectory module (reference: COSMO14 without trajectory module) 0.00 3.16 13.2 with trajectory module (reference: COSMO14 without trajectory module) 0.264 3.45 13.6 with trajectory module (reference: simulation without trajectory module) 0.264 0.0681 0.0366 the trajectory module is below 30 % for the 14 km simulation, and the impact of the trajectory calculation on the runtime decreases to a few percent for the higher-resolution simulations.This is because the number of trajectories remains constant, while the number of grid points is much larger.The observed runtime increase is not due to the integration of the trajectory equation itself, which is computationally cheap to solve, but is mainly caused by additional writing of output, the additional communication between processors and probably in some cases also by the extensive interpolation of trace variables.Of course the expected runtime increase strongly depends on the number of trajectories, the number of traced variables and the number of processors.In general the observed runtime increase is satisfactorily small in this test case.

Comparison of online and offline trajectories
The most prominent mesoscale flow features identified by Buzzi and Alberoni (1992) -the flow splitting at the Alps with a strong low-level wind on either side of the mountain range and the north foehn flow with a particularly strong outflow from the Simplon-Gotthard region -are well captured by the online and offline trajectory calculations (Fig. 6; foehn trajectories are those passing close to the black cross in both panels).Another interesting feature revealed by the trajectory analysis is the strongly ascending branch of air over eastern Europe, which is associated with ascent ahead of the upperlevel trough.Its passage over central Europe is accompanied by trajectories suddenly changing direction from southeast to northeast and rising as, for instance, also observed over the eastern Alps (trajectories rising above 5 km in Fig. 6).A first qualitative impression of the differences between online and offline trajectories can be obtained from Fig. 6: while the flow patterns of the online (bottom panel) and the offline trajectories based on three-hourly output (top panel) agree quite well on first order, significant differences are observable if subsynoptic-scale flow patterns are considered.For instance the flow around Corsica shows a much more detailed flow structure in the online trajectories, the ratio of trajectories passing over and around the Massif Central, the Pyrenees and the Alps varies quite strongly, and over the Po Valley, the westward curvature of the online trajectories that crossed the Alps is much stronger than that of the equivalent offline trajectories.To obtain a more thorough assessment of the path of air parcels calculated from wind fields with different temporal resolutions, trajectories starting at the same location are compared by calculating the average horizontal and vertical transport deviation (AHTD and AVTD).These distance measures are frequently used in the literature to quantify the differences of trajectories in different data sets (e.g., Rolph and Draxler, 1990;Stohl, 1998).The AHTD describes the average over N trajectories of the horizontal distance between each trajectory calculated with two data sets as a function of time: where (x n , y n ) is the position of the nth reference trajectory and (X n , Y n ) the position of the nth test trajectory.The AVTD is calculated with a similar equation for the trajectory heights z n and Z n instead of the horizontal position.Note that AHTD and AVTD only indicate the mean deviation of all pairs at a certain time.The sensitivity of trajectories with the same starting point to the temporal resolution of the wind field data can vary substantially with the flow features they encounter during their path: as illustrated in Fig. 7, online and offline trajectories sometimes take almost the same path, while at other occasions they diverge strongly and spread over entire Europe.It seems that the closer to the starting point a sensitive flow situation is encountered, the larger the final deviation.For instance the strongly diverging trajectory bundle shown in Fig. 7 (solid lines) enters a convective region off the coast of northern France.In this case the representation of small-scale vertical velocity structures strongly influences the final three-dimensional path of the parcel.The other set of trajectories shown in Fig. 7 (dashed-dotted lines) does not encounter such a sensitive situation and remains fairly coherent.AHTD and AVTD were computed with the online trajectories as the reference data set and offline trajectories as the test data set for all spatial and temporal resolutions used in this case study (Fig. 8).In addition online trajectories for simulations with different spatial resolutions are compared to assess the influence of a changing horizontal model resolution (Fig. 8, red lines).In all cases the AHTD increases more or less steadily with increasing simulation time, which can be explained by the increasing divergence of the trajectories once they enter specific flow regions.The AVTD increase is much less steady, which is probably due to the more localized structure of strong vertical winds; the strongest increases in AVTD occur during times when many trajectories pass over steep topography.
Comparing the AHTD and AVTD evolution for the same spatial resolution, but different data input frequencies for the offline trajectories (same line style in cyan, green and blue in Fig. 8) indicates a weaker deviation between offline and online trajectories with increasing temporal resolution of the wind fields used for the offline trajectories: for instance, if COSMO7 results are considered, the AHTD after 24 h (48 h) is 127 km (393 km) for six-hourly offline trajectories, 97 km (329 km) for three-hourly offline trajectories and 61 km (256 km) for one-hourly offline trajectories.For the spatial resolution of the wind fields the opposite behavior is observed: AHTD and AVTD are smallest for the COSMO14 simulation and largest for the COSMO2.2simulation.For example, if one-hourly offline trajectories are used as reference, the AHTD after 24 h (48 h) is 50 km (214 km) for COSMO14-based trajectories, 61 km (256 km) for COSMO7 and 133 km (444 km) for COSMO2.2.This is most likely related to the differences in atmospheric dynamics depending on the spatial resolution: while for the coarsest resolution the flow should be largely hydrostatic, the flow in the COSMO7 simulation has a nonhydrostatic component, and in the COSMO2.2simulation even deep convective motion is explicitly resolved on the Eulerian grid.The offline trajectory method performs worse in finding an accurate numerical solution to the trajectory equation if the flow is less homogeneous in space and time.A comparison of online trajectories based on COSMO simulations with different spatial resolutions confirms that the sensitivity of trajectories to the temporal resolution becomes larger with increasing spatial resolution (red lines in Fig. 8): AHTD and AVTD are smaller if the  (AHTD and AVTD, respectively).Online trajectories are used as reference trajectories and offline trajectories with a data input interval of 1 h (blue), 3 h (green) and 6 h (cyan) as test trajectories.This comparison is done for horizontal resolutions of the COSMO simulation of x = 14 km (dashed-dotted lines), x = 7 km (dashed lines) and x = 2.2 km (solid lines).In addition online trajectories computed for different spatial resolutions are compared, i.e., x = 2.2 km vs. x = 7 km (solid red line), x = 2.2 km vs. x = 14 km (dashed red lines) and x = 7 km vs. x = 14 km (dashed-dotted red lines).The calculated AHTD(t) and AVTD(t) take into account all trajectories that are inside the model domain at time t in the test and reference data set.
COSMO7 online trajectories are used as reference instead of the COSMO2.2online trajectories.Furthermore the explicit representation of deep convective motion in COSMO2.2 has a stronger impact on trajectories than the flow differences between COSMO14 and COSMO7.We conclude this from the fact that the AHTD and AVTD evolution is very similar for COSMO14 and COSMO7 online trajectories as test and COSMO2.2online trajectories as reference.
The total error after four days of forward integration is between 600 and 900 km in the horizontal (extrapolating the COSMO2.2results) and between 700 and 1000 m in the vertical.The AHTD and AVTD values found in the comparison between online and offline trajectories are somewhat larger than those found in other studies trying to estimate the accuracy of trajectories by comparison of different offline trajectory data sets: Rolph and Draxler (1990) found an AHTD of about 400-500 km after 96 h integration for offline trajectories based on six-hourly input data, and Kröner (2011) found an AHTD between 300 and 600 km after 96 h integration for offline trajectories based on three-hourly and six-hourly input data.The discrepancy may have three reasons: first, both cited studies used wind fields with much coarser spatial resolutions (50 to 360 km) than applied here, and it is obvious from our results that the errors are smaller if wind fields with coarser spatial resolution are considered.Secondly, Rolph and Draxler (1990) and Kröner (2011) averaged trajectories from different synoptic conditions and over different regions (North America and the entire Northern Hemisphere, respectively) for their comparison.This is anticipated to reduce the average error, as many meteorological conditions are less complex and variable than the crossing of a cold front over the Alps.Finally, they used offline trajectories based on onehourly wind data as reference data set for their evaluation, which are potentially affected by significant errors.In our case study we found that using offline trajectories based on one-hourly wind data as reference decreases the final AHTD by about 50 to 100 km and the final AVTD by about 50 to 100 m (not shown).

Foehn flow over the Alps
The analysis of the online and offline trajectories indicates that the differences are particularly large for mesoscale flow features.Because trajectories have been used quite frequently to study orographic flows in the Alps (e.g., Kljun et al., 2001;Würsch, 2009;Roch, 2011), we decided to perform a more detailed analysis of the representation of the north foehn flow in the different trajectory data sets.In each data set, all trajectories that reached a minimum elevation below 1500 m over the Po Valley after crossing the Alps were selected as foehn trajectories.
One of the most interesting features of the trajectories is the change in elevation across the Alpine ridge, which has significant implications for the long-standing discussion about foehn mechanisms (e.g., Steinacker, 2006;Drobinski et al., 2007).The elevation change of trajectories across the Alpine ridge is computed by subtracting the minimum elevation of the foehn trajectories over the Po Valley from their minimum elevation over the Swiss Plateau.For most foehn trajectories in all of our trajectory data sets, this elevation change (Fig. 9) is positive, which means that the majority of air parcels contributing to the foehn event descended during the passage of the Alpine ridge.However, the distribution of the elevation change varies on the one hand with the horizontal resolution of the model, but on the other hand also with the data input frequency used for the trajectory calculation (Fig. 9, left panels): for the online trajectories based on the 14 km simulation, the distribution is strongly peaked with a maximum around 1500 m, but for the simulations with higher horizontal resolution, the distribution becomes flatter and the maximum shifts to about 800-1000 m.If different offline and online trajectory data sets are compared, the shape of the distribution changes strongly: offline trajectories based on six-and three-hourly output data show a rather flat distribution of the elevation change, but for offline trajectories based on one-hourly output data and online trajectories, there is a clear peak around 800-1000 m elevation change.As the number of foehn trajectories varies from data set to data set, normalized distributions were also analyzed (not shown).The normalization has no effect on the qualitative differences in the distribution between data sets.
The comparison of the distributions of elevation south of the Alps (Fig. 9, right panels) for different data input frequencies shows a similar pattern to the elevation change across the Alps.However, here the distribution for offline trajectories based on one-hourly output data and for online trajectories differs significantly for elevations below 500 m: while the online trajectories show a structure reminiscent of a low-level jet with a core just below 500 m, no such features is visible in the other data sets.This "low-level jet" is only captured in the COSMO2.2online trajectories; at lower temporal resolution the distribution is flatter or the maximum is shifted to higher elevations.The differences in the distribution for online trajectories from COSMO simulations with different spatial resolutions reflect to a large degree the representation of the low-level jet in the Eulerian model.It becomes clear that for high-resolution simulations online trajectories are very beneficial for capturing and illustrating the physical processes related to foehn flow.

Potential and challenges
A new module for the nonhydrostatic numerical weather prediction model COSMO has been developed, which calculates air mass trajectories using the grid-scale model wind field at every time step during the integration of the Eulerian model.With this method no temporal interpolation of the wind field data is required and the trajectory equation is integrated with a very small time step corresponding to the Eulerian model time step.Such a small time step, as well as the elimination of the temporal interpolation, should make the numerical solution of the trajectory equation more accurate.The new module was tested by simulating an Alpine north foehn event in July 1987, which exhibited a rich mesoscale phenomenology along the Alpine ridge.Although it is not possible in this study to objectively verify trajectories with measurements, we can conclude that the pattern of the online trajectories is physically meaningful, compares well with offline trajectories on the synoptic scale and resolves many important flow phenomena on the mesoscale.The latter are in general not well represented by offline trajectories, particularly if they are based on low-frequency output data.Capturing smaller scale fluctuations in the wind field does not only add additional details to the trajectories but also alters their path significantly over the entire simulated period: for our Alpine foehn case study, after 96 h forward integration, an offline trajectory is on average displaced by about 600-900 km in the horizontal (AHTD) and by about 700-1000 m in the vertical (AVTD) compared to the online trajectory with the same starting point.
Besides the clear advantages of the online approach, there are also several challenges that should be kept in mind: first of all the COSMO-model cannot be used in backward mode, which means that only forward trajectories can be calculated online.This may complicate studies seeking the explanation of a certain point observation with the help of the history of the sampled air parcel, as for instance frequently done in air pollution (e.g., Forrer et al., 2000) or cloud studies (e.g., Haag and Kärcher, 2004).A related difficulty is the choice of the starting points in general.Some a priori knowledge about interesting meteorological phenomena and their spatio-temporal occurrence is needed to define the starting points before starting the model simulation.This problem may be negligible for problems dealing for instance with the dispersion of a pollutant from a fixed point source, but is not always trivial for other studies.
An unexpected challenge tied to the numerical implementation relates to terrain intersecting trajectories: while one would expect a decrease of the number of terrain intersecting trajectories with a decreasing time step, we observed an increase in the online trajectory data set compared to the offline trajectories based on one-hourly model output.A composite analysis suggests that this is caused by the combination of the still-required horizontal interpolation, the decrease of the time step leading to more data points along the trajectory, and narrow topography peaks.Another possible reason is the formulation of the lower boundary condition, which at the moment is obtained by vertical extrapolation from the two lowest model levels.A no-slip lower boundary condition has also been tested, but the number of trajectories that either hit the topography or get stuck close to the surface due to the zero wind velocity at the surface remains almost unaltered.Other yet unexplored potential remedies for this problem are a split of the integration time step close to the surface or different starting values for the iteration in the Petterssen scheme.At the moment the unsatisfactory loss of trajectories due to ground intersection is avoided by artificially placing these trajectories again 10 m above the surface.This solution is convenient, but for the future we hope to find a more physically justified numerical solution to this problem.
An essential property of the described online trajectory module is the neglect of subgrid-scale processes in the solution of the trajectory equation, which impacts potential applications of the module.Lagrangian parcel models as our online trajectory module have different strengths and weaknesses compared to Lagrangian particle dispersion models (LPDM), which explicitly include diffusive processes: the Lagrangian parcel model represents the average properties of an air parcel with a typical volume of a grid cell.The motion of such an air parcel represent the mean of a particle plume starting within a grid box in a Lagrangian particle model.As noted for instance by Stevens et al. (1996) and utilized in trajectory-based moisture source diagnostics (Sodemann et al., 2008), the time-averaged result of mixing is represented on the scale of grid boxes along parcel trajectories.Therefore on temporal scales corresponding to the grid spacing and the advection velocity, the variation in a finite size box is well captured by the parcel model.If subgrid scale variations and according timescales are the focus of the study, then Lagrangian particle dispersion models are the tool of choice.Note that a much larger number of particles must then be calculated (compared to the number of parcels with our approach) in order to statistically sample the subgridscale variations.For instance, as illustrated by Stevens et al. (1996) in a study on timescales in nonprecipitating stratocumulus clouds, a microphysical box model driven with a Lagrangian parcel model may have problems at cloud edges as warming and drying rates of individual parcels may be too strong due to the neglect of subgrid-scale variations in humidity and temperature.Nevertheless parcel models are successfully used in the literature for Lagrangian analyses of LES simulations (Yeo and Romps, 2013).In contrast to LPDM they allow for the study of the influence of nonresolved mixing, be it from parameterizations or numerical diffusion on the mean properties of air parcels.In addition, as Yeo and Romps (2013) pointed out, air parcel trajectories ensure a constant mass of dry air associated with the trajectory, while this is not the case if subgrid-scale velocities are additionally taken into account.
In addition it is important to keep in mind that the represented processes using mean-wind trajectories strongly depend on the grid spacing and hence differ between LES and NWP applications: while it may be inappropriate (or impossible) to study deep convection with a Lagrangian parcel model in a NWP model that does not resolve convective processes, it is justified in convection-resolving models.Online trajectories aim to represent the motion of air parcels as accurately as possible according to the resolved-scale wind field.Thereby the air parcels are not regarded as "closed boxes" but are permeable for subgrid-scale motions, which the model does not aim to represent explicitly.If the latter is the objective of a study, then Lagrangian particle models must be used.The differentiation between subgrid-scale processes and the resolved wind is fundamental, although it relates to different scales and processes for different model resolutions.
Despite the challenges associated with the online computation of trajectories, this novel possibility for performing Lagrangian studies is supposed to be useful for highresolution simulations and even mandatory for studying atmospheric phenomena with short temporal and spatial scales such as orographic flows or deep convection.There the advent of convection-resolving weather and climate predictions in combination with the calculation of online trajectories can lead to novel insight into the evolution of convective weather systems.For instance, as illustrated in Fig. 10, online trajectories can capture the rapid ascent in cumulonimbus clouds from the boundary layer to the upper troposphere and even provide enough data points during the ascent to study the in-cloud processes.Whether such trajectories are realistic of course depends to a large degree on the quality of the underlying Eulerian model, but such trajectories may also help to validate the Eulerian model in a more process-oriented way.Future studies addressing this verification aspect more closely, probably also employing observational data, will be very helpful in assessing the accuracy of the online trajectories.

Fig. 1 .
Fig. 1.Flowchart of the online trajectory module for the COSMOmodel.

Fig. 2 .
Fig. 2. Left: distribution of the vertical Courant numbers during all integration and iteration steps during the simulation of our Alpine case study.Right: number of iterations required until the solution of the trajectory equation converges.The criterion for convergence is that the horizontal position change by less than one-tenth of the horizontal grid spacing and less then 1 m in the vertical.

Fig. 3 .
Fig.3.Illustration of the terrain intersection problem.For the calculation a zig-zag-shaped topography is assumed (black line).The wind is assumed to be terrain following and linearly increasing with elevation above ground: the horizontal velocity is 3 m s −1 at the first level, which is 10 m above the surface and 6 m s −1 40 m above the surface.The vertical velocity is calculated such that the wind is terrain following on all levels.The horizontal velocity at the surface is obtained by linear vertical extrapolation from the two levels above, and the vertical velocity at the surface is calculated from this extrapolated velocity.If the air parcel is advected below the topography, the wind components at the surface are used.The blue line shows the calculated online trajectory (time step of 1 s) and the cyan line an offline trajectory (time step of 5.3 min).

Fig. 4 .
Fig. 4. Composite of the surface elevation (green), trajectory height (blue) and lowest model level (grey) for all terrain intersecting trajectories in the COSMO7 simulation of our Alpine case study.

Fig. 5 .
Fig.5.The potential vorticity distribution on the 320 K isentropic surface is shown in colors (in pvu) and the sea-level pressure field in blue contours (every 2 hPa) at 18:00 UTC, 26 July 1987 based on the COSMO14 simulation.

Fig. 6 .Fig. 6 .
Fig. 6.Offline trajectories based on three-hourly output (top) and online trajectories (bottom) calculated for a COSMO2.2simulation.Only trajectories starting south of 8.5 • N (rotated coordinates) and between 1400 m and 1500 m altitude are shown.The colors denote the height of the trajectories above sea-level (in meters).The trajectories that pass close to the black cross are the foehn trajectories in both panels.

Fig. 7 .
Fig. 7. Example comparisons of online (red) and offline trajectories based on one-(blue), three-(green) and six-hourly (cyan) COSMOmodel output starting at the same position calculated for the simulation with x = 2.2 km.The two examples (solid and dashed-dotted lines) illustrate extreme cases of similar and divergent trajectories calculated from wind data available at different temporal frequencies.

Fig. 8 .
Fig. 8. Average horizontal (top) and vertical (bottom) transport deviation for different trajectory data set pairs(AHTD and AVTD,  respectively).Online trajectories are used as reference trajectories and offline trajectories with a data input interval of 1 h (blue), 3 h (green) and 6 h (cyan) as test trajectories.This comparison is done for horizontal resolutions of the COSMO simulation of x = 14 km (dashed-dotted lines), x = 7 km (dashed lines) and x = 2.2 km (solid lines).In addition online trajectories computed for different spatial resolutions are compared, i.e., x = 2.2 km vs. x = 7 km (solid red line), x = 2.2 km vs. x = 14 km (dashed red lines) and x = 7 km vs. x = 14 km (dashed-dotted red lines).The calculated AHTD(t) and AVTD(t) take into account all trajectories that are inside the model domain at time t in the test and reference data set.

Fig. 9 .
Fig. 9.Histogram of the height difference of the trajectories between the southern and the northern side of the Alps (left panels) and of the height of the trajectories on the southern side of the Alps (right panels).In the top panels the distribution of heights is compared for online trajectories based on simulations with different horizontal resolutions, and in the bottom panels online and offline trajectories based on COSMO2.2 are compared.

Fig. 10 .
Fig. 10.Online trajectory ascending in a deep convective cloud in the COSMO2.2simulation over the southern Po Valley in the evening of the 26 July 1987.The colors denote the position of the trajectory relative to the W-E-oriented vertical section at 45 • (yellow: in the plane; red: further north; magenta: further south).The color shading shows the total hydrometeor content in g kg −1 and the contours indicate vertical velocity (solid: upward motion; dashed: downward motion) at 20:00 UTC, 26 July 1987, which is the model output time step closest to the ascent of the parcel.