# The Lagrangian particle dispersion model FLEXPART version 10.4

Ignacio Pisso, Espen Sollum, Henrik Grythe, Nina Kristiansen, Massimo Cassiani, Sabine Eckhardt, Delia Arnold, Don Morton, Rona Thompson, Christine Groot Zwaaitink, et al.

► To cite this version:

Ignacio Pisso, Espen Sollum, Henrik Grythe, Nina Kristiansen, Massimo Cassiani, et al. The Lagrangian particle dispersion model FLEXPART version 10.4. Geoscientific Model Development Discussions, 2019, 12, pp.4955-4997. 10.5194/gmd-12-4955-2019. hal-03135256

HAL Id: hal-03135256

https://hal.univ-reunion.fr/hal-03135256

Submitted on 16 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Lagrangian particle dispersion model FLEXPART version 10.4

Ignacio Pisso<fnref n="1" />, Espen Sollum<fnref n="1" />, Henrik Grythe<fnref n="1" />, Nina I. Kristiansen<fnref n="1" />,ª, Massimo Cassiani<fnref n="1" />, Sabine Eckhardt<fnref n="1" />, Delia Arnold<fnref n="2" />,<fnref n="3" />, Don Morton<fnref n="4" />, Rona L. Thompson<fnref n="1" />, Christine D. Groot Zwaaitink<fnref n="1" />, Nikolaos Evangelou<fnref n="1" />, Harald Sodemann<fnref n="5" />, Leopold Haimberger<fnref n="6" />, Stephan Henne<fnref n="7" />, Dominik Brunner<fnref n="7" />, John F. Burkhart<fnref n="8" />, Anne Fouilloux<fnref n="8" />, Jerome Brioude<fnref n="9" />, Anne Philipp<fnref n="6" />,<fnref n="10" />, Petra Seibert<fnref n="11" />, and Andreas Stohl<fnref n="1" />

**Correspondence:** Ignacio Pisso (ip@nilu.no)

Received: 21 December 2018 – Discussion started: 28 January 2019

Revised: 25 July 2019 – Accepted: 7 August 2019 – Published: 2 December 2019

**Abstract.** The Lagrangian particle dispersion model FLEXPART in its original version in the mid-1990s was designed for calculating the long-range and mesoscale dispersion of hazardous substances from point sources, such as those released after an accident in a nuclear power plant. Over the past decades, the model has evolved into a comprehensive tool for multi-scale atmospheric transport modeling and analysis and has attracted a global user community. Its application fields have been extended to a large range of atmospheric gases and aerosols, e.g., greenhouse gases, short-lived climate forcers like black carbon and volcanic ash, and it has also been used to study the atmospheric branch of the water cycle. Given suitable meteorological input data, it can be used for scales from dozens of meters to global. In particular, inverse modeling based on source–receptor relationships from FLEXPART has become widely used. In this paper, we present FLEXPART version 10.4, which works with meteorological input data from the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecast System (IFS) and data from the United States National Cen-

ters of Environmental Prediction (NCEP) Global Forecast System (GFS). Since the last publication of a detailed FLEXPART description (version 6.2), the model has been improved in different aspects such as performance, physicochemical parameterizations, input/output formats, and available preprocessing and post-processing software. The model code has also been parallelized using the Message Passing Interface (MPI). We demonstrate that the model scales well up to using 256 processors, with a parallel efficiency greater than 75% for up to 64 processes on multiple nodes in runs with very large numbers of particles. The deviation from 100% efficiency is almost entirely due to the remaining non-parallelized parts of the code, suggesting large potential for further speedup. A new turbulence scheme for the convective boundary layer has been developed that considers the skewness in the vertical velocity distribution (updrafts and downdrafts) and vertical gradients in air density. FLEXPART is the only model available considering both effects, making it highly accurate for small-scale applications, e.g., to quantify dispersion in the vicinity of a point source. The wet deposi-

<footnote n="1">¹Norwegian Institute for Air Research (NILU), Kjeller, Norway</footnote>

<footnote n="2">²Central Institute for Meteorology and Geodynamics (ZAMG), Vienna, Austria</footnote>

<footnote n="3">³Arnold Scientific Consulting, Manresa, Spain</footnote>

<footnote n="4">⁴Boreal Scientific Computing, Fairbanks, Alaska, USA</footnote>

<footnote n="5">⁵Geophysical Institute, University of Bergen and Bjerknes Centre for Climate Research, Bergen, Norway</footnote>

<footnote n="6">⁶Department of Meteorology and Geophysics, University of Vienna, Vienna, Austria</footnote>

<footnote n="7">⁷Empa, Swiss Federal Laboratories for Materials Science and Technology, Dübendorf, Switzerland</footnote>

<footnote n="8">⁸Department of Geosciences, University of Oslo, Oslo, Norway</footnote>

<footnote n="9">⁹Laboratoire de l'Atmosphère et des Cyclones (LACy), UMR8105, Université de la Réunion – CNRS – Météo-France, Saint-Denis de La Réunion, France</footnote>

<footnote n="10">¹⁰Aerosol Physics & Environmental Physics, University of Vienna, Vienna, Austria</footnote>

<footnote n="11">¹¹Institute of Meteorology, University of Natural Resources and Life Sciences, Vienna, Austria</footnote>

<footnote>ªnow at: Met Office, FitzRoy Road, Exeter, EX1 3PB, UK</footnote>
tion scheme for aerosols has been completely rewritten and a new, more detailed gravitational settling parameterization for aerosols has also been implemented. FLEXPART has had the option of running backward in time from atmospheric concentrations at receptor locations for many years, but this has now been extended to also work for deposition values and may become useful, for instance, for the interpretation of ice core measurements. To our knowledge, to date FLEXPART is the only model with that capability. Furthermore, the temporal variation and temperature dependence of chemical reactions with the OH radical have been included, allowing for more accurate simulations for species with intermediate lifetimes against the reaction with OH, such as ethane. Finally, user settings can now be specified in a more flexible namelist format, and output files can be produced in NetCDF format instead of FLEXPART's customary binary format. In this paper, we describe these new developments. Moreover, we present some tools for the preparation of the meteorological input data and for processing FLEXPART output data, and we briefly report on alternative FLEXPART versions.

# 1 Introduction

Multi-scale offline Lagrangian particle dispersion models (LPDMs) are versatile tools for simulating the transport and turbulent mixing of gases and aerosols in the atmosphere. Examples of such models are the Numerical Atmospheric dispersion Modelling Environment (NAME) (Jones et al., 2007), the Stochastic Time-Inverted Lagrangian Transport (STILT) model (Lin et al., 2003), the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Stein et al., 2015) and the FLEXible PARTicle (FLEXPART) model (Stohl et al., 1998, 2005). LPDMs are stochastic models that compute trajectories for a large number of notional particles that do not represent real aerosol particles but points moving with the ambient flow. The trajectories represent the transport by mean flow as well as turbulent, diffusive transport by unresolved parameterized subgrid-scale transport processes (e.g., turbulence, meandering, deep convection, etc.) and can also include gravitational settling. Each particle carries a certain mass, which can be affected by loss processes such as radioactive decay, chemical loss, or dry and wet deposition.

The theoretical basis for most currently used atmospheric particle models was laid down by Thomson (1987). He introduced the criterion to formulate Lagrangian stochastic models that produce particle trajectories consistent with pre-defined Eulerian probability density functions in physical and velocity space. Rodean (1996) and Wilson and Sawford (1996) provided detailed descriptions of the theory and formulation of LPDMs in constant density flows and under different atmospheric stability conditions. Stohl and Thomson (1999) extended this to flows with vertically variable air den-

sity. An important characteristic of LPDMs is their ability to run backward in time in a framework that is theoretically consistent with both the Eulerian flow field and LPDM forward calculations. This was discussed by Thomson (1987, 1990), further developed by Flesch et al. (1995), and extended to global-scale dispersion by Stohl et al. (2003) and Seibert and Frank (2004). The more practical aspects and efficiency of LPDMs were discussed by Zannetti (1992) and Uliasz (1994). A history of their development was provided by Thomson and Wilson (2013).

Lagrangian models exhibit much less numerical diffusion than Eulerian or semi-Lagrangian models (e.g., Reithmeier and Sausen, 2002; Cassiani et al., 2016), even though some numerical errors also arise in the discretization of their stochastic differential equations (Ramli and Esler, 2016). Due to their low level of numerical diffusion, tracer filaments generated by dispersion in the atmosphere (Ottino, 1989) are much better captured in Lagrangian models than in Eulerian models. It has been noticed, for instance, that Eulerian models have difficulties simulating the fine tracer structures created by intercontinental pollution transport (Rastigejev et al., 2010), while these are well preserved in LPDMs (e.g., Stohl et al., 2003). 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 computational grid and can account for point or line sources with potentially infinitesimal spatial resolution. When combined with their capability to run backward in time, this means that LPDMs can also be used to investigate the history of air parcels affecting, for instance, an atmospheric measurement site (e.g., for in situ monitoring of atmospheric composition).

The computational efficiency of LPDMs depends on the type of application. One important aspect is that their computational cost does not increase substantially with the number of species transported (excluding aerosol particles with different gravitational settling, for which trajectories deviate from each other), making multispecies simulations efficient. On the other hand, the computational time scales linearly with the number of particles used, while the statistical error in the model output decreases only with the square root of the particle density. Thus, it can be computationally costly to reduce statistical errors, and data input/output can require substantial additional resources. Generally, a high particle density can be achieved with a small number of released particles in the vicinity of a release location, where statistical errors, relative to simulated concentrations, are typically small. However, particle density and thus the relative accuracy of the results decrease with distance from the source. Methods should therefore be used to reduce the statistical error (e.g., Heinz et al., 2003), such as kernels or particle splitting, and it is important to quantify the statistical error.
## 1.1 The Lagrangian particle dispersion model
FLEXPART

One of the most widely used LPDMs is the open-source model FLEXPART, which simulates the transport, diffusion, dry and wet deposition, radioactive decay, and 1st-order chemical reactions (e.g., OH oxidation) of tracers released from point, line, area or volume sources, or filling the whole atmosphere (Stohl et al., 1998, 2005). FLEXPART development started more than 2 decades ago (Stohl et al., 1998) and the model has been free software ever since it was first released. The status as a free software is formally established by releasing the code under the GNU General Public License (GPL) Version 3. However, the last peer-reviewed publication describing FLEXPART (version 6.2) was published as a technical note about 14 years ago (Stohl et al., 2005). Since then, while updates of FLEXPART's source code and a manual were made available from the web page at https://flexpart.eu/ (last access: 30 October 2019), no citable reference was provided. In this paper, we describe FLEXPART developments since Stohl et al. (2005), which led to the current version 10.4 (subsequently abbreviated as v10.4).

FLEXPART can be run either forward or backward in time. For forward simulations, particles are released from one or more sources and concentrations (or mixing ratios) are determined on a regular latitude–longitude–altitude grid. In backward mode, the location where particles are released represents a receptor (e.g., a measurement site). Like in the forward mode, particles are sampled on a latitude–longitude–altitude grid, which in this case corresponds to potential sources. The functional values obtained represent the source–receptor relationship (SRR) (Seibert and Frank, 2004), also called source–receptor sensitivity (Wotawa et al., 2003) or simply emission sensitivity, and are related to the particles’ residence time in the output grid cells. Backward modeling is more efficient than forward modeling for calculating SRRs if the number of receptors is smaller than the number of (potential) sources. Seibert and Frank (2004) explained in detail the theory of backward modeling, and Stohl et al. (2003) gave a concrete backward modeling example. FLEXPART can also be used in a domain-filling mode whereby the entire atmosphere is represented by (e.g., a few million) particles of equal mass (Stohl and James, 2004). The number of particles required for domain-filling simulations, not unlike those needed for other types of simulations, depends on the scientific question to be answered. For instance, a few million particles distributed globally are often enough to investigate the statistical properties of air mass transport (e.g., monthly average residence times in a particular area that is not too small) but would not be enough for a case study of airstreams related to a particular synoptic situation (e.g., describing flow in the warm conveyor belt of a particular cyclone).

FLEXPART is an offline model that uses meteorological fields (analyses or forecasts) as input. Such data are available from several different numerical weather prediction (NWP)

models. For the model version described here, v10.4, data from the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecast System (IFS) and data from the United States National Centers of Environmental Prediction (NCEP) Global Forecast System (GFS) can be used. Common spatial resolutions for IFS depending on the application include 1°×1° at 3 h (standard for older products, e.g., ERA-Interim), 0.5°×0.5° at 1 h (standard for newer products, e.g., ERA5) and 0.1°×0.1° at 1 h (current ECMWF operational data). The ECMWF IFS model currently has 137 vertical levels. NCEP GFS input files are usually used at 1°×1° horizontal resolution, with 64 vertical levels and 3 h time resolution. NCEP GFS input files are also available at 0.5°×0.5° and 0.25°×0.25° horizontal resolution. Other FLEXPART model branches have been developed for input data from various limited-area models, for example the Weather Research and Forecasting (WRF) meteorological model (Brioude et al., 2013) and the Consortium for Small-scale Modeling (COSMO) model (Oney, 2015), which extend the applicability of FLEXPART down to the mesogamma scale. Notice that the turbulence parameterizations of FLEXPART are valid at even smaller scales. Another FLEXPART model version, FLEXPART–NorESM/CAM (Cassiani et al., 2016), uses the meteorological output data generated by the Norwegian Earth System Model (NorESM1-M) with its atmospheric component CAM (Community Atmosphere Model). The current paper does not document these other model branches, but most share many features with FLEXPART v10.4 and some are briefly described in Appendix C. A key aspect of these model branches is the ability to read meteorological input other than that from ECMWF or NCEP.

## 1.2 FLEXPART and its history

FLEXPART's first version (v1) was a further development of the trajectory model FLEXTRA (Stohl et al., 1995) and was coded in Fortran 77. It provided gridded output of concentrations of chemical species and radionuclides. Its meteorological input data were based on the ECMWF-specific GRIB-1 (gridded binary) format. The model was first applied in an extensive validation study using measurements from three large-scale tracer experiments (Stohl et al., 1998). A deposition module was added in version 2. Version 3 saw improvements in performance and the addition of a subgrid-scale terrain effects parameterization. In v3.1 the output format was optimized (sparse matrix) and mixing ratio output could optionally be produced. It also allowed for the output of particle positions. Furthermore, a density correction was added to account for decreasing air density with height in the boundary layer (Stohl and Thomson, 1999). Further v3 releases included the addition of a convection scheme (Seibert et al., 2001) based on Emanuel and Živković-Rothman (1999), the option to calculate mass fluxes across grid cell faces and age spectra, and free format input (v3.2). The preliminary convection scheme of v3.2 was revised in v4 (see
Forster et al., 2007). In v5 the output unit of backward calculations was changed to seconds and improvements in the input/output handling were made. Comprehensive validation of these early FLEXPART versions was done during intercontinental air pollution transport studies at the end of the 1990s and early 2000s (Stohl and Trickl, 1999; Forster et al., 2001; Spichtinger et al., 2001; Stohl et al., 2002, 2003). Special developments were also made in order to extend FLEXPART's forecasting capabilities for large-scale field campaigns (Stohl et al., 2004). Version 6.0 saw corrections to the numerics in the convection scheme, the addition of a domain-filling option used, for instance, in water cycle studies (Stohl and James, 2004) and the possibility to use nested output. Version 6.2, which added the ability to model sources and receptors in both mass and mixing ratio units (Seibert and Frank, 2004), is currently the last version described in a publication (Stohl et al., 2005). A separate sister model branch (v6.4) was adapted to run with NCEP GFS meteorological input data. The current paper describes the most important model developments since v6.2 (for ECMWF) and v6.4 (for GFS).

Version 8.0 unified the model branches based on ECMWF IFS and NCEP GFS input data in one source package but still required the building of two different executables. Importantly, Fortran 90 constructs were introduced in parts of the code, such as initial support for dynamic memory allocation. Furthermore, a global land use inventory was added, allowing for more accurate dry deposition calculations everywhere on the globe (before, land use data were provided only for Europe). The reading of the – at the time – newly introduced GRIB-2 format with the ECMWF grib_api library was implemented in v8.2. An option to calculate the sensitivity to initial conditions in backward model runs (in addition to the emission sensitivity calculations) was also implemented in v8.2. Version 8 was also the first version that distinguished between in-cloud and below-cloud scavenging for washout, relying on simple diagnostics for clouds based on grid-scale relative humidity. With a growing number of parameters defining removal processes, each species was given its own definition file, whereas in previous versions the properties for all species were contained in a single file. The gravitational settling scheme was improved in v8.2.1 (Stohl et al., 2011), and this is briefly described in this paper in section 2.3.

For v9, the code was transformed to the Fortran 90 free-form source format. The option to read the run specifications from Fortran namelists instead of the standard input files was introduced, as described in Sect. 5 of this paper. This change was motivated by the resulting greater flexibility, in particular with regard to setting default values, optional arguments, when new developments require adding new parameters and when specifying parameter lists. In addition, an option to produce output in compressed NetCDF 4 format was provided (see Sect. 6.3). Another option to write some model output only for the first vertical layer to save storage

space for inverse modeling applications was also introduced (Thompson and Stohl, 2014) (see Sect. 2.6).

## 1.3 FLEXPART version 10.4

For v10.4 of FLEXPART, described in this paper, several more changes and improvements were made. First, an optional new scheme applying more realistic skewed rather than Gaussian turbulence statistics in the convective atmospheric boundary layer (ABL) was developed (Sect. 2.1). Second, the wet deposition scheme for aerosols was totally revised (Grythe et al., 2017), introducing dependencies on aerosol size, precipitation type (rain or snow), and distinguishing between in-cloud and below-cloud scavenging (see Sect. 2.4). The code now also allows for the reading of three-dimensional (3-D) cloud water fields from meteorological input files. Third, a method to calculate the sensitivity of deposited quantities to sources in backward mode was developed (Sect. 2.5). Fourth, chemical reactions with the hydroxyl radical (OH) are now made dependent on the temperature and vary sub-monthly (Sect. 2.7). Fifth, large parts of the code were parallelized using the Message Passing Interface (MPI), thus facilitating a substantial speedup for certain applications (see Sect. 3), and the code was unified so that a single executable can now use both ECMWF and GFS input data. Sixth, a dust mobilization scheme that can be run as a FLEXPART preprocessor was developed (Sect. 2.8). Seventh, the software used to retrieve data from the ECMWF has been modernized and can now also be used by scientists from non-ECMWF member states (Sect. 5.2.1). Finally, a testing environment was created that allows users to verify their FLEXPART installation and compare results (Sect. 7).

Despite the many changes and additions, in large part the operation of FLEXPART v10.4 still resembles the original version 1 design. Throughout this paper, we avoid repeating information on aspects of the model that have not changed since earlier model descriptions. The paper should therefore always be considered together with the publications of Stohl et al. (1998, 2005). To provide the necessary context for the rest of this paper, we provide a brief overview of the FLEXPART v10.4 directory structure in Table 1. The source code is contained in directory *src*. The pathnames of the input and output directories are stated in the file *pathnames* read by the FLEXPART executable. The directory *options* contains the parameters that define a run in files such as *COMMAND* (e.g., start and end times of the simulation, output frequency, etc.), *RELEASES* (definition of the particle releases), *OUTGRID* (output grid specifications) and others. All the output is written in a directory unique for each run. There are also other directories containing the model testing environment and example runs, as well as preprocessing and post-processing software (see Table 1).

Sensu stricto FLEXPART consists of the (Fortran) source files required to build an executable, not including external libraries such as those needed for input reading. The make-
files and the sample input as provided in the “options” may also be included under this term. However, in order to do real work with FLEXPART, one also needs to obtain meteorological input data (in the case of the ECMWF this is not trivial, so the *flex_extract* package is provided), and one needs to do post-processing. This is the reason why we include a selection of such tools here.

## 2 Updates of the model physics and chemistry

This section gives an overview of the main updates of the model physics and chemistry since the last published FLEXPART version, v6.2 (Stohl et al., 2005). Some developments have been published already separately, and in such cases we keep the description short, focusing on technical aspects of the implementation in FLEXPART that are important for model users or demonstrating applications not covered in the original papers.

### 2.1 Boundary layer turbulence

Subgrid-scale atmospheric motions unresolved by the meteorological input data need to be parameterized in FLEXPART. This is done by adding stochastic fluctuations based on Langevin equations for the particle velocity components (Stohl et al., 2005). In the ABL, the stochastic differential equations are formulated according to the well-mixed criteria proposed by Thomson (1987). Until FLEXPART version 9.2, the Eulerian probability density functions (PDFs) for the three velocity components were assumed to be three independent Gaussian PDFs. However, for the vertical velocity component, the Gaussian turbulence model is well suited only for stable and neutral conditions. In the convective ABL (CBL), turbulence is skewed since a larger area is occupied by downdrafts than by updrafts (e.g., Stull, 1988; Luhar and Britter, 1989). In such conditions, the Gaussian turbulence model is not appropriate for sources within the ABL, as it cannot reproduce the observed upward bending of plumes from near-ground sources or the rapid downward transport of plumes from elevated sources (Venkatram and Wyngaard, 1988). However, the Gaussian approximation has negligible influence once the tracer is mixed throughout the whole ABL.

Cassiani et al. (2015) developed an alternative Langevin equation model for the vertical particle velocity including both skewed turbulence and a vertical density gradient, which is now implemented in FLEXPART v10.4. This scheme can be activated by setting the switch CBL to 1 in the file COMMAND. In this case, the time step requirement for numerical stability is much more stringent than for the standard Gaussian turbulence model (typically, values of CTL=10 and IFINE=10 are required, also set in the file COMMAND). Therefore, also considering that the computation time required for each time step is about 2.5 times that of

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Comparison of FLEXPART results obtained with the skewed turbulence parameterization (a) and with the Gaussian turbulence parameterization (b). Shown are the tracer concentrations integrated over all latitudes as a function of altitude and longitude. The simulations used a point source emitting 100 kg of tracer per hour for a period starting on 1 July 2017 at 12:00 UTC and ending at 13:30 UTC. The source was located near Vienna (Austria) at 47.9157° N and 16.3274° E, 250 m above ground level. Results are averaged for the time period 12:40 to 13:00 UTC. Notice that the maximum ground-level concentration in panel (a) (ca. 7 mg m⁻²) is about 30 % higher than in panel (b) (5 mg m⁻²).</figcaption></figure>

