Index: documentation/flexpart9.2.tex
===================================================================
--- documentation/flexpart9.2.tex	(revision 6b2046dd26c178ee4ab9df46bd4f25d15464b14c)
+++ documentation/flexpart9.2.tex	(revision 6b2046dd26c178ee4ab9df46bd4f25d15464b14c)
@@ -0,0 +1,2835 @@
+%% Version of source file:
+%% Date: 2003 July
+%%  **************************************
+%   *  TEMPLATE FOR EGU STYLE        *
+%   * FOR ARTICLES IN JOURNALS AND BOOKS *
+%%  **************************************
+%  Various alternatives for the input are shown, commented out.
+%  E.g., for the documentstyle options
+%        for the bibliography
+%  Feel free to play around with these variations, especially with
+%    the style options ms and twocolumn and 11pt/12pt
+
+ %%%% LATEX2E: SELECT ONE OF THE NEXT LINES
+\documentclass{egu}            % MANUSCRIPT, 10PT TYPE SIZE
+%\documentclass[ms,11pt]{egu}       % MANUSCRIPT, 11PT TYPE SIZE
+%\documentclass{egu}                 % for EGU TWO-COLUMN REVISED COPY
+%% <<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+%%   REMOVE THIS LINE ONLY IF IT CAUSES PROBLEMS
+  \usepackage{times}                   % WITH TIMES ROMAN FONT
+% ^^^^^^^^^^^^^^^^^^
+%%   This is STRONGLY recommended for the revised version!!
+%% <<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+
+%% ADD THIS PACKAGE IF GRAPHICS ARE TO BE IMPORTED
+\usepackage{graphicx} 
+\usepackage{amsmath} 
+%\usepackage{psfig}
+
+\printfigures              % PRINTS OUT FIGURES AT END OF MANUSCRIPT 
+\newcommand{\degreee}{{$\, \rm ^{\circ}$}}
+\newcommand{\degreen}{{$\, \rm ^{\circ}$~}}
+
+\begin{document} 
+
+\title{The Lagrangian particle dispersion model FLEXPART version 9.2}
+
+\author[1]{A. Stohl}
+\author[1]{H. Sodemann}
+\author[1]{S. Eckhardt}
+\author[2]{A. Frank}
+\author[2]{P. Seibert}
+\author[3]{G. Wotawa}
+
+\affil[1]{Norwegian Institute of Air Research, Kjeller, Norway}
+\affil[2]{Institute of Meteorology, University of Natural Resources and Applied Life Sciences, Vienna, Austria}
+\affil[3]{Preparatory Commission for the Comprehensive Nuclear Test Ban Treaty Organization, Vienna, Austria}
+
+\runningtitle{FLEXPART description}
+\runningauthor{Stohl et al.}
+%\runninghead{} % This may also be used instead of \runningtitle and \runningauthor
+
+\correspondence{A. Stohl (ast@nilu.no)}
+
+
+%\journal{ACP}       % for NONLINEAR PROCESSES IN GEOPHYSICS
+
+\date{Manuscript version from 29 November 2010}
+
+% ADDITIONAL NON-STANDARD COMMANDS FOR TITLE BLOCK INFORMATON
+\firstauthor{Stohl}
+%\proofs{A. Stohl\\
+%Norwegian Institute for Air Research\\
+%Instituttveien 18, 2027 Kjeller, Norway}
+%\offsets{A. Stohl\\
+%Norwegian Institute for Air Research\\
+%Instituttveien 18, 2027 Kjeller, Norway}
+
+% The manuscript number is supplied by the editorial office
+\msnumber{12345}
+
+% The following commands will be activated by the EGU editorial/production office
+% after inserting appropriate values.
+
+
+\maketitle % YOU MUST USE THE \maketitle COMMAND
+
+\begin{abstract}
+The Lagrangian particle dispersion model FLEXPART was originally (in its first
+release in 1998) designed for calculating the long-range and mesoscale
+dispersion of air pollutants from point sources, such as after an accident in a
+nuclear power plant.  In the meantime FLEXPART has evolved into a comprehensive
+tool for atmospheric transport modeling and analysis.  Its application fields
+were extended from air pollution studies to other topics where atmospheric
+transport plays a role (e.g., exchange between the stratosphere and
+troposphere, or the global water cycle).  It has evolved into a true community
+model that is now being used by at least 35 groups from 14 different countries
+and is seeing both operational and research applications.  The last citable
+manuscript for FLEXPART is: \citep{stohl2005}
+\end{abstract}
+
+\section{Updates since FLEXPART version 8.0} 
+
+In this version 8.2 of FLEXPART the representation of physical processes was
+improved as well as a number of technical changes and bugfixes implemented.  In
+addition, the program is now released under a GNU GPL license.  For the first
+time, detailed installation instructions are provided to help users getting
+started with running FLEXPART.  A short section on the new Python routines
+pflexpart for reading FLEXPART output data has been included as well.
+\bigskip
+
+\noindent {\bf Technical Changes:}
+
+\begin{itemize} 
+
+\item grib2 compatibility for ECMWF data which will be introduced soon (new
+data retrieval routines available)
+
+\item AVAILABLE files can now contain up to 256 characters, and include the
+path directly with the input file name.  This is used to gather data files
+stored in different directories.  The third line in the \verb|pathnames| file
+should be left empty if long data file names are used.
+
+\item The code was updated to work with the ECMWF grib\_api V1.6.1 and above
+(tested up to 1.9.5)
+
+\item An important bug the concentration output routines was fixed.  A
+numerical problem lead in some circumstances to garbled data.  The header
+version identifier has been changed to version 8.2.  New routines for loading
+data available for download on the FLEXPART homepage.
+
+\item Several bugs in the wet deposition parameterization were fixed.
+
+\item A bug which lead to ignoring landuses that was introduced in version 8.0
+has been fixed.
+
+\item Several bugs concerning nested grids have been fixed.
+
+\end{itemize}
+
+\noindent {\bf Algorithm Changes:}
+
+\begin{itemize}
+
+\item A new, more detailed settling parameterisation for aerosols was
+implemented.  The dynamic viscosity of air is now calculated as a function of
+temperature, and the calculation of settling velocities is now more accurate
+for higher Reynolds numbers.
+
+\item It is now possible to produce output files for backward runs that can be used
+to interface FLEXPART with model output from another model or from a FLEXPART
+forward run.  A gridded output file containing the exact sensitivities to these
+initial conditions can produced.  To this end, a new switch has been added in
+the COMMAND file.
+
+\end{itemize}
+
+\section{Introduction}
+
+Lagrangian particle models compute trajectories of a large number of so-called
+particles (not necessarily representing real particles, but infinitesimally
+small air parcels) to describe the transport and diffusion of tracers in the
+atmosphere.  The main advantage of Lagrangian models is that, unlike in
+Eulerian models, there is no numerical diffusion.  Furthermore, in Eulerian
+models a tracer released from a point source is instantaneously mixed within a
+grid box, whereas Lagrangian models are independent of a computational grid and
+have, in principle, infinitesimally small resolution.
+
+The basis for current atmospheric particle models was laid by
+\citet{thomson1987}, who stated the criteria that must be fulfilled in order
+for a model to be theoretically correct.  A monograph on the theory of
+stochastic Lagrangian models was published by \citet{rodean1996} and another
+good review was written by \citet{wilson1996}.  The theory of modeling
+dispersion backward in time with Lagrangian particle models was developed by
+\citet{flesch1995} and \citet{seibert2004}.  Reviews of the more practical
+aspects of particle modeling were provided by \citet{zannetti1992} and
+\citet{uliasz1994}.
+
+This note describes FLEXPART, a Lagrangian particle dispersion model that
+simulates the long-range and mesoscale transport, diffusion, dry and wet
+deposition, and radioactive decay of tracers released from point, line, area or
+volume sources.  It can also be used in a domain-filling mode where the entire
+atmosphere is represented by particles of equal mass.  FLEXPART can be used
+forward in time to simulate the dispersion of tracers from their sources, or
+backward in time to determine potential source contributions for given
+receptors.  The management of input data was largely taken from FLEXTRA, a
+kinematic trajectory model \citep{stohl1995}.  FLEXPART's first version was
+developed during the first author's military service at the
+nuclear-biological-chemical school of the Austrian Forces, the deposition code
+was added soon later (version 2), and this version was validated using data
+from three large tracer experiments \citep{stohletal1998}.  Version 3 saw
+performance optimizations and the development of a density correction
+\citep{stohlthomson1999}.  Further updates included the addition of a
+convection scheme \citep{seibertetal2001} (version 4), better backward
+calculation capabilities \citep{seibert2004}, and improvements in the
+input/output handling (version 5).  Validation was done during intercontinental
+air pollution transport studies \citep{stohl1999, forster2001, spichtinger2001,
+stohl2002, stohl2003}, for which also special developments for FLEXPART were
+made in order to extend the forecasting capabilities \citep{stohl2004}.
+Version 6.2 saw corrections to the numerics in the convection scheme, the
+addition of a domain-filling option, and the possibility to use output nests.
+Version 6.4 runs with NCEP GFS model data.  Version 7.0 was a transition
+version that was not publicly released.  Version 8.0 was a major release.  It
+unifies the ECMWF and GFS model versions in one source package. GFS data in
+GRIB2 format can be read, using ECMWF's grib\_api library. Each species
+got its own definition file.  Output can be written individually for multiple
+species in backward runs.  The output format was changed to a compressed sparse
+matrix format.  Memory is partly allocated dynamically.  Furthermore, dry and
+wet deposition algorithms were updated and new, global landuse inventory
+introduced.  OH reaction based on a monthly averaged 3 dimensional OH-field is
+available as an option.
+
+FLEXPART is coded following the Fortran 95 standard and tested with several
+compilers (gfortran, Absoft, Portland Group) under a number of operating
+systems (Linux, Solaris, Mac~OS~X, etc.).  The code is carefully documented and
+optimized for run-time performance.  No attempts have been made to parallelize
+the code because the model is strictly linear and, therefore, it is most
+effective to partition problems such that they run on single processors and to
+combine the results if needed.
+
+FLEXPART's source code and a manual are freely available from the internet page
+http://transport.nilu.no/flexpart.  According to a recent user survey, at least
+34 groups from 17 countries are currently using FLEXPART.  The user community
+maintains discussion by a mailing list to which one can subscribe on the
+flexpart home page.  The version of FLEXPART described here is based on model
+level data of the numerical weather prediction model of the European Centre for
+Medium-Range Weather Forecasts (ECMWF).  The standard source code distribution
+contains also the source files for a version of FLEXPART using the global
+National Centers of Environmental Prediction (NCEP) model data on pressure
+levels.  Other users have developed FLEXPART versions using input data from a
+suite of different and meso-scale (e.g.  MM5, WRF, COSMO) models, some of which
+are available from the FLEXPART website but are not described here.
+
+\section{License}
+
+FLEXPART has been free software ever since it first was released. The status as
+a free software is now more formally established by releasing the code under
+the GNU General Public License (GPL) Version 3.  The text of the license is
+included in the file COPYING in the source code archive.
+
+\section{Installation}\label{sec:installation}
+
+Getting FLEXPART up and running on their systems has been a major hurdle to many
+new users of the software. Since the inclusion of the grib\_api library has further
+complicated the compilation of FLEXPART, we include a step-by-step installation
+instruction of FLEXPART.  The steps below are described for an Ubunto 10.4 LT
+UNIX release.
+
+\subsection{Installing required libraries}
+
+In order to process input files in GRIB2 format, FLEXPART V8.2 needs to be
+compiled with ECMWF's grib\_api library version 1.6.1 and above.  Since the API
+of the grib\_api library may change in the future, upward compatibility cannot
+be guaranteed.  The input files in GRIB2 format can be compressed to save
+bandwitdth and storage space.  In order to make use of this feature, the jasper
+library needs to be installed on the system.  Optionally, the emos library can
+be used to run FLEXPART with input files in GRIB1 format.
+
+\noindent The following steps should be executed in sequence:
+
+\subsubsection{Install jasper library (Version 1.900.1)}
+
+Download the jasper library from the jasper project page\footnote{http://www.ece.uvic.ca/~mdadams/jasper/} 
+
+\begin{small}
+\begin{verbatim}
+unzip jasper-1.900.1.zip
+cd jasper-1.900.1
+./configure [--prefix=<installation path>]
+make
+make check
+make install
+\end{verbatim}
+\end{small}
+
+Parameters in brackets are optional but may require root privileges.  
+
+\subsubsection{Install grib\_api library (Version 1.6.1 or later)}
+  
+Download the grib\_api library from the ECMWF
+website\footnote{http://www.ecmwf.int/products/data/software/download/grib\_api.html}
+
+\begin{small}
+\begin{verbatim}
+tar -xvf grib_api-1.12.3.tar.gz
+./configure [--with-jasper=<jasper path>]
+make
+make check
+make install
+\end{verbatim}
+\end{small}
+
+\subsubsection{Optional: Install emos library (Version 000372)}
+
+Download the emos library from the ECMWF
+website\footnote{http://www.ecmwf.int/products/data/software/interpolation.html}
+
+\begin{small}
+\begin{verbatim}
+tar -xvf emos_000372.tar.gz
+./build_library
+./install
+\end{verbatim}
+\end{small}
+
+\subsection{Compiling FLEXPART V8.2}
+
+Download the FLEXPART source code archive from the FLEXPART
+homepage\footnote{http://transport.nilu.no/flexpart/flexpart/view}
+
+\begin{small}
+\begin{verbatim}
+tar -xvf flexpart82.tar.gz
+\end{verbatim}
+\end{small}
+optionally edit the file includepar to set parameters for the data center, grid
+dimension, and particle number edit the LIBRARY path variable in the makefiles
+according to the position of libgrib\_api and libjasper
+\begin{small}
+\begin{verbatim}
+make -f makefile.<center>_<compiler>_<system>
+\end{verbatim}
+\end{small}
+In the above statement, center can be one of: gfs, ecmwf, ecmwf\_emos,
+gfs\_emos compiler can be one of absoft or gfortran (emos library only with
+absoft) system can be one of 32, 64 (emos only 32).  The system parameter must
+match that of the compiled libraries. See also Table~\ref{tab:makefile} for all 
+available makefiles.
+
+When recompiling after making changes, all object and module files can be
+removed safely by using
+\begin{small}
+\begin{verbatim}
+make -f makefile.<xxx> clean
+\end{verbatim}
+\end{small}
+
+\begin{table*}
+\setlength{\tabcolsep}{1.1mm}
+\caption{\label{tab:makefile}
+List of available makefiles }
+\vspace{3mm}
+{\centerline{
+\begin{tabular}{llc} \hline
+Makefile GFS & Makefile ECMWF & GRIB version \\ 
+\hline
+makefile.gfs\_gfortran\_32 & makefile.ecmwf\_gfortran\_32 & 1/2 \\
+makefile.gfs\_gfortran\_64 & makefile.ecmwf\_gfortran\_64 & 1/2 \\
+makefile.gfs\_absoft\_32 & makefile.ecmwf\_absoft\_32 & 1/2 \\
+makefile.gfs\_absoft\_64 & makefile.ecmwf\_absoft\_64 & 1/2 \\
+makefile.gfs\_emos\_absoft\_32 & makefile.ecmwf\_emos\_absoft\_32 & 1 \\
+\hline
+\end{tabular}}}
+\end{table*}
+
+\section{Input data and grid definitions}
+
+FLEXPART is an off-line model that uses meteorological fields (analyses or
+forecasts) in Gridded Binary (GRIB) format in version 1 or 2 from the ECMWF
+numerical weather prediction model \citep{ecmwf1995} on a latitude/longitude
+grid and on native ECMWF model levels as input.  Optionally, GRIB data from 
+NCEP's GFS model, available on pressure levels, can be used.  The ECMWF data
+can be retrieved from the ECMWF archives using a pre-processor that is also
+available from the FLEXPART website but not described here.  The GRIB decoding
+libraries is {\it not} provided with the FLEXPART source codes but is publicly
+available (see Sec.~\ref{sec:installation}.  The data can be global or only
+cover a limited area.  Furthermore, higher-resolution domains can be nested
+into a mother domain.
+
+The file \verb|includepar| contains all relevant FLEXPART parameter settings,
+both physical constants and maximum field dimensions.  As the memory required
+by FLEXPART is determined by the various field dimensions, it is recommended
+that they are adjusted to actual needs before compilation.  The file
+\verb|includecom| defines all FLEXPART global variables and fields, i.e., those
+shared between most subroutines.
+
+\subsection{Input data organisation}
+
+A file \verb|pathnames| must exist in the directory where FLEXPART is started.
+It must contain at least four lines:\\ 1.  line: Directory where all the
+FLEXPART command files are stored.\\ 2.  line: Directory to which the model
+output is written.\\ 3.  line: Directory where the GRIB input fields are
+located.\\ 4.  line: Path name of the \verb|AVAILABLE| file (see below).\\ If
+nests with higher-resolution input data shall also be used, lines 3 and 4 must
+be repeated for every nest, thus also specifying the nesting level order.  Any
+number of nesting levels can be used up to a maximum (parameter
+\verb|maxnests|).
+
+The meteorological input data must be organised such that all data for a domain
+and a certain date must be contained in a single GRIB file.  The
+\verb|AVAILABLE| file lists all available dates and the corresponding file
+names.  For each nesting level, the input files must be stored in a different
+directory and the \verb|AVAILABLE| file must contain the same dates as for the
+mother domain.  Given a certain particle position, the last (i.e., innermost)
+nest is checked first whether it contains the particle or not.  If the particle
+resides in this nest, the meteorological data from this nest is interpolated
+linearly to the particle position.  If not, the next nest is checked, and so
+forth until the mother domain is reached.  There is no nesting in the vertical
+direction and the poles must not be contained in any nest.
+
+The maximum dimensions of the meteorological fields are specified by the
+parameters \verb|nxmax, nymax, nuvzmax, nwzmax, nzmax| in file
+\verb|includepar|, for x, y, and three z dimensions, respectively.  The three z
+dimensions are for the original ECMWF data (\verb|nuvzmax, nwzmax| for model
+half levels and model levels, respectively) and transformed data (\verb|nzmax|,
+see below), respectively.  The horizontal dimensions of the nests must be
+smaller than the parameters \verb|nxmaxn, nymaxn|.  Grid dimensions and other
+basic things are checked in routine \verb|gridcheck.f| (\verb|gridcheck_gfs.f|
+for the GFS version), and error messages are issued if necessary.
+
+The longitude/latitude range of the mother grid is also used as the
+computational domain.  All internal FLEXPART coordinates run from the
+western/southern domain boundary with coordinates (0,0) to the eastern/northern
+boundary with coordinates (nx-1,ny-1), where (nx,ny) are the mother grid
+dimensions.  For global input data, FLEXPART repeats the westernmost grid cells
+at the easternmost domain "boundary", in order to facilitate interpolation on
+all locations of the globe (e.g., if input data run from 0 to 359\degreen with
+1\degreen resolution, 0\degreen data are repeated at 360\degreee).  A global
+mother domain can be shifted by \verb|nxshift| (file \verb|includepar|) data
+columns (subroutines \verb|shift_field.f| and \verb|shift_field_0.f|) if nested
+input fields would otherwise overlap the "boundaries".  For instance, a domain
+stretching from 320\degreen to 30\degreen can be nested into the mother grid of
+the above example by shifting the mother grid by 30\degreee. The default setting
+for global ECMWF fields is \verb|nxshift=359|, while GFS fields the default
+value is \verb|nxshift=0|.
+
+\subsection{Vertical model structure and required data}
+
+FLEXPART needs five three-dimensional fields: horizontal and vertical wind
+components, temperature and specific humidity.  Input data must be on ECMWF
+model (i.e.  $\eta$) levels which are defined by a hybrid coordinate system.
+The conversion from $\eta$ to pressure coordinates is given by $p_k=A_k+B_kp_s$
+and the heights of the $\eta$ surfaces are defined by $\eta_k=A_k/p_0+B_k$,
+where $\eta_k$ is the value of $\eta$ at the $k^{\rm_{th}}$ model level, $p_s$
+is the surface pressure and $p_0=$101325 Pa.  $A_k$ and $B_k$ are coefficients,
+chosen such that the levels closest to the ground follow the topography, while
+the highest levels coincide with pressure surfaces; intermediate levels
+transition between the two.  The vertical wind in hybrid coordinates is
+calculated mass-consistently from spectral data by the pre-processor.  A
+surface level is defined in addition to the regular $\eta$ levels.  2~m
+temperature, 10~m winds and specific humidity from the first regular model
+level are assigned to this level, to represent "surface" values.
+
+Parameterized random velocities in the atmospheric boundary layer (ABL, see
+sections~\ref{PBLparameterization} and \ref{diffusion}) are calculated in units
+of m s$^{-1}$, and not in $\eta$ coordinates.  Therefore, in order to avoid
+time-consuming coordinate transformations every time step, all
+three-dimensional data are interpolated linearly from the ECMWF model levels to
+terrain-following Cartesian coordinates $\tilde{z}=z-z_t$, where $z_t$ is the
+height of the topography (subroutine \verb|verttransform.f|).  The conversion
+of vertical wind speeds from the eta coordinate system into the
+terrain-following co-ordinate system follows as
+
+\begin{equation}
+\tilde{w}=\dot{\tilde{z}}=\dot{\tilde{\eta}} \left( \frac{\partial p}{\partial z} \right)^{-1} + \left.\frac{\partial \tilde{z}}{\partial t}\right|_\eta +  \vec{v}_h \cdot \nabla_\eta \tilde{z}  \,
+\end{equation}
+
+where $\dot{\tilde{\eta}}=\dot{\eta} \partial p / \partial \eta$.  The second
+term on the right hand side is missing in the FLEXPART transformation because
+it is much smaller than the others.  One colleague has implemented this term in
+his version of FLEXPART and found virtually no differences (B.  Legras,
+personal communication).
+
+FLEXPART also needs the two-dimensional fields: surface pressure, total cloud
+cover, 10~m horizontal wind components, 2~m temperature and dew point
+temperature, large scale and convective precipitation, sensible heat flux,
+east/west and north/south surface stress, topography, land-sea-mask and subgrid
+standard deviation of topography.  The landuse inventory of \citet{belward1999}
+is provided in an extra file in the options directory (\verb|IGBP_int1.dat|).
+
+Note that GFS data is provided on pressure levels (currently 26) which requires
+that the FLEXPART code is compiled with a makefile for the GFS model. GFS data
+are freely available from the NCEP data archives. An example for a retrieval
+of GFS data files via ftp is given below:
+
+\begin{small}
+\begin{verbatim}
+cd $DIR_GFS_05TEMPORARY
+ftp ftpprd.ncep.noaa.gov << EOF11
+binary
+cd /pub/data/nccf/com/gfs/prod/gfs.2009111006
+get gfs.t06z.pgrb2f00 GF09111006
+quit
+EOF11
+\end{verbatim}
+\end{small}
+
+\subsection{Meteorological input data formats}
+
+While NCEP has already switched to the more flexible and compact GRIB2 format,
+the ECMWF is only gradually transitioning from GRIB1 to the new GRIB2 format.
+Transition to GRIB2 becomes necessary, at least for 3D fields, when the ECMWF
+increases the vertical resolution beyond 100 model levels in 2012 (which cannot
+be handled by the GRIB1 format).  Currently, GRIB2 codes have been defined for
+all 3D-fields required by FLEXPART, while definitions are still missing for a
+part of the 2D-fields.
+
+According to the ECMWF it may take another year to define the GRIB2 codes for
+all 2D-fields used by FLEXPART.  Therefore, in version 8.2, it is now possible
+to read in GRIB files that contain only GRIB1 coded fields, files that contain
+only GRIB2 coded fields, and files that contain mixed GRIB1/2 fields.  Such
+mixed GRIB1/2 files are produced by the current version 4 of the FLEXPART data
+retrieval routines (available from the FLEXPART hompage).
+
+\section[Physical parameterizations]{Physical parameterization of boundary layer parameters}
+
+\label{PBLparameterization} Accumulated surface sensible heat fluxes and
+surface stresses are available from ECMWF forecasts.  The pre-processor selects
+the shortest forecasts available for that date from the ECMWF archives and
+deaccumulates the flux data.  The total surface stress is computed from
+
+\begin{equation}
+\tau=\sqrt{\tau_1^2+\tau_2^2} \;,
+\end{equation}
+
+where $\tau_1$ and $\tau_2$ are the surface stresses in east/west and
+north/south direction, respectively.  Friction velocity is then calculated in
+subroutine \verb|scalev.f| as
+
+\begin{equation}
+u_*=\sqrt{\tau/\rho} \;,
+\end{equation}
+
+where $\rho$ is the air density \citep{wotawa1996}.  Friction velocities and
+heat fluxes calculated using this method are most accurate \citep{wotawa1997}.
+However, if deaccumulated surface stresses and surface sensible heat fluxes are
+not available, the profile method after \citet{berkowicz1982} (subroutine
+\verb|pbl_profile.f|) is applied to wind and temperature data at the second
+model level and at 10~m (for wind) and 2~m (for temperature) (note that
+previously the first model level was used; as ECMWF has its first model level
+now close to 10~m, the second level is used instead).  The following three
+equations are solved iteratively:
+
+\begin{equation}
+u_*=\frac{\kappa\Delta u}{\ln {\frac{z_l}{10}} -\Psi_m(\frac{z_l}{L}) +\Psi_m(\frac{10}{L})} \; ,
+\end{equation}
+
+\begin{equation}
+\Theta_*=\frac{\kappa\Delta \Theta}{0.74 \left[{\ln {\frac{z_l}{2}} -\Psi_h(\frac{z_l}{L}) +\Psi_h(\frac{2}{L})}\right]} \; ,
+\end{equation}
+
+\begin{equation}
+L= \frac{\overline{T}u_*^2}{g \kappa \Theta_*} \;,
+\end{equation}
+
+where $\kappa$ is the von K\'{a}rm\'{a}n constant (0.4), $z_l$ is the height of
+the second model level, $\Delta u$ is the difference between wind speed at the
+second model level and at 10~m, $\Delta \Theta$ is the difference between
+potential temperature at the second model level and at 2~m, $\Psi_m$ and
+$\Psi_h$ are the stability correction functions for momentum and heat
+\citep{businger1971, beljaars1991}, $g$ is the acceleration due to gravity,
+$\Theta_*$ is the temperature scale and $\overline{T}$ is the average surface
+layer temperature (taken as $T$ at the first model level).  The heat flux is
+then computed by
+
+\begin{equation}
+(\overline{w'\Theta_v'})_0=-\rho c_p u_* \Theta_* \;,
+\end{equation}
+
+where $\rho c_p$ is the specific heat capacity of air at constant pressure.
+
+ABL heights are calculated according to \citet{vogelezang1996} using the
+critical Richardson number concept (subroutine \verb|richardson.f|).  The ABL
+height $h_{mix}$ is set to the height of the first model level $l$ for which
+the Richardson number
+
+\begin{equation}
+Ri_l=\frac{(g/\Theta_{v1})(\Theta_{vl}-\Theta_{v1})(z_l-z_1)}{(u_l-u_1)^2+(v_l-v_1)^2+100u_*^2}\;,
+\label{richardson}
+\end{equation}
+
+exceeds the critical value of 0.25.  $\Theta_{v1}$ and $\Theta_{vl}$ are the
+virtual potential temperatures, $z_1$ and $z_l$ are the heights of, and
+$(u_1,v_1)$, and $(u_l,v_l)$ are the wind components at the 1$^{\rm st}$ and
+$l^{\rm th}$ model level, respectively.  The formulation of
+Eq.~\ref{richardson} can be improved for convective situations by replacing
+$\Theta_{v1}$ with
+
+\begin{equation}
+\Theta_{v1}'=\Theta_{v1}+8.5\frac{(\overline{w'\Theta_v'})_0}{w_* c_p}\;,
+\label{excess}
+\end{equation}
+
+where
+
+\begin{equation}
+w_*=\left[{ \frac{(\overline{w'\Theta_v'})_0 g h_{mix}}{\Theta_{v1} c_p} }\right]^{1/3}
+\end{equation}
+
+is the convective velocity scale.  The second term on the right hand side of
+Eq.~\ref{excess} represents a temperature excess of rising thermals.  As $w_*$
+is unknown beforehand, $h_{mix}$ and $w_*$ are calculated iteratively.
+
+Spatial and temporal variations of ABL heights on scales not resolved by the
+ECMWF model play an important role in determining the thickness of the layer
+over which tracer is effectively mixed.  The height of the convective ABL
+reaches its maximum value (say 1500~m) in the afternoon (say, at 1700 local
+time (LT)), before a much shallower stable ABL forms.  Now, if meteorological
+data are available only at 1200 and 1800 LT and the ABL heights at those times
+are, say, 1200~m and 200~m, and linear interpolation is used, the ABL height at
+1700 LT is significantly underestimated (370~m instead of 1500~m).  If tracer
+is released at the surface shortly before the breakdown of the convective ABL,
+this would lead to a serious overestimation of the surface concentrations (a
+factor of four in the above example).  Similar arguments hold for spatial
+variations of ABL heights due to complex topography and variability in landuse
+or soil wetness \citep{hubbe1997}.  The thickness of a tracer cloud traveling
+over such a patchy surface would be determined by the maximum rather than by
+the average ABL height.
+
+In FLEXPART a somewhat arbitrary parameterization is used to avoid a
+significant bias in the tracer cloud thickness and the surface tracer
+concentrations.  To account for spatial variations induced by topography, we
+use an "envelope" ABL height
+
+\begin{equation}
+H_{env}=h_{mix}+\min \left[\sigma_Z, c \frac{V}{N} \right]\; .
+\end{equation}
+
+Here, $\sigma_Z$ is the standard deviation of the ECMWF model subgrid
+topography, $c$ is a constant (here: 2.0), $V$ is the wind speed at height
+$h_{mix}$, and $N$ is the Brunt-Vaisala frequency.  Under convective
+conditions, the envelope ABL height is, thus, the diagnosed ABL height plus the
+subgrid topography (assuming that the ABL height over the hill tops effectively
+determines the dilution of a tracer cloud located in a convective ABL).  Under
+stable conditions, air tends to flow around topographic obstacles rather than
+above it, but some lifting is possible due to the available kinetic energy.
+$\frac{V}{N}$ is the local Froude number (i.e., the ratio of inertial to
+buoyant forces) times the length scale of the sub-grid topographic obstacle.
+The factor $c \frac{V}{N}$, thus, limits the effect of the subgrid topography
+under stable conditions, with $c=2$ being a subjective scaling factor.
+$H_{env}$ rather than $h_{mix}$ is used for all subsequent calculations.  In
+addition, $H_{env}$ is not interpolated to the particle position, but the
+maximum $H_{env}$ of the grid points surrounding a particle's position in space
+and time is used.
+
+\section{\label{diffusion}Particle transport and diffusion}
+
+\subsection{Particle trajectory calculations}
+
+FLEXPART generally uses the simple ``zero acceleration'' scheme 
+
+\begin{equation}
+{\vec{X}}(t+\Delta t)={\vec{X}}(t)+{\vec{v}}({\vec{X}},t) \Delta t\,,
+\label{firstorder}
+\end{equation}
+
+which is accurate to the first order, to integrate the trajectory equation \citep{stohl1998}
+
+\begin{equation}
+\frac{d {\vec{X}}}{dt}={\vec{v}}[{\vec{X}}(t)] \,,
+\end{equation}
+
+with $t$ being time, $\Delta t$ the time increment, ${\vec{X}}$ the position
+vector, and ${\vec{v}}=\overline{\vec{v}}+{\vec{v}}_t+{\vec{v}}_m$ the wind
+vector that is composed of the grid scale wind $\overline{\vec{v}}$, the
+turbulent wind fluctuations ${\vec{v}}_t$ and the mesoscale wind fluctuations
+${\vec{v}}_m$.
+
+Since FLEXPART version 5.0, numerical accuracy has been improved by making one
+iteration of the \citet{petterssen1940} scheme (which is accurate to the second
+order) whenever this is possible, but only for the grid-scale winds.  It is
+implemented as a correction applied to the position obtained with the ``zero
+acceleration'' scheme.  In three cases it cannot be applied.  First, the
+Petterssen scheme needs winds at a second time which may be outside the time
+interval of the two wind fields kept in memory.  Second, if a particle crosses
+the boundaries of nested domains, and third in the ABL if \verb|ctl|$>0$ (see
+below).
+
+Particle transport and turbulent dispersion are handled by the subroutine
+\verb|advance.f| where calls are issued to procedures that interpolate winds
+and other data to the particle position and the Langevin equations (see below)
+are solved.  The poles are singularities on a latitude/longitude grid.  Thus,
+horizontal winds (variables \verb|uu,vv|) poleward of latitudes
+(\verb|switchnorth, switchsouth|) are transformed to a polar stereographic
+projection (variables \verb|uupol,vvpol|) on which particle advection is
+calculated.  As \verb|uupol,vvpol| are also stored on the latitude/longitude
+grid, no additional interpolation is made.
+
+\subsection{The Langevin equation}
+
+Turbulent motions ${\vec{v}}_{t}$ for wind components $i$ are parameterized
+assuming a Markov process based on the Langevin equation \citep{thomson1987}
+
+\begin{equation}
+dv_{t_i}=a_i({\vec{x}},{\vec{v}}_t,t)dt+b_{ij}({\vec{x}},{\vec{v}}_t,t)dW_
+j \,,
+\label{langevin}
+\end{equation}
+
+where the drift term $a$ and the diffusion term $b$ are functions of the
+position, the turbulent velocity and time.  $dW_j$ are incremental components
+of a Wiener process with mean zero and variance $dt$, which are uncorrelated in
+time \citep{legg1982}.  Cross-correlations between the different wind
+components are also not taken into account, since they have little effect for
+long-range dispersion \citep{uliasz1994}.
+
+Gaussian turbulence is assumed in FLEXPART, which is strictly valid only for
+stable and neutral conditions.  Under convective conditions, when turbulence is
+skewed and larger areas are occupied by downdrafts than by updrafts, this
+assumption is violated, but for transport distances where particles are rather
+well mixed throughout the ABL, the error is minor.
+
+With the above assumptions, the Langevin equation for the vertical wind
+component $w$ can be written as
+
+\begin{multline}
+d w = \\
+- w \frac{dt}{\tau_{L_w}}
++ \frac{\partial \sigma_w^2}{\partial z} dt
++ \frac{\sigma_w^2}{\rho} \frac{\partial \rho}{\partial z} dt
++ \left( \frac{2}{\tau_{L_w}} \right)^{1/2} \sigma_w\; dW \;,
+\label{legg}
+\end{multline}
+
+where $w$ and $\sigma_w$ are the turbulent vertical wind component and its
+standard deviation, $\tau_{L_w}$ is the Lagrangian timescale for the vertical
+velocity autocorrelation and $\rho$ is density.  The second and the third term
+on the right hand side are the drift correction \citep{mcnider1988} and the
+density correction \citep{stohlthomson1999}, respectively.  This Langevin
+equation is identical to the one described by \citet{legg1982}, except for the
+term from \citet{stohlthomson1999} which accounts for the decrease of air
+density with height.
+
+Alternatively, the Langevin equation can be re-expressed in terms of
+$w/\sigma_w$ instead of $w$ \citep{wilson1983}:
+
+\begin{multline}
+d \left( \frac{w}{\sigma_w} \right) =\\
+- \frac{w}{\sigma_w} \frac{dt}{\tau_{L_w}}
++ \frac{\partial \sigma_w}{\partial z} dt
++ \frac{\sigma_w}{\rho} \frac{\partial \rho}{\partial z} dt
++ \left( \frac{2}{\tau_{L_w}} \right)^{1/2} dW \;,
+\label{wilson}
+\end{multline}
+
+This form was shown by \citet{thomson1987} to fulfill the well-mixed criterion
+which states that ``if a species of passive marked particles is initially mixed
+uniformly in position and velocity space in a turbulent flow, it will stay that
+way'' \citep{rodean1996}.  Although the method proposed by \citet{legg1982}
+violates this criterion in strongly inhomogeneous turbulence, their formulation
+was found to be practical, as numerical experiments have shown that it is more
+robust against an increase in the integration time step.  Therefore,
+Eq.~\ref{legg} is used with long time steps (see section~\ref{timestep});
+otherwise, Eq.~\ref{wilson} is used.  For the horizontal wind components, the
+Langevin equation is identical to Eq.~\ref{legg}, with no drift and density
+correction terms.
+
+For the discrete time step implementation of the above Langevin equations (at
+the example of Eq.~\ref{wilson}), two different methods are used.  When $\left
+(\Delta t/\tau_{L_w} \right) \ge 0.5$,
+
+\begin{multline}
+\left (\frac{w}{\sigma_w} \right)_{k+1}=
+r_w \left (\frac{w}{\sigma_w} \right)_k
++ \frac{\partial \sigma_w}{\partial z} \tau_{L_w} \left (1-r_w \right )\\
++ \frac{\sigma_w}{\rho} \frac{\partial \rho}{\partial z} \tau_{L_w} \left (1-r_w \right )
++ \left (1-r_w^2 \right )^{1/2} \zeta \;,
+\end{multline}
+
+where $r_w=\exp(-\Delta t/ \tau_{L_w})$ is the autocorrelation of the vertical
+wind, and $\zeta$ is a normally distributed random number with mean zero and
+unit standard deviation.  The subscripts $k$ and $k+1$ refer to subsequent
+times separated by $\Delta t$.
+
+To save computation time for cases when $\left (\Delta t/\tau_{L_w} \right) <
+0.5$, the following first order approximation is used in order to avoid the
+computation of the exponential function:
+
+\begin{multline}
+\left (\frac{w}{\sigma_w} \right)_{k+1}=
+\left (1-\frac{\Delta t}{\tau_{L_w}} \right )\left (\frac{w}{\sigma_w} \right)_k\\
++ \frac{\partial \sigma_w}{\partial z} \Delta t
++ \frac{\sigma_w}{\rho} \frac{\partial \rho}{\partial z} \Delta t
++ \left( \frac{2 \Delta t}{\tau_{L_w}} \right)^{1/2} \zeta \;.
+\end{multline}
+
+When a particle reaches the surface or the top of the ABL, it is reflected and
+the sign of the turbulent velocity is changed \citep{wilson1993}.
+
+\subsection{\label{timestep}Determination of the time step}
+
+FLEXPART can be used in two different modes.  The computationally faster one
+(\verb|ctl|$<$0 in file \verb|COMMAND|) does not adapt the computation time
+step to the Lagrangian timescales $\tau_{L_i}$ (where $i$ is one of the three
+wind components) and FLEXPART uses constant time steps of one synchronisation
+time interval (\verb|lsynctime|, specified in file \verb|COMMAND|, typically
+900 seconds).  Usually, autocorrelations are very low in this mode and
+turbulence is not described well.  Nevertheless, for large scale applications
+FLEXPART works very well with this option \citep{stohletal1998}.  If turbulence
+shall be described more accurately, the time steps must be limited by $\tau_L$.
+Since the vertical wind is most important, only $\tau_{L_w}$ is used for this.
+The user must specify two constants, \verb|ctl| and \verb|ifine| in file
+\verb|COMMAND|.  The first one determines the time step $\Delta t_i$ according
+to
+
+\begin{equation}
+\Delta t_i=\frac{1}{c_{tl}} \min \left ({\tau_{L_w}\, , \; \frac{h}{2w}\, , \; \frac{0.5}{\partial \sigma_w / \partial z}} \right ) \;.
+\end{equation}
+
+The minimum value of $\Delta t_i$ is 1~second.
+$\Delta t_i$ is used for solving the Langevin equations for the horizontal turbulent wind components.
+
+For solving the Langevin equation for the vertical wind component, a shorter
+time step $\Delta t_w=\Delta t_i / \verb|ifine|$ is used.  However, note that
+there is no interaction between horizontal and vertical wind components on
+timescales less than $\Delta t_i$.  This strategy (given sufficiently large
+values for \verb|ctl| and \verb|ifine|) ensures that the particles stay
+vertically well-mixed also in very inhomogeneous turbulence, while keeping the
+computational cost at a minimum.
+
+\subsection{Parameterization of the wind fluctuations}
+
+For $\sigma_{v_i}$ and $\tau_{L_i}$ \citet{hanna1982} proposed a
+parameterization scheme based on the boundary layer parameters $h$, $L$, $w_*$,
+$z_0$ and $u_*$, i.e.  ABL height, Monin-Obukhov length, convective velocity
+scale, roughness length and friction velocity, respectively.  It is used in
+subroutines \verb|hanna.f, hanna1.f, hanna_short.f| with a modification taken
+from \citet{ryall1997} for $\sigma_w$, as Hanna's scheme does not always yield
+smooth profiles of $\sigma_w$ throughout the whole convective ABL.  In the
+following, subscripts $u$ and $v$ refer to the along-wind and the cross-wind
+components (transformed to grid coordinates in subroutine \verb|windalign.f|),
+respectively, and $w$ to the vertical component of the turbulent velocities;
+$f$ is the Coriolis parameter.  The minimum $\tau_{L_u}$, $\tau_{L_v}$ and
+$\tau_{L_w}$ used are 10~s, 10~s and 30~s, respectively, in order to avoid
+excessive computation times for particles close to the surface.\\[0.3cm]
+
+{\bf Unstable conditions:}
+
+\begin{equation}
+\frac{\sigma_u}{u_*}=\frac{\sigma_v}{u_*}=\left(12+\frac{h}{2|L|}\right)^{1/3}
+\end{equation}
+\begin{equation}
+\tau_{L_u}=\tau_{L_v}=0.15\frac{h}{\sigma_u}
+\end{equation}
+\begin{multline}
+\sigma_w=\\
+\left[ 1.2 w_*^2 \left (1-0.9\frac{z}{h} \right )
+\left ( \frac{z}{h} \right )^{2/3}
++\left (1.8-1.4 \frac{z}{h} \right ) u_*^2 \right]^{1/2}
+\end{multline}
+For $z/h<0.1$ and $z-z_0>-L$:
+\begin{equation}
+\tau_{L_w}=0.1\frac{z}{\sigma_w\left[0.55-0.38\left(z-z_0\right)/L\right]}
+\end{equation}
+For $z/h<0.1$ and $z-z_0<-L$:
+\begin{equation}
+\tau_{L_w}=0.59\frac{z}{\sigma_w}
+\end{equation}
+For $z/h>0.1$:
+\begin{equation}
+\tau_{L_w}=0.15\frac{h}{\sigma_w}\left[1-\exp\left(\frac{-5z}{h}\right)\right]
+\end{equation}
+
+{\bf Neutral conditions:}
+
+\begin{equation}
+\frac{\sigma_u}{u_*}=2.0\exp(-3fz/u_*)
+\end{equation}
+\begin{equation}
+\frac{\sigma_v}{u_*}=\frac{\sigma_w}{u_*}=1.3\exp(-2fz/u_*)
+\end{equation}
+\begin{equation}
+\tau_{L_u}=\tau_{L_v}=\tau_{L_w}=\frac{0.5z/\sigma_w}{1+15fz/u_*}
+\end{equation}
+
+{\bf Stable conditions:}
+
+\begin{equation}
+\frac{\sigma_u}{u_*}=2.0\left(1-\frac{z}{h}\right)
+\end{equation}
+\begin{equation}
+\frac{\sigma_v}{u_*}=\frac{\sigma_w}{u_*}=1.3\left(1-\frac{z}{h}\right)
+\end{equation}
+\begin{equation}
+\tau_{L_u}=0.15\frac{h}{\sigma_u}\left(\frac{z}{h}\right)^{0.5}
+\end{equation}
+\begin{equation}
+\tau_{L_v}=0.07\frac{h}{\sigma_v}\left(\frac{z}{h}\right)^{0.5}
+\end{equation}
+\begin{equation}
+\tau_{L_w}=0.1\frac{h}{\sigma_w}\left(\frac{z}{h}\right)^{0.5}
+\end{equation}
+
+Lacking suitable turbulence parameterizations above the ABL ($z>h$), a constant
+vertical diffusivity $D_z$=0.1~m$^2$s$^{-1}$ is used in the stratosphere,
+following recent work of \citet{legras2003}, whereas a horizontal diffusivity
+$D_h$=50~m$^2$s$^{-1}$ is used in the free troposphere.  Stratosphere and
+troposphere are distinguished based on a threshold of 2~pvu (potential
+vorticity units).  Diffusivities are converted into velocity scales using
+$\sigma_{v_i}=\sqrt{D_i/dt}$.
+
+\subsection{Mesoscale velocity fluctuations}
+
+Mesoscale motions are neither resolved by the ECMWF data nor covered by the
+turbulence parameterization.  This unresolved spectral interval needs to be
+taken into account at least in an approximate way, since mesoscale motions can
+significantly accelerate the growth of a dispersing plume \citep{gupta1997}.
+For this, we use a similar method as \citet{maryon1998}, namely to solve an
+independent Langevin equation for the mesoscale wind velocity fluctuations
+(``meandering'' in Maryon's terms).  Assuming that the variance of the wind at
+the grid scale provides some information on its subgrid variance, the wind
+velocity standard deviation used for the mesoscale Langevin equation is set to
+\verb|turbmesoscale| (set in file \verb|includepar|) times the standard
+deviation of the grid points surrounding the particle's position.  The
+corresponding time scale is taken as half the interval at which wind fields are
+available, assuming that the linear interpolation between the grid points can
+recover half the subgrid variability, not an unlikely assumption
+\citep{stohl1995}.  This empirical approach does not describe actual mesoscale
+phenomena, but it is similar to the ensemble methods used to assess trajectory
+accuracy \citep{kahl1996, baumann1997, stohl1998}.
+
+\subsection{Moist convection}
+
+An important transport mechanism are the updrafts in convective clouds.  They
+occur in conjunction with downdrafts within the clouds and compensating
+subsidence in the cloud-free surroundings.  These convective transports are
+grid-scale in the vertical, but sub-grid scale in the horizontal, and are not
+represented by the ECMWF vertical velocity.
+
+To represent convective transport in a particle dispersion model, it is
+necessary to redistribute particles in the entire vertical column.  For
+FLEXPART we chose the convective parameterization scheme by
+\citet{emanuel1999}, as it relies on the grid-scale temperature and humidity
+fields and calculates a displacement matrix providing the necessary mass flux
+information for the particle redistribution.  The convective parameterization
+is switched on using \verb|lconvection| in file \verb|COMMAND|.  It's
+computation time scales to the square of the number of vertical model levels
+and may account for up to 70\% of FLEXPART's computation time using current
+60-level ECMWF data.
+
+The convection is computed within the subroutines {\tt convmix.f}, {\tt
+calcmatrix.f}, {\tt convect43c.f}, and {\tt redist.f}.  It is called every
+FLEXPART \verb|lsynctime| time step (typically 900~s) with time-interpolated
+temperature and specific humiditiy profiles from the ECMWF data.  Note that the
+original ECMWF model levels, not the Cartesian coordinates, are used in the
+convection scheme.  For efficiency reasons, particles are sorted according to
+their horizontal grid positions ({\tt sort2.f}) before calling the convection
+scheme once per grid column.
+
+In the Emanuel scheme ({\tt convect43c.f}), convection is triggered whenever
+
+\begin{equation}
+T_{vp}^{LCL+1} \ge T_{v}^{LCL+1} + T_{thres}
+\end{equation} 
+
+with $T_{vp}^{LCL+1}$ the virtual temperature of a surface air parcel lifted to
+the level above the lifting condensation level $LCL$, $T_{v}^{LCL+1}$ the
+virtual temperature of the environment there, and $T_{thres}=0.9$~K a threshold
+temperature value.  Based on the buoyancy sorting principle \citep{emanuel1991,
+telford1975}, a matrix $MA$ of the saturated upward and downward mass fluxes
+within clouds is calculated by accounting for entrainment and detrainment:
+
+\begin{multline}
+MA^{i,j}=\\
+\frac{M^i(|\sigma^{i,j+1}- \sigma^{i,j}|+ |\sigma^{i,j} - \sigma^{i,j-1}|)}{(1-\sigma^{i,j}) \displaystyle \sum_{j=LCL}^{LNB} [|\sigma^{i,j+1}- \sigma^{i,j}|+ |\sigma^{i,j} - \sigma^{i,j-1}|]}
+\label{matrix}
+\end{multline}
+
+Here $MA^{i,j}$ are the mass fractions displaced from level $i$ to level $j$,
+$M^{i}$ the mass fraction displaced from the surface to level $i$, $LNB$ the
+level of neutral buoyancy of a surface air parcel and $0 < \sigma^{i,j} < 1 $
+the mixing fraction between level $i$ and level $j$.  The fraction
+$\sigma^{i,j}$ is determined by the environmental potential temperature
+$\theta^{j}$, the liquid potential temperature $\theta_{l}^{i,j}$ of air
+displaced adiabatically from $i$ to $j$, and the liquid potential temperature
+$\theta_{lp}^{i,j}$ of an air parcel first lifted adiabatically to level $i$
+and further to level $j$:
+
+\begin{equation}
+\sigma^{i,j}=\frac{\theta^{j}-\theta_{lp}^{i,j}}{\theta_{l}^{i,j}-\theta_{lp}^{i,j}}
+\end{equation}
+
+By summing up over all levels $j$, we then calculate the saturated up- and
+downdrafts at each level $i$ from Eq.~\ref{matrix} and assume that these fluxes
+are balanced by a subsidence mass flux in the environment.
+
+The particles in each convectively active box are then redistributed ({\tt
+redist.f}) according to the matrix $MA$.  If the mass of an ECMWF model layer
+$i$ is $m^i$ and the mass flux from layer $i$ to layer $j$ accumulated over one
+time step is $\Delta MA^{i,j}$, then the probability of a particle to be moved
+from layer $i$ to layer $j$ is $\Delta MA^{i,j}/m^i$.  Whether a given particle
+is displaced or not is determined by drawing a random number between [0,1],
+which also determines the position of the particle within the destination layer
+$j$.  After the convective redistribution of the particles, the compensating
+subsidence mass fluxes are converted to a vertical velocity acting on those
+particles in the grid box that are not displaced by convective drafts.  By
+calculating a subsidence velocity rather than displacing particles randomly
+between layers the scheme's numerical diffusion in the cloud-free environment
+is eliminated.  The scheme was tested and fulfills the well-mixed criterion,
+i.e., if a tracer is well mixed in the whole atmospheric column, it remains so
+after the convection.
+
+\subsection{Particle splitting}
+
+During the initial phase of dispersion from a point source in the atmosphere,
+particles normally form a compact cloud.  Relatively few particles suffice to
+simulate this initial phase correctly.  After some time, however, the particle
+cloud gets distorted and particles spread over a much larger area.  More
+particles are now needed.  FLEXPART allows the user to specify a time constant
+$\Delta t_s$ (file \verb|COMMAND|).  Particles are split into two (each of
+which receives half of the mass of the original particle) after travel times of
+$\Delta t_s$, $2\Delta t_s$, $4\Delta t_s$, $8\Delta t_s$, and so on
+(subroutine \verb|timemanager.f|).\par
+
+\section{\label{backward}Forward and backward modeling}
+
+Normally, when FLEXPART is run forward in time (\verb|ldirect=1| in file
+\verb|COMMAND|), particles are released from one or a number of sources and
+concentrations are determined downwind on a grid.  However, FLEXPART can also
+run backward in time (\verb|ldirect=-1|), which is more efficient than forward
+modeling for calculating source-receptor relationships if the number of
+receptors is smaller than the number of (potential) sources.  In the backward
+mode, particles are released from a receptor location (e.g., a measurement
+site) and a four-dimensional (3 space dimensions plus time) response function
+(sensitivity) to emission input is calculated.
+
+Since this version (6.2) of FLEXPART, the calculation of the source-receptor
+relationships is generalized for both forward and backward runs, allowing much
+greater flexibility regarding the input and output units than before.
+\verb|ind_source| and \verb|ind_receptor| in file \verb|COMMAND| switch between
+mass and mass mixing ratio units at the source and at the receptor,
+respectively.  Note that source always stands for the physical source and not
+the location of the particle release, which is done at the source in forward
+mode but at the receptor in backward mode.  Table ~\ref{units} gives an
+overview of the units used in forward and backward modeling for different
+settings of the above switches.  A "normal" forward simulation which specifies
+the release in mass units (i.e., kg) and also samples the output in mass units
+(i.e., a concentration in ng m$^{-3}$) requires both switches to be set to 1.
+
+\begin{table}
+\setlength{\tabcolsep}{1.1mm}
+\caption{\label{units}
+Physical units of the input (in file RELEASES) and output data for forward
+(files grid\_conc\_date) and backward (files grid\_time\_date) runs for the
+various settings of the unit switches ind\_source and ind\_receptor (in both
+switches 1 refers to mass units, 2 to mass mixing ratio units).} \vspace{3mm}
+{\centerline{
+\begin{tabular}{ccccc} \hline
+Direction & ind\_source & ind\_receptor & input unit & output unit \\ \hline
+Forward  & 1 & 1 & kg & ng m$^{-3}$ \\
+Forward  & 1 & 2 & kg & ppt by mass \\
+Forward  & 2 & 1 & 1 & ng m$^{-3}$ \\
+Forward  & 2 & 2 & 1 & ppt by mass \\ 
+Backward & 1 & 1 & 1 & s \\
+Backward & 1 & 2 & 1 & s m$^3\,$kg$^{-1}$ \\
+Backward & 2 & 1 & 1 & s kg m$^{-3}$ \\
+Backward & 2 & 2 & 1 & s \\ \hline
+\end{tabular}}}
+\end{table}
+
+In the backward mode, any value not equal zero can be entered as the release
+"mass" in file \verb|RELEASES| because the output is normalized by this value.
+The calculated response function is related to the particles' residence time in
+the output grid cells.  The unit of the output response function varies,
+depending on how the switches are set.  If \verb|ind_source=1| and
+\verb|ind_receptor=1|, the response function has the unit s.  If this response
+function is folded (i.e., multiplied) with a 3-d field of emission mass fluxes
+into the output grid boxes (in kg~m$^{-3}$s$^{-1}$), a concentration at the
+receptor (kg~m$^{-3}$) is obtained.  If \verb|ind_source=1| and
+\verb|ind_receptor=2|, the response function has the unit s~m$^3\,$kg$^{-1}$
+and if it is folded with the emission mass flux (again in kg~m$^{-3}$s$^{-1}$),
+a mass mixing ratio at the receptor is obtained.  The units of the response
+function for \verb|ind_source=2| can be understood in analogy.
+
+In the case of loss processes (dry or wet deposition, decay) the response
+function is ``corrected'' for these loss processes.  See \citet{seibert2001}
+and, particularly, \citet{seibert2004} for a description of these generalized
+in- and output options and the implementation of backward modeling in FLEXPART.
+\citet{seibert2004} also describe the theory of backward modeling and give some
+examples, and \citet{stohl2003} presents an application.
+
+\section{\label{plumetraj}Plume trajectories}
+
+In a recent paper, \citet{stohl2002} proposed a method to condense the complex
+and large FLEXPART output using a cluster analysis \citep{dorling1992}.  The
+idea behind this is to cluster, at every output time, the positions of all
+particles originating from a release point, and write out only clustered
+particle positions, along with additional information (e.g., fraction of
+particles in the ABL and in the stratosphere).  This creates information that
+is almost as compact as traditional trajectories but accounts for turbulence
+and convection.  This option can be activated by setting \verb|iout| to 4 or 5
+in file \verb|COMMAND|.  The number of clusters can be set with the parameter
+\verb|ncluster| in file \verb|includepar|.  The clustering is handled and
+output is produced by subroutine \verb|plumetraj.f|.
+
+\section{\label{removal}Removal processes}
+
+FLEXPART takes into account radioactive (or other) decay, wet deposition, and
+dry deposition by reducing a particle's mass.  However, as atmospheric
+transport is the same for all chemical species, a single particle can represent
+several (up to \verb|maxspec|) chemical species, each affected differently by
+the removal processes.
+
+\subsection{\label{radioactive}Radioactive decay}
+
+Radioactive decay is accounted for by reducing the particle mass according to
+
+\begin{equation}
+m(t+\Delta t)=m(t) \exp (-\Delta t /\beta) \;,
+\end{equation}
+
+where $m$ is particle mass, and the time constant $\beta=T_{1/2}/\ln(2)$ is
+determined from the half life $T_{1/2}$ specified in file \verb|SPECIES_nnn|.
+Deposited pollutant mass decays at the same rate.
+
+\subsection{OH Reaction}
+
+If a positive value for the OH reaction rate is given in the file
+\verb|SPECIES_nnn|, OH reaction is performed in the model simulation and tracer
+mass is lost by this reaction.  A monthly averaged 3\degreee$\times$ 5\degreen resolution OH
+field averaged to 7 atmospheric levels is used.  The fields have been obtained
+from the GEOS-CHEM model \citep{bey2001}.  The reaction rate is temperature
+corrected and an activation rate of 1000 J\,mol$^{-1}$ is assumed.
+
+\subsection{\label{wetdepo}Wet deposition}
+
+Since version 8.0, in-cloud and below-cloud scavenging are treated differently.
+Based on the humidity and temperature from the meteorological input data, the
+occurence of clouds is calculated.  If relative humidity exceeds 80\% the
+occurence of a cloud is assumed.
+
+{\bf In cloud scavenging} is treated differently for gases and particles.  
+The implementation follows the scheme of \citet{hertel1995}.
+
+In general, the scavenging coefficient $\Lambda$ ((s$^{-1}$) depends on the precipitation
+rate $I$ (mm\,h$^{-1}$) and the height $H_i$ over which scavenging takes place:
+
+\begin{equation}
+\Lambda=\frac{S_i I}{H_i}
+\end{equation}
+
+$S_i$ is different for gases and particles. For particles,
+
+\begin{equation}
+S_i=0.9/cl
+\end{equation}
+
+where $cl$ is the cloud liquid water content
+
+\begin{equation}
+cl=2\times 10^{-7}\cdot I^{0.36} \;.
+\end{equation}
+
+For gases,
+
+\begin{equation}
+S_i=1/cl_\text{eff} \;,
+\end{equation}
+
+where $cl_\text{eff}$ is an effective cloud liquid water content:
+
+\begin{equation}
+cl_\text{eff}=\frac{(1-cl)}{H_\text{eff}RT}+cl \; .
+\end{equation}
+
+{\bf Below cloud scavenging} takes the form of an exponential decay process
+\citep{mcmahon1979}.  With the scavenging coefficient $\Lambda$, wet
+deposition is described as
+
+\begin{equation}
+m(t+\Delta t)=m(t) \exp (- \Lambda \Delta t) \;,
+\end{equation}
+
+where $m$ is the particle mass.  $\Lambda$ increases with the precipitation
+rate $I$ according to
+
+\begin{equation}
+\Lambda=AI^{B} \;,
+\label{wetscav}
+\end{equation}
+
+where $A$ [s$^{-1}$] is the scavenging coefficient at $I=$1~mm/hour and $B$
+gives the dependency on precipitation rate.  Both $A$ and $B$ must be specified
+in file \verb|SPECIES_nnn|.  FLEXPART uses the same scavenging coefficients for
+snow and rain.
+
+As wet deposition depends nonlinearly on precipitation rate, subgrid
+variability of precipitation must be accounted for \citep{hertel1995}.  The
+area fraction which experiences precipitation given a certain grid-scale
+precipitation rate is calculated by
+
+\begin{equation}
+F=\max\left[0.05,CC\frac{I_l fr_l(I_l) + I_c fr_c(I_c)}{I_l+I_c}\right] \;,
+\end{equation}
+
+where $CC$ is the total cloud cover, $I_l$ and $I_c$ are the large scale and
+convective precipitation rates, respectively, and $fr_l$ and $fr_c$ are
+correction factors that depend on $I_l$ and $I_c$ (see Table~\ref{corrfacts}).
+The subgrid scale precipitation rate is then $I_s=(I_l+I_c)/F \;$.
+
+\begin{table}
+\setlength{\tabcolsep}{1.1mm}
+\caption{\label{corrfacts}
+Correction factors used for the calculation of the area fraction that
+experiences precipitation.  Precipitation rates are in mm/hour.} \vspace{3mm}
+{\centerline{
+\begin{tabular}{cccccc} \hline
+~&\multicolumn{5}{c|}{$I_l$ and $I_c$} \\ \hline
+Factor&$I \le 1$&$1 < I \le 3$&$3 < I\le8$&$8 < I \le 20$&$20 < I$ \\ \hline
+$fr_l$&0.50&0.65&0.80&0.90&0.95 \\
+$fr_c$&0.40&0.55&0.70&0.80&0.90 \\ \hline
+\end{tabular}}}
+\end{table}
+
+\subsection{Dry deposition}
+
+Dry deposition is described in FLEXPART by a deposition velocity 
+
+\begin{equation}
+v_d(z)=-F_C/C(z) \;,
+\end{equation}
+
+where $F_C$ and $C$ are the flux and the concentration of a species at height
+$z$ within the constant flux layer.  A constant deposition velocity $v_d$ can
+be set (file \verb|SPECIES_nnn|).  Alternatively, if the physical and chemical
+properties of a substance are known (file \verb|SPECIES_nnn|), more complex
+parameterizations for gases and particles are also available.
+
+\subsubsection{Dry deposition of gases}
+
+The deposition velocity of a gas is calculated with the {\em resistance method}
+\citep{wesely1977} in subroutine \verb|getvdep.f| according to
+
+\begin{equation}
+|v_d(z)|=\left[r_a(z)+r_b+r_c\right]^{-1} \;,
+\end{equation}
+
+where $r_a$ is the aerodynamic resistance between $z$ and the surface, $r_b$ is
+the quasilaminar sublayer resistance, and $r_c$ is the bulk surface resistance.
+
+The aerodynamic resistance $r_a$ is calculated in function \verb|raerod.f|
+using the flux-profile relationship based on Monin-Obukhov similarity theory
+\citep{stull1988}
+
+\begin{equation}
+r_a(z)=\frac{1}{\kappa u_*}[\ln(z/z_0)- \Psi _h(z/L)+ \Psi _h(z_0/L)] \;.
+\label{ra}
+\end{equation}
+
+Following \citet{erisman1994}, the quasilaminar sublayer resistance is
+
+\begin{equation}
+r_b=\frac{2}{\kappa u_*} \left(\frac{Sc}{Pr}\right)^{2/3} \;,
+\end{equation}
+
+where $Sc$ and $Pr$ are the Schmidt and Prandtl numbers, respectively.  $Pr$ is
+0.72 and $Sc=\upsilon/D_i$, with $\upsilon$ being the kinematic viscosity of
+air and $D_i$ being the molecular diffusivity of species $i$ in air.  The
+slight dependency of $\upsilon$ on air temperature is formulated in accordance
+with \citet{pruppacher1978}.  $r_b$ is calculated in function
+\verb|getrb.f|.\par
+
+The surface resistance is calculated in function \verb|getrc.f| following
+\citet{wesely1989} as
+
+\begin{equation}
+\frac{1}{r_c}=
+\frac{1}{r_s+r_m}+\frac{1}{r_{lu}}+\frac{1}{r_{dc}+r_{cl}}+\frac{1}{r_{ac}+r_{gs}} \;,
+\end{equation}
+
+where $r_s$, $r_m$ and $r_{lu}$ represent the bulk values for leaf stomatal,
+leaf mesophyll and leaf cuticle surface resistances (alltogether the upper
+canopy resistance) , $r_{dc}$ represents the gas-phase transfer affected by
+buoyant convection in canopies, $r_{cl}$ the resistance of leaves, twig, bark
+and other exposed surfaces in the lower canopy, $r_{ac}$ the resistance for
+transfer that depends only on canopy height and density, and $r_{gs}$ the
+resistance for the soil, leaf litter, etc., at the ground.  Each of these
+resistances is parameterized according to the species' chemical reactivity and
+solubility, the landuse type, and the meteorological conditions.  The IGBP
+landuse inventory \citep{belward1999} provides the area fractions of 13 landuse
+classes for which roughness lengths $z_0$ are estimated, on a grid with 0.3
+resolution (Table~\ref{landuses}).  Charnock's relationship \citep{stull1988}
+$z_0=0.016u_*^2/g$ is used to calculate $z_0$ for the classes ``Ocean'' and
+``Inland water'', because of its dependence on wave height.  Deposition
+velocities are calculated for all landuse classes and weighted with their
+respective areas.  The original resolution of the IGBP land cover
+classification is 1 km x 1 km.  To save storage space, this was regridded to a
+0.3 x 0.3 degree resolution.  Table~\ref{conversion} shows how the initial 17
+classse were transfered in the 13 classes used in the deposition scheme of
+\cite{wesely1989}.  According to \cite{helmig2007} a resistance of snow to
+ozone deposition of 10000 s m$^{-1}$ was used.  For SO$_2$ a value of 100
+according to \citet{zhang2002} was used in the model.  As the Wesely scheme was
+initially developed for North America, the rain forest was not represented
+well.  Therefore a new category was introduced and resistance values according
+to \cite{jacob1990} were used.  The resistance values are dependent on 5
+different seasons.  As on the southern hemisphere they are asynchronous to the
+northern hemisphere half a year was added when choosing the appropriate
+seasonal category.  The snow depth is included in the original
+\cite{wesely1989} parameterization to some extend, as at a certain time of the
+year snow cover is assumed.  In order to have more accurate information we use
+the snow cover based on the snow depth in the ECMWF fields.  If the snow cover
+gets over 1 mm of water equivalent, the grid cell is considered as snow covered
+and category 12 (snow and ice) is applied.
+
+\begin{table}
+\setlength{\tabcolsep}{1.1mm}
+\caption{\label{landuses}
+List of the landuse classes and roughness lengths used by FLEXPART.
+``Charnock'' indicates that Charnock's relationship is used to calculate the
+roughness length.} \vspace{3mm}
+{\centerline{
+\begin{tabular}{lc} \hline
+Urban land & 0.7	\\
+Agricultural land	&  0.1 \\
+Range land & 0.1 \\
+Deciduous forest &  1.0 \\
+Coniferous forest & 1.0 \\
+Mixed forest including wetland &	0.7	\\
+Water, both salt and fresh		&   Charnock \\
+Barren land mostly desert &  0.01 \\
+Nonforested wetland &	0.1 \\
+Mixed agricultural and range land &	0.1 \\
+Rocky open areas with low growing shrubs  & 	0.05 \\
+Snow and ice & 0.001 \\
+Rainforest &   1.0 \\ \hline
+\end{tabular}}}
+\end{table}
+
+\begin{table*}
+\setlength{\tabcolsep}{1.1mm}
+\caption{\label{conversion}
+Conversion from the IGBP land cover legend to the landuse categories used by Wesely}
+\vspace{3mm}
+{\centerline{
+\begin{tabular}{clcl} \hline
+\multicolumn{2}{c}{IGBP Land Cover Legend} & \multicolumn{2}{c}{Wesely} \\
+Value &	Description	&		Value &	Description \\
+1& 	Evergreen Needleleaf Forest&			5&	Coniferous forest\\
+2& 	Evergreen Broadleaf Forest	&	     13& Rainforest\\
+3& 	Deciduous Needleleaf Forest&			4&	Deciduous forest\\
+4& 	Deciduous Broadleaf Forest	&			4&	Deciduous forest\\
+5& 	Mixed Forest&6&	Mixed forest including wetland\\
+6& 	Closed Shrublands&	      11& 	rocky open areas with low growing shrubs\\
+7& 	Open Shrublands&		   11& 	rocky open areas with low growing shrubs\\
+8& 	Woody Savannas&11& 	rocky open areas with low growing shrubs\\
+9& 	Savannas&11&rocky open areas with low growing shrubs\\
+10& 	Grasslands			         	&		3	& Range land\\
+11& 	Permanent Wetlands			   &      9	& nonforested wetland	\\
+12& 	Croplands						   &   	2	& Agricultural land\\
+13& 	Urban and Built-Up				&		1	& Urban land\\
+14& 	Cropland/Natural Vegetation Mosaic&	10	& mixed agricultural and range land		\\
+15& 	Snow and Ice							 &  12 & snow and ice\\
+16& 	Barren or Sparsely Vegetated		&	8& 	barren land mostly desert\\
+17& 	Water Bodies							&   7	& water, both salt and fresh\\
+\end{tabular}}}
+\end{table*}
+
+\subsubsection{Dry deposition of particulate matter}
+
+The deposition of particulates is calculated in subroutine \verb|partdep.f|
+according to
+
+\begin{equation}
+v_d(z)=\left[r_a(z)+r_b+r_a(z)r_bv_g\right]^{-1}+v_g \;,
+\end{equation}
+
+where $v_g$ is the gravitational settling velocity calculated from
+\citep{slinn1982}
+
+\begin{equation}
+v_g=\frac{g \rho_p d_p^2C_{cun}}{18 \mu} \;,
+\end{equation}
+
+where $\rho_p$ and $d_p$ are the particle density and diameter, $\mu$ the
+dynamic viscosity of air (0.000018 kg m$^{-1}$s$^{-1}$) and $C_{cun}$ the
+Cunningham slip-flow correction.  The quasilaminar sublayer resistance is
+calculated from the same relationship as for gases, with an additional
+impaction term.  For further details see \citet{slinn1982}.\par
+
+Settling and dry deposition velocities are strongly dependent on particulate
+size.  FLEXPART assumes a logarithmic normal size distribution of the
+particulate mass.  The user must specify the mean particulate diameter
+$\overline{d_p}$ and a measure of the variation around $\overline{d_p}$,
+$\sigma_p$.  Then, the settling and deposition velocities are calculated for
+several particle diameters and are weighted with their respective particulate
+mass fractions.
+
+Gravitational settling is important not only for the computation of the dry
+deposition velocity, but also affects the particle's trajectory.  As a FLEXPART
+particle can normally represent several species, gravitational settling can
+only be taken into account correctly (i.e., influence particle trajectories) in
+single-species simulations.  Pay attention that gravitational settling is
+simply switched off in multi-species simulations, without warning.
+
+With FLEXPART version 8.2, the temperature dependence of the dynamic viscosity
+is taken into account.  Furthermore, we have extended the calculation of
+settling velocities to higher Reynolds numbers.  For this, an iterative
+procedure has been introduced in subroutine \verb|get_settling.f|, where the
+iteration is started with Stokes’ law settling \citep{naeslund1991}, and then
+Reynolds number and settling velocity are calculated until convergence is
+achieved.
+
+\subsubsection{Loss of particle mass due to dry deposition}
+
+The depositon velocity is calculated for a reference height (parameter
+\verb|href| in file \verb|includepar|) of 15~m.  For all particles below
+$2h_{ref}$, the mass lost by deposition is calculated by
+
+\begin{equation}
+\Delta m(t)=m(t)\left[{1-\exp\left({\frac{-v_d(h_{ref})\Delta t}{2h_{ref}}}\right)}\right] \;.
+\end{equation}
+
+\section{\label{conccalc}Calculation of concentrations, uncertainties, age spectra, and mass fluxes}
+
+Output quantities $C_{T_c}$ at time $T_c$ (output interval \verb|loutstep| is
+set in file \verb|COMMAND|) are calculated as time-averages over period
+$[T_c-\Delta T_c/2,T_C+\Delta T_c/2]$.  $\Delta T_c$ must be specified
+(\verb|loutaver|) in file \verb|COMMAND|.  To calculate the time-averages,
+concentrations $C_{T_s}$ at times $T_s$ within $[T_c-\Delta T_c/2,T_C+\Delta
+T_c/2]$ are sampled at shorter intervals $\Delta T_s$ (\verb|loutsample| in
+file \verb|COMMAND|) and are then divided by the number $N=\frac{\Delta
+T_c}{\Delta T_s}$ of samples taken:
+
+\begin{equation}
+C_{T_c}= \frac{1}{N} \sum_{i=1}^N {C_{T_s}} \;.
+\end{equation}
+
+Both $\Delta T_c$ and $\Delta T_s$ must be multiples of the FLEXPART
+synchronisation interval (\verb|lsynctime| in file \verb|COMMAND|).  The
+shorter the sampling interval $\Delta T_s$, the more samples are taken and the
+more accurate are thus the time-averaged concentrations.
+
+\subsection{Concentrations, mixing ratios, and emission response functions}
+
+The user can choose (\verb|iout| in file \verb|COMMAND|, which must be set to 1
+for backward runs) whether concentrations, volume mixing ratios or both shall
+be produced.  We shall use the term "concentration" and particle mass here, but
+note that the actual units are determined by the settings of \verb|ind_source|
+and \verb|ind_receptor|, according to Table~\ref{units}.  The concentration in
+a grid cell is calculated in subroutine \verb|conccalc.f| by sampling the
+tracer mass fractions of all particles within the grid cell and dividing by the
+grid cell volume
+
+\begin{equation}
+C_{T_s} =\frac{1}{V}\sum_{i=1}^N (m_i f_i) \;,
+\end{equation}
+
+with $V$ being the grid cell volume, $m_i$ particle mass, $N$ the total number
+of particles, and $f_i$ the fraction of the mass of particle $i$ attributed to
+the respective grid cell.  This mass fraction is calculated by a uniform kernel
+with bandwidths $(\Delta x,\Delta y)$, where $\Delta x$ and $\Delta y$ are the
+grid distances on the longitude-latitude output grid.  Figure~\ref{kernel}
+illustrates this: The particle is located at the center of the shaded rectangle
+with side lengths $(\Delta x, \Delta y)$.  Generally, the shaded area stretches
+over four grid cells, each of which receives a fraction of the particle's mass
+equal to the fraction of the shaded area falling within this cell.  The uniform
+kernel is not used during the first 3 hours after a particle's release (when
+the mass is attributed only to the grid cell it resides in), in order to avoid
+smoothing close to the source.
+
+\begin{figure}[htb]
+\begin{minipage}[t]{2.8cm}
+\end{minipage}\hfill
+{\begin{minipage}[t]{12.5cm}
+
+\setlength{\unitlength}{2.5cm}
+\begin{picture}(3.0,3.0)
+
+\thicklines
+\multiput(0.5,0)(1,0){3}{\line(0,1){3.0}}
+\multiput(0,0.5)(0,1){3}{\line(1,0){3.0}}
+
+\thinlines
+\put(0.5,0.35){\vector(1,0){1.0}}
+\put(1.5,0.35){\vector(-1,0){1.0}}
+\put(0.35,0.5){\vector(0,1){1.0}}
+\put(0.35,1.5){\vector(0,-1){1.0}}
+
+\put(0.9,0.15){$\Delta x$}
+\put(0.05,0.92){$\Delta y$}
+
+\multiput(1.2909,1.25)(0.0909,0){10}{\line(0,1){1.0}}
+\multiput(1.2,1.34909)(0,0.0909){10}{\line(1,0){1.0}}
+
+\thicklines
+\put(1.67,1.76){\line(1,0){0.06}}
+\put(1.70,1.73){\line(0,1){0.06}}
+
+\put(1.2,1.25){\line(1,0){1.0}}
+\put(2.2,1.25){\line(0,1){1.0}}
+\put(2.2,2.25){\line(-1,0){1.0}}
+\put(1.2,2.25){\line(0,-1){1.0}}
+
+
+\end{picture}
+\end{minipage}}
+\caption{\label{kernel} Illustration of the uniform kernel used to calculate
+gridded concentration and deposition fields.  The particle position is marked
+by ``{\rm +}''}
+\end{figure}
+
+Wet and dry deposition fields are calculated on the same output grid
+(subroutines \verb|wetdepokernel.f| and \verb|drydepokernel.f|) and are written
+to all output grid files.  The deposited matter is accumulated over the course
+of a model run, i.e.  it generally increases with model time.  However,
+radioactive decay is calculated also for the deposited matter.
+
+\subsection{Uncertainties}
+
+The uncertainty of the output is estimated by carrying \verb|nclassunc| classes
+of particles in the model simulation, and determining the concentration
+separately for each class (subroutine \verb|conccalc.f|).  The standard
+deviation, calculated from \verb|nclassunc| concentration estimates and divided
+by $\sqrt{\tt{nclassunc}}$, is the standard deviation of the mean concentration
+(subroutine \verb|concoutput.f|), which is also written to the output files for
+every grid cell.  Note that the memory needed for some auxiliary fields
+increases with \verb|nclassunc| {ít and} the number of age classes (see below).
+It may, thus, be necessary to reduce \verb|nclassunc| for runs with large
+output grids and age spectra calculations or in the backward mode.
+
+\subsection{Age spectra}
+
+The age spectra option is switched on using \verb|lagespectra| in file
+\verb|COMMAND|, with the age classes specified in seconds in file
+\verb|AGECLASSES|.  Concentrations are split into contributions from particles
+of different age, defined as the time passed since their release.  Particles
+are terminated once they are older than the oldest age class and their storage
+space is made available to new particles.  Therefore, the age spectra option
+can be used also with a single age class for defining a maximum particle age.
+
+\subsection{Parabolic kernel}
+
+In addition to the simple uniform kernel method, a computationally demanding
+parabolic kernel as described in \citep{uliasz1994} can be used to calculate
+surface concentrations for a limited number of receptor points (age spectra are
+not available in this case):
+
+\begin{equation}
+C_{T_s}(x,y,z=0)=\sum_{i=1}^N \left[ \frac{2m_iK(r_x,r_y,r_z)}{h_{x_i}h_{y_i}h_{z_i}} \right]\;,
+\end{equation}
+
+where $h_{x_i}$, $h_{y_i}$ and $h_{z_i}$ are the kernel bandwidths which
+determine the degree of smoothing, $r_x=(X_i-x)/h_{x_i}$,
+$r_y=(Y_i-y)/h_{y_i}$, $r_z=Z_i/h_{z_i}$ with $X_i$, $Y_i$ and $Z_i$ being the
+position of particle $i$.  The kernel bandwidths are a function of the
+particles' age.
+
+\subsection{Mass fluxes}
+
+Mass flux calculations can be switched on using \verb|iflux| in file
+\verb|COMMAND|.  Mass fluxes are calculated separately for eastward, westward,
+northward, southward, upward and downward directions and contain both
+grid-scale and subgrid-scale motions.  Mass fluxes are determined for the
+centerlines of the output grid cells, e.g.  vertical fluxes are calculated for
+motions across the half level of each output cell.
+
+\section{Domain-filling option}
+
+\subsection{General}
+
+If \verb|mdomainfill=1| in file \verb|COMMAND| particles are not released at
+specific locations.  Instead, the longitudes and latitudes specified for the
+first release in the \verb|RELEASES| file are used to set up a global or
+limited model domain.  The particles (number is also taken from
+\verb|RELEASES|) are then distributed in the model domain proportionally to air
+density (subroutine \verb|init_domainfill.f|).  Each particle receives the same
+mass, altogether accounting for the total atmospheric mass.  Subsequently,
+particles move freely in the atmosphere.
+
+If a limited domain is chosen, mass fluxes are determined in small grid boxes
+at the boundary of this domain (boundaries must be at least one grid box away
+from the boundaries of the meteorological input data).  In the grid cells with
+air flowing into the model domain, mass fluxes are accumulated over time and
+whenever the accumulated mass exceeds the mass of a particle, a new particle
+(or more, if required) is released at a randomly chosen position at the
+boundary of the box (subroutine \verb|boundcond_domainfill.f|).  At the
+outflowing boundaries particles are terminated.  Note that, due to the change
+of mass of the atmosphere in the model domain and due to numerical effects, the
+number of particles used is not exactly constant throughout the simulation.
+
+\subsection{Stratospheric ozone tracer}
+
+If \verb|mdomainfill=2|, the domain-filling option is used to simulate a
+stratospheric ozone tracer.  Upon particle creation, the potential vorticity
+(PV) at its position is determined by interpolation from the ECMWF data.
+Particles initially located in the troposphere (PV$<$\verb|pvcrit| potential
+vorticity units (pvu), default 2 pvu) are not used.  In contrast, stratospheric
+particles (PV$>$\verb|pvcrit|) are given a mass according to:
+
+\begin{equation}
+M_{O_3}=M_{air}\, P \, C \, 48/29
+\end{equation}
+
+where $M_{air}$ is the mass of air a particle represents, $P$ is PV in pvu,
+$C=60\times$10$^{-9}$ pvu$^{-1}$ is the ozone/PV relationship \citep{stohl2000}
+(parameter \verb|ozonescale|), and the factor 48/29 converts from volume to
+mass mixing ratio.  Particles are then allowed to advect through the
+stratosphere and into the troposphere according to the winds.
+
+\section{Model output}
+
+Tracer concentrations and/or mixing ratios (for forward runs), or emission
+sensitivity response functions (for backward runs) are calculated on a
+three-dimensional longitude-latitude grid, defined in file \verb|OUTGRID|,
+whose domain and resolution can differ from the grid on which meteorological
+input data are given.  Two-dimensional wet and dry deposition fields are
+calculated over the same spatial domain, and tracer mass fluxes can also be
+determined on the 3-d grid.  Except for the mass fluxes, output can also be
+produced on one nested output grid with higher horizontal but the same vertical
+resolution, defined in file \verb|OUTGRID_NEST|.  For certain locations,
+specified in file \verb|RECEPTORS|, concentrations can also be calculated
+independently from a grid (see below).  The time interval (variable
+\verb|loutstep|) at which output is produced is read in from file
+\verb|COMMAND|.  For every output time and for every species (\verb|nnn|),
+files are created, with file names ending with date, time and species number in
+the format \verb|yyyymmddhhmmss_nnn|.  A list of all these output times is
+written to the formatted file \verb|dates|.  The dates indicate the ending time
+of an output sampling interval (see section~\ref{conccalc}).
+
+\subsection{Gridded output}
+
+There are several output options in FLEXPART, which can all be selected in file
+\verb|COMMAND|.  Gridded output fields can be concentrations (files
+\verb|grid_conc_date_nnn|), volume mixing ratios (files
+\verb|grid_pptv_date_nnn|), emission response sensitivity in backward
+simulations (files \verb|grid_time_date_nnn|), or fluxes (files
+\verb|grid_flux_date|, unit 10$^{-12}$ kg m$^{-2}$ s$^{-1}$ for forward runs),
+or, in backward mode, sensitivities to initial conditions taken from another
+model (\verb|grid\_initial\_nnn|).
+
+The species number identifier \verb|nnn| starts at one (with leading zeros) and
+increases to the maximum number of species used in the simulation.  Files
+\verb|grid_conc_date_nnn| are created only in forward runs, whereas files
+\verb|grid_time_date_nnn| are only created in backward runs.  Note that the
+units of the files \verb|grid_conc_date_nnn| and \verb|grid_time_date_nnn|
+depend on the settings of the switches \verb|ind_source| and
+\verb|ind_receptor|, following Table~\ref{units}.  In particular, the units of
+\verb|grid_conc_date| can also be mass mixing ratios.  For forward runs,
+additional files \verb|grid_pptv_date_nnn| can be created, which contain volume
+mixing ratios for gases.  Output files \verb|grid_conc*|, \verb|grid_pptv*|,
+and \verb|grid_time*| also contain wet and dry deposition fields (unit
+10$^{-12}$ kg m$^{-2}$ in forward mode), and all files contain, for each grid
+cell, corresponding uncertainties.  All these file types share a common header,
+file \verb|header| produced by subroutine \verb|writeheader.f|, where important
+information on the model run (start of simulation, grid domain, number and
+position of vertical levels, age classes, release points, etc.) is stored.  In
+all postprocessing programs, the header must be read in before the actual data
+files.  File names for the output nests follow the same nomenclature as
+described above, but with \verb|_nest| added (e.g., \verb|header_nest|, or
+\verb|grid_conc_nest_date_nnn|).  The output files are written with subroutines
+\verb|concoutput.f| and \verb|fluxoutput.f|.
+
+FLEXPART output typically contains many grid cells with zero values.  It would
+be inefficient to write out all these zeroes.  Therefore, a special format has
+been designed that compresses the information to the relevant information.  In
+previous versions (up to version 7), the output consisted of a mixture of a
+full grid dump (including zeroes) and a sparse matrix format - output was
+switched between these two formats based on which one was smaller.
+
+In version 8.0, the output has been redesigned completely for yet more efficiency
+such that output file sizes are only about 40\% of what they used to be.  The
+output grid is searched for consecutive sequences of non-zero values.  The
+variable \verb|sp_count_i| gives the number of such sequences (n), and the
+integer field \verb|sparse_dump_i(n)| contains the field positions of the first
+non-zero element for every sequence.  The variable \verb|sp_count_r| gives the
+total number (k) of non-zero values written out.  They are contained in the
+real field \verb|sparse_dump_r(n)|.  Since all physical output quantities of
+FLEXPART are non-zero, sequences are written out alternatingly as positive or
+negative values.  Every switch between positive and negative values indicates
+that a new sequence with non-zero values starts.  The field position for that
+start is contained in \verb|sparse_dump_i(k)| for the k-th switch between
+positive and negative.  Zero values in between sequences are ignored and not
+written out.  Field positions within the 3-d output field are coded such that a
+single integer value is sufficient.  It can later be converted back to give
+positions in all three coordinates.
+
+The possibility to output the sensitivities to initial conditions for backward
+runs has been introduced with FLEXPART version 8.2. This option can be used to
+calculate the exact sensitivities of the concentrations (or mixing ratios,
+depending on what is simulated) to initial conditions either in concentration
+or mixing ratio units taken from a gridded data set, which can be produced
+either by another model or by a FLEXPART forward simulation.  This option has
+been introduced to allow interfacing FLEXPART with other models.  Multiplying
+the sensitivities in files \verb|grid_initial_nnn| with the corresponding
+concentrations (or mixing ratios, depending on option chosen for switch
+\verb|LINIT_COND| in file \verb|COMMAND|) from the forward simulation at the
+interfacing time gives the concentration (or mixing ratio) response at the
+receptor (taking into account possible loss processes during the FLEXPART
+simulation time).
+
+\subsection{Receptor point output}
+
+For a list of points at the surface, concentrations or mixing ratios in forward
+simulations can be determined with a grid-independent method.  This information
+is written to files \verb|receptor_conc| and \verb|receptor_pptv|,
+respectively, for all dates of a simulation.
+
+\subsection{Particle dump and warm start option}
+
+Particle information (3-d position, release time, release point, and release
+masses for all species) can be written out to files (subroutine
+\verb|partoutput.f|) either continuously (binary files \verb|partposit_date|),
+or `only at the end' of a simulation (file \verb|partposit_end|).  In both
+cases output is written every output interval but file \verb|partposit_end| is
+overwritten upon each new output.  If FLEXPART must be terminated, it can be
+continued later on by reading in files \verb|header| and \verb|partposit_end|
+produced by the previous run (subroutine \verb|readpartpositions.f|).  Such a
+warm start is done if variable \verb|ipin| is set to 1 in file \verb|COMMAND|.
+
+If option \verb|mquasilag| is chosen in file \verb|COMMAND|, particle dumps
+every output interval are produced in a very compact format by converting the
+positions to an \verb|integer*2| format (subroutine \verb|partoutput_short.f|).
+As some accuracy is lost in the conversion, this output is not used for the
+warm start option.  Another difference to the normal particle dump is that
+every particle gets a unique number, thus allowing postprocessing routines to
+identify continuous particle trajectories.
+
+\subsection{Clustered plume trajectories}
+
+Condensed particle output using the clustering algorithm described in
+section~\ref{plumetraj} is written to the formatted file
+\verb|trajectories.txt|.  Information on the release points (coordinates,
+release start and end, number of particles) is written by subroutine
+\verb|openouttraj.f| to the beginning of file \verb|trajectories.txt|.
+Subsequently, \verb|plumetraj.f| writes out a time sequence of the clustering
+results for each release point: release point number, time in seconds elapsed
+since the middle of the release interval, plume centroid position coordinates,
+various overall statistics (e.g., fraction of particles residing in the ABL and
+troposphere), and then for each cluster the cluster centroid position, the
+fraction of particles belonging to the cluster, and the root-mean-square
+distance of cluster member particles from the cluster centroid.
+
+\section{pflexpart: a Python routine for the analysis of FLEXPART output}
+
+To assist in the usage and analysis of FLEXPART data we have created a Python
+module that is available with this new release.  The Python module 'pflexpart'
+enables the user to easily read and access the header and grid output data of
+the FLEXPART model runs.  Furthermore, we provide some basic classes that
+assist in conducted standard analysis of backward and forward model runs.
+
+The module is released under the same GPL license as FLEXPART.  As open source
+code it is constantly undergoing revision and updates from the community.
+Thus, the core functionality of the module is described online.  See:
+http://transport.nilu.no/pflexpart
+
+The module is largely based on the well known Numpy and matplotlib Python
+modules as well as the {``basemap``} tool box of matplotlib.  Additionally, for
+efficiency in reading the data, there are several custom build modules that use
+the \verb|readgrid.f| FORTRAN routines described above.  Note, that some of
+these modules will require the user to compile them for their system to achieve
+maximum benefit.
+
+The central framework to the module is the pflexpart \verb|Header| class.  This
+class will read a FLEXPART header file and provide some functionality toward
+data analysis.  For instance, the \verb|Header| class has a
+\verb|fill_backward| method for backward runs that will calculate the Total
+Column and Footprint sensitivities from the N ageclasses used in the run.  From
+this method, plotting of residence times and emission sensitivities is
+relatively simple.
+
+Other features include wrapper functions around the matplotlib toolboxes that
+allow for plotting of the data easily.  The functions \verb|plot_totalcolumn|
+and \verb|plot_footprint| for example are customized wrappers of the
+\verb|plot_sensitivity| function.  These functions take data objects that are
+created by the \verb|fill_backward| method of the \verb|Header| class and allow
+the user to create plots quickly.  For forward runs, the \verb|Header.readgrid|
+method can be used, again with the \verb|plot_sensitivity| function of the
+pflexpart module.
+
+For more details, the reader is referred again to the more frequently updated
+module home page.
+
+\section{Final remark}
+
+In this note, we have described the Lagrangian particle dispersion model FLEXPART in version
+8.2.  As FLEXPART is developed further this text will be kept up to date and will be accessible
+from the FLEXPART home page at http://transport.nilu.no/flexpart.
+
+\begin{thebibliography}{}
+
+\bibitem[Asman(1995)]{asman1995}
+Asman, W. A. H.:
+Parameterization of below-cloud scavenging of highly soluble gases under convective conditions.
+Atmos. Environ., 29, 1359--1368, 1995.
+\bibitem[Baumann and Stohl(1997)]{baumann1997}
+Baumann, K. and Stohl, A.:
+Validation of a long-range trajectory model using gas balloon tracks from the Gordon Bennett Cup 95.
+J. Appl. Meteor., 36, 711--720, 1997.
+\bibitem[Beljaars and Holtslag(1991)]{beljaars1991}
+Beljaars, A. C. M. and Holtslag, A. A. M.:
+Flux parameterization over land surfaces for atmospheric models.
+J. Appl. Meteor., 30, 327--341, 1991.
+\bibitem[Belward et al.(1999)]{belward1999}
+Belward, A.S., Estes, J.E., and Kline, K.D.: 
+The IGBP-DIS 1-Km Land-Cover Data Set DISCover: A Project Overview. 
+Photogrammetric Engineering and Remote Sensing , v. 65, no. 9, p. 1013--1020, 1999.
+\bibitem[Berkowicz and Prahm(1982)]{berkowicz1982}
+Berkowicz, R. and Prahm, L. P.:
+Evaluation of the profile method for estimation of surface fluxes of momentum and heat.
+Atmos. Environ., 16, 2809--2819, 1982.
+\bibitem[Bey et al.(2001)]{bey2001}
+Bey I, Jacob DJ, Yantosca RM, et al.:
+Asian chemical outflow to the Pacific in spring: Origins, pathways, and budgets.
+J. Geophys. Res., 106, 19, 23073-23095, 2001.
+\bibitem[Businger et al.(1971)]{businger1971}
+Businger, J. A., Wyngaard, J. C., Izumi, Y. and Bradley, E. F.:
+Flux-profile relationships in the atmospheric surface layer.
+J. Atmos. Sci., 28, 181--189, 1971. 
+\bibitem[Dorling et al.(1992)]{dorling1992}
+ Dorling, S. R., Davies, T. D. and Pierce, C.E.:
+Cluster analysis: a technique for estimating the synoptic meteorological controls on air and precipitation chemistry - method and applications.
+Atmos. Environ., 26A, 2575--2581, 1992.
+\bibitem[ECMWF(1995)]{ecmwf1995}
+ECMWF,
+User Guide to ECMWF Products 2.1. Meteorological Bulletin M3.2. Reading, UK, 1995.
+\bibitem[Emanuel(1991)]{emanuel1991}
+Emanuel, K. A.:
+A scheme for representing cumulus convection in large-scale models,
+J. Atmos. Sci., 48, 2313--2335, 1991.
+\bibitem[Emanuel and \v{Z}ivkovi\'{c}-Rothman(1999)]{emanuel1999}
+Emanuel, K. A., and \v{Z}ivkovi\'{c}-Rothman, M.:
+Development and evaluation of a convection scheme for use in climate models.
+J. Atmos. Sci., 56, 1766--1782, 1999.
+\bibitem[Erisman et al.(1994)]{erisman1994}
+Erisman, J. W., Van Pul, A. and Wyers, P.:
+Parametrization of surface resistance for the quantification of atmospheric deposition of acidifying pollutants and ozone.
+Atmos. Environ., 28, 2595--2607, 1994.
+\bibitem[Flesch et al.(1995)]{flesch1995}
+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.
+\bibitem[Forster et al.(2001)]{forster2001}
+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\"ager, H., Stohl, A.:
+Transport of boreal forest fire emissions from Canada to Europe.
+J. Geophys. Res., 106, 22,887--22,906, 2001.
+\bibitem[Gupta et al.(1997)]{gupta1997}
+Gupta, S., McNider, R. T., Trainer, M., Zamora, R. J., Knupp, K. and Singh, M.P.:
+Nocturnal wind structure and plume growth rates due to inertial oscillations.
+J. Appl. Meteor., 36, 1050--1063, 1997.
+\bibitem[Hanna(1982)]{hanna1982}
+Hanna, S.R.:
+Applications in air pollution modeling. In: Nieuwstadt F.T.M. and H. van Dop (ed.): {\em Atmospheric Turbulence and Air Pollution Modelling}.
+D. Reidel Publishing Company, Dordrecht, Holland, 1982.
+\bibitem[Helmig et al.(2007)]{helmig2007}
+Helmig D., Ganzeveld, L., Butler, T. and Oltmans S.J.: 
+The role of ozone atmosphere-snow gas exchange on polar, boundary-layer tropospheric ozone - a review and sensitivity analysis 
+Atmos. Chem. and Physics, 7, 15--30, 2007.
+\bibitem[Hertel et al.(1995)]{hertel1995}
+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.
+\bibitem[Hubbe et al.(1997)]{hubbe1997}
+Hubbe, J. M., Doran, J. C., Liljegren, J.C. and Shaw, W. J.:
+Observations of spatial variations of boundary layer structure over the southern Great Plains cloud and radiation testbed.
+J. Appl. Meteor., 36, 1221--1231, 1997.
+\bibitem[Jacob and Wofsy(1990)]{jacob1990}
+Jacob D. J., Wofsy W. C.:
+Budgets of reactive nitrogen, hydrocarbons, and ozone over the amazon-forest during the wet season. 
+J. Geophys. Res., 95, 16737--16754, 1990. 
+\bibitem[Kahl (1996)]{kahl1996}
+Kahl, J. D. W.:
+On the prediction of trajectory model error.
+Atmos. Environ., 30, 2945--2957, 1996.
+\bibitem[Legg and Raupach(1982)]{legg1982}
+Legg, B. J., and Raupach, M. R.:
+Markov-chain simulation of particle dispersion in inhomogeneous flows: the mean drift velocity induced by a gradient in Eulerian velocity variance.
+Bound.-Layer Met., 24, 3--13, 1982.
+\bibitem[Legras et al.(2003)]{legras2003}
+Legras, B., Joseph, B. and Lefevre, F.:
+Vertical diffusivity in the lower stratosphere from Lagrangian back-trajectory reconstructions of ozone profiles.
+J. Geophys. Res., 108, 4562, doi:10.1029/2002JD003045, 2003.
+\bibitem[Maryon(1998)]{maryon1998}
+Maryon, R. H.:
+Determining cross-wind variance for low frequency wind meander.
+Atmos. Environ., 32, 115--121, 1998.
+\bibitem[McMahon(1979)]{mcmahon1979}
+McMahon, T. A., and Denison, P. J.:
+Empirical atmospheric deposition parameters - a survey.
+Atmos. Environ., 13, 571--585, 1979.
+\bibitem[McNider et al.(1988)]{mcnider1988}
+McNider, R. T., Moran, M. D. and Pielke, R. A.:
+Influence of diurnal and inertial boundary layer oscillations on long-range dispersion.
+Atmos. Environ., 22, 2445--2462, 1988.
+\bibitem[Naeslund and Thaning(1991)]{naeslund1991}
+Naeslund, E., and Thaning, L.:
+On the settling velocity in a nonstationary atmosphere.
+Aerosol Sci. Technol., 14, 247-256, 1991.
+\bibitem[Petterssen(1940)]{petterssen1940}
+Petterssen, S.:
+Weather Analysis and Forecasting.
+McGraw-Hill Book Company, New York. pp. 221--223, 1940.
+\bibitem[Pruppacher and Klett(1978)]{pruppacher1978}
+Pruppacher, H. R. and Klett, J. D.:
+Microphysics of clouds and precipitation.
+D. Reidel Publishing Company, Dordrecht, 714p., 1978.
+\bibitem[Rodean(1996)]{rodean1996}
+Rodean, H.:
+Stochastic Lagrangian models of turbulent diffusion.
+Meteorological Monographs, 26 (48). American Meteorological Society, Boston, USA, 1996.
+\bibitem[Ryall et al.(1997)]{ryall1997}
+Ryall, D. B. and Maryon, R. H.:
+Validation of the UK Met Office's NAME model against the ETEX dataset. In: Nodop, K. (editor): ETEX Symposium on Long-Range Atmospheric Transport, Model Verification and Emergency Response, European Commission EUR 17346, 151--154, 1997.
+\bibitem[Seibert(2001)]{seibert2001}
+Seibert, P.:
+Inverse modelling with a Lagrangian particle dispersion model:  application to point releases over limited time intervals. In: Gryning, S. E., Schiermeier, F.A.  (eds.): Air Pollution Modeling and its Application XIV. Proc. of ITM Boulder. New York: Plenum Press, 381--389, 2001.
+\bibitem[Seibert and Frank(2004)]{seibert2004}
+Seibert, P. and Frank, A.:
+Source-receptor matrix calculation with a Lagrangian particle dispersion model in backward mode.
+Atmos. Chem. Phys., 4, 51--63, 2004.
+\bibitem[Seibert et al.(2001)]{seibertetal2001}
+Seibert, P., Kr\"uger, B. and Frank, A.:
+Parametrisation of convective mixing in a Lagrangian particle dispersion model, Proceedings of the 5th GLOREAM Workshop, Wengen, Switzerland, September 24--26, 2001.
+\bibitem[Slinn(1982)]{slinn1982}
+Slinn, W. G. N.:
+Predictions for particle deposition to vegetative canopies.
+Atmos. Environ., 16, 1785--1794, 1982.
+\bibitem[Spichtinger et al.(2001)]{spichtinger2001}
+Spichtinger, N., Wenig, M., James, P., Wagner, T., Platt, U. and Stohl, A.:
+Satellite detection of a continental-scale plume of nitrogen oxides from boreal forest fires,
+Geophys. Res. Lett., 28, 4579--4582, 2001.
+\bibitem[Stohl(1998)]{stohl1998}
+Stohl, A.:
+Computation, accuracy and applications of trajectories -- a review and bibliography.
+Atmos. Environ., 32, 947--966, 1998.
+\bibitem[Stohl et al.(2005)]{stohl2005} 
+Stohl, A., Forster, C., Frank, A., Seibert, P., Wotawa, G.: 
+Technical note: The Lagrangian particle dispersion model FLEXPART version 6.2., 
+Atmos. Chem. Phys., 5, 2461-2474, 2005
+\bibitem[Stohl et al.(2004)]{stohl2004}
+Stohl, A., Cooper, O. R., Damoah, R., Fehsenfeld, F. C., Forster, C., Hsie, E.-Y., Hübler, G., Parrish,
+D. D. and Trainer, M.:
+Forecasting for a Lagrangian aircraft campaign.
+Atmos. Chem. Phys., 4, 1113--1124, 2004.
+\bibitem[Stohl et al.(2002)]{stohl2002}
+Stohl, A., Eckhardt, S., Forster, C., James, P., Spichtinger, N. and Seibert, P.:
+A replacement for simple back trajectory calculations in the interpretation of atmospheric trace substance measurements.
+Atmos. Environ., 36, 4635--4648, 2002.
+\bibitem[Stohl et al.(2003)]{stohl2003}
+Stohl, A., Forster, C., Eckhardt, S., Spichtinger, N., Huntrieser, H., Heland, J., Schlager, H., Wilhelm, S., Arnold, F. and Cooper, O.:
+A backward modeling study of intercontinental pollution transport using aircraft measurements.
+J. Geophys. Res., 108, 4370, doi:10.1029/2002JD002862, 2003.
+\bibitem[Stohl et al.(1998)]{stohletal1998}
+Stohl, A., Hittenberger, M. and Wotawa, G.:
+Validation of the Lagrangian particle dispersion model FLEXPART against large scale tracer experiment data.
+Atmos. Environ., 32, 4245--4264, 1998.
+\bibitem[Stohl et al.(2000)]{stohl2000}
+Stohl, A., Spichtinger-Rakowsky, N., Bonasoni, P., Feldmann, H., Memmesheimer, M., Scheel, H. E., Trickl, T., H\"ubener, S. H., Ringer, W. and Mandl, M.:
+The influence of stratospheric intrusions on alpine ozone concentrations.
+Atmos. Environ., 34, 1323--1354, 2000. 
+\bibitem[Stohl and Thomson(1999)]{stohlthomson1999}
+Stohl, A. and Thomson, D. J.:
+A density correction for Lagrangian particle dispersion models.
+Bound.-Layer Met., 90, 155--167, 1999.
+\bibitem[Stohl and Trickl(1999)]{stohl1999}
+Stohl, A. and Trickl, T.:
+A textbook example of long-range transport: Simultaneous observation of ozone maxima of stratospheric and North American origin in the free troposphere over Europe,
+J. Geophys. Res., 104, 30,445--30,462, 1999.
+\bibitem[Walmsley and Wesely(1996]{walmsley1996}
+Walmsley J.L., Wesely M.L.: 
+Modification of coded parametrizations of surface resistances to gaseous dry deposition.  
+Atmos. Environ., 30, 1181--1188, 1996. 
+\bibitem[Stohl et al.(1995)]{stohl1995} 
+Stohl, A., Wotawa, G., Seibert, P. and Kromp-Kolb, H.:
+Interpolation errors in wind fields as a function of spatial and temporal resolution and their impact on different types of kinematic trajectories.
+J. Appl. Meteor., 34, 2149--2165, 1995.
+\bibitem[Stull(1988)]{stull1988}
+Stull, R. B.:
+An Introduction to Boundary Layer Meteorology. Kluwer Academic Publishers, Dordrecht, 1988.
+\bibitem[Telford(1975)]{telford1975}
+Telford, J. W.:
+Turbulence, entrainment and mixing in cloud dynamics,
+Pure Appl. Geophys., 113, 1067--1084, 1975.
+\bibitem[Thomson(1987)]{thomson1987}
+Thomson D. J.:
+Criteria for the selection of stochastic models of particle trajectories in turbulent flows.
+J. Fluid Mech., 180, 529--556, 1987.
+\bibitem[Uliasz(1994)]{uliasz1994}
+Uliasz, M.:
+Lagrangian particle dispersion modeling in mesoscale applications. In: Zannetti, P. (ed.): {\it Environmental Modeling, Vol. II}. Computational Mechanics Publications, Southampton, UK, 1994.
+\bibitem[Velde et al.(1994)]{velde1994}
+Velde van de, R. J., Faber, W. S., Katwijk van, V. F., Kuylenstierna, J. C. I., Scholten., H. J., Thewessen, T. J. M., Verspuij, M. and Zevenbergen, M.:
+The Preparation of a European Land Use Database. National Institute of Public Health and Environmental Protection, Report nr 712401001, Bilthoven, The Netherlands, 1994.
+\bibitem[Vogelezang and Holtslag(1996)]{vogelezang1996}
+Vogelezang, D. H. P. and Holtslag, A. A. M.:
+Evaluation and model impacts of alternative boundary-layer height formulations.
+Bound.-Layer Met., 81, 245--269, 1996.
+\bibitem[Wesely and Hicks(1977)]{wesely1977}
+Wesely, M. L. and Hicks, B. B.:
+Some factors that affect the deposition rates of sulfur dioxide and similar gases on vegetation.
+J. Air Poll. Contr. Assoc., 27, 1110--1116, 1977.
+\bibitem[Wesely(1989)]{wesely1989}
+Wesely, M. L.:
+Parameterization of surface resistances to gaseous dry deposition in regional-scale numerical models.
+Atmos. Environ., 23, 1293--1304, 1989.
+\bibitem[Wilson and Flesch(1993)]{wilson1993}
+Wilson, D. J. and Flesch, T. K.:
+Flow boundaries in random-flight dispersion models: enforcing the well-mixed condition.
+J. Appl. Meteor., 32, 1695--1707, 1993.
+\bibitem[Wilson et al.(1983)]{wilson1983}
+Wilson, J. D., Legg, B. J. and Thomson, D. J.:
+Calculation of particle trajectories in the presence of a gradient in turbulent-velocity scale.
+Bound.-Layer Met., 27, 163--169, 1983.
+\bibitem[Wilson and Sawford(1996)]{wilson1996}
+Wilson, J. D. and Sawford, B. L.:
+Review of Lagrangian stochastic models for trajectories in the turbulent atmosphere.
+Bound.-Layer Met., 78, 191--210, 1996.
+\bibitem[Wotawa et al.(1996)]{wotawa1996}
+Wotawa, G., Stohl, A. and Kromp-Kolb, H.:
+Parameterization of the planetary boundary layer over Europe - a data comparison between the observation based OML preprocessor and ECMWF model data.
+Contr. Atmos. Phys., 69, 273--284, 1996.
+\bibitem[Wotawa and Stohl(1997)]{wotawa1997}
+Wotawa, G. and Stohl, A.:
+Boundary layer heights and surface fluxes of momentum and heat derived from ECMWF data for use in pollutant dispersion models -- problems with data accuracy.
+In: Gryning, S.-E., Beyrich, F. and E. Batchvarova (editors): The Determination of the Mixing Height -- Current Progress and Problems. EURASAP Workshop Proceedings, 1-3 October 1997, Ris{\o} National Laboratory, Denmark, 1997.
+\bibitem[Zannetti(1992)]{zannetti1992}
+Zannetti, P.:
+Particle Modeling and Its Application for Simulating Air Pollution Phenomena. In: Melli P. and P. Zannetti (ed.): Environmental Modelling. Computational Mechanics Publications, Southampton, UK, 1992.
+\bibitem[Zhang et al.(2002)]{zhang2002}
+Zhang L.M., Moran M.D., Makar P.A., Brook J.R. and S. Gong
+Modelling gaseous dry deposition in AURAMS: a unified regional air-quality modelling system 
+Atmos. Environ., 36, 537-560, 2002.
+
+\end{thebibliography}
+
+\newpage
+
+\appendix
+
+\onecolumn
+
+\section{FLEXPART sample input files}
+\subsection{The pathnames file}
+A file \verb|pathnames| must exist in the directory where FLEXPART is started.
+It states the pathnames (absolute or relative) of input and output files:
+\begin{footnotesize}\begin{verbatim}
+/home/as/FLEXPART50/options/
+/volc/as/contrace/modelresults/forward/
+/volc/windcontrace/
+/volc/windcontrace/AVAILABLE
+/volc/nested/
+/volc/nested/AVAILABLE
+============================================
+
+Line 1: path where control files "COMMAND" and "RELEASES" are available
+Line 2: name of directory where output files are generated
+Line 3: path where meteorological fields are available (mother grid)
+Line 4: full filename of "AVAILABLE"-file (mother grid)
+
+Subsequent lines:
+Line 2n+3: path where meteorological fields are available (nested grid n)
+Line 2n+4: full filename of "AVAILABLE"-file (nested grid n)
+
+Line below last pathname must be:
+============================================
+
+The grids must be arranged such as that the coarse-scale nests
+come before the fine-scale nests. Multiple nests of the same
+nesting level are allowed. In that case, the order is arbitrary.
+\end{verbatim}\end{footnotesize}
+
+\newpage
+
+\subsection{Files in directory windfields}
+
+The directory where the meteorological input data are stored, here called
+\verb|windfields| (\verb|/volc/windcontrace/| in the above example
+\verb|pathnames| file), contains grib-code files containing the ECMWF data.
+All meteorological fields must have the same structure, i.e.  the same
+computational domain and the same resolution.  An example listing of this
+directory is given below.
+\begin{footnotesize}\begin{verbatim}
+AVAILABLE       EN01102806      EN01102815
+EN01102800      EN01102809      EN01102818
+EN01102803      EN01102812      EN01102821
+\end{verbatim}\end{footnotesize}
+
+The file names of the grib-code files and their validation dates and times (in UTC) must be listed in the file \verb|AVAILABLE|. While it is practical to have this file reside in the same directory as the wind fields, this is no necessity and it can also be located elsewhere, as its file name is also given in the \verb|pathnames| file.
+\begin{footnotesize}\begin{verbatim}
+DATE     TIME         FILENAME     SPECIFICATIONS
+YYYYMMDD HHMISS
+________ ______      __________      __________
+20011028 000000      EN01102800      ON DISC
+20011028 030000      EN01102803      ON DISC
+20011028 060000      EN01102806      ON DISC
+20011028 090000      EN01102809      ON DISC
+20011028 120000      EN01102812      ON DISC
+20011028 150000      EN01102815      ON DISC
+20011028 180000      EN01102818      ON DISC
+20011028 210000      EN01102821      ON DISC
+\end{verbatim}\end{footnotesize}
+
+Since version 7.0, FILENAME can be up to 256 characters long.  Thus, the path
+to a file can be directly placed in the AVAILABLE file.  If this feature is
+used, line 3 in the \verb|pathnames| file should be left empty.
+\begin{footnotesize}\begin{verbatim}
+DATE     TIME         FILENAME     SPECIFICATIONS
+YYYYMMDD HHMISS
+________ ______      __________      __________
+20011028 000000      /work_disk/data/EN01102800      ON DISC
+20011028 030000      /work_disk/data/EN01102803      ON DISC
+20011028 060000      /work_disk/data/EN01102806      ON DISC
+20011028 090000      /work_disk/data/EN01102809      ON DISC
+20011028 120000      /work_disk/data/EN01102812      ON DISC
+20011028 150000      /work_disk/data/EN01102815      ON DISC
+20011028 180000      /work_disk/data/EN01102818      ON DISC
+20011028 210000      /work_disk/data/EN01102821      ON DISC
+\end{verbatim}\end{footnotesize}
+
+Nested wind fields must be stored in one or more different directory/ies, as specified in the \verb|pathnames| file.
+
+\newpage
+
+\subsection{Files in directory options}
+
+The files in directory \begin{footnotesize}\verb|options|\end{footnotesize} are used to specify the model run.
+An example listing of \begin{footnotesize}\verb|options|\end{footnotesize} is given below.
+\begin{footnotesize}\begin{verbatim}
+AGECLASSES           IGBP_int1.dat       RECEPTORS             SPECIES
+COMMAND              OH_7lev_agl.dat     RELEASES              surfdata.t
+COMMAND.alternative  OUTGRID             RELEASES.alternative  surfdepo.t
+COMMAND.reference    OUTGRID_NEST        RELEASES.reference              
+\end{verbatim}\end{footnotesize}
+
+Here, SPECIES is a subdirectory containing a number of files (this feature was introduced in version 8.0):
+\begin{footnotesize}\begin{verbatim}
+SPECIES:
+SPECIES_001  SPECIES_004  SPECIES_007
+SPECIES_002  SPECIES_005  SPECIES_008
+SPECIES_003  SPECIES_006  SPECIES_009
+\end{verbatim}\end{footnotesize}
+
+\subsubsection{File COMMAND}
+
+The most important configuration file for setting up a FLEXPART simulation is
+the \verb|COMMAND| file which specifies (1) the simulation direction (either
+forward or backward), (2) the start and (3) the end time of the simulation, (4)
+the frequency $T_c$ of the model output, (5) the averaging time $\Delta T_c$ of
+model output, and (6) the intervals $\Delta T_s$ at which concentrations are
+sampled, (7) the time constant for particle splitting $\Delta t_s$, (8) the
+synchronisation interval of FLEXPART, (9) the factor $c_{tl}$ by which the time
+steps must be smaller than the Lagrangian time scale, and (10) the refinement
+factor for the time step used for solving the Langevin equation of the vertical
+component of the turbulent wind.  If (9) ($c_{tl}$) is negative, the Langevin
+equations are solved with constant time steps according to the synchronisation
+interval.  In that case, the value of (10) is arbitrary.  The synchronisation
+interval is the minimum time interval used by the model for all activities
+(such as concentration calculations, wet deposition calculations, interpolation
+of data, mesoscale wind fluctuations or output of data) other than the
+simulation of turbulent transport and dry deposition (if ($c_{tl}>0$).  Further
+switches determine (11) whether concentrations, mixing ratios, residence times
+or plume trajectories (or combinations thereof) are to be calculated, (12) the
+option of particle position dump either at the end of or continuously during
+the simulation, (13) on/off of subgrid terrain effect parameterization, (14)
+on/off of deep convection parameterization, (15) on/off calculation of age
+spectra, (16) continuation of simulation from previous particle dump, (17)
+write output for each RELEASE location on/off, (18) on/off for mass flux
+calculations and output, (19) on/off for the domain-filling option of FLEXPART,
+(20) an indicator that determines whether mass or mass mixing ratio units are
+to be used at the source, (21) an indicator that determines whether mass or
+mass mixing ratio units are to be used at the receptor, (22) on/off of
+additional compact dump of the positions of numbered particles, (23) on/off for
+the use of nested output fields, (24) write sensitivity to initial conditions
+in backward mode off/mass units/mass mixing ratio units.
+
+Two versions of \verb|COMMAND| may be used, which both can be read in by
+FLEXPART: the first contains formatted input (i.e., a mask to be filled for the
+various input options that must be filled in), the second contains largely
+unformatted input and is recommended for the more experienced FLEXPART user.
+The following example is for formatted input.
+
+\begin{scriptsize}\begin{verbatim}
+********************************************************************************
+*                                                                              *
+*      Input file for the Lagrangian particle dispersion model FLEXPART        *
+*                           Please select your options                         *
+*                                                                              *
+********************************************************************************
+
+1. __                3X, I2
+    1       
+   LDIRECT           1 FOR FORWARD SIMULATION, -1 FOR BACKWARD SIMULATION
+
+2. ________ ______   3X, I8, 1X, I6
+   20040626 000000
+   YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION
+
+3. ________ ______   3X, I8, 1X, I6
+   20040816 120000
+   YYYYMMDD HHMISS   ENDING DATE OF SIMULATION
+
+4. _____             3X, I5
+    7200 
+   SSSSS             OUTPUT EVERY SSSSS SECONDS
+
+5. _____             3X, I5
+    7200 
+   SSSSS             TIME AVERAGE OF OUTPUT (IN SSSSS SECONDS)
+
+6. _____             3X, I5
+     900 
+   SSSSS             SAMPLING RATE OF OUTPUT (IN SSSSS SECONDS)
+
+7. _________         3X, I9
+   999999999
+   SSSSSSSSS         TIME CONSTANT FOR PARTICLE SPLITTING (IN SECONDS)
+
+8. _____             3X, I5
+     900 
+   SSSSS             SYNCHRONISATION INTERVAL OF FLEXPART (IN SECONDS)
+
+9.  ---.--           4X, F6.4
+     -5.0
+    CTL              FACTOR, BY WHICH TIME STEP MUST BE SMALLER THAN TL
+
+10. ---              4X, I3
+      4
+    IFINE            DECREASE OF TIME STEP FOR VERTICAL MOTION BY FACTOR IFINE
+
+11. -                4X, I1
+    3  
+    IOUT             1 CONCENTRATION (RESIDENCE TIME FOR BACKWARD RUNS) OUTPUT, 2 MIXING RATIO OUTPUT, 3 BOTH,4 PLUME TRAJECT., 5=1+4
+
+12. -                4X, I1
+    2  
+    IPOUT            PARTICLE DUMP: 0 NO, 1 EVERY OUTPUT INTERVAL, 2 ONLY AT END
+
+13. _                4X, I1
+    1
+    LSUBGRID         SUBGRID TERRAIN EFFECT PARAMETERIZATION: 1 YES, 0 NO
+
+14. _                4X, I1
+    1
+    LCONVECTION      CONVECTION: 1 YES, 0 NO
+
+15. _                4X, I1
+    0
+    LAGESPECTRA      AGE SPECTRA: 1 YES, 0 NO
+
+16. _                4X, I1
+    0
+    IPIN             CONTINUE SIMULATION WITH DUMPED PARTICLE DATA: 1 YES, 0 NO
+
+17. _
+    0                4X,I1
+    IOFR             IOUTPUTFOREACHREL CREATE AN OUTPUT FILE FOR EACH RELEASE LOCATION: 1 YES, 0 NO
+
+18. _                4X, I1
+    0
+    IFLUX            CALCULATE FLUXES: 1 YES, 0 NO
+
+19. _                4X, I1
+    0
+    MDOMAINFILL      DOMAIN-FILLING TRAJECTORY OPTION: 1 YES, 0 NO, 2 STRAT. O3 TRACER
+
+20. _                4X, I1
+    1
+    IND_SOURCE       1=MASS UNIT , 2=MASS MIXING RATIO UNIT
+
+21. _                4X, I1
+    1
+    IND_RECEPTOR     1=MASS UNIT , 2=MASS MIXING RATIO UNIT
+
+22. _                4X, I1
+    0
+    MQUASILAG        QUASILAGRANGIAN MODE TO TRACK INDIVIDUAL PARTICLES: 1 YES, 0 NO
+
+23. _                4X, I1
+    0
+    NESTED_OUTPUT    SHALL NESTED OUTPUT BE USED? 1 YES, 0 NO
+24. _                4X, I1
+    2
+    LINIT_COND       INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
+    
+
+1. Simulation direction, 1 for forward, -1 for backward in time
+        (consult Seibert and Frank, 2004 for backward runs)
+
+2. Beginning date and time of simulation. Must be given in format
+   YYYYMMDD HHMISS, where YYYY is YEAR, MM is MONTH, DD is DAY, HH is HOUR,
+   MI is MINUTE and SS is SECOND. Current version utilizes UTC.
+
+3. Ending date and time of simulation. Same format as 3.
+
+4. Average concentrations are calculated every SSSSS seconds.
+
+5. The average concentrations are time averages of SSSSS seconds
+   duration. If SSSSS is 0, instantaneous concentrations are outputted.
+
+6. The concentrations are sampled every SSSSS seconds to calculate the time
+   average concentration. This period must be shorter than the averaging time.
+
+7. Time constant for particle splitting. Particles are split into two
+   after SSSSS seconds, 2xSSSSS seconds, 4xSSSSS seconds, and so on.
+
+8. All processes are synchronized with this time interval (lsynctime).
+   Therefore, all other time constants must be multiples of this value.
+   Output interval and time average of output must be at least twice lsynctime.
+
+9. CTL must be >1 for time steps shorter than the  Lagrangian time scale
+   If CTL<0, a purely random walk simulation is done
+
+10.IFINE=Reduction factor for time step used for vertical wind
+
+11.IOUT determines how the output shall be made: concentration
+   (ng/m3, Bq/m3), mixing ratio (pptv), or both, or plume trajectory mode,
+   or concentration + plume trajectory mode.
+   In plume trajectory mode, output is in the form of average trajectories.
+
+12.IPOUT determines whether particle positions are outputted (in addition
+   to the gridded concentrations or mixing ratios) or not.
+   0=no output, 1 output every output interval, 2 only at end of the
+   simulation
+
+13.Switch on/off subgridscale terrain parameterization (increase of
+   mixing heights due to subgridscale orographic variations)
+
+14.Switch on/off the convection parameterization
+
+15.Switch on/off the calculation of age spectra: if yes, the file AGECLASSES
+   must be available
+
+16. If IPIN=1, a file "partposit_end" from a previous run must be available in
+    the output directory. Particle positions are read in and previous simulation
+    is continued. If IPIN=0, no particles from a previous run are used
+
+17. IF IOUTPUTFOREACHRELEASE is set to 1, one output field for each location
+    in the RELEASES file is created. For backward calculation this should be
+    set to 1. For forward calculation both possibilities are applicable.
+
+18. If IFLUX is set to 1, fluxes of each species through each of the output
+    boxes are calculated. Six fluxes, corresponding to northward, southward,
+    eastward, westward, upward and downward are calculated for each grid cell of
+    the output grid. The control surfaces are placed in the middle of each
+    output grid cell. If IFLUX is set to 0, no fluxes are determined.
+
+19. If MDOMAINFILL is set to 1, the first box specified in file RELEASES is used
+    as the domain where domain-filling trajectory calculations are to be done.
+    Particles are initialized uniformly distributed (according to the air mass
+    distribution) in that domain at the beginning of the simulation, and are
+    created at the boundaries throughout the simulation period.
+
+20. IND_SOURCE switches between different units for concentrations at the source
+    NOTE that in backward simulations the release of computational particles
+    takes place at the "receptor" and the sampling of particles at the "source".
+          1=mass units (for bwd-runs = concentration)
+          2=mass mixing ratio units
+21. IND_RECEPTOR switches between different units for concentrations at the receptor
+          1=mass units (concentrations)
+          2=mass mixing ratio units
+
+22. MQUASILAG indicates whether particles shall be numbered consecutively (1) or
+    with their release location number (0). The first option allows tracking of
+    individual particles using the partposit output files
+
+23. NESTED_OUTPUT decides whether model output shall be made also for a nested
+    output field (normally with higher resolution)
+
+24. LINIT_COND determines whether, for backward runs only, the sensitivity to initial
+    conditions shall be calculated and written to output files
+    0=no output, 1 or 2 determines in which units the initial conditions are provided.
+
+\end{verbatim}\end{scriptsize}
+
+\newpage
+
+\subsubsection{File OUTGRID}
+The file \verb|OUTGRID| specifies the output grid.
+\begin{scriptsize}\begin{verbatim}
+********************************************************************************
+*                                                                              *
+*      Input file for the Lagrangian particle dispersion model FLEXPART        *
+*                       Please specify your output grid                        *
+*                                                                              *
+********************************************************************************
+
+1.  ------.----       4X,F11.4
+       -10.0000       GEOGRAFICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
+    OUTLONLEFT        (left boundary of the first grid cell - not its centre)
+
+
+2.  ------.----       4X,F11.4
+        40.0000       GEOGRAFICAL LATITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
+    OUTLATLOWER       (lower boundary of the first grid cell - not its centre)
+
+3.  -----             4X,I5
+      101             NUMBER OF GRID POINTS IN X DIRECTION (= No. of cells + 1)
+    NUMXGRID
+
+4.  -----             4X,I5
+       47             NUMBER OF GRID POINTS IN Y DIRECTION (= No. of cells + 1)
+    NUMYGRID
+
+5.  ------.---        4X,F10.3
+         0.500        GRID DISTANCE IN X DIRECTION
+    DXOUTLON
+
+6.  ------.---        4X,F10.3
+         0.500        GRID DISTANCE IN Y DIRECTION
+    DYOUTLAT
+
+7.  -----.-           4X, F7.1
+      100.0
+    LEVEL 1           HEIGHT OF LEVEL (UPPER BOUNDARY)
+
+8.  -----.-           4X, F7.1
+      300.0
+    LEVEL 2           HEIGHT OF LEVEL (UPPER BOUNDARY)
+
+9.  -----.-           4X, F7.1
+      600.0
+    LEVEL 3           HEIGHT OF LEVEL (UPPER BOUNDARY)
+
+10. -----.-           4X, F7.1
+     1000.0
+    LEVEL 4           HEIGHT OF LEVEL (UPPER BOUNDARY)
+
+11. -----.-           4X, F7.1
+     2000.0
+    LEVEL 5           HEIGHT OF LEVEL (UPPER BOUNDARY)
+
+12. -----.-           4X, F7.1
+     3000.0
+    LEVEL 6           HEIGHT OF LEVEL (UPPER BOUNDARY)
+\end{verbatim}\end{scriptsize}
+In order to define the grid for a nested output field, the file
+\verb|OUTGRID_NEST| must exist.  It has the same format as file \verb|OUTGRID|,
+but does not contain the vertical level information:
+
+\begin{scriptsize}\begin{verbatim}
+********************************************************************************
+*                                                                              *
+*      Input file for the Lagrangian particle dispersion model FLEXPART        *
+*                       Please specify your output grid                        *
+*                                                                              *
+********************************************************************************
+
+1.  ------.----       4X,F11.4
+      -125.0000       GEOGRAFICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
+    OUTLONLEFT        (left boundary of the first grid cell - not its centre)
+
+2.  ------.----       4X,F11.4
+        25.0000       GEOGRAFICAL LATITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
+    OUTLATLOWER       (lower boundary of the first grid cell - not its centre)
+
+3.  -----             4X,I5
+        1             NUMBER OF GRID POINTS IN X DIRECTION (= No. of cells + 1)
+    NUMXGRID
+
+4.  -----             4X,I5
+        1             NUMBER OF GRID POINTS IN Y DIRECTION (= No. of cells + 1)
+    NUMYGRID
+
+5.  ------.-----      4X,F12.5
+         0.33333      GRID DISTANCE IN X DIRECTION
+    DXOUTLON
+
+6.  ------.-----      4X,F12.5
+         0.25000      GRID DISTANCE IN Y DIRECTION
+    DYOUTLAT
+\end{verbatim}\end{scriptsize}
+
+\newpage
+
+\subsubsection{File RECEPTORS}
+
+\verb|RECEPTORS| specifies the receptor locations for which the parabolic
+kernel method shall be applied to calculate air concentrations.  The maximum
+number of receptor sites is set by parameter \verb|maxreceptor| in file
+\verb|includepar|.
+\begin{scriptsize}\begin{verbatim}
+********************************************************************************
+*                                                                              *
+*      Input file for the Lagrangian particle dispersion model FLEXPART        * 
+*                       Please specify your receptor points                    * 
+*     For the receptor points, ground level concentrations are calculated      * 
+*                                                                              *
+********************************************************************************
+1.  ----------------  4X,A16
+    F15               NAME OF RECEPTOR POINT
+    RECEPTORNAME
+
+2.  ------.----       4X,F11.4
+         6.1333       GEOGRAFICAL LONGITUDE
+    XRECEPTOR
+
+3.  ------.----       4X,F11.4
+        49.0833       GEOGRAFICAL LATITUDE
+    YRECEPTOR
+================================================================================
+1.  ----------------  4X,A16
+    NL01              NAME OF RECEPTOR POINT
+    RECEPTORNAME
+
+2.  ------.----       4X,F11.4
+         5.7833       GEOGRAFICAL LONGITUDE
+    XRECEPTOR
+
+3.  ------.----       4X,F11.4
+        50.9167       GEOGRAFICAL LATITUDE
+    YRECEPTOR
+================================================================================
+\end{verbatim}\end{scriptsize}
+
+\newpage
+
+\subsubsection{File RELEASES}
+
+\begin{footnotesize}\verb|RELEASES|\end{footnotesize} defines the release
+specifications.  In the first input line, the number $N$ of emitted species is
+defined (1 in the example below).  At all locations, the same species must be
+released.  The next $N$ input lines give a cross-reference to the respective file
+\begin{footnotesize}\verb|SPECIES_nnn|\end{footnotesize}, where the physical and
+chemical properties of the released species are given (also the temporal
+variations of emissions is defined for each species).  Then follows a list of
+release sites, for each of which the release characteristics must be entered:
+the beginning and the ending time of the release, geographical coordinates of
+the lower left and upper right corners of the release location, type of
+vertical coordinate (above ground level, or above sea level), lower level and
+upper level of source box, the number of particles to be used, and the total
+mass emitted.  Note that the mass entry must be repeated $N$ times, one mass
+per species released.  Finally, a name is assigned to each release point.\par
+
+The particles are released from random locations within a four-dimensional box
+extending from the lower to the upper level above a rectangle (on a lat/lon
+grid) defined by the geographical coordinates, and between the release's start
+and end.  With some of the coordinates set identically, line or point sources
+can be specified.
+
+As for \verb|COMMAND|, the \verb|RELEASES| file can be provided formatted or
+unformatted.  The example below shows the formatted version.
+
+\begin{scriptsize}\begin{verbatim}
+*************************************************************************
+*                                                                       *
+*                                                                       *
+*                                                                       *
+*   Input file for the Lagrangian particle dispersion model FLEXPART    *
+*                        Please select your options                     *
+*                                                                       *
+*                                                                       *
+*                                                                       *
+*************************************************************************
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  1      
+___                        i3    Total number of species emitted
+
+ 24
+___                        i3    Index of species in file SPECIES
+
+=========================================================================
+20011028  150007
+________ ______            i8,1x,i6 Beginning date and time of release
+
+20011028  150046
+________ ______            i8,1x,i6 Ending date and time of release
+
+   9.4048 
+____.____                  f9.4  Longitude [DEG] of lower left corner
+
+  48.5060
+____.____                  f9.4  Latitude [DEG] of lower left corner
+
+   9.5067
+____.____                  f9.4  Longitude [DEG] of upper right corner
+
+  48.5158
+____.____                  f9.4  Latitude [DEG] of upper right corner
+
+        2
+_________                  i9    1 for m above ground, 2 for m above sea level, 3 for pressure in hPa
+
+ 6933.60
+_____.___                  f10.3 Lower z-level (in m agl or m asl)
+ 
+ 6950.40
+_____.___                  f10.3 Upper z-level (in m agl or m asl)
+ 
+    20000                
+_________                  i9    Total number of particles to be released
+
+1.0000E00
+_.____E__                  e9.4  Total mass emitted
+
+FLIGHT_11242
+________________________________________   character*40 comment
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+20011028  150047
+________ ______            i8,1x,i6 Beginning date and time of release
+
+20011028  150107
+________ ______            i8,1x,i6 Ending date and time of release
+
+   9.3038 
+____.____                  f9.4  Longitude [DEG] of lower left corner
+
+  48.5158
+____.____                  f9.4  Latitude [DEG] of lower left corner
+
+   9.4048
+____.____                  f9.4  Longitude [DEG] of upper right corner
+
+  48.5906
+____.____                  f9.4  Latitude [DEG] of upper right corner
+
+        2
+_________                  i9    1 for m above ground, 2 for m above sea level, 3 for pressure in hPa
+
+ 6833.50
+_____.___                  f10.3 Lower z-level (in m agl or m asl)
+ 
+ 6950.40
+_____.___                  f10.3 Upper z-level (in m agl or m asl)
+ 
+    20000                
+_________                  i9    Total number of particles to be released
+
+1.0000E00
+_.____E__                  e9.4  Total mass emitted
+
+FLIGHT_11185
+________________________________________   character*40 comment
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+\end{verbatim}\end{scriptsize}
+
+\newpage
+
+\subsubsection{File AGECLASSES}
+
+\verb|AGECLASSES| provides the times for the age class calculation.  In the
+first data line, the number $n$ of age classes is set, and ages are listed in
+the following $n$ lines.  The entries specify the end times (in seconds) of the
+respective intervals to be used, the first one starting at zero seconds.
+Particles are dropped from the simulation once they exceed the maximum age.
+Even if no age classes are needed, this option (with the number of age classes
+set to 1) can be useful to determine the age at which particles are removed
+from the simulation.
+
+\begin{footnotesize}\begin{verbatim}
+************************************************
+*                                              *
+*Lagrangian particle dispersion model FLEXPART *
+*         Please select your options           *
+*                                              *
+*This file determines the ageclasses to be used*
+*                                              *
+*Ages are given in seconds. The first class    *
+*starts at age zero and goes up to the first   *
+*age specified. The last age gives the maximum *
+*time a particle is carried in the simulation. *
+*                                              *
+************************************************
+ 6          Integer        Number of age classes
+43200       Integer        Age class 1
+86400       Integer        Age class 2
+129600
+172800
+259200
+345600
+\end{verbatim}\end{footnotesize}
+
+\newpage
+
+\subsubsection{Files SPECIES\_nnn and surfdata.t}
+
+\begin{footnotesize}\verb|SPECIES_nnn|\end{footnotesize} (where \verb|nnn| is 
+a zero-padded identifier of the species number) specifies all physico-chemical
+properties for the given species.  Entries are the half life (due to
+radioactive or chemical decay), wet deposition information (\verb|A| and \verb|B| are the
+factors defined by Eq.~\ref{wetscav}), dry deposition information for
+gases ($D=D_{H_2O}/D_i$, $D_{H_2O}$ is the diffusivity of water vapor and \verb|D_i|
+is the diffusivity of the species, \verb|H| is the effective Henry's constant, and
+\verb|f0| varies between 0 and 1 and gives the reactivity of a species relative to
+that of ozone.  For nonreactive species \verb|f0| is 0, for slightly reactive it is
+0.1 and for highly reactive it is 1.), dry deposition information for
+particulates (\verb|rho| specifies the density of the substance, \verb|dquer| its mean
+diameter $\overline{d_p}$, and \verb|dsig| the measure of variation $\sigma_p$).
+Radioactive decay is switched off by specifying a negative half life, wet
+deposition is switched off by specifying negative \verb|A|, dry deposition of gases
+is switched off by negative \verb|D|, dry deposition of particles is switched off by
+negative \verb|rho|.  If no detailed information for deposition velocity calculation
+is available, a constant deposition velocity \verb|vd| (cm s$^{-1}$) can be used.
+\verb|molweight| gives the molecular weight of the species, which is needed for
+mixing ratio output.  For degradation in an monthly averaged OH field
+\citep{bey2001} the OH Reaction ratio at 25\degreee C for the species can be given,
+units are (cm$^3$/s).
+
+Optionally an emission variation information can be added at the end of the
+file, defined as following: Since FLEXPART version 6.0, emission factors can be
+defined that change the temporal variation of particle releases.  This is
+useful, for instance, to simulate the typical daily and weekly cycle of
+anthropogenic emissions.  The emission factors are given in the file of the
+corresponding species \verb|SPECIES_nnn|, where \verb|nnn| is the species
+number defined in file \verb|RELEASES|.  If no emission variation information
+is given, emission rates for species \verb|nnn| are taken as constant.  Release
+rates can vary with the hour of the day and with the day of the week, according
+to the local time at the release location.  Emission factors must be 1 on
+average.  24 hourly as well as 7 daily values must be specified.  Furthermore,
+different disaggregation factors must be given for area sources and for point
+sources.  FLEXPART distinguishes between the two using the lower altitude of
+the release box: area sources are assumed to start below 0.5 m above the
+ground, whereas point sources are assumed to be higher.  Please note that when
+this option is used, it is not so easy to determine the maximum number of
+particles present at a particular time of the model run.  It might then be
+necessary to increase the parameter \verb|maxpart| to a higher value than what
+would otherwise be needed.  The following is an example for an
+\verb|SPECIES_nnn| file including emission variation.
+
+\begin{scriptsize}\begin{verbatim}
+****************************************************************************
+*                                                                          *
+*     Input file for the Lagrangian particle dispersion model FLEXPART     *
+*               Definition file of chemical species/radionuclides          *
+*                                                                          *
+****************************************************************************
+EU-CO              Tracer name
+-999.9             Species half life
+-9.9E-09           Wet deposition - A
+                   Wet deposition - B
+-9.9               Dry deposition (gases) - D
+                   Dry deposition (gases) - Henrys const.
+                   Dry deposition (gases) - f0 (reactivity)
+-9.9E09            Dry deposition (particles) - rho
+                   Dry deposition (particles) - dquer
+                   Dry deposition (particles) - dsig
+-9.99              Alternative: dry deposition velocity
+ 28.00             molweight
+-9.9E-09           OH Reaction rate at 25 deg, [cm^3/s]
+-9                 number of associated specias (neg. none) - not implemented yet
+-99.99             KOA - organic matter air partitioning
+hr_start co_area  co_point
+ 0       0.535       0.932      0-1 local time
+ 1       0.405       0.931      1-2 local time
+ 2       0.317       0.927
+ 3       0.265       0.926
+ 4       0.259       0.928
+ 5       0.367       0.936
+ 6       0.668       0.952
+ 7       1.039       0.975
+ 8       1.015       1.046
+ 9       0.965       1.055
+10       1.016       1.061
+11       1.133       1.064
+12       1.269       1.067
+13       1.368       1.068
+14       1.516       1.069
+15       1.681       1.068
+16       1.777       1.024
+17       1.827       1.017
+18       1.538       1.008
+19       1.282       1.007
+20       1.136       1.004
+21       1.020       0.996
+22       0.879       0.981
+23       0.723       0.958      23-24 local time
+week_day co_area  co_point
+1        1.060       1.000      Monday
+2        1.060       1.000      Tuesday
+3        1.060       1.000      Wednesday
+4        1.060       1.000      Thursday
+5        1.060       1.000      Friday
+6        0.900       1.000      Saturday
+7        0.800       1.000      Sunday
+\end{verbatim}\end{scriptsize}
+
+\begin{footnotesize}\verb|IGBP_int1.dat|\end{footnotesize} contains the landuse inventory in binary format, and
+\begin{footnotesize}\verb|surfdata.t|\end{footnotesize}, shown below, gives the roughness lengths for each landuse class:
+
+\begin{footnotesize}
+\begin{verbatim}
+13 landuse categories are related roughness length
+--------------------------------------------------------
+landuse   comment                               z0
+--------------------------------------------------------
+ 1	Urban land                                   0.7	
+ 2	Agricultural land                            0.1
+ 3	Range land                                   0.1
+ 4	Deciduous forest                             1.
+ 5	Coniferous forest                            1.
+ 6	Mixed forest including wetland               0.7	
+ 7	Water, both salt and fresh                   0.001
+ 8	Barren land mostly desert                    0.01
+ 9	Nonforested wetland                          0.1
+10	Mixed agricultural and range land            0.1
+11	Rocky open areas with low growing shrubs     0.05
+12 Snow and ice                                 0.001
+13 Rainforest                                   1.
+\end{verbatim}
+\end{footnotesize}
+
+\newpage
+
+\subsubsection{File surfdepo.t}
+\verb|surfdepo.t| gives the resistances needed for the parameterization of dry deposition of gases for the 13 landuse classes and five seasonal categories.
+This file must not be changed by the user.\par
+ 
+\begin{scriptsize}\begin{verbatim}
+DRY DEPOSITION
+==============================================================================
+AFTER WESELY, 1989
+==============================================================================
+1 to 11: Landuse types after Wesely; 12 .. snow, 13 .. rainforest
+==============================================================================
+Values are tabulated for 5 seasonal categories:
+1     Midsummer with lush vegetation
+2     Autumn with unharvested cropland
+3     Late autumn after frost, no snow
+4     Winter, snow on ground and subfreezing
+5     Transitional spring with partially green short annuals
+==============================================================================
+               1       2       3       4       5       6       7       8       9     10       11      12      13  
+________________________________________________________________________________________________________________
+ri         9999.     60.    120.     70.    130.    100.   9999.   9999.     80.    100.    150.    9999.    200.  1
+rlu        9999.   2000.   2000.   2000.   2000.   2000.   9999.   9999.   2500.   2000.   4000.    9999.   1000.
+rac         100.    200.    100.   2000.   2000.   2000.      0.      0.    300.    150.    200.       0.   2000.
+rgss        400.    150.    350.    500.    500.    100.      0.   1000.      0.    220.    400.     100.    200.
+rgso        300.    150.    200.    200.    200.    300.   2000.    400.   1000.    180.    200.   10000.    200.
+rcls       9999.   2000.   2000.   2000.   2000.   2000.   9999.   9999.   2500.   2000.   4000.    9999.   9999.
+rclo       9999.   1000.   1000.   1000.   1000.   1000.   9999.   9999.   1000.   1000.   1000.    9999.   9999.
+_________________________________________________________________________________________________________________
+ri         9999.   9999.   9999.   9999.    250.    500.   9999.   9999.   9999.   9999.   9999.    9999.    200. 2
+rlu        9999.   9000.   9000.   9000.   4000.   8000.   9999.   9999.   9000.   9000.   9000.    9999.   1000.
+rac         100.    150.    100.   1500.   2000.   1700.      0.      0.    200.    120.    140.       0.   2000.
+rgss        400.    200.    350.    500.    500.    100.      0.   1000.      0.    300.    400.     100.    200.
+rgso        300.    150.    200.    200.    200.    300.   2000.    400.    800.    180.    200.   10000.    200.
+rcls       9999.   9000.   9000.   9000.   2000.   4000.   9999.   9999.   9000.   9000.   9000.    9999.   9999.
+rclo       9999.    400.    400.    400.   1000.    600.   9999.   9999.    400.    400.    400.    9999.   9999.
+_________________________________________________________________________________________________________________
+ri         9999.   9999.   9999.   9999.    250.    500.   9999.   9999.   9999.   9999.   9999.    9999.    200. 3
+rlu        9999.   9999.   9000.   9000.   4000.   8000.   9999.   9999.   9000.   9000.   9000.    9999.   1000.
+rac         100.     10.    100.   1000.   2000.   1500.      0.      0.    100.     50.    120.       0.   2000.
+rgss        400.    150.    350.    500.    500.    200.      0.   1000.      0.    200.    400.     100.    200.
+rgso        300.    150.    200.    200.    200.    300.   2000.    400.   1000.    180.    200.   10000.    200.
+rcls       9999.   9999.   9000.   9000.   3000.   6000.   9999.   9999.   9000.   9000.   9000.    9999.   9999.
+rclo       9999.   1000.    400.    400.   1000.    600.   9999.   9999.    800.    600.    600.    9999.   9999.
+_________________________________________________________________________________________________________________
+ri         9999.   9999.   9999.   9999.    400.    800.   9999.   9999.   9999.   9999.   9999.    9999.    200. 4
+rlu        9999.   9999.   9999.   9999.   6000.   9000.   9999.   9999.   9000.   9000.   9000.    9999.   1000.
+rac         100.     10.     10.   1000.   2000.   1500.      0.      0.     50.     10.     50.       0.   2000.
+rgss        100.    100.    100.    100.    100.    100.      0.   1000.    100.    100.     50.     100.    200.
+rgso        600.   3500.   3500.   3500.   3500.   3500.   2000.    400.   3500.   3500.   3500.   10000.    200. 
+rcls       9999.   9999.   9999.   9000.    200.    400.   9999.   9999.   9000.   9999.   9000.    9999.   9999.
+rclo       9999.   1000.   1000.    400.   1500.    600.   9999.   9999.    800.   1000.    800.    9999.   9999.
+_________________________________________________________________________________________________________________
+ri         9999.    120.    240.    140.    250.    190.   9999.   9999.    160.    200.    300.    9999.    200.  5
+rlu        9999.   4000.   4000.   4000.   2000.   3000.   9999.   9999.   4000.   4000.   8000.    9999.   1000.
+rac         100.     50.     80.   1200.   2000.   1500.      0.      0.    200.     60.    120.       0.   2000.
+rgss        500.    150.    350.    500     500.    200.      0.   1000.      0.    250.    400.     100.    200.
+rgso        300.    150.    200.    200.    200.    300.   2000.    400.   1000.    180.    200.   10000.    200.
+rcls       9999.   4000.   4000.   4000.   2000.   3000.   9999.   9999.   4000.   4000.   8000.    9999.   9999.
+rclo       9999.   1000.    500.    500.   1500.    700.   9999.   9999.    600.    800.    800.    9999.   9999.
+_________________________________________________________________________________________________________________
+\end{verbatim}\end{scriptsize}
+
+\end{document}
