Changeset 64cf353 in flex_extract.git for python


Ignore:
Timestamp:
Feb 8, 2018, 9:54:05 PM (6 years ago)
Author:
Anne Philipp <bscannephilipp@…>
Branches:
master, ctbto, dev
Children:
02c8c50
Parents:
6180177
Message:

pep8 changes + documentation added + minor code style changes

Location:
python
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • python/FlexpartTools.py

    rd69b677 r64cf353  
    1 #
    2 # This software is licensed under the terms of the Apache Licence Version 2.0
    3 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
    4 #
    5 # Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs
    6 #
    7 # Creation: October  2014 - Anne Fouilloux - University of Oslo
    8 # Extension November 2015 - Leopold Haimberger - University of Vienna for:
    9 # - using the WebAPI also for general MARS retrievals
    10 # - job submission on ecgate and cca
    11 # - job templates suitable for twice daily operational dissemination
    12 # - dividing retrievals of longer periods into digestable chunks
    13 # - retrieve also longer term forecasts, not only analyses and short term forecast data
    14 # - conversion into GRIB2
    15 # - conversion into .fp format for faster execution of FLEXPART
    16 #
    17 # Requirements:
    18 # in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed
    19 # ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/
    20 # dateutils
    21 # matplotlib (optional, for debugging)
    22 #
    23 #
     1#!/usr/bin/env python
     2# -*- coding: utf-8 -*-
     3#************************************************************************
     4# TODO AP
     5#AP
     6# - Functionality Provided is not correct for this file
     7# - localpythonpath should not be set in module load section!
     8# - Change History ist nicht angepasst ans File! Original geben lassen
     9# - def myerror muss angepasst werden da derzeit manuelle modifikation notwendig
     10# - def --init-- marsretrieval class, remove dataset and expect?!
     11# - the EIFlexpart class is jsut for EI ? So we should maybe create classes
     12#   for each possible data set so the boundarys of the attributes can be
     13#   validated!
     14#************************************************************************
     15"""
     16@Author: Anne Fouilloux (University of Oslo)
     17
     18@Date: October 2014
     19
     20@ChangeHistory:
     21    November 2015 - Leopold Haimberger (University of Vienna):
     22        - using the WebAPI also for general MARS retrievals
     23        - job submission on ecgate and cca
     24        - job templates suitable for twice daily operational dissemination
     25        - dividing retrievals of longer periods into digestable chunks
     26        - retrieve also longer term forecasts, not only analyses and
     27          short term forecast data
     28        - conversion into GRIB2
     29        - conversion into .fp format for faster execution of FLEXPART
     30
     31    February 2018 - Anne Philipp (University of Vienna):
     32        - applied PEP8 style guide
     33        - added documentation
     34
     35@License:
     36    (C) Copyright 2014 UIO.
     37
     38    This software is licensed under the terms of the Apache Licence Version 2.0
     39    which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
     40
     41@Requirements:
     42    - A standard python 2.6 or 2.7 installation
     43    - dateutils
     44    - matplotlib (optional, for debugging)
     45    - ECMWF specific packages, all available from https://software.ecmwf.int/
     46        ECMWF WebMARS, gribAPI with python enabled, emoslib and
     47        ecaccess web toolkit
     48
     49@Description:
     50    Further documentation may be obtained from www.flexpart.eu.
     51
     52    Functionality provided:
     53        Prepare input 3D-wind fields in hybrid coordinates +
     54        surface fields for FLEXPART runs
     55"""
     56# ------------------------------------------------------------------------------
     57# MODULES
     58# ------------------------------------------------------------------------------
    2459import subprocess
    25 from argparse import ArgumentParser,ArgumentDefaultsHelpFormatter
     60from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
    2661import traceback
    2762import shutil
    28 import os, errno,sys,inspect,glob
     63import os
     64import errno
     65import sys
     66import inspect
     67import glob
    2968import datetime
    30 #import re
    31 #import copy
    32 
    3369from string import atoi
    34 from numpy  import *
    35 
    36 #sys.stdout=open('ECMWFDATA_'+str(os.getpid()),'w')
    37 #from ecmwfapi import ECMWFDataServer
    38 ecapi=True
     70from numpy import *
     71ecapi = True
    3972try:
    4073    import ecmwfapi
    4174except ImportError:
    42     ecapi=False
     75    ecapi = False
    4376from gribapi import *
    4477from GribTools import GribTools
    4578
    46 def interpret_args_and_control(*args,**kwargs):
    47 
    48     parser = ArgumentParser(description='Retrieve FLEXPART input from ECMWF MARS archive',
     79def interpret_args_and_control():
     80    '''
     81    @Description:
     82        Assigns the command line arguments and reads control file
     83        content. Apply default values for non mentioned arguments.
     84
     85    @Input:
     86        <nothing>
     87
     88    @Return:
     89        args: instance of ArgumentParser
     90            Contains the commandline arguments from script/program call.
     91
     92        c: instance of class Control
     93            Contains all necessary information of a control file. The parameters
     94            are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     95            NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST,
     96            RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA,
     97            SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS,
     98            ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR
     99            For more information about format and content of the parameter see
     100            documentation.
     101
     102    '''
     103    parser = ArgumentParser(description='Retrieve FLEXPART input from \
     104                            ECMWF MARS archive',
    49105                            formatter_class=ArgumentDefaultsHelpFormatter)
    50106
    51 # the most important arguments 
     107# the most important arguments
    52108    parser.add_argument("--start_date", dest="start_date",
    53109                        help="start date YYYYMMDD")
    54     parser.add_argument( "--end_date", dest="end_date",
    55                          help="end_date YYYYMMDD")
    56     parser.add_argument( "--date_chunk", dest="date_chunk",default=None,
    57                          help="# of days to be retrieved at once")
    58    
     110    parser.add_argument("--end_date", dest="end_date",
     111                        help="end_date YYYYMMDD")
     112    parser.add_argument("--date_chunk", dest="date_chunk", default=None,
     113                        help="# of days to be retrieved at once")
     114
    59115# some arguments that override the default in the control file
    60     parser.add_argument("--basetime", dest="basetime", help="base such as 00/12 (for half day retrievals)")
    61     parser.add_argument("--step", dest="step", help="steps such as 00/to/48")
    62     parser.add_argument("--levelist", dest="levelist",help="Vertical levels to be retrieved, e.g. 30/to/60")
    63     parser.add_argument("--area", dest="area", help="area defined as north/west/south/east")
    64    
     116    parser.add_argument("--basetime", dest="basetime",
     117                        help="base such as 00/12 (for half day retrievals)")
     118    parser.add_argument("--step", dest="step",
     119                        help="steps such as 00/to/48")
     120    parser.add_argument("--levelist", dest="levelist",
     121                        help="Vertical levels to be retrieved, e.g. 30/to/60")
     122    parser.add_argument("--area", dest="area",
     123                        help="area defined as north/west/south/east")
     124
    65125# set the working directories
    66     parser.add_argument("--inputdir", dest="inputdir",default=None,
     126    parser.add_argument("--inputdir", dest="inputdir", default=None,
    67127                        help="root directory for storing intermediate files")
    68     parser.add_argument("--outputdir", dest="outputdir",default=None,
     128    parser.add_argument("--outputdir", dest="outputdir", default=None,
    69129                        help="root directory for storing output files")
    70130    parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
    71                         help="FLEXPART root directory (to find grib2flexpart and COMMAND file)\n\
    72                               Normally ECMWFDATA resides in the scripts directory of the FLEXPART distribution")
    73    
     131                        help="FLEXPART root directory (to find grib2flexpart \
     132                        and COMMAND file)\n\ Normally ECMWFDATA resides in \
     133                        the scripts directory of the FLEXPART distribution")
     134
    74135# this is only used by prepareFLEXPART.py to rerun a postprocessing step
    75136    parser.add_argument("--ppid", dest="ppid",
    76                         help="Specify parent process id for rerun of prepareFLEXPART")
     137                        help="Specify parent process id for \
     138                        rerun of prepareFLEXPART")
    77139
    78140# arguments for job submission to ECMWF, only needed by submit.py
    79     parser.add_argument("--job_template", dest='job_template',default="job.temp",
     141    parser.add_argument("--job_template", dest='job_template',
     142                        default="job.temp",
    80143                        help="job template file for submission to ECMWF")
    81144    #parser.add_argument("--remote", dest="remote",
    82                         #help="target for submission to ECMWF (ecgate or cca etc.)")
     145                        #help="target for submission to ECMWF \
     146                        #(ecgate or cca etc.)")
    83147    parser.add_argument("--queue", dest="queue",
    84                         help="queue for submission to ECMWF (e.g. ecgate or cca )")
    85 
    86 
    87     parser.add_argument("--controlfile", dest="controlfile",default='CONTROL.temp',
     148                        help="queue for submission to ECMWF \
     149                        (e.g. ecgate or cca )")
     150    parser.add_argument("--controlfile", dest="controlfile",
     151                        default='CONTROL.temp',
    88152                        help="file with control parameters")
    89     parser.add_argument("--debug", dest="debug",default=0,
     153    parser.add_argument("--debug", dest="debug", default=0,
    90154                        help="Debug mode - leave temporary files intact")
     155
    91156    args = parser.parse_args()
    92157
    93158    try:
    94         c=Control(args.controlfile)
     159        c = Control(args.controlfile)
    95160    except IOError:
    96         try:
    97             c=Control(localpythonpath+args.controlfile)
    98    
    99         except:
    100             print 'Could not read control file "'+args.controlfile+'"'
    101             print 'Either it does not exist or its syntax is wrong.'
    102             print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information'
    103             exit(1)
    104 
    105     if  args.start_date==None and getattr(c,'start_date')==None:
    106         print 'start_date specified neither in command line nor in control file '+args.controlfile
    107         print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information'
    108         exit(1)
    109 
    110     if args.start_date!=None:
    111         c.start_date=args.start_date 
    112     if args.end_date!=None:
    113         c.end_date=args.end_date
    114     if c.end_date==None:
    115         c.end_date=c.start_date
    116     if args.date_chunk!=None:
    117         c.date_chunk=args.date_chunk 
    118    
    119     if not hasattr(c,'debug'):
    120         c.debug=args.debug 
    121        
    122 
    123     if args.inputdir==None and args.outputdir==None:
    124         c.inputdir='../work'
    125         c.outputdir='../work'
     161        try:
     162            c = Control(localpythonpath + args.controlfile)
     163        except:
     164            print('Could not read control file "' + args.controlfile + '"')
     165            print('Either it does not exist or its syntax is wrong.')
     166            print('Try "' + sys.argv[0].split('/')[-1] +
     167                  ' -h" to print usage information')
     168            exit(1)
     169
     170    if  args.start_date is None and getattr(c, 'start_date') is None:
     171        print('start_date specified neither in command line nor \
     172               in control file ' + args.controlfile)
     173        print('Try "' + sys.argv[0].split('/')[-1] +
     174              ' -h" to print usage information')
     175        exit(1)
     176
     177    if args.start_date is not None:
     178        c.start_date = args.start_date
     179    if args.end_date is not None:
     180        c.end_date = args.end_date
     181    if c.end_date is None:
     182        c.end_date = c.start_date
     183    if args.date_chunk is not None:
     184        c.date_chunk = args.date_chunk
     185
     186    if not hasattr(c, 'debug'):
     187        c.debug = args.debug
     188
     189    if args.inputdir is None and args.outputdir is None:
     190        c.inputdir = '../work'
     191        c.outputdir = '../work'
    126192    else:
    127         if args.inputdir!=None:
    128             c.inputdir=args.inputdir
    129             if args.outputdir==None:
    130                 c.outputdir=args.inputdir
    131         if args.outputdir!=None:
    132             c.outputdir=args.outputdir
    133             if args.inputdir==None:
    134                 c.inputdir=args.outputdir
    135    
    136     if hasattr(c,'outputdir')==False and args.outputdir==None:
    137         c.outputdir=c.inputdir
     193        if args.inputdir is not None:
     194            c.inputdir = args.inputdir
     195        if args.outputdir is None:
     196            c.outputdir = args.inputdir
     197        if args.outputdir is not None:
     198            c.outputdir = args.outputdir
     199        if args.inputdir is None:
     200            c.inputdir = args.outputdir
     201
     202    if hasattr(c, 'outputdir') is False and args.outputdir is None:
     203        c.outputdir = c.inputdir
    138204    else:
    139         if args.outputdir!=None:
    140             c.outputdir=args.outputdir
    141            
    142     if args.area!=None:
    143         afloat='.' in args.area
    144         l=args.area.split('/')
    145         if afloat:
    146             for i in range(len(l)):
    147                 l[i]=str(int(float(l[i])*1000))
    148         c.upper,c.left,c.lower,c.right=l
    149 
    150 # basetime aktiviert den ''operational mode''   
    151     if args.basetime!=None:
    152         c.basetime=args.basetime
    153     #if int(c.basetime)==0:
    154         #c.start_date=datetime.datetime.strftime(datetime.datetime.strptime(c.start_date,'%Y%m%d')-datetime.timedelta(days=1),'%Y%m%d')
    155         #c.end_date=c.start_date
    156        
    157     if args.step!=None:
    158         l=args.step.split('/')
    159         if 'to' in args.step.lower():
    160             if 'by' in args.step.lower():
    161                 ilist=arange(int(l[0]),int(l[2])+1,int(l[4]))
    162                 c.step=['{:0>3}'.format(i) for i in ilist]
    163             else:
    164                 myerror(None,args.step+':\n'+'please use "by" as well if "to" is used')
    165         else:
    166             c.step=l
    167        
    168     if args.levelist!=None:
    169         c.levelist=args.levelist
    170         if 'to' in c.levelist:
    171             c.level=c.levelist.split('/')[2]
    172         else:
    173             c.level=c.levelist.split('/')[-1]
    174                    
    175    
    176     if args.flexpart_root_scripts!=None:
    177         c.flexpart_root_scripts=args.flexpart_root_scripts
    178        
    179     return args,c
     205        if args.outputdir is not None:
     206            c.outputdir = args.outputdir
     207
     208    if args.area is not None:
     209        afloat = '.' in args.area
     210        l = args.area.split('/')
     211        if afloat:
     212            for i in range(len(l)):
     213                l[i] = str(int(float(l[i]) * 1000))
     214        c.upper, c.left, c.lower, c.right = l
     215
     216# basetime aktiviert den ''operational mode''
     217    if args.basetime is not None:
     218        c.basetime = args.basetime
     219
     220    if args.step is not None:
     221        l = args.step.split('/')
     222        if 'to' in args.step.lower():
     223            if 'by' in args.step.lower():
     224                ilist = arange(int(l[0]), int(l[2]) + 1, int(l[4]))
     225                c.step = ['{:0>3}'.format(i) for i in ilist]
     226            else:
     227                myerror(None, args.step + ':\n' +
     228                        'please use "by" as well if "to" is used')
     229        else:
     230            c.step = l
     231
     232    if args.levelist is not None:
     233        c.levelist = args.levelist
     234        if 'to' in c.levelist:
     235            c.level = c.levelist.split('/')[2]
     236        else:
     237            c.level = c.levelist.split('/')[-1]
     238
     239    if args.flexpart_root_scripts is not None:
     240        c.flexpart_root_scripts = args.flexpart_root_scripts
     241
     242    return args, c
     243
    180244
    181245def install_args_and_control():
    182 
    183     parser = ArgumentParser(description='Install ECMWFDATA software locally or on ECMWF machines',
     246    '''
     247    @Description:
     248        Assigns the command line arguments for installation and reads
     249        control file content. Apply default values for non mentioned arguments.
     250
     251    @Input:
     252        <nothing>
     253
     254    @Return:
     255        args: instance of ArgumentParser
     256            Contains the commandline arguments from script/program call.
     257
     258        c: instance of class Control
     259            Contains all necessary information of a control file. The parameters
     260            are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     261            NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST,
     262            RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA,
     263            SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS,
     264            ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR
     265            For more information about format and content of the parameter see
     266            documentation.
     267    '''
     268    parser = ArgumentParser(description='Install ECMWFDATA software locally or \
     269                            on ECMWF machines',
    184270                            formatter_class=ArgumentDefaultsHelpFormatter)
    185     parser.add_argument('--target',dest='install_target',
    186                     help="Valid targets: local | ecgate | cca , the latter two are at ECMWF")
    187     parser.add_argument("--makefile", dest="makefile",help='Name of Makefile to use for compiling CONVERT2')
    188     parser.add_argument("--ecuid", dest="ecuid",help='user id at ECMWF')
    189     parser.add_argument("--ecgid", dest="ecgid",help='group id at ECMWF')
    190     parser.add_argument("--gateway", dest="gateway",help='name of local gateway server')
    191     parser.add_argument("--destination", dest="destination",help='ecaccess destination, e.g. leo@genericSftp')
     271
     272    parser.add_argument('--target', dest='install_target',
     273                        help="Valid targets: local | ecgate | cca , \
     274                        the latter two are at ECMWF")
     275    parser.add_argument("--makefile", dest="makefile",
     276                        help='Name of Makefile to use for compiling CONVERT2')
     277    parser.add_argument("--ecuid", dest="ecuid",
     278                        help='user id at ECMWF')
     279    parser.add_argument("--ecgid", dest="ecgid",
     280                        help='group id at ECMWF')
     281    parser.add_argument("--gateway", dest="gateway",
     282                        help='name of local gateway server')
     283    parser.add_argument("--destination", dest="destination",
     284                        help='ecaccess destination, e.g. leo@genericSftp')
    192285
    193286    parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
    194                         help="FLEXPART root directory on ECMWF servers (to find grib2flexpart and COMMAND file)\n\
    195                               Normally ECMWFDATA resides in the scripts directory of the FLEXPART distribution, thus the:")
    196    
     287                        help="FLEXPART root directory on ECMWF servers \
     288                        (to find grib2flexpart and COMMAND file)\n\
     289                        Normally ECMWFDATA resides in the scripts directory \
     290                        of the FLEXPART distribution, thus the:")
     291
    197292# arguments for job submission to ECMWF, only needed by submit.py
    198     parser.add_argument("--job_template", dest='job_template',default="job.temp.o",
     293    parser.add_argument("--job_template", dest='job_template',
     294                        default="job.temp.o",
    199295                        help="job template file for submission to ECMWF")
    200296
    201     parser.add_argument("--controlfile", dest="controlfile",default='CONTROL.temp',
     297    parser.add_argument("--controlfile", dest="controlfile",
     298                        default='CONTROL.temp',
    202299                        help="file with control parameters")
     300
    203301    args = parser.parse_args()
    204302
    205303    try:
    206 #    if True:
    207         c=Control(args.controlfile)
    208 #    else:
     304        c = Control(args.controlfile)
    209305    except:
    210         print 'Could not read control file "'+args.controlfile+'"'
    211         print 'Either it does not exist or its syntax is wrong.'
    212         print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information'
    213         exit(1)
    214 
    215     if args.install_target!='local':
    216         if args.ecgid==None or args.ecuid==None or args.gateway==None or args.destination==None :
    217             print 'Please enter your ECMWF user id and group id as well as the \nname of the local gateway and the ectrans destination '
    218             print 'with command line options --ecuid --ecgid --gateway --destination'
    219             print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information'
    220             print 'Please consult ecaccess documentation or ECMWF user support for further details'
    221             sys.exit(1)
    222         else:
    223             c.ecuid=args.ecuid
    224             c.ecgid=args.ecgid
    225             c.gateway=args.gateway
    226             c.destination=args.destination
    227            
     306        print('Could not read control file "' + args.controlfile + '"')
     307        print('Either it does not exist or its syntax is wrong.')
     308        print('Try "' + sys.argv[0].split('/')[-1] +
     309              ' -h" to print usage information')
     310        exit(1)
     311
     312    if args.install_target != 'local':
     313        if (args.ecgid is None or args.ecuid is None or args.gateway is None
     314            or args.destination is None):
     315            print('Please enter your ECMWF user id and group id as well as \
     316                   the \nname of the local gateway and the ectrans \
     317                   destination ')
     318            print('with command line options --ecuid --ecgid \
     319                   --gateway --destination')
     320            print('Try "' + sys.argv[0].split('/')[-1] +
     321                  ' -h" to print usage information')
     322            print('Please consult ecaccess documentation or ECMWF user support \
     323                   for further details')
     324            sys.exit(1)
     325        else:
     326            c.ecuid = args.ecuid
     327            c.ecgid = args.ecgid
     328            c.gateway = args.gateway
     329            c.destination = args.destination
     330
    228331    try:
    229         c.makefile=args.makefile
     332        c.makefile = args.makefile
    230333    except:
    231         pass
    232    
    233     if args.install_target=='local':
    234         if args.flexpart_root_scripts==None:
    235             c.flexpart_root_scripts='../'
    236         else:
    237             c.flexpart_root_scripts=args.flexpart_root_scripts
    238        
    239     if args.install_target!='local':
    240         if args.flexpart_root_scripts==None:
    241             c.ec_flexpart_root_scripts='${HOME}'
    242         else:
    243             c.ec_flexpart_root_scripts=args.flexpart_root_scripts
    244        
    245     return args,c
     334        pass
     335
     336    if args.install_target == 'local':
     337        if args.flexpart_root_scripts is None:
     338            c.flexpart_root_scripts = '../'
     339        else:
     340            c.flexpart_root_scripts = args.flexpart_root_scripts
     341
     342    if args.install_target != 'local':
     343        if args.flexpart_root_scripts is None:
     344            c.ec_flexpart_root_scripts = '${HOME}'
     345        else:
     346            c.ec_flexpart_root_scripts = args.flexpart_root_scripts
     347
     348    return args, c
     349
    246350
    247351def cleanup(c):
    248     print "cleanup"
    249     cleanlist=glob.glob(c.inputdir+"/*")
    250 #    cleanlist[-1:-1]=glob.glob(c.inputdir+"/flux*")
    251 #    cleanlist[-1:-1]=glob.glob(c.inputdir+"/*.grb")
    252 #    if c.ecapi==False and (c.ectrans=='1' or c.ecstorage=='1'):
    253 #       cleanlist[-1:-1]=glob.glob(c.inputdir+"/"+c.prefix+"*")
    254 
    255    
    256 #    silentremove(c.inputdir+"/VERTICAL.EC")
    257 #    silentremove(c.inputdir+"/date_time_stepRange.idx")
     352    '''
     353    @Description:
     354        Remove all files from intermediate directory
     355        (inputdir from control file).
     356
     357    @Input:
     358        c: instance of class Control
     359            Contains all the parameters of control files, which are:
     360            DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     361            NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL,
     362            LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA,
     363            ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX,
     364            ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL,
     365            GRIB2FLEXPART, FLEXPARTDIR
     366            For more information about format and content of the parameter
     367            see documentation.
     368
     369    @Return:
     370        <nothing>
     371    '''
     372
     373    print("cleanup")
     374
     375    cleanlist = glob.glob(c.inputdir + "/*")
    258376    for cl in cleanlist:
    259         if c.prefix not in cl:
    260             silentremove(cl)
    261         if c.ecapi==False and (c.ectrans=='1' or c.ecstorage=='1'):
    262             silentremove(cl)
    263        
    264        
    265        
    266     print "Done"
    267 
    268 
    269 def myerror(c,message='ERROR'):
    270    
    271     try:
    272         target=c.mailfail
    273     except AttributeError:
    274         target=os.getenv('USER')
    275        
    276     if(type(target) is not list):
    277         target=[target]
    278     print message
     377        if c.prefix not in cl:
     378            silentremove(cl)
     379        if c.ecapi is False and (c.ectrans == '1' or c.ecstorage == '1'):
     380            silentremove(cl)
     381
     382    print("Done")
     383
     384    return
     385
     386
     387def myerror(c, message='ERROR'):
     388    '''
     389    @Description:
     390        Print error message.
     391
     392    @Input:
     393        c: instance of class Control
     394            Contains all the parameters of control files, which are:
     395            DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     396            NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL,
     397            LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA,
     398            ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX,
     399            ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL,
     400            GRIB2FLEXPART, FLEXPARTDIR
     401            For more information about format and content of the parameter
     402            see documentation.
     403
     404        message: string
     405            Error message. Default value is "ERROR".
     406
     407    @Return:
     408        <nothing>
     409    '''
     410    # uncomment if user wants email notification directly from python
     411    #try:
     412        #target = c.mailfail
     413    #except AttributeError:
     414        #target = os.getenv('USER')
     415
     416    #if(type(target) is not list):
     417        #target = [target]
     418
     419    print(message)
     420
    279421    # uncomment if user wants email notification directly from python
    280422    #for t in target:
    281         #p=subprocess.Popen(['mail','-s ECMWFDATA v7.0 ERROR', os.path.expandvars(t)], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,bufsize=1)
    282         #tr='\n'.join(traceback.format_stack())
    283         #pout=p.communicate(input=message+'\n\n'+tr)[0]
    284         #print 'Email sent to '+os.path.expandvars(t) # +' '+pout.decode()
    285    
     423    #p = subprocess.Popen(['mail','-s ECMWFDATA v7.0 ERROR', os.path.expandvars(t)],
     424    #                     stdin = subprocess.PIPE, stdout = subprocess.PIPE,
     425    #                     stderr = subprocess.PIPE, bufsize = 1)
     426    #tr = '\n'.join(traceback.format_stack())
     427    #pout = p.communicate(input = message+'\n\n'+tr)[0]
     428    #print 'Email sent to '+os.path.expandvars(t) # +' '+pout.decode()
     429
    286430    exit(1)
    287    
    288 def normalexit(c,message='Done!'):
    289    
    290     try:
    291         target=c.mailops
    292            
    293    
    294         if(type(target) is not list):
    295             target=[target]
    296         # Uncomment if user wants notification directly from python
    297         #for t in target:
    298             #p=subprocess.Popen(['mail','-s ECMWFDATA v7.0 normal exit', os.path.expandvars(t)], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,bufsize=1)
    299             #pout=p.communicate(input=message+'\n\n')[0]
    300             #print pout.decode()
    301                
    302     except:
    303         pass
    304    
    305     print message
    306431
    307432    return
    308433
     434
     435def normalexit(c, message='Done!'):
     436    '''
     437    @Description:
     438        Print exit message.
     439
     440    @Input:
     441        c: instance of class Control
     442            Contains all the parameters of control files, which are:
     443            DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     444            NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL,
     445            LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA,
     446            ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX,
     447            ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL,
     448            GRIB2FLEXPART, FLEXPARTDIR
     449            For more information about format and content of the parameter
     450            see documentation.
     451
     452        message: string
     453            Message for exiting program. Default value is "Done!".
     454
     455    @Return:
     456        <nothing>
     457
     458    '''
     459    # Uncomment if user wants notification directly from python
     460    #try:
     461        #target = c.mailops
     462        #if(type(target) is not list):
     463            #target = [target]
     464        #for t in target:
     465            #p = subprocess.Popen(['mail','-s ECMWFDATA v7.0 normal exit',
     466            #                      os.path.expandvars(t)],
     467            #                      stdin = subprocess.PIPE,
     468            #                      stdout = subprocess.PIPE,
     469            #                      stderr = subprocess.PIPE, bufsize = 1)
     470            #pout = p.communicate(input = message+'\n\n')[0]
     471            #print pout.decode()
     472    #except:
     473        #pass
     474
     475    print(message)
     476
     477    return
     478
     479
    309480def product(*args, **kwds):
     481    '''
     482    @Description:
     483
     484    @Input:
     485
     486    @Return:
     487
     488    '''
    310489    # product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
    311     # product(range(2), repeat=3) --> 000 001 010 011 100 101 110 111
     490    # product(range(2), repeat = 3) --> 000 001 010 011 100 101 110 111
    312491    pools = map(tuple, args) * kwds.get('repeat', 1)
    313492    result = [[]]
    314493    for pool in pools:
    315         result = [x+[y] for x in result for y in pool]
     494        result = [x + [y] for x in result for y in pool]
    316495    for prod in result:
    317496        yield tuple(prod)
    318497
    319 ###############################################################
    320 # utility to remove a file if it exists
    321 # it does not fail if the file does not exist
    322 ###############################################################
     498    return
     499
     500
    323501def silentremove(filename):
     502    '''
     503    @Description:
     504        Removes the file which name is passed to the function if
     505        it exists. The function does not fail if the file does not
     506        exist.
     507
     508    @Input:
     509        filename: string
     510            The name of the file to be removed without notification.
     511
     512    @Return:
     513        <nothing>
     514    '''
    324515    try:
    325516        os.remove(filename)
    326     except OSError as e: # this would be "except OSError, e:" before Python 2.6
    327         if e.errno != errno.ENOENT: # errno.ENOENT = no such file or directory
    328             raise # re-raise exception if a different error occured
     517    except OSError as e:
     518        # this would be "except OSError, e:" before Python 2.6
     519        if e.errno is not  errno.ENOENT:
     520            # errno.ENOENT  =  no such file or directory
     521            raise  # re-raise exception if a different error occured
     522
     523    return
    329524
    330525
    331526def init128(fn):
    332    
    333     table128=dict()
     527    '''
     528    @Description:
     529        Opens and reads the grib table 128 file.
     530
     531    @Input:
     532        fn: string
     533            Path to file of ECMWF grib table number 128.
     534
     535    @Return:
     536        table128: dictionary
     537            Contains the ECMWF grib table 128 information.
     538    '''
     539    table128 = dict()
    334540    with open(fn) as f:
    335         fdata = f.read().split('\n')
    336         for data in fdata:
    337             if data[0]!='!':
    338                 table128[data[0:3]]=data[59:64].strip()
    339                
     541        fdata = f.read().split('\n')
     542    for data in fdata:
     543        if data[0] != '!':
     544            table128[data[0:3]] = data[59:64].strip()
     545
    340546    return table128
    341547
    342 def toparamId(pars,table):
    343    
    344     cpar=pars.upper().split('/')
    345     ipar=[]
     548
     549def toparamId(pars, table):
     550    '''
     551    @Description:
     552        Transform parameter names to parameter ids
     553        with ECMWF grib table 128.
     554
     555    @Input:
     556        pars: string
     557            Addpar argument from control file in the format of
     558            parameter names instead of ids.
     559
     560        table: dictionary
     561            Contains the ECMWF grib table 128 information.
     562
     563    @Return:
     564        ipar: list of integer
     565            List of addpar parameters from control file transformed to
     566            parameter ids in the format of integer.
     567    '''
     568    cpar = pars.upper().split('/')
     569    ipar = []
    346570    for par in cpar:
    347         found=False
    348         for k,v in table.iteritems():
    349             if par==k or par==v:
    350                 ipar.append(int(k))
    351                 found=True
    352                 break
    353         if found==False:
    354             print 'Warning: par '+par+' not found in table 128'
     571        found = False
     572        for k, v in table.iteritems():
     573            if par == k or par == v:
     574                ipar.append(int(k))
     575                found = True
     576                break
     577        if found is False:
     578            print('Warning: par ' + par + ' not found in table 128')
     579
    355580    return ipar
    356581
     582
    357583def dapoly(alist):
    358 
    359    
    360     pya=(alist[3]-alist[0]+3.*(alist[1]-alist[2]))/6.
    361     pyb=(alist[2]+alist[0])/2.-alist[1]-9.*pya/2.
    362     pyc=alist[1]-alist[0]-7.*pya/2.-2.*pyb
    363     pyd=alist[0]-pya/4.-pyb/3.-pyc/2.
    364     sfeld=8.*pya+4.*pyb+2.*pyc+pyd
    365    
     584    '''
     585    @Description:
     586
     587    @Input:
     588
     589    @Return:
     590
     591    '''
     592    pya = (alist[3] - alist[0] + 3. * (alist[1] - alist[2])) / 6.
     593    pyb = (alist[2] + alist[0]) / 2. - alist[1] - 9. * pya / 2.
     594    pyc = alist[1] - alist[0] - 7. * pya / 2. - 2. * pyb
     595    pyd = alist[0] - pya / 4. - pyb / 3. - pyc / 2.
     596    sfeld = 8. * pya + 4. * pyb + 2. * pyc + pyd
     597
    366598    return sfeld
    367599
     600
    368601def darain(alist):
    369     xa=alist[0]
    370     xb=alist[1]
    371     xc=alist[2]
    372     xd=alist[3]
    373     xa[xa<0.]=0.
    374     xb[xb<0.]=0.
    375     xc[xc<0.]=0.
    376     xd[xd<0.]=0.
    377 
    378     xac=0.5*xb
    379     mask=xa+xc>0.
    380     xac[mask]=xb[mask]*xc[mask]/(xa[mask]+xc[mask])
    381     xbd=0.5*xc
    382     mask=xb+xd>0.
    383     xbd[mask]=xb[mask]*xc[mask]/(xb[mask]+xd[mask])
    384     sfeld=xac+xbd
    385    
     602    '''
     603    @Description:
     604        Disaggregate a list of 4 precipitation values to generate a new value
     605        for the second position of the 4 value list. This is used for
     606        precipitation fields.
     607
     608    @Input:
     609        alist: list of size 4, float
     610            List of 4 precipitation values.
     611
     612    @Return:
     613        sfeld: float
     614            New value for the second position of the disaggregated
     615            precipitation field.
     616    '''
     617    xa = alist[0]
     618    xb = alist[1]
     619    xc = alist[2]
     620    xd = alist[3]
     621    xa[xa < 0.] = 0.
     622    xb[xb < 0.] = 0.
     623    xc[xc < 0.] = 0.
     624    xd[xd < 0.] = 0.
     625
     626    xac = 0.5 * xb
     627    mask = xa + xc > 0.
     628    xac[mask] = xb[mask] * xc[mask] / (xa[mask] + xc[mask])
     629    xbd = 0.5 * xc
     630    mask = xb + xd > 0.
     631    xbd[mask] = xb[mask] * xc[mask] / (xb[mask] + xd[mask])
     632    sfeld = xac + xbd
     633
    386634    return sfeld
    387    
     635
     636
    388637class Control:
    389     'class containing the information of the ECMWFDATA control file'
    390    
    391     def __init__(self,filename):
    392        
    393         #try:
    394             with open(filename) as f:
    395                 fdata = f.read().split('\n')
    396                 for ldata in fdata:
    397                     data=ldata.split()
    398                     if len(data)>1:
    399                         if 'm_' in data[0].lower():
    400                             data[0]=data[0][2:]
    401                         if data[0].lower()=='class':
    402                             data[0]='marsclass'
    403                         if data[0].lower()=='day1':
    404                             data[0]='start_date'
    405                         if data[0].lower()=='day2':
    406                             data[0]='end_date'
    407                         if data[0].lower()=='addpar':
    408                             if '/' in data[1]:
    409                                 if data[1][0]=='/':
    410                                     data[1]=data[1][1:]
    411                                 dd=data[1].split('/')
    412                                 data=[data[0]]
    413                                 for d in dd:
    414                                     data.append(d)
    415                             pass
    416                         if len(data)==2:
    417                             if '$' in data[1]:
    418 #                               print data[1]
    419                                 setattr(self,data[0].lower(),data[1])                           
    420                                 while '$' in data[1]:
    421                                     i=data[1].index('$')
    422                                     j=data[1].find('{')
    423                                     k=data[1].find('}')
    424                                     var=os.getenv(data[1][j+1:k])
    425                                     if var is not None:
    426                                         data[1]=data[1][:i]+var+data[1][k+1:]
    427                                     else:
    428                                         myerror(None,'Could not find variable '+data[1][j+1:k]+' while reading '+filename)
    429                                        
    430 #                                   print data[0],data[1]
    431                                 setattr(self,data[0].lower()+'_expanded',data[1])                               
    432                             else:
    433                                 if data[1].lower()!='none':
    434                                     setattr(self,data[0].lower(),data[1])
    435                                 else:
    436                                     setattr(self,data[0].lower(),None)
    437                                    
    438                         elif len(data)>2:
    439                             setattr(self,data[0].lower(),(data[1:]))
    440                     else:
    441                         pass
    442                    
    443                 if not hasattr(self,'start_date'):
    444                     self.start_date=None
    445                 if not hasattr(self,'end_date'):
    446                     self.end_date=self.start_date
    447                 if not hasattr(self,'accuracy'):
    448                     self.accuracy=24
    449                 if not hasattr(self,'omega'):
    450                     self.omega='0'
    451                 if not hasattr(self,'cwc'):
    452                     self.cwc='0'
    453                 if not hasattr(self,'omegadiff'):
    454                     self.omegadiff='0'
    455                 if not hasattr(self,'etadiff'):
    456                     self.etadiff='0'
    457                 if not hasattr(self,'levelist'):
    458                     if not hasattr(self,'level'):
    459                         print 'Warning: neither levelist nor level specified in CONTROL file'   
    460                     else:
    461                         self.levelist='1/to/'+self.level
    462                 else:
    463                     if 'to' in self.levelist:
    464                         self.level=self.levelist.split('/')[2]
    465                     else:
    466                         self.level=self.levelist.split('/')[-1]
    467                    
    468                 if not hasattr(self,'maxstep'):
    469                     self.maxstep=0
    470                     for s in self.step:
    471                         if int(s)>self.maxstep:
    472                             self.maxstep=int(s)
    473                 else:
    474                     self.maxstep=int(self.maxstep)
    475                 if not hasattr(self,'prefix'):
    476                     self.prefix='EN'
    477                 if not hasattr(self,'makefile'):
    478                     self.makefile=None
    479                 if not hasattr(self,'basetime'):
    480                     self.basetime=None
    481                 if not hasattr(self,'date_chunk'):
    482                     self.date_chunk='3'
    483                    
    484                 self.ecmwfdatadir=os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))+'/../' # script directory
    485 #               if not hasattr(self,'exedir'):
    486                 self.exedir=self.ecmwfdatadir+'src/'
    487                 if not hasattr(self,'grib2flexpart'):
    488                     self.grib2flexpart='0'
    489                 if not hasattr(self,'flexpart_root_scripts'):
    490                     self.flexpart_root_scripts=self.ecmwfdatadir
    491                
    492            
    493 
    494             return
     638    '''
     639    Class containing the information of the ECMWFDATA control file
     640
     641    Contains all the parameters of control files, which are:
     642    DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     643    NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL,
     644    LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA,
     645    ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX,
     646    ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL,
     647    GRIB2FLEXPART, FLEXPARTDIR
     648    For more information about format and content of the parameter
     649    see documentation.
     650    '''
     651
     652    def __init__(self, filename):
     653        '''
     654        @Description:
     655            Initialises the instance of Control class and defines and
     656            assign all controlfile variables.
     657
     658        @Input:
     659            self: instance of Control class
     660                Description see class documentation.
     661
     662            filename: string
     663                Name of control file.
     664
     665        @Return:
     666            <nothing>
     667        '''
     668        with open(filename) as f:
     669            fdata = f.read().split('\n')
     670        for ldata in fdata:
     671            data = ldata.split()
     672            if len(data) > 1:
     673                if 'm_' in data[0].lower():
     674                    data[0] = data[0][2:]
     675                if data[0].lower() == 'class':
     676                    data[0] = 'marsclass'
     677                if data[0].lower() == 'day1':
     678                    data[0] = 'start_date'
     679                if data[0].lower() == 'day2':
     680                    data[0] = 'end_date'
     681                if data[0].lower() == 'addpar':
     682                    if '/' in data[1]:
     683                        # remove leading '/' sign from addpar content
     684                        if data[1][0] == '/':
     685                            data[1] = data[1][1:]
     686                        dd = data[1].split('/')
     687                        data = [data[0]]
     688                        for d in dd:
     689                            data.append(d)
     690                    pass
     691                if len(data) == 2:
     692                    if '$' in data[1]:
     693                        setattr(self, data[0].lower(), data[1])
     694                        while '$' in data[1]:
     695                            i = data[1].index('$')
     696                            j = data[1].find('{')
     697                            k = data[1].find('}')
     698                            var = os.getenv(data[1][j+1:k])
     699                            if var is not None:
     700                                data[1] = data[1][:i] + var + data[1][k+1:]
     701                            else:
     702                                myerror(None, 'Could not find variable ' +
     703                                        data[1][j+1:k] + ' while reading ' +
     704                                        filename)
     705                        setattr(self, data[0].lower()+'_expanded', data[1])
     706                    else:
     707                        if data[1].lower() != 'none':
     708                            setattr(self, data[0].lower(), data[1])
     709                        else:
     710                            setattr(self, data[0].lower(), None)
     711                elif len(data) > 2:
     712                    setattr(self, data[0].lower(), (data[1:]))
     713            else:
     714                pass
     715
     716        if not hasattr(self, 'start_date'):
     717            self.start_date = None
     718        if not hasattr(self, 'end_date'):
     719            self.end_date = self.start_date
     720        if not hasattr(self, 'accuracy'):
     721            self.accuracy = 24
     722        if not hasattr(self, 'omega'):
     723            self.omega = '0'
     724        if not hasattr(self, 'cwc'):
     725            self.cwc = '0'
     726        if not hasattr(self, 'omegadiff'):
     727            self.omegadiff = '0'
     728        if not hasattr(self, 'etadiff'):
     729            self.etadiff = '0'
     730        if not hasattr(self, 'levelist'):
     731            if not hasattr(self, 'level'):
     732                print(('Warning: neither levelist nor level \
     733                       specified in CONTROL file'))
     734            else:
     735                self.levelist = '1/to/' + self.level
     736        else:
     737            if 'to' in self.levelist:
     738                self.level = self.levelist.split('/')[2]
     739            else:
     740                self.level = self.levelist.split('/')[-1]
     741
     742        if not hasattr(self, 'maxstep'):
     743            self.maxstep = 0
     744            for s in self.step:
     745                if int(s) > self.maxstep:
     746                    self.maxstep = int(s)
     747        else:
     748            self.maxstep = int(self.maxstep)
     749
     750        if not hasattr(self, 'prefix'):
     751            self.prefix = 'EN'
     752        if not hasattr(self, 'makefile'):
     753            self.makefile = None
     754        if not hasattr(self, 'basetime'):
     755            self.basetime = None
     756        if not hasattr(self, 'date_chunk'):
     757            self.date_chunk = '3'
     758
     759        self.ecmwfdatadir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))+'/../' # script directory
     760#AP wieso nicht abfragen? Wieso immer setzen?
     761#        if not hasattr(self,'exedir'):
     762        self.exedir = self.ecmwfdatadir + 'src/'
     763        if not hasattr(self, 'grib2flexpart'):
     764            self.grib2flexpart = '0'
     765        if not hasattr(self, 'flexpart_root_scripts'):
     766            self.flexpart_root_scripts = self.ecmwfdatadir
     767
     768        return
     769
    495770    def __str__(self):
    496        
    497         attrs = vars(self)
    498         # {'kids': 0, 'name': 'Dog', 'color': 'Spotted', 'age': 10, 'legs': 2, 'smell': 'Alot'}
    499         # now dump this in some way or another
    500         return ', '.join("%s: %s" % item for item in attrs.items())
    501    
     771        '''
     772        @Description:
     773            Prepares a single string with all the comma seperated Control
     774            class attributes including their values. The string has the format
     775#AP            ????????????????? anschaun
     776
     777        @Input:
     778            self: instance of Control class
     779                Description see class documentation.
     780
     781        @Return:
     782            string of Control class attributes with their values
     783        '''
     784        attrs = vars(self)
     785        # {'kids': 0, 'name': 'Dog', 'color': 'Spotted', 'age': 10, 'legs': 2, 'smell': 'Alot'}
     786        # now dump this in some way or another
     787        return ', '.join("%s: %s" % item for item in attrs.items())
     788
    502789    def tolist(self):
    503         attrs=vars(self)
    504         l=list()
    505         for item in attrs.items():
    506             if '_expanded' in item[0]:
    507                 pass
    508             elif 'exedir' in item[0]:
    509                 pass
    510             elif 'flexpart_root_scripts' in item[0]:
    511                 pass
    512             elif 'ecmwfdatadir' in item[0]:
    513                 pass
    514            
    515             else:
    516                 if type(item[1]) is list:
    517                     stot=''
    518                     for s in item[1]:
    519                         stot+=s+' '
    520                        
    521                     l.append("%s %s" % (item[0],stot))
    522                 else:
    523                     l.append("%s %s" % item )
    524         return sorted(l)
    525    
    526    
    527        
    528 ##############################################################
    529 # MARSretrieval class
    530 ##############################################################
     790        '''
     791        @Description:
     792            Just generates a list of strings containing the attributes and
     793            assigned values except the attributes "_expanded", "exedir",
     794            "ecmwfdatadir" and "flexpart_root_scripts".
     795
     796        @Input:
     797            self: instance of Control class
     798                Description see class documentation.
     799
     800        @Return:
     801            l: list
     802                A sorted list of the all Control class attributes with
     803                their values except the attributes "_expanded", "exedir",
     804                "ecmwfdatadir" and "flexpart_root_scripts".
     805        '''
     806        attrs = vars(self)
     807        l = list()
     808        for item in attrs.items():
     809            if '_expanded' in item[0]:
     810                pass
     811            elif 'exedir' in item[0]:
     812                pass
     813            elif 'flexpart_root_scripts' in item[0]:
     814                pass
     815            elif 'ecmwfdatadir' in item[0]:
     816                pass
     817            else:
     818                if type(item[1]) is list:
     819                    stot = ''
     820                    for s in item[1]:
     821                        stot += s + ' '
     822
     823                    l.append("%s %s" % (item[0], stot))
     824                else:
     825#AP syntax error with doubled %s ???
     826                    l.append("%s %s" % item )
     827        return sorted(l)
     828
     829
    531830class MARSretrieval:
    532831    'class for MARS retrievals'
    533832
    534     def __init__(self,server, marsclass="ei",type="",levtype="",levelist="",
    535                  repres="", date="",resol="",stream="",area="",time="",step="",expver="1",
    536                  number="",accuracy="", grid="", gaussian="", target="",param="",dataset="",expect=""):
    537 #        self.dataset=dataset
    538         self.marsclass=marsclass
    539         self.type=type
    540         self.levtype=levtype
    541         self.levelist=levelist
    542         self.repres=repres
    543         self.date=date
    544         self.resol=resol
    545         self.stream=stream
    546         self.area=area
    547         self.time=time
    548         self.step=step
    549         self.expver=expver
    550         self.target=target
    551         self.param=param
    552         self.number=number
    553         self.accuracy=accuracy
    554         self.grid=grid
    555         self.gaussian=gaussian
    556 #       self.expect=expect
    557         self.server=server
     833    def __init__(self, server, marsclass = "ei", type = "", levtype = "",
     834                 levelist = "", repres = "", date = "", resol = "", stream = "",
     835                 area = "", time = "", step = "", expver = "1", number = "",
     836                 accuracy = "", grid = "", gaussian = "", target = "",
     837                 param = "", dataset = "", expect = ""):
     838        '''
     839        @Description:
     840
     841        @Input:
     842            self
     843            server
     844            marsclass = "ei"
     845            type = ""
     846            levtype = ""
     847            levelist = ""
     848            repres = ""
     849            date = ""
     850            resol = ""
     851            stream = ""
     852            area = ""
     853            time = ""
     854            step = ""
     855            expver = "1"
     856            number = ""
     857            accuracy = ""
     858            grid = ""
     859            gaussian = ""
     860            target = ""
     861            param = ""
     862            dataset = ""
     863            expect = ""
     864
     865        @Return:
     866            <nothing>
     867        '''
     868
     869#        self.dataset = dataset
     870        self.marsclass = marsclass
     871        self.type = type
     872        self.levtype = levtype
     873        self.levelist = levelist
     874        self.repres = repres
     875        self.date = date
     876        self.resol = resol
     877        self.stream = stream
     878        self.area = area
     879        self.time = time
     880        self.step = step
     881        self.expver = expver
     882        self.target = target
     883        self.param = param
     884        self.number = number
     885        self.accuracy = accuracy
     886        self.grid = grid
     887        self.gaussian = gaussian
     888#    self.expect = expect
     889        self.server = server
     890
     891        return
    558892
    559893    def displayInfo(self):
    560         attrs=vars(self)
    561         for item in attrs.items():
    562             if item[0] in ('server'):
    563                 pass
    564             else:
    565                 print item[0]+': '+str(item[1])
    566        
    567         return
     894        '''
     895        @Description:
     896            .
     897
     898        @Input:
     899            self:
     900
     901        @Return:
     902            <nothing>
     903        '''
     904        attrs = vars(self)
     905        for item in attrs.items():
     906            if item[0] in ('server'):
     907                pass
     908            else:
     909                print(item[0] + ': ' + str(item[1]))
     910
     911        return
    568912
    569913    def dataRetrieve(self):
    570         attrs=vars(self)
    571 
    572 #        self.server.retrieve(dicolist)
    573         s='ret'
    574         for k,v in attrs.iteritems():
    575             if k in ('server'):
    576                 continue
    577             if k=='marsclass':
    578                 k='class'
    579             if v=='':
    580                 continue
    581             if k.lower()=='target':
    582                 target=v
    583             else:
    584                 s=s+','+k+'='+str(v)
    585 
    586         if self.server != False:
    587             try:
    588                 self.server.execute(s,target)
    589             except:
    590                 print 'MARS Request failed, have you already registered at apps.ecmwf.int?'
    591                 raise IOError
    592             if os.stat(target).st_size==0:
    593                 print 'MARS Request returned no data - please check request'
    594                 raise IOError
     914        '''
     915        @Description:
     916            .
     917
     918        @Input:
     919            self:
     920
     921        @Return:
     922            <nothing>
     923        '''
     924        attrs = vars(self)
     925    #   self.server.retrieve(dicolist)
     926        s = 'ret'
     927        for k, v in attrs.iteritems():
     928            if k in ('server'):
     929                continue
     930            if k == 'marsclass':
     931                k = 'class'
     932            if v == '':
     933                continue
     934            if k.lower() == 'target':
     935                target = v
     936            else:
     937                s = s + ',' + k + '=' + str(v)
     938
     939        if self.server !=  False:
     940            try:
     941                self.server.execute(s,target)
     942            except:
     943                print('MARS Request failed, have you already registered at apps.ecmwf.int?')
     944                raise IOError
     945            if os.stat(target).st_size == 0:
     946                print('MARS Request returned no data - please check request')
     947                raise IOError
    595948        else:
    596            
    597             s+=',target="'+target+'"'
    598             p=subprocess.Popen(['mars'], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,bufsize=1)
    599             pout=p.communicate(input=s)[0]
    600             print pout.decode()
    601             if 'Some errors reported' in pout.decode():
    602                 print 'MARS Request failed - please check request'
    603                 raise IOError
    604            
    605             if os.stat(target).st_size==0:
    606                 print 'MARS Request returned no data - please check request'
    607                 raise IOError
    608            
    609         return
    610            
     949            s += ',target = "' + target + '"'
     950            p = subprocess.Popen(['mars'], stdin=subprocess.PIPE,
     951                                 stdout=subprocess.PIPE,
     952                                 stderr=subprocess.PIPE, bufsize = 1)
     953            pout = p.communicate(input=s)[0]
     954            print(pout.decode())
     955            if 'Some errors reported' in pout.decode():
     956                print('MARS Request failed - please check request')
     957                raise IOError
     958
     959            if os.stat(target).st_size == 0:
     960                print('MARS Request returned no data - please check request')
     961                raise IOError
     962
     963        return
     964
     965
     966##############################################################
    611967##############################################################
    612968class EIFlexpart:
    613     'class to retrieve Era Interim data for running FLEXPART'
    614 ##############################################################
    615 
    616     def __init__(self,c,fluxes=False):
    617 # different mars types for retrieving reanalysis data for flexpart
    618        
    619         self.types=dict()
    620         try:
    621             if c.maxstep>len(c.type):    # Pure forecast mode
    622                 c.type=[c.type[1]]
    623                 c.step=['{:0>3}'.format(int(c.step[0]))]
    624                 c.time=[c.time[0]]
    625                 for i in range(1,c.maxstep+1):
    626                     c.type.append(c.type[0])
    627                     c.step.append('{:0>3}'.format(i))
    628                     c.time.append(c.time[0])
    629         except:
    630             pass
    631        
    632         self.inputdir=c.inputdir
    633         self.basetime=c.basetime
    634         self.dtime=c.dtime
    635         self.mars={}                       
    636         i=0
    637         j=0
    638         if fluxes==True and c.maxstep<24:
    639             self.types[c.type[1]]={'times':'00/12','steps':'{}/to/12/by/{}'.format(c.dtime,c.dtime)}
    640             i=1
    641             for k in [0,12]:
    642                 for j in range(int(c.dtime),13,int(c.dtime)):
    643                     self.mars['{:0>3}'.format(i*int(c.dtime))]=[c.type[1],'{:0>2}'.format(k),'{:0>3}'.format(j)]
    644                     i+=1
    645         else:
    646             for ty,st,ti in zip(c.type,c.step,c.time):
    647                 btlist=range(24)
    648                 if c.basetime=='12':
    649                     btlist=[1,2,3,4,5,6,7,8,9,10,11,12]
    650 #                   btday =[0,0,0,0,0,0,0,0,0, 0, 0, 0]
    651                 if c.basetime=='00':
    652                     btlist=[13,14,15,16,17,18,19,20,21,22,23,0]
    653 #                   btday =[-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1]
    654 #                   if c.step[0]=='00':
    655 #                       btday[-1]=0
    656                
    657                 if mod(i,int(c.dtime))==0 and (i in btlist or c.maxstep>24):
    658                     if ty not in self.types.keys():
    659                         self.types[ty]={'times':'','steps':''}
    660                     if ti not in self.types[ty]['times']:
    661                         if len(self.types[ty]['times'])>0:
    662                             self.types[ty]['times']+='/'
    663                         self.types[ty]['times']+=ti
    664                     if st not in self.types[ty]['steps']:
    665                         if len(self.types[ty]['steps'])>0:
    666                             self.types[ty]['steps']+='/'
    667                         self.types[ty]['steps']+=st
    668                    
    669                     self.mars['{:0>3}'.format(j)]=[ty,'{:0>2}'.format(int(ti)),'{:0>3}'.format(int(st))]
    670                     j+=int(c.dtime)
    671                
    672                 i+=1
    673                
    674 # Different grids need different retrievals     
    675 # SH=Spherical Harmonics, GG=Gaussian Grid, OG=Output Grid, ML=MultiLevel, SL=SingleLevel
    676         self.params={'SH__ML':'','SH__SL':'',
    677                      'GG__ML':'','GG__SL':'',
    678                      'OG__ML':'','OG__SL':'',
    679                      'OG_OROLSM_SL':'','OG_acc_SL':''}
    680        
    681         self.marsclass=c.marsclass
    682         self.stream=c.stream
    683         self.number=c.number
    684         self.resol=c.resol
    685         if 'N' in c.grid:  # Gaussian output grid
    686             self.grid=c.grid
    687             self.area='G'
    688         else:
    689             self.grid='{}/{}'.format(int(c.grid)/1000.,int(c.grid)/1000.)
    690             self.area='{}/{}/{}/{}'.format(int(c.upper)/1000.,int(c.left)/1000.,int(c.lower)/1000.,int(c.right)/1000.)
    691 
    692         self.accuracy=c.accuracy
    693         self.level=c.level
    694         try:
    695             self.levelist=c.levelist
    696         except:
    697             self.levelist='1/to/'+c.level
    698         self.glevelist='1/to/'+c.level
    699         try:
    700             self.gaussian=c.gaussian
    701         except:
    702             self.gaussian=''
    703         try:
    704             self.dataset=c.dataset
    705         except:
    706             self.dataset=''
    707         try:
    708             self.expver=c.expver
    709         except:
    710             self.expver='1'
    711         try:
    712             self.number=c.number
    713         except:
    714             self.number='0'
    715        
    716         self.outputfilelist=[]
    717            
    718        
    719 # Now comes the nasty part that deals with the different scenarios we have:
    720 # 1) Calculation of etadot on
    721 #    a) Gaussian grid
    722 #    b) Output grid
    723 #    c) Output grid using parameter 77 retrieved from MARS
    724 # 3) Calculation/Retrieval of omega
    725 # 4) Download also data for WRF
    726 
    727         if fluxes==False:
    728             self.params['SH__SL']=['LNSP','ML','1','OFF']
    729 #           self.params['OG__SL']=["SD/MSL/TCC/10U/10V/2T/2D/129/172",'SFC','1',self.grid]
    730             self.params['OG__SL']=["141/151/164/165/166/167/168/129/172",'SFC','1',self.grid]
    731             self.params['OG_OROLSM__SL']=["160/27/28/173",'SFC','1',self.grid]
    732            
    733             if len(c.addpar)>0:
    734                 if c.addpar[0]=='/':
    735                     c.addpar=c.addpar[1:]
    736                 self.params['OG__SL'][0]+='/'+'/'.join(c.addpar)
    737             self.params['OG__ML']=['T/Q','ML',self.levelist,self.grid]
    738            
    739             if c.gauss=='0' and c.eta=='1': # the simplest case
    740                 self.params['OG__ML'][0]+='/U/V/77'
    741             elif c.gauss=='0' and c.eta=='0': # this is not recommended (inaccurate)
    742                 self.params['OG__ML'][0]+='/U/V'
    743             elif c.gauss=='1' and c.eta=='0': # this is needed for data before 2008, or for reanalysis data
    744                 self.params['GG__SL']=['Q','ML','1','{}'.format((int(self.resol)+1)/2)]
    745                 self.params['SH__ML']=['U/V/D','ML',self.glevelist,'OFF']
    746             else:
    747                 print 'Warning: This is a very costly parameter combination, use only for debugging!'
    748                 self.params['GG__SL']=['Q','ML','1','{}'.format((int(self.resol)+1)/2)]
    749                 self.params['GG__ML']=['U/V/D/77','ML',self.glevelist,'{}'.format((int(self.resol)+1)/2)]
    750                
    751             if c.omega=='1':
    752                 self.params['OG__ML'][0]+='/W'
    753                
    754             try:            # add cloud water content if necessary
    755                 if c.cwc=='1':
    756                     self.params['OG__ML'][0]+='/CLWC/CIWC'
    757             except:
    758                 pass
    759            
    760             try:            # add vorticity and geopotential height for WRF if necessary
    761                 if c.wrf=='1':
    762                     self.params['OG__ML'][0]+='/Z/VO'
    763                     if '/D' not in self.params['OG__ML'][0]:
    764                         self.params['OG__ML'][0]+='/D'
    765 #                   wrf_sfc='sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper()
    766 #                            134/151/235/167/165/166/168/129/172/034/031/141/139/170/183/236/039/040/041/042
    767                     wrf_sfc='134/235/167/165/166/168/129/172/34/31/141/139/170/183/236/39/40/41/42'.upper()
    768                     lwrt_sfc=wrf_sfc.split('/')
    769                     for par in lwrt_sfc:
    770                         if par not in self.params['OG__SL'][0]:
    771                             self.params['OG__SL'][0]+='/'+par
    772                    
    773             except:
    774                 pass
    775         else:
    776             self.params['OG_acc_SL']=["LSP/CP/SSHF/EWSS/NSSS/SSR",'SFC','1',self.grid]
    777                
    778                
    779         return
    780         # add additional WRF specific parameters here
    781                
    782     def create_namelist(self,c,filename):
    783        
    784 #       if inputdir=="":
    785 #           self.inputdir='.'
    786 #       else:
    787 #           self.inputdir=inputdir
    788        
    789         self.inputdir=c.inputdir
    790         area=asarray(self.area.split('/')).astype(float)
    791         grid=asarray(self.grid.split('/')).astype(float)
    792        
    793 #       if area[3]<0:
    794 #           area[3]+=360
    795         if area[1]>area[3]:
    796             area[1]-=360
    797         zyk=abs((area[3]-area[1]-360.)+grid[1])<1.e-6
    798         maxl=int((area[3]-area[1])/grid[1])+1
    799         maxb=int((area[0]-area[2])/grid[0])+1
    800 
    801         f=open(self.inputdir+'/'+filename,'w')
    802         f.write('&NAMGEN\n')
    803         f.write(',\n  '.join(['maxl='+str(maxl),'maxb='+str(maxb),
    804                     'mlevel='+self.level,'mlevelist='+'"'+self.levelist+'"',
    805                     'mnauf='+self.resol,'metapar='+'77',
    806                     'rlo0='+str(area[1]),'rlo1='+str(area[3]),'rla0='+str(area[2]),'rla1='+str(area[0]),
    807                     'momega='+c.omega,'momegadiff='+c.omegadiff,'mgauss='+c.gauss,
    808                     'msmooth='+c.smooth,'meta='+c.eta,'metadiff='+c.etadiff,'mdpdeta='+c.dpdeta]))
    809        
    810         f.write('\n/\n')
    811         f.close()
    812         return
    813        
    814              
    815     def retrieve(self, server, dates,times, inputdir=''):
    816         self.dates=dates
    817         self.server=server
    818 
    819         if inputdir=="":
    820             self.inputdir='.'
    821         else:
    822             self.inputdir=inputdir
    823 
    824           # Retrieve Q not for using Q but as a template for a reduced gaussian grid one date and time is enough
    825 # Take analysis at 00
    826         qdate=self.dates
    827         idx=qdate.find("/")
    828         if (idx >0):
    829           qdate=self.dates[:idx]
    830        
    831         #QG= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream, type="an", levtype="ML", levelist="1",
    832                          #gaussian="reduced",grid='{}'.format((int(self.resol)+1)/2), resol=self.resol,accuracy=self.accuracy,target=self.inputdir+"/"+"QG.grb",
    833                          #date=qdate, time="00",expver=self.expver, param="133.128")
     969    '''
     970    Class to retrieve Era Interim data for running FLEXPART
     971    '''
     972#AP change class name to ECFlexpart
     973    def __init__(self, c, fluxes=False):
     974        '''
     975        @Description:
     976            Creates an object/instance of EIFlexpart with the
     977            associated settings of its attributes for the retrieval.
     978
     979        @Input:
     980            self: instance of EIFlexpart
     981                The current object of the class.
     982
     983            c: instance of class Control
     984                Contains all the parameters of control files, which are:
     985                DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     986                NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL,
     987                LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA,
     988                ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX,
     989                ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL,
     990                GRIB2FLEXPART, FLEXPARTDIR
     991                For more information about format and content of the parameter
     992                see documentation.
     993
     994            fluxes: boolean, optional
     995                Decides if a the flux parameter settings are stored or
     996                the rest of the parameter list.
     997                Default value is False.
     998
     999        @Return:
     1000            <nothing>
     1001        '''
     1002        # different mars types for retrieving reanalysis data for flexpart
     1003
     1004        self.types = dict()
     1005        try:
     1006            if c.maxstep > len(c.type):    # Pure forecast mode
     1007                c.type = [c.type[1]]
     1008                c.step = ['{:0>3}'.format(int(c.step[0]))]
     1009                c.time = [c.time[0]]
     1010                for i in range(1, c.maxstep + 1):
     1011                    c.type.append(c.type[0])
     1012                    c.step.append('{:0>3}'.format(i))
     1013                    c.time.append(c.time[0])
     1014        except:
     1015            pass
     1016
     1017        self.inputdir = c.inputdir
     1018        self.basetime = c.basetime
     1019        self.dtime = c.dtime
     1020        self.mars = {}
     1021        i = 0
     1022        j = 0
     1023        if fluxes is True and c.maxstep < 24:
     1024            # only FC with start times at 00/12 (without 00/12)
     1025            self.types[c.type[1]] = {'times': '00/12',
     1026                                     'steps': '{}/to/12/by/{}'.format(c.dtime,
     1027                                                                      c.dtime)}
     1028            i = 1
     1029            for k in [0, 12]:
     1030                for j in range(int(c.dtime), 13, int(c.dtime)):
     1031                    self.mars['{:0>3}'.format(i * int(c.dtime))] = \
     1032                           [c.type[1], '{:0>2}'.format(k), '{:0>3}'.format(j)]
     1033                    i += 1
     1034        else:
     1035            for ty, st, ti in zip(c.type, c.step, c.time):
     1036                btlist = range(24)
     1037                if c.basetime == '12':
     1038                    btlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
     1039                if c.basetime == '00':
     1040                    btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0]
     1041
     1042                if mod(i, int(c.dtime)) == 0 and \
     1043                    (i in btlist or c.maxstep > 24):
     1044
     1045                    if ty not in self.types.keys():
     1046                        self.types[ty] = {'times': '', 'steps': ''}
     1047
     1048                    if ti not in self.types[ty]['times']:
     1049                        if len(self.types[ty]['times']) > 0:
     1050                            self.types[ty]['times'] += '/'
     1051                        self.types[ty]['times'] += ti
     1052
     1053                    if st not in self.types[ty]['steps']:
     1054                        if len(self.types[ty]['steps']) > 0:
     1055                            self.types[ty]['steps'] += '/'
     1056                        self.types[ty]['steps'] += st
     1057
     1058                    self.mars['{:0>3}'.format(j)] = [ty,
     1059                                                     '{:0>2}'.format(int(ti)),
     1060                                                     '{:0>3}'.format(int(st))]
     1061                    j += int(c.dtime)
     1062
     1063                i += 1
     1064
     1065        # Different grids need different retrievals
     1066        # SH = Spherical Harmonics, GG = Gaussian Grid,
     1067        # OG = Output Grid, ML = MultiLevel, SL = SingleLevel
     1068        self.params = {'SH__ML': '', 'SH__SL': '',
     1069                       'GG__ML': '', 'GG__SL': '',
     1070                       'OG__ML': '', 'OG__SL': '',
     1071                       'OG_OROLSM_SL': '', 'OG_acc_SL': ''}
     1072        self.marsclass = c.marsclass
     1073        self.stream = c.stream
     1074        self.number = c.number
     1075        self.resol = c.resol
     1076        if 'N' in c.grid:  # Gaussian output grid
     1077            self.grid = c.grid
     1078            self.area = 'G'
     1079        else:
     1080            self.grid = '{}/{}'.format(int(c.grid) / 1000., int(c.grid) / 1000.)
     1081            self.area = '{}/{}/{}/{}'.format(int(c.upper) / 1000.,
     1082                                             int(c.left) / 1000.,
     1083                                             int(c.lower) / 1000.,
     1084                                             int(c.right) / 1000.)
     1085
     1086        self.accuracy = c.accuracy
     1087        self.level = c.level
     1088        try:
     1089            self.levelist = c.levelist
     1090        except:
     1091            self.levelist = '1/to/' + c.level
     1092
     1093        self.glevelist = '1/to/' + c.level
     1094
     1095        try:
     1096            self.gaussian = c.gaussian
     1097        except:
     1098            self.gaussian = ''
     1099
     1100        try:
     1101            self.dataset = c.dataset
     1102        except:
     1103            self.dataset = ''
     1104
     1105        try:
     1106            self.expver = c.expver
     1107        except:
     1108            self.expver = '1'
     1109
     1110        try:
     1111            self.number = c.number
     1112        except:
     1113            self.number = '0'
     1114
     1115        self.outputfilelist = []
     1116
     1117
     1118        # Now comes the nasty part that deals with the different scenarios we have:
     1119        # 1) Calculation of etadot on
     1120        #    a) Gaussian grid
     1121        #    b) Output grid
     1122        #    c) Output grid using parameter 77 retrieved from MARS
     1123        # 3) Calculation/Retrieval of omega
     1124        # 4) Download also data for WRF
     1125
     1126        if fluxes is False:
     1127            self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF']
     1128    #        self.params['OG__SL'] = ["SD/MSL/TCC/10U/10V/2T/2D/129/172",
     1129    #                                 'SFC','1',self.grid]
     1130            self.params['OG__SL'] = ["141/151/164/165/166/167/168/129/172", \
     1131                                     'SFC', '1', self.grid]
     1132            self.params['OG_OROLSM__SL'] = ["160/27/28/173", \
     1133                                            'SFC', '1', self.grid]
     1134
     1135            if len(c.addpar) > 0:
     1136                if c.addpar[0] == '/':
     1137                    c.addpar = c.addpar[1:]
     1138                self.params['OG__SL'][0] += '/' + '/'.join(c.addpar)
     1139            self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]
     1140
     1141            if c.gauss == '0' and c.eta == '1':
     1142                # the simplest case
     1143                self.params['OG__ML'][0] += '/U/V/77'
     1144            elif c.gauss == '0' and c.eta == '0':
     1145                # this is not recommended (inaccurate)
     1146                self.params['OG__ML'][0] += '/U/V'
     1147            elif c.gauss == '1' and c.eta == '0':
     1148                # this is needed for data before 2008, or for reanalysis data
     1149                self.params['GG__SL'] = ['Q', 'ML', '1', \
     1150                                         '{}'.format((int(self.resol) + 1) / 2)]
     1151                self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF']
     1152            else:
     1153                print('Warning: This is a very costly parameter combination, \
     1154                       use only for debugging!')
     1155                self.params['GG__SL'] = ['Q', 'ML', '1', \
     1156                                         '{}'.format((int(self.resol) + 1) / 2)]
     1157                self.params['GG__ML'] = ['U/V/D/77', 'ML', self.glevelist, \
     1158                                         '{}'.format((int(self.resol) + 1) / 2)]
     1159
     1160            if c.omega == '1':
     1161                self.params['OG__ML'][0] += '/W'
     1162
     1163            try:
     1164                # add cloud water content if necessary
     1165                if c.cwc == '1':
     1166                    self.params['OG__ML'][0] += '/CLWC/CIWC'
     1167            except:
     1168                pass
     1169
     1170            try:
     1171                # add vorticity and geopotential height for WRF if necessary
     1172                if c.wrf == '1':
     1173                    self.params['OG__ML'][0] += '/Z/VO'
     1174                    if '/D' not in self.params['OG__ML'][0]:
     1175                        self.params['OG__ML'][0] += '/D'
     1176        #            wrf_sfc = 'sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/ /
     1177        #                       stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper()
     1178                    wrf_sfc = '134/235/167/165/166/168/129/172/34/31/141/ \
     1179                               139/170/183/236/39/40/41/42'.upper()
     1180                    lwrt_sfc = wrf_sfc.split('/')
     1181                    for par in lwrt_sfc:
     1182                        if par not in self.params['OG__SL'][0]:
     1183                            self.params['OG__SL'][0] += '/' + par
     1184            except:
     1185                pass
     1186        else:
     1187            self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR", \
     1188                                        'SFC', '1', self.grid]
     1189
     1190        # if needed, add additional WRF specific parameters here
     1191
     1192        return
     1193
     1194
     1195    def create_namelist(self, c, filename):
     1196        '''
     1197        @Description:
     1198            Creates a namelist file in the temporary directory.
     1199            It will contain values for maxl, maxb, mlevel,
     1200            mlevelist, mnauf, metapar, rlo0, rlo1, rla0, rla1,
     1201            momega, momegadiff, mgauss, msmooth, meta, metadiff, mdpdeta
     1202
     1203        @Input:
     1204            self: instance of EIFlexpart
     1205                The current object of the class.
     1206
     1207            c: instance of class Control
     1208                Contains all the parameters of control files, which are:
     1209                DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     1210                NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL,
     1211                LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA,
     1212                ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX,
     1213                ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL,
     1214                GRIB2FLEXPART, FLEXPARTDIR
     1215                For more information about format and content of the parameter
     1216                see documentation.
     1217
     1218            filename: string
     1219                Name of the namelist file.
     1220
     1221        @Return:
     1222            <nothing>
     1223        '''
     1224        self.inputdir = c.inputdir
     1225        area = asarray(self.area.split('/')).astype(float)
     1226        grid = asarray(self.grid.split('/')).astype(float)
     1227
     1228        if area[1] > area[3]:
     1229            area[1] -= 360
     1230        zyk = abs((area[3] - area[1] - 360.) + grid[1]) < 1.e-6
     1231        maxl = int((area[3] - area[1]) / grid[1]) + 1
     1232        maxb = int((area[0] - area[2]) / grid[0]) + 1
     1233
     1234        f = open(self.inputdir + '/' + filename, 'w')
     1235        f.write('&NAMGEN\n')
     1236        f.write(',\n  '.join(['maxl = ' + str(maxl), 'maxb = ' + str(maxb),
     1237                'mlevel = ' + self.level,
     1238                'mlevelist = ' + '"' + self.levelist + '"',
     1239                'mnauf = ' + self.resol, 'metapar = ' + '77',
     1240                'rlo0 = ' + str(area[1]), 'rlo1 = ' + str(area[3]),
     1241                'rla0 = ' + str(area[2]), 'rla1 = ' + str(area[0]),
     1242                'momega = ' + c.omega, 'momegadiff = ' + c.omegadiff,
     1243                'mgauss = ' + c.gauss, 'msmooth = ' + c.smooth,
     1244                'meta = ' + c.eta, 'metadiff = ' + c.etadiff,
     1245                'mdpdeta = ' + c.dpdeta]))
     1246
     1247        f.write('\n/\n')
     1248        f.close()
     1249        return
     1250
     1251    def retrieve(self, server, dates, times, inputdir=''):
     1252        '''
     1253        @Description:
     1254
     1255
     1256        @Input:
     1257            self: instance of EIFlexpart
     1258
     1259            server: instance of ECMWFService
     1260
     1261            dates:
     1262
     1263            times:
     1264
     1265            inputdir: string
     1266                Default string is empty ('').
     1267
     1268        @Return:
     1269            <nothing>
     1270        '''
     1271        self.dates = dates
     1272        self.server = server
     1273
     1274        if inputdir == "":
     1275            self.inputdir = '.'
     1276        else:
     1277            self.inputdir = inputdir
     1278
     1279        # Retrieve Q not for using Q but as a template for a reduced gaussian
     1280        # grid one date and time is enough
     1281        # Take analysis at 00
     1282        qdate = self.dates
     1283        idx = qdate.find("/")
     1284        if (idx > 0):
     1285            qdate = self.dates[:idx]
     1286
     1287        #QG =  MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream, type = "an", levtype = "ML", levelist = "1",
     1288                             #gaussian = "reduced",grid = '{}'.format((int(self.resol)+1)/2), resol = self.resol,accuracy = self.accuracy,target = self.inputdir+"/"+"QG.grb",
     1289                             #date = qdate, time = "00",expver = self.expver, param = "133.128")
    8341290        #QG.displayInfo()
    8351291        #QG.dataRetrieve()
    8361292
    837         oro=False
     1293        oro = False
    8381294        for ftype in self.types:
    839             for pk,pv in self.params.iteritems():
    840                 if isinstance(pv,str):     # or pk!='GG__SL' :
    841                     continue
    842                 mftype=''+ftype
    843                 mftime=self.types[ftype]['times']
    844                 mfstep=self.types[ftype]['steps']
    845                 mfdate=self.dates
    846                 mfstream=self.stream
    847                 mftarget=self.inputdir+"/"+ftype+pk+'.'+self.dates.split('/')[0]+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
    848                 if pk=='OG__SL':
    849                     pass
    850                 if pk=='OG_OROLSM__SL':
    851                     if oro==False:
    852                         mfstream='OPER'
    853                         mftype='AN'
    854                         mftime='00'
    855                         mfstep='000'
    856                         mfdate=self.dates.split('/')[0]
    857                         mftarget=self.inputdir+"/"+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
    858                         oro=True
    859                     else:
    860                         continue
    861                
    862                 if pk=='GG__SL' and pv[0]=='Q':
    863                     area=""
    864                     gaussian='reduced'
    865                 else:
    866                     area=self.area
    867                     gaussian=self.gaussian
    868 
    869                 if self.basetime==None:
    870                     MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=mfstream,
    871                                   type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
    872                                   accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
    873                                   date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
    874                
    875                     MR.displayInfo()
    876                     MR.dataRetrieve()
    877 # The whole else section is only necessary for operational scripts. It could be removed
    878                 else: # check if mars job requests fields beyond basetime. If yes eliminate those fields since they may not
    879                       # be accessible with user's credentials
    880                     sm1=-1
    881                     if 'by' in mfstep:
    882                         sm1=2
    883                     tm1=-1
    884                     if 'by' in mftime:
    885                         tm1=2
    886                     maxtime=datetime.datetime.strptime(mfdate.split('/')[-1]+mftime.split('/')[tm1],'%Y%m%d%H')+ \
    887                     datetime.timedelta(hours=int(mfstep.split('/')[sm1]))
    888                    
    889                     elimit=datetime.datetime.strptime(mfdate.split('/')[-1]+self.basetime,'%Y%m%d%H')
    890 
    891                     if self.basetime=='12':
    892                         if 'acc' in pk:
    893 
    894                     # Strategy: if maxtime-elimit>=24h reduce date by 1,
    895                     # if 12h<=maxtime-elimit<12h reduce time for last date
    896                     # if maxtime-elimit<12h reduce step for last time
    897                     # A split of the MARS job into 2 is likely necessary.
    898                             maxtime=elimit-datetime.timedelta(hours=24)
    899                             mfdate='/'.join(('/'.join(mfdate.split('/')[:-1]),datetime.datetime.strftime(maxtime,'%Y%m%d')))
    900                    
    901                             MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
    902                                                 type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
    903                                                 accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
    904                                                 date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
    905                
    906                             MR.displayInfo()
    907                             MR.dataRetrieve()
    908                    
    909                             maxtime=elimit-datetime.timedelta(hours=12)
    910                             mfdate=datetime.datetime.strftime(maxtime,'%Y%m%d')
    911                             mftime='00'
    912                             mftarget=self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
    913 
    914                             MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
    915                                                 type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
    916                                                 accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
    917                                                 date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
    918                
    919                             MR.displayInfo()
    920                             MR.dataRetrieve()
    921                         else:
    922                             MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
    923                                                 type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
    924                                                 accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
    925                                                 date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
    926                
    927                             MR.displayInfo()
    928                             MR.dataRetrieve()
    929                     else:
    930                         maxtime=elimit-datetime.timedelta(hours=24)
    931                         mfdate=datetime.datetime.strftime(maxtime,'%Y%m%d')
    932                        
    933                         mftimesave=''.join(mftime)
    934                        
    935                         if '/' in mftime:
    936                             times=mftime.split('/')
    937                             while int(times[0])+int(mfstep.split('/')[0])<=12 and pk!='OG_OROLSM__SL' and 'acc' not in pk:
    938                                 times=times[1:]
    939                                 if len(times)>1:
    940                                     mftime='/'.join(times)
    941                                 else:
    942                                     mftime=times[0]
    943                            
    944                         MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
    945                                       type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
    946                                       accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
    947                                       date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
    948                    
    949                         MR.displayInfo()
    950                         MR.dataRetrieve()
    951                        
    952                         if int(mftimesave.split('/')[0])==0 and int(mfstep.split('/')[0])==0 and pk!='OG_OROLSM__SL':
    953                             mfdate=datetime.datetime.strftime(elimit,'%Y%m%d')
    954                             mftime='00'
    955                             mfstep='000'
    956                             mftarget=self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
    957                    
    958                             MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
    959                                           type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
    960                                           accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
    961                                           date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
    962                        
    963                             MR.displayInfo()
    964                             MR.dataRetrieve()
    965                                                        
    966                    
    967                    
    968         print "MARS retrieve done... "
    969        
    970     def getFlexpartTime(self, type,step, time):
    971         cstep='{:0>3}'.format(step)
    972         ctime='{:0>2}'.format(int(time/100))
    973 
     1295            for pk, pv in self.params.iteritems():
     1296                if isinstance(pv, str):     # or pk != 'GG__SL' :
     1297                    continue
     1298                mftype = ''+ftype
     1299                mftime = self.types[ftype]['times']
     1300                mfstep = self.types[ftype]['steps']
     1301                mfdate = self.dates
     1302                mfstream = self.stream
     1303                mftarget = self.inputdir+"/"+ftype+pk+'.'+self.dates.split('/')[0]+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
     1304                if pk == 'OG__SL':
     1305                    pass
     1306                if pk == 'OG_OROLSM__SL':
     1307                    if oro is False:
     1308                        mfstream = 'OPER'
     1309                        mftype = 'AN'
     1310                        mftime = '00'
     1311                        mfstep = '000'
     1312                        mfdate = self.dates.split('/')[0]
     1313                        mftarget = self.inputdir + "/" + pk + '.' + mfdate + \
     1314                                   '.' + str(os.getppid()) + '.' + \
     1315                                   str(os.getpid()) + ".grb"
     1316                        oro = True
     1317                    else:
     1318                        continue
     1319
     1320                if pk == 'GG__SL' and pv[0] == 'Q':
     1321                    area = ""
     1322                    gaussian = 'reduced'
     1323                else:
     1324                    area = self.area
     1325                    gaussian = self.gaussian
     1326
     1327                if self.basetime is None:
     1328                    MR =  MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = mfstream,
     1329                              type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian,
     1330                              accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area,
     1331                              date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0])
     1332
     1333                    MR.displayInfo()
     1334                    MR.dataRetrieve()
     1335    # The whole else section is only necessary for operational scripts.
     1336    # It could be removed
     1337                else:
     1338                    # check if mars job requests fields beyond basetime.
     1339                    # If yes eliminate those fields since they may not
     1340                    # be accessible with user's credentials
     1341                    sm1 = -1
     1342                    if 'by' in mfstep:
     1343                        sm1 = 2
     1344                    tm1 = -1
     1345                    if 'by' in mftime:
     1346                        tm1 = 2
     1347                    maxtime = datetime.datetime.strptime(mfdate.split('/')[-1]+mftime.split('/')[tm1],'%Y%m%d%H')+ \
     1348                    datetime.timedelta(hours = int(mfstep.split('/')[sm1]))
     1349
     1350                    elimit = datetime.datetime.strptime(mfdate.split('/')[-1]+self.basetime,'%Y%m%d%H')
     1351
     1352                    if self.basetime == '12':
     1353                        if 'acc' in pk:
     1354
     1355                # Strategy: if maxtime-elimit> = 24h reduce date by 1,
     1356                # if 12h< = maxtime-elimit<12h reduce time for last date
     1357                # if maxtime-elimit<12h reduce step for last time
     1358                # A split of the MARS job into 2 is likely necessary.
     1359                            maxtime = elimit-datetime.timedelta(hours = 24)
     1360                            mfdate = '/'.join(('/'.join(mfdate.split('/')[:-1]),datetime.datetime.strftime(maxtime,'%Y%m%d')))
     1361
     1362                            MR =  MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream,
     1363                                                type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian,
     1364                                                accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area,
     1365                                                date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0])
     1366
     1367                            MR.displayInfo()
     1368                            MR.dataRetrieve()
     1369
     1370                            maxtime = elimit-datetime.timedelta(hours = 12)
     1371                            mfdate = datetime.datetime.strftime(maxtime,'%Y%m%d')
     1372                            mftime = '00'
     1373                            mftarget = self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
     1374
     1375                            MR =  MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream,
     1376                                                type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian,
     1377                                                accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area,
     1378                                                date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0])
     1379
     1380                            MR.displayInfo()
     1381                            MR.dataRetrieve()
     1382                        else:
     1383                            MR =  MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream,
     1384                                        type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian,
     1385                                        accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area,
     1386                                        date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0])
     1387
     1388                            MR.displayInfo()
     1389                            MR.dataRetrieve()
     1390                    else:
     1391                        maxtime = elimit-datetime.timedelta(hours = 24)
     1392                        mfdate = datetime.datetime.strftime(maxtime,'%Y%m%d')
     1393
     1394                        mftimesave = ''.join(mftime)
     1395
     1396                        if '/' in mftime:
     1397                            times = mftime.split('/')
     1398                            while int(times[0])+int(mfstep.split('/')[0])< = 12 and pk! = 'OG_OROLSM__SL' and 'acc' not in pk:
     1399                                times = times[1:]
     1400                            if len(times)>1:
     1401                                mftime = '/'.join(times)
     1402                            else:
     1403                                mftime = times[0]
     1404
     1405                        MR =  MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream,
     1406                                  type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian,
     1407                                  accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area,
     1408                                  date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0])
     1409
     1410                        MR.displayInfo()
     1411                        MR.dataRetrieve()
     1412
     1413                        if int(mftimesave.split('/')[0]) == 0 and int(mfstep.split('/')[0]) == 0 and pk! = 'OG_OROLSM__SL':
     1414                            mfdate = datetime.datetime.strftime(elimit,'%Y%m%d')
     1415                            mftime = '00'
     1416                            mfstep = '000'
     1417                            mftarget = self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
     1418
     1419                            MR =  MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream,
     1420                                          type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian,
     1421                                          accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area,
     1422                                          date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0])
     1423
     1424                            MR.displayInfo()
     1425                            MR.dataRetrieve()
     1426
     1427            print("MARS retrieve done... ")
     1428
     1429            return
     1430
     1431
     1432    def getFlexpartTime(self, type, step, time):
     1433#AP remove this function, no longer needed
     1434        '''
     1435        @Description:
     1436            ????
     1437#AP wozu? wird das überhaupt noch gebraucht? in create auskommentiert
     1438
     1439        @Input:
     1440            self: instance of EIFlexpart
     1441                The current object of the class.
     1442
     1443            type: string
     1444                Type of the data field. E.g. FC - forecast or
     1445                AN - analysis
     1446
     1447            step: integer
     1448                Forecast step in hours.
     1449
     1450            time: integer
     1451                Time in hours.
     1452
     1453        @Return:
     1454            cflextime:
     1455
     1456        '''
     1457        cstep = '{:0>3}'.format(step)
     1458        ctime = '{:0>2}'.format(int(time / 100))
    9741459        ctype = str(type).upper()
    975         myinfo = [ctype,ctime, cstep]
     1460
     1461        myinfo = [ctype, ctime, cstep]
    9761462        cflextime = None
    9771463        for t, marsinfo in self.mars.items():
    9781464            if myinfo == marsinfo:
    979                 cflextime=t
     1465                cflextime = t
     1466
    9801467        return cflextime
    9811468
    9821469    def process_output(self, c):
    983        
    984         print 'Postprocessing:\n Format: {}\n'.format(c.format)
    985         if c.ecapi==False:
    986             print 'ecstorage: {}\n ecfsdir: {}\n'.format(c.ecstorage,c.ecfsdir)
    987             if not hasattr(c,'gateway'):
    988                 c.gateway=os.getenv('GATEWAY')
    989             if not hasattr(c,'destination'):
    990                 c.destination=os.getenv('DESTINATION')
    991             print 'ectrans: {}\n gateway: {}\n destination: {}\n '.format(c.ectrans, c.gateway,c.destination)
    992         print 'Output filelist: \n',self.outputfilelist
    993        
    994         if c.format.lower()=='grib2':
    995             for ofile in self.outputfilelist:
    996                 p=subprocess.check_call(['grib_set','-s','edition=2,productDefinitionTemplateNumber=8',ofile,ofile+'_2'])
    997                 p=subprocess.check_call(['mv',ofile+'_2',ofile])
    998        
    999         if int(c.ectrans)==1 and c.ecapi==False:
    1000             for ofile in self.outputfilelist:
    1001                 p=subprocess.check_call(['ectrans','-overwrite','-gateway',c.gateway,
    1002                                          '-remote',c.destination,'-source',ofile])
    1003                 print 'ectrans:',p
    1004         if int(c.ecstorage)==1 and c.ecapi==False:
    1005             for ofile in self.outputfilelist:
    1006                 p=subprocess.check_call(['ecp','-o',ofile,os.path.expandvars(c.ecfsdir)])
    1007 
    1008 #20131107 000000      EN13110700              ON DISC
    1009         if c.outputdir!=c.inputdir:
    1010             for ofile in self.outputfilelist:
    1011                 p=subprocess.check_call(['mv',ofile,c.outputdir])
    1012            
    1013         if c.grib2flexpart=='1':
    1014             f=open(c.outputdir+'/'+'AVAILABLE','w')
    1015             clist=[]
    1016             for ofile in self.outputfilelist: # generate AVAILABLE file
    1017                 fname=ofile.split('/')
    1018                 if '.' in fname[-1]:
    1019                     l=fname[-1].split('.')
    1020                     timestamp=datetime.datetime.strptime(l[0][-6:]+l[1],'%y%m%d%H')
    1021                     timestamp+=datetime.timedelta(hours=int(l[2]))
    1022                     cdate=datetime.datetime.strftime(timestamp,'%Y%m%d')
    1023                     chms=datetime.datetime.strftime(timestamp,'%H%M%S')
    1024                    
    1025                 else:
    1026                     cdate='20'+fname[-1][-8:-2]
    1027                     chms=fname[-1][-2:]+'0000'
    1028                 clist.append(cdate+' '+chms+' '*6+fname[-1]+' '*14+'ON DISC')
    1029             clist.sort()
    1030             f.write('\n'.join(clist)+'\n')
    1031             f.close()
    1032            
    1033             pwd=os.path.abspath(c.outputdir)
    1034             f=open(pwd+'/pathnames','w')
    1035             f.write(pwd+'/Options/\n')
    1036             f.write(pwd+'/\n')
    1037             f.write(pwd+'/\n')
    1038             f.write(pwd+'/AVAILABLE\n')
    1039             f.write('==================\n')
    1040             f.close()
    1041            
    1042             if not os.path.exists(pwd+'/Options'):
    1043                 os.makedirs(pwd+'/Options')
    1044             f=open(os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))+'/../Options/COMMAND','r')
    1045             lflist=f.read().split('\n')
    1046             i=0
    1047             for l in lflist:
    1048                 if 'LDIRECT' in l.upper():
    1049                     break
    1050                 i+=1
    1051                
    1052             clist.sort()
    1053             lflist=lflist[:i+1]+[clist[0][:16],clist[-1][:16]]+lflist[i+3:]
    1054             g=open(pwd+'/Options/COMMAND','w')
    1055             g.write('\n'.join(lflist)+'\n')
    1056             g.close()
    1057            
    1058             os.chdir(c.outputdir)
    1059             p=subprocess.check_call([os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))+'/../FLEXPART_PROGRAM/grib2flexpart',
    1060                                      'useAvailable','.'])
    1061             os.chdir(pwd)
    1062            
    1063     def create(self, inputfiles, c,):
    1064        
    1065         table128=init128(c.ecmwfdatadir+'/grib_templates/ecmwf_grib1_table_128')
    1066         wrfpars=toparamId('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',table128)
    1067 #        index_keys=["date","time","stepRange"]
    1068         index_keys=["date","time","step"]
    1069         indexfile=c.inputdir+"/date_time_stepRange.idx"
     1470        '''
     1471        @Description:
     1472
     1473
     1474        @Input:
     1475            self: instance of EIFlexpart
     1476                The current object of the class.
     1477
     1478            c: instance of class Control
     1479                Contains all the parameters of control files, which are:
     1480                DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     1481                NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL,
     1482                LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA,
     1483                ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX,
     1484                ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL,
     1485                GRIB2FLEXPART, FLEXPARTDIR
     1486                For more information about format and content of the parameter
     1487                see documentation.
     1488
     1489        @Return:
     1490            <nothing>
     1491
     1492        '''
     1493
     1494        print 'Postprocessing:\n Format: {}\n'.format(c.format)
     1495        if c.ecapi is False:
     1496            print 'ecstorage: {}\n ecfsdir: {}\n'.format(c.ecstorage, c.ecfsdir)
     1497            if not hasattr(c, 'gateway'):
     1498                c.gateway = os.getenv('GATEWAY')
     1499            if not hasattr(c, 'destination'):
     1500                c.destination = os.getenv('DESTINATION')
     1501            print 'ectrans: {}\n gateway: {}\n destination: {}\n '}
     1502                    .format(c.ectrans, c.gateway, c.destination)
     1503        print 'Output filelist: \n', self.outputfilelist
     1504
     1505        if c.format.lower() == 'grib2':
     1506            for ofile in self.outputfilelist:
     1507                p = subprocess.check_call(['grib_set', '-s', 'edition = 2, \
     1508                                        productDefinitionTemplateNumber = 8',
     1509                                        ofile, ofile + '_2'])
     1510                p = subprocess.check_call(['mv', ofile + '_2', ofile])
     1511
     1512        if int(c.ectrans) == 1 and c.ecapi is False:
     1513            for ofile in self.outputfilelist:
     1514                p = subprocess.check_call(['ectrans', '-overwrite', '-gateway',
     1515                                           c.gateway, '-remote', c.destination,
     1516                                           '-source', ofile])
     1517                print 'ectrans:',p
     1518        if int(c.ecstorage) == 1 and c.ecapi is False:
     1519            for ofile in self.outputfilelist:
     1520                p = subprocess.check_call(['ecp', '-o', ofile,
     1521                                           os.path.expandvars(c.ecfsdir)])
     1522
     1523#    20131107 000000      EN13110700              ON DISC
     1524        if c.outputdir != c.inputdir:
     1525            for ofile in self.outputfilelist:
     1526                p = subprocess.check_call(['mv', ofile, c.outputdir])
     1527
     1528        if c.grib2flexpart == '1':
     1529            f = open(c.outputdir + '/' + 'AVAILABLE', 'w')
     1530            clist = []
     1531            for ofile in self.outputfilelist:  # generate AVAILABLE file
     1532                fname = ofile.split('/')
     1533                if '.' in fname[-1]:
     1534                    l = fname[-1].split('.')
     1535                    timestamp = datetime.datetime.strptime(l[0][-6:] + l[1],
     1536                                                           '%y%m%d%H')
     1537                    timestamp += datetime.timedelta(hours=int(l[2]))
     1538                    cdate = datetime.datetime.strftime(timestamp, '%Y%m%d')
     1539                    chms = datetime.datetime.strftime(timestamp, '%H%M%S')
     1540                else:
     1541                    cdate = '20' + fname[-1][-8:-2]
     1542                    chms = fname[-1][-2:] + '0000'
     1543                clist.append(cdate + ' ' + chms + ' '*6 +
     1544                             fname[-1] + ' '*14 + 'ON DISC')
     1545            clist.sort()
     1546            f.write('\n'.join(clist) + '\n')
     1547            f.close()
     1548
     1549            pwd = os.path.abspath(c.outputdir)
     1550            f = open(pwd + '/pathnames','w')
     1551            f.write(pwd + '/Options/\n')
     1552            f.write(pwd + '/\n')
     1553            f.write(pwd + '/\n')
     1554            f.write(pwd + '/AVAILABLE\n')
     1555            f.write(' = == = == = == = == = == ==  = \n')
     1556            f.close()
     1557
     1558            if not os.path.exists(pwd + '/Options'):
     1559                os.makedirs(pwd+'/Options')
     1560            f = open(os.path.expandvars(
     1561                     os.path.expanduser(c.flexpart_root_scripts)) +
     1562                     '/../Options/COMMAND', 'r')
     1563            lflist = f.read().split('\n')
     1564            i = 0
     1565            for l in lflist:
     1566                if 'LDIRECT' in l.upper():
     1567                    break
     1568                i += 1
     1569
     1570            clist.sort()
     1571            lflist = lflist[:i+1] + \
     1572                     [clist[0][:16], clist[-1][:16]] + \
     1573                     lflist[i+3:]
     1574
     1575            with open(pwd + '/Options/COMMAND', 'w') as g:
     1576                g.write('\n'.join(lflist) + '\n')
     1577
     1578            os.chdir(c.outputdir)
     1579            p = subprocess.check_call([os.path.expandvars(
     1580                        os.path.expanduser(c.flexpart_root_scripts)) +
     1581                        '/../FLEXPART_PROGRAM/grib2flexpart',
     1582                        'useAvailable', '.'])
     1583            os.chdir(pwd)
     1584
     1585        return
     1586
     1587    def create(self, inputfiles, c):
     1588        '''
     1589        @Description:
     1590#AP WOZU und ist einrueckung richtig
     1591
     1592        @Input:
     1593            self: instance of EIFlexpart
     1594                The current object of the class.
     1595
     1596            inputfiles: list of strings
     1597                A list of filenames.
     1598
     1599            c: instance of class Control
     1600                Contains all the parameters of control files, which are:
     1601                DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     1602                NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL,
     1603                LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA,
     1604                ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX,
     1605                ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL,
     1606                GRIB2FLEXPART, FLEXPARTDIR
     1607                For more information about format and content of the parameter
     1608                see documentation.
     1609
     1610        @Return:
     1611            <nothing>
     1612        '''
     1613
     1614        table128 = init128(c.ecmwfdatadir +
     1615                           '/grib_templates/ecmwf_grib1_table_128')
     1616#AP wieso wrf
     1617        wrfpars = toparamId('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/ \
     1618                            stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
     1619                            table128)
     1620        index_keys = ["date", "time", "step"]
     1621        indexfile = c.inputdir + "/date_time_stepRange.idx"
    10701622        silentremove(indexfile)
    1071         grib=GribTools(inputfiles.files)
    1072         iid=grib.index(index_keys=index_keys, index_file = indexfile)
    1073 
    1074         print 'index done...'
    1075         fdict={'10':None,'11':None,'12':None,'13':None,'16':None,'17':None,'19':None,'21':None,'22':None,'20':None}
    1076         for f in fdict.keys():
    1077             silentremove(c.inputdir+"/fort."+f)
    1078 
     1623        grib = GribTools(inputfiles.files)
     1624        iid = grib.index(index_keys=index_keys, index_file=indexfile)
     1625        print 'index done...'
     1626
     1627        fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
     1628                 '17':None, '19':None, '21':None, '22':None, '20':None}
     1629        for f in fdict.keys():
     1630            silentremove(c.inputdir + "/fort." + f)
     1631
     1632        index_vals = []
     1633        for key in index_keys:
     1634            key_vals = grib_index_get(iid, key)
     1635            print key_vals
     1636            index_vals.append(key_vals)
     1637
     1638        # creates index file
     1639        for prod in product(*index_vals):
     1640            for i in range(len(index_keys)):
     1641                grib_index_select(iid, index_keys[i],  prod[i])
     1642
     1643            gid = grib_new_from_index(iid)
     1644            # do convert2 program if gid at this time is not None,
     1645            # therefore save in hid
     1646            hid = gid
     1647            cflextime = None
     1648            for k, f in fdict.iteritems():
     1649                fdict[k] = open(c.inputdir + '/fort.' + k, 'w')
     1650            if gid is not None:
     1651                cdate = str(grib_get(gid, 'date'))
     1652                time = grib_get(gid, 'time')
     1653                type = grib_get(gid, 'type')
     1654                step = grib_get(gid, 'step')
     1655        #       cflextime  =  self.getFlexpartTime(type,step, time)
     1656                timestamp = datetime.datetime.strptime(
     1657                                cdate + '{:0>2}'.format(time/100), '%Y%m%d%H')
     1658                timestamp += datetime.timedelta(hours=int(step))
     1659        #       print gid,index_keys[i],prod[i],cdate,time,step,timestamp
     1660
     1661                cdateH = datetime.datetime.strftime(timestamp, '%Y%m%d%H')
     1662                chms = datetime.datetime.strftime(timestamp, '%H%M%S')
     1663
     1664                if c.basetime is not None:
     1665                    slimit = datetime.datetime.strptime(
     1666                                c.start_date + '00', '%Y%m%d%H')
     1667                    bt = '23'
     1668                    if c.basetime == '00':
     1669                        bt = '00'
     1670                        slimit = datetime.datetime.strptime(
     1671                                    c.end_date + bt, '%Y%m%d%H') - \
     1672                                 datetime.timedelta(hours=12-int(c.dtime))
     1673                    if c.basetime == '12':
     1674                        bt = '12'
     1675                        slimit = datetime.datetime.strptime(
     1676                                    c.end_date + bt, '%Y%m%d%H') - \
     1677                                 datetime.timedelta(hours=12-int(c.dtime))
     1678
     1679                    elimit = datetime.datetime.strptime(c.end_date + bt, '%Y%m%d%H')
     1680
     1681                    if timestamp < slimit or timestamp > elimit:
     1682                        continue
     1683
     1684            try:
     1685                if c.wrf == '1':
     1686                    if 'olddate' not in locals():
     1687                        fwrf = open(c.outputdir + '/WRF' + cdate +
     1688                                    '.{:0>2}'.format(time) + '.000.grb2', 'w')
     1689                        olddate = cdate[:]
     1690                    else:
     1691                        if cdate != olddate:
     1692                            fwrf = open(c.outputdir + '/WRF' + cdate +
     1693                                        '.{:0>2}'.format(time) + '.000.grb2', 'w')
     1694                            olddate = cdate[:]
     1695            except AttributeError:
     1696                pass
     1697
     1698#           print 'cyear '+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime
     1699
     1700            savedfields = []
     1701            while 1:
     1702                if gid is None:
     1703                    break
     1704                paramId = grib_get(gid, 'paramId')
     1705                gridtype = grib_get(gid, 'gridType')
     1706                datatype = grib_get(gid, 'dataType')
     1707                levtype = grib_get(gid, 'typeOfLevel')
     1708                if paramId == 133 and gridtype == 'reduced_gg':
     1709                # Relative humidity (Q.grb) is used as a template only
     1710                # so we need the first we "meet"
     1711                    fout = open(c.inputdir + '/fort.18', 'w')
     1712                    grib_write(gid, fout)
     1713                    fout.close()
     1714                elif paramId == 131 or paramId == 132:
     1715                    grib_write(gid, fdict['10'])
     1716                elif paramId == 130:
     1717                    grib_write(gid, fdict['11'])
     1718                elif paramId == 133 and gridtype != 'reduced_gg':
     1719                    grib_write(gid, fdict['17'])
     1720                elif paramId == 152:
     1721                    grib_write(gid, fdict['12'])
     1722                elif paramId == 155 and gridtype == 'sh':
     1723                    grib_write(gid, fdict['13'])
     1724                elif paramId in [129, 138, 155] and levtype == 'hybrid' \
     1725                                                and c.wrf == '1':
     1726        #            print paramId,'not written'
     1727                    pass
     1728                elif paramId == 246 or paramId == 247:
     1729                    # cloud liquid water and ice
     1730                    if paramId == 246:
     1731                        clwc = grib_get_values(gid)
     1732                    else:
     1733                        clwc += grib_get_values(gid)
     1734                        grib_set_values(gid, clwc)
     1735            #            grib_set(gid,'shortName','qc')
     1736                        grib_set(gid, 'paramId', 201031)
     1737                        grib_write(gid, fdict['22'])
     1738                elif paramId == 135:
     1739                    grib_write(gid, fdict['19'])
     1740                elif paramId == 77:
     1741                    grib_write(gid, fdict['21'])
     1742                else:
     1743                    if paramId not in savedfields:
     1744                        grib_write(gid, fdict['16'])
     1745                        savedfields.append(paramId)
     1746                    else:
     1747                        print 'duplicate ' + str(paramId) + ' not written'
     1748
     1749                try:
     1750                    if c.wrf == '1':
     1751                        if levtype == 'hybrid':
     1752                            if paramId in [129, 130, 131, 132, 133, 138, 155]:
     1753                                grib_write(gid, fwrf)
     1754                        else:
     1755                            if paramId in wrfpars:
     1756                                grib_write(gid, fwrf)
     1757                except AttributeError:
     1758                    pass
     1759
     1760                grib_release(gid)
     1761                gid  =  grib_new_from_index(iid)
     1762
     1763        for f in fdict.values():
     1764            f.close()
     1765
     1766        # call for CONVERT2
     1767
     1768        if hid is not None:
     1769            pwd = os.getcwd()
     1770            os.chdir(c.inputdir)
     1771            if os.stat('fort.21').st_size == 0 and int(c.eta) == 1:
     1772                print 'Parameter 77 (etadot) is missing, most likely it is \
     1773                       not available for this type or date/time\n'
     1774                print 'Check parameters CLASS, TYPE, STREAM, START_DATE\n'
     1775                myerror(c, 'fort.21 is empty while parameter eta is set \
     1776                            to 1 in CONTROL file')
     1777
     1778            p = subprocess.check_call([os.path.expandvars(
     1779                    os.path.expanduser(c.exedir)) + '/CONVERT2'], shell=True)
     1780            os.chdir(pwd)
     1781            # create the corresponding output file fort.15
     1782            # (generated by CONVERT2)
     1783            # + fort.16 (paramId 167 and paramId 168)
     1784            fnout = c.inputdir + '/' + c.prefix
     1785            if c.maxstep > 12:
     1786                suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \
     1787                         '.{:0>3}'.format(step)
     1788            else:
     1789                suffix = cdateH[2:10]
     1790
     1791            fnout += suffix
     1792            print "outputfile = " + fnout
     1793            self.outputfilelist.append(fnout) # needed for final processing
     1794            fout = open(fnout, 'wb')
     1795            shutil.copyfileobj(open(c.inputdir + '/fort.15', 'rb'), fout)
     1796            if c.cwc == '1':
     1797                shutil.copyfileobj(open(c.inputdir + '/fort.22', 'rb'), fout)
     1798            shutil.copyfileobj(open(c.inputdir + '/flux' + cdate[0:2] +
     1799                                    suffix, 'rb'), fout)
     1800            shutil.copyfileobj(open(c.inputdir + '/fort.16', 'rb'), fout)
     1801            orolsm = glob.glob(c.inputdir +
     1802                               '/OG_OROLSM__SL.*.' + c.ppid + '*')[0]
     1803            shutil.copyfileobj(open(orolsm, 'rb'), fout)
     1804            fout.close()
     1805            if c.omega == '1':
     1806                fnout = c.outputdir + '/OMEGA'
     1807                fout  =  open(fnout, 'wb')
     1808                shutil.copyfileobj(open(c.inputdir + '/fort.25', 'rb'), fout)
     1809
     1810        try:
     1811            if c.wrf == '1':
     1812                fwrf.close()
     1813        except:
     1814            pass
     1815
     1816        grib_index_release(iid)
     1817
     1818        return
     1819
     1820
     1821    def deacc_fluxes(self, inputfiles, c):
     1822        '''
     1823        @Description:
     1824
     1825
     1826        @Input:
     1827            self: instance of EIFlexpart
     1828                The current object of the class.
     1829
     1830            inputfiles: list of strings
     1831                A list of filenames.
     1832
     1833            c: instance of class Control
     1834                Contains all the parameters of control files, which are:
     1835                DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     1836                NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL,
     1837                LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA,
     1838                ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX,
     1839                ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL,
     1840                GRIB2FLEXPART, FLEXPARTDIR
     1841                For more information about format and content of the parameter
     1842                see documentation.
     1843
     1844        @Return:
     1845            <nothing>
     1846        '''
     1847
     1848        table128 = init128(c.ecmwfdatadir +
     1849                           '/grib_templates/ecmwf_grib1_table_128')
     1850        pars = toparamId(self.params['OG_acc_SL'][0], table128)
     1851        index_keys = ["date", "time", "step"]
     1852        indexfile = c.inputdir + "/date_time_stepRange.idx"
     1853        silentremove(indexfile)  # delete, if it exists already
     1854        grib = GribTools(inputfiles.files)
     1855        iid = grib.index(index_keys=index_keys, index_file=indexfile)
     1856        print 'index done...'
     1857
     1858        # get the date, time and step values from indexfile
    10791859        index_vals = []
    10801860        for key in index_keys:
    10811861            key_vals = grib_index_get(iid,key)
    1082             print key_vals
    1083 
     1862            print key_vals
     1863            # have to sort the steps, therefore convert to int first
     1864            if key == 'step':
     1865                l = []
     1866                for k in key_vals:
     1867                    l.append(int(k))
     1868                l.sort()
     1869                # now convert back to str and overwrite key_vals
     1870                key_vals = []
     1871                for k in l:
     1872                    key_vals.append(str(k))
    10841873            index_vals.append(key_vals)
    10851874
    1086         for prod in product(*index_vals):
    1087             for i in range(len(index_keys)):
    1088                 grib_index_select(iid,index_keys[i],prod[i])
    1089 
    1090                    
    1091             gid = grib_new_from_index(iid)
    1092             hid = gid
    1093             cflextime = None
    1094             for k,f in fdict.iteritems():
    1095                 fdict[k] = open(c.inputdir+'/fort.'+k,'w')
    1096             if gid is not None:
    1097                 cdate = str(grib_get(gid, 'date'))
    1098                 time = grib_get(gid, 'time')
    1099                 type = grib_get(gid, 'type')
    1100                 step = grib_get(gid, 'step')
    1101 #                step = grib_get(gid, 'stepRange')
    1102 #               cflextime = self.getFlexpartTime(type,step, time)
    1103                 timestamp=datetime.datetime.strptime(cdate+'{:0>2}'.format(time/100),'%Y%m%d%H')
    1104                 timestamp+=datetime.timedelta(hours=int(step))
    1105 #               print gid,index_keys[i],prod[i],cdate,time,step,timestamp
    1106                
    1107                 cdateH=datetime.datetime.strftime(timestamp,'%Y%m%d%H')
    1108                 chms=datetime.datetime.strftime(timestamp,'%H%M%S')
    1109 
    1110                 if c.basetime !=None:
    1111                     slimit=datetime.datetime.strptime(c.start_date+'00','%Y%m%d%H')
    1112                     bt='23'
    1113                     if c.basetime=='00':
    1114                         bt='00'
    1115                         slimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')-datetime.timedelta(hours=12-int(c.dtime))
    1116                        
    1117                     if c.basetime=='12':
    1118                         bt='12'
    1119                         slimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')-datetime.timedelta(hours=12-int(c.dtime))
    1120                        
    1121                     elimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')
    1122                    
    1123                     if timestamp<slimit or timestamp>elimit:
    1124                         continue
    1125                    
    1126                    
    1127            
    1128             try:
    1129                 if c.wrf=='1':
    1130                     if 'olddate' not in locals():
    1131                         fwrf = open(c.outputdir+'/WRF'+cdate+'.{:0>2}'.format(time)+'.000.grb2','w')
    1132                         olddate=cdate[:]
    1133                     else:
    1134                         if cdate!=olddate:
    1135                             fwrf = open(c.outputdir+'/WRF'+cdate+'.{:0>2}'.format(time)+'.000.grb2','w')
    1136                             olddate=cdate[:]
    1137             except AttributeError:
    1138                 pass
    1139                        
    1140 #                print 'cyear '+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime
    1141 
    1142             savedfields=[]
    1143             while 1:
    1144                 if gid is None: break
    1145                 paramId = grib_get(gid, 'paramId')
    1146                 gridtype = grib_get(gid, 'gridType')
    1147                 datatype = grib_get(gid, 'dataType')
    1148                 levtype = grib_get(gid, 'typeOfLevel')
    1149                 if paramId == 133 and gridtype=='reduced_gg':
    1150 # Relative humidity (Q.grb) is used as a template only so we need the first we "meet"
    1151                     fout=open(c.inputdir+'/fort.18','w')
    1152                     grib_write(gid,fout)
    1153                     fout.close()
    1154                 elif paramId == 131 or paramId == 132:
    1155                     grib_write(gid, fdict['10'])
    1156                 elif paramId == 130:
    1157                     grib_write(gid, fdict['11'])
    1158                 elif paramId == 133 and gridtype!='reduced_gg':
    1159                     grib_write(gid, fdict['17'])
    1160                 elif paramId == 152:
    1161                     grib_write(gid, fdict['12'])
    1162                 elif paramId == 155 and gridtype=='sh':
    1163                     grib_write(gid, fdict['13'])
    1164                 elif paramId in [129,138,155] and levtype=='hybrid' and c.wrf=='1':
    1165 #                   print paramId,'not written'
    1166                     pass
    1167                 elif paramId == 246 or paramId == 247:  # cloud liquid water and ice
    1168                     if paramId==246:
    1169                         clwc=grib_get_values(gid)
    1170                     else:
    1171                         clwc+=grib_get_values(gid)
    1172                         grib_set_values(gid,clwc)
    1173 #                       grib_set(gid,'shortName','qc')
    1174                         grib_set(gid,'paramId',201031)
    1175                         grib_write(gid, fdict['22'])
    1176                    
    1177                 elif paramId == 135:
    1178                     grib_write(gid, fdict['19'])
    1179                 elif paramId == 77:
    1180                     grib_write(gid, fdict['21'])
    1181                 else:
    1182                     if paramId not in savedfields:
    1183                         grib_write(gid, fdict['16'])
    1184                         savedfields.append(paramId)
    1185                     else:
    1186                         print 'duplicate '+str(paramId)+' not written'
    1187                
    1188                 try:
    1189                     if c.wrf=='1':
    1190                         if levtype=='hybrid':
    1191                             if paramId in [129,130,131,132,133,138,155]:
    1192                                 grib_write(gid,fwrf)
    1193                         else:
    1194                             if paramId in wrfpars:
    1195                                 grib_write(gid,fwrf)
    1196                 except AttributeError:
    1197                     pass
    1198                            
    1199                            
    1200                
    1201                 grib_release(gid)
    1202                 gid = grib_new_from_index(iid)
    1203 # call for CONVERT2
    1204            
    1205             for f in fdict.values():
    1206                 f.close()
    1207 
    1208             if hid is not None:
    1209                 pwd=os.getcwd()
    1210                 os.chdir(c.inputdir)
    1211                 if os.stat('fort.21').st_size==0 and int(c.eta)==1:
    1212                     print 'Parameter 77 (etadot) is missing, most likely it is not available for this type or date/time\n'
    1213                     print 'Check parameters CLASS, TYPE, STREAM, START_DATE\n'             
    1214                     myerror(c,'fort.21 is empty while parameter eta is set to 1 in CONTROL file')
    1215                
    1216                 p=subprocess.check_call([os.path.expandvars(os.path.expanduser(c.exedir))+'/CONVERT2'],shell=True)
    1217                 os.chdir(pwd)
    1218 # create the corresponding output file  fort.15 (generated by CONVERT2) + fort.16 (paramId 167 and paramId 168) 
    1219                 fnout=c.inputdir+'/'+c.prefix
    1220                 if c.maxstep>12:
    1221                     suffix=cdate[2:8]+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step)
    1222                 else:
    1223                     suffix=cdateH[2:10]
    1224                    
    1225                 fnout+=suffix
    1226                 print "outputfile = " + fnout
    1227                 self.outputfilelist.append(fnout) # needed for final processing
    1228                 fout = open(fnout,'wb')   
    1229                 shutil.copyfileobj(open(c.inputdir+'/fort.15','rb'), fout)     
    1230                 if c.cwc=='1':
    1231                     shutil.copyfileobj(open(c.inputdir+'/fort.22','rb'), fout)     
    1232                 shutil.copyfileobj(open(c.inputdir+'/flux'+cdate[0:2]+suffix,'rb'), fout)     
    1233                 shutil.copyfileobj(open(c.inputdir+'/fort.16','rb'), fout)
    1234                 orolsm=glob.glob(c.inputdir+'/OG_OROLSM__SL.*.'+c.ppid+'*')[0]
    1235                 shutil.copyfileobj(open(orolsm,'rb'), fout)
    1236                 fout.close()
    1237                 if c.omega=='1':
    1238                     fnout=c.outputdir+'/OMEGA'
    1239                     fout = open(fnout,'wb')   
    1240                     shutil.copyfileobj(open(c.inputdir+'/fort.25','rb'), fout)   
    1241                    
    1242                
    1243                    
    1244         try:   
    1245             if c.wrf=='1':
    1246                 fwrf.close()
    1247         except:
    1248             pass
    1249 
    1250         grib_index_release(iid)
    1251 
    1252     def deacc_fluxes(self, inputfiles, c):
    1253        
    1254         table128=init128(c.ecmwfdatadir+'/grib_templates/ecmwf_grib1_table_128')
    1255         pars=toparamId(self.params['OG_acc_SL'][0],table128)
    1256         index_keys=["date","time","step"]
    1257         indexfile=c.inputdir+"/date_time_stepRange.idx"
    1258         silentremove(indexfile)
    1259         grib=GribTools(inputfiles.files)
    1260         iid=grib.index(index_keys=index_keys, index_file = indexfile)
    1261 
    1262         print 'index done...'
    1263 
    1264         index_vals = []
    1265         for key in index_keys:
    1266             key_vals = grib_index_get(iid,key)
    1267             print key_vals
    1268             if key=='step':
    1269                 l=[]
    1270                 for k in key_vals:
    1271                     l.append(int(k))
    1272                 l.sort()
    1273                 key_vals=[]
    1274                 for k in l:
    1275                     key_vals.append(str(k))
    1276 
    1277             index_vals.append(key_vals)
    1278 
    1279 
    1280         valsdict={}
    1281         svalsdict={}
    1282         stepsdict={}
    1283         for p in pars:
    1284             valsdict[str(p)]=[]
    1285             svalsdict[str(p)]=[]
    1286             stepsdict[str(p)]=[]
    1287            
    1288         for prod in product(*index_vals):
    1289             for i in range(len(index_keys)):
    1290                 grib_index_select(iid,index_keys[i],prod[i])
    1291 
    1292             #for k,f in fdict.iteritems():
    1293                 #fdict[k] = open('fort.'+k,'w')
    1294             gid = grib_new_from_index(iid)
    1295             hid = gid
    1296             cflextime = None
    1297             if gid is not None:
    1298                 cdate = grib_get(gid, 'date')
    1299                 #cyear = cdate[:4]
    1300                 #cmonth = cdate[4:6]
    1301                 #cday = cdate[6:8]
    1302                 time = grib_get(gid, 'time')
    1303                 type = grib_get(gid, 'type')
    1304                 step = grib_get(gid, 'step')
    1305                 # date+time+step-2*dtime (since interpolated value valid for step-2*dtime)
    1306                 sdate=datetime.datetime(year=cdate/10000,month=mod(cdate,10000)/100,
    1307                                         day=mod(cdate,100),hour=time/100)
    1308                 fdate=sdate+datetime.timedelta(hours=step-2*int(c.dtime))
    1309                 sdates=sdate+datetime.timedelta(hours=step)
    1310             else:
    1311                 break
    1312            
    1313             fnout=c.inputdir+'/'
    1314             if c.maxstep>12:
    1315                 fnout+='flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step-2*int(c.dtime))
    1316                 gnout=c.inputdir+'/flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step-int(c.dtime))
    1317                 hnout=c.inputdir+'/flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step)
    1318                 g=open(gnout,'w')           
    1319                 h=open(hnout,'w')           
    1320             else:
    1321                 fnout+='flux'+fdate.strftime('%Y%m%d%H')
    1322                 gnout=c.inputdir+'/flux'+(fdate+datetime.timedelta(hours=int(c.dtime))).strftime('%Y%m%d%H')
    1323                 hnout=c.inputdir+'/flux'+sdates.strftime('%Y%m%d%H')
    1324                 g=open(gnout,'w')           
    1325                 h=open(hnout,'w')           
    1326             print "outputfile = " + fnout
    1327             f=open(fnout,'w')       
    1328                        
    1329 #                print 'cyear '+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime
    1330             while 1:
    1331                 if gid is None: break
     1875        valsdict = {}
     1876        svalsdict = {}
     1877        stepsdict = {}
     1878        for p in pars:
     1879            valsdict[str(p)] = []
     1880            svalsdict[str(p)] = []
     1881            stepsdict[str(p)] = []
     1882#AP muss das hier wirklich eingerückt sein?
     1883            for prod in product(*index_vals):
     1884                for i in range(len(index_keys)):
     1885                    grib_index_select(iid, index_keys[i], prod[i])
     1886
     1887                gid = grib_new_from_index(iid)
     1888                hid = gid
     1889                cflextime = None
     1890                if gid is not None:
     1891                    cdate = grib_get(gid, 'date')
     1892                    #cyear = cdate[:4]
     1893                    #cmonth = cdate[4:6]
     1894                    #cday  = cdate[6:8]
     1895                    time = grib_get(gid, 'time')
     1896                    type = grib_get(gid, 'type')
     1897                    step = grib_get(gid, 'step')
     1898                    # date+time+step-2*dtime (since interpolated value valid for step-2*dtime)
     1899                    sdate = datetime.datetime(year=cdate / 10000,
     1900                                              month=mod(cdate, 10000) / 100,
     1901                                              day=mod(cdate, 100),
     1902                                              hour=time / 100)
     1903                    fdate = sdate + datetime.timedelta(
     1904                                        hours=step - 2 * int(c.dtime))
     1905                    sdates = sdate + datetime.timedelta(hours=step)
     1906                else:
     1907                    break
     1908
     1909            if c.maxstep > 12:
     1910                fnout = c.inputdir + '/flux' + \
     1911                    sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
     1912                    '.{:0>3}'.format(step-2*int(c.dtime))
     1913                gnout = c.inputdir + '/flux' + \
     1914                    sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100) + \
     1915                    '.{:0>3}'.format(step-int(c.dtime))
     1916                hnout = c.inputdir + '/flux' + \
     1917                    sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100) + \
     1918                    '.{:0>3}'.format(step)
     1919                g = open(gnout, 'w')
     1920                h = open(hnout, 'w')
     1921            else:
     1922                fnout = c.inputdir + '/flux' + fdate.strftime('%Y%m%d%H')
     1923                gnout = c.inputdir + '/flux' + (fdate+datetime.timedelta(
     1924                    hours = int(c.dtime))).strftime('%Y%m%d%H')
     1925                hnout = c.inputdir + '/flux' + sdates.strftime('%Y%m%d%H')
     1926                g = open(gnout, 'w')
     1927                h = open(hnout, 'w')
     1928            print "outputfile = " + fnout
     1929            f = open(fnout, 'w')
     1930
     1931            while 1:
     1932                if gid is None:
     1933                    break
    13321934                cparamId = str(grib_get(gid, 'paramId'))
    13331935                step = grib_get(gid, 'step')
    13341936                atime = grib_get(gid, 'time')
    1335                 ni=grib_get(gid, 'Ni')
    1336                 nj=grib_get(gid, 'Nj')
    1337                 if cparamId in valsdict.keys():
    1338                     values = grib_get_values(gid)
    1339                     vdp=valsdict[cparamId]
    1340                     svdp=svalsdict[cparamId]
    1341                     sd=stepsdict[cparamId]
    1342                    
    1343                     if cparamId=='142' or cparamId=='143':
    1344                         fak=1./1000.
    1345                     else:
    1346                         fak=3600.
    1347                    
    1348                     values=(reshape(values,(nj,ni))).flatten()/fak
    1349                     vdp.append(values[:]) # save the accumulated values
    1350                     if step<=int(c.dtime):
    1351                         svdp.append(values[:]/int(c.dtime))
    1352                     else:
    1353                         svdp.append((vdp[-1]-vdp[-2])/int(c.dtime))
    1354                    
    1355                     print cparamId,atime,step,len(values),values[0],std(values)
    1356                     #svdp.append(values[:]) # save the 1-hourly or 3-hourly specific values
    1357                     sd.append(step)
    1358                     if len(svdp)>=3:
    1359                         if len(svdp)>3:
    1360                             if cparamId=='142' or cparamId=='143':
    1361                                 values=darain(svdp)
    1362                             else:
    1363                                 values=dapoly(svdp)
    1364                                                    
    1365                             if not (step==c.maxstep and c.maxstep>12 or sdates==elimit):
    1366                                 vdp.pop(0)
    1367                                 svdp.pop(0)
    1368                         else:
    1369                             if c.maxstep>12:
    1370                                 values=svdp[1]
    1371                             else:
    1372                                 values=svdp[0]
    1373                            
    1374                         grib_set_values(gid, values)
    1375                         if c.maxstep>12:
    1376                             grib_set(gid,'step',max(0,step-2*int(c.dtime)))
    1377                         else:
    1378                             grib_set(gid,'step',0)
    1379                             grib_set(gid,'time',fdate.hour*100)
    1380                             grib_set(gid,'date',fdate.year*10000+fdate.month*100+fdate.day)
    1381                         grib_write(gid, f)
    1382                        
    1383                         if c.basetime is not None:
    1384                             elimit=datetime.datetime.strptime(c.end_date+c.basetime,'%Y%m%d%H')
    1385                         else:
    1386                             elimit=sdate+datetime.timedelta(2*int(c.dtime))
    1387                        
    1388                         # squeeze out information of last two steps contained in svdp
    1389 #                       if step+int(c.dtime)==c.maxstep and c.maxstep>12 or sdates+datetime.timedelta(hours=int(c.dtime))>=elimit:
    1390 #                       Note that svdp[0] has not been popped in this case
    1391                            
    1392                    
    1393                         if step==c.maxstep and c.maxstep>12 or sdates==elimit:
    1394                             values=svdp[3]
    1395                             grib_set_values(gid, values)
    1396                             grib_set(gid,'step',0)
    1397                             truedatetime=fdate+datetime.timedelta(hours=2*int(c.dtime))
    1398                             grib_set(gid,'time',truedatetime.hour*100)
    1399                             grib_set(gid,'date',truedatetime.year*10000+truedatetime.month*100+truedatetime.day)
    1400                             grib_write(gid, h)
    1401                            
    1402                             #values=(svdp[1]+svdp[2])/2.
    1403                             if cparamId=='142' or cparamId=='143':
    1404                                 values=darain(list(reversed(svdp)))
    1405                             else:
    1406                                 values=dapoly(list(reversed(svdp)))
    1407                            
    1408                             grib_set(gid,'step',0)
    1409                             truedatetime=fdate+datetime.timedelta(hours=int(c.dtime))
    1410                             grib_set(gid,'time',truedatetime.hour*100)
    1411                             grib_set(gid,'date',truedatetime.year*10000+truedatetime.month*100+truedatetime.day)
    1412                             grib_set_values(gid, values)
    1413                             grib_write(gid, g)
    1414                        
    1415                     grib_release(gid)
    1416                                        
    1417                 gid = grib_new_from_index(iid)
    1418                
    1419             f.close()
    1420             g.close()
    1421             h.close()
    1422 
    1423 
    1424         grib_index_release(iid)
    1425 
     1937                ni = grib_get(gid, 'Ni')
     1938                nj = grib_get(gid, 'Nj')
     1939                if cparamId in valsdict.keys():
     1940                    values = grib_get_values(gid)
     1941                    vdp = valsdict[cparamId]
     1942                    svdp = svalsdict[cparamId]
     1943                    sd = stepsdict[cparamId]
     1944
     1945                    if cparamId == '142' or cparamId == '143':
     1946                        fak = 1./1000.
     1947                    else:
     1948                        fak = 3600.
     1949
     1950                    values = (reshape(values, (nj, ni))).flatten() / fak
     1951                    vdp.append(values[:])  # save the accumulated values
     1952                    if step <= int(c.dtime):
     1953                        svdp.append(values[:] / int(c.dtime))
     1954                    else:
     1955                        svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
     1956
     1957                    print cparamId, atime, step, len(values), \
     1958                          values[0], std(values)
     1959                    # save the 1/3-hourly or specific values
     1960                    # svdp.append(values[:])
     1961                    sd.append(step)
     1962                    if len(svdp) >= 3:
     1963                        if len(svdp) > 3:
     1964                            if cparamId == '142' or cparamId == '143':
     1965                                values = darain(svdp)
     1966                            else:
     1967                                values = dapoly(svdp)
     1968
     1969                            if not (step == c.maxstep and c.maxstep > 12 \
     1970                                    or sdates == elimit):
     1971                                vdp.pop(0)
     1972                                svdp.pop(0)
     1973                        else:
     1974                            if c.maxstep > 12:
     1975                                values = svdp[1]
     1976                            else:
     1977                                values = svdp[0]
     1978
     1979                        grib_set_values(gid, values)
     1980                        if c.maxstep > 12:
     1981                            grib_set(gid, 'step', max(0, step-2*int(c.dtime)))
     1982                        else:
     1983                            grib_set(gid, 'step', 0)
     1984                            grib_set(gid, 'time', fdate.hour*100)
     1985                            grib_set(gid, 'date', fdate.year*10000 +
     1986                                     fdate.month*100+fdate.day)
     1987                        grib_write(gid, f)
     1988
     1989                        if c.basetime is not None:
     1990                            elimit = datetime.datetime.strptime(c.end_date +
     1991                                                                c.basetime,
     1992                                                                '%Y%m%d%H')
     1993                        else:
     1994                            elimit = sdate + datetime.timedelta(2*int(c.dtime))
     1995
     1996                        # squeeze out information of last two steps contained
     1997                        # in svdp
     1998                        # if step+int(c.dtime) == c.maxstep and c.maxstep>12
     1999                        # or sdates+datetime.timedelta(hours = int(c.dtime))
     2000                        # >= elimit:
     2001                        # Note that svdp[0] has not been popped in this case
     2002
     2003                        if (step == c.maxstep and c.maxstep > 12
     2004                            or sdates == elimit):
     2005                            values = svdp[3]
     2006                            grib_set_values(gid, values)
     2007                            grib_set(gid, 'step', 0)
     2008                            truedatetime = fdate + datetime.timedelta(
     2009                                     hours=2*int(c.dtime))
     2010                            grib_set(gid, 'time', truedatetime.hour * 100)
     2011                            grib_set(gid, 'date', truedatetime.year * 10000 +
     2012                                     truedatetime.month * 100 +
     2013                                     truedatetime.day)
     2014                            grib_write(gid, h)
     2015
     2016                            #values = (svdp[1]+svdp[2])/2.
     2017                            if cparamId == '142' or cparamId == '143':
     2018                                values = darain(list(reversed(svdp)))
     2019                            else:
     2020                                values = dapoly(list(reversed(svdp)))
     2021
     2022                            grib_set(gid, 'step',0)
     2023                            truedatetime = fdate + datetime.timedelta(
     2024                                     hours=int(c.dtime))
     2025                            grib_set(gid, 'time', truedatetime.hour * 100)
     2026                            grib_set(gid, 'date', truedatetime.year * 10000 +
     2027                                     truedatetime.month * 100 +
     2028                                     truedatetime.day)
     2029                            grib_set_values(gid, values)
     2030                            grib_write(gid, g)
     2031
     2032                    grib_release(gid)
     2033
     2034            gid = grib_new_from_index(iid)
     2035
     2036            f.close()
     2037            g.close()
     2038            h.close()
     2039
     2040            grib_index_release(iid)
     2041
     2042        return
     2043
  • python/GribTools.py

    rd69b677 r64cf353  
    1 #
    2 # (C) Copyright 2014 UIO.
    3 #
    4 # This software is licensed under the terms of the Apache Licence Version 2.0
    5 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
    6 #
    7 #
    8 # Creation: July 2014 - Anne Fouilloux - University of Oslo
    9 #
    10 #
    11 
     1#!/usr/bin/env python
     2# -*- coding: utf-8 -*-
     3#************************************************************************
     4# TODO AP
     5#AP
     6# - GribTools name möglicherweise etwas verwirrend.
     7# - change self.filename in self.filenames!!!
     8# -
     9#************************************************************************
     10"""
     11@Author: Anne Fouilloux (University of Oslo)
     12
     13@Date: July 2014
     14
     15@ChangeHistory:
     16   February 2018 - Anne Philipp (University of Vienna):
     17        - applied PEP8 style guide
     18        - added documentation
     19        - changed some naming
     20
     21@License:
     22    (C) Copyright 2014 UIO.
     23
     24    This software is licensed under the terms of the Apache Licence Version 2.0
     25    which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
     26
     27@Requirements:
     28    - A standard python 2.6 or 2.7 installation
     29    - dateutils
     30    - matplotlib (optional, for debugging)
     31    - ECMWF specific packages, all available from https://software.ecmwf.int/
     32        ECMWF WebMARS, gribAPI with python enabled, emoslib and
     33        ecaccess web toolkit
     34
     35@Description:
     36    Further documentation may be obtained from www.flexpart.eu.
     37
     38    Functionality provided:
     39        Prepare input 3D-wind fields in hybrid coordinates +
     40        surface fields for FLEXPART runs
     41"""
     42# ------------------------------------------------------------------------------
     43# MODULES
     44# ------------------------------------------------------------------------------
    1245from gribapi import *
    1346import traceback
    14 import sys,os
    15 
    16 
    17 ##############################################################
    18 #  GRIB utilities
    19 ##############################################################
     47import sys, os
     48
     49# ------------------------------------------------------------------------------
     50# CLASS
     51# ------------------------------------------------------------------------------
    2052class GribTools:
    21     'class for GRIB API with new methods'
    22     def __init__(self,filename):
    23         self.filename=filename
    24    
    25 # get keyvalues for a given list of keynames
    26 # a where statment can be given (list of key and list of values)
    27 
    28     def getkeys(self,keynames,wherekeynames=[],wherekeyvalues=[]):
    29         fileid=open(self.filename,'r')
    30 
    31         return_list=[]
    32        
     53    '''
     54    Class for GRIB utilities (new methods) based on GRIB API
     55    '''
     56    # --------------------------------------------------------------------------
     57    # FUNCTIONS
     58    # --------------------------------------------------------------------------
     59    def __init__(self, filenames):
     60        '''
     61        @Description:
     62            Initialise an object of GribTools and assign a list
     63            of filenames.
     64
     65        @Input:
     66            filenames: list of strings
     67                A list of filenames.
     68
     69        @Return:
     70            <nothing>
     71        '''
     72
     73        self.filename = filenames
     74
     75        return
     76
     77
     78    def getkeys(self, keynames, wherekeynames=[], wherekeyvalues=[]):
     79        '''
     80        @Description:
     81            get keyvalues for a given list of keynames
     82            a where statement can be given (list of key and list of values)
     83
     84        @Input:
     85            keynames: list of strings
     86                List of keynames.
     87
     88            wherekeynames: list of ???
     89                Default value is an empty list.
     90
     91            wherekeyvalues: list of ???
     92                Default value is an empty list.
     93
     94        @Return:
     95            return_list: list of strings
     96                List of keyvalues for given keynames.
     97        '''
     98
     99        fileid = open(self.filename, 'r')
     100
     101        return_list = []
     102
    33103        while 1:
    34             gid_in = grib_new_from_file(fileid) 
    35             if gid_in is None: break
    36             select=True
    37             i=0
    38             if len(wherekeynames) != len(wherekeyvalues): raise Exception("Give a value for each keyname!")
    39 
     104            gid_in = grib_new_from_file(fileid)
     105
     106            if gid_in is None:
     107                break
     108
     109            if len(wherekeynames) != len(wherekeyvalues):
     110                raise Exception("Number of key values and key names must be \
     111                                 the same. Give a value for each keyname!")
     112
     113            select = True
     114            i = 0
    40115            for wherekey in wherekeynames:
    41                 if not grib_is_defined(gid_in, wherekey): raise Exception("where Key was not defined")
    42                 select=select and (str(wherekeyvalues[i])==str(grib_get(gid_in, wherekey)))
    43                 i=i+1
     116                if not grib_is_defined(gid_in, wherekey):
     117                    raise Exception("where key was not defined")
     118
     119                select = (select and (str(wherekeyvalues[i]) ==
     120                                      str(grib_get(gid_in, wherekey))))
     121                i += 1
     122
    44123            if select:
    45124                llist = []
     
    47126                    llist.extend([str(grib_get(gid_in, key))])
    48127                return_list.append(llist)
     128
    49129            grib_release(gid_in)
     130
    50131        fileid.close()
     132
    51133        return return_list
    52      
    53 # set keyvalues for a given list of keynames
    54 # a where statment can be given (list of key and list of values)
    55 # an input file must be given as an input for reading grib messages
    56 # note that by default all messages are written out
    57 # if you want to get only those meeting the where statement, use
    58 # strict=true
    59     def setkeys(self,fromfile,keynames,keyvalues, wherekeynames=[],wherekeyvalues=[], strict=False, filemode='w'):
    60         fout=open(self.filename,filemode)
    61         fin=open(fromfile)
    62        
     134
     135
     136    def setkeys(self, fromfile, keynames, keyvalues, wherekeynames=[],
     137                wherekeyvalues=[], strict=False, filemode='w'):
     138        '''
     139        @Description:
     140            Opens the file to read the grib messages and then write
     141            them to a new output file. By default all messages are
     142            written out. Also, the keyvalues of the passed list of
     143            keynames are set or only those meeting the where statement.
     144            (list of key and list of values).
     145
     146        @Input:
     147            fromfile: string
     148                Filename of the input file to read the grib messages from.
     149
     150            keynames: list of ???
     151                List of keynames. Default is an empty list.
     152
     153            keyvalues: list of ???
     154                List of keynames. Default is an empty list.
     155
     156            wherekeynames: list of ???
     157                Default value is an empty list.
     158
     159            wherekeyvalues: list of ???
     160                Default value is an empty list.
     161
     162            strict: boolean
     163                Decides if everything from keynames and keyvalues
     164                is written out the grib file (False) or only those
     165                meeting the where statement (True). Default is False.
     166
     167            filemode:
     168                Sets the mode for the output file. Default is "w".
     169
     170        @Return:
     171            <nothing>
     172
     173        '''
     174        fout = open(self.filename, filemode)
     175        fin = open(fromfile)
     176
    63177        while 1:
    64             gid_in = grib_new_from_file(fin) 
    65             if gid_in is None: break
    66            
    67             select=True
    68             i=0
    69             if len(wherekeynames) != len(wherekeyvalues): raise Exception("Give a value for each keyname!")
    70 
     178            gid_in = grib_new_from_file(fin)
     179
     180            if gid_in is None:
     181                break
     182
     183            if len(wherekeynames) != len(wherekeyvalues):
     184                raise Exception("Give a value for each keyname!")
     185
     186            select = True
     187            i = 0
    71188            for wherekey in wherekeynames:
    72                 if not grib_is_defined(gid_in, wherekey): raise Exception("where Key was not defined")
    73                 select=select and (str(wherekeyvalues[i])==str(grib_get(gid_in, wherekey)))
    74                 i=i+1
     189                if not grib_is_defined(gid_in, wherekey):
     190                    raise Exception("where Key was not defined")
     191
     192                select = (select and (str(wherekeyvalues[i]) ==
     193                                      str(grib_get(gid_in, wherekey))))
     194                i += 1
     195
     196#AP is it secured that the order of keynames is equal to keyvalues?
    75197            if select:
    76                 i=0
     198                i = 0
    77199                for key in keynames:
    78200                    grib_set(gid_in, key, keyvalues[i])
    79                     i=i+1
     201                    i += 1
     202
     203#AP this is a redundant code section
     204# delete the if/else :
     205#
     206#           grib_write(gid_in, fout)
     207#
    80208            if strict:
    81209                if select:
    82                     grib_write(gid_in,fout)
     210                    grib_write(gid_in, fout)
    83211            else:
    84                 grib_write(gid_in,fout)
     212                grib_write(gid_in, fout)
     213#AP end
    85214            grib_release(gid_in)
     215
    86216        fin.close()
    87217        fout.close()
    88218
    89 # Add the content of a grib file but only messages
    90 # corresponding to keys/values       
    91 # if selectWhere is False select fields that are different from keynames/keyvalues
    92     def copy(self,filename_in, selectWhere=True, keynames=[], keyvalues=[],filemode='w'):
    93         fin=open(filename_in)
    94         fout=open(self.filename,filemode)
    95        
     219        return
     220
     221    def copy(self, filename_in, selectWhere=True,
     222             keynames=[], keyvalues=[], filemode='w'):
     223        '''
     224        Add the content of another input grib file to the objects file but
     225        only messages corresponding to keys/values passed to the function.
     226        The selectWhere switch decides if to copy the keys equal to (True) or
     227        different to (False) the keynames/keyvalues list passed to the function.
     228
     229        @Input:
     230            filename_in: string
     231                Filename of the input file to read the grib messages from.
     232
     233            selectWhere: boolean
     234                Decides if to copy the keynames and values equal to (True) or
     235                different to (False) the keynames/keyvalues list passed to the
     236                function. Default is True.
     237
     238            keynames: list of ???
     239                List of keynames. Default is an empty list.
     240
     241            keyvalues: list of ???
     242                List of keynames. Default is an empty list.
     243
     244            filemode:
     245                Sets the mode for the output file. Default is "w".
     246
     247        @Return:
     248            <nothing>
     249        '''
     250
     251        fin = open(filename_in)
     252        fout = open(self.filename, filemode)
     253
    96254        while 1:
    97             gid_in = grib_new_from_file(fin) 
    98             if gid_in is None: break
    99            
    100             select=True
    101             i=0
    102             if len(keynames) != len(keyvalues): raise Exception("Give a value for each keyname!")
    103 
     255            gid_in = grib_new_from_file(fin)
     256
     257            if gid_in is None:
     258                break
     259
     260            if len(keynames) != len(keyvalues):
     261                raise Exception("Give a value for each keyname!")
     262
     263            select = True
     264            i = 0
    104265            for key in keynames:
    105                 if not grib_is_defined(gid_in, key): raise Exception("Key was not defined")
    106                
     266                if not grib_is_defined(gid_in, key):
     267                    raise Exception("Key was not defined")
     268
    107269                if selectWhere:
    108                     select=select and (str(keyvalues[i])==str(grib_get(gid_in, key)))
     270                    select = (select and (str(keyvalues[i]) ==
     271                                          str(grib_get(gid_in, key))))
    109272                else:
    110                     select=select and (str(keyvalues[i])!=str(grib_get(gid_in, key)))
    111                 i=i+1
     273                    select = (select and (str(keyvalues[i]) !=
     274                                          str(grib_get(gid_in, key))))
     275                i += 1
     276
    112277            if select:
    113                 grib_write(gid_in,fout)
     278                grib_write(gid_in, fout)
     279
    114280            grib_release(gid_in)
     281
    115282        fin.close()
    116283        fout.close()
    117284
    118 # Create index from a list of files if it does not exist or read it
    119     def index(self,index_keys=["mars"], index_file = "my.idx"):
    120         print "index to be done"
     285        return
     286
     287    def index(self, index_keys=["mars"], index_file="my.idx"):
     288        '''
     289        @Description:
     290            Create index from a list of files if it does not exist or
     291            read an index file.
     292
     293        @Input:
     294            index_keys: list of ???
     295                Default is a list with a single entry string "mars".
     296
     297            index_file: string
     298                Filename where the indices are stored.
     299                Default is "my.idx".
     300
     301        @Return:
     302            iid: integer
     303                Grib index id.
     304        '''
     305        print("... index will be done")
    121306        self.iid = None
    122  
     307
    123308        if (os.path.exists(index_file)):
    124309            self.iid = grib_index_read(index_file)
    125             print "Use existing index file: %s "%(index_file)
     310            print("Use existing index file: %s " % (index_file))
    126311        else:
     312#AP does the for loop overwrite the iid all the time?
    127313            for file in self.filename:
    128                 print "Inputfile: %s "%(file)
    129                 if self.iid == None:
    130                     self.iid = grib_index_new_from_file(file,index_keys)
     314                print("Inputfile: %s " % (file))
     315                if self.iid is None:
     316                    self.iid = grib_index_new_from_file(file, index_keys)
    131317                else:
    132                      grib_index_add_file(self.iid,file)
    133 
    134             if self.iid != None:
    135               grib_index_write(self.iid,index_file)
    136         return self.iid
    137 
    138 
    139        
    140 
    141        
     318                    grib_index_add_file(self.iid, file)
     319#AP or does the if has to be in the for loop?
     320#AP would make more sense?
     321            if self.iid is not None:
     322                grib_index_write(self.iid, index_file)
     323
     324        return self.iid
     325
     326
     327
     328
     329
  • python/UIOTools.py

    r6180177 r64cf353  
     1#!/usr/bin/env python
     2# -*- coding: utf-8 -*-
    13#************************************************************************
    24# TODO AP
    3 #
     5#AP
    46# - File name und Klassenname gleichsetzen?
    57# - checken welche regelmässigen methoden auf diese Files noch angewendet werden
     
    1315
    1416@ChangeHistory:
    15    Anne Philipp (University of Vienna) - February 2018:
    16        Added documentation and applied pep8 style guides
     17   February 2018 - Anne Philipp (University of Vienna):
     18        - applied PEP8 style guide
     19        - added documentation
    1720
    1821@License:
     
    2124    This software is licensed under the terms of the Apache Licence Version 2.0
    2225    which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
     26
     27@Requirements:
     28    - A standard python 2.6 or 2.7 installation
     29    - dateutils
     30    - matplotlib (optional, for debugging)
     31    - ECMWF specific packages, all available from https://software.ecmwf.int/
     32        ECMWF WebMARS, gribAPI with python enabled, emoslib and
     33        ecaccess web toolkit
     34
     35@Description:
     36    Further documentation may be obtained from www.flexpart.eu.
     37
     38    Functionality provided:
     39        Prepare input 3D-wind fields in hybrid coordinates +
     40        surface fields for FLEXPART runs
    2341"""
    24 
    25 
    2642# ------------------------------------------------------------------------------
    2743# MODULES
     
    5470            suffix: list of strings
    5571                Types of files to manipulate such as
    56                 ['.grib', 'grb', 'grib1', 'grib2', 'grb1','grb2']
     72                ['grib', 'grb', 'grib1', 'grib2', 'grb1', 'grb2']
    5773
    5874        @Return:
     
    8399            <nothing>
    84100        '''
     101#AP pathname zu path ändern
     102#AP is it possible for each possible file extension ? mabye regexx?
    85103        # Get the absolute path of the pathname parameter
    86104        pathname = os.path.abspath(pathname)
  • python/getMARSdata.py

    rd69b677 r64cf353  
    22#
    33# This software is licensed under the terms of the Apache Licence Version 2.0
    4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 
    5 # 
     4# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
     5#
    66# Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs
    77#
    88# Creation: October  2014 - Anne Fouilloux - University of Oslo
    9 # Extension November 2015 - Leopold Haimberger - University of Vienna for: 
     9# Extension November 2015 - Leopold Haimberger - University of Vienna for:
    1010# - using the WebAPI also for general MARS retrievals
    1111# - job submission on ecgate and cca
     
    1515# - conversion into GRIB2
    1616# - conversion into .fp format for faster execution of FLEXPART
    17 # 
     17#
    1818#
    1919# Further documentation may be obtained from www.flexpart.eu
    20 # 
    21 # Requirements: 
     20#
     21# Requirements:
    2222# in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed
    2323# ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/
     
    2626#
    2727# Get MARS GRIB fields from ECMWF for FLEXPART
    28 # 
     28#
    2929
    3030#import socket
     
    5151
    5252from FlexpartTools import MARSretrieval, EIFlexpart, silentremove, \
    53                           Control,myerror,normalexit, interpret_args_and_control
     53                          Control, myerror, normalexit, \
     54                          interpret_args_and_control
    5455
    5556
    56 def getMARSdata(args,c):
    57    
    58        
     57def getMARSdata(args, c):
     58
     59
    5960    if not os.path.exists(c.inputdir):
    6061        os.makedirs(c.inputdir)
    61     print "start date %s "%(c.start_date)
    62     print "end date %s "%(c.end_date)
     62    print("start date %s " % (c.start_date))
     63    print("end date %s " % (c.end_date))
    6364
    64    
    65 #    server = ECMWFDataServer()
    6665    if ecapi:
    6766        server = ecmwfapi.ECMWFService("mars")
     
    6968        server = False
    7069
    71     c.ecapi=ecapi
    72     print 'ecapi:',c.ecapi
     70    c.ecapi = ecapi
     71    print 'ecapi:', c.ecapi
     72
    7373# Retrieve ERA interim data for running flexpart
     74#AP change this variant to correct format conversion with datetime
     75#AP import datetime and timedelta explicitly
     76    syear = int(c.start_date[:4])
     77    smonth = int(c.start_date[4:6])
     78    sday = int(c.start_date[6:])
     79    start = datetime.date(year=syear, month=smonth, day=sday)
     80    startm1 = start - datetime.timedelta(days=1)
     81    if c.basetime == '00':
     82        start = startm1
     83    eyear = int(c.end_date[:4])
     84    emonth = int(c.end_date[4:6])
     85    eday = int(c.end_date[6:])
     86    end = datetime.date(year=eyear, month=emonth, day=eday)
     87    if c.basetime == '00' or c.basetime == '12':
     88        endp1 = end + datetime.timedelta(days=1)
     89    else:
     90        endp1 = end + datetime.timedelta(days=2)
    7491
    75     syear=int(c.start_date[:4])
    76     smonth=int(c.start_date[4:6])
    77     sday=int(c.start_date[6:])
    78     start = datetime.date( year = syear, month = smonth, day = sday )
    79     startm1=start- datetime.timedelta(days=1)
    80     if c.basetime=='00':
    81         start=startm1
    82     eyear=int(c.end_date[:4])
    83     emonth=int(c.end_date[4:6])
    84     eday=int(c.end_date[6:])
    85 
    86     end = datetime.date( year = eyear, month = emonth, day = eday )
    87     if c.basetime=='00' or c.basetime=='12':
    88         endp1=end+ datetime.timedelta(days=1)
    89     else:
    90         endp1=end+ datetime.timedelta(days=2)
    91    
    92     datechunk=datetime.timedelta(days=int(c.date_chunk))
    93     print 'removing content of '+c.inputdir
    94     tobecleaned=glob.glob(c.inputdir+'/*_acc_*.'+str(os.getppid())+'.*.grb')
     92    datechunk = datetime.timedelta(days=int(c.date_chunk))
     93    print 'removing content of ' + c.inputdir
     94    tobecleaned = glob.glob(c.inputdir + '/*_acc_*.' + str(os.getppid()) + '.*.grb')
    9595    for f in tobecleaned:
    9696        os.remove(f)
    97        
     97
    9898    times=None
    9999    if c.maxstep<24:
     
    103103                flexpart = EIFlexpart(c,fluxes=True)
    104104                if day+datechunk-datetime.timedelta(days=1)<endp1:
    105                     dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d") 
     105                    dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d")
    106106                else:
    107                     dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d") 
    108                    
     107                    dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d")
     108
    109109                print "retrieve " + dates + " in dir " + c.inputdir
    110110                try:
     
    112112                except IOError:
    113113                    myerror(c,'MARS request failed')
    114                    
     114
    115115                day+=datechunk
    116116    else:
     
    120120                flexpart = EIFlexpart(c,fluxes=True)
    121121                if day+datechunk-datetime.timedelta(days=1)<end:
    122                     dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d") 
     122                    dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d")
    123123                else:
    124                     dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d") 
    125                    
     124                    dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d")
     125
    126126                print "retrieve " + dates + " in dir " + c.inputdir
    127127                flexpart.retrieve(server, dates, times, c.inputdir)
    128128                day+=datechunk
    129        
     129
    130130
    131131
     
    136136    times=None
    137137    while day<=end:
    138        
     138
    139139               # we need to retrieve MARS data for this period (maximum one month)
    140140            flexpart = EIFlexpart(c)
    141141            if day+datechunk-datetime.timedelta(days=1)<end:
    142                 dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d") 
     142                dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d")
    143143            else:
    144                 dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d") 
     144                dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d")
    145145            print "retrieve " + dates + " in dir " + c.inputdir
    146            
     146
    147147            flexpart.retrieve(server, dates, times, c.inputdir)
    148148            day+=datechunk
    149    
     149
    150150
    151151if __name__ == "__main__":
    152    
     152
    153153    args,c=interpret_args_and_control()
    154154    getMARSdata(args,c)
  • python/install.py

    ra4b6cef r64cf353  
    11#!/usr/bin/env python
    2 #
    3 # This software is licensed under the terms of the Apache Licence Version 2.0
    4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
    5 #
    6 # Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs
    7 #
    8 # Creation: October  2014 - Anne Fouilloux - University of Oslo
    9 # Extension November 2015 - Leopold Haimberger - University of Vienna for:
    10 # - using the WebAPI also for general MARS retrievals
    11 # - job submission on ecgate and cca
    12 # - job templates suitable for twice daily operational dissemination
    13 # - dividing retrievals of longer periods into digestable chunks
    14 # - retrieve also longer term forecasts, not only analyses and short term forecast data
    15 # - conversion into GRIB2
    16 # - conversion into .fp format for faster execution of FLEXPART
    17 #
    18 #
    19 # Further documentation may be obtained from www.flexpart.eu
    20 #
    21 # Requirements:
    22 # in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed
    23 # ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/
    24 # dateutils
    25 # matplotlib (optional, for debugging)
    26 #
    27 #
     2# -*- coding: utf-8 -*-
     3#************************************************************************
     4# TODO AP
     5#AP
     6# - Functionality Provided is not correct for this file
     7# - localpythonpath should not be set in module load section!
     8# - create a class Installation and divide installation in 3 subdefs for
     9#   ecgate, local and cca seperatly
     10# - Change History ist nicht angepasst ans File! Original geben lassen
     11#************************************************************************
     12"""
     13@Author: Anne Fouilloux (University of Oslo)
     14
     15@Date: October 2014
     16
     17@ChangeHistory:
     18    November 2015 - Leopold Haimberger (University of Vienna):
     19        - using the WebAPI also for general MARS retrievals
     20        - job submission on ecgate and cca
     21        - job templates suitable for twice daily operational dissemination
     22        - dividing retrievals of longer periods into digestable chunks
     23        - retrieve also longer term forecasts, not only analyses and
     24          short term forecast data
     25        - conversion into GRIB2
     26        - conversion into .fp format for faster execution of FLEXPART
     27
     28    February 2018 - Anne Philipp (University of Vienna):
     29        - applied PEP8 style guide
     30        - added documentation
     31
     32@License:
     33    (C) Copyright 2014 UIO.
     34
     35    This software is licensed under the terms of the Apache Licence Version 2.0
     36    which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
     37
     38@Requirements:
     39    - A standard python 2.6 or 2.7 installation
     40    - dateutils
     41    - matplotlib (optional, for debugging)
     42    - ECMWF specific packages, all available from https://software.ecmwf.int/
     43        ECMWF WebMARS, gribAPI with python enabled, emoslib and
     44        ecaccess web toolkit
     45
     46@Description:
     47    Further documentation may be obtained from www.flexpart.eu.
     48
     49    Functionality provided:
     50        Prepare input 3D-wind fields in hybrid coordinates +
     51        surface fields for FLEXPART runs
     52"""
     53# ------------------------------------------------------------------------------
     54# MODULES
     55# ------------------------------------------------------------------------------
    2856import calendar
    2957import shutil
     
    4472from prepareFLEXPART import prepareFLEXPART
    4573
    46 
     74# ------------------------------------------------------------------------------
     75# FUNCTIONS
     76# ------------------------------------------------------------------------------
    4777def main():
    4878    '''
    4979    '''
    50 #    calledfromdir = os.getcwd()
    5180    os.chdir(localpythonpath)
    5281    args, c = install_args_and_control()
    53 #    if c.outputdir[0]!='/':
    54 #    c.outputdir=os.path.join(calledfromdir,c.outputdir)
    55 #    c.inputdir=c.outputdir
    5682    if args.install_target is not None:
    5783        install_via_gateway(c, args.install_target)
     
    230256        print('You should get an email with subject flexcompile \
    231257                within the next few minutes')
     258
    232259    else:
    233260        print('ERROR: unknown installation target ', target)
  • python/plot_retrieved.py

    rd69b677 r64cf353  
    11#!/usr/bin/env python
    22# This software is licensed under the terms of the Apache Licence Version 2.0
    3 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 
    4 # 
     3# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
     4#
    55# Functionality provided: Simple tool for creating maps and time series of retrieved fields.
    6 # 
    7 # Requirements: 
     6#
     7# Requirements:
    88# in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed
    99# ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/
     
    1919if localpythonpath not in sys.path:
    2020    sys.path.append(localpythonpath)
    21    
     21
    2222from matplotlib.pylab import *
    2323import matplotlib.patches as mpatches
     
    4343    c.levels=asarray(c.levels,dtype='int')
    4444    c.area=asarray(c.area)
    45    
     45
    4646    index_keys=["date","time","step"]
    4747    indexfile=c.inputdir+"/date_time_stepRange.idx"
     
    6767
    6868        index_vals.append(key_vals)
    69    
     69
    7070    fdict=dict()
    7171    fmeta=dict()
     
    101101                fdatetimestep=fdatetime+datetime.timedelta(hours=step)
    102102                if len(fstamp)==0:
    103                     fstamp[key].append(fdatetimestamp)   
     103                    fstamp[key].append(fdatetimestamp)
    104104                    fmeta[key].append((paramId,parameterName,gtype,fdatetime,gtime,step,level))
    105105                    fdict[key].append(flipud(reshape(grib_get_values(gid),[gdict['Nj'],gdict['Ni']])))
     
    118118                        fmeta[key].append((paramId,parameterName,gtype,fdatetime,gtime,step,level))
    119119                        fdict[key].append(flipud(reshape(grib_get_values(gid),[gdict['Nj'],gdict['Ni']])))
    120                    
     120
    121121
    122122            grib_release(gid)
    123123            gid = grib_new_from_index(iid)
    124            
     124
    125125    for k in fdict.keys():
    126126        fml=fmeta[k]
    127127        fdl=fdict[k]
    128        
     128
    129129        for fd,fm in zip(fdl,fml):
    130130            ftitle=fm[1]+' {} '.format(fm[-1])+datetime.datetime.strftime(fm[3],'%Y%m%d%H')+' '+stats(fd)
    131131            pname='_'.join(fm[1].split())+'_{}_'.format(fm[-1])+datetime.datetime.strftime(fm[3],'%Y%m%d%H')+'.{:0>3}'.format(fm[5])
    132132            plotmap(fd, fm,gdict,ftitle,pname+'.eps')
    133        
     133
    134134    for k in fdict.keys():
    135135        fml=fmeta[k]
     
    143143            lat=-20
    144144            lon=20
    145             plottimeseries(fdl,fml,fsl,lat,lon, gdict, ftitle, pname+'.eps')   
    146    
     145            plottimeseries(fdl,fml,fsl,lat,lon, gdict, ftitle, pname+'.eps')
     146
    147147def plottimeseries(flist,fmetalist,ftimestamps,lat,lon,gdict,ftitle,filename):
    148148    t1=time.time()
     
    158158    plt.close(f)
    159159    print time.time()-t1,'s'
    160    
     160
    161161def plotmap(flist,fmetalist,gdict,ftitle,filename):
    162162    t1=time.time()
    163163    f=plt.figure(figsize=(12,6.7))
    164     mbaxes = f.add_axes([0.05, 0.15, 0.8, 0.7]) 
     164    mbaxes = f.add_axes([0.05, 0.15, 0.8, 0.7])
    165165    m =Basemap(llcrnrlon=-180.,llcrnrlat=-90.,urcrnrlon=180,urcrnrlat=90.)
    166166    #if bw==0 :
     
    168168    #else:
    169169        #fill_color=rgb(0.85,0.85,0.85)
    170    
     170
    171171    lw=0.3
    172172    m.drawmapboundary()
     
    179179    xleft=gdict['longitudeOfFirstGridPointInDegrees']
    180180    if xleft>180.0:
    181         xleft-=360. 
     181        xleft-=360.
    182182    x=linspace(xleft,gdict['longitudeOfLastGridPointInDegrees'],gdict['Ni'])
    183183    y=linspace(gdict['latitudeOfLastGridPointInDegrees'],gdict['latitudeOfFirstGridPointInDegrees'],gdict['Nj'])
     
    186186    s=m.contourf(xx,yy, flist)
    187187    title(ftitle,y=1.1)
    188     cbaxes = f.add_axes([0.9, 0.2, 0.04, 0.6]) 
     188    cbaxes = f.add_axes([0.9, 0.2, 0.04, 0.6])
    189189    cb=colorbar(cax=cbaxes)
    190    
     190
    191191    savefig(c.outputdir+'/'+filename)
    192192    print 'created ',c.outputdir+'/'+filename
     
    200200                            formatter_class=ArgumentDefaultsHelpFormatter)
    201201
    202 # the most important arguments 
     202# the most important arguments
    203203    parser.add_argument("--start_date", dest="start_date",
    204204                        help="start date YYYYMMDD")
     
    241241            print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information'
    242242            exit(1)
    243    
    244     if args.levelist:       
     243
     244    if args.levelist:
    245245        c.levels=args.levelist.split('/')
    246246    else:
     
    250250    else:
    251251        c.area='[0,0]'
    252    
     252
    253253    c.paramIds=args.paramIds.split('/')
    254254    if args.start_step:
     
    269269    else:
    270270        c.outputdir=c.inputdir
    271    
     271
    272272        return args,c
    273273
    274274if __name__ == "__main__":
    275    
     275
    276276    args,c=interpret_plotargs()
    277277    plot_retrieved(args,c)
  • python/prepareFLEXPART.py

    rd69b677 r64cf353  
    11#!/usr/bin/env python
    2 #             
     2#
    33# This software is licensed under the terms of the Apache Licence Version 2.0
    4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 
    5 # 
     4# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
     5#
    66# Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs
    77#
    88# Creation: October  2014 - Anne Fouilloux - University of Oslo
    9 # Extension November 2015 - Leopold Haimberger - University of Vienna for: 
     9# Extension November 2015 - Leopold Haimberger - University of Vienna for:
    1010# - using the WebAPI also for general MARS retrievals
    1111# - job submission on ecgate and cca
     
    1515# - conversion into GRIB2
    1616# - conversion into .fp format for faster execution of FLEXPART
    17 # 
    18 # Requirements: 
     17#
     18# Requirements:
    1919# in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed
    20 # ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/
     20# ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit,
     21# all available from https://software.ecmwf.int/
    2122# dateutils
    2223# matplotlib (optional, for debugging)
    23 # 
     24#
    2425import calendar
    2526import shutil
     
    3637#from string import strip
    3738from GribTools import GribTools
    38 from FlexpartTools import EIFlexpart, Control,interpret_args_and_control, cleanup
     39from FlexpartTools import EIFlexpart, Control, interpret_args_and_control, cleanup
    3940
    4041hostname=socket.gethostname()
     
    4546        import ecmwfapi
    4647except ImportError:
    47     ecapi=False
    48    
     48    ecapi = False
     49
    4950
    5051def prepareFLEXPART(args,c):
     
    5354
    5455    namelist='fort.4'
    55    
     56
    5657    if not args.ppid:
    5758        c.ppid=str(os.getppid())
     
    7879    inputfiles.listFiles(c.inputdir, '*OG_acc_SL*.'+c.ppid+'.*')
    7980    if not os.path.exists(c.outputdir):
    80         os.makedirs(c.outputdir)
    81        
     81    os.makedirs(c.outputdir)
     82
    8283    flexpart = EIFlexpart(c,fluxes=True)
    8384    flexpart.create_namelist(c,'fort.4')
     
    9798
    9899    if int(c.debug)!=0:
    99         print('Temporary files left intact')
     100    print('Temporary files left intact')
    100101    else:
    101         cleanup(c)
     102    cleanup(c)
    102103
    103104if __name__ == "__main__":
    104     args,c=interpret_args_and_control()
    105     prepareFLEXPART(args,c)
     105    args, c = interpret_args_and_control()
     106    prepareFLEXPART(args, c)
    106107    cleanup(c)
    107    
     108
  • python/submit.py

    rd69b677 r64cf353  
    11#!/usr/bin/env python
    2 #             
    3 # This software is licensed under the terms of the Apache Licence Version 2.0
    4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
    5 #
    6 # Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs
    7 #
    8 # Creation: October  2014 - Anne Fouilloux - University of Oslo
    9 # Extension November 2015 - Leopold Haimberger - University of Vienna for:
    10 # - using the WebAPI also for general MARS retrievals
    11 # - job submission on ecgate and cca
    12 # - job templates suitable for twice daily operational dissemination
    13 # - dividing retrievals of longer periods into digestable chunks
    14 # - retrieve also longer term forecasts, not only analyses and short term forecast data
    15 # - conversion into GRIB2
    16 # - conversion into .fp format for faster execution of FLEXPART
    17 #
    18 # Further documentation may be obtained from www.flexpart.eu
    19 #
    20 # Requirements:
    21 # in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed
    22 # ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/
    23 # dateutils
    24 # matplotlib (optional, for debugging)
    25 #
    26 #
     2# -*- coding: utf-8 -*-
     3#************************************************************************
     4# TODO AP
     5#AP
     6# - Change History ist nicht angepasst ans File! Original geben lassen
     7# - dead code ? what to do?
     8# - seperate operational and reanlysis for clarification
     9#************************************************************************
     10"""
     11@Author: Anne Fouilloux (University of Oslo)
     12
     13@Date: October 2014
     14
     15@ChangeHistory:
     16    November 2015 - Leopold Haimberger (University of Vienna):
     17        - using the WebAPI also for general MARS retrievals
     18        - job submission on ecgate and cca
     19        - job templates suitable for twice daily operational dissemination
     20        - dividing retrievals of longer periods into digestable chunks
     21        - retrieve also longer term forecasts, not only analyses and
     22          short term forecast data
     23        - conversion into GRIB2
     24        - conversion into .fp format for faster execution of FLEXPART
     25
     26    February 2018 - Anne Philipp (University of Vienna):
     27        - applied PEP8 style guide
     28        - added documentation
     29        - minor changes in programming style for consistence
     30
     31@License:
     32    (C) Copyright 2014 UIO.
     33
     34    This software is licensed under the terms of the Apache Licence Version 2.0
     35    which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
     36
     37@Requirements:
     38    - A standard python 2.6 or 2.7 installation
     39    - dateutils
     40    - matplotlib (optional, for debugging)
     41    - ECMWF specific packages, all available from https://software.ecmwf.int/
     42        ECMWF WebMARS, gribAPI with python enabled, emoslib and
     43        ecaccess web toolkit
     44
     45@Description:
     46    Further documentation may be obtained from www.flexpart.eu.
     47
     48    Functionality provided:
     49        Prepare input 3D-wind fields in hybrid coordinates +
     50        surface fields for FLEXPART runs
     51"""
     52# ------------------------------------------------------------------------------
     53# MODULES
     54# ------------------------------------------------------------------------------
    2755import calendar
    2856import shutil
    2957import datetime
    3058import time
    31 import os,sys,glob
     59import os, sys, glob
    3260import subprocess
    3361import inspect
    3462# add path to submit.py to pythonpath so that python finds its buddies
    35 localpythonpath=os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
     63localpythonpath = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
    3664sys.path.append(localpythonpath)
    3765from UIOTools import UIOFiles
    3866from string import strip
    39 from argparse import ArgumentParser,ArgumentDefaultsHelpFormatter
     67from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
    4068from GribTools import GribTools
    41 from FlexpartTools import EIFlexpart, Control, interpret_args_and_control, normalexit,myerror
     69from FlexpartTools import EIFlexpart, Control, interpret_args_and_control, normalexit, myerror
    4270from getMARSdata import getMARSdata
    4371from prepareFLEXPART import prepareFLEXPART
     72# ------------------------------------------------------------------------------
     73# FUNCTIONS
     74# ------------------------------------------------------------------------------
     75def main():
     76    '''
     77    @Description:
     78        Get the arguments from script call and initialize an object from
     79        Control class. Decides from the argument "queue" if the local version
     80        is done "queue=None" or the gateway version "queue=ecgate".
    4481
     82    @Input:
     83        <nothing>
    4584
    46 
    47 def main():
    48 
    49     calledfromdir=os.getcwd()
    50 #    os.chdir(localpythonpath)
    51     args,c=interpret_args_and_control()
    52     if args.queue==None:
    53         if c.inputdir[0]!='/':
    54             c.inputdir=os.path.join(calledfromdir,c.inputdir)
    55         if c.outputdir[0]!='/':
    56             c.outputdir=os.path.join(calledfromdir,c.outputdir)
    57         getMARSdata(args,c)
    58         prepareFLEXPART(args,c)
     85    @Return:
     86        <nothing>
     87    '''
     88    calledfromdir = os.getcwd()
     89    args, c = interpret_args_and_control()
     90    if args.queue is None:
     91        if c.inputdir[0] != '/':
     92            c.inputdir = os.path.join(calledfromdir, c.inputdir)
     93        if c.outputdir[0] != '/':
     94            c.outputdir = os.path.join(calledfromdir, c.outputdir)
     95        getMARSdata(args, c)
     96        prepareFLEXPART(args, c)
    5997        normalexit(c)
    6098    else:
    61         submit(args.job_template,c,args.queue)
    62        
    63        
    64 def submit(jtemplate,c,queue):
    65    
    66     f=open(jtemplate)
    67     lftext=f.read().split('\n')
    68     insert_point=lftext.index('EOF')
    69     f.close()
    70    
    71     clist=c.tolist()
    72     colist=[]
    73     mt=0
     99        submit(args.job_template, c, args.queue)
     100
     101    return
     102
     103def submit(jtemplate, c, queue):
     104    #AP divide in two submits , ondemand und operational
     105    '''
     106    @Description:
     107        Prepares the job script and submit it to the specified queue.
     108
     109    @Input:
     110        jtemplate: string
     111            Job template file for submission to ECMWF. It contains all necessary
     112            module and variable settings for the ECMWF environment as well as
     113            the job call and mail report instructions.
     114            Default is "job.temp".
     115
     116        c: instance of class Control
     117            Contains all necessary information of a control file. The parameters
     118            are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     119            NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST,
     120            RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA,
     121            SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
     122            MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR
     123            For more information about format and content of the parameter see
     124            documentation.
     125
     126        queue: string
     127            Name of queue for submission to ECMWF (e.g. ecgate or cca )
     128
     129    @Return:
     130        <nothing>
     131    '''
     132
     133    # read template file and split from newline signs
     134    with open(jtemplate) as f:
     135        lftext = f.read().split('\n')
     136        insert_point = lftext.index('EOF')
     137#AP es gibt mehrere EOFs überprüfen!
     138
     139    # put all parameters of control instance into a list
     140    clist = c.tolist()  # reanalysis (EI)
     141    colist = []  # operational
     142    mt = 0
     143
     144#AP wieso 2 for loops?
     145#AP dieser part ist für die CTBTO Operational retrieves bis zum aktuellen Tag.
    74146    for elem in clist:
    75147        if 'maxstep' in elem:
    76             mt=int(elem.split(' ')[1])
    77            
     148            mt = int(elem.split(' ')[1])
     149
    78150    for elem in clist:
    79151        if 'start_date' in elem:
    80             elem='start_date '+'${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}'
     152            elem = 'start_date ' + '${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}'
    81153        if 'end_date' in elem:
    82             elem='start_date '+'${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}'
     154#AP Fehler?! Muss end_date heissen
     155            elem = 'start_date ' + '${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}'
    83156        if 'base_time' in elem:
    84             elem='base_time '+'${MSJ_BASETIME}'
    85         if 'time' in elem and mt>24:
    86             elem='time '+'${MSJ_BASETIME} {MSJ_BASETIME}'
     157            elem = 'base_time ' + '${MSJ_BASETIME}'
     158        if 'time' in elem and mt > 24:
     159            elem = 'time ' + '${MSJ_BASETIME} {MSJ_BASETIME}'
    87160        colist.append(elem)
    88        
    89     lftextondemand=lftext[:insert_point]+clist+lftext[insert_point+2:]
    90     lftextoper=lftext[:insert_point]+colist+lftext[insert_point+2:]
    91    
    92     h=open('job.ksh','w')
    93     h.write('\n'.join(lftextondemand)) 
     161#AP end
     162
     163#AP whats the difference between clist and colist ?! What is MSJ?
     164
     165    lftextondemand = lftext[:insert_point] + clist + lftext[insert_point + 2:]
     166    lftextoper = lftext[:insert_point] + colist + lftext[insert_point + 2:]
     167
     168    with open('job.ksh', 'w') as h:
     169        h.write('\n'.join(lftextondemand))
     170
     171#AP this is not used ?! what is it for?
     172#maybe a differentiation is needed
     173    h = open('joboper.ksh', 'w')
     174    h.write('\n'.join(lftextoper))
    94175    h.close()
    95    
    96     h=open('joboper.ksh','w')
    97     h.write('\n'.join(lftextoper)) 
    98     h.close()
    99    
    100     try:   
    101         p=subprocess.check_call(['ecaccess-job-submit','-queueName',queue,'job.ksh'])
     176#AP end
     177
     178    # submit job script to queue
     179    try:
     180        p = subprocess.check_call(['ecaccess-job-submit', '-queueName',
     181                                   queue,' job.ksh'])
    102182    except:
    103         print 'ecaccess-job-submit failed, probably eccert has expired'
     183        print('ecaccess-job-submit failed, probably eccert has expired')
    104184        exit(1)
    105     #pout=p.communicate(input=s)[0]
    106     print 'You should get an email with subject flex.hostname.pid'
    107185
     186    print('You should get an email with subject flex.hostname.pid')
    108187
     188    return
    109189
    110            
    111190if __name__ == "__main__":
    112191    main()
  • python/testsuite.py

    rd69b677 r64cf353  
    22
    33# This software is licensed under the terms of the Apache Licence Version 2.0
    4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 
    5 # 
     4# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
     5#
    66# Leopold Haimberger, Dec 2015
    77#
    88# Functionality provided: This script triggers the ECMWFDATA test suite. Call with
    99# testsuite.py [test group]
    10 # 
     10#
    1111#
    1212# Further documentation may be obtained from www.flexpart.eu
    13 # 
     13#
    1414# Test groups are specified in testsuite.json
    1515# in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed
     
    2121import json
    2222import subprocess
     23
    2324try:
    2425    taskfile=open('testsuite.json')
     
    4142    os.makedirs('../test')
    4243if len(sys.argv)>1:
    43     groups=sys.argv[1:] 
     44    groups=sys.argv[1:]
    4445else:
    4546    groups=['xinstall','default','ops','work','cv','fc']#,'hires']
     
    5354    garglist=[]
    5455    for ttk,ttv in tv.iteritems():
    55         if isinstance(ttv,basestring):   
     56        if isinstance(ttv,basestring):
    5657            if ttk!='script':
    5758                garglist.append('--'+ttk)
     
    6465            arglist=[]
    6566            for tttk,tttv in ttv.iteritems():
    66                 if isinstance(tttv,basestring):   
     67                if isinstance(tttv,basestring):
    6768                        arglist.append('--'+tttk)
    6869                        if '$' in tttv[0]:
    69                             arglist.append(os.path.expandvars(tttv))                       
     70                            arglist.append(os.path.expandvars(tttv))
    7071                        else:
    7172                            arglist.append(tttv)
     
    8687print str(jobcounter-jobfailed)+' successful, '+str(jobfailed)+' failed'
    8788print 'If tasks have been submitted via ECACCESS please check emails'
    88                                        
     89
    8990#print tasks
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG