source: flex_extract.git/python/EcFlexpart.py @ e1228f3

devfeature/makefilestask/language-editing
Last change on this file since e1228f3 was e1228f3, checked in by Anne Philipp <anne.philipp@…>, 2 years ago

resolved loop import between controlfile and tools

  • Property mode set to 100644
File size: 56.3 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3#************************************************************************
4# ToDo AP
5# - specifiy file header documentation
6# - add class description in header information
7# - apply classtests
8# - add references to ECMWF specific software packages
9# - add describtion of deacc_fluxes
10# - change name of func deacc ( weil disagg )
11# - add desc of retrieve function
12#************************************************************************
13#*******************************************************************************
14# @Author: Anne Fouilloux (University of Oslo)
15#
16# @Date: October 2014
17#
18# @Change History:
19#
20#    November 2015 - Leopold Haimberger (University of Vienna):
21#        - extended with class Control
22#        - removed functions mkdir_p, daterange, years_between, months_between
23#        - added functions darain, dapoly, to_param_id, init128, normal_exit,
24#          my_error, clean_up, install_args_and_control,
25#          interpret_args_and_control,
26#        - removed function __del__ in class EIFLexpart
27#        - added the following functions in EIFlexpart:
28#            - create_namelist
29#            - process_output
30#            - deacc_fluxes
31#        - modified existing EIFlexpart - functions for the use in
32#          flex_extract
33#        - retrieve also longer term forecasts, not only analyses and
34#          short term forecast data
35#        - added conversion into GRIB2
36#        - added conversion into .fp format for faster execution of FLEXPART
37#          (see https://www.flexpart.eu/wiki/FpCtbtoWo4FpFormat)
38#
39#    February 2018 - Anne Philipp (University of Vienna):
40#        - applied PEP8 style guide
41#        - added documentation
42#        - removed function getFlexpartTime in class EcFlexpart
43#        - outsourced class ControlFile
44#        - outsourced class MarsRetrieval
45#        - changed class name from EIFlexpart to EcFlexpart
46#        - applied minor code changes (style)
47#        - removed "dead code" , e.g. retrieval of Q since it is not needed
48#        - removed "times" parameter from retrieve-method since it is not used
49#
50# @License:
51#    (C) Copyright 2014-2018.
52#
53#    This software is licensed under the terms of the Apache Licence Version 2.0
54#    which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
55#
56# @Class Description:
57#    FLEXPART needs grib files in a specifc format. All necessary data fields
58#    for one time step are stored in a single file. The class represents an
59#    instance with all the parameter and settings necessary for retrieving
60#    MARS data and modifing them so they are fitting FLEXPART need. The class
61#    is able to disaggregate the fluxes and convert grid types to the one needed
62#    by FLEXPART, therefore using the FORTRAN program.
63#
64# @Class Content:
65#    - __init__
66#    - write_namelist
67#    - retrieve
68#    - process_output
69#    - create
70#    - deacc_fluxes
71#
72# @Class Attributes:
73#    - dtime
74#    - basetime
75#    - server
76#    - marsclass
77#    - stream
78#    - resol
79#    - accuracy
80#    - number
81#    - expver
82#    - glevelist
83#    - area
84#    - grid
85#    - level
86#    - levelist
87#    - types
88#    - dates
89#    - area
90#    - gaussian
91#    - params
92#    - inputdir
93#    - outputfilelist
94#
95#*******************************************************************************
96#pylint: disable=unsupported-assignment-operation
97# this is disabled because its an error in pylint for this specific case
98#pylint: disable=consider-using-enumerate
99# this is not useful in this case
100# ------------------------------------------------------------------------------
101# MODULES
102# ------------------------------------------------------------------------------
103import subprocess
104import shutil
105import os
106import glob
107from datetime import datetime, timedelta
108import numpy as np
109from gribapi import grib_set, grib_index_select, grib_new_from_index, grib_get,\
110                    grib_write, grib_get_values, grib_set_values, grib_release,\
111                    grib_index_release, grib_index_get
112
113# software specific classes and modules from flex_extract
114from GribTools import GribTools
115from tools import init128, to_param_id, silent_remove, product, my_error
116from MarsRetrieval import MarsRetrieval
117import disaggregation
118
119# ------------------------------------------------------------------------------
120# CLASS
121# ------------------------------------------------------------------------------
122class EcFlexpart(object):
123    '''
124    Class to retrieve FLEXPART specific ECMWF data.
125    '''
126    # --------------------------------------------------------------------------
127    # CLASS FUNCTIONS
128    # --------------------------------------------------------------------------
129    def __init__(self, c, fluxes=False):
130        '''
131        @Description:
132            Creates an object/instance of EcFlexpart with the
133            associated settings of its attributes for the retrieval.
134
135        @Input:
136            self: instance of EcFlexpart
137                The current object of the class.
138
139            c: instance of class ControlFile
140                Contains all the parameters of CONTROL file, which are e.g.:
141                DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
142                STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
143                LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
144                OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
145                ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
146                MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
147                DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
148
149                For more information about format and content of the parameter
150                see documentation.
151
152            fluxes: boolean, optional
153                Decides if a the flux parameter settings are stored or
154                the rest of the parameter list.
155                Default value is False.
156
157        @Return:
158            <nothing>
159        '''
160
161        # different mars types for retrieving data for flexpart
162        self.types = dict()
163
164        if c.maxstep > len(c.type):    # Pure forecast mode
165            c.type = [c.type[1]]
166            c.step = ['{:0>3}'.format(int(c.step[0]))]
167            c.time = [c.time[0]]
168            for i in range(1, c.maxstep + 1):
169                c.type.append(c.type[0])
170                c.step.append('{:0>3}'.format(i))
171                c.time.append(c.time[0])
172
173        self.inputdir = c.inputdir
174        self.basetime = c.basetime
175        self.dtime = c.dtime
176        i = 0
177        if fluxes is True and c.maxstep < 24:
178            # no forecast beyond one day is needed!
179            # Thus, prepare flux data manually as usual
180            # with only forecast fields with start times at 00/12
181            # (but without 00/12 fields since these are
182            # the initialisation times of the flux fields
183            # and therefore are zero all the time)
184            self.types[c.type[1]] = {'times': '00/12', 'steps':
185                                     '{}/to/12/by/{}'.format(c.dtime, c.dtime)}
186        else:
187            for ty, st, ti in zip(c.type, c.step, c.time):
188                btlist = range(24)
189                if c.basetime == '12':
190                    btlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
191                if c.basetime == '00':
192                    btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0]
193
194                if i % int(c.dtime) == 0 and c.maxstep > 24 or i in btlist:
195
196                    if ty not in self.types.keys():
197                        self.types[ty] = {'times': '', 'steps': ''}
198
199                    if ti not in self.types[ty]['times']:
200                        if self.types[ty]['times']:
201                            self.types[ty]['times'] += '/'
202                        self.types[ty]['times'] += ti
203
204                    if st not in self.types[ty]['steps']:
205                        if self.types[ty]['steps']:
206                            self.types[ty]['steps'] += '/'
207                        self.types[ty]['steps'] += st
208                i += 1
209
210
211        self.marsclass = c.marsclass
212        self.stream = c.stream
213        self.number = c.number
214        self.resol = c.resol
215        self.accuracy = c.accuracy
216        self.level = c.level
217
218        if c.levelist:
219            self.levelist = c.levelist
220        else:
221            self.levelist = '1/to/' + c.level
222
223        # for gaussian grid retrieval
224        self.glevelist = '1/to/' + c.level
225
226        if hasattr(c, 'gaussian') and c.gaussian:
227            self.gaussian = c.gaussian
228        else:
229            self.gaussian = ''
230
231        if hasattr(c, 'expver') and c.expver:
232            self.expver = c.expver
233        else:
234            self.expver = '1'
235
236        if hasattr(c, 'number') and c.number:
237            self.number = c.number
238        else:
239            self.number = '0'
240
241        if 'N' in c.grid:  # Gaussian output grid
242            self.grid = c.grid
243            self.area = 'G'
244        else:
245            self.grid = '{}/{}'.format(int(c.grid) / 1000., int(c.grid) / 1000.)
246            self.area = '{}/{}/{}/{}'.format(int(c.upper) / 1000.,
247                                             int(c.left) / 1000.,
248                                             int(c.lower) / 1000.,
249                                             int(c.right) / 1000.)
250
251        self.outputfilelist = []
252
253
254        # Now comes the nasty part that deals with the different
255        # scenarios we have:
256        # 1) Calculation of etadot on
257        #    a) Gaussian grid
258        #    b) Output grid
259        #    c) Output grid using parameter 77 retrieved from MARS
260        # 3) Calculation/Retrieval of omega
261        # 4) Download also data for WRF
262
263
264        # Different grids need different retrievals
265        # SH = Spherical Harmonics, GG = Gaussian Grid,
266        # OG = Output Grid, ML = MultiLevel, SL = SingleLevel
267        self.params = {'SH__ML': '', 'SH__SL': '',
268                       'GG__ML': '', 'GG__SL': '',
269                       'OG__ML': '', 'OG__SL': '',
270                       'OG_OROLSM_SL': '', 'OG_acc_SL': ''}
271
272        if fluxes is False:
273            self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF']
274            #                        "SD/MSL/TCC/10U/10V/2T/2D/129/172"
275            self.params['OG__SL'] = ["141/151/164/165/166/167/168/129/172", \
276                                     'SFC', '1', self.grid]
277            if c.addpar:
278                if c.addpar[0] == '/':
279                    c.addpar = c.addpar[1:]
280                self.params['OG__SL'][0] += '/' + '/'.join(c.addpar)
281
282            self.params['OG_OROLSM__SL'] = ["160/27/28/173", \
283                                            'SFC', '1', self.grid]
284
285            self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]
286
287            if c.gauss == '0' and c.eta == '1':
288                # the simplest case
289                self.params['OG__ML'][0] += '/U/V/77'
290            elif c.gauss == '0' and c.eta == '0':
291            # this is not recommended (inaccurate)
292                self.params['OG__ML'][0] += '/U/V'
293            elif c.gauss == '1' and c.eta == '0':
294                # this is needed for data before 2008, or for reanalysis data
295                self.params['GG__SL'] = ['Q', 'ML', '1', \
296                                         '{}'.format((int(self.resol) + 1) / 2)]
297                self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF']
298            else:
299                print 'Warning: This is a very costly parameter combination, \
300                       use only for debugging!'
301                self.params['GG__SL'] = ['Q', 'ML', '1', \
302                                         '{}'.format((int(self.resol) + 1) / 2)]
303                self.params['GG__ML'] = ['U/V/D/77', 'ML', self.glevelist, \
304                                         '{}'.format((int(self.resol) + 1) / 2)]
305
306            if c.omega == '1':
307                self.params['OG__ML'][0] += '/W'
308
309            # add cloud water content if necessary
310            if c.cwc == '1':
311                self.params['OG__ML'][0] += '/CLWC/CIWC'
312
313            # add vorticity and geopotential height for WRF if necessary
314            if c.wrf == '1':
315                self.params['OG__ML'][0] += '/Z/VO'
316                if '/D' not in self.params['OG__ML'][0]:
317                    self.params['OG__ML'][0] += '/D'
318                #wrf_sfc = 'sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/ /
319                #           stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper()
320                wrf_sfc = '134/235/167/165/166/168/129/172/34/31/141/ \
321                           139/170/183/236/39/40/41/42'.upper()
322                lwrt_sfc = wrf_sfc.split('/')
323                for par in lwrt_sfc:
324                    if par not in self.params['OG__SL'][0]:
325                        self.params['OG__SL'][0] += '/' + par
326
327        else:
328            self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR", \
329                                        'SFC', '1', self.grid]
330
331        # if needed, add additional WRF specific parameters here
332
333        return
334
335
336    def write_namelist(self, c, filename):
337        '''
338        @Description:
339            Creates a namelist file in the temporary directory and writes
340            the following values to it: maxl, maxb, mlevel,
341            mlevelist, mnauf, metapar, rlo0, rlo1, rla0, rla1,
342            momega, momegadiff, mgauss, msmooth, meta, metadiff, mdpdeta
343
344        @Input:
345            self: instance of EcFlexpart
346                The current object of the class.
347
348            c: instance of class ControlFile
349                Contains all the parameters of CONTROL files, which are e.g.:
350                DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
351                STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
352                LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
353                OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
354                ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
355                MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
356                DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
357
358                For more information about format and content of the parameter
359                see documentation.
360
361            filename: string
362                Name of the namelist file.
363
364        @Return:
365            <nothing>
366        '''
367
368        self.inputdir = c.inputdir
369        area = np.asarray(self.area.split('/')).astype(float)
370        grid = np.asarray(self.grid.split('/')).astype(float)
371
372        if area[1] > area[3]:
373            area[1] -= 360
374        maxl = int((area[3] - area[1]) / grid[1]) + 1
375        maxb = int((area[0] - area[2]) / grid[0]) + 1
376
377        with open(self.inputdir + '/' + filename, 'w') as f:
378            f.write('&NAMGEN\n')
379            f.write(',\n  '.join(['maxl = ' + str(maxl), 'maxb = ' + str(maxb),
380                                  'mlevel = ' + self.level,
381                                  'mlevelist = ' + '"' + self.levelist + '"',
382                                  'mnauf = ' + self.resol,
383                                  'metapar = ' + '77',
384                                  'rlo0 = ' + str(area[1]),
385                                  'rlo1 = ' + str(area[3]),
386                                  'rla0 = ' + str(area[2]),
387                                  'rla1 = ' + str(area[0]),
388                                  'momega = ' + c.omega,
389                                  'momegadiff = ' + c.omegadiff,
390                                  'mgauss = ' + c.gauss,
391                                  'msmooth = ' + c.smooth,
392                                  'meta = ' + c.eta,
393                                  'metadiff = ' + c.etadiff,
394                                  'mdpdeta = ' + c.dpdeta]))
395
396            f.write('\n/\n')
397
398        return
399
400    def retrieve(self, server, dates, inputdir='.'):
401        '''
402        @Description:
403            Finalizing the retrieval information by setting final details
404            depending on grid type.
405            Prepares MARS retrievals per grid type and submits them.
406
407        @Input:
408            self: instance of EcFlexpart
409                The current object of the class.
410
411            server: instance of ECMWFService or ECMWFDataServer
412                The connection to the ECMWF server. This is different
413                for member state users which have full access and non
414                member state users which have only access to the public
415                data sets. The decision is made from command line argument
416                "public"; for public access its True (ECMWFDataServer)
417                for member state users its False (ECMWFService)
418
419            dates: string
420                Contains start and end date of the retrieval in the format
421                "YYYYMMDD/to/YYYYMMDD"
422
423            inputdir: string, optional
424                Path to the directory where the retrieved data is about
425                to be stored. The default is the current directory ('.').
426
427        @Return:
428            <nothing>
429        '''
430        self.dates = dates
431        self.server = server
432        self.inputdir = inputdir
433        oro = False
434
435        for ftype in self.types:
436            for pk, pv in self.params.iteritems():
437                if isinstance(pv, str):
438                    continue
439                mftype = '' + ftype
440                mftime = self.types[ftype]['times']
441                mfstep = self.types[ftype]['steps']
442                mfdate = self.dates
443                mfstream = self.stream
444                mftarget = self.inputdir + "/" + ftype + pk + '.' + \
445                           self.dates.split('/')[0] + '.' + str(os.getppid()) +\
446                           '.' + str(os.getpid()) + ".grb"
447                if pk == 'OG__SL':
448                    pass
449                if pk == 'OG_OROLSM__SL':
450                    if oro is False:
451                        mfstream = 'OPER'
452                        mftype = 'AN'
453                        mftime = '00'
454                        mfstep = '000'
455                        mfdate = self.dates.split('/')[0]
456                        mftarget = self.inputdir + "/" + pk + '.' + mfdate + \
457                                   '.' + str(os.getppid()) + '.' + \
458                                   str(os.getpid()) + ".grb"
459                        oro = True
460                    else:
461                        continue
462                if pk == 'GG__SL' and pv[0] == 'Q':
463                    area = ""
464                    gaussian = 'reduced'
465                else:
466                    area = self.area
467                    gaussian = self.gaussian
468
469    # ------  on demand path  --------------------------------------------------
470                if self.basetime is None:
471                    MR = MarsRetrieval(self.server,
472                                       marsclass=self.marsclass,
473                                       stream=mfstream,
474                                       type=mftype,
475                                       levtype=pv[1],
476                                       levelist=pv[2],
477                                       resol=self.resol,
478                                       gaussian=gaussian,
479                                       accuracy=self.accuracy,
480                                       grid=pv[3],
481                                       target=mftarget,
482                                       area=area,
483                                       date=mfdate,
484                                       time=mftime,
485                                       number=self.number,
486                                       step=mfstep,
487                                       expver=self.expver,
488                                       param=pv[0])
489
490                    MR.display_info()
491                    MR.data_retrieve()
492    # ------  operational path  ------------------------------------------------
493                else:
494                    # check if mars job requests fields beyond basetime.
495                    # If yes eliminate those fields since they may not
496                    # be accessible with user's credentials
497                    if 'by' in mfstep:
498                        sm1 = 2
499                    else:
500                        sm1 = -1
501
502                    if 'by' in mftime:
503                        tm1 = 2
504                    else:
505                        tm1 = -1
506
507                    maxdate = datetime.strptime(mfdate.split('/')[-1] +
508                                                mftime.split('/')[tm1],
509                                                '%Y%m%d%H')
510                    istep = int(mfstep.split('/')[sm1])
511                    maxtime = maxdate + timedelta(hours=istep)
512
513                    elimit = datetime.strptime(mfdate.split('/')[-1] +
514                                               self.basetime, '%Y%m%d%H')
515
516                    if self.basetime == '12':
517                        # --------------  flux data ----------------------------
518                        if 'acc' in pk:
519
520                        # Strategy:
521                        # if maxtime-elimit >= 24h reduce date by 1,
522                        # if 12h <= maxtime-elimit<12h reduce time for last date
523                        # if maxtime-elimit<12h reduce step for last time
524                        # A split of the MARS job into 2 is likely necessary.
525                            maxtime = elimit - timedelta(hours=24)
526                            mfdate = '/'.join(['/'.join(mfdate.split('/')[:-1]),
527                                               datetime.strftime(maxtime,
528                                                                 '%Y%m%d')])
529
530                            MR = MarsRetrieval(self.server,
531                                               marsclass=self.marsclass,
532                                               stream=self.stream,
533                                               type=mftype,
534                                               levtype=pv[1],
535                                               levelist=pv[2],
536                                               resol=self.resol,
537                                               gaussian=gaussian,
538                                               accuracy=self.accuracy,
539                                               grid=pv[3],
540                                               target=mftarget,
541                                               area=area,
542                                               date=mfdate,
543                                               time=mftime,
544                                               number=self.number,
545                                               step=mfstep,
546                                               expver=self.expver,
547                                               param=pv[0])
548
549                            MR.display_info()
550                            MR.data_retrieve()
551
552                            maxtime = elimit - timedelta(hours=12)
553                            mfdate = datetime.strftime(maxtime, '%Y%m%d')
554                            mftime = '00'
555                            mftarget = self.inputdir + "/" + ftype + pk + \
556                                       '.' + mfdate + '.' + str(os.getppid()) +\
557                                       '.' + str(os.getpid()) + ".grb"
558
559                            MR = MarsRetrieval(self.server,
560                                               marsclass=self.marsclass,
561                                               stream=self.stream,
562                                               type=mftype,
563                                               levtype=pv[1],
564                                               levelist=pv[2],
565                                               resol=self.resol,
566                                               gaussian=gaussian,
567                                               accuracy=self.accuracy,
568                                               grid=pv[3],
569                                               target=mftarget,
570                                               area=area,
571                                               date=mfdate,
572                                               time=mftime,
573                                               number=self.number,
574                                               step=mfstep,
575                                               expver=self.expver,
576                                               param=pv[0])
577
578                            MR.display_info()
579                            MR.data_retrieve()
580                        # --------------  non flux data ------------------------
581                        else:
582                            MR = MarsRetrieval(self.server,
583                                               marsclass=self.marsclass,
584                                               stream=self.stream,
585                                               type=mftype,
586                                               levtype=pv[1],
587                                               levelist=pv[2],
588                                               resol=self.resol,
589                                               gaussian=gaussian,
590                                               accuracy=self.accuracy,
591                                               grid=pv[3],
592                                               target=mftarget,
593                                               area=area,
594                                               date=mfdate,
595                                               time=mftime,
596                                               number=self.number,
597                                               step=mfstep,
598                                               expver=self.expver,
599                                               param=pv[0])
600
601                            MR.display_info()
602                            MR.data_retrieve()
603                    else: # basetime == 0 ??? #AP
604
605                        maxtime = elimit - timedelta(hours=24)
606                        mfdate = datetime.strftime(maxtime, '%Y%m%d')
607                        mftimesave = ''.join(mftime)
608
609                        if '/' in mftime:
610                            times = mftime.split('/')
611                            while ((int(times[0]) +
612                                    int(mfstep.split('/')[0]) <= 12) and
613                                   (pk != 'OG_OROLSM__SL') and 'acc' not in pk):
614                                times = times[1:]
615                            if len(times) > 1:
616                                mftime = '/'.join(times)
617                            else:
618                                mftime = times[0]
619
620                        MR = MarsRetrieval(self.server,
621                                           marsclass=self.marsclass,
622                                           stream=self.stream,
623                                           type=mftype,
624                                           levtype=pv[1],
625                                           levelist=pv[2],
626                                           resol=self.resol,
627                                           gaussian=gaussian,
628                                           accuracy=self.accuracy,
629                                           grid=pv[3],
630                                           target=mftarget,
631                                           area=area,
632                                           date=mfdate,
633                                           time=mftime,
634                                           number=self.number,
635                                           step=mfstep,
636                                           expver=self.expver,
637                                           param=pv[0])
638
639                        MR.display_info()
640                        MR.data_retrieve()
641
642                        if (int(mftimesave.split('/')[0]) == 0 and
643                                int(mfstep.split('/')[0]) == 0 and
644                                pk != 'OG_OROLSM__SL'):
645
646                            mfdate = datetime.strftime(elimit, '%Y%m%d')
647                            mftime = '00'
648                            mfstep = '000'
649                            mftarget = self.inputdir + "/" + ftype + pk + \
650                                       '.' + mfdate + '.' + str(os.getppid()) +\
651                                       '.' + str(os.getpid()) + ".grb"
652
653                            MR = MarsRetrieval(self.server,
654                                               marsclass=self.marsclass,
655                                               stream=self.stream,
656                                               type=mftype,
657                                               levtype=pv[1],
658                                               levelist=pv[2],
659                                               resol=self.resol,
660                                               gaussian=gaussian,
661                                               accuracy=self.accuracy,
662                                               grid=pv[3],
663                                               target=mftarget,
664                                               area=area,
665                                               date=mfdate,
666                                               time=mftime,
667                                               number=self.number,
668                                               step=mfstep,
669                                               expver=self.expver,
670                                               param=pv[0])
671
672                            MR.display_info()
673                            MR.data_retrieve()
674
675        print "MARS retrieve done... "
676
677        return
678
679
680    def process_output(self, c):
681        '''
682        @Description:
683            The grib files are postprocessed depending on the selection in
684            CONTROL file. The resulting files are moved to the output
685            directory if its not equla to the input directory.
686            The following modifications might be done if
687            properly switched in CONTROL file:
688            GRIB2 - Conversion to GRIB2
689            ECTRANS - Transfer of files to gateway server
690            ECSTORAGE - Storage at ECMWF server
691            GRIB2FLEXPART - Conversion of GRIB files to FLEXPART binary format
692
693        @Input:
694            self: instance of EcFlexpart
695                The current object of the class.
696
697            c: instance of class ControlFile
698                Contains all the parameters of CONTROL file, which are e.g.:
699                DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
700                STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
701                LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
702                OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
703                ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
704                MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
705                DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
706
707                For more information about format and content of the parameter
708                see documentation.
709
710        @Return:
711            <nothing>
712
713        '''
714
715        print '\n\nPostprocessing:\n Format: {}\n'.format(c.format)
716
717        if c.ecapi is False:
718            print('ecstorage: {}\n ecfsdir: {}\n'.
719                  format(c.ecstorage, c.ecfsdir))
720            if not hasattr(c, 'gateway'):
721                c.gateway = os.getenv('GATEWAY')
722            if not hasattr(c, 'destination'):
723                c.destination = os.getenv('DESTINATION')
724            print('ectrans: {}\n gateway: {}\n destination: {}\n '
725                  .format(c.ectrans, c.gateway, c.destination))
726
727        print 'Output filelist: \n'
728        print self.outputfilelist
729
730        if c.format.lower() == 'grib2':
731            for ofile in self.outputfilelist:
732                p = subprocess.check_call(['grib_set', '-s', 'edition=2, \
733                                           productDefinitionTemplateNumber=8',
734                                           ofile, ofile + '_2'])
735                p = subprocess.check_call(['mv', ofile + '_2', ofile])
736
737        if int(c.ectrans) == 1 and c.ecapi is False:
738            for ofile in self.outputfilelist:
739                p = subprocess.check_call(['ectrans', '-overwrite', '-gateway',
740                                           c.gateway, '-remote', c.destination,
741                                           '-source', ofile])
742                print('ectrans:', p)
743
744        if int(c.ecstorage) == 1 and c.ecapi is False:
745            for ofile in self.outputfilelist:
746                p = subprocess.check_call(['ecp', '-o', ofile,
747                                           os.path.expandvars(c.ecfsdir)])
748
749        if c.outputdir != c.inputdir:
750            for ofile in self.outputfilelist:
751                p = subprocess.check_call(['mv', ofile, c.outputdir])
752
753        # prepare environment for the grib2flexpart run
754        # to convert grib to flexpart binary
755        if c.grib2flexpart == '1':
756
757            # generate AVAILABLE file
758            # Example of AVAILABLE file data:
759            # 20131107 000000      EN13110700              ON DISC
760            clist = []
761            for ofile in self.outputfilelist:
762                fname = ofile.split('/')
763                if '.' in fname[-1]:
764                    l = fname[-1].split('.')
765                    timestamp = datetime.strptime(l[0][-6:] + l[1],
766                                                  '%y%m%d%H')
767                    timestamp += timedelta(hours=int(l[2]))
768                    cdate = datetime.strftime(timestamp, '%Y%m%d')
769                    chms = datetime.strftime(timestamp, '%H%M%S')
770                else:
771                    cdate = '20' + fname[-1][-8:-2]
772                    chms = fname[-1][-2:] + '0000'
773                clist.append(cdate + ' ' + chms + ' '*6 +
774                             fname[-1] + ' '*14 + 'ON DISC')
775            clist.sort()
776            with open(c.outputdir + '/' + 'AVAILABLE', 'w') as f:
777                f.write('\n'.join(clist) + '\n')
778
779            # generate pathnames file
780            pwd = os.path.abspath(c.outputdir)
781            with open(pwd + '/pathnames', 'w') as f:
782                f.write(pwd + '/Options/\n')
783                f.write(pwd + '/\n')
784                f.write(pwd + '/\n')
785                f.write(pwd + '/AVAILABLE\n')
786                f.write(' = == = == = == = == = == ==  = \n')
787
788            # create Options dir if necessary
789            if not os.path.exists(pwd + '/Options'):
790                os.makedirs(pwd+'/Options')
791
792            # read template COMMAND file
793            with open(os.path.expandvars(os.path.expanduser(
794                c.flexpart_root_scripts)) + '/../Options/COMMAND', 'r') as f:
795                lflist = f.read().split('\n')
796
797            # find index of list where to put in the
798            # date and time information
799            # usually after the LDIRECT parameter
800            i = 0
801            for l in lflist:
802                if 'LDIRECT' in l.upper():
803                    break
804                i += 1
805
806            # insert the date and time information of run start and end
807            # into the list of lines of COMMAND file
808            lflist = lflist[:i+1] + \
809                     [clist[0][:16], clist[-1][:16]] + \
810                     lflist[i+3:]
811
812            # write the new COMMAND file
813            with open(pwd + '/Options/COMMAND', 'w') as g:
814                g.write('\n'.join(lflist) + '\n')
815
816            # change to outputdir and start the grib2flexpart run
817            # afterwards switch back to the working dir
818            os.chdir(c.outputdir)
819            p = subprocess.check_call([
820                os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))
821                + '/../FLEXPART_PROGRAM/grib2flexpart', 'useAvailable', '.'])
822            os.chdir(pwd)
823
824        return
825
826    def create(self, inputfiles, c):
827        '''
828        @Description:
829            This method is based on the ECMWF example index.py
830            https://software.ecmwf.int/wiki/display/GRIB/index.py
831
832            An index file will be created which depends on the combination
833            of "date", "time" and "stepRange" values. This is used to iterate
834            over all messages in each grib file which were passed through the
835            parameter "inputfiles" to seperate specific parameters into fort.*
836            files. Afterwards the FORTRAN program Convert2 is called to convert
837            the data fields all to the same grid and put them in one file
838            per unique time step (combination of "date", "time" and
839            "stepRange").
840
841        @Input:
842            self: instance of EcFlexpart
843                The current object of the class.
844
845            inputfiles: instance of UioFiles
846                Contains a list of files.
847
848            c: instance of class ControlFile
849                Contains all the parameters of CONTROL files, which are e.g.:
850                DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
851                STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
852                LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
853                OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
854                ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
855                MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
856                DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
857
858                For more information about format and content of the parameter
859                see documentation.
860
861        @Return:
862            <nothing>
863        '''
864
865        table128 = init128(c.ecmwfdatadir +
866                           '/grib_templates/ecmwf_grib1_table_128')
867        wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
868                            stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
869                            table128)
870
871        index_keys = ["date", "time", "step"]
872        indexfile = c.inputdir + "/date_time_stepRange.idx"
873        silent_remove(indexfile)
874        grib = GribTools(inputfiles.files)
875        # creates new index file
876        iid = grib.index(index_keys=index_keys, index_file=indexfile)
877
878        # read values of index keys
879        index_vals = []
880        for key in index_keys:
881            index_vals.append(grib_index_get(iid, key))
882            print index_vals[-1]
883            # index_vals looks for example like:
884            # index_vals[0]: ('20171106', '20171107', '20171108') ; date
885            # index_vals[1]: ('0', '1200', '1800', '600') ; time
886            # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
887
888        fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
889                 '17':None, '19':None, '21':None, '22':None, '20':None}
890
891        for prod in product(*index_vals):
892            # flag for Fortran program CONVERT2, initially False
893            convertFlag = False
894            print 'current prod: ', prod
895            # e.g. prod = ('20170505', '0', '12')
896            #             (  date    ,time, step)
897            # per date e.g. time = 0, 600, 1200, 1800
898            # per time e.g. step = 0, 3, 6, 9, 12
899            for i in range(len(index_keys)):
900                grib_index_select(iid, index_keys[i], prod[i])
901
902            # get first id from current product
903            gid = grib_new_from_index(iid)
904
905            # if there is data for this product combination
906            # prepare some date and time parameter before reading the data
907            if gid is not None:
908                # Fortran program CONVERT2 is only done if gid at this time is
909                # not None, therefore save information in convertFlag
910                convertFlag = True
911                # remove old fort.* files and open new ones
912                for k, f in fdict.iteritems():
913                    silent_remove(c.inputdir + "/fort." + k)
914                    fdict[k] = open(c.inputdir + '/fort.' + k, 'w')
915
916                cdate = str(grib_get(gid, 'date'))
917                time = grib_get(gid, 'time')
918                step = grib_get(gid, 'step')
919                # create correct timestamp from the three time informations
920                # date, time, step
921                timestamp = datetime.strptime(cdate + '{:0>2}'.format(time/100),
922                                              '%Y%m%d%H')
923                timestamp += timedelta(hours=int(step))
924
925                cdateH = datetime.strftime(timestamp, '%Y%m%d%H')
926
927                if c.basetime is not None:
928                    slimit = datetime.strptime(c.start_date + '00', '%Y%m%d%H')
929                    bt = '23'
930                    if c.basetime == '00':
931                        bt = '00'
932                        slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
933                            - timedelta(hours=12-int(c.dtime))
934                    if c.basetime == '12':
935                        bt = '12'
936                        slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
937                            - timedelta(hours=12-int(c.dtime))
938
939                    elimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')
940
941                    if timestamp < slimit or timestamp > elimit:
942                        continue
943
944            try:
945                if c.wrf == '1':
946                    if 'olddate' not in locals():
947                        fwrf = open(c.outputdir + '/WRF' + cdate +
948                                    '.{:0>2}'.format(time) + '.000.grb2', 'w')
949                        olddate = cdate[:]
950                    else:
951                        if cdate != olddate:
952                            fwrf = open(c.outputdir + '/WRF' + cdate +
953                                        '.{:0>2}'.format(time) + '.000.grb2',
954                                        'w')
955                            olddate = cdate[:]
956            except AttributeError:
957                pass
958
959            # helper variable to remember which fields are already used.
960            savedfields = []
961            while 1:
962                if gid is None:
963                    break
964                paramId = grib_get(gid, 'paramId')
965                gridtype = grib_get(gid, 'gridType')
966                levtype = grib_get(gid, 'typeOfLevel')
967                if paramId == 133 and gridtype == 'reduced_gg':
968                # Relative humidity (Q.grb) is used as a template only
969                # so we need the first we "meet"
970                    with open(c.inputdir + '/fort.18', 'w') as fout:
971                        grib_write(gid, fout)
972                elif paramId == 131 or paramId == 132:
973                    grib_write(gid, fdict['10'])
974                elif paramId == 130:
975                    grib_write(gid, fdict['11'])
976                elif paramId == 133 and gridtype != 'reduced_gg':
977                    grib_write(gid, fdict['17'])
978                elif paramId == 152:
979                    grib_write(gid, fdict['12'])
980                elif paramId == 155 and gridtype == 'sh':
981                    grib_write(gid, fdict['13'])
982                elif  paramId in [129, 138, 155] and levtype == 'hybrid' \
983                        and c.wrf == '1':
984                    pass
985                elif paramId == 246 or paramId == 247:
986                    # cloud liquid water and ice
987                    if paramId == 246:
988                        clwc = grib_get_values(gid)
989                    else:
990                        clwc += grib_get_values(gid)
991                        grib_set_values(gid, clwc)
992                        grib_set(gid, 'paramId', 201031)
993                        grib_write(gid, fdict['22'])
994                elif paramId == 135:
995                    grib_write(gid, fdict['19'])
996                elif paramId == 77:
997                    grib_write(gid, fdict['21'])
998                else:
999                    if paramId not in savedfields:
1000                        grib_write(gid, fdict['16'])
1001                        savedfields.append(paramId)
1002                    else:
1003                        print 'duplicate ' + str(paramId) + ' not written'
1004
1005                try:
1006                    if c.wrf == '1':
1007                        if levtype == 'hybrid': # model layer
1008                            if paramId in [129, 130, 131, 132, 133, 138, 155]:
1009                                grib_write(gid, fwrf)
1010                        else: # sfc layer
1011                            if paramId in wrfpars:
1012                                grib_write(gid, fwrf)
1013                except AttributeError:
1014                    pass
1015
1016                grib_release(gid)
1017                gid = grib_new_from_index(iid)
1018
1019            for f in fdict.values():
1020                f.close()
1021
1022            # call for CONVERT2 if flag is True
1023            if convertFlag:
1024                pwd = os.getcwd()
1025                os.chdir(c.inputdir)
1026                if os.stat('fort.21').st_size == 0 and int(c.eta) == 1:
1027                    print 'Parameter 77 (etadot) is missing, most likely it is \
1028                           not available for this type or date/time\n'
1029                    print 'Check parameters CLASS, TYPE, STREAM, START_DATE\n'
1030                    my_error(c, 'fort.21 is empty while parameter eta is set \
1031                                to 1 in CONTROL file')
1032
1033                # create the corresponding output file fort.15
1034                # (generated by CONVERT2) + fort.16 (paramId 167 and 168)
1035                p = subprocess.check_call(
1036                    [os.path.expandvars(os.path.expanduser(c.exedir)) +
1037                     '/CONVERT2'], shell=True)
1038                os.chdir(pwd)
1039
1040                # create final output filename, e.g. EN13040500 (ENYYMMDDHH)
1041                fnout = c.inputdir + '/' + c.prefix
1042                if c.maxstep > 12:
1043                    suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \
1044                             '.{:0>3}'.format(step)
1045                else:
1046                    suffix = cdateH[2:10]
1047                fnout += suffix
1048                print "outputfile = " + fnout
1049                self.outputfilelist.append(fnout) # needed for final processing
1050
1051                # create outputfile and copy all data from intermediate files
1052                # to the outputfile (final GRIB files)
1053                orolsm = os.path.basename(glob.glob(
1054                    c.inputdir + '/OG_OROLSM__SL.*.' + c.ppid + '*')[0])
1055                fluxfile = 'flux' + cdate[0:2] + suffix
1056                if c.cwc != '1':
1057                    flist = ['fort.15', fluxfile, 'fort.16', orolsm]
1058                else:
1059                    flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]
1060
1061                with open(fnout, 'wb') as fout:
1062                    for f in flist:
1063                        shutil.copyfileobj(open(c.inputdir + '/' + f, 'rb'), fout)
1064
1065                if c.omega == '1':
1066                    with open(c.outputdir + '/OMEGA', 'wb') as fout:
1067                        shutil.copyfileobj(
1068                            open(c.inputdir + '/fort.25', 'rb'), fout)
1069
1070        if c.wrf == '1':
1071            fwrf.close()
1072
1073        grib_index_release(iid)
1074
1075        return
1076
1077    def deacc_fluxes(self, inputfiles, c):
1078        '''
1079        @Description:
1080            Goes through all flux fields in ordered time and de-accumulate
1081            the fields. Afterwards the fields are disaggregated in time.
1082            Different versions of disaggregation is provided for rainfall
1083            data (darain, modified linear) and the surface fluxes and
1084            stress data (dapoly, cubic polynomial).
1085
1086        @Input:
1087            self: instance of EcFlexpart
1088                The current object of the class.
1089
1090            inputfiles: instance of UioFiles
1091                Contains a list of files.
1092
1093            c: instance of class ControlFile
1094                Contains all the parameters of CONTROL file, which are e.g.:
1095                DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
1096                STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
1097                LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
1098                OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
1099                ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
1100                MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
1101                DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
1102
1103                For more information about format and content of the parameter
1104                see documentation.
1105
1106        @Return:
1107            <nothing>
1108        '''
1109
1110        table128 = init128(c.ecmwfdatadir +
1111                           '/grib_templates/ecmwf_grib1_table_128')
1112        pars = to_param_id(self.params['OG_acc_SL'][0], table128)
1113        index_keys = ["date", "time", "step"]
1114        indexfile = c.inputdir + "/date_time_stepRange.idx"
1115        silent_remove(indexfile)
1116        grib = GribTools(inputfiles.files)
1117        # creates new index file
1118        iid = grib.index(index_keys=index_keys, index_file=indexfile)
1119
1120        # read values of index keys
1121        index_vals = []
1122        for key in index_keys:
1123            key_vals = grib_index_get(iid, key)
1124            print key_vals
1125            # have to sort the steps for disaggregation,
1126            # therefore convert to int first
1127            if key == 'step':
1128                key_vals = [int(k) for k in key_vals]
1129                key_vals.sort()
1130                key_vals = [str(k) for k in key_vals]
1131            index_vals.append(key_vals)
1132            # index_vals looks for example like:
1133            # index_vals[0]: ('20171106', '20171107', '20171108') ; date
1134            # index_vals[1]: ('0', '1200') ; time
1135            # index_vals[2]: (3', '6', '9', '12') ; stepRange
1136
1137        valsdict = {}
1138        svalsdict = {}
1139        stepsdict = {}
1140        for p in pars:
1141            valsdict[str(p)] = []
1142            svalsdict[str(p)] = []
1143            stepsdict[str(p)] = []
1144
1145        print 'maxstep: ', c.maxstep
1146
1147        for prod in product(*index_vals):
1148            # e.g. prod = ('20170505', '0', '12')
1149            #             (  date    ,time, step)
1150            # per date e.g. time = 0, 1200
1151            # per time e.g. step = 3, 6, 9, 12
1152            for i in range(len(index_keys)):
1153                grib_index_select(iid, index_keys[i], prod[i])
1154
1155            gid = grib_new_from_index(iid)
1156            if gid is not None:
1157                cdate = grib_get(gid, 'date')
1158                time = grib_get(gid, 'time')
1159                step = grib_get(gid, 'step')
1160                # date+time+step-2*dtime
1161                # (since interpolated value valid for step-2*dtime)
1162                sdate = datetime(year=cdate/10000,
1163                                 month=(cdate % 10000)/100,
1164                                 day=(cdate % 100),
1165                                 hour=time/100)
1166                fdate = sdate + timedelta(hours=step-2*int(c.dtime))
1167                sdates = sdate + timedelta(hours=step)
1168                elimit = None
1169            else:
1170                break
1171
1172            if c.maxstep > 12:
1173                fnout = c.inputdir + '/flux' + \
1174                    sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
1175                    '.{:0>3}'.format(step-2*int(c.dtime))
1176                gnout = c.inputdir + '/flux' + \
1177                    sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
1178                    '.{:0>3}'.format(step-int(c.dtime))
1179                hnout = c.inputdir + '/flux' + \
1180                    sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
1181                    '.{:0>3}'.format(step)
1182                g = open(gnout, 'w')
1183                h = open(hnout, 'w')
1184            else:
1185                fnout = c.inputdir + '/flux' + fdate.strftime('%Y%m%d%H')
1186                gnout = c.inputdir + '/flux' + (fdate +
1187                                                timedelta(hours=int(c.dtime))
1188                                               ).strftime('%Y%m%d%H')
1189                hnout = c.inputdir + '/flux' + sdates.strftime('%Y%m%d%H')
1190                g = open(gnout, 'w')
1191                h = open(hnout, 'w')
1192
1193            print "outputfile = " + fnout
1194            f = open(fnout, 'w')
1195
1196            # read message for message and store relevant data fields
1197            # data keywords are stored in pars
1198            while 1:
1199                if gid is None:
1200                    break
1201                cparamId = str(grib_get(gid, 'paramId'))
1202                step = grib_get(gid, 'step')
1203                atime = grib_get(gid, 'time')
1204                ni = grib_get(gid, 'Ni')
1205                nj = grib_get(gid, 'Nj')
1206                if cparamId in valsdict.keys():
1207                    values = grib_get_values(gid)
1208                    vdp = valsdict[cparamId]
1209                    svdp = svalsdict[cparamId]
1210                    sd = stepsdict[cparamId]
1211
1212                    if cparamId == '142' or cparamId == '143':
1213                        fak = 1. / 1000.
1214                    else:
1215                        fak = 3600.
1216
1217                    values = (np.reshape(values, (nj, ni))).flatten() / fak
1218                    vdp.append(values[:])  # save the accumulated values
1219                    if step <= int(c.dtime):
1220                        svdp.append(values[:] / int(c.dtime))
1221                    else:  # deaccumulate values
1222                        svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
1223
1224                    print(cparamId, atime, step, len(values),
1225                          values[0], np.std(values))
1226                    # save the 1/3-hourly or specific values
1227                    # svdp.append(values[:])
1228                    sd.append(step)
1229                    # len(svdp) correspond to the time
1230                    if len(svdp) >= 3:
1231                        if len(svdp) > 3:
1232                            if cparamId == '142' or cparamId == '143':
1233                                values = disaggregation.darain(svdp)
1234                            else:
1235                                values = disaggregation.dapoly(svdp)
1236
1237                            if not (step == c.maxstep and c.maxstep > 12 \
1238                                    or sdates == elimit):
1239                                vdp.pop(0)
1240                                svdp.pop(0)
1241                        else:
1242                            if c.maxstep > 12:
1243                                values = svdp[1]
1244                            else:
1245                                values = svdp[0]
1246
1247                        grib_set_values(gid, values)
1248                        if c.maxstep > 12:
1249                            grib_set(gid, 'step', max(0, step-2*int(c.dtime)))
1250                        else:
1251                            grib_set(gid, 'step', 0)
1252                            grib_set(gid, 'time', fdate.hour*100)
1253                            grib_set(gid, 'date', fdate.year*10000 +
1254                                     fdate.month*100+fdate.day)
1255                        grib_write(gid, f)
1256
1257                        if c.basetime is not None:
1258                            elimit = datetime.strptime(c.end_date +
1259                                                       c.basetime, '%Y%m%d%H')
1260                        else:
1261                            elimit = sdate + timedelta(2*int(c.dtime))
1262
1263                        # squeeze out information of last two steps contained
1264                        # in svdp
1265                        # if step+int(c.dtime) == c.maxstep and c.maxstep>12
1266                        # or sdates+timedelta(hours = int(c.dtime))
1267                        # >= elimit:
1268                        # Note that svdp[0] has not been popped in this case
1269
1270                        if step == c.maxstep and c.maxstep > 12 or \
1271                           sdates == elimit:
1272
1273                            values = svdp[3]
1274                            grib_set_values(gid, values)
1275                            grib_set(gid, 'step', 0)
1276                            truedatetime = fdate + timedelta(hours=
1277                                                             2*int(c.dtime))
1278                            grib_set(gid, 'time', truedatetime.hour * 100)
1279                            grib_set(gid, 'date', truedatetime.year * 10000 +
1280                                     truedatetime.month * 100 +
1281                                     truedatetime.day)
1282                            grib_write(gid, h)
1283
1284                            #values = (svdp[1]+svdp[2])/2.
1285                            if cparamId == '142' or cparamId == '143':
1286                                values = disaggregation.darain(list(reversed(svdp)))
1287                            else:
1288                                values = disaggregation.dapoly(list(reversed(svdp)))
1289
1290                            grib_set(gid, 'step', 0)
1291                            truedatetime = fdate + timedelta(hours=int(c.dtime))
1292                            grib_set(gid, 'time', truedatetime.hour * 100)
1293                            grib_set(gid, 'date', truedatetime.year * 10000 +
1294                                     truedatetime.month * 100 +
1295                                     truedatetime.day)
1296                            grib_set_values(gid, values)
1297                            grib_write(gid, g)
1298
1299                    grib_release(gid)
1300
1301                    gid = grib_new_from_index(iid)
1302
1303            f.close()
1304            g.close()
1305            h.close()
1306
1307        grib_index_release(iid)
1308
1309        exit()
1310        return
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG