1  Introduction
Lagrangian particle models compute trajectories of a large number of so-called particles (not necessarily representing real particles, but infinitesimally small air parcels) to de- scribe the transport and diffusion of tracers in the atmo- sphere.  The main advantage of Lagrangian models is that, unlike in Eulerian models, there is no numerical diffusion. Furthermore, in Eulerian models a tracer released from a point source is instantaneously mixed within a grid box, whereas Lagrangian models are independent of a computa- tional grid and have, in principle, infinitesimally small reso- lution. The basis for current atmospheric particle models was laid by Thomson (1987), who stated the criteria that must be ful- Correspondence to: A. Stohl (ast@nilu.no) filled in order for a model to be theoretically correct.  A monograph on the theory of stochastic Lagrangian models was published by Rodean (1996) and another good review was written by Wilson and Sawford (1996).  The theory of modeling dispersion backward in time with Lagrangian par- ticle models was developed by Flesch et al. (1995) and Seib- ert and Frank (2004). Reviews of the more practical aspects of particle modeling were provided by Zannetti (1992) and Uliasz (1994). This note describes FLEXPART, a Lagrangian particle dis- persion model that simulates the long-range and mesoscale transport, diffusion, dry and wet deposition, and radioactive decay of tracers released from point, line, area or volume sources. It can also be used in a domain-filling mode where the entire atmosphere is represented by particles of equal mass.  FLEXPART can be used forward in time to simulate the dispersion of tracers from their sources, or backward in time to determine potential source contributions for given re- ceptors.  The management of input data was largely taken from FLEXTRA, a kinematic trajectory model (Stohl et al., FLEXPART’s first version was developed during 1995). the first author’s military service at t the nuclear-biological- chemical school of the Austrian Forces, the deposition code was added soon later (version 2), and this version was val- idated using data from three large tracer experiments (Stohl et al., 1998).  Version 3 saw performance optimizations and the development of a density correction (Stohl and Thomson, 1999). Further updates included the addition of a convection scheme (Seibert et al., 2001) (version 4), better backward cal- culation capabilities (Seibert and Frank, 2004), and improve- ments in the input/output handling (version 5).  Validation was done during intercontinental air pollution transport stud- ies (Stohl and Trickl, 1999; Forster et al., 2001; Spichtinger et al., 2001; Stohl et al., 2002, 2003), for which also special developments for FLEXPART were made in order to extend the forecasting capabilities (Stohl et al., 2004).  The most recent version described here is 6.2, which saw corrections © 2005 Author(s). This work is licensed under a Creative Commons License. to the numerics in the convection scheme, the addition of a domain-filling option, and the possibility to use output nests. FLEXPART is coded following the Fortran 77 standard and tested with several compilers (gnu, Absoft, Portland Group) under a number of operating systems (Linux, So- laris, etc.). The code is carefully documented and optimized for run-time performance.  No attempts have been made to parallelize the code because the model is strictly linear and, therefore, it is most effective to partition problems such that they run on single processors and to combine the results if needed. FLEXPART's source code and a manual are freely avail- able from the internet page http://zardoz.nilu.no/~~andreas/ flextra+flexpart.html.  According to a recent user survey, at least 25 groups from 14 countries are currently using FLEX- PART. The version of FLEXPART described here is based on model level data of the numerical weather prediction model of the European Centre for Medium-Range Weather Fore- casts (ECMWF). Other users have developed FLEXPART versions using input data from a suite of different global (e.g., from the National Centers of Environmental Prediction) and meso-scale (e.g., MM5) models, some of which are available from the FLEXPART website but are not described here.
2  Input data and grid definitions
FLEXPART is an off-line model that uses meteorological fields (analyses or forecasts) in Gridded Binary (GRIB) for- mat from the ECMWF numerical weather prediction model (ECMWF, 1995) on a latitude/longitude grid and on na- tive ECMWF model levels as input.  The data can be re- trieved from the ECMWF archives using a pre-processor that is also available from the FLEXPART website but not de- scribed here.  The GRIB decoding software is not provided with FLEXPART but is publicly available (see links on the FLEXPART website). The data can be global or only cover a limited area. Furthermore, higher-resolution domains can be nested into a mother domain. The file includepar contains all relevant FLEXPART parameter settings, both physical constants and maximum field dimensions.  As the memory required by FLEXPART is determined by the various field dimensions, it is recom- mended that they are adjusted to actual needs before compi- lation. The file includecom defines all FLEXPART global variables and fields, i.e., those shared between most subrou- tines.
2.1
Input data organisation
A file pathnames must exist in the directory where FLEX- PART is started. It must contain at least four lines: 1. line: Directory where all the FLEXPART command files are stored. 2. line: Directory to which the model output is written. 3. line: Directory where the GRIB input fields are located. 4. line: Path name of the AVAILABLE file (see below). If nests with higher-resolution input data shall also be used, lines 3 and 4 must be repeated for every nest, thus also spec- ifying the nesting level order. Any number of nesting levels can be used up to a maximum (parameter maxnests). The meteorological input data must be organised such that all data for a domain and a certain date must be contained in a single GRIB file. The AVAILABLE file lists all avail- able dates and the corresponding file names.  For each nest- ing level, the input files must be stored in a different directory and the AVAILABLE file must contain the same dates as for the mother domain. Given a certain particle position, the last (i.e., innermost) nest is checked first whether it contains the particle or not.  If the particle resides in this nest, the mete- orological data from this nest is interpolated linearly to the particle position. If not, the next nes is checked, and so forth until the mother domain is reached.  There is no nesting in the vertical direction and the poles must not be contained in any nest. The maximum dimensions of the meteorological fields are specified by the parameters nxmax, nymax, nuvz- max, nwzmax, nzmax in file includepar, for x, y, and three z dimensions, respectively. The three z dimensions are for the original ECMWF data (nuvzmax, nwzmax for model half levels and model levels, respectively) and trans- formed data (nzmax, see below), respectively. The horizon- tal dimensions of the nests must be smaller than the param- eters nxmaxn, nymaxn. Grid dimensions and other ba- sic things are checked in routine gridcheck.f , and error messages are issued if necessary. The longitude/latitude range of the mother grid is also used as the computational domain.  All internal FLEXPART coordinates run from the western/southern domain bound- ary with coordinates (0,0) to the eastern/northern bound- ary with coordinates (nx-1,ny-1), where (nx,ny) are the mother grid dimensions.  For global input data, FLEX- PART repeats the westernmost grid cells at the eastern- most domain "boundary", in order to facilitate interpola- tion on all locations of the globe (e.g., if input data run from 0 to 359° with 1° resolution, 0° data are repeated at 360°). A global mother domain can be shifted by nxshift (file includepar) data columns (subroutines shift_field.f and shift_field_0.f) if nested in- put fields would otherwise overlap the "'boundaries"'. For in- stance, a domain stretching from 320° to 30° can be nested into the mother grid of the above example by shifting the mother grid by 30°.
2.2
Vertical model structure and required data
FLEXPART needs five three-dimensional fields: horizon- tal and vertical wind components, temperature and specific humidity.  Input data must be on ECMWF model (i.e., η) levels which are defined by a hybrid coordinate system. The conversion from η to pressure coordinates is given by p k = A k + B k p s  and the heights of the η surfaces are defined by η k = A k / P 0 + B k , where η k is the value of η at the k model level, p s is the surface pressure and p 0 =101325 Pa. A k and B k are coefficients, chosen such that the levels closest to the ground follow the topography, while the highest levels co- incide with pressure surfaces; intermediate levels transition between the two.  The vertical wind in hybrid coordinates is calculated mass-consistently from spectral data by the pre- processor. A surface level is defined in addition to the regular η levels. 2 m temperature, 10 m winds and specific humidity from the first regular model level are assigned to this level, to represent “surface” values. Parameterized random velocities in the atmospheric boundary layer (ABL, see Sects. 3 and 4) are calculated in units of m s−1, and not in η coordinates. Therefore, in order to avoid time-consuming coordinate transformations every time step, all three-dimensional data are interpolated linearly from the ECMWF model levels to terrain-following Carte- sian coordinates  ξ=z−z r , where z r is the height of the to- pography (subroutine verttransform.f). The conver- sion of vertical wind speeds from the eta coordinate system into the terrain-following co-ordinate system follows as
where ň=ἡ∂p/∂η. The second term on the right hand side is missing in the FLEXPART transformation because it is much smaller than the others. One colleague has implemented this term in his version of FLEXPART and found virtually no differences (B. Legras, personal communication). FLEXPART also needs the two-dimensional fields:  sur- face pressure, total cloud cover, 10 m horizontal wind components, 2 m temperature and dew point temperature, large scale and convective precipitation, sensible heat flux, east/west and north/south surface stress, topography, land- sea-mask and subgrid standard deviation of topography. A landuse inventory (for Europe, the data set of Velde et al. (1994) is used) must be provided in an extra file (landuse.asc).
Physical parameterization of boundary layer param-
3
eters
Accumulated surface sensible heat fluxes and surface stresses are available from ECMWF forecasts. The pre-processor se- lects the shortest forecasts available for that date from the ECMWF archives and deaccumulates the flux data. The total surface stress is computed from where τ 1 and τ 2 are the surface stresses in east/west and north/south direction, respectively.  Friction velocity is then calculated in subroutine scalev.f as
where ρ is the air density (Wotawa et al., 1996).  Friction velocities and heat fluxes calculated using this method are most accurate (Wotawa and Stohl, 1997).  However, if deac- cumulated surface stresses and surface sensible heat fluxes are not available, the profile method after Berkowicz and Prahm (1982) (subroutine pbl_profile.f) is applied to wind and temperature data at the second model level and at 10 m (for wind) and 2 m (for temperature) (note that previ- ously the first model level was used; as ECMWF has its first model level now close to 10 m, the second level is used in- stead). The following three equations are solved iteratively:
where κ is the von Kármán constant (0.4), z l is the height of the second model level, ∆ u is the difference between wind speed at t the second model level and at 10 m, ΔΘ is the dif- ference between potential temperature at the second model level and at 2 m, Ψ m and Ψ h are the stability correction func- tions for momentum and heat (Businger et al., 1971; Beljaars and Holtslag, 1991), g is the acceleration due to gravity, Θ* is the temperature scale and T is the average surface layer temperature (taken as T at the first model level).  The heat flux is then computed by
where ρc p  is the specific heat capacity of air at constant pres- sure. ABL heights are calculated according to Vogelezang and Holtslag (1996) using the critical Richardson number con- cept (subroutine richardson.f). The ABL height h mix is set to the height of the first model level l for which the Richardson number
exceeds the critical value of 0.25. Θ v 1 and Θ v 1 are the vir- tual potential temperatures, z 1 and z 1 are the heights of, and ( u 1 , v 1 ), and ( u 1 , v 1 ) are the wind components at the 1 st and l th model level, respectively. The formulation of Eq. 8 can be improved for convective situations by replacing Θ v 1 with
www.atmos-chem-phys.org/acp/5/2461/ Atmos. Chem. Phys., 5, 2461–2474, 2005 where is the convective velocity scale. The second term on the right hand side of Eq. 9 represents a temperature excess of rising thermals.  As w ∗ is unknown beforehand, h mix and w ∗ are calculated iteratively. Spatial and temporal variations of ABL heights on scales not resolved by the ECMWF model play an important role in determining the thickness of the layer over which tracer is effectively mixed. The height of the convective ABL reaches its maximum value (say 1500 m) in the afternoon (say, at 1700 local time (LT)), before a much shallower sta- ble ABL forms.  Now, if meteorological data are available only at 12:00 and 18:00 LT and the ABL heights at those times are, say, 1200 m and 200 m, and linear interpolation is used, the ABL height at 17:00 LT is significantly underes- timated (370 m instead of 1500 m).  If tracer is released at the surface shortly before the breakdown of the convective ABL, this would lead to a serious overestimation of the sur- face concentrations (a factor of four in the above example). Similar arguments hold for spatial variations of ABL heights due to complex topography and variability in landuse or soil wetness (Hubbe et al., 1997). The thickness of a tracer cloud traveling over such a patchy surface would be determined by the maximum rather than by the average ABL height. In FLEXPART a somewhat arbitrary parameterization is used to avoid a significant bias in the tracer cloud thickness and the surface tracer concentrations.  To account for spa- tial variations induced by topography, we use an "envelope" ABL height
Here, σ z  is the standard deviation of the ECMWF model sub- grid topography, c is a constant (here: 2.0), V is the wind speed at height h mix , and N is the Brunt-Vaisala frequency. Under convective conditions, the envelope ABL height is, thus, the diagnosed ABL height plus the subgrid topogra- phy (assuming that the ABL height over the hill tops effec- tively determines the dilution of a tracer cloud located in a convective ABL). Under stable conditions, air tends to flow around topographic obstacles rather than above it, but some lifting is possible due to the available kinetic energy.  ♕ is the local Froude number (i.e., the ratio of inertial to buoyant forces) times the length scale of the sub-grid topographic ob- stacle.  The factor c H , thus, limits the effect of the subgrid topography under stable conditions, with c=2 being a sub- jective scaling factor.  Henv rather than h mix is used for all subsequent calculations. In addition, Henv is not interpolated to the particle position, but the maximum Henv of the grid points surrounding a particle's position in space and time is used.
Particle transport and diffusion
4
Particle trajectory calculations 4.1
FLEXPART generally uses the simple "'zero acceleration'' scheme
which is accurate to the first order, to integrate the trajectory equation (Stohl, 1998)
with t being time, Δt t the time increment, X the position vec- tor, and v=v +v t +v m  the wind vector that is composed of the grid scale windข, the turbulent wind fluctuations v t and the mesoscale wind fluctuations vm. Since FLEXPART version 5.0, numerical accuracy has been improved by making one iteration of the Petterssen (1940) scheme (which is accurate to the second order) when- ever this is possible, but only for the grid-scale winds. It is implemented as a correction applied to the position obtained with the "'zero acceleration'' scheme.  In three cases it can- not be applied.  First, the Petterssen scheme needs winds at a second time which may be outside the time interval of the two wind fields kept in memory. Second, if a particle crosses the boundaries of nested domains, and third in the ABL if ctl>0 (see below). Particle transport and turbulent dispersion are handled by the subroutine advance.f where calls are issued to pro- cedures that interpolate winds and other data to the par- ticle position and the Langevin equations (see below) are solved.  The poles are singularities on a latitude/longitude Thus, horizontal winds (variables uu,vv) pole- grid. ward of latitudes ( switchnorth, switchsouth ) are transformed to a polar stereographic projection (variables uupol,vvpol) on which particle advection is calculated. As uupol , vvpol are also stored on the latitude/longitude grid, no additional interpolation is made.
The Langevin equation
4.2
Turbulent motions v t for wind components i are parame- terized assuming a Markov process based on the Langevin equation (Thomson, 1987)
where the drift term a and the diffusion term b are functions of the position, the turbulent velocity and time. dW j are in- cremental components of a Wiener process with mean zero and variance dt , which are uncorrelated in time (Legg and Raupach, 1982).  Cross-correlations between the different wind components are also not taken into account, since they have little effect for long-range dispersion (Uliasz, 1994). Atmos. Chem. Phys., 5, 2461–2474, 2005 Gaussian turbulence is assumed in FLEXPART, which is strictly valid only for stable and neutral conditions.  Under convective conditions, when turbulence is skewed and larger areas are occupied by downdrafts than by updrafts, this as- sumption is violated, but for transport distances where parti- cles are rather well mixed throughout the ABL, the error is minor. With the above assumptions, the Langevin equation for the
where w and σ w are the turbulent vertical wind component and its standard deviation, τ L w  is the Lagrangian timescale for the vertical velocity autocorrelation and ρ is density. The second and the third term on the right hand side are the drift correction (McNider et al., 1988) and the density correc- tion (Stohl and Thomson, 1999), respectively. This Langevin equation is identical to the one described by Legg and Rau- pach (1982), except for the term from Stohl and Thomson (1999) which accounts for the decrease of air density with height. Alternatively, the Langevin equation can be re-expressed
This form was shown by Thomson (1987) to fulfill the well- mixed criterion which states that "if a species of passive marked particles is initially mixed uniformly in position and velocity space in a turbulent flow, it will stay that way" (Rodean, 1996). Although the method proposed by Legg and Raupach (1982) violates this criterion in strongly inhomoge- neous turbulence, their formulation was found to be practical, as numerical experiments have shown that it is more robust against an increase in the integration time step.  Therefore, Eq. (15) is used with long time steps (see Sect. 4.3); other- wise, Eq. (16) is used. For the horizontal wind components, the Langevin equation is identical to Eq. (15), with no drift and density correction terms. For the discrete time step implementation of the above Langevin equations (at the example of Eq. 16), two differ- ent methods are used. When ( ∆ t / τ L w ) ≥ 0.5,
where r w = exp( − Δt / τ L w ) is the autocorrelation of the ver- tical wind, and ζ is a normally distributed random number with mean zero and unit standard deviation. The subscripts k and k+1 refer to subsequent times separated by ∆ t . To save computation time for cases when (∆t /τ L … ) <0.5, the following first order approximation is used in order to avoid the computation of the exponential function:
When a particle reaches the surface or the top of the ABL, it is reflected and the sign of the turbulent velocity is changed (Wilson and Flesch, 1993). 4.3
Determination of the time step
FLEXPART can be used in two different modes.  The com- putationally faster one (ctl<0 in file COMMAND) does not adapt the computation time step to the Lagrangian timescales TL; (where i  is one of the three wind components) and FLEX- PART uses constant time steps of one synchronisation time interval (lsynctime, specified in file COMMAND, typically 900 seconds).  Usually, autocorrelations are very low in this mode and turbulence is not described well. Nevertheless, for large scale applications FLEXPART works very well with this option (Stohl et al., 1998).  If turbulence shall be de- scribed more accurately, the time steps must be limited by TL . Since the vertical wind is most important, only TL , , is used for this. The user must specify two constants, ctl and ifine in file COMMAND. The first one determines the time step Δt  according to
The minimum value of Δ t is 1 second. Δ t is used for solv- ing the Langevin equations for the horizontal turbulent wind components. For solving the Langevin equation for the vertical wind component, a shorter time step Δt w =Δt i /ifine is used. However, note that there is no interaction between horizontal and vertical wind components on timescales less than ∆ t i . This strategy (given sufficiently large values for ctl and ifine) ensures that the particles stay vertically well-mixed also in very inhomogeneous turbulence, while keeping the computational cost at a minimum.
4.4
Parameterization of the wind fluctuations
For σ vi and τ Lj Hanna (1982) proposed a parameteri- zation scheme based on the boundary layer parameters h, L, w*, z0 and u*, i.e. ABL height, Monin-Obukhov length, convective velocity scale, roughness length and friction velocity, respectively. It is used in subroutines hanna.f, hannal.f, hanna_short.f with a mod- ification taken from Ryall et al. (1997) for σ w , as Hanna's scheme does not always yield smooth profiles of σ w through- out the whole convective ABL. In the following, sub- scripts u and v refer to the along-wind and the cross-wind components (transformed to grid coordinates in subroutine windalign.f ), respectively, and w to the vertical compo- nent of the turbulent velocities; f is the Coriolis parameter. The minimum τ L α , τ L ν and τ L w used are 10 s, 10 s and 30 s, respectively, in order to avoid excessive computation times for particles close to the surface.
Unstable conditions:

