The Particle Dispersion Model
Version 3.1
User Guide

ANDREAS STOHL
Lehrstuhl für Bioklimatologie und Immissionsforschung,
University of Munich
Am Hochanger 13, 85354 Freising, Germany

Contents

Acknowledgments: The new subgrid topography scheme has been developed by Gerhard Wotawa and Petra Seibert, and the new sparse matrix output format has been suggested by Petra Seibert. Also, discussions on other subjects with Gerhard and Petra continue to be valuable for the further development of FLEXPART.

Chapter 1
Introduction

Particle modeling is the Lagrangian approach towards the simulation of atmospheric dispersion. Lagrangian particle models compute trajectories of a large number of individual so-called particles (not necessarily representing real particles, but small air parcels) to describe the transport and diffusion of substances in the atmosphere. Usually they are used to simulate the dispersion of air pollutants, released for instance due to an accident in a nuclear power plant.

A main advantage of Lagrangian models is that there is no artificial numerical diffusion as occurring in Eulerian models. This is of special importance in the vicinity of the source of an air pollutant, where in Eulerian models the pollutant is instantaneously mixed over at least one grid box, which can cause large subsequent transport errors. The description of turbulence with particle models is of similar quality as large eddy simulations, and it is more accurate than with Eulerian models. Particle models are independent of a computational grid and, in principle, provide infinitesimally small resolution. Practically, however, their effective resolution is determined by the number of particles released and by the resolution of meteorological input data.

Tens or hundreds of thousand particles may be required to accurately describe a multiday transport episode. This posts high demands on computer resources and explains why particle models have only become popular in recent years since powerful computers became available. The basis for most of nowadays particle models was developed in a paper by Thomson (1987). A very thorough monograph on the theoretical background of stochastic Lagrangian models has been published by Rodean (1996). Another interesting paper is the one by Wilson and Sawford (1996). Reviews of particle modeling, which also describe practical aspects, were provided by Zannetti (1992) and Uliasz (1994).

is a newly developed Lagrangian particle dispersion model designed for emergency response and research applications. It simulates the long-range transport, diffusion, dry and wet deposition, and radioactive decay of air pollutants released from point, line, area or volume sources. The management of input data was largely taken from FLEXTRA, a kinematic trajectory model (Stohl et al., 1995) developed at the Institute of Meteorology and Physics of the University of Agricultural Sciences in Vienna. The basic transport and diffusion code of the model was developed during the authors military service at the NBC school of the Austrian Forces (version 1.x). The deposition code was developed later on as a contract work for the Central Institute for Meteorology and Geodynamics in Vienna (version 2.x). Recently, the code underwent a major revision (current version 3.0), which optimized the runtime performance, especially when using short time steps. Some minor bugs were found and removed, the numerics of the model has been improved, and a density correction was introduced.

is based on model level data of the T213 L31 numerical weather prediction model of the European Centre for Medium-Range Weather Forecasts (ECMWF). Gridded data of this model are available at a maximum resolution of 0.5 ^° longitude and latitude, 31 vertical levels, and with a frequency of 3 hours. Usually, however, data with lower horizontal and temporal resolution are used.

may be used in forward mode to simulate the dispersion of pollutants released at a location, or in backward mode to determine the ``influence area'' of a given receptor point. Forward trajectories are calculated in the forward mode, backward trajectories in the backward mode. Deposition processes are only simulated in forward mode.

was evaluated using data from practically all large-scale tracer experiments available, namely from the Cross-Appalachian Tracer Experiment (CAPTEX), the Across North America Tracer Experiment (ANATEX) and the European Tracer Experiment (ETEX), comprising a total of 40 usable releases (Stohl et al., 1998). The model performed rather well compared to other models. The new version 3.1 was again validated with the CAPTEX and ETEX data. The performance was similar to previous versions.

Chapter 2
New features

1. Global datasets allowed: Use of cyclic boundary conditions; polar stereographic projection near the poles.
2. Nesting: An unlimited number of nests at various nesting levels may be used.
3. Uncertainty output: For every grid cell, FLEXPART calculates the uncertainty caused by the limited number of particles being used for the simulation. Summary uncertainty statistics are continuously created during the model run, so that one can re-start the model with more particles at an early stage, if necessary.
4. Sparse matrix format output: For every output field, it is checked whether it is more efficient to save the field in sparse matrix format, which only contains the concentrations in grid boxes with non-zero concentrations, rather than to dump the whole field. This saves a lot of mass storage space, especially for rather short model runs (small tracer clouds).
5. Alternative output of mixing ratio: one can choose whether the output shall be in concentration (ng m^-3, Bq m^-3) or in mixing ratio (pptv) units or in both.
6. Additional optional direct output of particle positions, which may be useful for visualization. Also, particles may be dumped at the end of the simulation, and may be used for re-initialization of a subsequent model run.
7. Tested for inverse modeling: FLEXPART 3.1 now definitely also works backward in time, which is useful for inverse modeling purposes. A first attempt at the reconstruction of the ETEX source term currently undertaken by Seibert and Stohl shows promising results.
8. Subgrid variability of mixing heights is treated somewhat differently than before.

Chapter 3
Input data and vertical model structure

expects the input data on ECMWF model (i.e. h) levels which are defined by a hybrid coordinate system. The conversion from h to pressure coordinates is given by p_k = A_k+B_kp_s and the heights of the h surfaces are defined by h_k = A_k/p₀+B_k, where h_k is the value of h at the k^_th model level, p_s is the surface pressure and p₀ is a pressure constant (101325 Pa). A_k and B_k are coefficients, chosen in such a way that the levels closest to the ground follow the topography, while the highest levels coincide with pressure surfaces. Intermediate levels show a gradual transition between topography-following and pressure levels. The vertical wind in hybrid coordinates is not directly available from the ECMWF archives, but is calculated by the pre-processor that extracts the data. A surface level is defined in addition to the regular ECMWF model h levels. 2 m temperature, 10 m winds and specific humidity from the first regular model level are attributed to this level. The inaccuracy introduced by the non-vanishing wind speed at the surface is marginal in most situations.

As mixing heights and parameterized random velocities (see sections 4 and 5) are given in meters and meters per second, respectively, the ECMWF model levels are transformed to terrain-following Cartesian z coordinates to avoid multiple coordinate transformations during the trajectory computations. The heights of the z levels equal the heights of the h levels at the lower left grid point of the computational domain and at the initial time. Thus, the vertical resolution is identical to that in the h coordinate system and the levels are close to the original levels in order to minimize interpolation errors. Interpolation from the h levels to the z levels is done with linear interpolation.

Surface parameters are needed as two-dimensional fields: surface pressure, total cloud cover, 10 m horizontal wind components, 2 m temperature and dew point temperature, large scale and convective precipitation, sensible heat flux, solar radiation, and east/west and north/south surface stress. If surface stress and heat flux are not available, friction velocity and surface sensible heat flux are computed from 10 m wind, 2 m temperature and their respective values at the first model level using the profile method (Berkowicz and Prahm, 1982). All other parameters must be available.

Topography, land-sea-mask and subgrid standard deviation of topography must be also available to the model. The former two may be included in the meteorological GRIB files, the latter, at the moment, must be provided in an additional file. If dry deposition of gases is to be calculated, a landuse inventory should be provided. For Europe, the inventory of Velde et al. (1994) is used.

Nesting is facilitated in FLEXPART 3.1. The input files for each nesting level must be stored in a different directory, and each directory must contain an AVAILABLE file, pointing towards the filenames. ECMWF input files must be available for ALL nests for ALL times, i.e., the AVAILABLE files of the nested domains must be consistent to the AVAILABLE file of the mother domain. The horizontal field dimensions of all nests can differ from each other, but the vertical discretization must be identical (i.e., no nesting in the vertical direction). However, the memory occupied by FLEXPART is determined by the nest with the largest number of grid points, since all nests are set to the same maximum field size (both in x and y direction) which must accomodate the largest nest. Therefore, it is most economical to have nests with comparable dimensions. Only the mother domain is initialized differently, and thus, is not subject to the above rule.

Chapter 4
Physical parameterization of boundary layer parameters

If deaccumulated surface stresses and surface sensible heat fluxes are not available, the profile method after Berkowicz and Prahm (1982) is applied to wind and temperature data provided at the first model level and at 10 m (for wind) and 2 m (for temperature). This is an iterative procedure. For its initialization, Obukhov length L is set to a large value. The following three equations are solved in five iterations which is sufficient to achieve convergence:

The heat flux is then computed by

Whenever possible, deaccumulated flux data from ECMWF should be provided since they are consistent with the ECMWF model calculations. Wotawa and Stohl (1997) have shown that friction velocities and heat fluxes calculated using the profile method are much less accurate.

Planetary boundary layer (PBL) heights are calculated according to Vogelezang and Holtslag (1996) using the critical Richardson number concept. The PBL height h_mix is set to the height of the first model level l for which the Richardson number

Spatial and temporal variations of PBL heights on scales not resolved by the ECMWF model play an important role in determining the thickness of the layer over which the tracer material is distributed. To elucidate this, we first consider unresolved temporal variations. Given a typical evolution of a convective mixed layer, the PBL height reaches its maximum value (say 1500 m) in the afternoon (for instance at 1700 LST), before a much shallower stable PBL forms. Now, if the meteorological data are available only at 1200 LST and at 1800 LST and the PBL heights at those times are, say, 1200 m and 200 m, and some standard interpolation procedure is used to determine the PBL height at 1700 LST, the PBL height is significantly underestimated. In the above example, using linear interpolation, a value of 370 m would have been obtained instead of 1500 m. If a tracer release takes place shortly before the breakdown of the convective mixed layer, this leads to a serious overestimation of the surface concentrations (a factor of four in the above example). Even if the release occurs already in the morning hours, the thickness of the tracer cloud in the evening would be restricted to 1200 m in the above example, whereas the correct thickness would be 1500 m.

Similar arguments are valid for spatial variations of PBL heights. A major reason for spatial variability of PBL heights is complex topography. It increases mechanical mixing due to enhanced shear stress, induces gravity waves during stable conditions, and it provides elevated sources of sensible heat during daytime. All these complex effects can produce turbulence and can thereby locally increase the PBL heights above the elevated areas. In addition to topographical effects, significant variations of the PBL heights can also be caused by differences in landuse or soil wetness (Hubbe et al., 1997). If a tracer cloud travels over such a patchy surface, its thickness would be determined by the maximum PBL height experienced along its path rather than by the average PBL height, which is obtained at the grid points of a coarse resolution meteorological model.

The best solution to this problem obviously is to use meteorological data with higher space and time resolution than the data from ECMWF. Unfortunately, such data are not as readily available, and if they are, it is not a priori clear whether the quality of these data (especially of the wind fields) is as high as the quality of the ECMWF data. Another solution would be to use a simple prognostic boundary layer model as an interpolation tool.

However, for the current version, we did not follow one of these solutions. Instead, we used a somewhat arbitrary parameterization to avoid a significant bias in the tracer cloud thickness which would result in strong overestimations of the surface tracer concentrations. To account for spatial variations induced by topography, we calculate an ënvelope" PBL height

Chapter 5
Particle transport and diffusion

5.1  Trajectory calculations

In order to save CPU time, , contrary to FLEXTRA, uses the simple scheme

5.2  The Langevin equation

Gaussian turbulence is assumed in , which is strictly valid only during stable and neutral conditions. During convective condtions, when turbulence is skewed (i.e. larger areas are occupied by downdrafts than by updrafts), this assumption is violated, but for long-range transport the error due to this is not so large. Cross-correlations between the different wind components are not taken into account, since Uliasz (1994) has shown that they have little effect for long-range dispersion.

Although density decreases with height, most particle models assume constant density flows. Especially for deep boundary layers, this is not justified. Stohl and Thomson (1998) have developed a density correction for Gaussian turbulence that takes the density effect into account.

With the above assumptions, the Langevin equation for the vertical wind component can be written as

Alternatively, the Langevin equation can be re-expressed in terms of w/s_w instead of w (Wilson et al., 1983):

For the discrete time step implementation of the above Langevin equations (at the example of equation 5.5), two different methods are used. When (Dt/t_Lw ) < 0.5,

If a particle reaches the surface or the top of the PBL layer, it is reflected according to Wilson and Flesch (1993), and the sign of the turbulent velocity is changed. If a particle reaches the top or the lateral boundary of the model domain, the computation of its trajectory is terminated.

5.3   Determination of the time steps

If the transport in the PBL shall be described more accurately, the time steps must be limited by the Lagrangian time scales. Since the vertical wind component is the most important one, only t_Lw is used to limit the time step. The user must specify two constants, c_tl and ifine. The first one determines the time step Dt_i according to

For solving the Langevin equation for the vertical wind component, a shorter time step Dt_w = Dt_i / ifine is used. This strategy (given sufficiently large values for c_tl and ifine) ensures that the particles stay vertically well-mixed also in very inhomogeneous turbulence, while keeping the computational cost at a minimum (i.e. only solution of the Langevin equation for w at the shortest time steps Dt_w).

5.4  The parameterization of the wind fluctuations

Neutral conditions:

Stable conditions:

Along-wind and cross-wind components are transformed to Cartesian wind components. In the free troposphere (z > h), s_u = s_v = 0.3 m/s and s_w = 0 m/s are used and Lagrangian time scales are set to 1 hour. The minimum t_Lu, t_Lv and t_Lw used are 10 seconds, 10 seconds and 30 seconds, respectively, to avoid excessive computation times for particles close to the surface.

5.5  Mesoscale velocity fluctuations

5.6  Particle splitting

Chapter 6
Radioactive decay

Chapter 7
Wet deposition

The scavenging coefficients L_i are species-dependent and increase with precipitation rate according to

Precipitation can vary considerably over the area covered by a single grid cell. As wet deposition depends nonlinearly on precipitation rate, this subgrid variability must be accounted for (Hertel et al., 1995). The area fraction which experiences precipitation given a certain grid-scale precipitation rate is calculated by

To save computing time, wet deposition is not calculated at each integration time step, but only once during each concentration sampling interval. With sampling intervals being shorter than the intervals between available meteorological fields, the inaccuracy introduced by this simplification is marginal.

Chapter 8
Dry deposition

8.1  Dry deposition of gases

The aerodynamic resistance r_a can be expressed with the flux-profile relationship which is based on Monin-Obukhov similarity theory (see, e.g., Stull, 1988)

Following Erisman et al. (1994), the quasilaminar sublayer resistance is

The surface resistance r_c is given by (Wesely, 1989)

The required data on landuse is taken from a European database developed by van de Velde et al. (1994). It provides the area fractions of eight landuse classes on a grid with 10' resolution. The roughness lengths z₀ needed for parameterization of r_a are also estimated from this landuse inventory. Table 8.1 shows the landuse classes and associated values of z₀. Charnock's relationship (see Stull, 1988) is used to calculate z₀ for the classes ``Ocean'' and ``Inland water'', because of its dependence on wave height

Deposition velocities in a grid cell are calculated for all landuse classes which cover a non-vanishing area of the cell. The average deposition velocity in the cell is then computed by weighting the individual deposition velocities with the respective area fractions of the landuse classes.

Outside the area of the European landuse database, the landuse classes are determined only from the ECMWF land-sea-mask, attributing the landuse classes ``Ocean'' to the sea surfaces and ``Grasslands'' to the land surfaces.

8.2  Dry deposition of particulate matter

The settling velocity is calculated from

Settling and dry deposition velocities are strongly dependent on the particulates sizes. Usually, their sizes cover some range (polydisperse aerosol) and it is not justified to describe particulate matter by a single average diameter. assumes a logarithmic normal distribution of the particulate mass over the range of particulate diameters. The user must specify the mean particulate diameter [`(d<sub>p</sub>)] and a measure of the variation around [`(d_p)], s_p. s_p = 0.1 or s_p = 10 both indicate that 68% of the mass are between 0.1[`(d<sub>p</sub>)] and 10[`(d_p)]. Then, the settling and deposition velocities are calculated for several particle diameters within the range of ±3 standard deviations of the logarithmic normal distribution. The particulate mass fractions represented by the diameter intervals are used to weight the respective velocities to yield average settling and deposition velocities. A value of s_p = 1 indicates monodisperse aerosol, i.e. aerosol with a uniform diameter.

Gravitational settling is important not only for the computation of the dry deposition velocity, but also affects the particles trajectory. For practical reasons, however, gravitational settling velocities are only taken into account for trajectory calculations when the particles represent a single species. If simulates the transport of multiple species, gravitational settling is neglected, because otherwise each species would follow its own trajectory. This would violate the model concept with several species being described by a single particle. Thus, if the dispersion of large particles with significant settling velocities is to be calculated, a single species calculation is suggested.

8.3  Loss of particle mass due to dry deposition

Chapter 9
Concentration and deposition calculation

The concentration of a pollutant at time t in a grid cell is calculated by sampling the mass fractions of all particles within the grid cell and dividing by the grid cell volume

Picture Omitted

The user can choose to output mixing ratios rather than concentrations. In this case, the density of the air is calculated from the data given on the mother domain and at some instant of time close to the actual time. The above approximations are sufficiently accurate.

The uncertainty of the gridded concentration output is estimated by carrying several classes of particles in the model simulation, and determining the concentration (or mixing ratio) at every grid point separately for every class. Then, the standard deviation of the individual classes is calculated, which is subsequently divided by the square root of the number of classes to yield the standard deviation of the mean of all classes. This number is reported in the output files. An overall uncertainty for the whole field is also written periodically to the standard output during the model run. This is reported as relative uncertainty (standard deviation divided by mean concentration).

In addition to the simple uniform kernel method, a computationally much more demanding parabolic kernel method as described in Uliasz (1994) can be used to calculate surface concentrations for a limited number of user-specified receptor points. It should be used for a few receptor sites only, since it can increase the computation time of FLEXPART considerably. The concentration is estimated from

The horizontal kernel bandwidths for each particle are defined according to

The vertical kernel bandwidth for each particle is given by

The kernel is defined as

With both the uniform and the parabolic kernel, concentrations C_Tc at time T_c are calculated as time-averages over the period [T_c-DT_c/2,T_C+DT_c/2]. DT_c must be specified by the user. To calculate the time-averages, concentrations C_Ts at times T_s within [T_c-DT_c/2,T_C+DT_c/2] are sampled at shorter intervals DT_s and are then divided by the number of samples taken:

Chapter 10
References

Asman, W.A.H. (1995): Parameterization of below-cloud scavenging of highly soluble gases under convective conditions. Atmos. Environ., 29, 1359-1368.

Baumann, K., and A. Stohl (1997): Validation of a long-range trajectory model using gas balloon tracks from the Gordon Bennett Cup 95. J. Appl. Meteor., 36, 711-720.

Businger, J.A., J.C. Wyngaard, Y. Izumi, and E.F. Bradley (1971): Flux-profile relationships in the atmospheric surface layer. J. Atmos. Sci., 28, 181-189.

De Baas, A.F., H. Van Dop, F.T.M. Nieuwstadt (1986): An application of the Langevin equation for inhomogeneous conditions to dispersion in a convective boundary layer. Q. J. R. Meteorol. Soc., 112, 165-180.

Beljaars, A.C.M., and A.A.M. Holtslag (1991): Flux parameterization over land surfaces for atmospheric models. J. Appl. Meteor., 30, 327-341.

Berkowicz, R., and L.P. Prahm (1982): Evaluation of the profile method for estimation of surface fluxes of momentum and heat. Atmos. Environ., 16, 2809-2819.

Businger, J.A., J.C. Wyngaard, Y. Izumi, and E.F. Bradley (1971): Flux-profile relationships in the atmospheric surface layer. J. Atmos. Sci., 28, 181-189.

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

Erisman, J.W., A. Van Pul, and P. Wyers (1994): Parametrization of surface resistance for the quantification of atmospheric deposition of acidifying pollutants and ozone. Atmos. Environ., 28, 2595-2607.

Gupta, S., McNider, R.T., Trainer, M., Zamora, R.J., Knupp, K. and Singh, M.P. (1997): Nocturnal wind structure and plume growth rates due to inertial oscillations. J. Appl. Meteor., 36, 1050-1063.

Hanna, S.R. (1982): Applications in air pollution modeling. In: Nieuwstadt F.T.M. and H. van Dop (ed.): Atmospheric Turbulence and Air Pollution Modelling. D. Reidel Publishing Company, Dordrecht, Holland.

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

Hubbe, J.M., Doran, J.C., Liljegren, J.C. and Shaw, W.J. (1997) Observations of spatial variations of boundary layer structure over the southern Great Plains cloud and radiation testbed. J. Appl. Meteor., 36, 1221-1231.

Kahl, J.D.W. (1996): On the prediction of trajectory model error. Atmos. Environ., 30, 2945-2957.

Legg, B.J., and M.R. Raupach (1982): Markov-chain simulation of particle dispersion in inhomogeneous flows: the mean drift velocity induced by a gradient in Eulerian velocity variance. Boundary-Layer Meteorol., 24, 3-13.

Maryon, R.H. (1998) Determining cross-wind variance for low frequency wind meander. Atmos. Environ., 32, 115-121.

McMahon, T.A., and P.J. Denison (1979): Empirical atmospheric deposition parameters - a survey. Atmos. Environ., 13, 571-585.

McNider, R.T., M.D. Moran, and R.A. Pielke (1988): Influence of diurnal and inertial boundary layer oscillations on long-range dispersion. Atmos. Environ., 22, 2445-2462.

Pruppacher, H.R., and J.D. Klett (1978): Microphysics of clouds and precipitation. D. Reidel Publishing Company, Dordrecht.

Rodean, H. (1996): Stochastic Lagrangian models of turbulent diffusion. Meteorological Monographs, 26 (48). American Meteorological Society, Boston, USA.

Ryall, D.B. and Maryon, R.H. (1997): Validation of the UK Met Office's NAME model against the ETEX dataset. In: Nodop, K. (editor): ETEX Symposium on Long-Range Atmospheric Transport, Model Verification and Emergency Response, European Commission EUR 17346, 151-154.

Slinn (1982): Predictions for particle deposition to vegetative canopies. Atmos. Environ., 16, 1785-1794.

Stohl, A. (1998) Computation, accuracy and applications of trajectories - a review and bibliography. Atmos. Environ.In press.

Stohl, A., Hittenberger, M., Wotawa, G. (1998): Validation of the Lagrangian particle dispersion model FLEXPART against large scale tracer experiment data. Atmos. Environ., in press.

Stohl, A., and P. Seibert (1998): Accuracy of trajectories as determined from the conservation of meteorological tracers. Q. J. R. Meteorol. Soc., in press.

Stohl, A. and Thomson, D.J. (1998): A density correction for Lagrangian particle dispersion models. Submitted to Boundary-Layer Meteorol.

Stohl, A., G. Wotawa, P. Seibert, and H. Kromp-Kolb (1995): Interpolation errors in wind fields as a function of spatial and temporal resolution and their impact on different types of kinematic trajectories. J. Appl. Meteor., 34, 2149-2165.

Stull, R.B. (1988): An Introduction to Boundary Layer Meteorology. Kluwer Academic Publishers, Dordrecht.

Thomson D.J. (1984): Random walk modelling of diffusion in inhomogeneous turbulence. Q. J. R. Meteorol. Soc., 110, 1107-1120.

Thomson D.J. (1987): Criteria for the selection of stochastic models of particle trajectories in turbulent flows. J. Fluid Mech., 180, 529-556.

Uliasz, M. (1994): Lagrangian particle dispersion modeling in mesoscale applications. In: Zannetti, P. (ed.): Environmental Modeling, Vol. II. Computational Mechanics Publications, Southampton, UK.

Velde van de, R.J., W.S. Faber, V.F. Katwijk van,, J.C.I. Kuylenstierna, H.J. Scholten., T.J.M. Thewessen, M. Verspuij, and M. Zevenbergen (1994): The Preparation of a European Land Use Database. National Institute of Public Health and Environmental Protection, Report nr 712401001, Bilthoven, The Netherlands.

Vogelezang, D.H.P., and A.A.M. Holtslag (1996): Evaluation and model impacts of alternative boundary-layer height formulations. Boundary-Layer Meteorol., in press.

Wesely, M.L., and B.B. Hicks (1977): Some factors that affect the deposition rates of sulfur dioxide and similar gases on vegetation. J. Air Poll. Contr. Assoc., 27, 1110-1116.

Wesely, M.L. (1989): Parameterization of surface resistances to gaseous dry deposition in regional-scale numerical models. Atmos. Environ., 23, 1293-1304.

Wilson, D.J., and Flesch, T.K. (1993): Flow boundaries in random-flight dispersion models: enforcing the well-mixed condition. J. Appl. Meteor., 32, 1695-1707.

Wilson, J.D., B.J. Legg, and D.J. Thomson (1983): Calculation of particle trajectories in the presence of a gradient in turbulent-velocity scale. Boundary-Layer Meteorol., 27, 163-169.

Wilson, J.D., and Sawford, B.L. (1996): Review of Lagrangian stochastic models for trajectories in the turbulent atmosphere. Boundary-Layer Meteorol., 78, 191-210.

Wilson, J.D., G.W. Thurtell, and G.E. Kidd (1981): Numerical simulation of particle trajectories in inhomogeneous turbulence, II: Systems with variable turbulent velocity scale. Boundary-Layer Meteorol., 21, 423-441.

Wotawa, G., A. Stohl, and H. Kromp-Kolb (1996): Parameterization of the planetary boundary layer over Europe - a data comparison between the observation based OML preprocessor and ECMWF model data. Contr. Atmos. Phys., 69, 273-284.

Wotawa, G. and Stohl, A. (1997) Boundary layer heights and surface fluxes of momentum and heat derived from ECMWF data for use in pollutant dispersion models - problems with data accuracy. In: Gryning, S.-E., Beyrich, F. and E. Batchvarova (editors): The Determination of the Mixing Height - Current Progress and Problems. EURASAP Workshop Proceedings, 1-3 October 1997, Risø National Laboratory, Denmark.

Zannetti, P. (1992): Particle Modeling and Its Application for Simulating Air Pollution Phenomena. In: Melli P. and P. Zannetti (ed.): Environmental Modelling. Computational Mechanics Publications, Southampton, UK.

Chapter 11
Input files of the FLEXPART model

11.1  The pathnames file

/home/atmos1/as/flextra3.0/options/
/home/atmos1/as/flextra3.0/output/
/daten2/ecfields/2.0/
/daten2/ecfields/2.0/AVAILABLE
/daten2/ecfields/1.0/
/daten2/ecfields/1.0/AVAILABLE
/daten2/ecfields/0.5/
/daten2/ecfields/0.5/AVAILABLE
============================================

Line 1: path where control files "COMMAND" and "STARTPOINTS" are available
Line 2: name of directory where output files are generated
Line 3: path where meteorological fields are available (mother grid)
Line 4: full filename of "AVAILABLE"-file (mother grid)
Subsequent lines:
Line 2n+3: path where meteorological fields are available (nested grid n)
Line 2n+4: full filename of "AVAILABLE"-file (nested grid n)

Line below last pathname must be:
============================================

The grids must be arranged such as that the coarse-scale nests
come before the fine-scale nests. Multiple nests of the same
nesting level are allowed. In that case, the order is arbitrary.

11.2  Files in directory windfields

AVAILABLE       EN94102300      EN94102303      EN94102306
EN94102309      EN94102312      EN94102315      LANDSEAMASK
OROGRAPHY       EXCESSORO

The names of the grib code-files and the times (in UTC) for which their data are valid must be indicated in the file AVAILABLE:

DATE     TIME         FILENAME     SPECIFICATIONS
YYYYMMDD HHMISS
________ ______      __________      __________
19941023 000000      EN94102300      ON DISC
19941023 030000      EN94102303      ON DISC
19941023 060000      EN94102306      ON DISC
19941023 090000      EN94102309      ON DISC
19941023 120000      EN94102312      ON DISC
19941023 150000      EN94102315      ON DISC

LANDSEAMASK, OROGRAPHY and SUBGRIDORO contain the ECMWF land-sea-mask, the topography and the subgrid topography standard deviation of the model domain. All three files are in ASCII. Only the file SUBGRIDORO is necessary when the topography and land-sea-mask information is contained in the GRIB files. The subgrid topography is, at the moment, read from the ASCII file in any case.

11.3  Files in directory options

COMMAND         RELEASES        surfdata.t
OUTGRID         SPECIES         surfdepo.t
RECEPTORS       landuse.t       

COMMAND specifies (1) the simulation mode (either forward or backward), (2) the start and (3) the end time of the simulation, (4) the frequency T_c, (5) the averaging time DT_c and (6) the sampling intervals DT_s of the output fields, (7) the time constant for particle splitting Dt_s, (8) the synchronisation interval of , (9) the factor by which the time steps must be below the Lagrangian time scale c_tl, and (10) the refinement factor for the Langevin equation for the vertical turbulent wind component ifine. If (9) c_tl is negative, the Langevin equations are solved with constant time steps according to the synchronisation interval. The synchronisation interval is the minimum time interval used by the model for all activities (such as concentration calculations, wet deposition calculations, interpolation of data, mesoscale wind fluctuations or output of data) other than the simulation of turbulent transport and dry deposition.

********************************************************************************
*                                                                              *
*      Input file for the Lagrangian particle dispersion model FLEXPART        *
*                           Please select your options                         *
*                                                                              *
********************************************************************************

1. __                3X, I2
    1       
   DIRECTION         1 FOR FORWARD SIMULATION, -1 FOR BACKWARD SIMULATION

2. ________ ______   3X, I8, 1X, I6
   19831014 150000
   YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION

3. ________ ______   3X, I8, 1X, I6
   19831016 210000
   YYYYMMDD HHMISS   ENDING DATE OF SIMULATION

4. _____             3X, I5
   21600 
   SSSSS             OUTPUT EVERY SSSSS SECONDS

5. _____             3X, I5
   21600 
   SSSSS             TIME AVERAGE OF OUTPUT (IN SSSSS SECONDS)

6. _____             3X, I5
     900 
   SSSSS             SAMPLING RATE OF OUTPUT (IN SSSSS SECONDS)

7. ______            3X, I6
   999999 
   SSSSSS            TIME CONSTANT FOR PARTICLE SPLITTING (IN SECONDS)

8. _____             3X, I5
     900 
   SSSSS             SYNCHRONISATION INTERVAL OF FLEXPART (IN SECONDS)

9.  ---.--           4X, F6.4
      5.0
    CTL              FACTOR, BY WHICH TIME STEP MUST BE SMALLER THAN TL

10. ---              4X, I3
      4
    IFINE            DECREASE OF TIME STEP FOR VERTICAL MOTION BY FACTOR IFINE

11. -                4X, I1
    2  
    IOUT             1 CONCENTRATION OUTPUT, 2 MIXING RATIO OUTPUT, 3 BOTH

12. -                4X, I1
    0  
    POUT             PARTICLE DUMP: 0 NO, 1 EVERY OUTPUT INTERVAL, 2 ONLY AT END

13. ___              4X, A3
    YES
  TSUBGRID           SHALL SUBGRID SCALE TERRAIN EFFECTS BE PARAMETERIZED (YES or NO)?


1. Simulation direction, 1 for forward, -1 for backward in time

2. Beginning date and time of simulation. Must be given in format
   YYYYMMDD HHMISS, where YYYY is YEAR, MM is MONTH, DD is DAY, HH is HOUR,
   MI is MINUTE and SS is SECOND. Current version utilizes UTC.

3. Ending date and time of simulation. Same format as 3.

4. Average concentrations are calculated every SSSSS seconds.

5. The average concentrations are time averages of SSSSS seconds
   duration. If SSSSS is 0, instantaneous concentrations are oututted.

6. The concentrations are sampled every SSSSS seconds to calculate the time
   average concentration. This period must be shorter than the averaging time.

7. Time constant for particle splitting. Particles are split into two
   after SSSSS seconds, 2SSSSS seconds, 4SSSSS seconds, and so on.

8. All processes are synchronized with this time interval. Therefore,
   all other time constants must be multiples of this value.

9. CTL=factor by which time step must be shorter than the Lagrangian
   time scale
   CTL should be greater 1 if account of Lagrangian time scale is to
   be taken.
   If CTL<0, a purely random walk simulation is done

10.IFINE=Reduction factor for time step used for vertical wind

11 IOUT determines how the output shall be made: concentration
   (ng/m3, Bq/m3), mixing ratio (pptv), or both

12.POUT determines whether particle positions are outputted (in addition
   to the gridded concentrations or mixing ratios) or not.
   0=no output, 1 output every output interval, 2 only at end of the
   simulation

13.Switch on/off subgridscale terrain parameterization (increase of
   mixing heights due to subgridscale orographic variations)

OUTGRID specifies the output grid. Up to 29 vertical levels may be chosen. This can, however, also be changed in file includepar.

********************************************************************************
*                                                                              *
*      Input file for the Lagrangian particle dispersion model FLEXPART        *
*                       Please specify your output grid                        *
*                                                                              *
********************************************************************************

1.  ------.----       4X,F11.4
       -10.0000       GEOGRAFICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
    OUTLONLEFT        (left boundary of the first grid cell - not its centre)


2.  ------.----       4X,F11.4
        41.0000       GEOGRAFICAL LATITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
    OUTLATLOWER       (lower boundary of the first grid cell - not its centre)

3.  -----             4X,I5
      101             NUMBER OF GRID POINTS IN X DIRECTION (= No. of cells + 1)
    NUMXGRID

4.  -----             4X,I5
       47             NUMBER OF GRID POINTS IN Y DIRECTION (= No. of cells + 1)
    NUMYGRID

5.  ------.---        4X,F10.3
         0.500        GRID DISTANCE IN X DIRECTION
    DXOUTLON

6.  ------.---        4X,F10.3
         0.500        GRID DISTANCE IN Y DIRECTION
    DYOUTLAT

7.  -----.-           4X, F7.1
      100.0
    LEVEL 1           HEIGHT OF LEVEL (UPPER BOUNDARY)

8.  -----.-           4X, F7.1
      300.0
    LEVEL 2           HEIGHT OF LEVEL (UPPER BOUNDARY)

9.  -----.-           4X, F7.1
      600.0
    LEVEL 3           HEIGHT OF LEVEL (UPPER BOUNDARY)

10. -----.-           4X, F7.1
     1000.0
    LEVEL 4           HEIGHT OF LEVEL (UPPER BOUNDARY)

11. -----.-           4X, F7.1
     2000.0
    LEVEL 5           HEIGHT OF LEVEL (UPPER BOUNDARY)

12. -----.-           4X, F7.1
     3000.0
    LEVEL 6           HEIGHT OF LEVEL (UPPER BOUNDARY)

RECEPTORS specifies the receptor locations for which the kernel method shall be applied to calculate air concentrations. Up to 200 receptors may be specified:

********************************************************************************
*                                                                              *
*      Input file for the Lagrangian particle dispersion model FLEXPART        * 
*                       Please specify your receptor points                    * 
*     For the receptor points, ground level concentrations are calculated      * 
*                                                                              *
********************************************************************************
1.  ----------------  4X,A16
    F15               NAME OF RECEPTOR POINT
    RECEPTORNAME

2.  ------.----       4X,F11.4
         6.1333       GEOGRAFICAL LONGITUDE
    XRECEPTOR

3.  ------.----       4X,F11.4
        49.0833       GEOGRAFICAL LATITUDE
    YRECEPTOR
================================================================================
1.  ----------------  4X,A16
    NL01              NAME OF RECEPTOR POINT
    RECEPTORNAME

2.  ------.----       4X,F11.4
         5.7833       GEOGRAFICAL LONGITUDE
    XRECEPTOR

3.  ------.----       4X,F11.4
        50.9167       GEOGRAFICAL LATITUDE
    YRECEPTOR
================================================================================

RELEASES defines the release specifications. In the first input line, the number N of emitted species is defined (two in the example below). At all locations, the same species must be released. The next N input lines give a cross-reference to file SPECIES, where the features of the species that are released are specified. Then follow up to 10 releases specified by the trajectory type to be used, the beginning and the ending time of the release, geographical coordinates of the lower left and upper right corners of the release location, lower level and upper level of source height, the number of particles to be used, and the total mass emitted. Please note that the mass specification must be made for as many species (N) as have been specified above. Finally, the name of the release point must be given.

The particles are released evenly distributed within a volume stretching from the lower to the upper level (both given in m above ground) above a rectangular area (in grid coordinates) defined by the geographical coordinates of the lower left and upper right corners. In the following example, two release points are specified, the first a volume source, the second a ground level point source.

*************************************************************************
*                                                                       *
*   Input file for the Lagrangian particle dispersion model FLEXPART    *
*                        Please select your options                     *
*                                                                       *
*   Kind of trajectory: 1 = 3-dimensional                               *
*                       2 = 2-dimensional (terrain-following)           *
*                                                                       *
*                                                                       *
*************************************************************************
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2
___                        i3    Total number of species emitted

  3
___                        i3    Index of species in file SPECIES

  5
___                        i3    Index of species in file SPECIES

=========================================================================
 1
 _                         1X,I1 Kind of trajectory (see file header)

19941023 160000
________ ______            i8,1x,i6 Beginning date and time of release

19941024 035000
________ ______            i8,1x,i6 Ending date and time of release

  -2.50
____.____                  f9.4  Longitude [DEG] of lower left corner

  48.00
____.____                  f9.4  Latitude [DEG] of lower left corner

  -2.20
____.____                  f9.4  Longitude [DEG] of upper right corner

  48.40
____.____                  f9.4  Latitude [DEG] of upper right corner

   50.0
_____.___                  f10.3 Lower z-level (in m above ground)
 
  150.0
_____.___                  f10.3 Upper z-level (in m above ground)

    30000
_________                  i9    Total number of particles to be released

3.0000E00
_.____E__                  e9.4  Total mass emitted

5.0000E00
_.____E__                  e9.4  Total mass emitted

TEST1
________________________________________   character*40 comment
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1
 _                         1X,I1 Kind of trajectory (see file header)

19941023 180000
________ ______            i8,1x,i6 Beginning date and time of release

19941024 030000
________ ______            i8,1x,i6 Ending date and time of release

  -2.50
____.____                  f9.4  Longitude [DEG] of lower left corner

  48.00
____.____                  f9.4  Latitude [DEG] of lower left corner

  -2.50
____.____                  f9.4  Longitude [DEG] of upper right corner

  48.00
____.____                  f9.4  Latitude [DEG] of upper right corner

    0.0
_____.___                  f10.3 Lower z-level (in m above ground)
 
    0.0
_____.___                  f10.3 Upper z-level (in m above ground)

    10000
_________                  i9    Total number of particles to be released

1.0000E00
_.____E__                  e9.4  Total mass emitted
 
1.0000E00
_.____E__                  e9.4  Total mass emitted
 
TEST2
________________________________________   character*40 comment
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

SPECIES gives a list of species from which the user can select the ones to be simulated. The file contains the physical properties of these species used to simulate radioactive decay, and wet and dry deposition. For the wet deposition, A and B specify the factors defined by equation 7.3.

For dry deposition of gases, D = D_H2O/D_i, where D_H2O is the diffusivity of water vapor and D_i is the diffusivity of the species, H is the effective Henry's constant, and f0 varies between 0 and 1 and gives the reactivity of a species relative to that of ozone. For nonreactive species f0 is 0, for slightly reactive it is 0.1 and for highly reactive it is 1.