the standard Gaussian formulation, the CBL option is much more computationally demanding and not recommended for large-scale applications. However, for studies of tracer dispersion in the vicinity of individual point sources, the CBL option is essential to reproduce the characteristic features of CBL dispersion (Weil, 1985), while the additional computational burden remains tolerable.

Figure 1 shows a comparison between two simulations of dispersion from an elevated source, with the skewed and with the Gaussian turbulence model. It can be seen that the maximum time-averaged ground-level concentration is about 30 % higher for the skewed turbulence parameterization. This is the result of the plume centerline tilting downward to the surface in the vicinity of the source for the skewed turbulence case due to downdrafts being more frequent than updrafts. The plume also spreads faster in this case. These results are similar to those obtained by others (e.g., Luhar and Britter, 1989).

It is important to note that the CBL formulation smoothly transits to a standard Gaussian formulation when the stabil-
<table><caption>**Table 1.** Directory structure overview of the FLEXPART v10.4 software distribution. All listed directories are subdirectories of the installation root directory `$flexhome/`.</caption><thead><tr><th></th><th>Subdirectory or file</th><th>Contents</th><th>Description and comments</th></tr></thead><tbody><tr><td rowspan="10">$flexhome/</td><td>pathnames<br>(file)</td><td>$options/<br>$output/<br>$flex_winds/<br>$AVAILABLE</td><td>path to options directory<br>path to output directory<br>path to meteorological input files<br>path to AVAILABLE file</td></tr><tr><td>src/</td><td>*.f90<br>makefile<br>FLEXPART</td><td>Fortran source files<br>see Sect. 4 and Appendix A<br>executable file (see Sect. 4)</td></tr><tr><td>options/</td><td>COMMAND, RELEASES, OUTGRID, SPECIES, AGECLASSES, OUTGRID_NEST, RECEPTORS, IGBP_int1.dat, surfdata.t, surfdepo.t, OH_variables.bin</td><td>user input files<br>see Sect. 5 and Table 7</td></tr><tr><td>AVAILABLE</td><td>list of meteorological input data files</td><td>file containing list, see Sect. 5</td></tr><tr><td>output/</td><td>FLEXPART output files</td><td>see Sect. 6 and Table 11</td></tr><tr><td>preprocess/</td><td>flex_extract/</td><td>see Sect. 5.2</td></tr><tr><td>postprocess/</td><td>read_flex_fortran/, read_flex matlab/,</td><td>see Sect. 6.4</td></tr><tr><td>tests/<br>tests/examples/</td><td>development tests for FLEXPART and ancillary software<br>example runs illustrating various FLEXPART functionalities</td><td>see Sect. 7<br>and Appendix C</td></tr></tbody></table>

ity changes towards neutral (Cassiani et al., 2015). However, the actual equation solved inside the model for the Gaussian condition is still different from the standard version as actual particle velocities rather than the scaled ones are advanced (see, e.g., Wilson et al., 1981; Rodean, 1996). Full details of the CBL implementation can be found in Cassiani et al. (2015).

To date, FLEXPART has mainly been used for large-scale applications. With this new CBL option, FLEXPART is now also well suited for the simulation of small-scale tracer dispersion or for the inverse modeling of point source emissions from near-field measurements – at least if the resolved winds are representative of the situation considered. In fact, to our knowledge FLEXPART is the only particle model considering both skewness in the vertical velocity distribution and vertical gradients in air density. Both these effects are particularly important in deep CBLs and can be additive with respect to simulated ground-level concentrations.

## 2.2 Turbulent diffusivity above the boundary layer

Above the atmospheric boundary layer, turbulent fluctuations can be represented with a turbulent diffusivity. The value of the diffusivity tensor controls the size and lifetimes of the filamentary structures caused by advection. Diffusivities are converted into velocity scales using $\sigma_{v_i} = \sqrt{2D_i/dt}$, where i is the direction. This corresponds to a mean diffusive displacement of $\sigma_{x_i} = \sqrt{2D_i dt}$, characteristic of Brownian motion. For i, only the vertical (v) and horizontal (h) directions are considered. The value of the vertical diffusivity $D_z$ is related

to the value of the horizontal diffusivity $D_h$ by the square of the typical atmospheric aspect ratio for tracer structures $\kappa \approx 100-250$ (Haynes and Anglade, 1997).

FLEXPART uses by default a constant vertical diffusivity $D_z = 0.1 \text{ m}^2 \text{ s}^{-1}$ in the stratosphere, following Legras et al. (2003), whereas a horizontal diffusivity $D_h = 50 \text{ m}^2 \text{ s}^{-1}$ is used in the free troposphere. In general in the atmosphere, the values of the turbulent diffusivity depend on time and location, showing in particular seasonal, latitudinal and altitude variability: e.g., $D_v = 10^{-2} \text{ m}^2 \text{ s}^{-1}$ in the stratosphere (Balluch and Haynes, 1997) and $D_h = 10^4 \text{ m}^2 \text{ s}^{-1}$ (Pisso et al., 2009) in the troposphere. The values can be modified by the user in the `COMMAND` file (namelist variables `d_trop` and `d_strat`, in $\text{m}^2 \text{ s}^{-1}$). As mentioned above, $D_h \approx \kappa^2 D_z$, and therefore both values can be used interchangeably.

In FLEXPART version 6.2, the stratosphere and troposphere were distinguished based on a threshold of 2 pvu (potential vorticity units), with a maximal height of 18 km in the tropics and a minimal height of 5 km elsewhere. Such a threshold is well suited to midlatitudes, but it can differ from the thermal tropopause in the polar regions and close to the Equator. In FLEXPART 10.4, the thermal tropopause definition is used and is calculated using the lapse rate definition (Hoinka, 1997).

## 2.3 Gravitational settling

Gravitational settling of aerosols is implemented in FLEXPART as a process that changes the particle trajectories. The settling velocity is determined at each time step and added
to the vertical wind velocity. In simulations in which a particle represents several species, all species are transported with the settling velocity of the first species. If this is not intended, simulations for the different species must be run separately. Gravitational settling velocities are also used in the calculation of dry deposition.

In older FLEXPART versions, gravitational settling was calculated using a single dynamic viscosity of air. With FLEXPART 8.2.1, the gravitational settling calculation was generalized to higher Reynolds numbers and it takes into account the temperature dependence of dynamic viscosity. This is done in subroutine `get_settling.f90` in an iterative loop, wherein both the Reynolds number and the settling velocity are determined (Naeslund and Thaning, 1991). For initialization of the loop, Stokes' law and a constant viscosity estimate is used. The dynamic viscosity is calculated as a function of temperature using the formula of Sutherland (1893). A spherical shape of the particles is assumed in the settling calculation, which could be further extended in the future to allow for more complex particle shapes. For particle sizes of about 10 µm, settling velocities in the new FLEXPART version are not much different from earlier versions using the old settling calculation, typically less than 20%. However, the differences are largest in the cold upper troposphere, implying, for instance, changes in the residence time of volcanic ash particles at heights relevant for aviation. The residence times in the upper troposphere are increased with the new scheme, but the effect is not particularly large, typically on the order of 20%.

## 2.4 Wet deposition

In FLEXPART, the calculation of wet scavenging is divided into three parts. First, where scavenging occurs and which form it takes is determined (e.g., below- or within-cloud scavenging). Second, the scavenging coefficient is determined. Third, the actual removal of particle mass is calculated.

With respect to the first part, it is important to understand that wet scavenging occurs only in the presence of clouds and where precipitation occurs. In air columns without clouds, above the top of the clouds, and where neither the large-scale nor the convective precipitation rate exceeds 0.01 mm h⁻¹, no scavenging occurs. To quickly know where a particle is located relative to the clouds, in subroutines `verttransform_ecmwf.f90` and `verttransform_gfs.f90` each grid cell is categorized as being in a cloud-free column, above a cloud, inside a cloud or below a cloud. This cloud definition has been completely revised compared to earlier versions and is described in Sect. 2.4.1.

With respect to the second step, the scavenging coefficient Λ (s⁻¹) is determined in subroutine `get_wetscav.f90`. After a series of updates, in particular Grythe et al. (2017), FLEXPART now distinguishes between below-cloud and in-

cloud scavenging and also has different parameterizations of Λ for gases and particles. For the latter, it also distinguishes between liquid-phase and ice-phase states. This yields in total six different parameterizations for Λ, described in Sect. 2.4.2 and 2.4.3.

In the third step, the removal of particle mass due to wet deposition is calculated. It takes the form of an exponential decay process (McMahon, 1979),

$$m(t + \Delta t) = m(t) \exp(-\Lambda \Delta t), \quad (1)$$

where m is the particle mass (kg) (it can also be a mass mixing ratio, depending on settings in file `COMMAND`). This removal of particle mass and corresponding accumulation of deposition at the surface is calculated in subroutine `wetdepo.f90` and has not been changed since earlier versions.

### 2.4.1 Definition of clouds, cloud water content and precipitation

The location of clouds, the total cloud column water content and phase, and precipitation intensity and phase are needed in the calculation of the wet scavenging. Therefore, a three-dimensional cloud mask is defined in subroutine `verttransform_ecmwf.f90` (or `verttransform_gfs.f90` for GFS data). In previous FLEXPART versions, the cloud definition scheme was very simple and based on relative humidity (RH). In grid columns with precipitation, grid cells with RH > 80% were defined as in-cloud, and those with RH < 80% were set as below-cloud up to the bottom of the uppermost layer with RH > 80%. This was appropriate for the first version of FLEXPART, as the ECMWF had a similarly simple definition of clouds and more detailed information was not available from the ECMWF archives at the time.

If no cloud information is available from the meteorological data, the old RH-based scheme is still used in FLEXPART. However, nowadays, specific cloud liquid water content (CLWC; kg kg⁻¹) and specific cloud ice water content (CIWC; kg kg⁻¹) are available as 3-D fields in meteorological analyses from the ECMWF, and NCEP also provides the 3-D cloud water (liquid plus ice) mixing ratio (CLWMR; kg kg⁻¹), from here on referred to as $q_c$. A cloudy grid cell is defined when $q_c > 0$. FLEXPART v10.4 can ingest the ECMWF CLWC and CIWC either separately or as a sum ($q_c$=CLWC+CIWC). However, to save storage space, we recommend retrieving only the sum, $q_c$, from the ECMWF, as the relative fractions of ice and liquid water can be parameterized quite accurately using Eq. (4).

The column cloud water ($c_l$; kg m⁻²), which is needed for the in-cloud scavenging parameterization, is calculated by integrating $q_c$ over all vertical z levels:

$$c_l = \sum_z q_c(z) \rho_{\text{air}}(z) \Delta z, \quad (2)$$
where $\rho_{\text{air}}(z)$ is the density of the air in the grid cell, and $\Delta z$ is the vertical extent of the grid cell. In older FLEXPART versions, $c_l$ was parameterized based on an empirical equation given in Hertel et al. (1995) using the subgrid (see below for a description of how subgrid is defined) surface precipitation rate $I_s$ (mm h$^{-1}$). While such a parameterization is not needed anymore if $q_c$ is available, it is still activated in the case that cloud water input data are missing. However, in order to ensure that $c_l$ from the parameterization is statistically consistent with the cloud data, we derived the modified expression

$$c_l = 0.5 \times I_s^{0.36} \qquad (3)$$

using a regression analysis between existing cloud and precipitation data.

Precipitation is not uniform within a grid cell. To account for subgrid variability, it is assumed that precipitation is enhanced within a subgrid area and that no precipitation (and thus no scavenging) occurs outside this subgrid area. The subgrid area fraction and precipitation rate ($I_s$) are estimated from the grid-scale precipitation rate ($I_t$) based on values tabulated in Hertel et al. (1995). This parameterization of subgrid variability is used for all scavenging processes in FLEXPART and maintained from previous FLEXPART versions as described in Stohl et al. (2005).

The precipitation phase needed for the below-cloud scavenging scheme is simply based on ambient grid-scale temperature, with snow occurring below 0°C and rain above. For cloud water, $c_l$, we assume a temperature-dependent mixture of liquid and solid particles, for which the liquid fraction ($\alpha_L$) is calculated based on the local temperature $T$,

$$\alpha_L = \left( \frac{T - T_l}{T_L - T_l} \right)^2, \qquad (4)$$

where $T_L = 0$°C and $T_l = -20$°C. For $T > T_L$, $\alpha_L = 1$ and for $T < T_l$, $\alpha_L = 0$. Even when CLWC and CIWC are available as separate fields, we derive the liquid fraction ($\alpha_L$) of cloud water from the local temperature. The ice fraction $\alpha_I$ is $1 - \alpha_L$. Comparisons have shown that CLWC is very accurately reproduced by $\alpha_L q_c$.

The cloud information should be linearly interpolated like the other variables, and situations in which the diagnosed cloud is incompatible with the precipitation rate (be it because of interpolation or because of convective precipitation accompanied by too shallow or lacking grid-scale clouds) need to receive special treatment. This is planned for a version upgrade in the near future in conjunction with a better interpolation scheme for precipitation (see Hittmeir et al., 2018). In certain cases, the deposition calculation of FLEXPART might be improved by using higher-resolution precipitation data from other sources such as from radar observations (Arnold et al., 2015); however, as the precipitation and ECMWF cloud data, and also the precipitation and wind fields, may not match, this does not guarantee better results.

## 2.4.2 Below-cloud scavenging

For gases, the scavenging coefficient, $\Lambda$, for below-cloud scavenging is calculated as described in Asman (1995),

$$\Lambda = A I_s^B, \qquad (5)$$

where the scavenging parameters $A$ and $B$ depend on the chemical properties of the gas and are specified in the SPECIES\_nnn file as described in Sect. 5.1.3 (nnn represents the species number (0–999) used for the simulation). In older FLEXPART versions, this scheme was used also for aerosols; however, Grythe et al. (2017) developed a new aerosol scavenging scheme that is implemented in FLEXPART v10.4 and briefly summarized here.

The relevant processes of collision and attachment of ambient aerosol particles to falling precipitation depend mainly on the relationship between the aerosol and hydrometeor size and type (rain or snow) as well as to a lesser degree on the density and hygroscopicity of the aerosol. In FLEXPART v10.4, the dependence of scavenging on the sizes of both the aerosol and falling hydrometeors are taken into account by the schemes of Laakso et al. (2003) for rain and Kyrö et al. (2009) for snow. Both schemes follow the equation

$$\log_{10}\left(\frac{\Lambda}{\Lambda_0}\right) = C_* (a + b D_p^{-4} + c D_p^{-3} + d D_p^{-2} + e D_p^{-1} + f \left(\frac{I_s}{I_0}\right)^{0.5}), \qquad (6)$$

where $D_p = \log_{10} \frac{d_p}{d_{p0}}$, $d_p$ is the particle dry diameter provided in the SPECIES\_nnn file, $d_{p0} = 1$ m, $\Lambda_0 = 1$ s$^{-1}$ and $I_0 = 1$ mm h$^{-1}$. Coefficients for factors $a-f$ are different for rain and snow scavenging and are given in Table 2. The $C_*$ values are collection efficiencies that reflect the properties of the aerosol and must be given for both rain ($C_* = C_{\text{rain}}$) and snow scavenging ($C_* = C_{\text{snow}}$) in the SPECIES\_nnn file. Notice that by setting $C_{\text{snow}} = 0$, below-cloud scavenging by snowfall is switched off (similarly, $C_{\text{rain}} = 0$ for rain).

## 2.4.3 In-cloud scavenging

For in-cloud scavenging of both aerosols and gases, $\Lambda$ is calculated as described in Grythe et al. (2017):

$$\Lambda = i_{c_r} S_i I_s, \qquad (7)$$

where $i_{c_r} = 6.2$ is the cloud water replenishment factor, which was determined empirically in Grythe et al. (2017) (there it was named $i_{c_r}$), and $S_i$ is proportional to the in-cloud scavenging ratio, which is derived differently for gases and aerosols.

For gases, $S_i = \frac{1}{\frac{1-c_l}{HRT} + c_l}$, where $H$ is Henry's constant (describing the solubility of the gas and specified in the SPECIES\_nnn file), $R$ is the gas constant and $T$ is temperature. Notice that this is applied for both liquid-phase and ice clouds.
<table><caption>**Table 2.** Parameters used in Eq. (6) for below-cloud scavenging.</caption><thead><tr><th></th><th>C<sub>*</sub></th><th>a</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th><th>Reference</th></tr></thead><tbody><tr><td>Rain scavenging</td><td>C<sub>rain</sub></td><td>274.36</td><td>332 839.6</td><td>226656</td><td>58 005.9</td><td>6588.38</td><td>0.24498</td><td>Laakso et al. (2003)</td></tr><tr><td>Snow scavenging</td><td>C<sub>snow</sub></td><td>22.7</td><td>0</td><td>0</td><td>1321</td><td>381</td><td>0</td><td>Kyrö et al. (2009)</td></tr></tbody></table>

For aerosols, the in-cloud scavenging is dominated by activated particles forming cloud droplets or ice nuclei. Those may eventually combine to form a hydrometeor that falls out of the cloud, thus removing all aerosol particles contained in it. Therefore, $S_i$ depends on the nucleation efficiency ($F_{nuc}$) and $c_i$:

$$S_i = \frac{F_{nuc}}{c_i} \qquad (8)$$

$F_{nuc}$ describes how efficiently the aerosol particles are activated as cloud droplet condensation nuclei (CCN) or ice nuclei (IN):

$$F_{nuc} = \alpha_L \text{CCN}_{\text{eff}} + \alpha_I \text{IN}_{\text{eff}}, \qquad (9)$$

where the relative abundances of the liquid and ice phase are accounted for by the factor $\alpha_L$. Values for the efficiencies, CCN<sub>eff</sub> and IN<sub>eff</sub>, are available from the literature for many different types of aerosols (e.g., black carbon, mineral dust particles or soluble particles) and some have been collected in SPECIES_nnn files distributed with FLEXPART (see Sect. 5.1.3). The CCN<sub>eff</sub> and IN<sub>eff</sub> values are set for an aerodynamic particle radius of 1 µm, but CCN and IN efficiencies increase with increasing particle size. The in-cloud parameterization takes this into account. For further details on the wet scavenging scheme used in FLEXPART, see Grythe et al. (2017).

### 2.4.4 Influence of wet scavenging on the aerosol lifetime

Aerosol wet scavenging controls the lifetime of most aerosols. In Fig. 2, we compare modeled e-folding lifetimes from a number of FLEXPART simulations using different model versions and switching off in-cloud and below-cloud scavenging in FLEXPART v10.4 with measured lifetimes. The parameter settings in FLEXPART used for these comparisons were the same as used by Grythe et al. (2017). To derive aerosol lifetimes in a consistent way from both measurements and model simulations, a radionuclide attached to ambient aerosol and a noble gas radionuclide were used. Kristiansen et al. (2016) used the same method to compare many different aerosol models, and we refer to their paper for more details on the method. For our model simulations, several size bins of aerosols were used, though total concentrations and lifetimes are largely controlled by 0.4 and 0.6 µm particles (see Grythe et al., 2017). E-folding lifetimes increase from 5.8 to 10.0 d between FLEXPART v9

and v10.4. A simulation performed with v10.4 but which emulated the in-cloud scavenging of v9 showed that the difference is mainly due to the decreased in-cloud scavenging in the new removal scheme compared to the old one. Notice that the lifetime obtained with v10.4 is much closer to the observation-based lifetimes. Turning off the below-cloud removal has a relatively small effect, increasing the lifetime to 11 d, whereas turning off the in-cloud removal extends the lifetime to the unrealistic value of 66 d (see bottom two panels in Fig. 2). This highlights the dominant role of in-cloud removal for accumulation-mode particles in FLEXPART.

Notice that compared to older versions of FLEXPART, the SPECIES_nnn files now include additional parameters related to the wet deposition scheme. Old input files, therefore, need to be updated for use with FLEXPART v10.4. The required format changes are detailed in Sect. 5.1.3.

## 2.5 Source-receptor matrix calculation of deposited mass backward in time

When running FLEXPART forward in time for a depositing species with a given emission flux (kilograms per release as specified in file RELEASES), the accumulated wet and dry deposition fluxes (ng m⁻²) are appended to the FLEXPART output files (grid_conc_date and/or grid_pptv_date, for which date represents the date and time in format YYYYMMDDhhmmss; see Sect. 6) containing the atmospheric concentration and/or volume mixing ratio output. The deposition is always given in mass units, even if atmospheric values are given in mixing ratio units. In contrast to concentration values, deposition quantities are accumulated over the time of the simulation, so the deposited quantities generally increase during a simulation (except when radioactive decay is activated, which also affects deposited quantities and can decrease them).

As discussed in Sect. 1, running FLEXPART backward in time for calculating SRRs is more efficient than running it forward if the number of (potential) sources is larger than the number of receptors. For atmospheric concentrations (or mixing ratios), the backward mode has been available from the very beginning and in an improved form since FLEXPART v5 (Stohl et al., 2003; Seibert and Frank, 2004). This has proven very useful for the interpretation of ground-based, shipborne or airborne observations (e.g., to characterize sources contributing to pollution plumes). Furthermore, the inversion scheme FLEXINVERT (Thompson and Stohl, 2014) that is used to determine the fluxes of greenhouse gases
<figure><img src="image_4.png" /></figure>

**Figure 2.** Aerosol lifetimes estimated from the decrease in radionuclide ratios (aerosol-bound <sup>137</sup>Cs and noble gas <sup>133</sup>Xe as a passive tracer) with time after the Fukushima nuclear accident, as measured and modeled at a number of global measurement stations. For details on the method, see Kristiansen et al. (2016). *E*-folding lifetimes, τ<sub>e</sub>, are estimated based on fits to the data and reported in each panel. In the top panel, the observed values are shown and in subsequent panels from the top, modeled values are given for FLEXPART v9, FLEXPART v10.4, FLEXPART v10.4 with parameter settings to emulate removal as in v9, FLEXPART v10.4 with no below-cloud removal and FLEXPART v10.4 with no in-cloud removal.

