source: flex_extract.git/python/EcFlexpart.py @ 54a8a01

ctbtodev
Last change on this file since 54a8a01 was 54a8a01, checked in by Anne Philipp <anne.philipp@…>, 6 years ago

restructuring, documentations and bug fixes

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