source: flex_extract.git/python/FlexpartTools.py @ 64cf353

ctbtodev
Last change on this file since 64cf353 was 64cf353, checked in by Anne Philipp <bscannephilipp@…>, 6 years ago

pep8 changes + documentation added + minor code style changes
  • Property mode set to 100644
File size: 77.2 KB
Line 
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# ------------------------------------------------------------------------------
59import subprocess
60from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
61import traceback
62import shutil
63import os
64import errno
65import sys
66import inspect
67import glob
68import datetime
69from string import atoi
70from numpy import *
71ecapi = True
72try:
73    import ecmwfapi
74except ImportError:
75    ecapi = False
76from gribapi import *
77from GribTools import GribTools
78
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',
105                            formatter_class=ArgumentDefaultsHelpFormatter)
106
107# the most important arguments
108    parser.add_argument("--start_date", dest="start_date",
109                        help="start date YYYYMMDD")
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
115# some arguments that override the default in the control file
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
125# set the working directories
126    parser.add_argument("--inputdir", dest="inputdir", default=None,
127                        help="root directory for storing intermediate files")
128    parser.add_argument("--outputdir", dest="outputdir", default=None,
129                        help="root directory for storing output files")
130    parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
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
135# this is only used by prepareFLEXPART.py to rerun a postprocessing step
136    parser.add_argument("--ppid", dest="ppid",
137                        help="Specify parent process id for \
138                        rerun of prepareFLEXPART")
139
140# arguments for job submission to ECMWF, only needed by submit.py
141    parser.add_argument("--job_template", dest='job_template',
142                        default="job.temp",
143                        help="job template file for submission to ECMWF")
144    #parser.add_argument("--remote", dest="remote",
145                        #help="target for submission to ECMWF \
146                        #(ecgate or cca etc.)")
147    parser.add_argument("--queue", dest="queue",
148                        help="queue for submission to ECMWF \
149                        (e.g. ecgate or cca )")
150    parser.add_argument("--controlfile", dest="controlfile",
151                        default='CONTROL.temp',
152                        help="file with control parameters")
153    parser.add_argument("--debug", dest="debug", default=0,
154                        help="Debug mode - leave temporary files intact")
155
156    args = parser.parse_args()
157
158    try:
159        c = Control(args.controlfile)
160    except IOError:
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'
192    else:
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
204    else:
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
244
245def install_args_and_control():
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',
270                            formatter_class=ArgumentDefaultsHelpFormatter)
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')
285
286    parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
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
292# arguments for job submission to ECMWF, only needed by submit.py
293    parser.add_argument("--job_template", dest='job_template',
294                        default="job.temp.o",
295                        help="job template file for submission to ECMWF")
296
297    parser.add_argument("--controlfile", dest="controlfile",
298                        default='CONTROL.temp',
299                        help="file with control parameters")
300
301    args = parser.parse_args()
302
303    try:
304        c = Control(args.controlfile)
305    except:
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
331    try:
332        c.makefile = args.makefile
333    except:
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
350
351def cleanup(c):
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 + "/*")
376    for cl in cleanlist:
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
421    # uncomment if user wants email notification directly from python
422    #for t in target:
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
430    exit(1)
431
432    return
433
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
480def product(*args, **kwds):
481    '''
482    @Description:
483
484    @Input:
485
486    @Return:
487
488    '''
489    # product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
490    # product(range(2), repeat = 3) --> 000 001 010 011 100 101 110 111
491    pools = map(tuple, args) * kwds.get('repeat', 1)
492    result = [[]]
493    for pool in pools:
494        result = [x + [y] for x in result for y in pool]
495    for prod in result:
496        yield tuple(prod)
497
498    return
499
500
501def 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    '''
515    try:
516        os.remove(filename)
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
524
525
526def init128(fn):
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()
540    with open(fn) as f:
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
546    return table128
547
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 = []
570    for par in cpar:
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
580    return ipar
581
582
583def dapoly(alist):
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
598    return sfeld
599
600
601def darain(alist):
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
634    return sfeld
635
636
637class Control:
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
770    def __str__(self):
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
789    def tolist(self):
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
830class MARSretrieval:
831    'class for MARS retrievals'
832
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
892
893    def displayInfo(self):
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
912
913    def dataRetrieve(self):
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
948        else:
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##############################################################
967##############################################################
968class EIFlexpart:
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")
1290        #QG.displayInfo()
1291        #QG.dataRetrieve()
1292
1293        oro = False
1294        for ftype in self.types:
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))
1459        ctype = str(type).upper()
1460
1461        myinfo = [ctype, ctime, cstep]
1462        cflextime = None
1463        for t, marsinfo in self.mars.items():
1464            if myinfo == marsinfo:
1465                cflextime = t
1466
1467        return cflextime
1468
1469    def process_output(self, c):
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"
1622        silentremove(indexfile)
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
1859        index_vals = []
1860        for key in index_keys:
1861            key_vals = grib_index_get(iid,key)
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))
1873            index_vals.append(key_vals)
1874
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
1934                cparamId = str(grib_get(gid, 'paramId'))
1935                step = grib_get(gid, 'step')
1936                atime = grib_get(gid, 'time')
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
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG