ECMWF DATA EXTRACTION FOR FLEXTRA/FLEXPART

IMP	Institute of Meteorology and Physics University of Agricultural Sciences, Vienna, Austria
	*Environmental Meteorology Group*

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

ATTENTION:
The version presented here is obsolete.
It has been superseded by a new version
prepared by Paul James in Munich. The new version takes care of the following aspects:

Improved deaccumulation of fluxes

Replacement of FG (first guess) fields by 4V (4DVAR analysis) and FC (short-term forecasts, up to +9h), reflecting changes in the operational procedures at ECMWF

Get this new version here.

Petra Seibert, 31 Jan 2001

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 Gerhard Wotawa

Description of tar file

All routines needed are enclosed in the file ecmwf.tar:

fluxjobs/                       (Retrieval/de-accumulation of flux data)
fluxjobs/job.flux
fluxjobs/caldate.f
fluxjobs/commonvar
fluxjobs/deaccumulate.f
fluxjobs/fluxinterpol.f
fluxjobs/juldate.f
fluxjobs/parameter
fluxjobs/readflux.f
fluxjobs/surfjob.new
fluxjobs/writeout.f

job.tra/                        (Retrieval of met. fields, calculation of d(eta)/dt
job.tra/caldate.f
job.tra/create_job.f
job.tra/execute_job.f
job.tra/extract.f
job.tra/juldate.f
job.tra/makeinclude.f
job.tra/CONTROL
job.tra/includepar

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

Step 1: Retrieval and de-accumulation of flux data

The programs can be found in the directory ./fluxjobs.

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 an ecfile-directory to save your flux data (ecfile: mass storage/data handling system established at ECMWF

e.g. emkdir ec:/uid/flux (with uid being your personal user id)

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
Plese notice: Do not retrieve more than 24 days within one job
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 job.flux (interactive execution)
Modify job.flux before first usage:

Enter ecfile directory where the input to FLUXACC is stored (see file surfjob.new).
Enter ecfile directory where the output of FLUXACC shall be stored.
In the given example file, input and output of FLUXACC are stored in different directories. But, these guys could also be stored in the same.

After succesful execution of job.flux, the data will be available in the selected ecfile directory in three 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 directory ./job.tra

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, ecfile 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 ecfile 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.

A sample CONTROL file is printed in Appendix 2.

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 20 April 1983 onwards. Please modify the spectral truncation and the maximum number of layers (currently: 31) 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.tra, and if you need flux data, you have to retrieve new flux data for the modified domain.
Do not retrieve more than one month of data in one job.
If you retrieve longer periods, 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 the previous version of job.tra.

If you want to use these routines, just download them. (NEW: 12.1.1999 -> Bug reports!!!)

The file ecmwf2.tar has the following structure:

job.tra1995+/CONTROL
job.tra1995+/CONTROL_0.5
job.tra1995+/CONTROL_1.0
job.tra1995+/CONTROL_2.0
job.tra1995+/caldate.f
job.tra1995+/create_job.f
job.tra1995+/execute_job.f
job.tra1995+/extract.f
job.tra1995+/juldate.f
job.tra1995+/makeinclude.f
fluxjobs_0.5/caldate.f
fluxjobs_0.5/deaccumulate.f
fluxjobs_0.5/fluxinterpol.f
fluxjobs_0.5/juldate.f
fluxjobs_0.5/readflux.f
fluxjobs_0.5/writeout.f
fluxjobs_0.5/parameter
fluxjobs_0.5/commonvar
fluxjobs_0.5/surfjob.new
fluxjobs_0.5/job.flux
fluxjobs_1.0/caldate.f
fluxjobs_1.0/deaccumulate.f
fluxjobs_1.0/fluxinterpol.f
fluxjobs_1.0/juldate.f
fluxjobs_1.0/readflux.f
fluxjobs_1.0/writeout.f
fluxjobs_1.0/parameter
fluxjobs_1.0/commonvar
fluxjobs_1.0/surfjob.new
fluxjobs_1.0/job.flux
fluxjobs_2.0/caldate.f
fluxjobs_2.0/deaccumulate.f
fluxjobs_2.0/fluxinterpol.f
fluxjobs_2.0/juldate.f
fluxjobs_2.0/readflux.f
fluxjobs_2.0/writeout.f
fluxjobs_2.0/parameter
fluxjobs_2.0/commonvar
fluxjobs_2.0/surfjob.new
fluxjobs_2.0/job.flux

Three different domains have been prepared for the retrieval of meteorological data:

DOMAIN 1: Global data 2.0x2.0 degrees (178 W - 180 E; 90 S - 90 N), all levels
DOMAIN 2: Europe 1.0x1.0 degrees (50 W - 50 E; 25 N - 80 N), all levels
DOMAIN 3: Central Europe 0.4x0.4 degrees (6 W - 26 E; 40 N - 60 N), all levels

If you would like to retrieve data, e.g., for domain 2, you have to proceed as follows:

cd ./fluxjobs_1.0

retrieve flux data (surfjob.new)
execute job.flux

cd ./job.tra1995+

cp CONTROL_1.0 CONTROL
a.out

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 1: Retrieval and de-accumulation of flux data from ECMWF: Sample files surfjob.new and job.flux

surfjob.new:

#QSUB -r SURFdata
#QSUB -s /bin/ksh
#QSUB -eo
#QSUB -lM 200Mb
#QSUB -lm 200Mb
#QSUB -lT 400
#QSUB -lt 400
#QSUB -mb
#QSUB -me
#QSUB

cd $TMPDIR

PATH=$PATH:.
DAY=19970526/to/19970528
HOUR1=00/06/12/18
HOUR2=03/09/15/21
AREALF=80.0/-50.0/25.0/50.0
GRIDLF=1./1.

export PATH
export DAY
export HOUR1
export HOUR2
export AREALF
export GRIDLF
cd $SCRATCH

cat >zrequest<<EOF
RETRIEVE,
TYPE=FG,
PARAM=LSP/CP/SSHF/EWSS/NSSS/SSR,
AREA=$AREALF,
GRID=$GRIDLF,
LEVTYPE=SFC,
LEVELIST=OFF,
REPRES=GG,
DATE=$DAY,
TIME=$HOUR1,
AC=N,
STEP=06,
TARGET="surf_06_ub"
END
EOF
mars zrequest

cat >zrequest<<EOF
RETRIEVE,
TYPE=FG,
PARAM=LSP/CP/SSHF/EWSS/NSSS/SSR,
AREA=$AREALF,
GRID=$GRIDLF,
LEVTYPE=SFC,
LEVELIST=OFF,
REPRES=GG,
DATE=$DAY,
TIME=$HOUR2,
AC=N,
STEP=03,
TARGET="surf_03_ub"
END
EOF
mars zrequest

ecp -o surf_* ec:/vo3/metdat/

job.flux:

#! /bin/ksh

#COPY INPUT OF FLUXACC TO /tmp
ecp -o ec:/vo3/metdat/surf_* /tmp/

#EXECUTE FLUXACC
FLUXACC

#SAVE OUTPUT OF FLUXACC IN ECFILE DIRECTORY
ecp /tmp/flux* ec:/vo3/flux/

#DELETE INPUT OF FLUXACC
erm ec:/vo3/metdat/surf_*

#CLEAN UP
rm /tmp/surf_*
rm /tmp/flux*

APPENDIX 2: Sample CONTROL file for retrieval of meteorological fields and calculation of vertical velocity component (directory job.tra)

CONTROL:

1995031415
1995031415
3
1.0
-50.0 50.0 25.0 80.0
28 29
31
319
ka7
/ws/home/ms/at/ka7/job.tra/
ec:/ka7/flux/
ec:/ka7/metdat/

* CONTROL FILE *

* 1. BEGINNING DATE [yyyymmddhh]
* 2. ENDING DATE [yyyymmddhh]
* 3. TIME INTERVAL [hr]
* 4. GRID RESOLUTION [DEG]
* 5. GRID DOMAIN - LLON/RLON/BLAT/TLAT
* 6. VERTICAL LAYERS TO EXTRACT (UVTQ and W)
* 7. NUMBER OF LAYERS AVAILABLE FROM ECMWF
* 8. SPECTRAL TRUNCATION NEEDED (<= TRUNCATION ARCHIVED)
* 9. USER NAME
* 10. DIRECTORY WHERE JOB.TRA IS INSTALLED
* 11. DIRECTORY WHERE DEACCUMULATED FLUX DATA ARE AVAILABLE
* (EC FILE SYSTEM)
* 12. DIRECTORY WHERE ENYYMMDD OUTPUT FILES ARE SAVED
* (EC FILE SYSTEM)

**********************************************************

DOMAIN:

 TLAT
 + + + + + + + + + RLON
 + + + + + + + + +
 + + + + + + + + +
 + + + + + + + + +
LLON BLAT

*****************************************************

SPECTRAL TRUNCATION/LAYERS AVAILABLE FROM ECMWF MODEL

T63 L16 Beginning: 20/04/1983 1800
T106 L16 Beginning: 01/05/1985 0000
T106 L19 Beginning: 03/04/1987 0000
T213 L31 Beginning: 17/09/1991 0000
T319 L31 Beginning: 01/04/1998 0000
T639 future

**********************************************

DEPENDENCE OF TRUNCATION TO NEEDED LAM-PHI OUTPUT


Gaussian Lat/long grid increment truncation
 ---- -------------------------------- ----
 N36 2.5 and greater T63
 N60 1.5 up to, but not including 2.5 T106
 N150 0.6 up to, but not including 1.5 T213
 N225 0.4 up to, but not including 0.6 T319
 N225 Up to, but not including 0.4 T639

Appendix 3: Plot of sample fields

Last update: 12 Jan 1999 by Gerhard Wotawa | URL of this page: http://www.boku.ac.at/imp/envmet/ecmwf_extr.html