source: flex_extract.git/python/EcFlexpart.py @ 2fb99de

devfeature/makefilesorigin/task/language-editingtask/language-editing
Last change on this file since 2fb99de was 2fb99de, checked in by Anne Philipp <anne.philipp@…>, 23 months ago

introduced config with path definitions and changed py files accordingly; Installation works; some tests were added for tarball making; Problems in submission to ecgate

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