ECMWF DATA EXTRACTION FOR FLEXTRA/FLEXPART

Extraction of ECMWF data for the trajectory model FLEXTRA and the particle dispersion model FLEXPART

NEW UPDATE

Step 1: Retrieval and de-accumulation of flux data (precipitation, radiation, sensible heat flux,...)
Step 2: Retrieval of other meteorological fields, calculation of vertical velocity component d(eta)/dt

If you need input data for FLEXTRA only, you can leave away step 1

Warning: In order to use these routines, it is necessary to have a user account on ecgate1.ecmwf.int

For more information: send email to Paul James

The tar file ecdocs.tar

All routines and documentation needed are enclosed in the file ecdocs.tar:

Download ecdocs.tar (379 kB)
Transfer to ecgate1
Extract directories and files: tar -xf ecdocs.tar

Step 1: Retrieval and de-accumulation of flux data

The programs can be found in the directory ./fluxjobs_gen and ./fluxjobs_era (for ERA-15 data) .

The following flux data needed by FLEXPART are retrieved on a user-selectable geographical grid (lamda-phi):

Abreviation	Full name	Unit	Field code
LSP	Large Scale Precipitation	mm/hr	142
CP	Convective Precipitation	mm/hr	143
SSHF	Surf. Sens. Heat Flux	W/m2	146
EWSS	East-West Surf. Stress	N/m2	180
NSSS	North-South Surf. Stress	N/m2	181
SSR	Surf. Solar Radiation	W/m2	176

Flux data are not available in the ECMWF analyses, but they can be obtained (in accumulated form) from the ECMWF short-term forecasts (first guess) stored in three-hourly intervals (see ECMWF MARS database manual).

For de-accumulation, the following strategy is followed:

Precipitation (LSP, CP)

Accumulated values for 3/6/9/... UTC are just divided by the number of hours, namely 3
Therefore, Precipitation at X UTC means: Mean precipitation per hour from X-3 to X UTC

Other flux quantities

Accumulated values are not only divided by the number of hours, but also interpolated to 3/6/9... UTC (for a three-hourly first guess with evaluation time X UTC, a flux value has been accumulated from X-3 UTC to X UTC and is therefore valid for X-1.5 UTC)
Therefore, e.g., Radiation (X UTC) means: Radiation valid for X UTC (three-hourly average)

If you want to work with these routines, you should read the following instructions:

Create a scratch-directory to save your flux data

e.g. mkdir $SCRATCH/flux
or for ERA-15 data mkdir $SCRATCH/flux_era

After each modification of the grid:

Modify the FORTRAN include file parameter according to your new grid specifications

nx, ny: number of grid points
xlon0, ylat0: longitude and latitude of lower left grid point
dx, dy: grid distance

Compile and link de-accumulation routines (program FLUXACC)

f77 *.f $EMOSLIB
mv a.out FLUXACC

Retrieve the accumulated flux data from the MARS archive in batch mode (surfjob.new)

Only first time: Modify the last line of surfjob.new according to your own ecfile structure created

e.g.: ecp -o surf_* ec:/uid/metdat

Modify beginning and ending day

e.g. DAY=19970526/to/19970528
In the given example, the following data will be retrieved: 26 May 1997 03 UTC to 28 May 1997 18 UTC

Put your job in batch queue

e.g.: qsub -q ecgate1.normal surfjob.new
You will be notified via e-mail if your job starts and after it is finished
Look at the job log file in ~uid/waitqueue named SURFdataXXXXXX with XXXXXX being a system-generated number

Deaccumulation of flux data, storage of results in ecfile system

For that purpose, use the script deacc.flux (interactive execution)
Modify deacc.flux before first usage:

Enter ecfile directory where the input to FLUXACC is stored (see file surfjob.new).
Enter scratch directory where the output of FLUXACC shall be stored.

After succesful execution of deacc.flux, the data will be available in the selected scratch directory in three (six for ERA data) hourly intervals. The input to FLUXACC is deleted. The filenames are fluxYYYYMMDDHH with YYYYMMDDHH being the date. These data are used in step 2 of the data retrieval procedure for FLEXTRA/FLEXPART

Step 2: Retrieval of other meteorological fields, calculation of the vertical velocity component d(eta)/dt

The programs can be found in the directories ./job_era ./job_pre99 ./job_y2k

This step has been fully automated. The user just has to fill in a CONTROL file, providing information like beginning date, ending date, time interval, number of levels to extract, scratch directory for (optional) flux data input, ecfile directory for output. Afterwards, the execution of a.out creates a batch script and sends it to batch queue ecgate1.normal.

The batch file is saved under /tmp/job.prepo.bak.uid. After execution, a job log file in ~uid/waitqueue named METEprepXXXXXX with XXXXXX being a system-generated, unique number, is created. Beginning and end of execution are anounced to the user via e-mail.

The batch-job works as follows:

It creates an executable CONVERT_PRE from source-files made available (you do not need to have them in your own directories, and they are therefore not included in the tar file)
It extracts ECMWF model data (U, V, T, Q, LNSP, ...) from the MARS archive
If available, it extracts your prepared flux data from the specified scratch directory
It executes the routine CONVERT_PRE, which serves the folowing purposes:

It reads in all ECMWF model data
It calculates gridded values of the vertical velocity component in the ECMWF eta vertical coordinate system (d(eta)/dt) applying the equation of continuity, making the wind fields 3D mass consistent.
It writes out all meteorological data including d(eta)/dt

It saves the resulting output files in the user-specified ecfile directory. The files are ECMWF grib format (for grib decoding routines, ask ECMWF user supports or people from your national weather service). The files are named ENYYMMDDHH.

To use these routines, please follow the instructions step by step:

Before first execution, create an ecfile directory for your output

e.g.: emkdir ec:/uid/metdat

Fill in the file CONTROL
execute a.out
After the batch job has finished, have a look at your log file
At last, you can transfer the data to your local machine. Please mind that transfer of grib encoded files via ftp requires the option binary. Otherwise, the files get corrupt and have to be transfered again.

Technical comments:

The routine has been tested for many grid domains in the northern hemisphere. It can also be used to retrieve hemispheric and global data.
The following important restrictions concerning the selection of the grid domain do apply:

Specify domain boundaries only as whole numbers (e.g. -5.0 ; Never: -5.5).
If the upper right boundary longitude is a positive number, the lower left boundary longitude has to be negative.

For technical reasons, you can only chose grid configurations with dx=dy.
Please mind that d(eta)/dt can not be calculated for the Poles. It is set to the values from the previous latitude.
This routine should work for all dates from 4 April 1995 onwards (ERA-15 version also from 1979 onwards). Please modify the spectral truncation and the maximum number of layers in the CONTROL file according to the information provided there (or have a closer look into your MARS manual).
You can (and should) select a lower spectral truncation than possible in the CONTROL file if you retrieve data on a coarser grid. Please look at the comments in the CONTROL file concerning the relationship between grid resolution and truncation.
Please mind: If you change the grid domain in job_xxx, and if you need flux data, you have to retrieve new flux data for the modified domain.
If you retrieve long periods at once (more than a month), check the status of your job at times. If it has, e.g., produced no output for more than 24 hours, you should terminate it (command qdel) and start it again for the time period not yet finished.

Update August 1998 - Retrieval of additional parameters, 3 prepared domains

In August 1998, the retrieval routines (Step 2) were revised to get additional quantities potentially useful for particle modelling, air quality modelling and mesoscale meteorological modelling. The following additional data are now provided:

Abreviation	Field code	Name
SDOR	160	Std. Deviaton of Orography
VEG	199	Percentage of vegetation
SR	173	Surface roughness
LCC	186	Low cloud cover
MCC	187	Medium cloud cover
HCC	188	High cloud cover
SKT	235	Skin temperature
STL1	139	Soil temperature level 1
SWL1	140	Soil wetness level 1

Please notice: The new routines just work for ECMWF data from 04/04/1995 0 UTC onwards. If you want to extract older ECMWF data, use older routines written by Gerhard Wotawa.

If you want to use these old routines, just download them.

Handling of wind components at the poles

At the south Pole (phi=-90 DEG) and the north Pole (phi=90 DEG), all scalar quantities (e.g., temperature) are set to the same value for -178<= lambda <= 180 DEG (i=1,NX) on the lambda-phi grid. The wind components u, v are handled according to the following convention:

North Pole:

u(lambda)=-ff sin(lambda+dd)

v(lambda)=-ff cos(lambda+dd)

where lambda are the geographical longitude (-pi to + pi), ff the wind speed at the pole, dd the wind direction at the pole. Both, ff and dd, are constant for the Pole.

The wind direction is calculated as follows:

dd(lambda)=atan(u(lambda)/v(lambda))-lambda [if v(lambda)<0]

dd(lambda)=pi+atan(u(lambda)/v(lambda))-lambda [else]

For example: A wind direction of pi (180 DEG) means that the wind at the pole blows along the Greenwich meridian towards the pole. If the wind direction is 0, the wind blows the other way round.

South Pole:

u(lambda)=+ff sin(lambda-dd)

v(lambda)=-ff cos(lambda-dd)

The wind direction is calculated as follows:

dd(lambda)=atan(u(lambda)/v(lambda))+lambda [if v(lambda)<0]

dd(lambda)=pi+atan(u(lambda)/v(lambda))+lambda [else]

For a better understandig, I provide a little program where wind at the south and north pole is transformed to speed and direction and subsequently transformed back according to the pole convention used by ECMWF:

* SOUTH POLE 
 
      DO 666 I=1,NX
         LLON=-178+FLOAT(I-1)*DX
         LLON1=LLON*PI/180.
         WSPEED=SQRT((U3D(I,1,1)**2)+(V3D(I,1,1)**2))
         IF(V3D(I,1,1).lt.0.) THEN
           DDPOL=ATAN(U3D(I,1,1)/V3D(I,1,1))+LLON1
         ELSE
           DDPOL=PI+ATAN(U3D(I,1,1)/V3D(I,1,1))+LLON1
         ENDIF
         IF(DDPOL.LT.0.) DDPOL=2.0*pi+DDPOL
         IF(DDPOL.GT.2.0*PI) DDPOL=DDPOL-2.0*pi
         UNEW=+WSPEED*SIN(LLON1-DDPOL)
         VNEW=-WSPEED*COS(LLON1-DDPOL)
666      WRITE(999,778) I,LLON,U3D(I,1,1),V3D(I,1,1),WSPEED,
     &         UNEW,VNEW

*NORTH POLE 

            DO 777 I=1,NX
         LLON=-178+FLOAT(I-1)*DX
         LLON1=LLON*PI/180.
         WSPEED=SQRT((U3D(I,NY,1)**2)+(V3D(I,NY,1)**2))
         IF(V3D(I,NY,1).lt.0.) THEN
           DDPOL=ATAN(U3D(I,NY,1)/V3D(I,NY,1))-LLON1
         ELSE
           DDPOL=PI+ATAN(U3D(I,NY,1)/V3D(I,NY,1))-LLON1
         ENDIF
         IF(DDPOL.LT.0.) DDPOL=2.0*pi+DDPOL
         IF(DDPOL.GT.2.0*PI) DDPOL=DDPOL-2.0*pi
         UNEW=-WSPEED*SIN(LLON1+DDPOL)
         VNEW=-WSPEED*COS(LLON1+DDPOL)
777      WRITE(999,778) I,LLON,U3D(I,NY,1),V3D(I,NY,1),WSPEED,
     &         UNEW,VNEW
778      FORMAT(i3,1x,f5.0,1x,5f7.2)

Important remark: A spectral truncation of T106 is default in the ECMWF mars system for the retrieval of 2x2 DEG lambda-phi data out of the Sperical Harmonics (SH) data. However, as I have tested, this truncation is insufficient to obtain accurate wind data at the poles. If you use the latest (6.1.1999) version of my routines, this problem is fixed. If not, make sure that U,V wind components are retrieved at least with T213 or T319. This problem is serious. Wind data at the poles are definitely unusable with the lower truncation.

Bug reports:

(1) Date: 12.1.1999

Modification: File makeinclude.f in directory job.tra1995+ modified

Please note: Values of etapoint and log surface pressure are wrong in the global 2x2 grid as well as in the 0.4 degrees grid prior to this modification.

Restrictions

This software can be used "as is" by all people having user accounts at the ECMWF workstation ecgate1.ecmwf.int, without any fee or cost. Any warranty is disclaimed and no support is given. This software is part of the FLEXTRA and FLEXPART program package. The same regulations as accepted by using FLEXTRA and FLEXPART do apply. Please keep in mind the data usage regulations of the European Centre for Medium Range Weather Forecasts and of your national weather service. Any use of this program, which is against these regulations, is prohibited.

You are kindly invited to report problems and bugs to me.

Appendix: Plot of sample fields

Last update: 7 August 2000 by Paul James