For the dry deposition of particulates, rho specifies the density of the substance, dquer its mean diameter [`(d_p)], and dsig the measure of variation s_p.

Radioactive decay is switched off by specifying a negative half life, wet deposition is switched off by specifying negative A, dry deposition of gases is switched off by negative D, dry deposition of particles is switched off by negative rho.

molweight gives the molecular weight of the species, which is needed for mixing ratio output.

****************************************************************************
*                                                                          *
*     Input file for the Lagrangian particle dispersion model FLEXPART     *
*               Definition file of chemical species/radionuclides          *
*                                                                          *
****************************************************************************
           Radioactivity    Wet depo         Dry depo (gases)      Dry depo (particles)  Dry depo
     SPECIES HALF LIFE [s]     A      B       D      H     f0      rho    dquer    dsig      vd   molweight
  1 TRACER        -999.9   -9.9E-09         -9.9                 -9.9E09                   -9.99  350.00 
  2 O3            -999.9   -9.9E-09          1.5  1.0e-02  1.0   -9.9E09                   -9.99   48.00
  3 NO            -999.9    8.0E-06  0.62    1.2  2.0e-03  0.0   -9.9E09                   -9.99   30.00
  4 NO2           -999.9    1.0E-05  0.62    1.6  1.0e-02  0.1   -9.9E09                   -9.99   46.00
  5 HNO3          -999.9    8.0E-04  0.62    1.9  1.0e+14  0.0   -9.9E09                   -9.99   63.00
  6 HNO2          -999.9   -9.9E-09          1.6  1.0e+05  0.1   -9.9E09                   -9.99   47.00
  7 H2O2          -999.9    1.0E-04  0.62    1.4  1.0e+05  1.0   -9.9E09                   -9.99   34.00       
  8 SO2           -999.9   -9.9E-09  0.62    2.0  1.0e+05  0.0   -9.9E09                   -9.99   64.00
  9 HCHO          -999.9   -9.9E-09          1.3  6.0e+03  0.0   -9.9E09                   -9.99   30.00
 10 PAN           -999.9   -9.9E-09          2.6  3.6e+00  0.1   -9.9E09                   -9.99  121.00
 11 NH3           -999.9   -9.9E-09          1.1  2.0e+14  0.0   -9.9E09                   -9.99   17.00
 12 SO4-aero      -999.9    1.0E-04  0.80   -9.9                  2.0E03  4.0E-7  3.0E-1   -9.99   -9.99
 13 NO3-aero      -999.9    1.0E-04  0.80   -9.9                  2.0E03  4.0E-7  3.0E-1   -9.99   -9.99
 14 I2-131      691200.0    8.0E-05  0.62    2.7  1.0e+05  0.1   -9.9E09                   -9.99   -9.99
 15 I-131       691200.0    1.0E-04  0.80   -9.9                  2.5E03  6.0E-7  3.0E-1   -9.99   -9.99
 16 Cs-137        -999.9    1.0E-04  0.80   -9.9                  2.5E03  6.0E-7  3.0E-1   -9.99   -9.99
 17 Y-91       5037120.0    1.0E-04  0.80   -9.9                  2.5E03  6.0E-7  3.0E-1   -9.99   -9.99
 18 Ru-106    31536000.0    1.0E-04  0.80   -9.9                  2.5E03  6.0E-7  3.0E-1   -9.99   -9.99
 19 Kr-85         -999.9   -9.9E-09         -9.9                 -9.9E09                   -9.99   -9.99
 20 Sr-90         -999.9    1.0E-04  0.80   -9.9                  2.5E03  6.0E-7  3.0E-1   -9.99   -9.99
 21 Xe-133      198720.0   -9.9E-09         -9.9                 -9.9E09                   -9.99   -9.99

landuse.t is a binary file that contains the landuse inventory. surfdata.t gives the roughness lengths for each landuse class:

8 landuse categories are related with Leaf Area Index and roughness length
--------------------------------------------------------
landuse    LAI         z0(m)       Albedo     comment
--------------------------------------------------------
 1         0.50       0.10          0.18      Grassland for agricultural use
 2         2.00       0.15          0.10      Arable land
 3         3.00       0.30          0.18      Permanent crops
 4         7.00       0.60          0.15      Forest
 5         0.00       0.10          0.12      Inland water
 6         0.20       0.70          0.20      Urban areas
 7         1.00       0.10          0.15      Other
 8         0.00       0.10          0.12      Ocean

surfdepo.t gives the resistances needed for the parameterization of dry deposition of gases for the eight landuse classes and five seasonal categories. This file must not be changed by the user.

==============================================================================
INPUT RESISTANCES (s/m) FOR THE COMPUTATION OF SURFACE RESISTANCES TO
DRY DEPOSITION
==============================================================================
AFTER WESELY, 1989
==============================================================================
1 to 8: Landuse types
==============================================================================
Values are tabulated for 5 seasonal categories:
1     Midsummer with lush vegetation
2     Autumn with unharvested cropland
3     Late autumn after frost, no snow
4     Winter, snow on ground and subfreezing
5     Transitional spring with partially green short annuals
==============================================================================
               1       2       3       4       5       6       7       8 
______________________________________________________________________________
ri          120.     60.     65.     90.   9999.   9999.    150.   9999.     1   
rlu        2000.   2000.   2000.   2000.   9999.   9999.   4000.   9999.
rac         100.    200.    365.   2000.      0.    100.    200.      0. 
rgss        350.    150.    230.    500.      0.    400.    400.      0. 
rgso        200.    150.    170.    200.   2000.    300.    200.   2000.
rcls       2000.   2000.   2000.   2000.   9999.   9999.   4000.   9999.
rclo       1000.   1000.   1000.   1000.   9999.   9999.   1000.   9999.
______________________________________________________________________________
ri         9999.   9999.   9999.    500.   9999.   9999.   9999.   9999.     2   
rlu        9000.   9000.   9000.   5500.   9999.   9999.   9000.   9999.
rac         100.    150.    270.   1710.      0.    100.    140.      0. 
rgss        350.    200.    285.    500.      0.    400.    400.      0. 
rgso        200.    150.    170.    200.   2000.    300.    200.   2000.
rcls       9000.   9000.   9000.   3270.   9999.   9999.   9000.   9999.
rclo        400.    400.    400.    570.   9999.   9999.    400.   9999.
______________________________________________________________________________
ri         9999.   9999.   9999.    500.   9999.   9999.   9999.   9999.     3   
rlu        9000.   9999.   9470.   5500.   9999.   9999.   9000.   9999.
rac         100.     10.     20.   1330.      0.    100.    120.      0. 
rgss        350.    150.    230.    500.      0.    400.    400.      0. 
rgso        200.    150.    170.    200.   2000.    300.    200.   2000.
rcls       9000.   9999.   9470.   4500.   9999.   9999.   9000.   9999.
rclo        400.   1000.    570.    570.   9999.   9999.    600.   9999.
______________________________________________________________________________
ri         9999.   9999.   9999.    800.   9999.   9999.   9999.   9999.     4   
rlu        9999.   9999.   9999.   1200.   9999.   9999.   9000.   9999.
rac          10.     10.     20.   1330.      0.    100.     50.      0. 
rgss        100.    100.    100.    100.      0.    100.     50.      0. 
rgso       3500.   3500.   3500.   3500.   2000.    600.   3500.   2000.
rcls       9999.   9999.   9470.    390.   9999.   9999.   9000.   9999.
rclo       1000.   1000.    570.    632.   9999.   9999.    800.   9999.
______________________________________________________________________________
ri          240.    120.    130.    180.   9999.   9999.    300.   9999.     5   
rlu        4000.   4000.   4000.   2670.   9999.   9999.   8000.   9999.
rac          80.     50.    100.   1500.      0.    100.    120.      0. 
rgss        350.    150.    230.    500.      0.    500.    400.      0. 
rgso        200.    150.    170.    200.   2000.    300.    200.   2000.
rcls       4000.   4000.   4000.   2670.   9999.   9999.   8000.   9999.
rclo        500.   1000.    670.    750.   9999.   9999.    800.   9999.
______________________________________________________________________________

11.4  Output files

For every field and at every output time it is checked whether it requires less mass storage space to dump only the values for grid points containing non-zero concentrations than to dump the whole field at once. In the former case, in addition to the concentrations, also the grid indices have to be dumped (obviously, otherwise, this method would always be more efficient). For this purpose, all three indices are combined into a single integer number to save storage space. As long as the tracer cloud is small compared to the output domain, this sparse matrix format is more efficient than the full field dump. In the output directory, there is a short program (flextrans.f), which converts the binary output to ASCII output. From this example, it can be seen how the FLEXPART output can be read into a FORTRAN program.

The concentrations and depositions are multiplied by 10¹² for output. This means that if the released mass is given in kg, concentrations are given in ngm^-3 and deposition in ngm^-2. If the release is given in TBq, then the output is in Bqm^-3 and Bqm^-2, respectively. For mixing ratio, this means that if the released mass is given in kg, the mixing ratios are given in pptv.

Chapter 12
Source code description

SHORT DESCRIPTION OF THE SOURCE CODE FILES USED BY THE FLEXPART MODEL
_____________________________________________________________________
advance.f         advection of the trajectory using a first order scheme for one
assignland.f      assigns the fractions of 8 landuse classes to each ECMWF grid
                  point
calcpar.f         calculates boundary layer parameters: friction velocity,
                  convective velocity scale, boundary layer height, Obukhov length
calcpar_nests.f   same as calcpar.f, but for the nested grids
caldate.f         compute calendar date from julian date
cmapf1.0.f        Program package for coordinate transformations between
                  latitude/longitude grid and polar stereographic projection
coordtrafo.f      transformation of release point coordinates from geografical
                  to grid coordinates
dist2.f           calculation of distance (m) between two points given in 
                  geografical coordinates
drydepokernel.f   calculates gridded dry deposition field using a uniform kernel
                  of bandwidths dx, dy
erf.f             error function
ew.f              calculation of saturation vapor pressure at given temperature
                  synchronisation interval; uses mean wind and Langevin equations
                  for parameterizing turbulence and mesoscale wind fluctuations
fixparticles.f    determine heights, release times and masses of all particles
FLEXPART.f        main program, manages reading of run specifications, calls
                  timemanager to do the actual calculations
getfields.f       management of wind fields in memory
getrb.f           computation of quasilaminar sublayer resistance to dry
                  deposition of gases
getrc.f           calculation of surface resistance to dry deposition of gases
getvdep.f         calculation of dry deposition velocities
gridcheck.f       determines grid specifications (grid area, etc.) from GRIB file
gridcheck_nests.f same as gridcheck.f, but for the nested grids
gridconc.f        calculation of concentrations on a grid using volume sampling
                  and at receptor points using kernel method
gridout.f         output of concentrations on grid and for receptor points
hanna.f           parameterization of the turbulent velocity variance and
                  the respective timescales; Langevin equation for w'/sigw
hanna1.f          alternative model to hanna.f: Langevin equation for w'
                  drift correction velocity after Wilson
hanna_short.f     same as hanna.f, but only for the vertical wind
includecom        include file containing all global variables
includeinterpol   include file containing variables used for interpolation
includehanna      include file used for turbulence parameterization
includepar        include file containing all globally used parameter statements
initialize.f      initialization of the turbulence statistics for each particle
interpol_all.f    interpolation of all data needed by advance for two levels within
                  the boundary layer and for the surface
interpol_all_nests.f  same as interpol_all.f, but for the nested grids
interpol_misslev.f interpolates all levels not yet interpolated by interpol_all.f
interpol_misslev_nests.f same as interpol_misslev.f, but for the nested grids
interpol_uvrho.f  interpolation of u, v and rho for an additional level
interpol_wrhograd.f interpolation of w and rhograd for an additional level
interpol_vdep.f   interpolation of the deposition velocity
interpol_vdep_nests.f  same as interpol_vdep.f, but for the nested grids
interpol_wind.f   interpolation of wind data above the PBL
interpol_wind_nests.f   same as interpol_wind.f, but for the nested grids
juldate.f         compute julian date from calendar date
levlin3interpol.f linear interpolation of three fields to a point on a model layer
                  and in time
mean.f            calculation of mean and standard deviation
obukhov.f         calculation of Obukhov length from surface heat flux
openoutput.f      opens the output files for gridded and receptor point output
openoutputppt.f   same as openoutput.f, but for the mixing ratio output rather
                  than the concentration output
openpartoutput.f  opens the output files for the particle position dumps
part0.f           calculation of time independent factors for the dry deposition
                  of particles
partdep.f         calculation of dry deposition velocities for particles
partout.f         dumps the particle positions
pbl_profile.f     calculation of surface stress and sensible heat flux using
                  the profile method
psih.f            stability correction term for heat
psim.f            stability correction term for momentum
raerod.f          calculation of the aerodynamical resistance to dry deposition
random.f          subroutines for generation of normally distriubted
                  random numbers
readavailable.f   read, which wind fields are available within the modelling period
readcommand.f     read user commands, e.g. start of release, duration of release
readdepo.f        read parameters needed for dry deposition of gases
readlanduse.f     read landuse inventory and the respective roughness lenghts
readoro.f         read ecmwf model orography and land sea mask
readoutgrid.f     read output grid specifications
readpaths.f       read pathnames, where input/output files are stored
readreceptors.f   read names and coordinates of receptor points; for receptor
                  points, concentrations are calculated with the kernel method 
readreleases.f    read release point specifications (e.g., coordinates, type of
                  trajectories to be used, number of particles, mass, etc.)
readspecies.f     read physical and chemical properties of species
readwind.f        reads the ECMWF meteorological data fields; GRIB decoding
readwind_nests.f  same as readwind.f, but for the nested grids
richardson.f      calculation of convective velocity scale and boundary layer
                  height using critical Richardson number concept
scalev.f          computation of friction velocity from surface stress
sigma.f           calculation of standard deviation of a 1-D field
timemanager.f     time management of FLEXPART; determines when a integration time
                  step of a trajectory is due; determines when a concentration
                  calculation is due, etc, and calls the respective subprograms
verttransform.f   transformation of three-dimensional fields from eta coordinates
                  to meter above ground coordinates
verttransform_nests.f  same as verttransform.f, but for the nested grids
wetdepo.f         calculation of wet deposition
wetdepokernel.f   calculates gridded wet deposition field using a uniform kernel
                  of bandwidths dx, dy
windalign.f       transformation of random velocities from along- and cross-wind
                  components to Cartesian u and v components

******************************************************************************
*            FLEXPART model structure for ECMWF wind fields.                 *
******************************************************************************

FLEXPART --> readpaths

         --> readcommand --> juldate

         --> readavailable --> juldate

         --> gridcheck

         --> gridcheck_nests

         --> readoutgrid

         --> readreceptors

         --> readspecies

         --> readoro

         --> readlanduse

         --> assignland

         --> readreleases --> part0

         --> readdepo

         --> coordtrafo

         --> fixparticles

         --> openoutput --> caldate

                        --> dist2

         --> openoutputppt --> caldate

                           --> dist2

         --> openpartoutput --> caldate

         --> gasdev1 --> ran3

         --> timemanager --> getfields --> %1

                         --> wetdepo --> levlin3interpol
 
                                     --> wetdepokernel
                                      
                         --> gridconc

                         --> gridout

                         --> partout

                         --> advance --> ran3

                                     --> interpol_all --> sigma

                                     --> interpol_all_nests --> sigma

                                     --> interpol_misslev --> sigma

                                     --> interpol_misslev_nests --> sigma

                                     --> hanna or hanna1

                                     --> hanna_short

                                     --> interpol_vdep

                                     --> interpol_vdep_nests

                                     --> interpol_wind --> sigma

                                     --> interpol_wind_nests --> sigma

                                     --> windalign
                                      
                         --> drydepokernel

______________________________________________________________________________
%1

 --> readwind (--> pbl_profile --> psim)

 --> readwind_nests (--> pbl_profile --> psim)

 --> calcpar --> scalev --> ew

             --> obukhov

             --> richardson

             --> getvdep --> caldate

                         --> getrb

                         --> raerod --> psih

                         --> getrc

                         --> partdep

 --> calcpar_nests --> scalev --> ew

             --> obukhov

             --> richardson

             --> getvdep --> caldate

                         --> getrb

                         --> raerod --> psih

                         --> getrc

                         --> partdep

 --> verttransform

 --> verttransform_nests

______________________________________________________________________________