Stable conditions: Lacking suitable turbulence parameterizations above the ABL ( z > h ), a constant vertical diffusivity D z =0.1 m 2 s − 1 is used in the stratosphere, following recent work of Legras et al. (2003), whereas a horizontal diffusivity D h =50 m² s−1 is used in the free troposphere.  Stratosphere and tropo- sphere are distinguished based on a threshold of 2 pvu (po- tential vorticity units). Diffusivities are converted into veloc- ity scales using σν i =√Di√dt .
4.5
Mesoscale velocity fluctuations
Mesoscale motions are neither resolved by the ECMWF data nor covered by the turbulence parameterization. This unre- solved spectral interval needs to be taken into account at least in an approximate way, since mesoscale motions can signif- icantly accelerate the growth of a dispersing plume (Gupta et al., 1997). For this, we use a similar method as Maryon (1998), namely to solve an independent Langevin equation for the mesoscale wind velocity fluctuations ("meandering" in Maryon's terms). Assuming that the variance of the wind at the grid scale provides some information on its subgrid variance, the wind velocity standard deviation used for the mesoscale Langevin equation is set to turbmesoscale (set in file includepar ) times the standard deviation of the grid points surrounding the particle's position. The cor- responding time scale is taken as half the interval at which wind fields are available, assuming that the linear interpo- lation between the grid points can recover half the subgrid variability, not an unlikely assumption (Stohl et al., 1995). This empirical approach does not describe actual mesoscale phenomena, but it is similar to the ensemble methods used to assess trajectory accuracy (Kahl , 1996; Baumann and Stohl, 1997; Stohl, 1998).
Moist convection
4.6

(27)
(33)
An important transport mechanism are the updrafts in con- vective clouds. They occur in conjunction with downdrafts within the clouds and compensating subsidence in the cloud- free surroundings. These convective transports are grid-scale in the vertical, but sub-grid scale in the horizontal, and are not represented by the ECMWF vertical velocity. To represent convective transport in a particle dispersion model, it is necessary to redistribute particles in the en- tire vertical column.  For FLEXPART we chose the con- vective parameterization scheme by Emanuel and Živković- Rothman (1999), as it relies on the grid-scale temperature and humidity fields and calculates a displacement matrix pro- viding the necessary mass flux information for the particle redistribution. The convective parameterization is switched on using lconvection in file COMMAND . It's computation time scales to the square of the number of vertical model levels and may account for up to 70% of FLEXPART's com- putation time using current 60-level ECMWF data. Atmos. Chem. Phys., 5, 2461–2474, 2005. The convection is computed within the subroutines con- vmix.f, calcmatrix.f, convect43c.f, and re- dist.f. It is called every FLEXPART lsynctime time step (typically 900 s) with time-interpolated temperature and specific humiditiy profiles from the ECMWF data. Note that the original ECMWF model levels, not the Cartesian coordi- nates, are used in the convection scheme. For efficiency rea- sons, particles are sorted according to their horizontal grid positions ( sort2.f ) before calling the convection scheme once per grid column. In the Emanuel scheme ( convect43c.f ), convection is triggered whenever
with TLCL+1 the virtual temperature of a surface air parcel lifted to the level above the lifting condensation level LCL , TLCL+1 the virtual temperature of the environment there, and Tthres =0.9 K a threshold temperature value. Based on the buoyancy sorting principle (Emanuel, 1991; Telford, 1975), a matrix M A of the saturated upward and downward mass fluxes within clouds is calculated by accounting for en-
Here M A i , j  are the mass fractions displaced from level i to level j , M i the mass fraction displaced from the surface to level i , LNB the level of neutral buoyancy of a surface air parcel and 0<σ i,j <1 the mixing fraction between level i  and level j . The fractionلاءُis determined by the environmen- tal potential temperature θ ্ , the liquid potential temperature θ f ; } of air displaced adiabatically from i to j , and the liquid potential temperature θ i , j  of an air parcel first lifted adiabat- ically to level i  and further to level  j :
By summing up over all levels j , we then calculate the sat- urated up- and downdrafts at each level i from Eq. (35) and assume that these fluxes are balanced by a subsidence mass flux in the environment. The particles in each convectively active box are then re- distributed ( redist.f ) according to the matrix MA . If the mass of an ECMWF model layer i is m 1 and the mass flux from layer i to layer j accumulated over one time step is △M A 1.1 , then the probability of a particle to be moved from layer i to layer j is ΔΜΑ i j /m i . Whether a given parti- cle is displaced or not is determined by drawing a random number between [0,1], which also determines the position of the particle within the destination layer j . After the con- vective redistribution of the particles, the compensating sub- sidence mass fluxes are converted to a vertical velocity act- ing on those particles in the grid box that are not displaced by convective drafts. By calculating a subsidence velocity rather than displacing particles randomly between layers the scheme's numerical diffusion in the cloud-free environment is eliminated. The scheme was tested and fulfills the well- mixed criterion, i.e., if a tracer is well mixed in the whole atmospheric column, it remains so after the convection.
Particle splitting
4.7
During the initial phase of dispersion from a point source in the atmosphere, particles normally form a compact cloud. Relatively few particles suffice to simulate this initial phase correctly.  After some time, however, the particle cloud gets distorted and particles spread over a much larger area. More particles are now needed.  FLEXPART allows the user to specify a time constant Δ t s (file COMMAND). Particles are split into two (each of which receives half of the mass of the original particle) after travel times of Δt s , 2Δt s , 4Δt s , 8Δt s , and so on (subroutine timemanager.f ).
Forward and backward modeling
5
Normally, when FLEXPART is run forward in time (ldirect=1 in file COMMAND), particles are released from one or a number of sources and concentrations are deter- mined downwind on a grid. However, FLEXPART can also run backward in time ( ldirect=-1 ), which is more effi- cient than forward modeling for calculating source-receptor relationships if the number of receptors is smaller than the number of (potential) sources.  In the backward mode, par- ticles are released from a receptor location (e.g., a measure- ment site) and a four-dimensional (3 space dimensions plus time) response function (sensitivity) to emission input is cal- culated. Since this version (6.2) of FLEXPART, the calcula- tion of the source-receptor relationships is generalized for both forward and backward runs, allowing much greater flexibility regarding the input and output units than be- fore. ind source and ind receptor in file COMMAND switch between mass and mass mixing ratio units at the source and at the receptor, respectively. Note that source al- ways stands for the physical source and not the location of the particle release, which is done at the source in forward mode but at the receptor in backward mode. Table 1 gives an overview of the units used in forward and backward model- ing for different settings of the above switches.  A "normal" forward simulation which specifies the release in mass units (i.e., kg) and also samples the output in mass units (i.e., a concentration in ng m−3) requires both switches to be set to 1.
In the backward mode, any value not equal zero can be entered as the release “mass” in file RELEASES because the output is normalized by this value. The calculated response function is related to the particles' residence time in the out- put grid cells.  The unit of the output response function varies, depending on how the switches are set. If ind_source=1 and ind_receptor=1 , the response function has the unit s. If this response function is folded (i.e., multiplied) with a 3-d field of emission mass fluxes into the output grid boxes (in kg m− 3 s− 1 ), a concentration at the receptor (kg m− 3 ) is obtained. If ind_source=1 and ind_receptor=2, the response function has the unit s m 3 kg− 1 and if it is folded with the emission mass flux (again in kg m− 3 s− 1 ), a mass mixing ratio at the receptor is obtained. The units of the re- sponse function for ind_source=2 can be understood in analogy. In the case of loss processes (dry or wet deposition, de- cay) the response function is "corrected" for these loss pro- See Seibert (2001) and, particularly, Seibert and cesses. Frank (2004) for a description of these generalized in- and output options and the implementation of backward model- ing in FLEXPART. Seibert and Frank (2004) also describe the theory of backward modeling and give some examples, and Stohl et al. (2003) presents an application. There are differences between the dimensions of the out- put fields needed to run FLEXPART forward and back- ward in time.  For backward simulations, the output for all the receptors (particle release locations) must be kept sepa- rate, contrary to the forward runs where releases from sev- eral source locations can be combined.  Therefore, for the backward runs the output fields must contain a dimension maxpoint , the maximum number of release points, that is not needed in forward runs. In order to avoid creating a fur- ther dimension for the output fields (thus increasing FLEX- PART's memory demands), maxspec , the maximum num- ber of chemical species used in forward runs, is replaced with maxpoint in the backward runs as a dimension of the output fields.  This has the disadvantage that only one species (with certain properties regarding the removal pro- cesses, see below) can be calculated in a backward run.  To switch between forward and backward runs, the parameter maxpointspec is used. It must be set (in includecom) to maxspec for forward runs and to maxpoint for back- ward runs.  FLEXPART must be recompiled upon changing this.
Plume trajectories
6
In a recent paper, Stohl et al. (2002) proposed a method to condense the complex and large FLEXPART output using a cluster analysis (Dorling et al., 1992).  The idea behind this is to cluster, at every output time, the positions of all par- ticles originating from a release point, and write out only clustered particle positions, along with additional informa- tion (e.g., fraction of particles in the ABL and in the strato- sphere).  This creates information that is almost as compact as traditional trajectories but accounts for turbulence and con- vection. This option can be activated by setting iout to 4 or 5 in file COMMAND. The number of clusters can be set with the parameter ncluster in file includepar. The clustering is handled and output is produced by subroutine plumetraj.f.
Removal processes
7
FLEXPART takes into account radioactive (or other) decay, wet deposition, and dry deposition by reducing a particle's mass. However, as atmospheric transport is the same for all chemical species, a single particle can represent several (up to maxspec) chemical species, each affected differently by the removal processes.
7.1
Radioactive decay
Radioactive decay is accounted for by reducing the particle mass according to
where m is particle mass, and the time constant β=T 1/2/ ln(2) is determined from the half life T 1/2 speci- fied in file SPECIES. Deposited pollutant mass decays at the same rate.
7.2
Wet deposition
Wet deposition removes aerosols and gases from the atmo- sphere.  In principle, in-cloud and below-cloud scavenging must be separated (Asman, 1995). However, as data on cloud base height and depth are not available, in-cloud and below- cloud scavenging are treated jointly in FLEXPART. Using mm/hour.
scavenging coefficients, wet deposition takes the form of an exponential decay process (McMahon, 1979)
where m and Λ are the particle mass and the scavenging coef- ficient, respectively. The scavenging coefficient A increases with precipitation rate according to
where I is the precipitation rate in mm/hour, A [s- 1 ] is the scavenging coefficient at I =1 mm/hour and B gives the de- pendency on precipitation rate. Both A and B must be speci- fied in file SPECIES . FLEXPART uses the same scavenging coefficients for snow and rain. As wet deposition depends nonlinearly on precipitation rate, subgrid variability of precipitation must be accounted for (Hertel et al., 1995). The area fraction which experiences precipitation given a certain grid-scale precipitation rate is calculated by
where CC is the total cloud cover, I t and I c are the large scale and convective precipitation rates, respectively, and f r 1 and f r c are correction factors that depend on I 1 and I c (see Table 2).  The subgrid scale precipitation rate is then I s =( I l + I c ) / F .
7.3
Dry deposition
Dry deposition is described in FLEXPART by a deposition velocity
where F C and C are the flux and the concentration of a species at height z within the constant flux layer. A constant deposition velocity vd can be set (file SPECIES ). Alterna- tively, if the physical and chemical properties of a substance are known (file SPECIES ), more complex parameterizations for gases and particles are also available. used to calculate the roughness length.