is based on backward simulations. However, there are also
measurements of deposition on the ground, e.g., in precipita-
tion samples or ice cores, and for this type of measurement
no backward simulations were possible until recently. There-
fore, Eckhardt et al. (2017) introduced the option to calculate
SRR values in backward mode also for wet and dry deposi-
tion, and a first application to ice core data was presented by
McConnell et al. (2018). It is anticipated that quantitative in-
terpretation of ice core data will be a major application of the
new backward mode, which is efficient enough to allow for

the calculation of, for example, 100 years of seasonally re-
solved deposition data in less than 24 h of computation time.

We illustrate the different backward modes and explain
the required settings with an example. The calculations were
run for a single receptor location, Ny-Ålesund in Spitsber-
gen (78.93° N, 11.92° E) and for the 24 h period from 18 Au-
gust 2013 at 20:00 UTC to 19 August 2013 at 20:00 UTC.
SRR values are calculated for the atmospheric concentration
averaged over the layer 0–100 m a.g.l., as well as for wet
and dry deposition. The substance transported is black car-
bon (BC), which is subject to both dry and wet deposition.
Backward simulations for wet and dry deposition must al-
ways be run separately. In order to obtain SRR values for
total deposition, results for wet and dry deposition need to be
summed.

The backward mode is activated by setting the simulation
direction, LDIRECT in file COMMAND (see Sect. 5), to -1.
The three simulations are obtained by consecutively setting
IND_RECEPTOR to 1 (for concentration), 3 (wet deposi-
tion) and 4 (dry deposition). IND_SOURCE is always set to
1, meaning that the sensitivities (SRR values) are calculated
with respect to physical emissions in mass units. A complete
list of possible options is reported in Table 1 of Eckhardt et
al. (2017).

Figure 3 shows the resulting SRR (i.e., emission sensitiv-
ity) fields for the concentration as well as dry and wet de-
position at the receptor. Dry deposition occurs on the Earth’s
surface, and therefore particles are released in a shallow layer
adjacent to the surface. Its height is consistent with the shal-
low depth over which dry deposition is calculated in forward
mode (user settings for the release height are ignored for dry
deposition backward calculations). Dry deposition rates are
the product of the surface concentration and the deposition
velocity. Therefore, the SRR fields for surface concentration
(Fig. 3a) and dry deposition (Fig. 3b) show similar patterns,
in this case indicating high sensitivity for sources over Scan-
dinavia and northwestern Russia. The differences in the spa-
tial patterns are mainly due to temporal variability in the dry
deposition velocity at the receptor caused by varying mete-
orological conditions (e.g., stability) and surface conditions
during the 24 h release interval.

Wet deposition, on the other hand, can occur anywhere in
the atmospheric column from the surface to the top of the
precipitating cloud. FLEXPART automatically releases par-
ticles in the whole atmospheric column (again, user settings
for the release height are ignored), but particles for which no
scavenging occurs (e.g., those above the cloud top or when
no precipitation occurs) are immediately terminated. There-
fore, and because of the vertical variability of the scaveng-
ing process, the sensitivity for the deposited mass can de-
viate significantly from the sensitivity corresponding to sur-
face concentration. Here (Fig. 3c), the sensitivity is high over
Scandinavia and northwestern Russia, as was already seen
for surface air concentrations and dry deposition. However,
in addition, sources located in North America and eastern
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Source–receptor relationships (for emissions occurring in the lowest 100 m a.g.l.) for black carbon observed at Ny-Ålesund in Svalbard for a 24 h period starting on 18 August 2013 at 20:00 UTC. The sensitivities were calculated for (a) concentrations (s) in the layer 0–100 m a.g.l., (b) dry deposition (mm) and (c) wet deposition (mm).</figcaption></figure>

Siberia also contribute strongly to wet deposition. The maximum over the ocean close to the North American east coast is likely due to lifting in a warm conveyor belt, followed by fast transport at high altitude.

Concentration, dry deposition and wet deposition at the receptor can be calculated from the SRR fields shown in Fig. 3 as follows.

$$ 
\begin{aligned}
c &= m_c \cdot q \\
d_d &= m_d \cdot q \\
d_w &= m_w \cdot q
\end{aligned} 
\qquad (10) $$

Here, c is the modeled concentration (in kg m⁻³), $d_d$ the dry deposition rate and $d_w$ the wet deposition rate (both in kg m⁻² s⁻¹). In this specific case with only a single scalar receptor, the source–receptor matrix degenerates to a vector of the SRR values, one for each of the three types of receptor ($m_c$ for concentration in units of seconds, $m_d$ for dry deposition and $m_w$ for wet deposition, both in units of meters). In order to obtain the concentration or the deposition rates, these vectors need to be multiplied with the vector of emissions $q$ (kg m⁻³ s⁻¹). If the total deposition is desired, the deposition rates $d_d$ and $d_w$ can be multiplied with the receptor time interval $\Delta T_r$, in our case 86 400 s (24 h). Note that this is the period during which particles are released according to

the specification of the RELEASES file. The emission fluxes must be volume averages over the output grid cells specified in the OUTGRID file, typically surface emission fluxes (in kg m⁻² s⁻¹) divided by the height of the lowermost model layer.

## 2.6 Sensitivity to initial conditions

Backward simulations with FLEXPART in the context of inverse modeling problems typically track particles for several days up to a few weeks. This is sufficient to estimate concentrations at the receptor only for species with atmospheric lifetimes shorter than this period. Many important species (e.g., greenhouse gases such as methane) have considerably longer lifetimes. For such long-lived species, most of the atmospheric concentration variability is still caused by emission and loss processes occurring within the last few days before a measurement because the impact of processes occurring at earlier times is smoothed out by atmospheric mixing. This leads to a relatively smooth “background” (in time series analyses sometimes also called a baseline) that is often a dominant fraction of the total concentration but that does not vary much with time, with short-term fluctuations on top of it. The signal of the regional emissions around the measurement site is mostly contained in the short-term concentration fluctuations but in order to use it in inverse modeling, the background still needs to be accounted for, as otherwise no direct comparison to measurements is possible.

One simple method is to estimate the background from the measurements as, e.g., in Stohl et al. (2009). A better approach is to use a concentration field taken from a long-term forward simulation with an Eulerian model or with FLEXPART itself, especially if nudged to observations (Groot Zwaaftink et al., 2018), as an initial condition for the backward simulation. This field needs to be interfaced with the FLEXPART backward simulation by calculating the receptor sensitivity to the initial conditions (see Eqs. 2–6 in Seibert and Frank, 2004). For instance, for a 10 d backward simulation, the concentration field needs to be sampled at those points in time and space when and where each particle trajectory terminates 10 d back in time. Furthermore, it is necessary to quantify the effects of deposition or chemical loss during the backward simulation on this background (the factor $p(0)$ in Seibert and Frank, 2004). For example, chemical reactions with hydroxyl radicals will reduce initial concentrations of methane en route to the receptor, even though this is not much during a 10 d period.

Since version 8.2, FLEXPART has provided an option to quantify the influence of initial conditions on the receptor in backward simulations, which is activated with the switch `LINIT_COND` in file `COMMAND`. Then, gridded fields containing the sensitivities to background mixing ratios (or concentrations, depending on user settings for the switch `LINIT_COND` in file `COMMAND`) are produced and stored in the output files `grid_initial_nnn` (nnn stands for the
species number) on the same 3-D grid as the regular output, defined in the files OUTGRID and OUTGRID_NEST. In this case, a concentration would be calculated as

$$c = m_i \cdot c_b + m_c \cdot q, \qquad (11)$$

where $m_i$ denotes the sensitivity to the initial condition and $c_b$ the background concentration when and where particles are terminated.

Figure 4 shows an example of the use of the sensitivities of receptor mixing ratios (here, of methane) to both surface emissions and initial conditions. The panel (b) shows the sensitivity to surface emissions on one particular day, and the panels (c) and (d) show the sensitivity to initial conditions below and above 3000 m for the same day. Both results are from an 8 d backward simulation from one receptor site in Sweden. It can be seen that the sensitivity to emissions is highest close to the station, but there is also substantial sensitivity to emission uptake over large parts of central and eastern Europe. The particles terminated 8 d before arrival at the receptor in a roughly croissant-shaped area covering large parts of Europe and the North Atlantic, as indicated by the sensitivity to initial conditions. Most of the sensitivity is located below 3000 m but there is also some influence from higher levels. Notice that only two layers are shown in Fig. 4, whereas the real model output has much higher vertical resolution.

The sensitivity to initial conditions was interfaced with a domain-filling methane forward simulation as described in Groot Zwaafink et al. (2018) (not shown), while the emission sensitivity was interfaced with an emission inventory for methane (not shown), as given by Eq. (11). This was done for daily simulations throughout 1 month, thus generating a time series of background mixing ratios (from the first term in Eq. 11 only) and total mixing ratios (Fig. 4a). The latter include the contributions from emissions during the 8 d backward simulation. It can be seen that the methane background advected from 8 d back varies relatively little between about 1910 and 1940 ppbv, while the emission contributions vary from 0 (on 29 October) to about 200 ppbv (on 19 October, the date for which the sensitivity plots are shown).

In practical applications for inverse modeling, source-receptor sensitivities are often only needed at the surface (as most emissions occur there), while sensitivities to the background are needed in 3-D. By setting the option SURF_ONLY to 1 in the COMMAND file, the regular output files grid_time_date_nnn containing the source-receptor sensitivities will include only the first vertical level as defined in the OUTGRID, while the full vertical resolution is retained in grid_initial_nnn files containing the sensitivities to the initial conditions. Since the data amounts stored in the grid_time_date_nnn files can be much larger than in the grid_initial_nnn files, this is a highly efficient way to save storage space. This setup also interfaces directly with the inverse modeling package FLEXINVERT (Thompson and Stohl, 2014). An application can be found in Thompson et al. (2017)

wherein initial conditions were taken from a gridded observation product. A further output option, which was also introduced for practical considerations of inverse modeling, is the LINVERSIONOUT option in the file COMMAND. If LINVERSIONOUT is set to 1, then the grid_time_date_nnn and grid_initial_nnn files are written per release with a time dimension of footprints instead of the default per footprint with a time dimension of releases. Since inverse modeling assimilates atmospheric observations and each observation is represented by a single release, it is computationally more efficient to read in the grid files separated by release. This output format also interfaces directly with FLEXINVERT.

## 2.7 Chemical reactions with the hydroxyl radical (OH)

The hydroxyl (OH) radical reacts with many gases and is the main cleansing agent in the atmosphere. While it is involved in highly nonlinear atmospheric chemistry, for many substances (e.g., methane) a simplified linear treatment of loss by OH is possible using prescribed OH fields. For this, monthly averaged 3° × 5° resolution OH fields for 17 atmospheric layers are used in FLEXPART. The fields were obtained from simulations with the GEOS-Chem model (Bey et al., 2001) and are read from the file OH_variables.bin by the subroutine readOHfield.f90.

Tracer mass is lost by reaction with OH if a positive value for the OH reaction rate is given in the file SPECIES_nnn. In FLEXPART v10.4, the OH reaction scheme was modified to account for (i) hourly variations in OH and (ii) the temperature dependence of the OH reaction rate (Thompson et al., 2015). This makes the chemical loss calculations more accurate, especially for substances with shorter lifetimes (of the order of weeks to months), for example ethane. Hourly OH fields are calculated from the stored monthly fields by correcting them with the photolysis rate of ozone calculated with a simple parameterization for cloud-free conditions based on the solar zenith angle (gethourlyOH.f90):

$$\mathrm{OH} = \frac{j}{j^*} \mathrm{OH}^*, \qquad (12)$$

where j represents the hourly photolysis rates calculated for all 3-D locations in the field, while $j^*$ represents the corresponding monthly mean rates, precalculated and stored in file OH_variables.bin together with the monthly mean fields $OH^*$ (see Sect. 5.1.8). The motivation for this is that OH production closely follows the production of O(1D) by the photolysis of ozone, allowing for this simple parameterization of OH variability. At any time, two hourly OH fields are in memory and are interpolated to the current time step. Figure 5 shows the annual and daily variation of OH for two locations as obtained with this simple parameterization.

The OH reaction rate $\kappa$ ($s^{-1}$) is calculated in ohreaction.f90 using the temperature-dependent
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Example of FLEXPART 8 d backward runs for methane from a site in southern Sweden (Hyltemossa) demonstrating the combined use of sensitivities to emissions and initial conditions. (a) Time series of methane background mixing ratios and total mixing ratios in October 2016. (b) Sensitivity of the methane mixing ratio at Hyltemossa on 19 October 2016 to methane emissions at the surface. (c) Sensitivity of the methane mixing ratio at Hyltemossa on 19 October 2016 to methane initial conditions below 3000 m. (d) Sensitivity of the methane mixing ratio at Hyltemossa on 19 October 2016 to methane initial conditions above 3000 m. Blue asterisks on the maps mark the receptor location.</figcaption></figure>

formulation

$$ \kappa = C T^N e^{-D/T} [\text{OH}], \qquad (13) $$

where C, N and D are species-specific constants (assigned in the SPECIES_nnn files), T is the absolute temperature, and [OH] the OH concentration (Atkinson, 1997). As the OH concentration in file OH_variables.bin is given in units of molecules per cubic centimeter, the unit of C needs to be in cubic centimeters per molecule per second (cm³ molecule⁻¹ s⁻¹). The mass m of a given species after reaction with OH is determined as

$$ m(t + \Delta t') = m(t) e^{-\kappa \Delta t'}, \qquad (14) $$

where Δt' is the reaction time step (given by lsynctime).

Backwards compatibility with the former temperature-independent specification of the OH reaction (version 9 and before) can be achieved by setting the constant N in the SPECIES_nnn file to zero. The constants C and D can be derived from the former parameters as follows:

$$ C = \kappa_r e^{D/T_r} \qquad (15) $$

and

$$ D = A/R, \qquad (16) $$

where A is the activation energy, R is the gas constant and κr is the former OH reaction rate (referring to T_r = 298 K), which were specified in the SPECIES_nnn file for earlier versions.

OH fields other than those provided with the model code have been tested in FLEXPART. These fields may have higher spatial and temporal resolution (e.g., Fang et al., 2016), which is important for chemical species with short lifetimes. Users are required to modify readOHfield.f90 and gethourlyOH.f90 to read in other OH fields and be aware that expressions of the OH reaction rate or reaction with OH might differ from those in the above equations. If this is the case users need to modify ohreaction.f90, too.

## 2.8 Dust mobilization scheme

Desert dust is a key natural aerosol with relevance for both climate and air quality. FLEXPART has been used earlier with preprocessors to initialize dust amounts from wind speed and surface properties following Tegen and Fung (1994) (Sodemann et al., 2015). Now a dust mobilization routine has been included as a preprocessing tool in FLEXPART v10.4. The scheme, called FLEXDUST, was devel-
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Annual and daily OH concentration variation as obtained with the simple parameterization based on photolysis rates of ozone for two locations, one in the Northern Hemisphere (a) and one in the Southern Hemisphere (b). Line labels correspond to the time of day.</figcaption></figure>

oped to simulate mineral dust transport with FLEXPART in forward or backward simulations (Groot Zwaastink et al., 2016). This module runs independently from FLEXPART and produces gridded output of mineral dust emissions as well as input files (RELEASES) that can be used for FLEXPART simulations of atmospheric transport. It can thus be considered a preprocessing (for forward simulations) or post-processing tool (for backward simulations) for FLEXPART v10.4.

In FLEXDUST, emission rates are estimated according to the emission scheme proposed by Marticorena and Bergametti (1995). We thereby assume that sandblasting occurs in the case that sand is present and a minimum threshold based on the size-dependent threshold friction velocity following Shao and Lu (2000) can be applied. The following are used as input for the model: ECMWF operational analysis or

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Average annual dust emission for the period 1990–2012 estimated with FLEXDUST driven with ERA-Interim meteorology.</figcaption></figure>

ERA-Interim reanalysis data, Global Land Cover by National Mapping Organizations version 2 (Tateishi et al., 2014), and sand and clay fractions from the Global Soil Data Task (2014). Erodibility is enhanced in topographic depressions, and dust emissions are modified by soil moisture and snow cover. The module includes high-latitude dust sources in the Northern Hemisphere. These sources are rarely included in global dust models, even though they appear important for the climate system and substantially contribute to dust in the Arctic (Bullard et al., 2016; Groot Zwaastink et al., 2016). Icelandic deserts in particular are known to be highly active, and a high-resolution surface type map for Iceland can therefore be included in FLEXDUST simulations (Arnalds et al., 2016; Groot Zwaastink et al., 2017). Like in FLEXPART, nested meteorological fields can be used for specific regions of interest. The size distribution of emitted dust follows Kok (2011), is independent of friction velocity and is by default represented by 10 size bins. This can be changed depending on known properties or assumptions of dust sources. The dust particles are assumed to be spherical in FLEXPART. An example of annual mean dust emissions from 1990 to 2012 calculated with FLEXDUST driven with ERA-Interim meteorology is shown in Fig. 6. Further details on FLEXDUST, including model evaluation, are given by Groot Zwaastink et al. (2016). The source code is available from the git repository: https://gitnilu.no/christine/flexdust.git (last access: 30 October 2019).

## 3 Parallelization

In a Lagrangian model like FLEXPART, particles move totally independently of each other. This facilitates efficient parallelization of the code. The most simple and often most effective way is running several instances of the model in parallel. For example, if the model is to be run backwards (for 10 d, for example) at regular intervals from a measurement site for a year, one could run the model separately, in par-
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Computational time (a, b) and speedup (c, d) for up to 16 processes on a single node. In panels (a, c), all processes read meteorological input data, whereas in panels (b, d), a dedicated process reads and distributes input data for $N_p \ge 4$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Computational time (a) and speedup (b) for up to 256 processes on 16 nodes. Logarithmic scaling along both axes. For $n \ge 4$ a dedicated process reads and distributes input data.</figcaption></figure>

allel, for monthly subperiods. The total computation time of the 12 monthly processes together is nearly the same as if the model is run as one process for the whole year. Some overhead in processing input data occurs because, in the above example, 10 extra days of data per process are needed to calculate trajectories 10 d back into the preceding month. One disadvantage of that approach is that the memory needed for holding the meteorological input data and the model output fields is multiplied. However, this overhead is often small; thus, this approach has been used very often by FLEXPART users in the past.

Even if a task cannot easily be decomposed into runs for different periods or sources, trivial parallelization is still possible if a large number of particles is desired, for example in a domain-filling simulation for which tens of millions of particles may be used. The strategy in this case would be to assign a fraction of the particles to each run. Note that different random seeds should be used for each run, which requires a manual change and recompilation of the code.

As a user-friendly alternative, FLEXPART v10.4 has been parallelized using standard parallelization libraries. Common parallelization libraries are Open Multi-Processing (OpenMP; http://www.openmp.org/, last access: 30 Octo-
ber 2019), which is designed for multicore processors with
shared memory, and Message Passing Interface (MPI, 2015)
for distributed memory environments. Examples of other
Lagrangian particle models that have been parallelized are
NAME (Jones et al., 2007) and FLEXPART-WRF (Brioude
et al., 2013), which use a hybrid approach (OpenMP+MPI).
For FLEXPART v10.3 we decided to use a pure MPI ap-
proach for the following reasons.

- It is simpler to program than a hybrid model and more flexible than a pure OpenMP model.

- While OpenMP in principle may be more effective in a shared memory environment, MPI can often perform equally well or better provided there is not excessive communication between the processes.

- MPI offers good scalability and potentially low overhead when running with many processes.

3.1 Implementation

The FLEXPART code contains several computational loops
over all the particles in the simulation, which is where
most of the computational time is spent for simulations with
many particles. The basic concept behind our parallel code
closely resembles the “trivial parallelization” concept de-
scribed above. When launched with a number of processes,
Np, each process will separately calculate how many parti-
cles to release per location, attempting to achieve an approx-
imately even distribution of particles among the processes
while keeping the total number of particles the same as for a
simulation with the serial version. Each running process will
generate an independent series of random numbers and sepa-
rately calculate trajectories and output data for its set of par-
ticles. Explicit communications between processes are only
used when the output fields are combined at the master pro-
cess (MPI rank 0) using MPI_Reduce operations, before
writing the output. Also, in the case in which the output of
all individual particle properties is desired (option IPOUT1
= 1 or 2 in file COMMAND), we let each process append its
data to the same file. We thus avoid the costly operation of
transferring particle properties between processes. The per-
formance of the implementation is discussed in Sect. 3.2 (see
Fig. 7).

Some parts of the code are not simply loops over parti-
cles, most notably the routines for reading and transforming
the input meteorological data. It follows that the performance
gain of using parallel FLEXPART in general is better for sim-
ulations with a larger number of particles. We have, however,
implemented a feature whereby instead of having each MPI
process read and process the same input data, one dedicated
MPI process is set aside for this purpose. When the simula-
tion time *t* lies in the interval between wind field time *T*<sub>*i*</sub> and
*T*<sub>*i*+1</sub>, all other processes calculate particle trajectories, while
this dedicated process ingests input fields from time *T*<sub>*i*+2</sub>. At

simulation time *t* = *T*<sub>*i*+1</sub> the dedicated “reader process” will
distribute the newest data to the other processes and imme-
diately start reading fields for time *T*<sub>*i*+3</sub>, while the other pro-
cesses continue doing trajectory calculations. A hard-coded
integer (*read_grp_min* in file *mpi_mod.f90*) is used to
set the minimum number of total MPI processes for which
this separate process will be reserved for reading input data.
For the examples shown in Sect. 3.2 a value of 4 was used
(Figs. 7 and 8).

3.2 Performance aspects

To assess the performance of the parallel code we performed
three scaling experiments of various size on different com-
putational platforms.

3.2.1 40 million particles, single 16-core node

In the following we present the results from running the code
on a machine equipped with an Opteron 6174 processor with
16 cores. Compilation was done using gfortran version 4.9.1
and OpenMPI version 1.8.3. For the experiment, 40 million
particles were released and propagated 48 h forward in time.
We ran with this setup with an increasing number of pro-
cesses, from 1 to 16. All time measurements in the code were
made with the MPI_wtime() subroutine.

For the first experiment, every process separately pro-
cessed the meteorological input data. Figure 7a and c show
the CPU time Tn used in the case of n processes and the rela-
tive speedup factor S(n) = T1/Tn. Time and speedup shown
for “particle loops” includes the three most computationally
demanding particle loops (integration of the Langevin equa-
tion, wet deposition and concentration calculations), but, in
addition, FLEXPART contains a few smaller loops over par-
ticles that exhibit similar performance improvements. We see
that for 40 million particles, the loops over particles take the
largest share, at least 87% of the total time when run with
one process. Close-to-perfect speedup is expected and ob-
served for these loops (compare results for “particle loops”
and “ideal (particle loops)” in Fig. 7a and c). The major bot-
tleneck for overall performance in this case is that each pro-
cess reads the same input files from disk, thus forcing the
others to wait. This bottleneck causes the speedup to deviate
substantially from the ideal situation when more than a few
processes are used (compare results for “total” and “ideal” in
Fig. 7a and c).

Next we repeated the experiment above but set aside a ded-
icated process for reading the meteorological data whenever
*n* ≥ 4. The results are shown in Fig. 7b and d. Numerical val-
ues for the speedup factors for selected numbers of processes
are given in Table 3. We observed that with *n* ≥ 7 there was
consistently a benefit to setting aside the dedicated reader
process, whereas for *n* < 7 it was more effective to have all
processes read data and thus an extra process available for
doing the trajectory calculations. These results will of course
<table><caption>**Table 3.** Computational speedup *S* for up to 16 processes (single-node experiment) for the two different MPI modes, with 40 million released particles.</caption><thead><tr><th>Number of processes</th><th>1</th><th>2</th><th>4</th><th>8</th><th>16</th></tr></thead><tbody><tr><td>All processes read</td><td>1</td><td>1.89</td><td>3.30</td><td>5.76</td><td>9.16</td></tr><tr><td>Dedicated read process</td><td>1</td><td>1.91</td><td>3.09</td><td>6.72</td><td>13.10</td></tr></tbody></table>

vary with the resolution of the input data, the number of particles and the system on which the program is run.

### 3.2.2 500 million particles, multiple 16-core nodes

We performed a larger-scale experiment at the Abel computer cluster<fnref n="1" /> using up to 256 cores on 16 nodes with Intel Xeon E5-2670 CPUs. For each node, up to 16 cores were used, and then the number of nodes was determined by the total number of processes launched. The FLEXPART setup was similar to the previous single-node experiment, but we increased the number of particles to 500 million and reduced the simulated time to 12 h. Compilation was done with Intel Fortran v16.0.1 and OpenMPI v1.10.2.

Run time and speedup factors are shown in Table 4 and Fig. 8. As before we see essentially perfect speedup of the computationally intensive parts (the particle loops), which is expected. Table 4 also gives the parallel efficiency, which is seen decreasing for larger $N_p$. This is partly due to the increased cost of MPI communications and also because the nonparallel parts of the code have relatively higher impact. With 256 processes there are only about 2 million particles per process and the CPU time is not as clearly dominated by the particle loops as when 500 million particles all run in one process. In addition, the initialization of the code (allocation of arrays, reading configuration files) takes around 20 s for this run, which is significant given a total run time of 299 s. Thus, parallel efficiency would increase for longer simulation times and/or for simulations with more particles per process, i.e., realistic cases that are more likely to be run with such a large number of processes.

### 3.2.3 900 000 particles, laptop and single 16-core node

Finally, we examined a small-scale experiment in which we released 900 000 particles and simulated 15 d of transport. The performance was tested on two systems; a ThinkPad P52s laptop (Intel i7-8550U CPU with four cores; results in Table 5) and a machine equipped with an AMD Opteron 6386 SE processor (16 cores; results in Table 6). With this relatively lower number of particles it is not surprising to see that the parallel efficiency is lower than in the preceding ex-

amples. Still, we see that a speedup of 2.38 on a 4-core laptop and 5.25 on a 16-core machine is attainable. We also note that for practical applications, users would likely use the serial version for applications with so few particles and, if there are many such runs to be done, use trivial parallelization by submitting many separate serial runs in parallel. The parallelization feature is most useful for cases with a very large number of particles that cannot so easily be split in many separate runs, such as domain-filling simulations.

## 3.3 Validation

In order to ensure that the parallel version produces results with the same accuracy as the serial version, we have performed a set of tests and validation experiments. A direct comparison between the versions can only be performed in statistical terms because FLEXPART uses Gaussian-distributed random numbers for calculating the turbulent velocities of the particles. For the parallel version we let each process independently calculate a set of random numbers, which leads to small numeric differences (arising from the random “noise”) between the versions.

To confirm that the only source of differences between the serial and parallel code is in the random number generation, we first observe that when the parallel executable is run using only one process, it produces results identical to the serial version. This is as expected, as the first MPI process (rank 0) always uses the same random number seeds as the serial version.

Next, we have done tests in which all random numbers are set to zero in both codes, corresponding to switching off the turbulent displacements, and we run the parallel version using multiple processes. The outputs from the serial and parallel versions of the code when run this way are identical except for small differences due to round-off errors (e.g., in concentration calculations – these round-off errors are typically larger in the serial version due to the larger number of particles).

# 4 Installation, compilation and execution

FLEXPART is usually used in a Linux environment, which we also assume for the following instructions. However, the model has also been implemented successfully under MacOS and MS Windows. The default Fortran compiler for FLEXPART v10.4 is gfortran, but ifort, Absoft and PGI compilers have been used as well.

## 4.1 Required libraries and FLEXPART download

