source: flexpart.git/documentation/flexpart9.2.tex @ 6b2046d

10.4.1_peseiFPv9.3.1FPv9.3.1b_testingFPv9.3.2GFS_025bugfixes+enhancementsdevfp9.3.1-20161214-nc4grib2nc4_repairrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 6b2046d was 6b2046d, checked in by flexpart <>, 9 years ago

update documentation: flexpart9.2.tex

  • Property mode set to 100644
File size: 136.2 KB
Line 
1%% Version of source file:
2%% Date: 2003 July
3%%  **************************************
4%   *  TEMPLATE FOR EGU STYLE        *
5%   * FOR ARTICLES IN JOURNALS AND BOOKS *
6%%  **************************************
7%  Various alternatives for the input are shown, commented out.
8%  E.g., for the documentstyle options
9%        for the bibliography
10%  Feel free to play around with these variations, especially with
11%    the style options ms and twocolumn and 11pt/12pt
12
13 %%%% LATEX2E: SELECT ONE OF THE NEXT LINES
14\documentclass{egu}            % MANUSCRIPT, 10PT TYPE SIZE
15%\documentclass[ms,11pt]{egu}       % MANUSCRIPT, 11PT TYPE SIZE
16%\documentclass{egu}                 % for EGU TWO-COLUMN REVISED COPY
17%% <<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
18%%   REMOVE THIS LINE ONLY IF IT CAUSES PROBLEMS
19  \usepackage{times}                   % WITH TIMES ROMAN FONT
20% ^^^^^^^^^^^^^^^^^^
21%%   This is STRONGLY recommended for the revised version!!
22%% <<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
23
24
25%% ADD THIS PACKAGE IF GRAPHICS ARE TO BE IMPORTED
26\usepackage{graphicx} 
27\usepackage{amsmath} 
28%\usepackage{psfig}
29
30\printfigures              % PRINTS OUT FIGURES AT END OF MANUSCRIPT
31\newcommand{\degreee}{{$\, \rm ^{\circ}$}}
32\newcommand{\degreen}{{$\, \rm ^{\circ}$~}}
33
34\begin{document} 
35
36\title{The Lagrangian particle dispersion model FLEXPART version 9.2}
37
38\author[1]{A. Stohl}
39\author[1]{H. Sodemann}
40\author[1]{S. Eckhardt}
41\author[2]{A. Frank}
42\author[2]{P. Seibert}
43\author[3]{G. Wotawa}
44
45\affil[1]{Norwegian Institute of Air Research, Kjeller, Norway}
46\affil[2]{Institute of Meteorology, University of Natural Resources and Applied Life Sciences, Vienna, Austria}
47\affil[3]{Preparatory Commission for the Comprehensive Nuclear Test Ban Treaty Organization, Vienna, Austria}
48
49\runningtitle{FLEXPART description}
50\runningauthor{Stohl et al.}
51%\runninghead{} % This may also be used instead of \runningtitle and \runningauthor
52
53\correspondence{A. Stohl (ast@nilu.no)}
54
55
56%\journal{ACP}       % for NONLINEAR PROCESSES IN GEOPHYSICS
57
58\date{Manuscript version from 29 November 2010}
59
60% ADDITIONAL NON-STANDARD COMMANDS FOR TITLE BLOCK INFORMATON
61\firstauthor{Stohl}
62%\proofs{A. Stohl\\
63%Norwegian Institute for Air Research\\
64%Instituttveien 18, 2027 Kjeller, Norway}
65%\offsets{A. Stohl\\
66%Norwegian Institute for Air Research\\
67%Instituttveien 18, 2027 Kjeller, Norway}
68
69% The manuscript number is supplied by the editorial office
70\msnumber{12345}
71
72% The following commands will be activated by the EGU editorial/production office
73% after inserting appropriate values.
74
75
76\maketitle % YOU MUST USE THE \maketitle COMMAND
77
78\begin{abstract}
79The Lagrangian particle dispersion model FLEXPART was originally (in its first
80release in 1998) designed for calculating the long-range and mesoscale
81dispersion of air pollutants from point sources, such as after an accident in a
82nuclear power plant.  In the meantime FLEXPART has evolved into a comprehensive
83tool for atmospheric transport modeling and analysis.  Its application fields
84were extended from air pollution studies to other topics where atmospheric
85transport plays a role (e.g., exchange between the stratosphere and
86troposphere, or the global water cycle).  It has evolved into a true community
87model that is now being used by at least 35 groups from 14 different countries
88and is seeing both operational and research applications.  The last citable
89manuscript for FLEXPART is: \citep{stohl2005}
90\end{abstract}
91
92\section{Updates since FLEXPART version 8.0} 
93
94In this version 8.2 of FLEXPART the representation of physical processes was
95improved as well as a number of technical changes and bugfixes implemented.  In
96addition, the program is now released under a GNU GPL license.  For the first
97time, detailed installation instructions are provided to help users getting
98started with running FLEXPART.  A short section on the new Python routines
99pflexpart for reading FLEXPART output data has been included as well.
100\bigskip
101
102\noindent {\bf Technical Changes:}
103
104\begin{itemize} 
105
106\item grib2 compatibility for ECMWF data which will be introduced soon (new
107data retrieval routines available)
108
109\item AVAILABLE files can now contain up to 256 characters, and include the
110path directly with the input file name.  This is used to gather data files
111stored in different directories.  The third line in the \verb|pathnames| file
112should be left empty if long data file names are used.
113
114\item The code was updated to work with the ECMWF grib\_api V1.6.1 and above
115(tested up to 1.9.5)
116
117\item An important bug the concentration output routines was fixed.  A
118numerical problem lead in some circumstances to garbled data.  The header
119version identifier has been changed to version 8.2.  New routines for loading
120data available for download on the FLEXPART homepage.
121
122\item Several bugs in the wet deposition parameterization were fixed.
123
124\item A bug which lead to ignoring landuses that was introduced in version 8.0
125has been fixed.
126
127\item Several bugs concerning nested grids have been fixed.
128
129\end{itemize}
130
131\noindent {\bf Algorithm Changes:}
132
133\begin{itemize}
134
135\item A new, more detailed settling parameterisation for aerosols was
136implemented.  The dynamic viscosity of air is now calculated as a function of
137temperature, and the calculation of settling velocities is now more accurate
138for higher Reynolds numbers.
139
140\item It is now possible to produce output files for backward runs that can be used
141to interface FLEXPART with model output from another model or from a FLEXPART
142forward run.  A gridded output file containing the exact sensitivities to these
143initial conditions can produced.  To this end, a new switch has been added in
144the COMMAND file.
145
146\end{itemize}
147
148\section{Introduction}
149
150Lagrangian particle models compute trajectories of a large number of so-called
151particles (not necessarily representing real particles, but infinitesimally
152small air parcels) to describe the transport and diffusion of tracers in the
153atmosphere.  The main advantage of Lagrangian models is that, unlike in
154Eulerian models, there is no numerical diffusion.  Furthermore, in Eulerian
155models a tracer released from a point source is instantaneously mixed within a
156grid box, whereas Lagrangian models are independent of a computational grid and
157have, in principle, infinitesimally small resolution.
158
159The basis for current atmospheric particle models was laid by
160\citet{thomson1987}, who stated the criteria that must be fulfilled in order
161for a model to be theoretically correct.  A monograph on the theory of
162stochastic Lagrangian models was published by \citet{rodean1996} and another
163good review was written by \citet{wilson1996}.  The theory of modeling
164dispersion backward in time with Lagrangian particle models was developed by
165\citet{flesch1995} and \citet{seibert2004}.  Reviews of the more practical
166aspects of particle modeling were provided by \citet{zannetti1992} and
167\citet{uliasz1994}.
168
169This note describes FLEXPART, a Lagrangian particle dispersion model that
170simulates the long-range and mesoscale transport, diffusion, dry and wet
171deposition, and radioactive decay of tracers released from point, line, area or
172volume sources.  It can also be used in a domain-filling mode where the entire
173atmosphere is represented by particles of equal mass.  FLEXPART can be used
174forward in time to simulate the dispersion of tracers from their sources, or
175backward in time to determine potential source contributions for given
176receptors.  The management of input data was largely taken from FLEXTRA, a
177kinematic trajectory model \citep{stohl1995}.  FLEXPART's first version was
178developed during the first author's military service at the
179nuclear-biological-chemical school of the Austrian Forces, the deposition code
180was added soon later (version 2), and this version was validated using data
181from three large tracer experiments \citep{stohletal1998}.  Version 3 saw
182performance optimizations and the development of a density correction
183\citep{stohlthomson1999}.  Further updates included the addition of a
184convection scheme \citep{seibertetal2001} (version 4), better backward
185calculation capabilities \citep{seibert2004}, and improvements in the
186input/output handling (version 5).  Validation was done during intercontinental
187air pollution transport studies \citep{stohl1999, forster2001, spichtinger2001,
188stohl2002, stohl2003}, for which also special developments for FLEXPART were
189made in order to extend the forecasting capabilities \citep{stohl2004}.
190Version 6.2 saw corrections to the numerics in the convection scheme, the
191addition of a domain-filling option, and the possibility to use output nests.
192Version 6.4 runs with NCEP GFS model data.  Version 7.0 was a transition
193version that was not publicly released.  Version 8.0 was a major release.  It
194unifies the ECMWF and GFS model versions in one source package. GFS data in
195GRIB2 format can be read, using ECMWF's grib\_api library. Each species
196got its own definition file.  Output can be written individually for multiple
197species in backward runs.  The output format was changed to a compressed sparse
198matrix format.  Memory is partly allocated dynamically.  Furthermore, dry and
199wet deposition algorithms were updated and new, global landuse inventory
200introduced.  OH reaction based on a monthly averaged 3 dimensional OH-field is
201available as an option.
202
203FLEXPART is coded following the Fortran 95 standard and tested with several
204compilers (gfortran, Absoft, Portland Group) under a number of operating
205systems (Linux, Solaris, Mac~OS~X, etc.).  The code is carefully documented and
206optimized for run-time performance.  No attempts have been made to parallelize
207the code because the model is strictly linear and, therefore, it is most
208effective to partition problems such that they run on single processors and to
209combine the results if needed.
210
211FLEXPART's source code and a manual are freely available from the internet page
212http://transport.nilu.no/flexpart.  According to a recent user survey, at least
21334 groups from 17 countries are currently using FLEXPART.  The user community
214maintains discussion by a mailing list to which one can subscribe on the
215flexpart home page.  The version of FLEXPART described here is based on model
216level data of the numerical weather prediction model of the European Centre for
217Medium-Range Weather Forecasts (ECMWF).  The standard source code distribution
218contains also the source files for a version of FLEXPART using the global
219National Centers of Environmental Prediction (NCEP) model data on pressure
220levels.  Other users have developed FLEXPART versions using input data from a
221suite of different and meso-scale (e.g.  MM5, WRF, COSMO) models, some of which
222are available from the FLEXPART website but are not described here.
223
224\section{License}
225
226FLEXPART has been free software ever since it first was released. The status as
227a free software is now more formally established by releasing the code under
228the GNU General Public License (GPL) Version 3.  The text of the license is
229included in the file COPYING in the source code archive.
230
231\section{Installation}\label{sec:installation}
232
233Getting FLEXPART up and running on their systems has been a major hurdle to many
234new users of the software. Since the inclusion of the grib\_api library has further
235complicated the compilation of FLEXPART, we include a step-by-step installation
236instruction of FLEXPART.  The steps below are described for an Ubunto 10.4 LT
237UNIX release.
238
239\subsection{Installing required libraries}
240
241In order to process input files in GRIB2 format, FLEXPART V8.2 needs to be
242compiled with ECMWF's grib\_api library version 1.6.1 and above.  Since the API
243of the grib\_api library may change in the future, upward compatibility cannot
244be guaranteed.  The input files in GRIB2 format can be compressed to save
245bandwitdth and storage space.  In order to make use of this feature, the jasper
246library needs to be installed on the system.  Optionally, the emos library can
247be used to run FLEXPART with input files in GRIB1 format.
248
249\noindent The following steps should be executed in sequence:
250
251\subsubsection{Install jasper library (Version 1.900.1)}
252
253Download the jasper library from the jasper project page\footnote{http://www.ece.uvic.ca/~mdadams/jasper/} 
254
255\begin{small}
256\begin{verbatim}
257unzip jasper-1.900.1.zip
258cd jasper-1.900.1
259./configure [--prefix=<installation path>]
260make
261make check
262make install
263\end{verbatim}
264\end{small}
265
266Parameters in brackets are optional but may require root privileges. 
267
268\subsubsection{Install grib\_api library (Version 1.6.1 or later)}
269 
270Download the grib\_api library from the ECMWF
271website\footnote{http://www.ecmwf.int/products/data/software/download/grib\_api.html}
272
273\begin{small}
274\begin{verbatim}
275tar -xvf grib_api-1.12.3.tar.gz
276./configure [--with-jasper=<jasper path>]
277make
278make check
279make install
280\end{verbatim}
281\end{small}
282
283\subsubsection{Optional: Install emos library (Version 000372)}
284
285Download the emos library from the ECMWF
286website\footnote{http://www.ecmwf.int/products/data/software/interpolation.html}
287
288\begin{small}
289\begin{verbatim}
290tar -xvf emos_000372.tar.gz
291./build_library
292./install
293\end{verbatim}
294\end{small}
295
296\subsection{Compiling FLEXPART V8.2}
297
298Download the FLEXPART source code archive from the FLEXPART
299homepage\footnote{http://transport.nilu.no/flexpart/flexpart/view}
300
301\begin{small}
302\begin{verbatim}
303tar -xvf flexpart82.tar.gz
304\end{verbatim}
305\end{small}
306optionally edit the file includepar to set parameters for the data center, grid
307dimension, and particle number edit the LIBRARY path variable in the makefiles
308according to the position of libgrib\_api and libjasper
309\begin{small}
310\begin{verbatim}
311make -f makefile.<center>_<compiler>_<system>
312\end{verbatim}
313\end{small}
314In the above statement, center can be one of: gfs, ecmwf, ecmwf\_emos,
315gfs\_emos compiler can be one of absoft or gfortran (emos library only with
316absoft) system can be one of 32, 64 (emos only 32).  The system parameter must
317match that of the compiled libraries. See also Table~\ref{tab:makefile} for all
318available makefiles.
319
320When recompiling after making changes, all object and module files can be
321removed safely by using
322\begin{small}
323\begin{verbatim}
324make -f makefile.<xxx> clean
325\end{verbatim}
326\end{small}
327
328\begin{table*}
329\setlength{\tabcolsep}{1.1mm}
330\caption{\label{tab:makefile}
331List of available makefiles }
332\vspace{3mm}
333{\centerline{
334\begin{tabular}{llc} \hline
335Makefile GFS & Makefile ECMWF & GRIB version \\ 
336\hline
337makefile.gfs\_gfortran\_32 & makefile.ecmwf\_gfortran\_32 & 1/2 \\
338makefile.gfs\_gfortran\_64 & makefile.ecmwf\_gfortran\_64 & 1/2 \\
339makefile.gfs\_absoft\_32 & makefile.ecmwf\_absoft\_32 & 1/2 \\
340makefile.gfs\_absoft\_64 & makefile.ecmwf\_absoft\_64 & 1/2 \\
341makefile.gfs\_emos\_absoft\_32 & makefile.ecmwf\_emos\_absoft\_32 & 1 \\
342\hline
343\end{tabular}}}
344\end{table*}
345
346\section{Input data and grid definitions}
347
348FLEXPART is an off-line model that uses meteorological fields (analyses or
349forecasts) in Gridded Binary (GRIB) format in version 1 or 2 from the ECMWF
350numerical weather prediction model \citep{ecmwf1995} on a latitude/longitude
351grid and on native ECMWF model levels as input.  Optionally, GRIB data from
352NCEP's GFS model, available on pressure levels, can be used.  The ECMWF data
353can be retrieved from the ECMWF archives using a pre-processor that is also
354available from the FLEXPART website but not described here.  The GRIB decoding
355libraries is {\it not} provided with the FLEXPART source codes but is publicly
356available (see Sec.~\ref{sec:installation}.  The data can be global or only
357cover a limited area.  Furthermore, higher-resolution domains can be nested
358into a mother domain.
359
360The file \verb|includepar| contains all relevant FLEXPART parameter settings,
361both physical constants and maximum field dimensions.  As the memory required
362by FLEXPART is determined by the various field dimensions, it is recommended
363that they are adjusted to actual needs before compilation.  The file
364\verb|includecom| defines all FLEXPART global variables and fields, i.e., those
365shared between most subroutines.
366
367\subsection{Input data organisation}
368
369A file \verb|pathnames| must exist in the directory where FLEXPART is started.
370It must contain at least four lines:\\ 1.  line: Directory where all the
371FLEXPART command files are stored.\\ 2.  line: Directory to which the model
372output is written.\\ 3.  line: Directory where the GRIB input fields are
373located.\\ 4.  line: Path name of the \verb|AVAILABLE| file (see below).\\ If
374nests with higher-resolution input data shall also be used, lines 3 and 4 must
375be repeated for every nest, thus also specifying the nesting level order.  Any
376number of nesting levels can be used up to a maximum (parameter
377\verb|maxnests|).
378
379The meteorological input data must be organised such that all data for a domain
380and a certain date must be contained in a single GRIB file.  The
381\verb|AVAILABLE| file lists all available dates and the corresponding file
382names.  For each nesting level, the input files must be stored in a different
383directory and the \verb|AVAILABLE| file must contain the same dates as for the
384mother domain.  Given a certain particle position, the last (i.e., innermost)
385nest is checked first whether it contains the particle or not.  If the particle
386resides in this nest, the meteorological data from this nest is interpolated
387linearly to the particle position.  If not, the next nest is checked, and so
388forth until the mother domain is reached.  There is no nesting in the vertical
389direction and the poles must not be contained in any nest.
390
391The maximum dimensions of the meteorological fields are specified by the
392parameters \verb|nxmax, nymax, nuvzmax, nwzmax, nzmax| in file
393\verb|includepar|, for x, y, and three z dimensions, respectively.  The three z
394dimensions are for the original ECMWF data (\verb|nuvzmax, nwzmax| for model
395half levels and model levels, respectively) and transformed data (\verb|nzmax|,
396see below), respectively.  The horizontal dimensions of the nests must be
397smaller than the parameters \verb|nxmaxn, nymaxn|.  Grid dimensions and other
398basic things are checked in routine \verb|gridcheck.f| (\verb|gridcheck_gfs.f|
399for the GFS version), and error messages are issued if necessary.
400
401The longitude/latitude range of the mother grid is also used as the
402computational domain.  All internal FLEXPART coordinates run from the
403western/southern domain boundary with coordinates (0,0) to the eastern/northern
404boundary with coordinates (nx-1,ny-1), where (nx,ny) are the mother grid
405dimensions.  For global input data, FLEXPART repeats the westernmost grid cells
406at the easternmost domain "boundary", in order to facilitate interpolation on
407all locations of the globe (e.g., if input data run from 0 to 359\degreen with
4081\degreen resolution, 0\degreen data are repeated at 360\degreee).  A global
409mother domain can be shifted by \verb|nxshift| (file \verb|includepar|) data
410columns (subroutines \verb|shift_field.f| and \verb|shift_field_0.f|) if nested
411input fields would otherwise overlap the "boundaries".  For instance, a domain
412stretching from 320\degreen to 30\degreen can be nested into the mother grid of
413the above example by shifting the mother grid by 30\degreee. The default setting
414for global ECMWF fields is \verb|nxshift=359|, while GFS fields the default
415value is \verb|nxshift=0|.
416
417\subsection{Vertical model structure and required data}
418
419FLEXPART needs five three-dimensional fields: horizontal and vertical wind
420components, temperature and specific humidity.  Input data must be on ECMWF
421model (i.e.  $\eta$) levels which are defined by a hybrid coordinate system.
422The conversion from $\eta$ to pressure coordinates is given by $p_k=A_k+B_kp_s$
423and the heights of the $\eta$ surfaces are defined by $\eta_k=A_k/p_0+B_k$,
424where $\eta_k$ is the value of $\eta$ at the $k^{\rm_{th}}$ model level, $p_s$
425is the surface pressure and $p_0=$101325 Pa.  $A_k$ and $B_k$ are coefficients,
426chosen such that the levels closest to the ground follow the topography, while
427the highest levels coincide with pressure surfaces; intermediate levels
428transition between the two.  The vertical wind in hybrid coordinates is
429calculated mass-consistently from spectral data by the pre-processor.  A
430surface level is defined in addition to the regular $\eta$ levels.  2~m
431temperature, 10~m winds and specific humidity from the first regular model
432level are assigned to this level, to represent "surface" values.
433
434Parameterized random velocities in the atmospheric boundary layer (ABL, see
435sections~\ref{PBLparameterization} and \ref{diffusion}) are calculated in units
436of m s$^{-1}$, and not in $\eta$ coordinates.  Therefore, in order to avoid
437time-consuming coordinate transformations every time step, all
438three-dimensional data are interpolated linearly from the ECMWF model levels to
439terrain-following Cartesian coordinates $\tilde{z}=z-z_t$, where $z_t$ is the
440height of the topography (subroutine \verb|verttransform.f|).  The conversion
441of vertical wind speeds from the eta coordinate system into the
442terrain-following co-ordinate system follows as
443
444\begin{equation}
445\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}  \,
446\end{equation}
447
448where $\dot{\tilde{\eta}}=\dot{\eta} \partial p / \partial \eta$.  The second
449term on the right hand side is missing in the FLEXPART transformation because
450it is much smaller than the others.  One colleague has implemented this term in
451his version of FLEXPART and found virtually no differences (B.  Legras,
452personal communication).
453
454FLEXPART also needs the two-dimensional fields: surface pressure, total cloud
455cover, 10~m horizontal wind components, 2~m temperature and dew point
456temperature, large scale and convective precipitation, sensible heat flux,
457east/west and north/south surface stress, topography, land-sea-mask and subgrid
458standard deviation of topography.  The landuse inventory of \citet{belward1999}
459is provided in an extra file in the options directory (\verb|IGBP_int1.dat|).
460
461Note that GFS data is provided on pressure levels (currently 26) which requires
462that the FLEXPART code is compiled with a makefile for the GFS model. GFS data
463are freely available from the NCEP data archives. An example for a retrieval
464of GFS data files via ftp is given below:
465
466\begin{small}
467\begin{verbatim}
468cd $DIR_GFS_05TEMPORARY
469ftp ftpprd.ncep.noaa.gov << EOF11
470binary
471cd /pub/data/nccf/com/gfs/prod/gfs.2009111006
472get gfs.t06z.pgrb2f00 GF09111006
473quit
474EOF11
475\end{verbatim}
476\end{small}
477
478\subsection{Meteorological input data formats}
479
480While NCEP has already switched to the more flexible and compact GRIB2 format,
481the ECMWF is only gradually transitioning from GRIB1 to the new GRIB2 format.
482Transition to GRIB2 becomes necessary, at least for 3D fields, when the ECMWF
483increases the vertical resolution beyond 100 model levels in 2012 (which cannot
484be handled by the GRIB1 format).  Currently, GRIB2 codes have been defined for
485all 3D-fields required by FLEXPART, while definitions are still missing for a
486part of the 2D-fields.
487
488According to the ECMWF it may take another year to define the GRIB2 codes for
489all 2D-fields used by FLEXPART.  Therefore, in version 8.2, it is now possible
490to read in GRIB files that contain only GRIB1 coded fields, files that contain
491only GRIB2 coded fields, and files that contain mixed GRIB1/2 fields.  Such
492mixed GRIB1/2 files are produced by the current version 4 of the FLEXPART data
493retrieval routines (available from the FLEXPART hompage).
494
495\section[Physical parameterizations]{Physical parameterization of boundary layer parameters}
496
497\label{PBLparameterization} Accumulated surface sensible heat fluxes and
498surface stresses are available from ECMWF forecasts.  The pre-processor selects
499the shortest forecasts available for that date from the ECMWF archives and
500deaccumulates the flux data.  The total surface stress is computed from
501
502\begin{equation}
503\tau=\sqrt{\tau_1^2+\tau_2^2} \;,
504\end{equation}
505
506where $\tau_1$ and $\tau_2$ are the surface stresses in east/west and
507north/south direction, respectively.  Friction velocity is then calculated in
508subroutine \verb|scalev.f| as
509
510\begin{equation}
511u_*=\sqrt{\tau/\rho} \;,
512\end{equation}
513
514where $\rho$ is the air density \citep{wotawa1996}.  Friction velocities and
515heat fluxes calculated using this method are most accurate \citep{wotawa1997}.
516However, if deaccumulated surface stresses and surface sensible heat fluxes are
517not available, the profile method after \citet{berkowicz1982} (subroutine
518\verb|pbl_profile.f|) is applied to wind and temperature data at the second
519model level and at 10~m (for wind) and 2~m (for temperature) (note that
520previously the first model level was used; as ECMWF has its first model level
521now close to 10~m, the second level is used instead).  The following three
522equations are solved iteratively:
523
524\begin{equation}
525u_*=\frac{\kappa\Delta u}{\ln {\frac{z_l}{10}} -\Psi_m(\frac{z_l}{L}) +\Psi_m(\frac{10}{L})} \; ,
526\end{equation}
527
528\begin{equation}
529\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]} \; ,
530\end{equation}
531
532\begin{equation}
533L= \frac{\overline{T}u_*^2}{g \kappa \Theta_*} \;,
534\end{equation}
535
536where $\kappa$ is the von K\'{a}rm\'{a}n constant (0.4), $z_l$ is the height of
537the second model level, $\Delta u$ is the difference between wind speed at the
538second model level and at 10~m, $\Delta \Theta$ is the difference between
539potential temperature at the second model level and at 2~m, $\Psi_m$ and
540$\Psi_h$ are the stability correction functions for momentum and heat
541\citep{businger1971, beljaars1991}, $g$ is the acceleration due to gravity,
542$\Theta_*$ is the temperature scale and $\overline{T}$ is the average surface
543layer temperature (taken as $T$ at the first model level).  The heat flux is
544then computed by
545
546\begin{equation}
547(\overline{w'\Theta_v'})_0=-\rho c_p u_* \Theta_* \;,
548\end{equation}
549
550where $\rho c_p$ is the specific heat capacity of air at constant pressure.
551
552ABL heights are calculated according to \citet{vogelezang1996} using the
553critical Richardson number concept (subroutine \verb|richardson.f|).  The ABL
554height $h_{mix}$ is set to the height of the first model level $l$ for which
555the Richardson number
556
557\begin{equation}
558Ri_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}\;,
559\label{richardson}
560\end{equation}
561
562exceeds the critical value of 0.25$\Theta_{v1}$ and $\Theta_{vl}$ are the
563virtual potential temperatures, $z_1$ and $z_l$ are the heights of, and
564$(u_1,v_1)$, and $(u_l,v_l)$ are the wind components at the 1$^{\rm st}$ and
565$l^{\rm th}$ model level, respectively.  The formulation of
566Eq.~\ref{richardson} can be improved for convective situations by replacing
567$\Theta_{v1}$ with
568
569\begin{equation}
570\Theta_{v1}'=\Theta_{v1}+8.5\frac{(\overline{w'\Theta_v'})_0}{w_* c_p}\;,
571\label{excess}
572\end{equation}
573
574where
575
576\begin{equation}
577w_*=\left[{ \frac{(\overline{w'\Theta_v'})_0 g h_{mix}}{\Theta_{v1} c_p} }\right]^{1/3}
578\end{equation}
579
580is the convective velocity scale.  The second term on the right hand side of
581Eq.~\ref{excess} represents a temperature excess of rising thermals.  As $w_*$
582is unknown beforehand, $h_{mix}$ and $w_*$ are calculated iteratively.
583
584Spatial and temporal variations of ABL heights on scales not resolved by the
585ECMWF model play an important role in determining the thickness of the layer
586over which tracer is effectively mixed.  The height of the convective ABL
587reaches its maximum value (say 1500~m) in the afternoon (say, at 1700 local
588time (LT)), before a much shallower stable ABL forms.  Now, if meteorological
589data are available only at 1200 and 1800 LT and the ABL heights at those times
590are, say, 1200~m and 200~m, and linear interpolation is used, the ABL height at
5911700 LT is significantly underestimated (370~m instead of 1500~m).  If tracer
592is released at the surface shortly before the breakdown of the convective ABL,
593this would lead to a serious overestimation of the surface concentrations (a
594factor of four in the above example).  Similar arguments hold for spatial
595variations of ABL heights due to complex topography and variability in landuse
596or soil wetness \citep{hubbe1997}.  The thickness of a tracer cloud traveling
597over such a patchy surface would be determined by the maximum rather than by
598the average ABL height.
599
600In FLEXPART a somewhat arbitrary parameterization is used to avoid a
601significant bias in the tracer cloud thickness and the surface tracer
602concentrations.  To account for spatial variations induced by topography, we
603use an "envelope" ABL height
604
605\begin{equation}
606H_{env}=h_{mix}+\min \left[\sigma_Z, c \frac{V}{N} \right]\; .
607\end{equation}
608
609Here, $\sigma_Z$ is the standard deviation of the ECMWF model subgrid
610topography, $c$ is a constant (here: 2.0), $V$ is the wind speed at height
611$h_{mix}$, and $N$ is the Brunt-Vaisala frequency.  Under convective
612conditions, the envelope ABL height is, thus, the diagnosed ABL height plus the
613subgrid topography (assuming that the ABL height over the hill tops effectively
614determines the dilution of a tracer cloud located in a convective ABL).  Under
615stable conditions, air tends to flow around topographic obstacles rather than
616above it, but some lifting is possible due to the available kinetic energy.
617$\frac{V}{N}$ is the local Froude number (i.e., the ratio of inertial to
618buoyant forces) times the length scale of the sub-grid topographic obstacle.
619The factor $c \frac{V}{N}$, thus, limits the effect of the subgrid topography
620under stable conditions, with $c=2$ being a subjective scaling factor.
621$H_{env}$ rather than $h_{mix}$ is used for all subsequent calculations.  In
622addition, $H_{env}$ is not interpolated to the particle position, but the
623maximum $H_{env}$ of the grid points surrounding a particle's position in space
624and time is used.
625
626\section{\label{diffusion}Particle transport and diffusion}
627
628\subsection{Particle trajectory calculations}
629
630FLEXPART generally uses the simple ``zero acceleration'' scheme
631
632\begin{equation}
633{\vec{X}}(t+\Delta t)={\vec{X}}(t)+{\vec{v}}({\vec{X}},t) \Delta t\,,
634\label{firstorder}
635\end{equation}
636
637which is accurate to the first order, to integrate the trajectory equation \citep{stohl1998}
638
639\begin{equation}
640\frac{d {\vec{X}}}{dt}={\vec{v}}[{\vec{X}}(t)] \,,
641\end{equation}
642
643with $t$ being time, $\Delta t$ the time increment, ${\vec{X}}$ the position
644vector, and ${\vec{v}}=\overline{\vec{v}}+{\vec{v}}_t+{\vec{v}}_m$ the wind
645vector that is composed of the grid scale wind $\overline{\vec{v}}$, the
646turbulent wind fluctuations ${\vec{v}}_t$ and the mesoscale wind fluctuations
647${\vec{v}}_m$.
648
649Since FLEXPART version 5.0, numerical accuracy has been improved by making one
650iteration of the \citet{petterssen1940} scheme (which is accurate to the second
651order) whenever this is possible, but only for the grid-scale winds.  It is
652implemented as a correction applied to the position obtained with the ``zero
653acceleration'' scheme.  In three cases it cannot be applied.  First, the
654Petterssen scheme needs winds at a second time which may be outside the time
655interval of the two wind fields kept in memory.  Second, if a particle crosses
656the boundaries of nested domains, and third in the ABL if \verb|ctl|$>0$ (see
657below).
658
659Particle transport and turbulent dispersion are handled by the subroutine
660\verb|advance.f| where calls are issued to procedures that interpolate winds
661and other data to the particle position and the Langevin equations (see below)
662are solved.  The poles are singularities on a latitude/longitude grid.  Thus,
663horizontal winds (variables \verb|uu,vv|) poleward of latitudes
664(\verb|switchnorth, switchsouth|) are transformed to a polar stereographic
665projection (variables \verb|uupol,vvpol|) on which particle advection is
666calculated.  As \verb|uupol,vvpol| are also stored on the latitude/longitude
667grid, no additional interpolation is made.
668
669\subsection{The Langevin equation}
670
671Turbulent motions ${\vec{v}}_{t}$ for wind components $i$ are parameterized
672assuming a Markov process based on the Langevin equation \citep{thomson1987}
673
674\begin{equation}
675dv_{t_i}=a_i({\vec{x}},{\vec{v}}_t,t)dt+b_{ij}({\vec{x}},{\vec{v}}_t,t)dW_
676j \,,
677\label{langevin}
678\end{equation}
679
680where the drift term $a$ and the diffusion term $b$ are functions of the
681position, the turbulent velocity and time.  $dW_j$ are incremental components
682of a Wiener process with mean zero and variance $dt$, which are uncorrelated in
683time \citep{legg1982}.  Cross-correlations between the different wind
684components are also not taken into account, since they have little effect for
685long-range dispersion \citep{uliasz1994}.
686
687Gaussian turbulence is assumed in FLEXPART, which is strictly valid only for
688stable and neutral conditions.  Under convective conditions, when turbulence is
689skewed and larger areas are occupied by downdrafts than by updrafts, this
690assumption is violated, but for transport distances where particles are rather
691well mixed throughout the ABL, the error is minor.
692
693With the above assumptions, the Langevin equation for the vertical wind
694component $w$ can be written as
695
696\begin{multline}
697d w = \\
698- w \frac{dt}{\tau_{L_w}}
699+ \frac{\partial \sigma_w^2}{\partial z} dt
700+ \frac{\sigma_w^2}{\rho} \frac{\partial \rho}{\partial z} dt
701+ \left( \frac{2}{\tau_{L_w}} \right)^{1/2} \sigma_w\; dW \;,
702\label{legg}
703\end{multline}
704
705where $w$ and $\sigma_w$ are the turbulent vertical wind component and its
706standard deviation, $\tau_{L_w}$ is the Lagrangian timescale for the vertical
707velocity autocorrelation and $\rho$ is density.  The second and the third term
708on the right hand side are the drift correction \citep{mcnider1988} and the
709density correction \citep{stohlthomson1999}, respectively.  This Langevin
710equation is identical to the one described by \citet{legg1982}, except for the
711term from \citet{stohlthomson1999} which accounts for the decrease of air
712density with height.
713
714Alternatively, the Langevin equation can be re-expressed in terms of
715$w/\sigma_w$ instead of $w$ \citep{wilson1983}:
716
717\begin{multline}
718d \left( \frac{w}{\sigma_w} \right) =\\
719- \frac{w}{\sigma_w} \frac{dt}{\tau_{L_w}}
720+ \frac{\partial \sigma_w}{\partial z} dt
721+ \frac{\sigma_w}{\rho} \frac{\partial \rho}{\partial z} dt
722+ \left( \frac{2}{\tau_{L_w}} \right)^{1/2} dW \;,
723\label{wilson}
724\end{multline}
725
726This form was shown by \citet{thomson1987} to fulfill the well-mixed criterion
727which states that ``if a species of passive marked particles is initially mixed
728uniformly in position and velocity space in a turbulent flow, it will stay that
729way'' \citep{rodean1996}.  Although the method proposed by \citet{legg1982}
730violates this criterion in strongly inhomogeneous turbulence, their formulation
731was found to be practical, as numerical experiments have shown that it is more
732robust against an increase in the integration time step.  Therefore,
733Eq.~\ref{legg} is used with long time steps (see section~\ref{timestep});
734otherwise, Eq.~\ref{wilson} is used.  For the horizontal wind components, the
735Langevin equation is identical to Eq.~\ref{legg}, with no drift and density
736correction terms.
737
738For the discrete time step implementation of the above Langevin equations (at
739the example of Eq.~\ref{wilson}), two different methods are used.  When $\left
740(\Delta t/\tau_{L_w} \right) \ge 0.5$,
741
742\begin{multline}
743\left (\frac{w}{\sigma_w} \right)_{k+1}=
744r_w \left (\frac{w}{\sigma_w} \right)_k
745+ \frac{\partial \sigma_w}{\partial z} \tau_{L_w} \left (1-r_w \right )\\
746+ \frac{\sigma_w}{\rho} \frac{\partial \rho}{\partial z} \tau_{L_w} \left (1-r_w \right )
747+ \left (1-r_w^2 \right )^{1/2} \zeta \;,
748\end{multline}
749
750where $r_w=\exp(-\Delta t/ \tau_{L_w})$ is the autocorrelation of the vertical
751wind, and $\zeta$ is a normally distributed random number with mean zero and
752unit standard deviation.  The subscripts $k$ and $k+1$ refer to subsequent
753times separated by $\Delta t$.
754
755To save computation time for cases when $\left (\Delta t/\tau_{L_w} \right) <
7560.5$, the following first order approximation is used in order to avoid the
757computation of the exponential function:
758
759\begin{multline}
760\left (\frac{w}{\sigma_w} \right)_{k+1}=
761\left (1-\frac{\Delta t}{\tau_{L_w}} \right )\left (\frac{w}{\sigma_w} \right)_k\\
762+ \frac{\partial \sigma_w}{\partial z} \Delta t
763+ \frac{\sigma_w}{\rho} \frac{\partial \rho}{\partial z} \Delta t
764+ \left( \frac{2 \Delta t}{\tau_{L_w}} \right)^{1/2} \zeta \;.
765\end{multline}
766
767When a particle reaches the surface or the top of the ABL, it is reflected and
768the sign of the turbulent velocity is changed \citep{wilson1993}.
769
770\subsection{\label{timestep}Determination of the time step}
771
772FLEXPART can be used in two different modes.  The computationally faster one
773(\verb|ctl|$<$0 in file \verb|COMMAND|) does not adapt the computation time
774step to the Lagrangian timescales $\tau_{L_i}$ (where $i$ is one of the three
775wind components) and FLEXPART uses constant time steps of one synchronisation
776time interval (\verb|lsynctime|, specified in file \verb|COMMAND|, typically
777900 seconds).  Usually, autocorrelations are very low in this mode and
778turbulence is not described well.  Nevertheless, for large scale applications
779FLEXPART works very well with this option \citep{stohletal1998}.  If turbulence
780shall be described more accurately, the time steps must be limited by $\tau_L$.
781Since the vertical wind is most important, only $\tau_{L_w}$ is used for this.
782The user must specify two constants, \verb|ctl| and \verb|ifine| in file
783\verb|COMMAND|.  The first one determines the time step $\Delta t_i$ according
784to
785
786\begin{equation}
787\Delta t_i=\frac{1}{c_{tl}} \min \left ({\tau_{L_w}\, , \; \frac{h}{2w}\, , \; \frac{0.5}{\partial \sigma_w / \partial z}} \right ) \;.
788\end{equation}
789
790The minimum value of $\Delta t_i$ is 1~second.
791$\Delta t_i$ is used for solving the Langevin equations for the horizontal turbulent wind components.
792
793For solving the Langevin equation for the vertical wind component, a shorter
794time step $\Delta t_w=\Delta t_i / \verb|ifine|$ is used.  However, note that
795there is no interaction between horizontal and vertical wind components on
796timescales less than $\Delta t_i$.  This strategy (given sufficiently large
797values for \verb|ctl| and \verb|ifine|) ensures that the particles stay
798vertically well-mixed also in very inhomogeneous turbulence, while keeping the
799computational cost at a minimum.
800
801\subsection{Parameterization of the wind fluctuations}
802
803For $\sigma_{v_i}$ and $\tau_{L_i}$ \citet{hanna1982} proposed a
804parameterization scheme based on the boundary layer parameters $h$, $L$, $w_*$,
805$z_0$ and $u_*$, i.e.  ABL height, Monin-Obukhov length, convective velocity
806scale, roughness length and friction velocity, respectively.  It is used in
807subroutines \verb|hanna.f, hanna1.f, hanna_short.f| with a modification taken
808from \citet{ryall1997} for $\sigma_w$, as Hanna's scheme does not always yield
809smooth profiles of $\sigma_w$ throughout the whole convective ABL.  In the
810following, subscripts $u$ and $v$ refer to the along-wind and the cross-wind
811components (transformed to grid coordinates in subroutine \verb|windalign.f|),
812respectively, and $w$ to the vertical component of the turbulent velocities;
813$f$ is the Coriolis parameter.  The minimum $\tau_{L_u}$, $\tau_{L_v}$ and
814$\tau_{L_w}$ used are 10~s, 10~s and 30~s, respectively, in order to avoid
815excessive computation times for particles close to the surface.\\[0.3cm]
816
817{\bf Unstable conditions:}
818
819\begin{equation}
820\frac{\sigma_u}{u_*}=\frac{\sigma_v}{u_*}=\left(12+\frac{h}{2|L|}\right)^{1/3}
821\end{equation}
822\begin{equation}
823\tau_{L_u}=\tau_{L_v}=0.15\frac{h}{\sigma_u}
824\end{equation}
825\begin{multline}
826\sigma_w=\\
827\left[ 1.2 w_*^2 \left (1-0.9\frac{z}{h} \right )
828\left ( \frac{z}{h} \right )^{2/3}
829+\left (1.8-1.4 \frac{z}{h} \right ) u_*^2 \right]^{1/2}
830\end{multline}
831For $z/h<0.1$ and $z-z_0>-L$:
832\begin{equation}
833\tau_{L_w}=0.1\frac{z}{\sigma_w\left[0.55-0.38\left(z-z_0\right)/L\right]}
834\end{equation}
835For $z/h<0.1$ and $z-z_0<-L$:
836\begin{equation}
837\tau_{L_w}=0.59\frac{z}{\sigma_w}
838\end{equation}
839For $z/h>0.1$:
840\begin{equation}
841\tau_{L_w}=0.15\frac{h}{\sigma_w}\left[1-\exp\left(\frac{-5z}{h}\right)\right]
842\end{equation}
843
844{\bf Neutral conditions:}
845
846\begin{equation}
847\frac{\sigma_u}{u_*}=2.0\exp(-3fz/u_*)
848\end{equation}
849\begin{equation}
850\frac{\sigma_v}{u_*}=\frac{\sigma_w}{u_*}=1.3\exp(-2fz/u_*)
851\end{equation}
852\begin{equation}
853\tau_{L_u}=\tau_{L_v}=\tau_{L_w}=\frac{0.5z/\sigma_w}{1+15fz/u_*}
854\end{equation}
855
856{\bf Stable conditions:}
857
858\begin{equation}
859\frac{\sigma_u}{u_*}=2.0\left(1-\frac{z}{h}\right)
860\end{equation}
861\begin{equation}
862\frac{\sigma_v}{u_*}=\frac{\sigma_w}{u_*}=1.3\left(1-\frac{z}{h}\right)
863\end{equation}
864\begin{equation}
865\tau_{L_u}=0.15\frac{h}{\sigma_u}\left(\frac{z}{h}\right)^{0.5}
866\end{equation}
867\begin{equation}
868\tau_{L_v}=0.07\frac{h}{\sigma_v}\left(\frac{z}{h}\right)^{0.5}
869\end{equation}
870\begin{equation}
871\tau_{L_w}=0.1\frac{h}{\sigma_w}\left(\frac{z}{h}\right)^{0.5}
872\end{equation}
873
874Lacking suitable turbulence parameterizations above the ABL ($z>h$), a constant
875vertical diffusivity $D_z$=0.1~m$^2$s$^{-1}$ is used in the stratosphere,
876following recent work of \citet{legras2003}, whereas a horizontal diffusivity
877$D_h$=50~m$^2$s$^{-1}$ is used in the free troposphere.  Stratosphere and
878troposphere are distinguished based on a threshold of 2~pvu (potential
879vorticity units).  Diffusivities are converted into velocity scales using
880$\sigma_{v_i}=\sqrt{D_i/dt}$.
881
882\subsection{Mesoscale velocity fluctuations}
883
884Mesoscale motions are neither resolved by the ECMWF data nor covered by the
885turbulence parameterization.  This unresolved spectral interval needs to be
886taken into account at least in an approximate way, since mesoscale motions can
887significantly accelerate the growth of a dispersing plume \citep{gupta1997}.
888For this, we use a similar method as \citet{maryon1998}, namely to solve an
889independent Langevin equation for the mesoscale wind velocity fluctuations
890(``meandering'' in Maryon's terms).  Assuming that the variance of the wind at
891the grid scale provides some information on its subgrid variance, the wind
892velocity standard deviation used for the mesoscale Langevin equation is set to
893\verb|turbmesoscale| (set in file \verb|includepar|) times the standard
894deviation of the grid points surrounding the particle's position.  The
895corresponding time scale is taken as half the interval at which wind fields are
896available, assuming that the linear interpolation between the grid points can
897recover half the subgrid variability, not an unlikely assumption
898\citep{stohl1995}.  This empirical approach does not describe actual mesoscale
899phenomena, but it is similar to the ensemble methods used to assess trajectory
900accuracy \citep{kahl1996, baumann1997, stohl1998}.
901
902\subsection{Moist convection}
903
904An important transport mechanism are the updrafts in convective clouds.  They
905occur in conjunction with downdrafts within the clouds and compensating
906subsidence in the cloud-free surroundings.  These convective transports are
907grid-scale in the vertical, but sub-grid scale in the horizontal, and are not
908represented by the ECMWF vertical velocity.
909
910To represent convective transport in a particle dispersion model, it is
911necessary to redistribute particles in the entire vertical column.  For
912FLEXPART we chose the convective parameterization scheme by
913\citet{emanuel1999}, as it relies on the grid-scale temperature and humidity
914fields and calculates a displacement matrix providing the necessary mass flux
915information for the particle redistribution.  The convective parameterization
916is switched on using \verb|lconvection| in file \verb|COMMAND|.  It's
917computation time scales to the square of the number of vertical model levels
918and may account for up to 70\% of FLEXPART's computation time using current
91960-level ECMWF data.
920
921The convection is computed within the subroutines {\tt convmix.f}, {\tt
922calcmatrix.f}, {\tt convect43c.f}, and {\tt redist.f}.  It is called every
923FLEXPART \verb|lsynctime| time step (typically 900~s) with time-interpolated
924temperature and specific humiditiy profiles from the ECMWF data.  Note that the
925original ECMWF model levels, not the Cartesian coordinates, are used in the
926convection scheme.  For efficiency reasons, particles are sorted according to
927their horizontal grid positions ({\tt sort2.f}) before calling the convection
928scheme once per grid column.
929
930In the Emanuel scheme ({\tt convect43c.f}), convection is triggered whenever
931
932\begin{equation}
933T_{vp}^{LCL+1} \ge T_{v}^{LCL+1} + T_{thres}
934\end{equation}
935
936with $T_{vp}^{LCL+1}$ the virtual temperature of a surface air parcel lifted to
937the level above the lifting condensation level $LCL$, $T_{v}^{LCL+1}$ the
938virtual temperature of the environment there, and $T_{thres}=0.9$~K a threshold
939temperature value.  Based on the buoyancy sorting principle \citep{emanuel1991,
940telford1975}, a matrix $MA$ of the saturated upward and downward mass fluxes
941within clouds is calculated by accounting for entrainment and detrainment:
942
943\begin{multline}
944MA^{i,j}=\\
945\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}|]}
946\label{matrix}
947\end{multline}
948
949Here $MA^{i,j}$ are the mass fractions displaced from level $i$ to level $j$,
950$M^{i}$ the mass fraction displaced from the surface to level $i$, $LNB$ the
951level of neutral buoyancy of a surface air parcel and $0 < \sigma^{i,j} < 1 $
952the mixing fraction between level $i$ and level $j$.  The fraction
953$\sigma^{i,j}$ is determined by the environmental potential temperature
954$\theta^{j}$, the liquid potential temperature $\theta_{l}^{i,j}$ of air
955displaced adiabatically from $i$ to $j$, and the liquid potential temperature
956$\theta_{lp}^{i,j}$ of an air parcel first lifted adiabatically to level $i$
957and further to level $j$:
958
959\begin{equation}
960\sigma^{i,j}=\frac{\theta^{j}-\theta_{lp}^{i,j}}{\theta_{l}^{i,j}-\theta_{lp}^{i,j}}
961\end{equation}
962
963By summing up over all levels $j$, we then calculate the saturated up- and
964downdrafts at each level $i$ from Eq.~\ref{matrix} and assume that these fluxes
965are balanced by a subsidence mass flux in the environment.
966
967The particles in each convectively active box are then redistributed ({\tt
968redist.f}) according to the matrix $MA$.  If the mass of an ECMWF model layer
969$i$ is $m^i$ and the mass flux from layer $i$ to layer $j$ accumulated over one
970time step is $\Delta MA^{i,j}$, then the probability of a particle to be moved
971from layer $i$ to layer $j$ is $\Delta MA^{i,j}/m^i$.  Whether a given particle
972is displaced or not is determined by drawing a random number between [0,1],
973which also determines the position of the particle within the destination layer
974$j$.  After the convective redistribution of the particles, the compensating
975subsidence mass fluxes are converted to a vertical velocity acting on those
976particles in the grid box that are not displaced by convective drafts.  By
977calculating a subsidence velocity rather than displacing particles randomly
978between layers the scheme's numerical diffusion in the cloud-free environment
979is eliminated.  The scheme was tested and fulfills the well-mixed criterion,
980i.e., if a tracer is well mixed in the whole atmospheric column, it remains so
981after the convection.
982
983\subsection{Particle splitting}
984
985During the initial phase of dispersion from a point source in the atmosphere,
986particles normally form a compact cloud.  Relatively few particles suffice to
987simulate this initial phase correctly.  After some time, however, the particle
988cloud gets distorted and particles spread over a much larger area.  More
989particles are now needed.  FLEXPART allows the user to specify a time constant
990$\Delta t_s$ (file \verb|COMMAND|).  Particles are split into two (each of
991which receives half of the mass of the original particle) after travel times of
992$\Delta t_s$, $2\Delta t_s$, $4\Delta t_s$, $8\Delta t_s$, and so on
993(subroutine \verb|timemanager.f|).\par
994
995\section{\label{backward}Forward and backward modeling}
996
997Normally, when FLEXPART is run forward in time (\verb|ldirect=1| in file
998\verb|COMMAND|), particles are released from one or a number of sources and
999concentrations are determined downwind on a grid.  However, FLEXPART can also
1000run backward in time (\verb|ldirect=-1|), which is more efficient than forward
1001modeling for calculating source-receptor relationships if the number of
1002receptors is smaller than the number of (potential) sources.  In the backward
1003mode, particles are released from a receptor location (e.g., a measurement
1004site) and a four-dimensional (3 space dimensions plus time) response function
1005(sensitivity) to emission input is calculated.
1006
1007Since this version (6.2) of FLEXPART, the calculation of the source-receptor
1008relationships is generalized for both forward and backward runs, allowing much
1009greater flexibility regarding the input and output units than before.
1010\verb|ind_source| and \verb|ind_receptor| in file \verb|COMMAND| switch between
1011mass and mass mixing ratio units at the source and at the receptor,
1012respectively.  Note that source always stands for the physical source and not
1013the location of the particle release, which is done at the source in forward
1014mode but at the receptor in backward mode.  Table ~\ref{units} gives an
1015overview of the units used in forward and backward modeling for different
1016settings of the above switches.  A "normal" forward simulation which specifies
1017the release in mass units (i.e., kg) and also samples the output in mass units
1018(i.e., a concentration in ng m$^{-3}$) requires both switches to be set to 1.
1019
1020\begin{table}
1021\setlength{\tabcolsep}{1.1mm}
1022\caption{\label{units}
1023Physical units of the input (in file RELEASES) and output data for forward
1024(files grid\_conc\_date) and backward (files grid\_time\_date) runs for the
1025various settings of the unit switches ind\_source and ind\_receptor (in both
1026switches 1 refers to mass units, 2 to mass mixing ratio units).} \vspace{3mm}
1027{\centerline{
1028\begin{tabular}{ccccc} \hline
1029Direction & ind\_source & ind\_receptor & input unit & output unit \\ \hline
1030Forward  & 1 & 1 & kg & ng m$^{-3}$ \\
1031Forward  & 1 & 2 & kg & ppt by mass \\
1032Forward  & 2 & 1 & 1 & ng m$^{-3}$ \\
1033Forward  & 2 & 2 & 1 & ppt by mass \\ 
1034Backward & 1 & 1 & 1 & s \\
1035Backward & 1 & 2 & 1 & s m$^3\,$kg$^{-1}$ \\
1036Backward & 2 & 1 & 1 & s kg m$^{-3}$ \\
1037Backward & 2 & 2 & 1 & s \\ \hline
1038\end{tabular}}}
1039\end{table}
1040
1041In the backward mode, any value not equal zero can be entered as the release
1042"mass" in file \verb|RELEASES| because the output is normalized by this value.
1043The calculated response function is related to the particles' residence time in
1044the output grid cells.  The unit of the output response function varies,
1045depending on how the switches are set.  If \verb|ind_source=1| and
1046\verb|ind_receptor=1|, the response function has the unit s.  If this response
1047function is folded (i.e., multiplied) with a 3-d field of emission mass fluxes
1048into the output grid boxes (in kg~m$^{-3}$s$^{-1}$), a concentration at the
1049receptor (kg~m$^{-3}$) is obtained.  If \verb|ind_source=1| and
1050\verb|ind_receptor=2|, the response function has the unit s~m$^3\,$kg$^{-1}$
1051and if it is folded with the emission mass flux (again in kg~m$^{-3}$s$^{-1}$),
1052a mass mixing ratio at the receptor is obtained.  The units of the response
1053function for \verb|ind_source=2| can be understood in analogy.
1054
1055In the case of loss processes (dry or wet deposition, decay) the response
1056function is ``corrected'' for these loss processes.  See \citet{seibert2001}
1057and, particularly, \citet{seibert2004} for a description of these generalized
1058in- and output options and the implementation of backward modeling in FLEXPART.
1059\citet{seibert2004} also describe the theory of backward modeling and give some
1060examples, and \citet{stohl2003} presents an application.
1061
1062\section{\label{plumetraj}Plume trajectories}
1063
1064In a recent paper, \citet{stohl2002} proposed a method to condense the complex
1065and large FLEXPART output using a cluster analysis \citep{dorling1992}.  The
1066idea behind this is to cluster, at every output time, the positions of all
1067particles originating from a release point, and write out only clustered
1068particle positions, along with additional information (e.g., fraction of
1069particles in the ABL and in the stratosphere).  This creates information that
1070is almost as compact as traditional trajectories but accounts for turbulence
1071and convection.  This option can be activated by setting \verb|iout| to 4 or 5
1072in file \verb|COMMAND|.  The number of clusters can be set with the parameter
1073\verb|ncluster| in file \verb|includepar|.  The clustering is handled and
1074output is produced by subroutine \verb|plumetraj.f|.
1075
1076\section{\label{removal}Removal processes}
1077
1078FLEXPART takes into account radioactive (or other) decay, wet deposition, and
1079dry deposition by reducing a particle's mass.  However, as atmospheric
1080transport is the same for all chemical species, a single particle can represent
1081several (up to \verb|maxspec|) chemical species, each affected differently by
1082the removal processes.
1083
1084\subsection{\label{radioactive}Radioactive decay}
1085
1086Radioactive decay is accounted for by reducing the particle mass according to
1087
1088\begin{equation}
1089m(t+\Delta t)=m(t) \exp (-\Delta t /\beta) \;,
1090\end{equation}
1091
1092where $m$ is particle mass, and the time constant $\beta=T_{1/2}/\ln(2)$ is
1093determined from the half life $T_{1/2}$ specified in file \verb|SPECIES_nnn|.
1094Deposited pollutant mass decays at the same rate.
1095
1096\subsection{OH Reaction}
1097
1098If a positive value for the OH reaction rate is given in the file
1099\verb|SPECIES_nnn|, OH reaction is performed in the model simulation and tracer
1100mass is lost by this reaction.  A monthly averaged 3\degreee$\times$ 5\degreen resolution OH
1101field averaged to 7 atmospheric levels is used.  The fields have been obtained
1102from the GEOS-CHEM model \citep{bey2001}.  The reaction rate is temperature
1103corrected and an activation rate of 1000 J\,mol$^{-1}$ is assumed.
1104
1105\subsection{\label{wetdepo}Wet deposition}
1106
1107Since version 8.0, in-cloud and below-cloud scavenging are treated differently.
1108Based on the humidity and temperature from the meteorological input data, the
1109occurence of clouds is calculated.  If relative humidity exceeds 80\% the
1110occurence of a cloud is assumed.
1111
1112{\bf In cloud scavenging} is treated differently for gases and particles. 
1113The implementation follows the scheme of \citet{hertel1995}.
1114
1115In general, the scavenging coefficient $\Lambda$ ((s$^{-1}$) depends on the precipitation
1116rate $I$ (mm\,h$^{-1}$) and the height $H_i$ over which scavenging takes place:
1117
1118\begin{equation}
1119\Lambda=\frac{S_i I}{H_i}
1120\end{equation}
1121
1122$S_i$ is different for gases and particles. For particles,
1123
1124\begin{equation}
1125S_i=0.9/cl
1126\end{equation}
1127
1128where $cl$ is the cloud liquid water content
1129
1130\begin{equation}
1131cl=2\times 10^{-7}\cdot I^{0.36} \;.
1132\end{equation}
1133
1134For gases,
1135
1136\begin{equation}
1137S_i=1/cl_\text{eff} \;,
1138\end{equation}
1139
1140where $cl_\text{eff}$ is an effective cloud liquid water content:
1141
1142\begin{equation}
1143cl_\text{eff}=\frac{(1-cl)}{H_\text{eff}RT}+cl \; .
1144\end{equation}
1145
1146{\bf Below cloud scavenging} takes the form of an exponential decay process
1147\citep{mcmahon1979}.  With the scavenging coefficient $\Lambda$, wet
1148deposition is described as
1149
1150\begin{equation}
1151m(t+\Delta t)=m(t) \exp (- \Lambda \Delta t) \;,
1152\end{equation}
1153
1154where $m$ is the particle mass.  $\Lambda$ increases with the precipitation
1155rate $I$ according to
1156
1157\begin{equation}
1158\Lambda=AI^{B} \;,
1159\label{wetscav}
1160\end{equation}
1161
1162where $A$ [s$^{-1}$] is the scavenging coefficient at $I=$1~mm/hour and $B$
1163gives the dependency on precipitation rate.  Both $A$ and $B$ must be specified
1164in file \verb|SPECIES_nnn|.  FLEXPART uses the same scavenging coefficients for
1165snow and rain.
1166
1167As wet deposition depends nonlinearly on precipitation rate, subgrid
1168variability of precipitation must be accounted for \citep{hertel1995}.  The
1169area fraction which experiences precipitation given a certain grid-scale
1170precipitation rate is calculated by
1171
1172\begin{equation}
1173F=\max\left[0.05,CC\frac{I_l fr_l(I_l) + I_c fr_c(I_c)}{I_l+I_c}\right] \;,
1174\end{equation}
1175
1176where $CC$ is the total cloud cover, $I_l$ and $I_c$ are the large scale and
1177convective precipitation rates, respectively, and $fr_l$ and $fr_c$ are
1178correction factors that depend on $I_l$ and $I_c$ (see Table~\ref{corrfacts}).
1179The subgrid scale precipitation rate is then $I_s=(I_l+I_c)/F \;$.
1180
1181\begin{table}
1182\setlength{\tabcolsep}{1.1mm}
1183\caption{\label{corrfacts}
1184Correction factors used for the calculation of the area fraction that
1185experiences precipitation.  Precipitation rates are in mm/hour.} \vspace{3mm}
1186{\centerline{
1187\begin{tabular}{cccccc} \hline
1188~&\multicolumn{5}{c|}{$I_l$ and $I_c$} \\ \hline
1189Factor&$I \le 1$&$1 < I \le 3$&$3 < I\le8$&$8 < I \le 20$&$20 < I$ \\ \hline
1190$fr_l$&0.50&0.65&0.80&0.90&0.95 \\
1191$fr_c$&0.40&0.55&0.70&0.80&0.90 \\ \hline
1192\end{tabular}}}
1193\end{table}
1194
1195\subsection{Dry deposition}
1196
1197Dry deposition is described in FLEXPART by a deposition velocity
1198
1199\begin{equation}
1200v_d(z)=-F_C/C(z) \;,
1201\end{equation}
1202
1203where $F_C$ and $C$ are the flux and the concentration of a species at height
1204$z$ within the constant flux layer.  A constant deposition velocity $v_d$ can
1205be set (file \verb|SPECIES_nnn|).  Alternatively, if the physical and chemical
1206properties of a substance are known (file \verb|SPECIES_nnn|), more complex
1207parameterizations for gases and particles are also available.
1208
1209\subsubsection{Dry deposition of gases}
1210
1211The deposition velocity of a gas is calculated with the {\em resistance method}
1212\citep{wesely1977} in subroutine \verb|getvdep.f| according to
1213
1214\begin{equation}
1215|v_d(z)|=\left[r_a(z)+r_b+r_c\right]^{-1} \;,
1216\end{equation}
1217
1218where $r_a$ is the aerodynamic resistance between $z$ and the surface, $r_b$ is
1219the quasilaminar sublayer resistance, and $r_c$ is the bulk surface resistance.
1220
1221The aerodynamic resistance $r_a$ is calculated in function \verb|raerod.f|
1222using the flux-profile relationship based on Monin-Obukhov similarity theory
1223\citep{stull1988}
1224
1225\begin{equation}
1226r_a(z)=\frac{1}{\kappa u_*}[\ln(z/z_0)- \Psi _h(z/L)+ \Psi _h(z_0/L)] \;.
1227\label{ra}
1228\end{equation}
1229
1230Following \citet{erisman1994}, the quasilaminar sublayer resistance is
1231
1232\begin{equation}
1233r_b=\frac{2}{\kappa u_*} \left(\frac{Sc}{Pr}\right)^{2/3} \;,
1234\end{equation}
1235
1236where $Sc$ and $Pr$ are the Schmidt and Prandtl numbers, respectively.  $Pr$ is
12370.72 and $Sc=\upsilon/D_i$, with $\upsilon$ being the kinematic viscosity of
1238air and $D_i$ being the molecular diffusivity of species $i$ in air.  The
1239slight dependency of $\upsilon$ on air temperature is formulated in accordance
1240with \citet{pruppacher1978}.  $r_b$ is calculated in function
1241\verb|getrb.f|.\par
1242
1243The surface resistance is calculated in function \verb|getrc.f| following
1244\citet{wesely1989} as
1245
1246\begin{equation}
1247\frac{1}{r_c}=
1248\frac{1}{r_s+r_m}+\frac{1}{r_{lu}}+\frac{1}{r_{dc}+r_{cl}}+\frac{1}{r_{ac}+r_{gs}} \;,
1249\end{equation}
1250
1251where $r_s$, $r_m$ and $r_{lu}$ represent the bulk values for leaf stomatal,
1252leaf mesophyll and leaf cuticle surface resistances (alltogether the upper
1253canopy resistance) , $r_{dc}$ represents the gas-phase transfer affected by
1254buoyant convection in canopies, $r_{cl}$ the resistance of leaves, twig, bark
1255and other exposed surfaces in the lower canopy, $r_{ac}$ the resistance for
1256transfer that depends only on canopy height and density, and $r_{gs}$ the
1257resistance for the soil, leaf litter, etc., at the ground.  Each of these
1258resistances is parameterized according to the species' chemical reactivity and
1259solubility, the landuse type, and the meteorological conditions.  The IGBP
1260landuse inventory \citep{belward1999} provides the area fractions of 13 landuse
1261classes for which roughness lengths $z_0$ are estimated, on a grid with 0.3
1262resolution (Table~\ref{landuses}).  Charnock's relationship \citep{stull1988}
1263$z_0=0.016u_*^2/g$ is used to calculate $z_0$ for the classes ``Ocean'' and
1264``Inland water'', because of its dependence on wave height.  Deposition
1265velocities are calculated for all landuse classes and weighted with their
1266respective areas.  The original resolution of the IGBP land cover
1267classification is 1 km x 1 km.  To save storage space, this was regridded to a
12680.3 x 0.3 degree resolution.  Table~\ref{conversion} shows how the initial 17
1269classse were transfered in the 13 classes used in the deposition scheme of
1270\cite{wesely1989}.  According to \cite{helmig2007} a resistance of snow to
1271ozone deposition of 10000 s m$^{-1}$ was used.  For SO$_2$ a value of 100
1272according to \citet{zhang2002} was used in the model.  As the Wesely scheme was
1273initially developed for North America, the rain forest was not represented
1274well.  Therefore a new category was introduced and resistance values according
1275to \cite{jacob1990} were used.  The resistance values are dependent on 5
1276different seasons.  As on the southern hemisphere they are asynchronous to the
1277northern hemisphere half a year was added when choosing the appropriate
1278seasonal category.  The snow depth is included in the original
1279\cite{wesely1989} parameterization to some extend, as at a certain time of the
1280year snow cover is assumed.  In order to have more accurate information we use
1281the snow cover based on the snow depth in the ECMWF fields.  If the snow cover
1282gets over 1 mm of water equivalent, the grid cell is considered as snow covered
1283and category 12 (snow and ice) is applied.
1284
1285\begin{table}
1286\setlength{\tabcolsep}{1.1mm}
1287\caption{\label{landuses}
1288List of the landuse classes and roughness lengths used by FLEXPART.
1289``Charnock'' indicates that Charnock's relationship is used to calculate the
1290roughness length.} \vspace{3mm}
1291{\centerline{
1292\begin{tabular}{lc} \hline
1293Urban land & 0.7        \\
1294Agricultural land       &  0.1 \\
1295Range land & 0.1 \\
1296Deciduous forest &  1.0 \\
1297Coniferous forest & 1.0 \\
1298Mixed forest including wetland &        0.7     \\
1299Water, both salt and fresh              &   Charnock \\
1300Barren land mostly desert &  0.01 \\
1301Nonforested wetland &   0.1 \\
1302Mixed agricultural and range land &     0.1 \\
1303Rocky open areas with low growing shrubs  &     0.05 \\
1304Snow and ice & 0.001 \\
1305Rainforest &   1.0 \\ \hline
1306\end{tabular}}}
1307\end{table}
1308
1309\begin{table*}
1310\setlength{\tabcolsep}{1.1mm}
1311\caption{\label{conversion}
1312Conversion from the IGBP land cover legend to the landuse categories used by Wesely}
1313\vspace{3mm}
1314{\centerline{
1315\begin{tabular}{clcl} \hline
1316\multicolumn{2}{c}{IGBP Land Cover Legend} & \multicolumn{2}{c}{Wesely} \\
1317Value & Description     &               Value & Description \\
13181&      Evergreen Needleleaf Forest&                    5&      Coniferous forest\\
13192&      Evergreen Broadleaf Forest      &            13& Rainforest\\
13203&      Deciduous Needleleaf Forest&                    4&      Deciduous forest\\
13214&      Deciduous Broadleaf Forest      &                       4&      Deciduous forest\\
13225&      Mixed Forest&6& Mixed forest including wetland\\
13236&      Closed Shrublands&            11&       rocky open areas with low growing shrubs\\
13247&      Open Shrublands&                   11&  rocky open areas with low growing shrubs\\
13258&      Woody Savannas&11&      rocky open areas with low growing shrubs\\
13269&      Savannas&11&rocky open areas with low growing shrubs\\
132710&     Grasslands                                      &               3       & Range land\\
132811&     Permanent Wetlands                         &      9     & nonforested wetland   \\
132912&     Croplands                                                  &    2       & Agricultural land\\
133013&     Urban and Built-Up                              &               1       & Urban land\\
133114&     Cropland/Natural Vegetation Mosaic&     10      & mixed agricultural and range land             \\
133215&     Snow and Ice                                                     &  12 & snow and ice\\
133316&     Barren or Sparsely Vegetated            &       8&      barren land mostly desert\\
133417&     Water Bodies                                                    &   7   & water, both salt and fresh\\
1335\end{tabular}}}
1336\end{table*}
1337
1338\subsubsection{Dry deposition of particulate matter}
1339
1340The deposition of particulates is calculated in subroutine \verb|partdep.f|
1341according to
1342
1343\begin{equation}
1344v_d(z)=\left[r_a(z)+r_b+r_a(z)r_bv_g\right]^{-1}+v_g \;,
1345\end{equation}
1346
1347where $v_g$ is the gravitational settling velocity calculated from
1348\citep{slinn1982}
1349
1350\begin{equation}
1351v_g=\frac{g \rho_p d_p^2C_{cun}}{18 \mu} \;,
1352\end{equation}
1353
1354where $\rho_p$ and $d_p$ are the particle density and diameter, $\mu$ the
1355dynamic viscosity of air (0.000018 kg m$^{-1}$s$^{-1}$) and $C_{cun}$ the
1356Cunningham slip-flow correction.  The quasilaminar sublayer resistance is
1357calculated from the same relationship as for gases, with an additional
1358impaction term.  For further details see \citet{slinn1982}.\par
1359
1360Settling and dry deposition velocities are strongly dependent on particulate
1361size.  FLEXPART assumes a logarithmic normal size distribution of the
1362particulate mass.  The user must specify the mean particulate diameter
1363$\overline{d_p}$ and a measure of the variation around $\overline{d_p}$,
1364$\sigma_p$.  Then, the settling and deposition velocities are calculated for
1365several particle diameters and are weighted with their respective particulate
1366mass fractions.
1367
1368Gravitational settling is important not only for the computation of the dry
1369deposition velocity, but also affects the particle's trajectory.  As a FLEXPART
1370particle can normally represent several species, gravitational settling can
1371only be taken into account correctly (i.e., influence particle trajectories) in
1372single-species simulations.  Pay attention that gravitational settling is
1373simply switched off in multi-species simulations, without warning.
1374
1375With FLEXPART version 8.2, the temperature dependence of the dynamic viscosity
1376is taken into account.  Furthermore, we have extended the calculation of
1377settling velocities to higher Reynolds numbers.  For this, an iterative
1378procedure has been introduced in subroutine \verb|get_settling.f|, where the
1379iteration is started with Stokes’ law settling \citep{naeslund1991}, and then
1380Reynolds number and settling velocity are calculated until convergence is
1381achieved.
1382
1383\subsubsection{Loss of particle mass due to dry deposition}
1384
1385The depositon velocity is calculated for a reference height (parameter
1386\verb|href| in file \verb|includepar|) of 15~m.  For all particles below
1387$2h_{ref}$, the mass lost by deposition is calculated by
1388
1389\begin{equation}
1390\Delta m(t)=m(t)\left[{1-\exp\left({\frac{-v_d(h_{ref})\Delta t}{2h_{ref}}}\right)}\right] \;.
1391\end{equation}
1392
1393\section{\label{conccalc}Calculation of concentrations, uncertainties, age spectra, and mass fluxes}
1394
1395Output quantities $C_{T_c}$ at time $T_c$ (output interval \verb|loutstep| is
1396set in file \verb|COMMAND|) are calculated as time-averages over period
1397$[T_c-\Delta T_c/2,T_C+\Delta T_c/2]$$\Delta T_c$ must be specified
1398(\verb|loutaver|) in file \verb|COMMAND|.  To calculate the time-averages,
1399concentrations $C_{T_s}$ at times $T_s$ within $[T_c-\Delta T_c/2,T_C+\Delta
1400T_c/2]$ are sampled at shorter intervals $\Delta T_s$ (\verb|loutsample| in
1401file \verb|COMMAND|) and are then divided by the number $N=\frac{\Delta
1402T_c}{\Delta T_s}$ of samples taken:
1403
1404\begin{equation}
1405C_{T_c}= \frac{1}{N} \sum_{i=1}^N {C_{T_s}} \;.
1406\end{equation}
1407
1408Both $\Delta T_c$ and $\Delta T_s$ must be multiples of the FLEXPART
1409synchronisation interval (\verb|lsynctime| in file \verb|COMMAND|).  The
1410shorter the sampling interval $\Delta T_s$, the more samples are taken and the
1411more accurate are thus the time-averaged concentrations.
1412
1413\subsection{Concentrations, mixing ratios, and emission response functions}
1414
1415The user can choose (\verb|iout| in file \verb|COMMAND|, which must be set to 1
1416for backward runs) whether concentrations, volume mixing ratios or both shall
1417be produced.  We shall use the term "concentration" and particle mass here, but
1418note that the actual units are determined by the settings of \verb|ind_source|
1419and \verb|ind_receptor|, according to Table~\ref{units}.  The concentration in
1420a grid cell is calculated in subroutine \verb|conccalc.f| by sampling the
1421tracer mass fractions of all particles within the grid cell and dividing by the
1422grid cell volume
1423
1424\begin{equation}
1425C_{T_s} =\frac{1}{V}\sum_{i=1}^N (m_i f_i) \;,
1426\end{equation}
1427
1428with $V$ being the grid cell volume, $m_i$ particle mass, $N$ the total number
1429of particles, and $f_i$ the fraction of the mass of particle $i$ attributed to
1430the respective grid cell.  This mass fraction is calculated by a uniform kernel
1431with bandwidths $(\Delta x,\Delta y)$, where $\Delta x$ and $\Delta y$ are the
1432grid distances on the longitude-latitude output grid.  Figure~\ref{kernel}
1433illustrates this: The particle is located at the center of the shaded rectangle
1434with side lengths $(\Delta x, \Delta y)$.  Generally, the shaded area stretches
1435over four grid cells, each of which receives a fraction of the particle's mass
1436equal to the fraction of the shaded area falling within this cell.  The uniform
1437kernel is not used during the first 3 hours after a particle's release (when
1438the mass is attributed only to the grid cell it resides in), in order to avoid
1439smoothing close to the source.
1440
1441\begin{figure}[htb]
1442\begin{minipage}[t]{2.8cm}
1443\end{minipage}\hfill
1444{\begin{minipage}[t]{12.5cm}
1445
1446\setlength{\unitlength}{2.5cm}
1447\begin{picture}(3.0,3.0)
1448
1449\thicklines
1450\multiput(0.5,0)(1,0){3}{\line(0,1){3.0}}
1451\multiput(0,0.5)(0,1){3}{\line(1,0){3.0}}
1452
1453\thinlines
1454\put(0.5,0.35){\vector(1,0){1.0}}
1455\put(1.5,0.35){\vector(-1,0){1.0}}
1456\put(0.35,0.5){\vector(0,1){1.0}}
1457\put(0.35,1.5){\vector(0,-1){1.0}}
1458
1459\put(0.9,0.15){$\Delta x$}
1460\put(0.05,0.92){$\Delta y$}
1461
1462\multiput(1.2909,1.25)(0.0909,0){10}{\line(0,1){1.0}}
1463\multiput(1.2,1.34909)(0,0.0909){10}{\line(1,0){1.0}}
1464
1465\thicklines
1466\put(1.67,1.76){\line(1,0){0.06}}
1467\put(1.70,1.73){\line(0,1){0.06}}
1468
1469\put(1.2,1.25){\line(1,0){1.0}}
1470\put(2.2,1.25){\line(0,1){1.0}}
1471\put(2.2,2.25){\line(-1,0){1.0}}
1472\put(1.2,2.25){\line(0,-1){1.0}}
1473
1474
1475\end{picture}
1476\end{minipage}}
1477\caption{\label{kernel} Illustration of the uniform kernel used to calculate
1478gridded concentration and deposition fields.  The particle position is marked
1479by ``{\rm +}''}
1480\end{figure}
1481
1482Wet and dry deposition fields are calculated on the same output grid
1483(subroutines \verb|wetdepokernel.f| and \verb|drydepokernel.f|) and are written
1484to all output grid files.  The deposited matter is accumulated over the course
1485of a model run, i.e.  it generally increases with model time.  However,
1486radioactive decay is calculated also for the deposited matter.
1487
1488\subsection{Uncertainties}
1489
1490The uncertainty of the output is estimated by carrying \verb|nclassunc| classes
1491of particles in the model simulation, and determining the concentration
1492separately for each class (subroutine \verb|conccalc.f|).  The standard
1493deviation, calculated from \verb|nclassunc| concentration estimates and divided
1494by $\sqrt{\tt{nclassunc}}$, is the standard deviation of the mean concentration
1495(subroutine \verb|concoutput.f|), which is also written to the output files for
1496every grid cell.  Note that the memory needed for some auxiliary fields
1497increases with \verb|nclassunc| {ít and} the number of age classes (see below).
1498It may, thus, be necessary to reduce \verb|nclassunc| for runs with large
1499output grids and age spectra calculations or in the backward mode.
1500
1501\subsection{Age spectra}
1502
1503The age spectra option is switched on using \verb|lagespectra| in file
1504\verb|COMMAND|, with the age classes specified in seconds in file
1505\verb|AGECLASSES|.  Concentrations are split into contributions from particles
1506of different age, defined as the time passed since their release.  Particles
1507are terminated once they are older than the oldest age class and their storage
1508space is made available to new particles.  Therefore, the age spectra option
1509can be used also with a single age class for defining a maximum particle age.
1510
1511\subsection{Parabolic kernel}
1512
1513In addition to the simple uniform kernel method, a computationally demanding
1514parabolic kernel as described in \citep{uliasz1994} can be used to calculate
1515surface concentrations for a limited number of receptor points (age spectra are
1516not available in this case):
1517
1518\begin{equation}
1519C_{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]\;,
1520\end{equation}
1521
1522where $h_{x_i}$, $h_{y_i}$ and $h_{z_i}$ are the kernel bandwidths which
1523determine the degree of smoothing, $r_x=(X_i-x)/h_{x_i}$,
1524$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
1525position of particle $i$.  The kernel bandwidths are a function of the
1526particles' age.
1527
1528\subsection{Mass fluxes}
1529
1530Mass flux calculations can be switched on using \verb|iflux| in file
1531\verb|COMMAND|.  Mass fluxes are calculated separately for eastward, westward,
1532northward, southward, upward and downward directions and contain both
1533grid-scale and subgrid-scale motions.  Mass fluxes are determined for the
1534centerlines of the output grid cells, e.g.  vertical fluxes are calculated for
1535motions across the half level of each output cell.
1536
1537\section{Domain-filling option}
1538
1539\subsection{General}
1540
1541If \verb|mdomainfill=1| in file \verb|COMMAND| particles are not released at
1542specific locations.  Instead, the longitudes and latitudes specified for the
1543first release in the \verb|RELEASES| file are used to set up a global or
1544limited model domain.  The particles (number is also taken from
1545\verb|RELEASES|) are then distributed in the model domain proportionally to air
1546density (subroutine \verb|init_domainfill.f|).  Each particle receives the same
1547mass, altogether accounting for the total atmospheric mass.  Subsequently,
1548particles move freely in the atmosphere.
1549
1550If a limited domain is chosen, mass fluxes are determined in small grid boxes
1551at the boundary of this domain (boundaries must be at least one grid box away
1552from the boundaries of the meteorological input data).  In the grid cells with
1553air flowing into the model domain, mass fluxes are accumulated over time and
1554whenever the accumulated mass exceeds the mass of a particle, a new particle
1555(or more, if required) is released at a randomly chosen position at the
1556boundary of the box (subroutine \verb|boundcond_domainfill.f|).  At the
1557outflowing boundaries particles are terminated.  Note that, due to the change
1558of mass of the atmosphere in the model domain and due to numerical effects, the
1559number of particles used is not exactly constant throughout the simulation.
1560
1561\subsection{Stratospheric ozone tracer}
1562
1563If \verb|mdomainfill=2|, the domain-filling option is used to simulate a
1564stratospheric ozone tracer.  Upon particle creation, the potential vorticity
1565(PV) at its position is determined by interpolation from the ECMWF data.
1566Particles initially located in the troposphere (PV$<$\verb|pvcrit| potential
1567vorticity units (pvu), default 2 pvu) are not used.  In contrast, stratospheric
1568particles (PV$>$\verb|pvcrit|) are given a mass according to:
1569
1570\begin{equation}
1571M_{O_3}=M_{air}\, P \, C \, 48/29
1572\end{equation}
1573
1574where $M_{air}$ is the mass of air a particle represents, $P$ is PV in pvu,
1575$C=60\times$10$^{-9}$ pvu$^{-1}$ is the ozone/PV relationship \citep{stohl2000}
1576(parameter \verb|ozonescale|), and the factor 48/29 converts from volume to
1577mass mixing ratio.  Particles are then allowed to advect through the
1578stratosphere and into the troposphere according to the winds.
1579
1580\section{Model output}
1581
1582Tracer concentrations and/or mixing ratios (for forward runs), or emission
1583sensitivity response functions (for backward runs) are calculated on a
1584three-dimensional longitude-latitude grid, defined in file \verb|OUTGRID|,
1585whose domain and resolution can differ from the grid on which meteorological
1586input data are given.  Two-dimensional wet and dry deposition fields are
1587calculated over the same spatial domain, and tracer mass fluxes can also be
1588determined on the 3-d grid.  Except for the mass fluxes, output can also be
1589produced on one nested output grid with higher horizontal but the same vertical
1590resolution, defined in file \verb|OUTGRID_NEST|.  For certain locations,
1591specified in file \verb|RECEPTORS|, concentrations can also be calculated
1592independently from a grid (see below).  The time interval (variable
1593\verb|loutstep|) at which output is produced is read in from file
1594\verb|COMMAND|.  For every output time and for every species (\verb|nnn|),
1595files are created, with file names ending with date, time and species number in
1596the format \verb|yyyymmddhhmmss_nnn|.  A list of all these output times is
1597written to the formatted file \verb|dates|.  The dates indicate the ending time
1598of an output sampling interval (see section~\ref{conccalc}).
1599
1600\subsection{Gridded output}
1601
1602There are several output options in FLEXPART, which can all be selected in file
1603\verb|COMMAND|.  Gridded output fields can be concentrations (files
1604\verb|grid_conc_date_nnn|), volume mixing ratios (files
1605\verb|grid_pptv_date_nnn|), emission response sensitivity in backward
1606simulations (files \verb|grid_time_date_nnn|), or fluxes (files
1607\verb|grid_flux_date|, unit 10$^{-12}$ kg m$^{-2}$ s$^{-1}$ for forward runs),
1608or, in backward mode, sensitivities to initial conditions taken from another
1609model (\verb|grid\_initial\_nnn|).
1610
1611The species number identifier \verb|nnn| starts at one (with leading zeros) and
1612increases to the maximum number of species used in the simulation.  Files
1613\verb|grid_conc_date_nnn| are created only in forward runs, whereas files
1614\verb|grid_time_date_nnn| are only created in backward runs.  Note that the
1615units of the files \verb|grid_conc_date_nnn| and \verb|grid_time_date_nnn|
1616depend on the settings of the switches \verb|ind_source| and
1617\verb|ind_receptor|, following Table~\ref{units}.  In particular, the units of
1618\verb|grid_conc_date| can also be mass mixing ratios.  For forward runs,
1619additional files \verb|grid_pptv_date_nnn| can be created, which contain volume
1620mixing ratios for gases.  Output files \verb|grid_conc*|, \verb|grid_pptv*|,
1621and \verb|grid_time*| also contain wet and dry deposition fields (unit
162210$^{-12}$ kg m$^{-2}$ in forward mode), and all files contain, for each grid
1623cell, corresponding uncertainties.  All these file types share a common header,
1624file \verb|header| produced by subroutine \verb|writeheader.f|, where important
1625information on the model run (start of simulation, grid domain, number and
1626position of vertical levels, age classes, release points, etc.) is stored.  In
1627all postprocessing programs, the header must be read in before the actual data
1628files.  File names for the output nests follow the same nomenclature as
1629described above, but with \verb|_nest| added (e.g., \verb|header_nest|, or
1630\verb|grid_conc_nest_date_nnn|).  The output files are written with subroutines
1631\verb|concoutput.f| and \verb|fluxoutput.f|.
1632
1633FLEXPART output typically contains many grid cells with zero values.  It would
1634be inefficient to write out all these zeroes.  Therefore, a special format has
1635been designed that compresses the information to the relevant information.  In
1636previous versions (up to version 7), the output consisted of a mixture of a
1637full grid dump (including zeroes) and a sparse matrix format - output was
1638switched between these two formats based on which one was smaller.
1639
1640In version 8.0, the output has been redesigned completely for yet more efficiency
1641such that output file sizes are only about 40\% of what they used to be.  The
1642output grid is searched for consecutive sequences of non-zero values.  The
1643variable \verb|sp_count_i| gives the number of such sequences (n), and the
1644integer field \verb|sparse_dump_i(n)| contains the field positions of the first
1645non-zero element for every sequence.  The variable \verb|sp_count_r| gives the
1646total number (k) of non-zero values written out.  They are contained in the
1647real field \verb|sparse_dump_r(n)|.  Since all physical output quantities of
1648FLEXPART are non-zero, sequences are written out alternatingly as positive or
1649negative values.  Every switch between positive and negative values indicates
1650that a new sequence with non-zero values starts.  The field position for that
1651start is contained in \verb|sparse_dump_i(k)| for the k-th switch between
1652positive and negative.  Zero values in between sequences are ignored and not
1653written out.  Field positions within the 3-d output field are coded such that a
1654single integer value is sufficient.  It can later be converted back to give
1655positions in all three coordinates.
1656
1657The possibility to output the sensitivities to initial conditions for backward
1658runs has been introduced with FLEXPART version 8.2. This option can be used to
1659calculate the exact sensitivities of the concentrations (or mixing ratios,
1660depending on what is simulated) to initial conditions either in concentration
1661or mixing ratio units taken from a gridded data set, which can be produced
1662either by another model or by a FLEXPART forward simulation.  This option has
1663been introduced to allow interfacing FLEXPART with other models.  Multiplying
1664the sensitivities in files \verb|grid_initial_nnn| with the corresponding
1665concentrations (or mixing ratios, depending on option chosen for switch
1666\verb|LINIT_COND| in file \verb|COMMAND|) from the forward simulation at the
1667interfacing time gives the concentration (or mixing ratio) response at the
1668receptor (taking into account possible loss processes during the FLEXPART
1669simulation time).
1670
1671\subsection{Receptor point output}
1672
1673For a list of points at the surface, concentrations or mixing ratios in forward
1674simulations can be determined with a grid-independent method.  This information
1675is written to files \verb|receptor_conc| and \verb|receptor_pptv|,
1676respectively, for all dates of a simulation.
1677
1678\subsection{Particle dump and warm start option}
1679
1680Particle information (3-d position, release time, release point, and release
1681masses for all species) can be written out to files (subroutine
1682\verb|partoutput.f|) either continuously (binary files \verb|partposit_date|),
1683or `only at the end' of a simulation (file \verb|partposit_end|).  In both
1684cases output is written every output interval but file \verb|partposit_end| is
1685overwritten upon each new output.  If FLEXPART must be terminated, it can be
1686continued later on by reading in files \verb|header| and \verb|partposit_end|
1687produced by the previous run (subroutine \verb|readpartpositions.f|).  Such a
1688warm start is done if variable \verb|ipin| is set to 1 in file \verb|COMMAND|.
1689
1690If option \verb|mquasilag| is chosen in file \verb|COMMAND|, particle dumps
1691every output interval are produced in a very compact format by converting the
1692positions to an \verb|integer*2| format (subroutine \verb|partoutput_short.f|).
1693As some accuracy is lost in the conversion, this output is not used for the
1694warm start option.  Another difference to the normal particle dump is that
1695every particle gets a unique number, thus allowing postprocessing routines to
1696identify continuous particle trajectories.
1697
1698\subsection{Clustered plume trajectories}
1699
1700Condensed particle output using the clustering algorithm described in
1701section~\ref{plumetraj} is written to the formatted file
1702\verb|trajectories.txt|.  Information on the release points (coordinates,
1703release start and end, number of particles) is written by subroutine
1704\verb|openouttraj.f| to the beginning of file \verb|trajectories.txt|.
1705Subsequently, \verb|plumetraj.f| writes out a time sequence of the clustering
1706results for each release point: release point number, time in seconds elapsed
1707since the middle of the release interval, plume centroid position coordinates,
1708various overall statistics (e.g., fraction of particles residing in the ABL and
1709troposphere), and then for each cluster the cluster centroid position, the
1710fraction of particles belonging to the cluster, and the root-mean-square
1711distance of cluster member particles from the cluster centroid.
1712
1713\section{pflexpart: a Python routine for the analysis of FLEXPART output}
1714
1715To assist in the usage and analysis of FLEXPART data we have created a Python
1716module that is available with this new release.  The Python module 'pflexpart'
1717enables the user to easily read and access the header and grid output data of
1718the FLEXPART model runs.  Furthermore, we provide some basic classes that
1719assist in conducted standard analysis of backward and forward model runs.
1720
1721The module is released under the same GPL license as FLEXPART.  As open source
1722code it is constantly undergoing revision and updates from the community.
1723Thus, the core functionality of the module is described online.  See:
1724http://transport.nilu.no/pflexpart
1725
1726The module is largely based on the well known Numpy and matplotlib Python
1727modules as well as the {``basemap``} tool box of matplotlib.  Additionally, for
1728efficiency in reading the data, there are several custom build modules that use
1729the \verb|readgrid.f| FORTRAN routines described above.  Note, that some of
1730these modules will require the user to compile them for their system to achieve
1731maximum benefit.
1732
1733The central framework to the module is the pflexpart \verb|Header| class.  This
1734class will read a FLEXPART header file and provide some functionality toward
1735data analysis.  For instance, the \verb|Header| class has a
1736\verb|fill_backward| method for backward runs that will calculate the Total
1737Column and Footprint sensitivities from the N ageclasses used in the run.  From
1738this method, plotting of residence times and emission sensitivities is
1739relatively simple.
1740
1741Other features include wrapper functions around the matplotlib toolboxes that
1742allow for plotting of the data easily.  The functions \verb|plot_totalcolumn|
1743and \verb|plot_footprint| for example are customized wrappers of the
1744\verb|plot_sensitivity| function.  These functions take data objects that are
1745created by the \verb|fill_backward| method of the \verb|Header| class and allow
1746the user to create plots quickly.  For forward runs, the \verb|Header.readgrid|
1747method can be used, again with the \verb|plot_sensitivity| function of the
1748pflexpart module.
1749
1750For more details, the reader is referred again to the more frequently updated
1751module home page.
1752
1753\section{Final remark}
1754
1755In this note, we have described the Lagrangian particle dispersion model FLEXPART in version
17568.2.  As FLEXPART is developed further this text will be kept up to date and will be accessible
1757from the FLEXPART home page at http://transport.nilu.no/flexpart.
1758
1759\begin{thebibliography}{}
1760
1761\bibitem[Asman(1995)]{asman1995}
1762Asman, W. A. H.:
1763Parameterization of below-cloud scavenging of highly soluble gases under convective conditions.
1764Atmos. Environ., 29, 1359--1368, 1995.
1765\bibitem[Baumann and Stohl(1997)]{baumann1997}
1766Baumann, K. and Stohl, A.:
1767Validation of a long-range trajectory model using gas balloon tracks from the Gordon Bennett Cup 95.
1768J. Appl. Meteor., 36, 711--720, 1997.
1769\bibitem[Beljaars and Holtslag(1991)]{beljaars1991}
1770Beljaars, A. C. M. and Holtslag, A. A. M.:
1771Flux parameterization over land surfaces for atmospheric models.
1772J. Appl. Meteor., 30, 327--341, 1991.
1773\bibitem[Belward et al.(1999)]{belward1999}
1774Belward, A.S., Estes, J.E., and Kline, K.D.:
1775The IGBP-DIS 1-Km Land-Cover Data Set DISCover: A Project Overview.
1776Photogrammetric Engineering and Remote Sensing , v. 65, no. 9, p. 1013--1020, 1999.
1777\bibitem[Berkowicz and Prahm(1982)]{berkowicz1982}
1778Berkowicz, R. and Prahm, L. P.:
1779Evaluation of the profile method for estimation of surface fluxes of momentum and heat.
1780Atmos. Environ., 16, 2809--2819, 1982.
1781\bibitem[Bey et al.(2001)]{bey2001}
1782Bey I, Jacob DJ, Yantosca RM, et al.:
1783Asian chemical outflow to the Pacific in spring: Origins, pathways, and budgets.
1784J. Geophys. Res., 106, 19, 23073-23095, 2001.
1785\bibitem[Businger et al.(1971)]{businger1971}
1786Businger, J. A., Wyngaard, J. C., Izumi, Y. and Bradley, E. F.:
1787Flux-profile relationships in the atmospheric surface layer.
1788J. Atmos. Sci., 28, 181--189, 1971.
1789\bibitem[Dorling et al.(1992)]{dorling1992}
1790 Dorling, S. R., Davies, T. D. and Pierce, C.E.:
1791Cluster analysis: a technique for estimating the synoptic meteorological controls on air and precipitation chemistry - method and applications.
1792Atmos. Environ., 26A, 2575--2581, 1992.
1793\bibitem[ECMWF(1995)]{ecmwf1995}
1794ECMWF,
1795User Guide to ECMWF Products 2.1. Meteorological Bulletin M3.2. Reading, UK, 1995.
1796\bibitem[Emanuel(1991)]{emanuel1991}
1797Emanuel, K. A.:
1798A scheme for representing cumulus convection in large-scale models,
1799J. Atmos. Sci., 48, 2313--2335, 1991.
1800\bibitem[Emanuel and \v{Z}ivkovi\'{c}-Rothman(1999)]{emanuel1999}
1801Emanuel, K. A., and \v{Z}ivkovi\'{c}-Rothman, M.:
1802Development and evaluation of a convection scheme for use in climate models.
1803J. Atmos. Sci., 56, 1766--1782, 1999.
1804\bibitem[Erisman et al.(1994)]{erisman1994}
1805Erisman, J. W., Van Pul, A. and Wyers, P.:
1806Parametrization of surface resistance for the quantification of atmospheric deposition of acidifying pollutants and ozone.
1807Atmos. Environ., 28, 2595--2607, 1994.
1808\bibitem[Flesch et al.(1995)]{flesch1995}
1809Flesch, T. K., Wilson, J. D., and Lee, E.:
1810Backward-time Lagrangian stochastic dispersion models and their application to estimate gaseous emissions,
1811J. Appl. Meteorol., 34, 1320--1333, 1995.
1812\bibitem[Forster et al.(2001)]{forster2001}
1813Forster, 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.:
1814Transport of boreal forest fire emissions from Canada to Europe.
1815J. Geophys. Res., 106, 22,887--22,906, 2001.
1816\bibitem[Gupta et al.(1997)]{gupta1997}
1817Gupta, S., McNider, R. T., Trainer, M., Zamora, R. J., Knupp, K. and Singh, M.P.:
1818Nocturnal wind structure and plume growth rates due to inertial oscillations.
1819J. Appl. Meteor., 36, 1050--1063, 1997.
1820\bibitem[Hanna(1982)]{hanna1982}
1821Hanna, S.R.:
1822Applications in air pollution modeling. In: Nieuwstadt F.T.M. and H. van Dop (ed.): {\em Atmospheric Turbulence and Air Pollution Modelling}.
1823D. Reidel Publishing Company, Dordrecht, Holland, 1982.
1824\bibitem[Helmig et al.(2007)]{helmig2007}
1825Helmig D., Ganzeveld, L., Butler, T. and Oltmans S.J.:
1826The role of ozone atmosphere-snow gas exchange on polar, boundary-layer tropospheric ozone - a review and sensitivity analysis
1827Atmos. Chem. and Physics, 7, 15--30, 2007.
1828\bibitem[Hertel et al.(1995)]{hertel1995}
1829Hertel, O., Christensen, J. Runge, E. H., Asman, W. A. H., Berkowicz, R., Hovmand, M. F. and Hov, O.:
1830Development and testing of a new variable scale air pollution model - ACDEP.
1831Atmos. Environ., 29, 1267--1290, 1995.
1832\bibitem[Hubbe et al.(1997)]{hubbe1997}
1833Hubbe, J. M., Doran, J. C., Liljegren, J.C. and Shaw, W. J.:
1834Observations of spatial variations of boundary layer structure over the southern Great Plains cloud and radiation testbed.
1835J. Appl. Meteor., 36, 1221--1231, 1997.
1836\bibitem[Jacob and Wofsy(1990)]{jacob1990}
1837Jacob D. J., Wofsy W. C.:
1838Budgets of reactive nitrogen, hydrocarbons, and ozone over the amazon-forest during the wet season.
1839J. Geophys. Res., 95, 16737--16754, 1990.
1840\bibitem[Kahl (1996)]{kahl1996}
1841Kahl, J. D. W.:
1842On the prediction of trajectory model error.
1843Atmos. Environ., 30, 2945--2957, 1996.
1844\bibitem[Legg and Raupach(1982)]{legg1982}
1845Legg, B. J., and Raupach, M. R.:
1846Markov-chain simulation of particle dispersion in inhomogeneous flows: the mean drift velocity induced by a gradient in Eulerian velocity variance.
1847Bound.-Layer Met., 24, 3--13, 1982.
1848\bibitem[Legras et al.(2003)]{legras2003}
1849Legras, B., Joseph, B. and Lefevre, F.:
1850Vertical diffusivity in the lower stratosphere from Lagrangian back-trajectory reconstructions of ozone profiles.
1851J. Geophys. Res., 108, 4562, doi:10.1029/2002JD003045, 2003.
1852\bibitem[Maryon(1998)]{maryon1998}
1853Maryon, R. H.:
1854Determining cross-wind variance for low frequency wind meander.
1855Atmos. Environ., 32, 115--121, 1998.
1856\bibitem[McMahon(1979)]{mcmahon1979}
1857McMahon, T. A., and Denison, P. J.:
1858Empirical atmospheric deposition parameters - a survey.
1859Atmos. Environ., 13, 571--585, 1979.
1860\bibitem[McNider et al.(1988)]{mcnider1988}
1861McNider, R. T., Moran, M. D. and Pielke, R. A.:
1862Influence of diurnal and inertial boundary layer oscillations on long-range dispersion.
1863Atmos. Environ., 22, 2445--2462, 1988.
1864\bibitem[Naeslund and Thaning(1991)]{naeslund1991}
1865Naeslund, E., and Thaning, L.:
1866On the settling velocity in a nonstationary atmosphere.
1867Aerosol Sci. Technol., 14, 247-256, 1991.
1868\bibitem[Petterssen(1940)]{petterssen1940}
1869Petterssen, S.:
1870Weather Analysis and Forecasting.
1871McGraw-Hill Book Company, New York. pp. 221--223, 1940.
1872\bibitem[Pruppacher and Klett(1978)]{pruppacher1978}
1873Pruppacher, H. R. and Klett, J. D.:
1874Microphysics of clouds and precipitation.
1875D. Reidel Publishing Company, Dordrecht, 714p., 1978.
1876\bibitem[Rodean(1996)]{rodean1996}
1877Rodean, H.:
1878Stochastic Lagrangian models of turbulent diffusion.
1879Meteorological Monographs, 26 (48). American Meteorological Society, Boston, USA, 1996.
1880\bibitem[Ryall et al.(1997)]{ryall1997}
1881Ryall, D. B. and Maryon, R. H.:
1882Validation 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.
1883\bibitem[Seibert(2001)]{seibert2001}
1884Seibert, P.:
1885Inverse 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.
1886\bibitem[Seibert and Frank(2004)]{seibert2004}
1887Seibert, P. and Frank, A.:
1888Source-receptor matrix calculation with a Lagrangian particle dispersion model in backward mode.
1889Atmos. Chem. Phys., 4, 51--63, 2004.
1890\bibitem[Seibert et al.(2001)]{seibertetal2001}
1891Seibert, P., Kr\"uger, B. and Frank, A.:
1892Parametrisation of convective mixing in a Lagrangian particle dispersion model, Proceedings of the 5th GLOREAM Workshop, Wengen, Switzerland, September 24--26, 2001.
1893\bibitem[Slinn(1982)]{slinn1982}
1894Slinn, W. G. N.:
1895Predictions for particle deposition to vegetative canopies.
1896Atmos. Environ., 16, 1785--1794, 1982.
1897\bibitem[Spichtinger et al.(2001)]{spichtinger2001}
1898Spichtinger, N., Wenig, M., James, P., Wagner, T., Platt, U. and Stohl, A.:
1899Satellite detection of a continental-scale plume of nitrogen oxides from boreal forest fires,
1900Geophys. Res. Lett., 28, 4579--4582, 2001.
1901\bibitem[Stohl(1998)]{stohl1998}
1902Stohl, A.:
1903Computation, accuracy and applications of trajectories -- a review and bibliography.
1904Atmos. Environ., 32, 947--966, 1998.
1905\bibitem[Stohl et al.(2005)]{stohl2005}
1906Stohl, A., Forster, C., Frank, A., Seibert, P., Wotawa, G.:
1907Technical note: The Lagrangian particle dispersion model FLEXPART version 6.2.,
1908Atmos. Chem. Phys., 5, 2461-2474, 2005
1909\bibitem[Stohl et al.(2004)]{stohl2004}
1910Stohl, A., Cooper, O. R., Damoah, R., Fehsenfeld, F. C., Forster, C., Hsie, E.-Y., Hübler, G., Parrish,
1911D. D. and Trainer, M.:
1912Forecasting for a Lagrangian aircraft campaign.
1913Atmos. Chem. Phys., 4, 1113--1124, 2004.
1914\bibitem[Stohl et al.(2002)]{stohl2002}
1915Stohl, A., Eckhardt, S., Forster, C., James, P., Spichtinger, N. and Seibert, P.:
1916A replacement for simple back trajectory calculations in the interpretation of atmospheric trace substance measurements.
1917Atmos. Environ., 36, 4635--4648, 2002.
1918\bibitem[Stohl et al.(2003)]{stohl2003}
1919Stohl, A., Forster, C., Eckhardt, S., Spichtinger, N., Huntrieser, H., Heland, J., Schlager, H., Wilhelm, S., Arnold, F. and Cooper, O.:
1920A backward modeling study of intercontinental pollution transport using aircraft measurements.
1921J. Geophys. Res., 108, 4370, doi:10.1029/2002JD002862, 2003.
1922\bibitem[Stohl et al.(1998)]{stohletal1998}
1923Stohl, A., Hittenberger, M. and Wotawa, G.:
1924Validation of the Lagrangian particle dispersion model FLEXPART against large scale tracer experiment data.
1925Atmos. Environ., 32, 4245--4264, 1998.
1926\bibitem[Stohl et al.(2000)]{stohl2000}
1927Stohl, A., Spichtinger-Rakowsky, N., Bonasoni, P., Feldmann, H., Memmesheimer, M., Scheel, H. E., Trickl, T., H\"ubener, S. H., Ringer, W. and Mandl, M.:
1928The influence of stratospheric intrusions on alpine ozone concentrations.
1929Atmos. Environ., 34, 1323--1354, 2000.
1930\bibitem[Stohl and Thomson(1999)]{stohlthomson1999}
1931Stohl, A. and Thomson, D. J.:
1932A density correction for Lagrangian particle dispersion models.
1933Bound.-Layer Met., 90, 155--167, 1999.
1934\bibitem[Stohl and Trickl(1999)]{stohl1999}
1935Stohl, A. and Trickl, T.:
1936A textbook example of long-range transport: Simultaneous observation of ozone maxima of stratospheric and North American origin in the free troposphere over Europe,
1937J. Geophys. Res., 104, 30,445--30,462, 1999.
1938\bibitem[Walmsley and Wesely(1996]{walmsley1996}
1939Walmsley J.L., Wesely M.L.:
1940Modification of coded parametrizations of surface resistances to gaseous dry deposition. 
1941Atmos. Environ., 30, 1181--1188, 1996.
1942\bibitem[Stohl et al.(1995)]{stohl1995}
1943Stohl, A., Wotawa, G., Seibert, P. and Kromp-Kolb, H.:
1944Interpolation errors in wind fields as a function of spatial and temporal resolution and their impact on different types of kinematic trajectories.
1945J. Appl. Meteor., 34, 2149--2165, 1995.
1946\bibitem[Stull(1988)]{stull1988}
1947Stull, R. B.:
1948An Introduction to Boundary Layer Meteorology. Kluwer Academic Publishers, Dordrecht, 1988.
1949\bibitem[Telford(1975)]{telford1975}
1950Telford, J. W.:
1951Turbulence, entrainment and mixing in cloud dynamics,
1952Pure Appl. Geophys., 113, 1067--1084, 1975.
1953\bibitem[Thomson(1987)]{thomson1987}
1954Thomson D. J.:
1955Criteria for the selection of stochastic models of particle trajectories in turbulent flows.
1956J. Fluid Mech., 180, 529--556, 1987.
1957\bibitem[Uliasz(1994)]{uliasz1994}
1958Uliasz, M.:
1959Lagrangian particle dispersion modeling in mesoscale applications. In: Zannetti, P. (ed.): {\it Environmental Modeling, Vol. II}. Computational Mechanics Publications, Southampton, UK, 1994.
1960\bibitem[Velde et al.(1994)]{velde1994}
1961Velde 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.:
1962The Preparation of a European Land Use Database. National Institute of Public Health and Environmental Protection, Report nr 712401001, Bilthoven, The Netherlands, 1994.
1963\bibitem[Vogelezang and Holtslag(1996)]{vogelezang1996}
1964Vogelezang, D. H. P. and Holtslag, A. A. M.:
1965Evaluation and model impacts of alternative boundary-layer height formulations.
1966Bound.-Layer Met., 81, 245--269, 1996.
1967\bibitem[Wesely and Hicks(1977)]{wesely1977}
1968Wesely, M. L. and Hicks, B. B.:
1969Some factors that affect the deposition rates of sulfur dioxide and similar gases on vegetation.
1970J. Air Poll. Contr. Assoc., 27, 1110--1116, 1977.
1971\bibitem[Wesely(1989)]{wesely1989}
1972Wesely, M. L.:
1973Parameterization of surface resistances to gaseous dry deposition in regional-scale numerical models.
1974Atmos. Environ., 23, 1293--1304, 1989.
1975\bibitem[Wilson and Flesch(1993)]{wilson1993}
1976Wilson, D. J. and Flesch, T. K.:
1977Flow boundaries in random-flight dispersion models: enforcing the well-mixed condition.
1978J. Appl. Meteor., 32, 1695--1707, 1993.
1979\bibitem[Wilson et al.(1983)]{wilson1983}
1980Wilson, J. D., Legg, B. J. and Thomson, D. J.:
1981Calculation of particle trajectories in the presence of a gradient in turbulent-velocity scale.
1982Bound.-Layer Met., 27, 163--169, 1983.
1983\bibitem[Wilson and Sawford(1996)]{wilson1996}
1984Wilson, J. D. and Sawford, B. L.:
1985Review of Lagrangian stochastic models for trajectories in the turbulent atmosphere.
1986Bound.-Layer Met., 78, 191--210, 1996.
1987\bibitem[Wotawa et al.(1996)]{wotawa1996}
1988Wotawa, G., Stohl, A. and Kromp-Kolb, H.:
1989Parameterization of the planetary boundary layer over Europe - a data comparison between the observation based OML preprocessor and ECMWF model data.
1990Contr. Atmos. Phys., 69, 273--284, 1996.
1991\bibitem[Wotawa and Stohl(1997)]{wotawa1997}
1992Wotawa, G. and Stohl, A.:
1993Boundary layer heights and surface fluxes of momentum and heat derived from ECMWF data for use in pollutant dispersion models -- problems with data accuracy.
1994In: 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.
1995\bibitem[Zannetti(1992)]{zannetti1992}
1996Zannetti, P.:
1997Particle Modeling and Its Application for Simulating Air Pollution Phenomena. In: Melli P. and P. Zannetti (ed.): Environmental Modelling. Computational Mechanics Publications, Southampton, UK, 1992.
1998\bibitem[Zhang et al.(2002)]{zhang2002}
1999Zhang L.M., Moran M.D., Makar P.A., Brook J.R. and S. Gong
2000Modelling gaseous dry deposition in AURAMS: a unified regional air-quality modelling system
2001Atmos. Environ., 36, 537-560, 2002.
2002
2003\end{thebibliography}
2004
2005\newpage
2006
2007\appendix
2008
2009\onecolumn
2010
2011\section{FLEXPART sample input files}
2012\subsection{The pathnames file}
2013A file \verb|pathnames| must exist in the directory where FLEXPART is started.
2014It states the pathnames (absolute or relative) of input and output files:
2015\begin{footnotesize}\begin{verbatim}
2016/home/as/FLEXPART50/options/
2017/volc/as/contrace/modelresults/forward/
2018/volc/windcontrace/
2019/volc/windcontrace/AVAILABLE
2020/volc/nested/
2021/volc/nested/AVAILABLE
2022============================================
2023
2024Line 1: path where control files "COMMAND" and "RELEASES" are available
2025Line 2: name of directory where output files are generated
2026Line 3: path where meteorological fields are available (mother grid)
2027Line 4: full filename of "AVAILABLE"-file (mother grid)
2028
2029Subsequent lines:
2030Line 2n+3: path where meteorological fields are available (nested grid n)
2031Line 2n+4: full filename of "AVAILABLE"-file (nested grid n)
2032
2033Line below last pathname must be:
2034============================================
2035
2036The grids must be arranged such as that the coarse-scale nests
2037come before the fine-scale nests. Multiple nests of the same
2038nesting level are allowed. In that case, the order is arbitrary.
2039\end{verbatim}\end{footnotesize}
2040
2041\newpage
2042
2043\subsection{Files in directory windfields}
2044
2045The directory where the meteorological input data are stored, here called
2046\verb|windfields| (\verb|/volc/windcontrace/| in the above example
2047\verb|pathnames| file), contains grib-code files containing the ECMWF data.
2048All meteorological fields must have the same structure, i.e.  the same
2049computational domain and the same resolution.  An example listing of this
2050directory is given below.
2051\begin{footnotesize}\begin{verbatim}
2052AVAILABLE       EN01102806      EN01102815
2053EN01102800      EN01102809      EN01102818
2054EN01102803      EN01102812      EN01102821
2055\end{verbatim}\end{footnotesize}
2056
2057The 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.
2058\begin{footnotesize}\begin{verbatim}
2059DATE     TIME         FILENAME     SPECIFICATIONS
2060YYYYMMDD HHMISS
2061________ ______      __________      __________
206220011028 000000      EN01102800      ON DISC
206320011028 030000      EN01102803      ON DISC
206420011028 060000      EN01102806      ON DISC
206520011028 090000      EN01102809      ON DISC
206620011028 120000      EN01102812      ON DISC
206720011028 150000      EN01102815      ON DISC
206820011028 180000      EN01102818      ON DISC
206920011028 210000      EN01102821      ON DISC
2070\end{verbatim}\end{footnotesize}
2071
2072Since version 7.0, FILENAME can be up to 256 characters long.  Thus, the path
2073to a file can be directly placed in the AVAILABLE file.  If this feature is
2074used, line 3 in the \verb|pathnames| file should be left empty.
2075\begin{footnotesize}\begin{verbatim}
2076DATE     TIME         FILENAME     SPECIFICATIONS
2077YYYYMMDD HHMISS
2078________ ______      __________      __________
207920011028 000000      /work_disk/data/EN01102800      ON DISC
208020011028 030000      /work_disk/data/EN01102803      ON DISC
208120011028 060000      /work_disk/data/EN01102806      ON DISC
208220011028 090000      /work_disk/data/EN01102809      ON DISC
208320011028 120000      /work_disk/data/EN01102812      ON DISC
208420011028 150000      /work_disk/data/EN01102815      ON DISC
208520011028 180000      /work_disk/data/EN01102818      ON DISC
208620011028 210000      /work_disk/data/EN01102821      ON DISC
2087\end{verbatim}\end{footnotesize}
2088
2089Nested wind fields must be stored in one or more different directory/ies, as specified in the \verb|pathnames| file.
2090
2091\newpage
2092
2093\subsection{Files in directory options}
2094
2095The files in directory \begin{footnotesize}\verb|options|\end{footnotesize} are used to specify the model run.
2096An example listing of \begin{footnotesize}\verb|options|\end{footnotesize} is given below.
2097\begin{footnotesize}\begin{verbatim}
2098AGECLASSES           IGBP_int1.dat       RECEPTORS             SPECIES
2099COMMAND              OH_7lev_agl.dat     RELEASES              surfdata.t
2100COMMAND.alternative  OUTGRID             RELEASES.alternative  surfdepo.t
2101COMMAND.reference    OUTGRID_NEST        RELEASES.reference             
2102\end{verbatim}\end{footnotesize}
2103
2104Here, SPECIES is a subdirectory containing a number of files (this feature was introduced in version 8.0):
2105\begin{footnotesize}\begin{verbatim}
2106SPECIES:
2107SPECIES_001  SPECIES_004  SPECIES_007
2108SPECIES_002  SPECIES_005  SPECIES_008
2109SPECIES_003  SPECIES_006  SPECIES_009
2110\end{verbatim}\end{footnotesize}
2111
2112\subsubsection{File COMMAND}
2113
2114The most important configuration file for setting up a FLEXPART simulation is
2115the \verb|COMMAND| file which specifies (1) the simulation direction (either
2116forward or backward), (2) the start and (3) the end time of the simulation, (4)
2117the frequency $T_c$ of the model output, (5) the averaging time $\Delta T_c$ of
2118model output, and (6) the intervals $\Delta T_s$ at which concentrations are
2119sampled, (7) the time constant for particle splitting $\Delta t_s$, (8) the
2120synchronisation interval of FLEXPART, (9) the factor $c_{tl}$ by which the time
2121steps must be smaller than the Lagrangian time scale, and (10) the refinement
2122factor for the time step used for solving the Langevin equation of the vertical
2123component of the turbulent wind.  If (9) ($c_{tl}$) is negative, the Langevin
2124equations are solved with constant time steps according to the synchronisation
2125interval.  In that case, the value of (10) is arbitrary.  The synchronisation
2126interval is the minimum time interval used by the model for all activities
2127(such as concentration calculations, wet deposition calculations, interpolation
2128of data, mesoscale wind fluctuations or output of data) other than the
2129simulation of turbulent transport and dry deposition (if ($c_{tl}>0$).  Further
2130switches determine (11) whether concentrations, mixing ratios, residence times
2131or plume trajectories (or combinations thereof) are to be calculated, (12) the
2132option of particle position dump either at the end of or continuously during
2133the simulation, (13) on/off of subgrid terrain effect parameterization, (14)
2134on/off of deep convection parameterization, (15) on/off calculation of age
2135spectra, (16) continuation of simulation from previous particle dump, (17)
2136write output for each RELEASE location on/off, (18) on/off for mass flux
2137calculations and output, (19) on/off for the domain-filling option of FLEXPART,
2138(20) an indicator that determines whether mass or mass mixing ratio units are
2139to be used at the source, (21) an indicator that determines whether mass or
2140mass mixing ratio units are to be used at the receptor, (22) on/off of
2141additional compact dump of the positions of numbered particles, (23) on/off for
2142the use of nested output fields, (24) write sensitivity to initial conditions
2143in backward mode off/mass units/mass mixing ratio units.
2144
2145Two versions of \verb|COMMAND| may be used, which both can be read in by
2146FLEXPART: the first contains formatted input (i.e., a mask to be filled for the
2147various input options that must be filled in), the second contains largely
2148unformatted input and is recommended for the more experienced FLEXPART user.
2149The following example is for formatted input.
2150
2151\begin{scriptsize}\begin{verbatim}
2152********************************************************************************
2153*                                                                              *
2154*      Input file for the Lagrangian particle dispersion model FLEXPART        *
2155*                           Please select your options                         *
2156*                                                                              *
2157********************************************************************************
2158
21591. __                3X, I2
2160    1       
2161   LDIRECT           1 FOR FORWARD SIMULATION, -1 FOR BACKWARD SIMULATION
2162
21632. ________ ______   3X, I8, 1X, I6
2164   20040626 000000
2165   YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION
2166
21673. ________ ______   3X, I8, 1X, I6
2168   20040816 120000
2169   YYYYMMDD HHMISS   ENDING DATE OF SIMULATION
2170
21714. _____             3X, I5
2172    7200 
2173   SSSSS             OUTPUT EVERY SSSSS SECONDS
2174
21755. _____             3X, I5
2176    7200 
2177   SSSSS             TIME AVERAGE OF OUTPUT (IN SSSSS SECONDS)
2178
21796. _____             3X, I5
2180     900 
2181   SSSSS             SAMPLING RATE OF OUTPUT (IN SSSSS SECONDS)
2182
21837. _________         3X, I9
2184   999999999
2185   SSSSSSSSS         TIME CONSTANT FOR PARTICLE SPLITTING (IN SECONDS)
2186
21878. _____             3X, I5
2188     900 
2189   SSSSS             SYNCHRONISATION INTERVAL OF FLEXPART (IN SECONDS)
2190
21919---.--           4X, F6.4
2192     -5.0
2193    CTL              FACTOR, BY WHICH TIME STEP MUST BE SMALLER THAN TL
2194
219510. ---              4X, I3
2196      4
2197    IFINE            DECREASE OF TIME STEP FOR VERTICAL MOTION BY FACTOR IFINE
2198
219911. -                4X, I1
2200    3 
2201    IOUT             1 CONCENTRATION (RESIDENCE TIME FOR BACKWARD RUNS) OUTPUT, 2 MIXING RATIO OUTPUT, 3 BOTH,4 PLUME TRAJECT., 5=1+4
2202
220312. -                4X, I1
2204    2 
2205    IPOUT            PARTICLE DUMP: 0 NO, 1 EVERY OUTPUT INTERVAL, 2 ONLY AT END
2206
220713. _                4X, I1
2208    1
2209    LSUBGRID         SUBGRID TERRAIN EFFECT PARAMETERIZATION: 1 YES, 0 NO
2210
221114. _                4X, I1
2212    1
2213    LCONVECTION      CONVECTION: 1 YES, 0 NO
2214
221515. _                4X, I1
2216    0
2217    LAGESPECTRA      AGE SPECTRA: 1 YES, 0 NO
2218
221916. _                4X, I1
2220    0
2221    IPIN             CONTINUE SIMULATION WITH DUMPED PARTICLE DATA: 1 YES, 0 NO
2222
222317. _
2224    0                4X,I1
2225    IOFR             IOUTPUTFOREACHREL CREATE AN OUTPUT FILE FOR EACH RELEASE LOCATION: 1 YES, 0 NO
2226
222718. _                4X, I1
2228    0
2229    IFLUX            CALCULATE FLUXES: 1 YES, 0 NO
2230
223119. _                4X, I1
2232    0
2233    MDOMAINFILL      DOMAIN-FILLING TRAJECTORY OPTION: 1 YES, 0 NO, 2 STRAT. O3 TRACER
2234
223520. _                4X, I1
2236    1
2237    IND_SOURCE       1=MASS UNIT , 2=MASS MIXING RATIO UNIT
2238
223921. _                4X, I1
2240    1
2241    IND_RECEPTOR     1=MASS UNIT , 2=MASS MIXING RATIO UNIT
2242
224322. _                4X, I1
2244    0
2245    MQUASILAG        QUASILAGRANGIAN MODE TO TRACK INDIVIDUAL PARTICLES: 1 YES, 0 NO
2246
224723. _                4X, I1
2248    0
2249    NESTED_OUTPUT    SHALL NESTED OUTPUT BE USED? 1 YES, 0 NO
225024. _                4X, I1
2251    2
2252    LINIT_COND       INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
2253   
2254
22551. Simulation direction, 1 for forward, -1 for backward in time
2256        (consult Seibert and Frank, 2004 for backward runs)
2257
22582. Beginning date and time of simulation. Must be given in format
2259   YYYYMMDD HHMISS, where YYYY is YEAR, MM is MONTH, DD is DAY, HH is HOUR,
2260   MI is MINUTE and SS is SECOND. Current version utilizes UTC.
2261
22623. Ending date and time of simulation. Same format as 3.
2263
22644. Average concentrations are calculated every SSSSS seconds.
2265
22665. The average concentrations are time averages of SSSSS seconds
2267   duration. If SSSSS is 0, instantaneous concentrations are outputted.
2268
22696. The concentrations are sampled every SSSSS seconds to calculate the time
2270   average concentration. This period must be shorter than the averaging time.
2271
22727. Time constant for particle splitting. Particles are split into two
2273   after SSSSS seconds, 2xSSSSS seconds, 4xSSSSS seconds, and so on.
2274
22758. All processes are synchronized with this time interval (lsynctime).
2276   Therefore, all other time constants must be multiples of this value.
2277   Output interval and time average of output must be at least twice lsynctime.
2278
22799. CTL must be >1 for time steps shorter than the  Lagrangian time scale
2280   If CTL<0, a purely random walk simulation is done
2281
228210.IFINE=Reduction factor for time step used for vertical wind
2283
228411.IOUT determines how the output shall be made: concentration
2285   (ng/m3, Bq/m3), mixing ratio (pptv), or both, or plume trajectory mode,
2286   or concentration + plume trajectory mode.
2287   In plume trajectory mode, output is in the form of average trajectories.
2288
228912.IPOUT determines whether particle positions are outputted (in addition
2290   to the gridded concentrations or mixing ratios) or not.
2291   0=no output, 1 output every output interval, 2 only at end of the
2292   simulation
2293
229413.Switch on/off subgridscale terrain parameterization (increase of
2295   mixing heights due to subgridscale orographic variations)
2296
229714.Switch on/off the convection parameterization
2298
229915.Switch on/off the calculation of age spectra: if yes, the file AGECLASSES
2300   must be available
2301
230216. If IPIN=1, a file "partposit_end" from a previous run must be available in
2303    the output directory. Particle positions are read in and previous simulation
2304    is continued. If IPIN=0, no particles from a previous run are used
2305
230617. IF IOUTPUTFOREACHRELEASE is set to 1, one output field for each location
2307    in the RELEASES file is created. For backward calculation this should be
2308    set to 1. For forward calculation both possibilities are applicable.
2309
231018. If IFLUX is set to 1, fluxes of each species through each of the output
2311    boxes are calculated. Six fluxes, corresponding to northward, southward,
2312    eastward, westward, upward and downward are calculated for each grid cell of
2313    the output grid. The control surfaces are placed in the middle of each
2314    output grid cell. If IFLUX is set to 0, no fluxes are determined.
2315
231619. If MDOMAINFILL is set to 1, the first box specified in file RELEASES is used
2317    as the domain where domain-filling trajectory calculations are to be done.
2318    Particles are initialized uniformly distributed (according to the air mass
2319    distribution) in that domain at the beginning of the simulation, and are
2320    created at the boundaries throughout the simulation period.
2321
232220. IND_SOURCE switches between different units for concentrations at the source
2323    NOTE that in backward simulations the release of computational particles
2324    takes place at the "receptor" and the sampling of particles at the "source".
2325          1=mass units (for bwd-runs = concentration)
2326          2=mass mixing ratio units
232721. IND_RECEPTOR switches between different units for concentrations at the receptor
2328          1=mass units (concentrations)
2329          2=mass mixing ratio units
2330
233122. MQUASILAG indicates whether particles shall be numbered consecutively (1) or
2332    with their release location number (0). The first option allows tracking of
2333    individual particles using the partposit output files
2334
233523. NESTED_OUTPUT decides whether model output shall be made also for a nested
2336    output field (normally with higher resolution)
2337
233824. LINIT_COND determines whether, for backward runs only, the sensitivity to initial
2339    conditions shall be calculated and written to output files
2340    0=no output, 1 or 2 determines in which units the initial conditions are provided.
2341
2342\end{verbatim}\end{scriptsize}
2343
2344\newpage
2345
2346\subsubsection{File OUTGRID}
2347The file \verb|OUTGRID| specifies the output grid.
2348\begin{scriptsize}\begin{verbatim}
2349********************************************************************************
2350*                                                                              *
2351*      Input file for the Lagrangian particle dispersion model FLEXPART        *
2352*                       Please specify your output grid                        *
2353*                                                                              *
2354********************************************************************************
2355
23561------.----       4X,F11.4
2357       -10.0000       GEOGRAFICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
2358    OUTLONLEFT        (left boundary of the first grid cell - not its centre)
2359
2360
23612------.----       4X,F11.4
2362        40.0000       GEOGRAFICAL LATITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
2363    OUTLATLOWER       (lower boundary of the first grid cell - not its centre)
2364
23653-----             4X,I5
2366      101             NUMBER OF GRID POINTS IN X DIRECTION (= No. of cells + 1)
2367    NUMXGRID
2368
23694-----             4X,I5
2370       47             NUMBER OF GRID POINTS IN Y DIRECTION (= No. of cells + 1)
2371    NUMYGRID
2372
23735------.---        4X,F10.3
2374         0.500        GRID DISTANCE IN X DIRECTION
2375    DXOUTLON
2376
23776------.---        4X,F10.3
2378         0.500        GRID DISTANCE IN Y DIRECTION
2379    DYOUTLAT
2380
23817-----.-           4X, F7.1
2382      100.0
2383    LEVEL 1           HEIGHT OF LEVEL (UPPER BOUNDARY)
2384
23858-----.-           4X, F7.1
2386      300.0
2387    LEVEL 2           HEIGHT OF LEVEL (UPPER BOUNDARY)
2388
23899-----.-           4X, F7.1
2390      600.0
2391    LEVEL 3           HEIGHT OF LEVEL (UPPER BOUNDARY)
2392
239310. -----.-           4X, F7.1
2394     1000.0
2395    LEVEL 4           HEIGHT OF LEVEL (UPPER BOUNDARY)
2396
239711. -----.-           4X, F7.1
2398     2000.0
2399    LEVEL 5           HEIGHT OF LEVEL (UPPER BOUNDARY)
2400
240112. -----.-           4X, F7.1
2402     3000.0
2403    LEVEL 6           HEIGHT OF LEVEL (UPPER BOUNDARY)
2404\end{verbatim}\end{scriptsize}
2405In order to define the grid for a nested output field, the file
2406\verb|OUTGRID_NEST| must exist.  It has the same format as file \verb|OUTGRID|,
2407but does not contain the vertical level information:
2408
2409\begin{scriptsize}\begin{verbatim}
2410********************************************************************************
2411*                                                                              *
2412*      Input file for the Lagrangian particle dispersion model FLEXPART        *
2413*                       Please specify your output grid                        *
2414*                                                                              *
2415********************************************************************************
2416
24171------.----       4X,F11.4
2418      -125.0000       GEOGRAFICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
2419    OUTLONLEFT        (left boundary of the first grid cell - not its centre)
2420
24212------.----       4X,F11.4
2422        25.0000       GEOGRAFICAL LATITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
2423    OUTLATLOWER       (lower boundary of the first grid cell - not its centre)
2424
24253-----             4X,I5
2426        1             NUMBER OF GRID POINTS IN X DIRECTION (= No. of cells + 1)
2427    NUMXGRID
2428
24294-----             4X,I5
2430        1             NUMBER OF GRID POINTS IN Y DIRECTION (= No. of cells + 1)
2431    NUMYGRID
2432
24335------.-----      4X,F12.5
2434         0.33333      GRID DISTANCE IN X DIRECTION
2435    DXOUTLON
2436
24376------.-----      4X,F12.5
2438         0.25000      GRID DISTANCE IN Y DIRECTION
2439    DYOUTLAT
2440\end{verbatim}\end{scriptsize}
2441
2442\newpage
2443
2444\subsubsection{File RECEPTORS}
2445
2446\verb|RECEPTORS| specifies the receptor locations for which the parabolic
2447kernel method shall be applied to calculate air concentrations.  The maximum
2448number of receptor sites is set by parameter \verb|maxreceptor| in file
2449\verb|includepar|.
2450\begin{scriptsize}\begin{verbatim}
2451********************************************************************************
2452*                                                                              *
2453*      Input file for the Lagrangian particle dispersion model FLEXPART        * 
2454*                       Please specify your receptor points                    * 
2455*     For the receptor points, ground level concentrations are calculated      * 
2456*                                                                              *
2457********************************************************************************
24581----------------  4X,A16
2459    F15               NAME OF RECEPTOR POINT
2460    RECEPTORNAME
2461
24622------.----       4X,F11.4
2463         6.1333       GEOGRAFICAL LONGITUDE
2464    XRECEPTOR
2465
24663------.----       4X,F11.4
2467        49.0833       GEOGRAFICAL LATITUDE
2468    YRECEPTOR
2469================================================================================
24701----------------  4X,A16
2471    NL01              NAME OF RECEPTOR POINT
2472    RECEPTORNAME
2473
24742------.----       4X,F11.4
2475         5.7833       GEOGRAFICAL LONGITUDE
2476    XRECEPTOR
2477
24783------.----       4X,F11.4
2479        50.9167       GEOGRAFICAL LATITUDE
2480    YRECEPTOR
2481================================================================================
2482\end{verbatim}\end{scriptsize}
2483
2484\newpage
2485
2486\subsubsection{File RELEASES}
2487
2488\begin{footnotesize}\verb|RELEASES|\end{footnotesize} defines the release
2489specifications.  In the first input line, the number $N$ of emitted species is
2490defined (1 in the example below).  At all locations, the same species must be
2491released.  The next $N$ input lines give a cross-reference to the respective file
2492\begin{footnotesize}\verb|SPECIES_nnn|\end{footnotesize}, where the physical and
2493chemical properties of the released species are given (also the temporal
2494variations of emissions is defined for each species).  Then follows a list of
2495release sites, for each of which the release characteristics must be entered:
2496the beginning and the ending time of the release, geographical coordinates of
2497the lower left and upper right corners of the release location, type of
2498vertical coordinate (above ground level, or above sea level), lower level and
2499upper level of source box, the number of particles to be used, and the total
2500mass emitted.  Note that the mass entry must be repeated $N$ times, one mass
2501per species released.  Finally, a name is assigned to each release point.\par
2502
2503The particles are released from random locations within a four-dimensional box
2504extending from the lower to the upper level above a rectangle (on a lat/lon
2505grid) defined by the geographical coordinates, and between the release's start
2506and end.  With some of the coordinates set identically, line or point sources
2507can be specified.
2508
2509As for \verb|COMMAND|, the \verb|RELEASES| file can be provided formatted or
2510unformatted.  The example below shows the formatted version.
2511
2512\begin{scriptsize}\begin{verbatim}
2513*************************************************************************
2514*                                                                       *
2515*                                                                       *
2516*                                                                       *
2517*   Input file for the Lagrangian particle dispersion model FLEXPART    *
2518*                        Please select your options                     *
2519*                                                                       *
2520*                                                                       *
2521*                                                                       *
2522*************************************************************************
2523+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2524  1     
2525___                        i3    Total number of species emitted
2526
2527 24
2528___                        i3    Index of species in file SPECIES
2529
2530=========================================================================
253120011028  150007
2532________ ______            i8,1x,i6 Beginning date and time of release
2533
253420011028  150046
2535________ ______            i8,1x,i6 Ending date and time of release
2536
2537   9.4048 
2538____.____                  f9.4  Longitude [DEG] of lower left corner
2539
2540  48.5060
2541____.____                  f9.4  Latitude [DEG] of lower left corner
2542
2543   9.5067
2544____.____                  f9.4  Longitude [DEG] of upper right corner
2545
2546  48.5158
2547____.____                  f9.4  Latitude [DEG] of upper right corner
2548
2549        2
2550_________                  i9    1 for m above ground, 2 for m above sea level, 3 for pressure in hPa
2551
2552 6933.60
2553_____.___                  f10.3 Lower z-level (in m agl or m asl)
2554 
2555 6950.40
2556_____.___                  f10.3 Upper z-level (in m agl or m asl)
2557 
2558    20000               
2559_________                  i9    Total number of particles to be released
2560
25611.0000E00
2562_.____E__                  e9.4  Total mass emitted
2563
2564FLIGHT_11242
2565________________________________________   character*40 comment
2566+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
256720011028  150047
2568________ ______            i8,1x,i6 Beginning date and time of release
2569
257020011028  150107
2571________ ______            i8,1x,i6 Ending date and time of release
2572
2573   9.3038 
2574____.____                  f9.4  Longitude [DEG] of lower left corner
2575
2576  48.5158
2577____.____                  f9.4  Latitude [DEG] of lower left corner
2578
2579   9.4048
2580____.____                  f9.4  Longitude [DEG] of upper right corner
2581
2582  48.5906
2583____.____                  f9.4  Latitude [DEG] of upper right corner
2584
2585        2
2586_________                  i9    1 for m above ground, 2 for m above sea level, 3 for pressure in hPa
2587
2588 6833.50
2589_____.___                  f10.3 Lower z-level (in m agl or m asl)
2590 
2591 6950.40
2592_____.___                  f10.3 Upper z-level (in m agl or m asl)
2593 
2594    20000               
2595_________                  i9    Total number of particles to be released
2596
25971.0000E00
2598_.____E__                  e9.4  Total mass emitted
2599
2600FLIGHT_11185
2601________________________________________   character*40 comment
2602+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2603
2604\end{verbatim}\end{scriptsize}
2605
2606\newpage
2607
2608\subsubsection{File AGECLASSES}
2609
2610\verb|AGECLASSES| provides the times for the age class calculation.  In the
2611first data line, the number $n$ of age classes is set, and ages are listed in
2612the following $n$ lines.  The entries specify the end times (in seconds) of the
2613respective intervals to be used, the first one starting at zero seconds.
2614Particles are dropped from the simulation once they exceed the maximum age.
2615Even if no age classes are needed, this option (with the number of age classes
2616set to 1) can be useful to determine the age at which particles are removed
2617from the simulation.
2618
2619\begin{footnotesize}\begin{verbatim}
2620************************************************
2621*                                              *
2622*Lagrangian particle dispersion model FLEXPART *
2623*         Please select your options           *
2624*                                              *
2625*This file determines the ageclasses to be used*
2626*                                              *
2627*Ages are given in seconds. The first class    *
2628*starts at age zero and goes up to the first   *
2629*age specified. The last age gives the maximum *
2630*time a particle is carried in the simulation. *
2631*                                              *
2632************************************************
2633 6          Integer        Number of age classes
263443200       Integer        Age class 1
263586400       Integer        Age class 2
2636129600
2637172800
2638259200
2639345600
2640\end{verbatim}\end{footnotesize}
2641
2642\newpage
2643
2644\subsubsection{Files SPECIES\_nnn and surfdata.t}
2645
2646\begin{footnotesize}\verb|SPECIES_nnn|\end{footnotesize} (where \verb|nnn| is
2647a zero-padded identifier of the species number) specifies all physico-chemical
2648properties for the given species.  Entries are the half life (due to
2649radioactive or chemical decay), wet deposition information (\verb|A| and \verb|B| are the
2650factors defined by Eq.~\ref{wetscav}), dry deposition information for
2651gases ($D=D_{H_2O}/D_i$, $D_{H_2O}$ is the diffusivity of water vapor and \verb|D_i|
2652is the diffusivity of the species, \verb|H| is the effective Henry's constant, and
2653\verb|f0| varies between 0 and 1 and gives the reactivity of a species relative to
2654that of ozone.  For nonreactive species \verb|f0| is 0, for slightly reactive it is
26550.1 and for highly reactive it is 1.), dry deposition information for
2656particulates (\verb|rho| specifies the density of the substance, \verb|dquer| its mean
2657diameter $\overline{d_p}$, and \verb|dsig| the measure of variation $\sigma_p$).
2658Radioactive decay is switched off by specifying a negative half life, wet
2659deposition is switched off by specifying negative \verb|A|, dry deposition of gases
2660is switched off by negative \verb|D|, dry deposition of particles is switched off by
2661negative \verb|rho|.  If no detailed information for deposition velocity calculation
2662is available, a constant deposition velocity \verb|vd| (cm s$^{-1}$) can be used.
2663\verb|molweight| gives the molecular weight of the species, which is needed for
2664mixing ratio output.  For degradation in an monthly averaged OH field
2665\citep{bey2001} the OH Reaction ratio at 25\degreee C for the species can be given,
2666units are (cm$^3$/s).
2667
2668Optionally an emission variation information can be added at the end of the
2669file, defined as following: Since FLEXPART version 6.0, emission factors can be
2670defined that change the temporal variation of particle releases.  This is
2671useful, for instance, to simulate the typical daily and weekly cycle of
2672anthropogenic emissions.  The emission factors are given in the file of the
2673corresponding species \verb|SPECIES_nnn|, where \verb|nnn| is the species
2674number defined in file \verb|RELEASES|.  If no emission variation information
2675is given, emission rates for species \verb|nnn| are taken as constant.  Release
2676rates can vary with the hour of the day and with the day of the week, according
2677to the local time at the release location.  Emission factors must be 1 on
2678average.  24 hourly as well as 7 daily values must be specified.  Furthermore,
2679different disaggregation factors must be given for area sources and for point
2680sources.  FLEXPART distinguishes between the two using the lower altitude of
2681the release box: area sources are assumed to start below 0.5 m above the
2682ground, whereas point sources are assumed to be higher.  Please note that when
2683this option is used, it is not so easy to determine the maximum number of
2684particles present at a particular time of the model run.  It might then be
2685necessary to increase the parameter \verb|maxpart| to a higher value than what
2686would otherwise be needed.  The following is an example for an
2687\verb|SPECIES_nnn| file including emission variation.
2688
2689\begin{scriptsize}\begin{verbatim}
2690****************************************************************************
2691*                                                                          *
2692*     Input file for the Lagrangian particle dispersion model FLEXPART     *
2693*               Definition file of chemical species/radionuclides          *
2694*                                                                          *
2695****************************************************************************
2696EU-CO              Tracer name
2697-999.9             Species half life
2698-9.9E-09           Wet deposition - A
2699                   Wet deposition - B
2700-9.9               Dry deposition (gases) - D
2701                   Dry deposition (gases) - Henrys const.
2702                   Dry deposition (gases) - f0 (reactivity)
2703-9.9E09            Dry deposition (particles) - rho
2704                   Dry deposition (particles) - dquer
2705                   Dry deposition (particles) - dsig
2706-9.99              Alternative: dry deposition velocity
2707 28.00             molweight
2708-9.9E-09           OH Reaction rate at 25 deg, [cm^3/s]
2709-9                 number of associated specias (neg. none) - not implemented yet
2710-99.99             KOA - organic matter air partitioning
2711hr_start co_area  co_point
2712 0       0.535       0.932      0-1 local time
2713 1       0.405       0.931      1-2 local time
2714 2       0.317       0.927
2715 3       0.265       0.926
2716 4       0.259       0.928
2717 5       0.367       0.936
2718 6       0.668       0.952
2719 7       1.039       0.975
2720 8       1.015       1.046
2721 9       0.965       1.055
272210       1.016       1.061
272311       1.133       1.064
272412       1.269       1.067
272513       1.368       1.068
272614       1.516       1.069
272715       1.681       1.068
272816       1.777       1.024
272917       1.827       1.017
273018       1.538       1.008
273119       1.282       1.007
273220       1.136       1.004
273321       1.020       0.996
273422       0.879       0.981
273523       0.723       0.958      23-24 local time
2736week_day co_area  co_point
27371        1.060       1.000      Monday
27382        1.060       1.000      Tuesday
27393        1.060       1.000      Wednesday
27404        1.060       1.000      Thursday
27415        1.060       1.000      Friday
27426        0.900       1.000      Saturday
27437        0.800       1.000      Sunday
2744\end{verbatim}\end{scriptsize}
2745
2746\begin{footnotesize}\verb|IGBP_int1.dat|\end{footnotesize} contains the landuse inventory in binary format, and
2747\begin{footnotesize}\verb|surfdata.t|\end{footnotesize}, shown below, gives the roughness lengths for each landuse class:
2748
2749\begin{footnotesize}
2750\begin{verbatim}
275113 landuse categories are related roughness length
2752--------------------------------------------------------
2753landuse   comment                               z0
2754--------------------------------------------------------
2755 1      Urban land                                   0.7       
2756 2      Agricultural land                            0.1
2757 3      Range land                                   0.1
2758 4      Deciduous forest                             1.
2759 5      Coniferous forest                            1.
2760 6      Mixed forest including wetland               0.7       
2761 7      Water, both salt and fresh                   0.001
2762 8      Barren land mostly desert                    0.01
2763 9      Nonforested wetland                          0.1
276410      Mixed agricultural and range land            0.1
276511      Rocky open areas with low growing shrubs     0.05
276612 Snow and ice                                 0.001
276713 Rainforest                                   1.
2768\end{verbatim}
2769\end{footnotesize}
2770
2771\newpage
2772
2773\subsubsection{File surfdepo.t}
2774\verb|surfdepo.t| gives the resistances needed for the parameterization of dry deposition of gases for the 13 landuse classes and five seasonal categories.
2775This file must not be changed by the user.\par
2776 
2777\begin{scriptsize}\begin{verbatim}
2778DRY DEPOSITION
2779==============================================================================
2780AFTER WESELY, 1989
2781==============================================================================
27821 to 11: Landuse types after Wesely; 12 .. snow, 13 .. rainforest
2783==============================================================================
2784Values are tabulated for 5 seasonal categories:
27851     Midsummer with lush vegetation
27862     Autumn with unharvested cropland
27873     Late autumn after frost, no snow
27884     Winter, snow on ground and subfreezing
27895     Transitional spring with partially green short annuals
2790==============================================================================
2791               1       2       3       4       5       6       7       8       9     10       11      12      13 
2792________________________________________________________________________________________________________________
2793ri         9999.     60.    120.     70.    130.    100.   9999.   9999.     80.    100.    150.    9999.    2001
2794rlu        9999.   2000.   2000.   2000.   2000.   2000.   9999.   9999.   2500.   2000.   4000.    9999.   1000.
2795rac         100.    200.    100.   2000.   2000.   2000.      0.      0.    300.    150.    200.       0.   2000.
2796rgss        400.    150.    350.    500.    500.    100.      0.   1000.      0.    220.    400.     100.    200.
2797rgso        300.    150.    200.    200.    200.    300.   2000.    400.   1000.    180.    200.   10000.    200.
2798rcls       9999.   2000.   2000.   2000.   2000.   2000.   9999.   9999.   2500.   2000.   4000.    9999.   9999.
2799rclo       9999.   1000.   1000.   1000.   1000.   1000.   9999.   9999.   1000.   1000.   1000.    9999.   9999.
2800_________________________________________________________________________________________________________________
2801ri         9999.   9999.   9999.   9999.    250.    500.   9999.   9999.   9999.   9999.   9999.    9999.    200. 2
2802rlu        9999.   9000.   9000.   9000.   4000.   8000.   9999.   9999.   9000.   9000.   9000.    9999.   1000.
2803rac         100.    150.    100.   1500.   2000.   1700.      0.      0.    200.    120.    140.       0.   2000.
2804rgss        400.    200.    350.    500.    500.    100.      0.   1000.      0.    300.    400.     100.    200.
2805rgso        300.    150.    200.    200.    200.    300.   2000.    400.    800.    180.    200.   10000.    200.
2806rcls       9999.   9000.   9000.   9000.   2000.   4000.   9999.   9999.   9000.   9000.   9000.    9999.   9999.
2807rclo       9999.    400.    400.    400.   1000.    600.   9999.   9999.    400.    400.    400.    9999.   9999.
2808_________________________________________________________________________________________________________________
2809ri         9999.   9999.   9999.   9999.    250.    500.   9999.   9999.   9999.   9999.   9999.    9999.    200. 3
2810rlu        9999.   9999.   9000.   9000.   4000.   8000.   9999.   9999.   9000.   9000.   9000.    9999.   1000.
2811rac         100.     10.    100.   1000.   2000.   1500.      0.      0.    100.     50.    120.       0.   2000.
2812rgss        400.    150.    350.    500.    500.    200.      0.   1000.      0.    200.    400.     100.    200.
2813rgso        300.    150.    200.    200.    200.    300.   2000.    400.   1000.    180.    200.   10000.    200.
2814rcls       9999.   9999.   9000.   9000.   3000.   6000.   9999.   9999.   9000.   9000.   9000.    9999.   9999.
2815rclo       9999.   1000.    400.    400.   1000.    600.   9999.   9999.    800.    600.    600.    9999.   9999.
2816_________________________________________________________________________________________________________________
2817ri         9999.   9999.   9999.   9999.    400.    800.   9999.   9999.   9999.   9999.   9999.    9999.    200. 4
2818rlu        9999.   9999.   9999.   9999.   6000.   9000.   9999.   9999.   9000.   9000.   9000.    9999.   1000.
2819rac         100.     10.     10.   1000.   2000.   1500.      0.      0.     50.     10.     50.       0.   2000.
2820rgss        100.    100.    100.    100.    100.    100.      0.   1000.    100.    100.     50.     100.    200.
2821rgso        600.   3500.   3500.   3500.   3500.   3500.   2000.    400.   3500.   3500.   3500.   10000.    200.
2822rcls       9999.   9999.   9999.   9000.    200.    400.   9999.   9999.   9000.   9999.   9000.    9999.   9999.
2823rclo       9999.   1000.   1000.    400.   1500.    600.   9999.   9999.    800.   1000.    800.    9999.   9999.
2824_________________________________________________________________________________________________________________
2825ri         9999.    120.    240.    140.    250.    190.   9999.   9999.    160.    200.    300.    9999.    2005
2826rlu        9999.   4000.   4000.   4000.   2000.   3000.   9999.   9999.   4000.   4000.   8000.    9999.   1000.
2827rac         100.     50.     80.   1200.   2000.   1500.      0.      0.    200.     60.    120.       0.   2000.
2828rgss        500.    150.    350.    500     500.    200.      0.   1000.      0.    250.    400.     100.    200.
2829rgso        300.    150.    200.    200.    200.    300.   2000.    400.   1000.    180.    200.   10000.    200.
2830rcls       9999.   4000.   4000.   4000.   2000.   3000.   9999.   9999.   4000.   4000.   8000.    9999.   9999.
2831rclo       9999.   1000.    500.    500.   1500.    700.   9999.   9999.    600.    800.    800.    9999.   9999.
2832_________________________________________________________________________________________________________________
2833\end{verbatim}\end{scriptsize}
2834
2835\end{document}
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG