IMP Institute of Meteorology and Physics
University of Agricultural Sciences, Vienna, Austria 
BOKU
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

    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

    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:

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

    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:

    To use these routines, please follow the instructions step by step: A sample CONTROL file is printed in Appendix 2.

    Technical comments:


    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: If you would like to retrieve data, e.g., for domain 2, you have to proceed as follows:

    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