As the meteorological data from numerical weather prediction models are usually distributed in GRIB format, a library for reading GRIB data is required. It is recommended to use `ecCodes` (https://software.ecmwf.int/wiki/display/ECC, last access: 30 October 2019), the primary GRIB

<footnote n="1">¹Owned by the University of Oslo and Uninett/Sigma2, operated by the Department for Research Computing at USIT, the University of Oslo IT department; https://www.uio.no/english/services/it/research/hpc/abel/ (last access: 30 October 2019).</footnote>
<table><caption>**Table 4.** Run time and speedup for the multi-node experiment with 500 million particles. Up to 16 nodes in the Abel cluster (University of Oslo).</caption><thead><tr><th>Number of processes <i>n</i></th><th>1</th><th>2</th><th>4</th><th>8</th><th>16</th><th>32</th><th>64</th><th>128</th><th>256</th></tr></thead><tbody><tr><th>Total run time (s)</th><td>39 536</td><td>19 681</td><td>14 123</td><td>6380</td><td>3061</td><td>1568</td><td>828</td><td>491</td><td>299</td></tr><tr><th>Speedup factor <i>S</i></th><td>1</td><td>2.01*</td><td>3.37</td><td>6.20</td><td>12.92</td><td>25.22</td><td>47.76</td><td>80.53</td><td>132.12</td></tr><tr><th>Parallel efficiency (<i>S</i>/<i>n</i>)</th><td>1</td><td>1.004</td><td>0.843</td><td>0.775</td><td>0.807</td><td>0.788</td><td>0.746</td><td>0.629</td><td>0.516</td></tr></tbody></table>

<table><caption>**Table 5.** Run time and speedup using up to four cores on a ThinkPad P52s laptop (900 000 particles).</caption><thead><tr><th>Number of processes <i>n</i></th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><th>Total run time (s)</th><td>3504</td><td>2035</td><td>1790</td><td>1470</td></tr><tr><th>Speedup factor <i>S</i></th><td>1</td><td>1.72</td><td>1.97</td><td>2.38</td></tr><tr><th>Parallel efficiency (<i>S</i>/<i>n</i>)</th><td>1</td><td>0.86</td><td>0.66</td><td>0.775</td></tr></tbody></table>

encoding-decoding package used at the ECMWF (recent enough versions of its predecessor `grib_api`, no longer supported after 2018, can also be used). Data in GRIB-2 format can be compressed. If this is the case for the input data, the `jasper` library is needed<fnref n="2" />. If it is desired to produce FLEXPART output in the NetCDF format, the NetCDF Fortran Library (https://www.unidata.ucar.edu/software/netcdf/, last access: 30 October 2019) is also required.

In order to obtain the FLEXPART source code, download the appropriate v10.4 tarball from the FLEXPART website<fnref n="3" /> and unpack it.

```
tar -xvf flexpart10.4.tar
```

To obtain the latest available model version, clone the FLEXPART git repository from the FLEXPART community website.

```
git clone https://www.flexpart.eu/gitmob/flexpart
```

This repository mirrors https://transport nilu.no (https://git Nilu.no/flexpart/flexpart/-/releases, last access: 18 November 2019). Additional mirrors exist, e.g., at GitHub (https://github.com/flexpart/flexpart, last access: 30 October 2019) and BitBucket (https://bitbucket.org/flexpart/flexpart, last access: 30 October 2019). See the “Code and data availability” section for additional information. After unpacking the tarball or cloning the repository, a local directory structure as shown in Table 1 is created.

The directory `src` contains the code and a makefile. The makefile needs to be adapted to the compiler and libraries present on the local system. Appendix A4 describes these steps in detail, including manual installation of the libraries. This was tested for Ubuntu 16.04.3 as well as 18.04.3 LTS Linux and MacOS (OS versions 14.5.0 and 18.6.0). Both a serial and a parallel executable can be built from the FLEXPART v10.4 source files.

## 4.2 Compiling and running the serial version

After correctly setting the library paths in the makefile, the command `make` produces the executable called `FLEXPART`. It can be executed from the command line by `./FLEXPART` and then expects a file `pathnames` to exist in the current working directory. This file contains information on where input data are located and where output data shall be produced (see Sect. 5). Note that `pathnames` is expected in the directory from which FLEXPART is started, which can be different from where the executable file is located. A different name of a `pathnames` file can be also given as an argument. FLEXPART can thus be invoked according to the following generic syntax:

```
path_to_flexpart/flexpart_executable
path_to_pathnames/pathnames_file.
```

Using an optional argument, `-v` (verbose mode), will display additional information during the run. Even more information, including clock time between different program units, will be printed with `-v2`. Invoking FLEXPART with the flags `-i` and `-i2` (info mode) will provide detailed run-specific information while reading input files. However, in this mode FLEXPART then stops before particle trajectories are calculated.

## 4.3 Compiling and running the parallel version

Most subroutines calling MPI functions are in a single module named `mpi_mod.f90`. Other FLEXPART source files that depend on this module are given the `_mpi.f90` suffix to distinguish them from the serial version. During compilation the makefile selects the source files automatically depending on whether the parallel or serial version is built.

In order to compile and run the parallel version, an MPI library must be installed, either as a package from

<footnote>* Superlinear speedup (efficiency greater than 100 %) as seen here is usually attributed to memory and/or cache effects.</footnote>

<footnote n="2">²The `jasper` package is available in Linux distributions; https://github.com/mdadams/jasper (last access: 30 October 2019).</footnote>

<footnote n="3">³The website, https://flexpart.eu/ (last access: 30 October 2019), provides additional information that can be used to supplement these instructions.</footnote>
<table><caption>**Table 6.** Run time and speedup using up to 16 cores on a machine equipped with an AMD Opteron 6386 SE processor (900 000 particles).</caption><thead><tr><th>Number of processes <i>n</i></th><th>1</th><th>2</th><th>4</th><th>6</th><th>8</th><th>10</th><th>12</th><th>16</th></tr></thead><tbody><tr><th>Total run time (s)</th><td>3788</td><td>2337</td><td>1376</td><td>976</td><td>840</td><td>765</td><td>747</td><td>717</td></tr><tr><th>Speedup factor <i>S</i></th><td>1</td><td>1.62</td><td>2.75</td><td>3.88</td><td>4.51</td><td>4.95</td><td>5.07</td><td>5.29</td></tr><tr><th>Parallel efficiency (<i>S</i>/<i>n</i>)</th><td>1</td><td>0.81</td><td>0.69</td><td>0.65</td><td>0.56</td><td>0.50</td><td>0.42</td><td>0.33</td></tr></tbody></table>

the distribution or built from source code. Both Open-MPI (http://www.open-mpi.org/software/ompi/v1.8/, last access: 30 October 2019) and mpich2 (http://www.mpich.org/downloads/, last access: 30 October 2019) work, but testing on some systems indicates slightly better performance with OpenMPI. As for the other libraries, the MPI library names and paths need to be adapted in the make-file. The `MPIF90` variable sets the Fortran compiler wrapper (usually `mpifort` or `mpif90`; in the case of coexisting OpenMPI and mpich installations, wrappers called `mpif90.openmpi` or `mpif90.mpich` may be defined). Compilation of the parallel version should then be done by the following.

### make mpi

This will produce an executable file `FLEXPART MPI`. If executed, this will run on a single processor and should produce results identical to the serial version. To activate the parallel features, the executable must be run through an MPI launcher (here it is important to use the launcher corresponding to the MPI library that was used for the compilation).

### mpirun -n <number> FLEXPART MPI

In this example, `<number>` specifies the number of processes one wishes to launch. For some installations, `mpirun` is called `mpiexec`, or, in the case of coexisting OpenMPI and mpich2 installations, `mpiexec.openmpi` or `mpiexec.mpich2`, respectively. Many command-line options exist for `mpirun` that can be helpful for improving performance, e.g., processor binding. For a list of these options, see `mpirun -help`. In practice, the optimal number of cores for a given simulation will depend on the size of the problem and the hardware availability among other factors.

## 5 FLEXPART input

In this section, we describe the different FLEXPART input files and, where appropriate, changes that have occurred since the last publication (Stohl et al., 2005). FLEXPART needs the following three types of input files.

1. The text file `pathnames` is located by default in the directory where FLEXPART is executed. It must contain at least four lines: first, the path to the directories where run-defining input files are located (the so-called `options` directory); second, the path where output

files are created; third, the path to the meteorological input GRIB files; and, fourth, the path to the so-called `AVAILABLE` file (see point 3). The last two lines can be repeated if nested input data shall be used. For each nesting level, one line for the GRIB data directory and one for the corresponding `AVAILABLE` file are needed.

2. The files containing the run-defining settings are located in a subdirectory (given in line 1 of `pathnames`) by default called `options` (see Table 1). The settings, which control FLEXPART's physics and program flow, are stored in different text files listed in Table 7 and described in Sect. 5.1. In addition, the `options` directory contains data files that are not usually changed by the user.

3. The meteorological input data, one file for each input time, are stored in GRIB format in a common directory (specified in line 3 of `pathnames`). To enable FLEXPART to find these files, a file usually named `AVAILABLE` (given in line 4 of `pathnames`) contains a list of all available meteorological input files and their corresponding time stamps. Additional files containing nested input data may also be provided. In this case, a separate file containing the input file names (e.g., named `AVAILABLE_NESTED`) must be given. Date and time entries in the `AVAILABLE*` files for mother and nested fields must be identical. Details on the meteorological input data are given in Sect. 5.2.

### 5.1 Run-defining settings: the options directory

Here, we give an overview of the information provided in the run-defining FLEXPART user input files listed in Table 7. In previous versions of FLEXPART, these files were formatted text files (coming alternatively in a long and a short format). For backward compatibility, these plain text formats are still supported. However, FLEXPART v10.4 also allows for the use of namelists, a standard Fortran feature whereby values are provided in a list with elements of the form `name=value`. When FLEXPART is started, it tries to open the files as namelists, and if this is not working, it expects the files to be in one of the two old plain text formats. We encourage users to update their input files to namelists for two reasons. Firstly, FLEXPART now has default user options for all input settings so that users only need to set those options that they want to deviate from the defaults. Secondly,
<table><caption>**Table 7.** Alphabetical list of the run-defining input files (upper part) and static input files (lower parts), usually contained in a directory called *options*. Processing of files marked with * depends on the run specifications. The other files are always read in.</caption><thead><tr><th>File name</th><th>Content</th></tr></thead><tbody><tr><td>AGECLASSES*</td><td>Age class definitions</td></tr><tr><td>COMMAND</td><td>Main control parameters</td></tr><tr><td>OUTGRID</td><td>Output grid definition</td></tr><tr><td>OUTGRID_NEST*</td><td>Nested output grid definition</td></tr><tr><td>RECEPTORS*</td><td>Receptor locations for receptor kernel output</td></tr><tr><td>RELEASES</td><td>Specification of the sources (forward run) or receptors (backward run)</td></tr><tr><td>SPECIES/</td><td>Directory containing files with definitions of physical and chemical parameters of species referenced in <i>RELEASES</i></td></tr><tr><td>IGBP_int1.dat</td><td>Land cover input data</td></tr><tr><td>surfdata.t</td><td>Roughness length, leaf area index for different land cover types</td></tr><tr><td>surfdepo.t</td><td>Seasonal surface resistances for different land cover types</td></tr><tr><td>OH_variables.bin*</td><td>OH field</td></tr></tbody></table>

namelists make it easier to add new user options, which may be required in future versions of FLEXPART. Thus, plain text input files may not be supported in future versions of FLEXPART. Examples for all formats of the user input files are contained in the FLEXPART distribution.

To convert user input files of any format to namelist format, the switch `nmlout=.TRUE.` (in file `com_mod.f90`) must be set before compilation. Then, run-defining user input files are written out in namelist format in the `output` directory, with the appendix `.namelist` added to the input file name (e.g., `COMMAND.namelist`). This feature also improves the traceability of FLEXPART model results and makes simulations easily reproducible by documenting the settings used for the model run.

In the following, we provide reference tables of the run-defining user input files including default settings (in bold) when using the namelist format. Notice that the default values are appropriate for regional-scale simulations, but simulations on smaller scales or with higher accuracy may need adjustments (in particular, shorter time steps and the use of the new CBL scheme).

### 5.1.1 File COMMAND

The `COMMAND` file contains the user settings controlling the simulation and the behavior of the run. The default `COMMAND` file contains a namelist `&COMMAND`, for which Table 8 provides a complete listing of all settings with their meaning and preset default values. It is important that users of previous FLEXPART versions who choose to use plain text input files update their `COMMAND` file, since new parameters have been added. However, the `cblflag` (and any option added afterwards) must be provided in namelist format in any case.

### 5.1.2 File RELEASES

The `RELEASES` file contains information related to when and where the particles are introduced in the simulation and other properties of the release points (e.g., the chemical species simulated). It consists of a namelist `&RELEASES_CTRL` that specifies header information. The header gives the total number of different species (i.e., different substances) to be released, followed by a corresponding list of the FLEXPART species numbers `nnn`, and `SPECIES_nnn` files define the species' physical properties (see Sect. 5.1.3). Following the header, there is an arbitrary number of namelists `&RELEASE` defining each release. For each such release, the following is given: the starting and ending time, the location and extension, the masses released (one value for each released species), and the number of particles to be released, as well as a comment string. The content of the `RELEASES` file is summarized in Table 9.

### 5.1.3 SPECIES files

The subdirectory `options/SPECIES/` needs to contain one or more files named `SPECIES_nnn`. For each species `nnn` listed in the header section of the `RELEASES` file, such a `SPECIES_nnn` file must exist. The parameters in the `SPECIES_nnn` file, contained in the namelist `&SPECIES_PARAMS`, set the species name and define the physicochemical properties of the species; they are described in Table 10. These are important for simulating radioactive or chemical decay, wet deposition (scavenging) for gases and aerosols, dry deposition for gases and aerosols, particle settling, and chemical reaction with the OH radical. Some parameters are only necessary for gas tracers and some are only necessary for aerosol tracers; thus, a namelist does not need to contain all parameters for both gases and particles. Optionally, since FLEXPART version 6.0, information about temporal emission variations can be added at the end of the file.
<table><caption>**Table 8.** Contents of the user input file COMMAND. Variable names with their meaning and all possible values are listed. Where appropriate, default values are given in bold. Note that not all input parameter combinations are allowed.</caption><thead><tr><th></th><th>Variable name</th><th>Description</th><th>Value (default)</th></tr></thead><tbody><tr><td>1</td><td>LDIRECT</td><td>Simulation direction in time</td><td>1 (forward) or -1 (backward)</td></tr><tr><td>2</td><td>IBDATE</td><td>Start date of the simulation</td><td>YYYYMMDD: YYYY=year, MM=month, DD=day</td></tr><tr><td>3</td><td>IBTIME</td><td>Start time of the simulation</td><td>HHMISS: HH hours, MI=minutes, SS=seconds. UTC zone.</td></tr><tr><td>4</td><td>IEDATE</td><td>End date of the simulation</td><td>Same format as IBDATE</td></tr><tr><td>5</td><td>IETIME</td><td>End time of the simulation</td><td>Same format as IBTIME</td></tr><tr><td>6</td><td>LOUTSTEP</td><td>Interval of model output</td><td>Average concentrations are calculated every LOUTSTEP (10 800 s)</td></tr><tr><td>7</td><td>LOUTAVER</td><td>Concentration averaging interval, instantaneous for value of zero</td><td>(10 800 s)</td></tr><tr><td>8</td><td>LOUTSAMPLE</td><td>Numerical sampling rate of output, higher statistical accuracy with shorter intervals</td><td>(900 s)</td></tr><tr><td>9</td><td>ITSPLIT</td><td>Time constant for particle splitting (particles are split in two after a given time)</td><td>(999 999 999 s)</td></tr><tr><td>10</td><td>LSYNCTIME</td><td>All processes are synchronized to this time interval; it has to divide all values above</td><td>(900 s)</td></tr><tr><td>11</td><td>CTL</td><td>Factor by which particle transport time step in the ABL must be smaller than the Lagrangian timescale t<sub>l</sub>; resulting time steps can be shorter than LSYNCTIME; LSYNCTIME is used if CTL &lt; 0</td><td>&gt; 1 for time steps shorter than t<sub>l</sub>; if CTL &lt; 0, a purely random walk simulation is done (-5.0)</td></tr><tr><td>12</td><td>IFINE</td><td>Additional reduction factor for time step used for vertical transport only considered if CTL &gt; 1</td><td>Positive integer (4)</td></tr><tr><td>13</td><td>IOUT</td><td>Switch determining the output type</td><td>(1) mass concentration (residence time backwards),<br>2 volume mixing ratio,<br>3 both 1 and 2,<br>4 plume trajectories,<br>5 both 1 and 4.<br>Add 8 for NetCDF output</td></tr><tr><td>14</td><td>IPOUT</td><td>Switch for particle position output</td><td>(0) no particle output,<br>1 particle output every output interval,<br>2 only at the end of the simulation (useful, e.g., for warm start)</td></tr><tr><td>15</td><td>LSUBGRID</td><td>Increase in ABL heights due to subgrid-scale orographic variations</td><td>(0)=off, 1=on</td></tr><tr><td>16</td><td>LCONVECTION</td><td>Switch for convection parameterization</td><td>0=off, (1)=on</td></tr><tr><td>17</td><td>LAGESPECTRA</td><td>Switch for calculation of age spectra (needs file AGECLASSES)</td><td>(0)=off, 1=on</td></tr><tr><td>18</td><td>IPIN</td><td>Warm start simulation, restarted from a particle dump (needs partposit_end file from previous simulation)</td><td>(0)=no, 1=yes</td></tr><tr><td>19</td><td>IOER</td><td>Separate output fields for each location in the RELEASE file</td><td>(0)=no, 1=yes</td></tr><tr><td>20</td><td>IFLUX</td><td>Output of mass fluxes through output grid box boundaries (northward, southward, eastward, westward, upward and downward)</td><td>(0)=off, 1=on</td></tr></tbody></table>
<table><caption>**Table 8. Continued.**</caption><tr><td>21</td><td>MDOMAINFILL</td><td>Switch for domain-filling calculations: particles are initialized to reproduce air density or stratospheric ozone density; for limited-area simulations, particles are generated at the domain boundaries</td><td>(0)=no, 1 like air density, 2 stratospheric ozone tracer</td></tr><tr><td>22</td><td>IND_SOURCE</td><td>Unit to be used at the source; see Seibert and Frank (2004); Eckhardt et al. (2017)</td><td>(1)=mass, 2=mass mixing ratio</td></tr><tr><td>23</td><td>IND_RECEPTOR</td><td>Unit to be used at the receptor; see Seibert and Frank (2004); Eckhardt et al. (2017)</td><td>(1)=mass, 2=mass mixing ratio, 3=bwd. wet. dep., 4=bwd. dry. dep.</td></tr><tr><td>24</td><td>MQUASILAG</td><td>Quasi-Lagrangian mode to track individual numbered particles</td><td>(0)=off, 1=on</td></tr><tr><td>25</td><td>NESTED_OUTPUT</td><td>Switch to produce output also for a nested domain</td><td>(0)=no, 1=yes</td></tr><tr><td>26</td><td>LINIT_COND</td><td>Switch to produce output sensitivity to initial conditions given in concentration or mixing ratio units (in backwards mode only)</td><td>(0)=no, 1=mass concentration, 2=mass mixing ratio</td></tr><tr><td>27</td><td>SURF_ONLY</td><td>Output of SRR for fluxes only for the lowest model layer, most useful for backward runs when LINIT_COND set to 1 or 2</td><td>(0)=no, 1=yes</td></tr><tr><td>28</td><td>CBLFLAG</td><td>Skewed rather than Gaussian turbulence in the convective ABL; when turned on, very short time steps should be used (see CTL and IFINE)</td><td>(0)=no, 1=yes</td></tr><tr><td>29</td><td>OHFIELDS_PATH</td><td>Default path for OH file</td><td>(50 m<sup>2</sup> s<sup>-1</sup>)</td></tr><tr><td>30</td><td>d_trop</td><td>Tropospheric horizontal turbulent diffusivity D<sub>h</sub></td><td>(0.1 m<sup>2</sup> s<sup>-1</sup>)</td></tr><tr><td>31</td><td>d_strat</td><td>Stratospheric vertical turbulent diffusivity D<sub>z</sub></td><td>(0.1 m<sup>2</sup> s<sup>-1</sup>)</td></tr></table>

<table><caption>**Table 9.** Contents of the user input file RELEASES.</caption><thead><tr><th>Variable name</th><th>Description</th><th>Format, valid values, variable type</th></tr></thead><tbody><tr><td colspan="3">Header (written only once and valid for all releases)</td></tr><tr><td>NSPEC</td><td>Total number of species</td><td>Integer number</td></tr><tr><td>SPECNUM_REL</td><td>Species numbers in dir. SPECIES</td><td>Integer array of size NSPEC</td></tr><tr><td colspan="3">For each release</td></tr><tr><td>IDATE1</td><td>Release start date</td><td>YYYYMMDD: YYYY=year, MM=month, DD=day</td></tr><tr><td>ITIME1</td><td>Release start time in UTC</td><td>HHMISS: HH hours, MI=minutes, SS=seconds; integer</td></tr><tr><td>IDATE2</td><td>Release end date</td><td>Same format as IDATE1</td></tr><tr><td>ITIME2</td><td>Release end time in UTC</td><td>Same format as ITIME1</td></tr><tr><td>LON1</td><td>Left longitude of release box</td><td>-180 &lt; LON1 &lt; 180, or according to input winds; real</td></tr><tr><td>LON2</td><td>Right longitude of release box</td><td>Same format as LON1; real</td></tr><tr><td>LAT1</td><td>Lower latitude of release box</td><td>-90 &lt; LAT1 &lt; 90, or according to input winds; real</td></tr><tr><td>LAT2</td><td>Upper latitude of release box</td><td>Same format as LAT1; real</td></tr><tr><td>ZKIND</td><td>Reference level</td><td>1: meters above ground, 2: meters above sea level, 3: pressure (hPa); integer</td></tr><tr><td>Z1</td><td>Lower height of release box</td><td>Meters above reference level; real</td></tr><tr><td>Z2</td><td>Upper height of release box</td><td>Meters above reference level; real</td></tr><tr><td>PARTS</td><td>Total number of particles to be released</td><td>Integer ≥ 1</td></tr><tr><td colspan="3">For each species (NSPEC times)</td></tr><tr><td>MASS</td><td>Total mass emitted</td><td>in, e.g., kilograms or unitless for mixing ratio; 1 in backward mode; real</td></tr><tr><td>COMMENT</td><td>Comment</td><td>40-character string (e.g., name of release point)</td></tr></tbody></table>
Notice that the format of the `SPECIES_nnn` files has changed from previous FLEXPART versions and users need to update their files accordingly. The use of `SPECIES_nnn` files from older FLEXPART versions may lead to run time errors or erroneous results.

The following specifies the parameters associated with each physicochemical process simulated.

- Radioactive or chemical decay: set with `pdecay`; off if `pdecay<0`.

- Wet deposition for gases: set with `pweta_gas`, `pwetb_gas` (for below-cloud) and `phenry` (for in-cloud). Switch off for both in- and below-cloud if either `pweta_gas` or `pwetb_gas` is negative.

- Wet deposition for aerosols: set with `pccn_aero`, `pin_aero` for in-cloud scavenging and `pcrain_aero`, `pcsnow_aero` and `pdquer` for below-cloud scavenging.

- Dry deposition for aerosols: set with `pdensity`, `pdquer` and `psigma`; off if `pdensity<0`.

- Dry deposition for gases: set with `phenry`, `pf0` and `preldiff`; off if `preldiff<0`. Alternatively, a constant dry deposition velocity `pdryvel` can be given.

- Settling of particles: set with `pdensity` and `pdquer`.

- OH reaction: chemical reaction with the OH radical can be turned on by giving parameter `pohconst` (cm³ molecule⁻¹ s⁻¹), `pohdconst` (K) and `pohnconst` (no unit) positive values; defined by Eq. (13).

- Emission variation: emission variation during the hours (local time) of the day and during the days of the week can be specified. Factors should be 1.0 on average to obtain unbiased emissions overall. The area source factors (useful, e.g., for traffic emissions) are applied to emissions with a lower release height below 0.5 m above ground level (a.g.l.) and the point source factors (useful, e.g., for power plant emissions) to emissions with a lower release height than 0.5 m a.g.l. Default values are 1.0.

### 5.1.4 File OUTGRID

The `OUTGRID` file specifies the domain and grid spacing of the three-dimensional output grid. Note that in a Lagrangian model, the domain and resolution of the gridded output are totally independent from those of the meteorological input (apart from the fact that the output domain must be contained within the computational domain). The `OUTGRID` file contains a namelist `&OUTGRID` specifying all parameters. The variables read in for this file and all the following input files have not changed in recent FLEXPART versions; thus, for

further explanation, see Stohl et al. (1995). Example files can be found in the `options` directory in the `FLEXPART` distribution.

### 5.1.5 File OUTGRID_NEST

Output can also be produced on one nested output grid with higher horizontal resolution, defined in the file `OUTGRID_NEST`, but with the same vertical resolution as given in `OUTGRID`. The `OUTGRID_NEST` file contains a namelist `&OUTGRIDN` specifying all parameters.

### 5.1.6 File AGECLASSES

The option to produce age class output can be activated in the `COMMAND` file. The file `AGECLASSES` then allows for the definition of a list of times (in seconds, in increasing order) that define the age classes used for model output. With this option, the model output (e.g., concentrations) is split into contributions from particles of different age, defined as the time passed since the particle release. Particles are dropped from the simulation once they exceed the maximum age, allowing their storage locations to be reused for new particles. This is an important technique to limit the memory usage for long-term simulations. Thus, even if the user is not interested in age information per se, it may often be useful to set one age class to define a maximum particle age.

### 5.1.7 File RECEPTORS

In addition to gridded model output, it is also possible to define receptor points. With this option output can be specifically produced for certain points at the surface in addition to gridded output. The `RECEPTORS` file contains a list with the definitions of the receptor name, longitude and latitude. If no such file is present, no receptors are written to `output`.

### 5.1.8 Static data input files

Several files contain static input data that are not usually modified by the user. These are (by default) also located in the `options` directory. If modeling a species requires calculating OH reactions, an OH field stored in file `OH_variables.bin` needs to be present. The file `IGBP_int1.dat` is a land cover inventory; file `surfdata.t` gives the roughness length and leaf area index of the different land cover types, and file `surfdepo.t` contains surface resistances for dry deposition calculations.

## 5.2 Meteorological data and preprocessing routines

FLEXPART can be run with meteorological input data for global domains or for smaller, limited-area domains. The FLEXPART computational domain always corresponds to this mother domain set by the input data, while the output domain can be smaller. FLEXPART can also ingest
<table><caption>**Table 10.** FLEXPART variables set in the user input file `SPECIES_nnn` for species number `nnn`. Note that the variable names given in the input namelist are the same as used subsequently in FLEXPART but with a prepended letter `p` (for parameter). For instance, `pspecies` corresponds to `species`.</caption><thead><tr><th>Variable name</th><th>Description</th><th>Unit</th><th>Type</th></tr></thead><tbody><tr><td>1 pspecies</td><td>Tracer name</td><td>no units</td><td>string</td></tr><tr><td>2 pdecay</td><td>Species half-life for radioactive or chemical decay; off if pdecay&lt;0</td><td>seconds</td><td>real</td></tr><tr><td>3 pweta_gas</td><td>Gases wet deposition: below-cloud scavenging parameter A (for precip. of 1 mm h<sup>-1</sup>)</td><td>1 s<sup>-1</sup></td><td>real</td></tr><tr><td>4 pwetb_gas</td><td>Gases wet deposition: below-cloud scavenging parameter B (dependency on precip. rate)</td><td>1</td><td>real</td></tr><tr><td>5 pcrain_aero</td><td>Aerosols wet deposition: below-cloud scavenging rain collection efficiency moderator for rain C<sub>rain</sub></td><td>1</td><td>real</td></tr><tr><td>6 pcsnow_aero</td><td>Aerosols wet deposition: below-cloud scavenging rain collection efficiency moderator for snow C<sub>snow</sub></td><td>1</td><td>real</td></tr><tr><td>7 pccn_aero</td><td>Aerosols wet deposition: in-cloud scavenging, cloud condensation nuclei efficiency CCN<sub>eff</sub></td><td>1</td><td>real</td></tr><tr><td>8 pin_aero</td><td>Aerosols wet deposition: in-cloud scavenging, ice nuclei efficiency IN<sub>eff</sub></td><td>1</td><td>real</td></tr><tr><td>9 preldiff</td><td>Gases dry deposition: ratio D = D<sub>H<sub>2</sub>O</sub>/D<sub>i</sub>, D<sub>H<sub>2</sub>O</sub> of the diffusivity of H<sub>2</sub>O to the diffusivity of the component D<sub>i</sub> (the diffusivity of the species in the SPECIES_nnn file); dry deposition of gases is switched off by negative D</td><td>1</td><td>real</td></tr><tr><td>10 phenry</td><td>Gases dry deposition and in-cloud scavenging: Henry's constant H</td><td>M atm<sup>-1</sup></td><td>real</td></tr><tr><td>11 pf0</td><td>Gases dry deposition: reactivity factor for oxidation of biological substances relative to that of ozone; (0 ≤ f<sub>0</sub> ≤ 1)</td><td>1</td><td>real</td></tr><tr><td></td><td>For nonreactive species f<sub>0</sub> is 0, for slightly reactive species it is 0.1 and for highly reactive species it is 1</td><td></td><td></td></tr><tr><td>12 pdensity</td><td>Aerosols dry deposition and settling: particle density ρ</td><td>kg m<sup>-3</sup></td><td>real</td></tr><tr><td>13 pdquer</td><td>Aerosol dry deposition, aerosol wet deposition: below-cloud scavenging: particle mean diameter d̄</td><td>m</td><td>real</td></tr><tr><td></td><td>Decides whether its gas(&lt;=0) or aerosol (&gt; 0)</td><td></td><td></td></tr><tr><td>14 psigma</td><td>Aerosol dry deposition: species diameter normalized deviation (σ &gt; 1)</td><td>1</td><td>real</td></tr><tr><td>15 pdryvel</td><td>Gases dry deposition: dry deposition velocity (only used if preldiff and pdensity &lt; 0)</td><td>m s<sup>-1</sup></td><td>real</td></tr><tr><td>16 pweightmolar</td><td>Gases: species molar weight, used for volume mixing ratio (pptv) output</td><td>g mol<sup>-1</sup></td><td>real</td></tr><tr><td>17 pohcconst</td><td>Gases OH reaction: C</td><td>cm<sup>3</sup> molecule s<sup>-1</sup></td><td>real</td></tr><tr><td>18 pohdconst</td><td>Gases OH reaction: D</td><td>K</td><td>real</td></tr><tr><td>19 pohnconst</td><td>Gases OH reaction: N</td><td>1</td><td>real</td></tr><tr><td>20 parea_hour</td><td>Emission variation factor (area source) for hour of the day, starting with 00:00–01:00 local time, 24 values</td><td>1</td><td>real</td></tr><tr><td></td><td>Local time from longitude, no correction for summer time</td><td></td><td></td></tr><tr><td>21 parea_dow</td><td>Emission variation factor (area source) for day of the week, starting with Monday, 7 values</td><td>1</td><td>real</td></tr><tr><td>22 ppoint_hour</td><td>Emission variation factor (point source) for hour of the day, starting with 00:00–01:00 local time, 24 values</td><td>1</td><td>real</td></tr><tr><td></td><td>Local time from longitude, no correction for summer time</td><td></td><td></td></tr><tr><td>23 ppoint_dow</td><td>Emission variation factor (point source) for day of the week, starting with Monday, 7 values</td><td>1</td><td>real</td></tr></tbody></table>

higher-resolution meteorological input data in subdomains of the mother domain. Such nested data must be available for the exact same times as those for the mother domain, checked by FLEXPART by comparing the time stamps in the two `AVAILABLE(_NESTED)` files. There is no nesting in the vertical direction and the poles must not be contained in any nest. To automatically produce the `AVAILABLE(_NESTED)` files, a Python script is available from the FLEXPART website (https://flexpart.eu/wiki/FpInputMetMkavail, last access: 25 June 2018) that checks which input files are present and then creates this file in the required format.

Compilation of FLEXPART v10.4 produces a single executable that automatically detects whether the meteorological input data come from the ECMWF IFS or NCEP GFS and whether they are in GRIB-1 or in GRIB-2 format. Nev-
ertheless, certain parameters may need to be adapted in `par_mod.f90` to the size of the meteorological input files (array dimensions), and the input grid may need to be shifted relative to the output grid (parameter `nxshift`). In the following, we describe how meteorological input data appropriate for FLEXPART can be retrieved from the ECMWF and NCEP.

### 5.2.1 ECMWF data retrieval

ECMWF data can be comprised of analysis and/or forecast data from the operational IFS data stream or specific re-analysis projects. For operational data, the meteorological fields can currently have a maximal temporal resolution of 1 h (more frequent data are not available), a vertical resolution of 137 model levels and $0.1° \times 0.1°$ horizontal resolution on a regular latitude-longitude grid. Other ECMWF datasets are not available at such high horizontal resolution. For example, ERA-Interim reanalysis data (Dee et al., 2011) with $1° \times 1°$ latitude-longitude resolution and 60 vertical levels can be retrieved 3-hourly by mixing 6-hourly analysis and 3 h forecast fields, but higher resolution is not available from the standard archive. The new Copernicus reanalysis ERA5 provides 1-hourly analysis fields with 137 model levels and a horizontal resolution of 31 km ($0.28125°$). Notice that access to some datasets, in particular the operational forecasts, is restricted and requires specific access (https://www.ecmwf.int/en/forecasts/accessing-forecasts, last access: 23 June 2018). However, reanalysis data (https://software.ecmwf.int/wiki/display/WEBAPI/Available+ECMWF+Public+Datasets, last access: 23 June 2018) are publicly available.

The IFS is a global model that uses spectral representation with spherical harmonics for the dynamical part and a grid-point representation on a reduced Gaussian grid for the physical part. However, FLEXPART needs the input data on a regular latitude-longitude grid, and thus IFS data have to be preprocessed. With respect to the vertical coordinate system, the data need to be on the native ECMWF model levels ($\eta$ levels), which are subsequently transformed within FLEXPART to a terrain-following vertical coordinate system.

As explained above, each ECMWF dataset has its own specific temporal and spatial resolution, and the meteorological parameters provided can be different from dataset to dataset. To produce meteorological GRIB files suitable for FLEXPART input from these different datasets, a software called `flex_extract` (current version 7.0.4) has been developed specifically for this purpose. In order to prepare the GRIB files from the ECMWF Meteorological Archival and Retrieval System (MARS; https://software.ecmwf.int/wiki/download/attachments/45759146/mars.pdf, last access: 24 June 2018), several retrieval requests using the MARS command language and some further processing steps are needed. Since all ECMWF datasets need to be handled differently and some may not even contain all information needed for FLEXPART, `flex_extract` has a focus on

some of the most important ones for driving FLEXPART. These are, in particular, the reanalysis datasets ERA-Interim (Dee et al., 2011), CERA-20C (the coupled climate reanalysis of the 20th century; Laloyaux et al., 2018) and the latest reanalysis ERA5, as well as data from the operational IFS stream. Each file (one for each time step) prepared by `flex_extract` for FLEXPART consists of a set of model-level and surface data as a combination of analysis and forecast fields depending on availability. For example, certain variables such as precipitation may only be available in forecast fields, whereas other data are also contained in analysis fields; `flex_extract` seeks an optimum combination of such data. Note that some parameters are stored as time-accumulated fields in the ECMWF archives and `flex_extract` calculates the instantaneous fluxes out of them (e.g., precipitation fluxes). For more details on this process of de-accumulation, see Hittmeir et al. (2018). Since FLEXPART needs the pressure hybrid coordinate vertical velocity as used in the ECMWF model, an important feature of `flex_extract` is the computation of this parameter from the horizontal wind field (see Stohl et al., 2001) for ERA-Interim and for the years when it was not operationally archived in MARS (before 2009).

The ECMWF is a European intergovernmental organization that grants full access to its multi-petabyte MARS archive for their member and cooperating states. Users with a full-access account can run `flex_extract` v7.0.4 directly on ECMWF servers or via a local gateway server. This mode is also required to retrieve the most recent operational data from the ECMWF. Users from member or cooperating states interested in this mode should contact the computing representative from their national meteorological service to obtain an account. Users from other countries worldwide can self-register at the ECMWF for a public account to be able to retrieve the public datasets (i.e., most reanalysis products); `flex_extract` v7.0.4 makes use of the WebAPI (https://software.ecmwf.int/wiki/display/WEBAPI/ECMWF+Web+API+Home, last access: 24 June 2018) tool provided by the ECMWF to access the data from outside their systems. This tool can distinguish between public and member state users. Therefore, it is also a convenient option for member and cooperating state users who only need data older than a few days from the operational stream or reanalysis data. A full-access account to ECMWF servers is no longer needed in this case.

The `flex_extract` software v7.0.4 is a set of Python routines combined with a Fortran program for faster computation of grid transformations and vertical velocity calculation. A Python 2.7 interpreter with several common modules, such as NumPy and date-time, are required and usually included in the Anaconda distribution (https://www.anaconda.com/download/, last access: 25 June 2018). Additionally, a Fortran compiler, the ECMWF WebAPI tool, the GRIB-API or ecCodes module, and the Emoslib interpolation library have to be available. Note that the GRIB-
API (or ecCodes) module has to be available for Python as
well as for Fortran. Installation instructions can be found
at ECMWF websites directly or in the Software Installation
Plan for flex_extract. Knowledge of Python, although
helpful, is not necessary for using the retrieval scripts. A
certain knowledge of the ECMWF dataset to be retrieved
is useful to understand the composition of retrievals, but
many basic examples of CONTROL files are provided in the
flex_extract distribution. These CONTROL files deter-
mine the key parameters for the flex_extract MARS
retrievals and can be adapted to change domain as well as
spatial and temporal resolution. Even for these few parame-
ters the user should check for availability upfront. For exam-
ple, ERA-Interim data have a maximum grid resolution of
0.75° × 0.75° and 6-hourly temporal resolution for the pub-
lic dataset.

The *flex_extract* v7.0.4 software is included in
the FLEXPART v10.4 file tree under the directory
*preprocess* (see Table 1). It can also be downloaded from
https://flexpart.eu/ (last access: 30 October 2019) (https:
//flexpart.eu/downloads/62, last access: 30 October 2019)
as a tarball. For more details the reader is referred to the
*flex_extract* v7.0.4 user documentation (e.g., the Soft-
ware Installation Plan - SIP.pdf - and the Software User Tu-
torial - SUT.pdf) in *preprocess/flex_extract*.

5.2.2 NCEP data retrieval

Meteorological data from the NCEP GFS are freely avail-
able, easily accessible and ingested by FLEXPART on pres-
sure levels, unlike ECMWF data. These pressure-level data
have lower resolution than model-level data but offer the
advantage of great consistency between different datasets.
Therefore, preprocessing NCEP data is much simpler than
ECMWF data and limited to precipitation data, which are
available only in forecast fields.

Both operational analysis data and several reanalysis
datasets are available. Notice that NCEP also provides fore-
cast data for free, which are not available from the ECMWF
even for member state users except for national meteorologi-
cal services or users with a special contract. The data retrieval
from NCEP is described in a wiki page on the FLEXPART
website (https://flexpart.eu/wiki/FpInputMetGfs, last access:
8 July 2018), where a script for downloading NCEP data can
also be found. Operational GFS data can be downloaded by
simple FTP or wget from a rolling archive of the meteoro-
logical forecast and analysis data (http://www.ftp.ncep.noaa.
gov/data/nccf/com/gfs/prod/, last access: 8 July 2018) un-
der the catalog gfs.YYYYMMDDHH, which contains fields in
GRIB-2 format. Six-hourly NCEP FNL (Final) Operational
Model Global Tropospheric Analyses (http://rda.ucar.edu/
datasets/ds083.2/, last access: 30 October 2019) are avail-
able in near-real time since July 1999. These data are sim-
ilar to the operational analyses, but NCEP also ingests late-
incoming observation data for their production. Archived re-

analysis datasets are also available from NCEP, e.g., the Cli-
mate Forecast System Reanalysis (CFSR) Selected Hourly
Time-Series Products (http://rda.ucar.edu/datasets/ds093.1/
last access: 30 October 2019) for the period January 1979
to March 2011.

6 FLEXPART output

6.1 Output files overview

In the following we describe the FLEXPART output files to-
gether with changes made since the last documented FLEX-
PART version (Stohl et al., 2005). An overview of all pos-
sible output files is provided in Table 11. Notice that not all
these files are written out in every model run; the user set-
tings control which files are produced. At the beginning of a
run, FLEXPART records descriptive metadata in the binary
file header. This information is also written into the plain
text files header_txt (with the exception of the orogra-
phy data and release information). The release information is
written in header_txt_releases. Corresponding files
header_nested are produced if nested output is selected.

At each output time, FLEXPART produces files con-
taining the gridded output. Separate files are created
for every species and domain (mother and, if re-
quested, nest). The naming convention for these files is
grid_type_date_nnn. For forward runs, type can be
conc or pptv for concentrations and mixing ratios or flux
for 3-D mass fluxes across the grid cell faces (Stohl et al.,
2005, Sect. 8.5). For backward runs, type can be time for
the sensitivity of receptor concentrations to emission fluxes,
drydep for the sensitivity of receptor dry deposition to
emissions or wetdep for the sensitivity of receptor wet de-
position to emissions. For backward runs, there can also be
an output file grid_initial_nnn, which gives the re-
ceptor sensitivity to initial conditions; date denotes the date
and time for which the output is valid, and nnn is the species
number as specified in RELEASES. The list of the output
times is progressively written to the text file dates. For the
nested output, grid is replaced by grid_nested.

Wet and dry deposition fields in forward runs are calcu-
lated on the same horizontal output grid and are appended to
grid_conc_date_nnn and grid_pptv_date_nnn
files. The deposited matter is accumulated over the course
of a model run. It generally increases with model time, but
for species with radioactive decay, losses are possible. As for
long simulations small deposition amounts may be added to
already large deposited quantities, the default precision of the
deposition fields was changed from single (in older FLEX-
PART versions) to double precision to avoid numerical in-
accuracies when deriving instantaneous fluxes from accumu-
lated quantities.

For a list of points at the surface, concentrations or mixing
ratios in forward simulations can also be calculated indepen-
<table><caption>**Table 11.** List of FLEXPART output files and for which user settings ("switches") they are produced. See Table 12 for details on the units in the gridded output.</caption><thead><tr><th>Name</th><th>Format</th><th>Switches</th><th>Description of contents</th></tr></thead><tbody><tr><td>header</td><td>binary</td><td>default output</td><td>run metadata + ancillary data</td></tr><tr><td>header_txt</td><td>text</td><td>default output</td><td>human-readable run metadata (from COMMAND)</td></tr><tr><td>header_txt_releases</td><td>text</td><td>default output</td><td>human-readable run metadata (from RELEASES)</td></tr><tr><td>dates</td><td>text</td><td>default output</td><td>time series: dates of output files</td></tr><tr><td>grid_conc_date_nnn</td><td>binary (sparse array)</td><td>LDIRECT=1 IOUT=1,3,5</td><td>3-D tracer mass density + 2-D deposition</td></tr><tr><td>grid_pptv_date_nnn</td><td>binary (sparse array)</td><td>LDIRECT=1 IOUT=2,3</td><td>3-D tracer volume mixing ratio + 2-D deposition</td></tr><tr><td>grid_time_date_nnn</td><td>binary (sparse array)</td><td>LDIRECT=-1 IOUT=1</td><td>3-D sensitivity of atmospheric receptor to emissions</td></tr><tr><td>grid_drydep_date_nnn</td><td>binary (sparse array)</td><td>LDIRECT=-1 IOUT=1, IND_RECEPTOR=3</td><td>3-D sensitivity of dry deposition receptor to emissions</td></tr><tr><td>grid_wetdep_date_nnn</td><td>binary (sparse array)</td><td>LDIRECT=-1 IOUT=1, IND_RECEPTOR=4</td><td>3-D sensitivity of wet deposition receptor to emissions</td></tr><tr><td>grid_conc_date_nnn.nc</td><td>binary (NetCDF)</td><td>LDIRECT=1 IOUT=9,11,13</td><td>3-D tracer + 2-D wet and dry deposition</td></tr><tr><td>grid_time_date_nnn.nc</td><td>binary (NetCDF)</td><td>LDIRECT=-1 IOUT=9</td><td>3-D sensitivity of atmospheric receptor to emissions</td></tr><tr><td>grid_drydep_date_nnn.nc</td><td>binary (NetCDF)</td><td>LDIRECT=-1 IOUT=9, IND_RECEPTOR=3</td><td>3-D sensitivity of dry deposition receptor to emissions</td></tr><tr><td>grid_wetdep_date_nnn.nc</td><td>binary (NetCDF)</td><td>LDIRECT=-1 IOUT=9, IND_RECEPTOR=4</td><td>3-D sensitivity of wet deposition receptor to emissions</td></tr><tr><td>grid_initial_nnn</td><td>binary (sparse array)</td><td>LDIRECT=-1, LINIT_COND&gt;0</td><td>3-D sensitivity of receptor concentrations and deposition to initial conditions</td></tr><tr><td>partposit_date</td><td>binary</td><td>IPOUT=1,2<br>IPOUT=1,2 MQUASILAG=1</td><td>particle positions and meteorological data<br>particle positions and meteorological data numbered<br>consecutively</td></tr><tr><td>partposit_average_date</td><td>binary</td><td>IPOUT=3</td><td>time-averaged particle positions and meteorological<br>data<br>clustered trajectories</td></tr><tr><td>trajectories.txt</td><td>text</td><td>IOUT=4, 5</td><td></td></tr><tr><td>receptor_conc</td><td>binary</td><td>LDIRECT=1 IOUT=1,3,5,9,11,13</td><td>mass density at receptors</td></tr><tr><td>receptor_pptv</td><td>binary</td><td>LDIRECT=1 IOUT=2,3,10,11</td><td>volume mixing ratio at receptors</td></tr><tr><td>header_nested</td><td>binary</td><td>NESTED_OUTPUT=1</td><td>nest metadata + ancillary data</td></tr><tr><td>grid_conc_nested_date_nnn</td><td>binary (sparse array)</td><td>as for mother grid</td><td>as for mother grid in a<br>higher-resolution latitude-longitude grid</td></tr><tr><td>grid_pptv_nested_date_nnn</td><td>binary (sparse array)</td><td>+NESTED_OUTPUT=1</td><td></td></tr><tr><td>grid_time_nested_date_nnn</td><td>binary (sparse array)</td><td></td><td></td></tr><tr><td>grid_drydep_nested_date_nnn</td><td>binary (sparse array)</td><td></td><td></td></tr><tr><td>grid_wetdep_nested_date_nnn.nc</td><td>binary (NetCDF)</td><td></td><td></td></tr><tr><td>grid_conc_nested_date_nnn.nc</td><td>binary (NetCDF)</td><td></td><td></td></tr><tr><td>grid_time_nested_date_nnn.nc</td><td>binary (NetCDF)</td><td></td><td></td></tr><tr><td>grid_drydep_nested_date_nnn.nc</td><td>binary (NetCDF)</td><td></td><td></td></tr><tr><td>grid_wetdep_nested_date_nnn.nc</td><td>binary (NetCDF)</td><td></td><td></td></tr></tbody></table>
dently from the grid using a kernel method and recorded in the files `receptor_conc` and/or `receptor_pptv`.

If the particle dump option is activated, in addition to the gridded output, the particle coordinates together with additional variables such as pressure, humidity, density, tropopause height, ABL height and orography height are recorded in the binary files `partposit_date`. These data can be useful for a variety of different purposes, for instance diagnostics of the water cycle (Stohl and James, 2004). FLEXPART version 10.4 also has the new option to write out time-averaged particle positions and meteorological data. These are recorded in the files `partposit_average_date`. Such output may be useful to obtain, for instance, more representative heights for particles in the boundary layer, where particle positions change rapidly and this is not sampled sufficiently with instantaneous output. If plume trajectory mode is activated, for every release the positions of trajectory clusters representing the centers of mass of all released particles are recorded in the file `trajectories.txt` (Stohl et al., 2002, 2005, Sect. 10).

The physical unit used for the output data in the files `grid_conc_date_nnn` and `grid_time_date_nnn` depends on the settings of the switches `ind_source` and `ind_receptor`, following Table 12. It is noteworthy that the unit of mass mixing ratio can also be used in `grid_conc_date_nnn`. For forward runs, additional files `grid_pptv_date_nnn` can be created (setting `IOUT` to values of 2 or 3), which contain data such as volume mixing ratios (requires molar weight in `SPECIES_nnn` file). Source-receptor relationships (i.e., emission sensitivities) in backward mode for atmospheric receptors are written out in `grid_time_date_nnn` files; those for deposited mass are recorded in files `grid_wetdep_date_nnn` and `grid_drydep_date_nnn` (see Seibert and Frank (2004), Eckhardt et al. (2017), Sect. 2.5, and Table 12 for output units). Notice that the user can also provide different input units. For instance, if emissions in a forward run are specified in Becquerels (Bq), the output would be in nanobecquerels per cubic meter (nBq m⁻³) with `ind_source=1` and `ind_receptor=1`. Notice further that all gridded output quantities in FLEXPART are grid cell averages, not point values.

## 6.2 Sparse matrix output

Depending on the type of model run, the gridded output can contain many grid cells with zero values (e.g., dispersion from a point source, backward run from a single receptor). The output is therefore written in a sparse matrix format, which is specific to FLEXPART. The array containing the data to be written out is scanned for sequences of nonzero values. The number of sequences found is stored in an integer variable `sp_count_i`, and the field positions at which each sequence begins are stored in a 1-D integer array,

`sparse_dump_i`, using a one-dimensional representation of the output field. The total number of nonzero values is stored in `sp_count_r` and the nonzero values themselves in the real vector `sparse_dump_r`. Since all physical output quantities of FLEXPART are greater than or equal to zero, nonzero sequences are stored in `sparse_dump_r` with alternating signs, which allows for the separation of different sequences upon reading. Finally, all four variables are written out to the unformatted output file. This format replaces the compression used up to version 7 (the smallest of a full dump and a simple sparse matrix format), saving up to 60% of disk space. The sparse matrix data can be read, for example, with the functions `readgrid.f` (Fortran) and `flex_read.m` (MATLAB) described in Sect. 6.4

## 6.3 NetCDF output

FLEXPART v10.4 can also support output in NetCDF format if the NetCDF libraries are available. To activate NetCDF support, append `ncf=yes` to the `make` command. If FLEXPART is compiled and linked to the NetCDF libraries, output files in NetCDF format can be produced by adding 8 to the `IOUT` parameter in the input file `COMMAND`, e.g., `IOUT=9` corresponds to `IOUT=1` with the standard binary output; see Table 11 and Sect. 5.1.1. In the NetCDF module `netcdf_output_mod.f90` a parameter `write_releases` determines at compile time if the information on the releases should also be written to the NetCDF file. Only one NetCDF file is written that contains all species and all time steps. Both mother and nested output (if present) are contained in that file. Since the NetCDF output is specified in the climate and forecast (CF) format, any standard software can be used for displaying and processing the output (e.g., panoply, ncview). NetCDF output data files are compressed.

The NetCDF output file contains information on the run settings and the simulation grid from the `COMMAND` and `OUTGRID*` files. It also contains additional information in the header on the producing center, as listed in Table 13. The content of these attributes can be adapted in the file `netcdf_output_mod.f90` before compilation.

## 6.4 Post-processing routines

For the NetCDF output of FLEXPART, standard visualization tools, for example Panoply, can be used. For the sparse matrix binary output, several post-processing routines (MATLAB, Fortran, R, Python and IDL) have been developed in order to assist in the usage and analysis of these data. A number of post-processing tools are available online (https://flexpart.eu/wiki/FpOutput, last access: 16 August 2018). Note that some of these tools require reading a text string containing the model version. Since the length of this string changed in FLEXPART v9.2, the post-processing routines now require the allocation of a longer string.
<table><caption>**Table 12.** Physical units of the input (in file `RELEASES`) and output data for forward (files `grid_conc_date_nnn`) and backward (files `grid_time_date_nnn`) runs for the various settings of the unit switches `ind_source` and `ind_receptor` (for both switches, 1 refers to mass units, 2 to mass mixing ratio units). `IOUT` is 1 (or 9 for NetCDF output) except where indicated; "(dep.)" in lines 5 and 6 of the table refer to the deposition output provided in addition to the atmospheric output in files `grid_conc_date_nnn`.</caption><thead><tr><th>Direction</th><th>File name</th><th>ind_source</th><th>ind_receptor</th><th>Input unit</th><th>Output unit</th></tr></thead><tbody><tr><td rowspan="7">Forward</td><td>grid_conc*</td><td>1</td><td>1</td><td>kg</td><td>ng m<sup>-3</sup></td></tr><tr><td>grid_conc*</td><td>1</td><td>2</td><td>kg</td><td>ppt by mass</td></tr><tr><td>grid_conc*</td><td>2</td><td>1</td><td>l</td><td>ng m<sup>-3</sup></td></tr><tr><td>grid_conc*</td><td>2</td><td>2</td><td>l</td><td>ppt by mass</td></tr><tr><td>grid_conc*</td><td>1</td><td>1 or 2 (dep.)</td><td>kg</td><td>ng m<sup>-2</sup></td></tr><tr><td>grid_conc*</td><td>2</td><td>1 or 2 (dep.)</td><td>l</td><td>ng m<sup>-2</sup></td></tr><tr><td>grid_pptv* (IOUT=2, 3)</td><td>1</td><td>1</td><td>l</td><td>ppt by volume</td></tr><tr><td rowspan="10">Backward</td><td>grid_time*</td><td>1</td><td>1</td><td>l</td><td>s</td></tr><tr><td>grid_time*</td><td>1</td><td>2</td><td>l</td><td>s m<sup>3</sup> kg<sup>-1</sup></td></tr><tr><td>grid_time*</td><td>2</td><td>1</td><td>l</td><td>s kg m<sup>-3</sup></td></tr><tr><td>grid_time*</td><td>2</td><td>2</td><td>l</td><td>s</td></tr><tr><td>grid_wetdep*</td><td>1</td><td>3 (wet dep.)</td><td>l</td><td>m</td></tr><tr><td>grid_drydep*</td><td>1</td><td>4 (dry dep.)</td><td>l</td><td>m</td></tr><tr><td>grid_wetdep*</td><td>2</td><td>3 (wet dep.)</td><td>l</td><td>kg m<sup>-2</sup></td></tr><tr><td>grid_drydep*</td><td>2</td><td>4 (dry dep.)</td><td>l</td><td>kg m<sup>-2</sup></td></tr><tr><td>grid_initial*</td><td>1</td><td>1</td><td>l</td><td>m<sup>3</sup> kg<sup>-1</sup></td></tr><tr><td>grid_initial*</td><td>1</td><td>2</td><td>l</td><td>m<sup>3</sup> kg<sup>-1</sup></td></tr></tbody></table>

<table><caption>**Table 13.** Additional information in the NetCDF output file as attributes.</caption><tr><td>Conventions</td><td>CF-1.6 (NetCDF CF convention identifier)</td></tr><tr><td>Title</td><td>FLEXPART model output (content title)</td></tr><tr><td>Institution</td><td>producer string “institution” set in netcdf_output_mod.f90</td></tr><tr><td>Source</td><td>creation string “flexversion” model output set in FLEXPART.f90</td></tr><tr><td>History</td><td>date string with login and host name</td></tr><tr><td>References</td><td>Stohl et al. (2005)</td></tr></table>

Fortran routines are available for download on the FLEXPART website with the subroutines `readheader.f` for reading the header and `readgrid.f` for reading the grided binary fields. Analysis or plotting programs written in Fortran can call these subroutines.

There are also MATLAB tools working in a similar way as the Fortran routines, with `flex_header.m` for reading the header and `flex_read.m` for reading the data fields. If particle dumps were made, the MATLAB function `readpart.m` reads the corresponding data files (a similar Fortran code is also available).

The R programs available for post-processing FLEXPART output include routines to read the binary output in the `grid_conc` (or `grid_pptv`) and `grid_time` files and to plot maps. Routines are also available to plot trajectories on a map from the file `trajectories.txt` and to plot time series of concentrations (or mixing ratios) from the file `receptor_conc` (or `receptor_pptv`).

Several Python tools are available for reading FLEXPART data from release 8.0 and above. The module `reflexible`, available from the FLEXPART website and also at https://github.com/spectraphilic/reflexible (last access: 6 August 2018), enables the user to easily read and access the header and grid output data of the FLEXPART model runs. It provides a simple tool that facilitates consistent reading of both the original sparse matrix output files and the NetCDF output. Some basic plotting functionality is provided to quickly assess and validate runs or to look at the input parameters. An alternative Python tool is `Quicklook` that can be also downloaded from the https://flexpart.eu website.

## 7 Application examples

In this section we provide 38 examples of the FLEXPART model that serve three purposes: (1) verification of a new FLEXPART installation; (2) demonstration of the model ca-
pabilities for new users; and (3) confirmation of consistency
in the model output when code changes are made that should
not change the results. These examples do not represent an
exhaustive set of all possible model uses, but they are de-
signed to demonstrate and test different widely used func-
tionalities of the model.

All examples are variations of a default example case,
which uses the settings in the user input files as distributed
with the FLEXPART v10.4 code package. These default in-
put files are located in the directory options (Sect. 5) and
are consistent with the default meteorological data retrieved
from the ECMWF by the flex_extract package (Ap-
pendices A5 and B1). An AVAILABLE file fitting with these
input data is also distributed with FLEXPART. These default
settings are described in detail in Appendix B2.

Using the default example as a basis, the different func-
tionalities of the model can be activated by adequately chang-
ing certain parameters in the user input files, thereby generat-
ing 36 other example runs. We have categorized these exam-
ples into 10 different groups; each group explores different
capabilities of the model. Table 14 lists all examples and the
parameter changes needed to produce them. The first group
includes the default example and explores the different op-
tions for producing gridded model output (e.g., output units,
output formats) for a simple forward model run with a single
starting point over the North Atlantic. The second group of
examples introduces FLEXPART's backward simulation ca-
pability. The third group demonstrates different usages of the
particle dump output. The fourth group gives examples for
the use of mass vs. mass mixing ratio units at both the source
and the receptor and for both forward and backward simula-
tions to establish source–receptor relationships as in Seibert
and Frank (2004). The fifth group produces output for differ-
ent chemical species and aerosols. The sixth group illustrates
the use of nested output fields. Group seven is constituted
by a single domain-filling run, as used, for instance, in Stohl
and James (2004). Group eight contains settings for a back-
ward run providing 2-D sensitivities to gridded surface fluxes
and 3-D sensitivities to initial conditions, as they are typi-
cally required for the inverse modeling of greenhouse gases
(e.g., Thompson et al., 2017). Group nine shows the use of
the new skewed turbulence parameterization (Cassiani et al.,
2015). Group 10 shows the use of the new backward deposi-
tion (Eckhardt et al., 2017). Group 11 contains a forward 2 d
run simulating instantaneous emissions from a hypothetical
Grímsvötn eruption (Fig. 9).

The list of examples may be extended in the future to allow
for the testing of even more model features and to provide a
reference archive to see how FLEXPART results may change
as the code is being developed further. The user can get these
reference results from https://flexpart.eu. A quick reference
containing mean and maximum grid values for every exam-
ple is also referenced in Appendix B5.

The directory tests/examples/ contains scripts that
generate all the files necessary to run the examples. These

scripts, described in Appendix B3, generate the input files by
modifying the namelists in the default options directory
provided with the distribution. This is done by the bash script
gen_options_all.sh. For instance, the example "bwd"
is generated by changing the line containing the parameter
LDIRECT to -1 in the file COMMAND.

After the input data files are generated, all examples can
be executed interactively from the command line. Alterna-
tively, the script gen_batch_jobs_cl.sh generates a
batch script for each case (to be run from the command line
or using a workload manager such as SLURM). This proce-
dure automates the sample output generation. Once the out-
put files are created, they can be read using the tools in the
directory postprocess. They can be plotted and analyzed
with, e.g., the reading routines described in Sect. 6.4. In ad-
dition, some testing capabilities have been added. These are
presented in Appendix B.

8 Final remarks, outlook and future code development

In this paper, we have described the Lagrangian particle dis-
persion model FLEXPART v10.4; 2 decades ago, the model
code was developed mainly by one person, with specific code
input from a few other researchers. At that time, no specific
measures were needed to ensure code consistency, track code
changes or identify coding bugs. However, as the number of
FLEXPART users has grown substantially in recent years,
more and more people have started to develop the code, con-
tributed code snippets, and reported or identified bugs. The
resulting code changes range from the adaptation of more
modern coding standards, parallelization and efficiency en-
hancements, the improvement of the model functionality, and

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Hypothetical Grímsvötn eruption on 1 April 2015 at 00:00 UTC (instantaneous release). Total column concentrations are shown (µg m⁻²) 18 h after the eruption.</figcaption></figure>
<table><caption>**Table 14.** List of the test cases for FLEXPART 10.4.</caption><thead><tr><th>Group</th><th>Test name</th><th>Change to default options</th><th>Description</th></tr></thead><tbody><tr><td rowspan="6">(1) Gridded output</td><td>1</td><td>default (IOUT=1)</td><td>default option: forward run with mass concentration output; see Table 12</td></tr><tr><td>2</td><td>IOUT=2</td><td>mixing ratio</td></tr><tr><td>3</td><td>IOUT=3</td><td>concentration and mixing ratio</td></tr><tr><td>5</td><td>IOUT=5</td><td>concentration and trajectory cluster</td></tr><tr><td>9</td><td>IOUT=9</td><td>mass concentration in NetCDF output format</td></tr><tr><td>10</td><td>IOUT=10</td><td>mixing ratio in NetCDF output format</td></tr><tr><td>11</td><td>IOUT=11</td><td>concentration and mixing ratio in NetCDF output format</td></tr><tr><td rowspan="3">(2) Backward</td><td>bwd</td><td>LDIRECT=-1</td><td>SRR</td></tr><tr><td>bwd5</td><td>LDIRECT=-1, IOUT=5</td><td>backward trajectory cluster</td></tr><tr><td>bwd_nc</td><td>LDIRECT=-1, IOUT=9</td><td>SRR in NetCDF output format</td></tr><tr><td rowspan="3">(3) Particles</td><td>part1</td><td>IPOUT=1</td><td>particle dump</td></tr><tr><td>part2</td><td>IPOUT=2</td><td>particle dump at end of simulation</td></tr><tr><td>part_bwd1</td><td>IPOUT=1, LDIRECT=-1</td><td>backward trajectories</td></tr><tr><td rowspan="6">(4) Units</td><td>ind_1_2</td><td>IND_RECEPTOR=2</td><td>receptor (gridded) in mass mixing ratio units</td></tr><tr><td>ind_2_1</td><td>IND_SOURCE=2</td><td>source in mass mixing ratio units</td></tr><tr><td>ind_2_2</td><td>IND_SOURCE=2, IND_RECEPTOR=2</td><td>source and receptor in mass mixing ratio units</td></tr><tr><td>bwd_ind_1_2</td><td>IND_RECEPTOR=2, LDIRECT=-1</td><td>receptor in mass mixing ratio units, backward</td></tr><tr><td>bwd_ind_2_1</td><td>IND_SOURCE=2, LDIRECT=-1</td><td>source in mass mixing ratio units, backward</td></tr><tr><td>bwd_ind_2_2</td><td>IND_SOURCE=2, IND_RECEPTOR=2, LDIRECT=-1</td><td>source and receptor in mixing ratio units, backward</td></tr><tr><td rowspan="5">(5) Species</td><td>specNO</td><td>SPECNUM_REL=3</td><td>nitric oxide species</td></tr><tr><td>specCO</td><td>SPECNUM_REL=22</td><td>carbon monoxide species</td></tr><tr><td>specAERO-TRACE</td><td>SPECNUM_REL=23</td><td>idealized aerosol simulation</td></tr><tr><td>specBC</td><td>SPECNUM_REL=40</td><td>black carbon simulation</td></tr><tr><td>nested</td><td>NESTED_OUTPUT=1</td><td>nested output</td></tr><tr><td rowspan="6">(6) Nests</td><td>nested mr</td><td>NESTED_OUTPUT=1, IOUT=2</td><td>volume mixing ratio nested output</td></tr><tr><td>nested_bwd</td><td>NESTED_OUTPUT=1, LDIRECT=-1</td><td>nested output backwards</td></tr><tr><td>nested nc</td><td>NESTED_OUTPUT=1, IOUT=9</td><td>nested output in NetCDF format</td></tr><tr><td>nested mr_nc</td><td>NESTED_OUTPUT=1, IOUT=10</td><td>volume mixing ratio nested output NetCDF</td></tr><tr><td>nested_bwd_nc</td><td>NESTED_OUTPUT=1, LDIRECT=-1, IOUT=9</td><td>nested output backwards in NetCDF output format</td></tr><tr><td>(7) DOMAINEFILL</td><td>MDOMAINEFILL=1</td><td>domain-filling run</td></tr><tr><td rowspan="3">(8) Init. cond.</td><td>init_cond</td><td>LINIT_COND=1, LDIRECT=-1</td><td>sensitivity to initial conditions</td></tr><tr><td>init_cond_ind_1_2</td><td>LINIT_COND=1, LDIRECT=-1, IND_RECEPTOR=2</td><td>sensitivity to initial conditions in mixing ratio</td></tr><tr><td>init_cond_surr</td><td>LINIT_COND=1, SURF_ONLY=1, LDIRECT=-1</td><td>sensitivity to initial conditions adapted to surface fluxes</td></tr><tr><td rowspan="2">(9) CBL</td><td>CBLFLAG</td><td>CBLFLAG=1, CTL=40, IFINE=5</td><td>skewed turbulence</td></tr><tr><td>CBLFLAG_bwd</td><td>CBLFLAG=1, LDIRECT=-1, CTL=40, IFINE=5</td><td>skewed turbulence backwards</td></tr><tr><td rowspan="2">(10) B.d.</td><td>bwd_ind_1_3</td><td>IND_RECEPTOR=3, LDIRECT=-1</td><td>backward wet deposition</td></tr><tr><td>bwd_ind_1_4</td><td>IND_RECEPTOR=4, LDIRECT=-1</td><td>backward dry deposition</td></tr><tr><td>(11) Volc</td><td>Modified RELEASES file</td><td>hypothetical volcanic eruption</td></tr></tbody></table>
the addition of output options, to revisions and extensions of
the model physics. All this has been documented in this pa-
per. Integration of all these changes into a single stable model
version represents a growing challenge in itself, and efforts
to address this challenge (e.g., model website and repository,
version control, testing environment) have also been docu-
mented here.

As FLEXPART is developed further, updates will con-
tinue to be made available on the FLEXPART website at
https://flexpart.eu. We encourage established and new users
to contribute to FLEXPART development by providing their
code changes, as well as a description of these changes, as
new feature branches of the latest commits in the FLEX-
PART git repository. New code should pass all test cases pro-
vided in the FLEXPART distribution and provide consistent
output, unless there are specific reasons why output should be
different, such as improvements in the model physics. This
will expedite the integration of important new model features
in the main development branch of the model.

*Code and data availability.* The code described in this work is archived as flexpart10.4.tar at https://doi.org/10.5281/zenodo.3542278 (Pisso et al., 2019). FLEXPART downloads are available at https://www.flexpart.eu/downloads (last access: 22 November 2019) and FLEXPART releases are available at https://gitnilu.no/flexpart/flexpart/-/releases (last access: 22 November 2019). The working git repository for this version 10.4 (branches master, dev and release) can be accessed at https://gitnilu.no/flexpart/flexpart (last access: 22 November 2019) as well as the mirrors https://www.flexpart.eu/gitmob/flexpart (last access: 22 November 2019), https://github.com/flexpart/flexpart (last access: 22 November 2019) and https://bitbucket.org/flexpart/flexpart (last access: 22 November 2019).
## Appendix A: Installing FLEXPART and flex_extract

Here, we provide step-by-step instructions on how to install FLEXPART on Linux from scratch. This has been tested on an Ubuntu 16.4 distribution running on a dedicated instance in the Amazon cloud. Notice that in most environments, some or all of the required libraries (e.g., a Fortran compiler) are already installed and an installation totally from scratch would thus not be needed. In such cases, we strongly recommend that these libraries are used instead of installing everything from scratch. However, sometimes it may be necessary to install them from source (e.g., to avoid incompatibilities between different compilers or different versions of the same compiler). In the following, we assume that the user has root privileges in the system, but it is also possible for normal users to install the libraries in nonstandard locations. It is possible to ask for help by writing to the FLEXPART user email list (registration needed) or by creating a ticket on the community website at https://flexpart.eu.

### A1 System requirements

To begin, ensure that the latest packages are being used. This section is for completeness only, and most users (if not starting from a new system installation) can skip it and jump to Sect. A2.

sudo apt-get update

FLEXPART is developed using gfortran.

sudo apt-get install g++ gfortran

Some libraries (e.g., `grib_api`, `jasper-1.900.1`) require the GNU autotools suite in order to configure, build and install.

sudo apt-get install autoconf libtool automake flex bison

Newer packages (e.g., `ecCodes`) use `CMake` instead.

sudo apt-get install cmake

Python is not required for FLEXPART itself but is necessary for some preprocessing and post-processing tools, in particular `flex_extract` for retrieving ECMWF wind fields. Git is recommended to access the code repositories. An editor (e.g., `vim`) is usually also necessary.

sudo apt-get install python-dev python-pip git-core vim

### A2 Installing GRIB libraries

If JPG compression is needed to decode the input meteorological winds, download the `jasper` library from the `jasper` project page (http://www.ece.uvic.ca/~mdadams/jasper/, last access: 30 October 2019) and install it.

curl https://www.ece.uvic.ca/~frodo/jasper/

software/jasper-1.900.1.zip
--output jasper-1.900.1.zip

sudo apt install unzip
unzip jasper-1.900.1.zip

cd jasper-1.900.1
./configure

make
make check
sudo make install

Download the `grib_api` library from the ECMWF website (https://software.ecmwf.int/wiki/display/GRIB/Home, last access: 30 October 2019) and install it.

gunzip grib_api-X.X.X.tar.gz

tar -xf grib_api-X.X.X.tar

./configure [--prefix=grib_api_dir]
[-with_jasper=<jasper installation path>]

make
make check
make install

If you have no root privileges in your system, give the full path of `grib_api_dir` to the prefix option. If `jasper` is in a nonstandard location, it has to be passed to the `grib_api` configuration script. Please note that `GRIB-API` is no longer maintained. The primary `GRIB` encoding-decoding package used at the ECMWF is currently `ecCodes`. Any new features for the handling of `GRIB` files will only be developed in `ecCodes`. However, for FLEXPART v10.4 `grib_api` is sufficient. We keep the `grib_api` instructions for backward consistency.

For future versions, ensure that the path `/usr/local/lib/` is in the environment variable `$PATH`; otherwise, `ecCodes` may not find it. Obtain and unpack `ecCodes`.

curl https://software.ecmwf.int/wiki/download/attachments/45757960/eccodes-2.7.3-Source.tar.gz \
--output eccodes-2.7.3-Source.tar.gz
tar -xvf eccodes-2.7.3-Source.tar.gz

The `ecCodes` requires `CMake`. The installation procedure is described on the ECMWF `ecCodes` web page.

### A3 Installing NetCDF libraries

NetCDF output is optional. In order to enable NetCDF output, the NetCDF library has to be available in the system. For building the NetCDF library it is recommended to first build `HDF5` with support for compression (`zlib`). For this, download `zlib` (version 1.2.8) from the `zlib` website (https://www.zlib.net/, last access: 30 October 2019) and install it.
```
tar -xzvf zlib-1.2.8.tar.gz
cd zlib-1.2.8/
./configure [--prefix=<installation path>]
```

```
make
make install
```

Download HDF5 from the HDF Group website (https://www.hdfgroup.org/HDF5/release/obtainsrc.html, last access: 30 October 2019) and install it.

```
tar -xzvf HDF5-1.8.17.tar.gz
cd HDF5-1.8.17/
./configure --with-zlib=<path to zlib>
 [--prefix=<installation path>]
```

```
make
make check
make install
```

Download the latest stable version of NetCDF-C from the Unidata website (https://www.unidata.ucar.edu/downloads/netcdf/, last access: 30 October 2019) and install it.

```
tar -xzvf netcdf-4.4.1.tar.gz
cd netcdf-4.4.1/
./configure --enable-netcdf-4
 [--prefix=<installation path>]
```

```
make
make check
make install
```

Download the latest stable version of NetCDF-Fortran from the Unidata website (https://www.unidata.ucar.edu/downloads/netcdf/, last access: 30 October 2019) and install it.

```
tar -xzvf netcdf-fortran-4.4.4.tar.gz
cd netcdf-fortran-4.4.4/
./configure [--prefix=<installation path>]
```

```
make
make check
make install
```

## A4 Installing FLEXPART

Download the latest release of the FLEXPART source from the FLEXPART community website (https://flexpart.eu/wiki/FpDownloads, last access: 14 November 2019) or from https://transport-nilu.no (last access: 14 November 2019) and untar it.

```
tar -xvf flexpart10.4.tar.gz
```

Alternatively, clone the FLEXPART repository directly from the FLEXPART community site git.

```
git clone https://www.flexpart.eu/gitmob/
flexpart
```

This mirrors https://git Nilu.no/flexpart/flexpart. Additional mirrors exist, e.g., at Bitbucket (https://bitbucket.org/flexpart/flexpart, last access 30 October 2019) and GitHub (https://github.com/flexpart/flexpart, last access 30 October 2019). Edit the library path variable in the makefile according to the position of libeccodes (or libgrib_api) and libjasper. Optionally, edit the file par_mod.f90 to set parameters for the meteorological data, grid dimension and maximum particle number (maxpart, maxspec, nxmax, nylimax, nuvzmax, nwzmax, nzmax, nxshift). The default values are set to work with the test cases in Sect. 7 but may be too small for large simulations or too large for the available system resources. Then type

```
make
```

in order to create the executable. Invoking the executable FLEXPART should now print in the standard output.

Welcome to FLEXPART Version 10.4
FLEXPART is free software released under
the GNU General Public License.

However, without access to valid input data, the program will issue an error. Appendix C explains how to generate valid output with the standard meteorological fields from the ECMWF that can be obtained following the procedure described in Sect. A5. The makefile also allows the following command.

```
make clean
```

This can be used to safely remove all object and module files, e.g., if one wants to recompile after compiler option changes.

## A5 Installing flex_extract

A short description of the installation steps for this software is given for the public user mode (other modes are described in the flex_extract documentation). For this mode, the user does not need to be a member state user (https://www.ecmwf.int/en/about/who-we-are/member-states, last access: 13 October 2018) but can simply register at the ECMWF website. For the other operating modes and a more detailed explanation, see the README.md file of the python directory in the flex_extract distribution or the documentation files SIP.pdf and SUT_ondemand.pdf.

First of all, the user should register at the ECMWF website (https://www.ecmwf.int/en/forecasts/accessing-forecasts/order-historical-datasets, last access: 13 October 2018). To access public datasets each dataset license has to be
accepted separately before the account can be used for
retrieval of these data. This can be done at the following
website: https://software.ecmwf.int/wiki/display/WEBAPI/
Available+ECMWF+Public+Datasets (last access: 13 Octo-
ber 2018).

A5.1 System preparation for flex_extract

`flex_extract` requires a Python environment and a Fortran compiler. See Sect. A1 for installation instructions. To prepare the environment for the `flex_extract` installation, it is advisable to consider the official documentation and information from the ECMWF websites. We recommend the following steps.

1. For important information read the Emoslib (https://software.ecmwf.int/wiki/display/EMOS/Emoslib, last access: 13 October 2018) installation instructions first.

2. Read the ECMWF blog about gfortran (https://software.ecmwf.int/wiki/display/SUP/2015/05/11/Building+ECMWF+software+with+gfortran, last access: 13 October 2018) for details on the installation process of the libraries.

3. Install FFTW (http://www.fftw.org, last access: 13 October 2018) for Fortran, which is a library for computing the discrete Fourier transformation. This library is necessary for Emoslib. (Note: apply make twice! Once without any options and once with single precision option; see the information on the Emoslib website).

4. Install the interpolation library `Emoslib` for Fortran.

5. Install `ecCodes` (https://software.ecmwf.int/wiki/display/ECC, last access: 13 October 2018) or `grib_api` (https://software.ecmwf.int/wiki/display/GRIB/Home, last access: 13 October 2018) (for Python and Fortran). The `grib_api` support will be discontinued at the end of 2018 but `ecCodes` is downward compatible with `grib_api`.

6. Install the ECMWF WebAPI (https://confluence.ecmwf.int/display/WEBAPI/Access+MARS, last access: 13 October 2018) client by following the instructions on the website. It is a Python library to provide external access to the ECMWF servers.

7. Check whether `LD_LIBRARY_PATH` and `PATH` environment variables contain all paths to the previously installed libraries. The user should modify the `.bashrc` or `.tcshrc` file to guarantee that the variables contain the paths every time a new console is used.

8. Install the python package ` numpy ` via pip (https://scipy.org/install.html, last access: 13 October 2018).

9. Check the availability of Python packages (e.g., check the Python console for the following commands: `import eccodes`, `import grib_api`, `import ecmwfapi`)

10. Start a simple test retrieval (following the instructions on the ECMWF WebAPI website).

11. Install `flex_extract` (see the next section).

It is important to use the same compiler and compiler version for all libraries and the Fortran program CONVERT2.

A5.2 Building flex_extract

To install `flex_extract` a script `install.py` was prepared. The user can find it in the `python` directory of the `flex_extract` distribution.

The public user mode requires a local installation of `flex_extract`. Hence, we recommend adapting the paths to `ecCodes`, `Emoslib` or `grib_api` in one of the prepared makefiles, such as `Makefile.local.gfortran`, which can be found in the `src` directory. If a different compiler is used, this must also be adapted in the makefile. Then the installation script can be called as follows.

./install.py --target=local
--makefile=Makefile.local.gfortran

With this setting `flex_extract` is installed within the
current `flex_extract` directory. To install it in a different
place, e. g. within a FLEXPART distribution, the user can set
the path with the parameter `flexpart_root scripts`.
The installation was successful if the compilation of the For-
tran program (CONVERT2) did not fail and is displayed at
the end in the terminal.

A5.3 Running flex_extract

`flex_extract` is controlled by providing `CONTROL` files that contain a list of parameter settings. These parameters are described in detail in the Software User Tutorial (SUT.pdf) in the docs directory. The `CONTROL` files specify which ECMWF dataset is to be retrieved, the time and spatial resolution, the format of the GRIB file, and other options. In the Python directory are some example `CONTROL` files for the different datasets and access modes. They can be used as templates. `CONTROL` files with a `.public` ending are usable for the public access mode. The main difference is the parameter `dataset`, which explicitly specifies the public datasets. Note that not all meteorological fields, times and parameters were archived in the public datasets. This is already considered in the public `CONTROL` files.

To run `flex_extract`, the main program `submit.py` must be called. It retrieves the ECMWF data and generates the FLEXPART input files. To show all possible parameter options one can use the `-h` option. The script must be called
from the `python` directory of the `flex_extract` distribution. From the `-h` output it is clear that most parameters have default values or were already set via a `CONTROL` file parameter, except for the date. To retrieve just one day, one only needs to provide the start date. The rest will be done by `flex_extract`. This leads to the following script call for an arbitrary date.

```
./submit.py --controlfile=CONTROL_EI public \
  --start_date=20120101 \
  --public=1
```

The program now displays each MARS request and some messages for the preparation of the `FLEXPART` input files. Eventually, the program will finish with a `Done!` message if there was no error. Output will be stored in the default directory `work`, which is a subdirectory of the distribution directory (`flex_extract_v7.0.4`). The produced files can serve as input to `FLEXPART`.

## Appendix B: Running and testing FLEXPART

After a working `FLEXPART` executable is built (Appendix A), the next step is running the model and generating valid output. This requires consistent meteorological input data and user input files. In this section we describe the following: how to obtain the necessary wind fields (1), how to test run the executable with a default example (2), how to generate other examples (3), and how to run these examples and compare them with a reference output (4). In the following, `$flexhome` indicates the path to the root `FLEXPART` directory (e.g., `$HOME/flexpart/`) and `$flex_extracthome` indicates the path to the `flex_extract` root directory (e.g., `$flexhome/preprocess/flex_extract/`).

### B1 Meteorological input for the examples

Appendix A describes how to build the `flex_extract` version included in the source code. Here, we describe the settings to produce the meteorological input data required for running the default (Sect. B2) and derived (Sect. B3) cases. The instructions are for ECMWF ERA5 reanalysis, which is a publicly available dataset (https://confluence.ecmwf.int/display/WEBAPI/Access+ECMWF+Public+Datasets, last access: 30 October 2019). Therefore, the data can be obtained via `ecmwfapi` and no special access rights to the ECMWF are needed. However, in order to retrieve the data the user needs to register, obtain a personal Secure Shell (SSH) key and properly configure the file `.ecmwfapirc`. The execution of the retrieval requires the Python packages `ecmwfapi` (for access) and `grib_api` or `eccodes` (for processing). To retrieve the data, execute the following commands.

```
export PYTHONPATH=path/to/ecmwfapi:/path/to/
grib_api
$flex_extracthome/Python>./submit.py
  --start_date 20170102 --controlfile
  CONTROL_EA5
```

This should generate the files `EA170102??` in the following directory.

`$flex_extracthome/work/`

An `AVAILABLE` consistent with these wind fields is shipped together with the `FLEXPART` distribution.

`$flexhome/AVAILABLE`

### B2 Running the default example: installation verification

With the input files, which are included in the `FLEXPART` distribution and described in Sect. 5, a first test case to verify that `FLEXPART` was installed correctly can be run. To start the model run, the meteorological data have to be in `$flex_extracthome/work/` (see Sect. B1), the file `pathnames` in `$flexhome`, and the executable in `$flexhome/src/` in the `$flexhome` directory type.

`$flexhome>./src/FLEXPART`

The results created by this run are stored, e.g., in a directory `$flexhome/output` (as defined in `pathnames`). The run should end with the following message.

CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!

If this message is received, the model has completed the simulation, which confirms that `FLEXPART` and all required libraries are installed correctly. However, it does not guarantee valid output. To verify that the results obtained are valid, see Sect. B5.

### B3 Generating variations of the default example

To demonstrate more functionalities, a set of shell scripts generating different `FLEXPART` setups are provided in `$flexhome/tests/examples`. The script `set_default_example.sh` takes the content of the `options` directory and `pathnames` file from Sect. B2 as a basis, and then `gen_options_all.sh` creates new `options_suffix` directories for all of the cases described in Table 14. Here, the suffix corresponds to the example name as given in column 2 in Table 14. Finally, the script `gen_pathnames.sh` generates corresponding `pathnames_suffix` files pointing to all the `options suffix` directories. With this, all example cases in Table 14 are ready to run.
## B4 Running the examples

The examples can be run interactively one by one by invoking FLEXPART with the corresponding `pathnames_suffix` file. Alternatively, the script `gen_batch_jobs_cl.sh` generates a one-line script for each example case containing a call of FLEXPART and the appropriate `pathnames_suffix` file as a command-line parameter. All example scripts can then be run sequentially with `run_batch_cl.sh`, which creates `output_suffix` directories with the results, as well as log files `batch_job_pathnames_suffix.stdout` for each run. The examples described above can now be read and plotted with the tools included in the distribution. These steps are also automated in a `makefile`. All of the files and directories created by executing the scripts from Sects. B2 to B4 can be removed again with the command `make clean`.

## B5 Comparing the results

To verify that FLEXPART is producing valid output, it is useful to compare the output of a new installation with existing model output. It is also useful to repeat such a comparison after code changes to make sure the output is not affected, except for model simulations in which changes in the results are intended. While comprehensive comparisons of model results are possible, here we provide only a very simple way of checking the model results. The directory included in the FLEXPART distribution `$flexhome/tests/examples_reference/` contains the output of the examples described in Table 14. The file `read_examples_output.txt` contains, for the relevant examples that produce gridded output, the mean and the maximum value that occurs in the gridded output files. This shall serve as a reference to which users can compare their results and thus verify that the model produces output as expected. In addition, the directory `compare_examples` contains the script `compare_grids.sh` that allows for the partial automation of this task (output in `compare_examples_output.txt`).

## Appendix C: FLEXPART model versions

In addition to the reference version of FLEXPART described in this paper, there are many different model branches that were developed either for special purposes or to ingest other meteorological input data. This Appendix provides an incomplete list and a short description of some of these other versions. Further contributions are welcome in order to keep this list up to date.

## C1 FLEXPART-NorESM/CAM

Recently, the FLEXPART model version FLEXPART-NorESM/CAM was developed, which is tailored to run with

the meteorological output data generated by the CMIP5 version of NorESM1-M (the Norwegian Earth System Model) with 1.89° × 2.5° horizontal resolution and 26 vertical levels. The standard time resolution of the NorESM/CAM meteorological data is 3 h. FLEXPART-NorESM/CAM is based on FLEXPART v9, and the atmospheric component of NorESM1-M is based on CAM4 (the Community Atmosphere Model). The adaptation of FLEXPART to NorESM required new routines to read meteorological fields, new post-processing routines to obtain the vertical velocity in the FLEXPART coordinate system and other changes, as detailed by Cassiani et al. (2016). The code can be downloaded from https://www.flexpart.eu/wiki/FpClimateNorESM (last access: 30 October 2019).

## C2 FLEXPART-WRF

This FLEXPART version uses output from the Weather Research and Forecasting (WRF) mesoscale meteorological model (Brioude et al., 2013). Originally it was developed at the PNNL (Pacific Northwest National Laboratory) and named PILT (PNNL Integrated Lagrangian Transport). Compared to PILT, the further developed FLEXPART-WRF can use both instantaneous and time-averaged meteorological output of the WRF model. The latest version also includes the skewed turbulence scheme that was subsequently ported to the standard FLEXPART version 10.4. FLEXPART-WRF output can either be in binary or Network Common Data Form (NetCDF) format, both of which have efficient data compression. FLEXPART-WRF also offers effective parallelization with OpenMP in shared memory and an MPI library in distributed memory. Released versions of the code can be downloaded from https://www.flexpart.eu/wiki/ (last access: 30 October 2019) or cloned from the open repository git@git Nilu.no: flexpart/flexpart-wrf.git.

## C3 FLEXPART-COSMO

In Europe several national weather services and research groups cooperate to develop and operate the non-hydrostatic limited-area atmospheric model COSMO (Consortium for Small-scale Modeling). At MeteoSwiss COSMO is operationally run with data assimilation on two grids with approximately 7 × 7 km² and 2 × 2 km² horizontal resolution centered over Switzerland. This enables the study of atmospheric transport over complex terrain on a long-term basis. To this end, we have developed a new version of FLEXPART that is offline coupled to COSMO output (FLEXPART-COSMO hereafter) and supports output from multiple COSMO nests. Particles are internally referenced against the native vertical coordinate system used in COSMO and not, as in standard FLEXPART, in a terrain-following z system. This eliminates the need for an additional interpolation step. A new flux accumulation scheme was introduced that removes the need for additional preprocessing of the input files. In addition to
the existing Emanuel-based convection parameterization, a convection parameterization based on the Tiedtke scheme, which is identical to the one implemented in COSMO itself, was introduced. A possibility for offline nesting of a FLEXPART-COSMO run into a FLEXPART-ECMWF run for backward simulations was developed that only requires minor modifications of the FLEXPART-ECMWF version and allows particles to leave the limited COSMO domain. The OpenMP shared memory parallelization to the model allows for asynchronous reading of input data. The code is available on request from dominik.brunner@empa.ch and stephan.henne@empa.ch.

## C4 FLEXPART-AROME

The Applications of Research to Operations at Mesoscale (AROME) numerical weather prediction model is run operationally by Météo-France at the mesoscale. AROME forecasts for Europe exist at a resolution ranging from 0.5 to 2.5 km. The standard time resolution of the AROME meteorological data is 1 h. Based on FLEXPART-WRF, a coupling between FLEXPART and AROME was developed at Laboratoire de l'Atmosphère et des Cyclones (LACy, a joint institute between CNRS, Météo-France and the University of Réunion Island) using AROME high-resolution (2.5 × 2.5 km²) forecasts over the southwest Indian Ocean. The FLEXPART-AROME branch (Verreyken et al., 2019b) simulates turbulent transport using the Thomson turbulent scheme (Thomson, 1987) already implemented by Lin et al. (2003) in the Stochastic Time-Inverted Lagrangian Transport (STILT) model. This method constrains mass transport between different turbulent regions to conserve mass locally for a passive well-mixed tracer. Turbulent kinetic energy profiles are taken directly from AROME model outputs. Such treatment of turbulent motion ensures consistency between the turbulence in the meteorological fields calculated by the NWP model and turbulence computed in the offline Lagrangian transport model. It has been noticed that the use of a dedicated ABL scheme such as Hanna in the FLEXPART model may generate inconsistency between the ABL turbulent domain and the resolved wind fields used to drive FLEXPART. Simulations using the Thomson scheme show a better representation of the turbulent mixing between boundary layer air and free tropospheric air.

## C5 TRACZILLA

This branch-off from FLEXPART version 5 was originally developed for studies of transport and mixing in the upper troposphere-lower stratosphere region (e.g., Legras et al., 2003; Pisso and Legras, 2008). The modifications from the FLEXPART advection scheme consist mainly of discarding the intermediate terrain-following coordinate system and performing a direct vertical interpolation of winds, linear in log pressure, from hybrid levels. The vertical velocities are computed by the FLEXPART preprocessor using a mass-conserving scheme in the hybrid ECMWF coordinates. Alternatively, the vertical velocities can be computed from the rates of diabatic heating from ECMWF winds. In addition to the reanalyses from the ECMWF, the current version can use MERRA (Modern-Era Retrospective analysis for Research and Applications) from NASA and JRA-55 (the Japanese 55-year Reanalysis) from the Japanese Meteorological Agency (JMA). The parallelization uses the OMP version of PGI. All arrays are allocated dynamically. The code can be obtained from https://github.com/bernard-legras/traczilla (last access: 30 October 2019).
**Author contributions.** IP coordinated the contributions to the paper and the code development since version 9, including I/O, updates to turbulent mixing, the implementation of the tests and the distributed version control. ES developed and wrote the description of the parallelized version of FLEXPART and led the assembling of the new code developments into the main model version 10.4. HG developed and tested the new wet deposition scheme for aerosols. NIK contributed to the new wet deposition scheme for aerosols by testing the new model version and coordinated the contributions to the first version of the paper. MC developed the optional new turbulence scheme and the NorESM version and contributed to the WRF version. SE developed and wrote the description of the backward deposition, performed the benchmark test case together with IP, and worked on ECMWF data retrieval and testing. DA and DM contributed to the CTBTO developments including the unified executable, the Vtables approach and testing environment. RLT developed the temporal variation and temperature dependence of the OH reaction. CDGZ developed the dust mobilization scheme around FLEXPART and performed testing of the new model version. NE tested the new model version for black carbon and radionuclide applications. HS implemented the namelist input file format and contributed to the implementation of the NetCDF output and GRIB input routines. LH developed versions 2.0–7.02 of the *flex_extract* retrieval software, in particular the codes for calculating the hybrid coordinate vertical velocity. He also supervised the most recent developments and wrote the description together with AP. SH and DB contributed to the implementation of the NetCDF output module. JB coordinated the development up to FLEXPART version 8.3 and contributed to post-processing the Python module. AF developed the new Phyton-based ECMWF data retrieval software. JB led the development of the WRF and AROME versions and contributed to the turbulence scheme. AP developed, maintained and wrote the description of the *flex_extract* retrieval routines and contributed to the testing environment. PS devised the community website http://flexpart.eu, takes care of tickets and the wiki, contributed to various parts of the code development, and contributed to editing the paper. AS developed the first version of the code in 1998 and supervised all recent developments, including the new settling parameterization for aerosols and time-averaged particle output, and worked on the writing and structuring of the paper.

**Competing interests.** The authors declare that they have no conflict of interest.

**Acknowledgements.** The work was performed at the Nordic Center of Excellence eSTICC. Pirmin Kaufmann and Martin Schraner (MeteoSwiss) are acknowledged for code reformatting from fixed Fortran 77 to free Fortran 90 format. We thank Mariëlle Mulder for comments on an early version of this paper. The resources for the numerical simulations were provided by UNINETT Sigma2 (the National Infrastructure for High Performance Computing and Data Storage in Norway) under projects NN9419K and NS9419K. Input wind fields were provided by the ECMWF.

**Financial support.** This research has been supported by NordForsk (the Nordic Center of Excellence eSTICC, grant no. 57001), the European Research Council (project COMTESSA, grant no. 670462), and the CTBTO (Enhancements of the FLEXPART software on a call-off basis).

**Review statement.** This paper was edited by Slimane Bekki and reviewed by two anonymous referees.

## References

Arnalds, O., Dagsson-Waldhauserova, P., and Olafsson, H.: The Icelandic volcanic aeolian environment: Processes and impacts – A review, Aeolian Res., 20, 176–195, https://doi.org/10.1016/j.aeolia.2016.01.004, 2016.

Arnold, D., Maurer, C., Wotawa, G., Draxler, R., Saito, K., and Seibert, P.: Influence of the meteorological input on the atmospheric transport modelling with FLEXPART of radionuclides from the Fukushima Daiichi nuclear accident, J. Environ. Radioactiv., 139, 212–225, https://doi.org/10.1016/j.jenvrad.2014.02.013, 2015.

Asman, W. A. H.: Parametrisation of below-cloud scavenging of highly soluble gases under convective conditions, Atmos. Environment, 29, 1359–1368, 1995.

Atkinson, R.: Gas-phase tropospheric chemistry of volatile organic compounds: 1. Alkanes and alkenes, J. Phys. Chem. Ref. Data, 26, 215–290, 1997.

Balluch, M., and Haynes, P.: Quantification of lower stratospheric mixing processes using aircraft data, J. Geophys. Res., 102, 23487–23504, 1997.

Bey, I., Jacob, D. J., Logan, J. A., and Yantosca, R. M.: Asian chemical outflow to the Pacific in spring: Origins, pathways, and budgets, J. Geophys. Res., 106, 23073–23095, https://doi.org/10.1029/2001jd000806, 2001.

Brioude, J., Arnold, D., Stohl, A., Cassiani, M., Morton, D., Seibert, P., Angevine, W., Evan, S., Dingwell, A., Fast, J. D., Easter, R. C., Pisso, I., Burkhart, J., and Wotawa, G.: The Lagrangian particle dispersion model FLEXPART-WRF version 3.1, Geosci. Model Dev., 6, 1889–1904, https://doi.org/10.5194/gmd-6-1889-2013, 2013.

Bullard, J. E., Baddock, M., Bradwell, T., Crusius, J., Darlington, E., Gaiero, D., Gasso, S., Gisladottir, G., Hodgkins, R., McCulloch, R., Neuman, C. M., Mockford, T., Stewart, H., and Thorsteinsson, T.: High latitude dust in the earth system, Rev. Geophys., 54, 447–485, https://doi.org/10.1002/2016RG000518, 2016.

Cassiani, M., Stohl, A., and Brioude, J.: Lagrangian stochastic modelling of dispersion in the convective boundary layer with skewed turbulence conditions and a vertical density gradient: Formulation and implementation in the FLEXPART Model, Bound.-Lay. Meteorol., 154, 367–390, https://doi.org/10.1007/s10546-014-9976-5, 2015.

Cassiani, M., Stohl, A., Olivié, D., Seland, Ø., Bethke, I., Pisso, I., and Iversen, T.: The offline Lagrangian particle model FLEXPART-NorESM/CAM (v1): model description and comparisons with the online NorESM transport scheme and with
the reference FLEXPART model, Geosci. Model Dev., 9, 4029–4048, https://doi.org/10.5194/gmd-9-4029-2016, 2016.

Cozic, J., Verheggen, B., Mertes, S., Connolly, P., Bower, K., Petzold, A., Baltensperger, U., and Weingartner, E.: Scavenging of black carbon in mixed phase clouds at the high alpine site Jungfraujoch, Atmos. Chem. Phys., 7, 1797–1807, https://doi.org/10.5194/acp-7-1797-2007, 2007.

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kallberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., and Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system. Q. J. Roy. Meteor. Soc., 137, 553–597, https://doi.org/10.1002/qj.828, 2011.

Eckhardt, S., Prata, A. J., Seibert, P., Stebel, K., and Stohl, A.: Estimation of the vertical profile of sulfur dioxide injection into the atmosphere by a volcanic eruption using satellite column measurements and inverse transport modeling, Atmos. Chem. Phys., 8, 3881–3897, https://doi.org/10.5194/acp-8-3881-2008, 2008.

Eckhardt, S., Cassiani, M., Evangelou, N., Sollum, E., Pisso, I., and Stohl, A.: Source-receptor matrix calculation for deposited mass with the Lagrangian particle dispersion model FLEXPART v10.2 in backward mode, Geosci. Model Dev., 10, 4605–4618, https://doi.org/10.5194/gmd-10-4605-2017, 2017.

ECMWF: User Guide to ECMWF Products 2.1. Meteorological Bulletin M3.2. Reading, UK, 1995.

Emanuel, K. A. and Živković-Rothman, M.: Development and evaluation of a convection scheme for use in climate models. J. Atmos. Sci., 56, 1766–1782, 1999.

Fang, X., Shao, M., Stohl, A., Zhang, Q., Zheng, J., Guo, H., Wang, C., Wang, M., Ou, J., Thompson, R. L., and Prinn, R. G.: Top-down estimates of benzene and toluene emissions in the Pearl River Delta and Hong Kong, China, Atmos. Chem. Phys., 16, 3369–3382, https://doi.org/10.5194/acp-16-3369-2016, 2016.

Flesch, T. K., Wilson, J. D., and Lee, E.: Backward-time Lagrangian stochastic dispersion models and their application to estimate gaseous emissions, J. Appl. Meteorol., 34, 1320–1333, 1995.

Forster, C., Wandinger, U., Wotawa, G., James, P., Mattis, I., Althausen, D., Simmonds, P., O'Doherty, S., Kleefeld, C., Jennings, S. G., Schneider, J., Trickl, T., Kreipl, S., Jäger, H., and Stohl, A.: Transport of boreal forest fire emissions from Canada to Europe, J. Geophys. Res., 106, 22887–22906, 2001.

Forster, C., Stohl, A., and Seibert, P.: Parameterization of Convective Transport in a Lagrangian Particle Dispersion Model and Its Evaluation, J. Appl. Meteorol. Clim., 46, 403–422, https://doi.org/10.1175/JAM2470.1, 2007.

Global Soil Data Task: Global soil data products CD-ROM contents (IGBP-DIS), Data Set, Oak Ridge Natl. Lab. Distrib. Active Arch. Cent., Oak Ridge, Tenn., https://doi.org/10.3334/ORNLDAAC/565, 2014.

Groot Zwaastink, C. D., Grythe, H., Skov, H., and Stohl, A.: Substantial contribution of northern high-latitude sources to mineral dust in the Arctic, J. Geophys. Res.-Atmos., 121, 13678–13697, https://doi.org/10.1002/2016JD025482, 2016.

Groot Zwaastink, C. D., Arnalds, Ó., Dagsson-Waldhauserova, P., Eckhardt, S., Prospero, J. M., and Stohl, A.: Temporal and spatial variability of Icelandic dust emissions and atmospheric transport, Atmos. Chem. Phys., 17, 10865–10878, https://doi.org/10.5194/acp-17-10865-2017, 2017.

Groot Zwaastink, C. D., Henne, S., Thompson, R. L., Dlugokencky, E. J., Machida, T., Paris, J.-D., Sasakawa, M., Segers, A., Sweeney, C., and Stohl, A.: Three-dimensional methane distribution simulated with FLEXPART 8-CTM-1.1 constrained with observation data, Geosci. Model Dev., 11, 4469–4487, https://doi.org/10.5194/gmd-11-4469-2018, 2018.

Grythe, H., Kristiansen, N. I., Groot Zwaastink, C. D., Eckhardt, S., Ström, J., Tunved, P., Krejci, R., and Stohl, A.: A new aerosol wet removal scheme for the Lagrangian particle model FLEXPART v10, Geosci. Model Dev., 10, 1447–1466, https://doi.org/10.5194/gmd-10-1447-2017, 2017.

Haynes, P. and Anglade, J.: The Vertical-Scale Cascade in Atmospheric Tracers due to Large-Scale Differential Advection, J. Atmos. Sci., 54, 1121–1136, 1997.

Heinz, S.: Statistical Mechanics of Turbulent Flows, Springer, Berlin, Heidelberg, Germany, 214 pp., 2003.

Henne, S., Brunner, D., Oney, B., Leuenberger, M., Eugster, W., Bamberger, I., Meinhardt, F., Steinbacher, M., and Emmengger, L.: Validation of the Swiss methane emission inventory by atmospheric observations and inverse modelling, Atmos. Chem. Phys., 16, 3683–3710, https://doi.org/10.5194/acp-16-3683-2016, 2016.

Henning, S., Bojinski, S., Diehl, K., Ghan, S., Nyeki, S., Weingartner, E., Wurzler, S., and Baltensperger, U.: Aerosol partitioning in natural mixed-phase clouds, Geophys. Res. Lett., 31, L06101, https://doi.org/10.1029/2003GL019025, 2004.

Hertel, O., Christensen, J. Runge, E. H., Asman, W. A. H., Berkowicz, R., Hovmand, M. F., and Hov, O.: Development and testing of a new variable scale air pollution model – ACDEP, Atmos. Environ., 29, 1267–1290, 1995.

Hittmeir, S., Philipp, A., and Seibert, P.: A conservative reconstruction scheme for the interpolation of extensive quantities in the Lagrangian particle dispersion model FLEXPART, Geosci. Model Dev., 11, 2503–2523, https://doi.org/10.5194/gmd-11-2503-2018, 2018.

Hoinka, K: The tropopause: discovery, definition and demarcation, Meteorol. Z., 6, 281–303, 1997.

Hoyer, S. and Hamman, J.: xarray: N-D labeled Arrays and Datasets in Python, Journal of Open Research Software, 5, p. 10, https://doi.org/http://doi.org/10.5334/jors.148, 2017.

Jones, A., Thomson, D., Hort, M., and Devenish, B.: The UK Met Office's next generation atmospheric dispersion model, NAME III, in Air Pollution Modelling and its Application XVII, edited by: Borrego, C. and Norman, A.-L., Springer, New York, 580–589, 2007.

Kok, J. F.: A scaling theory for the size distribution of emitted dust aerosols suggests climate models underestimate the size of the global dust cycle, P. Natl. Acad. Sci. USA, 108, 1016–1021, https://doi.org/10.1073/pnas.1014798108, 2011.

Kristiansen, N. I., Stohl, A., Prata, A. J., Richter, A., Eckhardt, S., Seibert, P., Hoffmann, A., Ritter, C., Bitar, L., Duck, T. J., and Stebel, K.: Remote sensing and inverse transport modelling of the Kasatochi eruption sulphur dioxide cloud, J. Geophys. Res., 115, D00L16, https://doi.org/10.1029/2009JD013286, 2010.