Dry deposition of gases
7.3.1
The deposition velocity of a gas is calculated with the re- sistance method (Wesely and Hicks, 1977) in subroutine getvdep.f according to
where r a is the aerodynamic resistance between z and the surface, r b is the quasilaminar sublayer resistance, and r c is the bulk surface resistance. The aerodynamic resistance r a is calculated in func- tion raerod.f using the flux-profile relationship based on Monin-Obukhov similarity theory (Stull, 1988)
Following Erisman et al. (1994), the quasilaminar sublayer resistance is
where Sc and Pr are the Schmidt and Prandtl numbers, re- spectively. Pr is 0.72 and Sc=υ / D i , with υ being the kine- matic viscosity of air and D i being the molecular diffusivity of species i  in air. The slight dependency of υ on air temper- ature is formulated in accordance with Pruppacher and Klett (1978). rb is calculated in function getrb.f. The surface resistance is calculated in function getrc.f following Wesely (1989) as
where r s , r m and r 1 μ represent the bulk values for leaf stom- atal, leaf mesophyll and leaf cuticle surface resistances (allto- gether the upper canopy resistance) , r dc represents the gas- phase transfer affected by buoyant convection in canopies, rcl the resistance of leaves, twig, bark and other exposed surfaces in the lower canopy, rac the resistance for trans- fer that depends only on canopy height and density, and r gs the resistance for the soil, leaf litter, etc., at the ground.
Each of these resistances is parameterized according to the species' chemical reactivity and solubility, the landuse type, and the meteorological conditions.  The landuse inventory (Velde et al., 1994) provides the area fractions of eight lan- duse classes for which roughness lengths zo are estimated, on a grid with 10' resolution (Table 3). Charnock's relationship (Stull, 1988) zo=0.016u 2 /g is used to calculate zo for the classes "Ocean" and "Inland water", because of its depen- dence on wave height.  Deposition velocities are calculated for all landuse classes and weighted with their respective ar- eas.  Outside of Europe, the landuse classes are determined only from the ECMWF land-sea-mask, attributing the lan- duse classes "Ocean" to the sea surfaces and "Grasslands" to the land surfaces. 7.3.2
Dry deposition of particulate matter
The deposition of particulates is calculated in subroutine partdep.f according to
where ν , is the gravitational settling velocity calculated from (Slinn, 1982)
where ρ p and d p are the particle density and diameter, μ the dynamic viscosity of air (0.000018 kg m − 1 s − 1 ) and Ccun the Cunningham slip-flow correction. The quasilaminar sub- layer resistance is calculated from the same relationship as for gases, with an additional impaction term. For further de- tails see Slinn (1982). Settling and dry deposition velocities are strongly depen- dent on particulate size. FLEXPART assumes a logarithmic normal size distribution of the particulate mass.  The user must specify the mean particulate diameter d p and a mea- sure of the variation around , σ ,. Then, the settling and deposition velocities are calculated for several particle diam- eters and are weighted with their respective particulate mass fractions. Gravitational settling is important not only for the com- putation of the dry deposition velocity, but also affects the particle's trajectory. As a FLEXPART particle can normally represent several species, gravitational settling can only be taken into account correctly (i.e., influence particle trajecto- ries) in single-species simulations. 7.3.3
Loss of particle mass due to dry deposition
The depositon velocity is calculated for a reference height (parameter href in file includepar) of 15 m. For all par- ticles below 2href , the mass lost by deposition is calculated by
Atmos. Chem. Phys., 5, 2461–2474, 2005 8
Calculation of concentrations, uncertainties, age
spectra, and mass fluxes
Output quantities C T at time T c (output interval loutstep is set in file COMMAND) are calculated as time-averages over period [Tc-ATc/2,Tc+ATc/2]. △T c must be specified ( loutaver ) in file COMMAND. To calculate the time-averages, concentrations C T s at times T s within [Tc-ATc/2,TC+ATc/2] are sampled at shorter intervals △Ts (loutsample in file COMMAND) and are then divided by the number N= ATs of samples taken:
Both △T c and △T s must be multiples of the FLEXPART synchronisation interval ( lsynctime in file COMMAND). The shorter the sampling interval △ T s , the more samples are taken and the more accurate are thus the time-averaged con- centrations. Concentrations, mixing ratios, and emission response 8.1 functions The user can choose ( iout in file COMMAND , which must be set to 1 for backward runs) whether concentrations, vol- ume mixing ratios or both shall be produced. We shall use the term "concentration" and particle mass here, but note that the actual units are determined by the settings of ind_source and ind_receptor , according to Table 1. The concentra- tion in a grid cell is calculated in subroutine conccalc.f by sampling the tracer mass fractions of all particles within the grid cell and dividing by the grid cell volume
with V being the grid cell volume, m i particle mass, N the total number of particles, and f i the fraction of the mass of particle i attributed to the respective grid cell. This mass fraction is calculated by a uniform kernel with bandwidths ( Ax, Ay), where Ax and Ay are the grid distances on the longitude-latitude output grid. Figure 1 illustrates this: The particle is located at the center of the shaded rectangle with side lengths ( ∆ x, ∆ y). Generally, the shaded area stretches over four grid cells, each of which receives a fraction of the particle's mass equal to the fraction of the shaded area falling within this cell. The uniform kernel is not used during the first 3 hours after a particle's release (when the mass is at- tributed only to the grid cell it resides in), in order to avoid smoothing close to the source. Wet and dry deposition fields are calculated on the same output grid (subroutines wetdepokernel.f and drydepokernel.f) and are written to all output grid files. The deposited matter is accumulated over the course of a model run, i.e. it generally increases with model time. However, radioactive decay is calculated also for the de- posited matter.
8.2
Uncertainties
The uncertainty of the output is estimated by carrying nclassunc classes of particles in the model simulation, and determining the concentration separately for each class (subroutine conccalc.f ). The standard deviation, cal- culated from nclassunc concentration estimates and di- vided by ✓nclassunc, is the standard deviation of the mean concentration (subroutine concoutput.f ), which is also written to the output files for every grid cell. Note that the memory needed for some auxiliary fields increases with nclassunc t and the number of age classes (see below). It may, thus, be necessary to reduce nclassunc for runs with large output grids and age spectra calculations or in the backward mode.
8.3
Age spectra
The age spectra option is switched on using lagespectra in file COMMAND, with the age classes specified in seconds in file AGECLASSES. Concentrations are split into contri- butions from particles of different age, defined as the time passed since their release. Particles are terminated once they are older than the oldest age class and their storage space is made available to new particles.  Therefore, the age spectra option can be used also with a single age class for defining a maximum particle age.
8.4
Parabolic kernel
In addition to the simple uniform kernel method, a computa- tionally demanding parabolic kernel as described in (Uliasz, 1994) can be used to calculate surface concentrations for a limited number of receptor points (age spectra are not avail- able in this case):
where h x i , h y i and h z i are the kernel bandwidths which determine the degree of smoothing, rx=(Xi-x)/hxi, ry=(Yi-y)/hyi, rz=Zi/hzi with Xi, Yi and Zi being the position of particle i .  The kernel bandwidths are a function of the particles' age.
8.5
Mass fluxes
Mass flux calculations can be switched on using iflux in file COMMAND. Mass fluxes are calculated separately for east- ward, westward, northward, southward, upward and down- ward directions and contain both grid-scale and subgrid-scale motions.  Mass fluxes are determined for the centerlines of the output grid cells, e.g. vertical fluxes are calculated for motions across the half level of each output cell.
Domain-filling option
9
9.1
General
If mdomainfill=1 in file COMMAND particles are not re- leased at specific locations. Instead, the longitudes and lat- itudes specified for the first release in the RELEASES file are used to set up a global or limited model domain.  The particles (number is also taken from RELEASES ) are then distributed in the model domain proportionally to air den- sity (subroutine init_domainfill.f ). Each particle re- ceives the same mass, altogether accounting for the total at- mospheric mass.  Subsequently, particles move freely in the atmosphere. If a limited domain is chosen, mass fluxes are determined in small grid boxes at the boundary of this domain (bound- aries must be at least one grid box away from the bound- aries of the meteorological input data). In the grid cells with air flowing into the model domain, mass fluxes are accumu- lated over time and whenever the accumulated mass exceeds the mass of a particle, a new particle (or more, if required) is released at a randomly chosen position at the boundary of the box (subroutine boundcond_domainfill.f). At the outflowing boundaries particles are terminated. Note that, due to the change of mass of the atmosphere in the model domain and due to numerical effects, the number of particles used is not exactly constant throughout the simulation.
9.2
Stratospheric ozone tracer
If mdomainfill=2 , the domain-filling option is used to simulate a stratospheric ozone tracer. Upon particle creation, the potential vorticity (PV) at its position is determined by interpolation from the ECMWF data.  Particles initially lo- cated in the troposphere (PV<pvcrit potential vorticity units (pvu), default 2 pvu) are not used.  In contrast, strato- spheric particles (PV>pvcrit) are given a mass according to:
Mo, = Mair P C 48/29
(52) where Mair is the mass of air a particle represents, P is PV in pvu, C=60×10−9 pvu−1 is the ozone/PV relationship (Stohl et al., 2000) (parameter ozonescale), and the factor 48/29 converts from volume to mass mixing ratio. Particles are then allowed to advect through the stratosphere and into the tro- posphere according to the winds.
Model output
10
Tracer concentrations and/or mixing ratios (for forward runs), or emission sensitivity response functions (for back- ward runs) are calculated on a three-dimensional longitude- latitude grid, defined in file OUTGRID , whose domain and resolution can differ from the grid on which meteorological input data are given. Two-dimensional wet and dry deposi- tion fields are calculated over the same spatial domain, and tracer mass fluxes can also be determined on the 3-d grid. Except for the mass fluxes, output can also be produced on one nested output grid with higher horizontal but the same vertical resolution, defined in file OUTGRID_NEST . For cer- tain locations, specified in file RECEPTORS , concentrations can also be calculated independently from a grid (see below). The time interval (variable loutstep) at which output is produced is read in from file COMMAND. For every output time, files are created, whose file name ends with the date and time in the format yyyymmddhhmmss.  A list of all these output times is written to the formatted file dates. The dates indicate the ending time of an output sampling in- terval (see Sect. 8).
10.1
Gridded output
There are several output options in FLEXPART, which can all be selected in file COMMAND. Gridded output fields can be concentrations (files grid_conc_date ), volume mixing ratios (files grid_pptv_date ), emission response sensi- tivity in backward simulations (files grid_time_date ), or fluxes (files grid_flux_date , unit 10-12kgm-2s-1 for forward runs). Files grid_conc_date are created only in forward runs, whereas files grid_time_date are only created in backward runs.  Note that the units of the files grid_conc_date and grid_time_date de- pend on the settings of the switches ind_source and ind_receptor , following Table 1.  In particular, the units of grid_conc_date can also be mass mixing ra- tios. For forward runs, additional files grid_pptv_date can be created, which contain volume mixing ratios for Output files grid_conc*, grid_pptv*, and gases. grid_time* also contain wet and dry deposition fields (unit 10-12 kg m-2 in forward mode), and all files contain, for each grid cell, corresponding uncertainties.  All these file types share a common header, file header produced by subroutine writeheader.f , where important information on the model run (start of simulation, grid domain, number and position of vertical levels, age classes, release points, etc.) is stored. In all postprocessing programs, the header must be read in before the actual data files.  File names for the output nests follow the same nomenclature as de- scribed above, but with _nest added (e.g., header_nest , or grid_conc_nest_date ). The output files are written with subroutines concoutput .f and fluxoutput .f . FLEXPART output files, except for dates, are all binary and often contain many grid cells with zero concentrations (or mixing ratios, fluxes, etc.). Writing out only those cells with non-zero values can produce smaller output than a full grid dump.  But in this case the grid indices (note that all three are combined into a single integer number) must also be written out and this produces bigger output than a full grid dump if most grid cells contain non-zero concentrations. Therefore, at every output time and for every output field the more efficient method is determined and used.
Receptor point output
10.2
For a list of points at the surface, concentrations or mix- ing ratios in forward simulations can be determined with a grid-independent method. This information is written to files receptor_conc and receptor_pptv , respec- tively, for all dates of a simulation.
10.3
Particle dump and warm start option
Particle information (3-d position, release time, release point, and release masses for all species) can be written out to files (subroutine partoutput.f) either continuously (binary files partposit_date ), or 'only at the end' of a simu- lation (file partposit_end ). In both cases output is writ- ten every output interval but file partposit_end is over- written upon each new output. If FLEXPART must be termi- nated, it can be continued later on by reading in files header and partposit_end produced by the previous run (sub- routine readpartpositions.f ). Such a warm start is done if variable ipin is set to 1 in file COMMAND. If option mquasilag is chosen in file COMMAND, particle dumps every output interval are produced in a very compact format by converting the positions to an integer*2 for- mat (subroutine partoutput_short.f ). As some accu- racy is lost in the conversion, this output is not used for the warm start option.  Another difference to the normal parti- cle dump is that every particle gets a unique number, thus allowing postprocessing routines to identify continuous par- ticle trajectories.
10.4
Clustered plume trajectories Condensed particle output using the clustering algorithm described in Sect. 6 is written to the formatted file trajectories.txt. Information on the release points (coordinates, release start and end, number of particles) is written by subroutine openouttraj.f to the beginning of file trajectories.txt. Subsequently, plumetraj.f writes out a time sequence of the clustering results for each release point: release point number, time in seconds elapsed since the middle of the release interval, plume centroid po- sition coordinates, various overall statistics (e.g., fraction of particles residing in the ABL and troposphere), and then for each cluster the cluster centroid position, the fraction of par- ticles belonging to the cluster, and the root-mean-square dis- tance of cluster member particles from the cluster centroid.
11
Final remark
In this note, we have described the particle dispersion model FLEXPART version 6.2 with the dual purpose of creating a citeable reference for FLEXPART and pro- viding an actual user manual.  In an appendix provided as an electronic supplement (http://www.atmos-chem- phys.org/acp/5/2461/acp-5-2461-sp.pdf), the various FLEXPART input files are briefly explained and exam- ples are given.  As FLEXPART develops this text will be kept actual and will be accessible from the internet site http://zardoz.nilu.no/~andreas/flextra+flexpart.html. Edited by: B. Kärcher