source: flex_extract.git/python/EcFlexpart.py @ 97e09f4

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

set path to grib1 table as global variable; test path

  • Property mode set to 100644
File size: 56.1 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(_config.PATH_GRIBTABLE)
855        wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
856                            stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
857                            table128)
858
859        index_keys = ["date", "time", "step"]
860        indexfile = c.inputdir + "/date_time_stepRange.idx"
861        silent_remove(indexfile)
862        grib = GribTools(inputfiles.files)
863        # creates new index file
864        iid = grib.index(index_keys=index_keys, index_file=indexfile)
865
866        # read values of index keys
867        index_vals = []
868        for key in index_keys:
869            index_vals.append(grib_index_get(iid, key))
870            print index_vals[-1]
871            # index_vals looks for example like:
872            # index_vals[0]: ('20171106', '20171107', '20171108') ; date
873            # index_vals[1]: ('0', '1200', '1800', '600') ; time
874            # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
875
876        fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
877                 '17':None, '19':None, '21':None, '22':None, '20':None}
878
879        for prod in product(*index_vals):
880            # flag for Fortran program CONVERT2 and file merging
881            convertFlag = False
882            print 'current prod: ', prod
883            # e.g. prod = ('20170505', '0', '12')
884            #             (  date    ,time, step)
885            # per date e.g. time = 0, 600, 1200, 1800
886            # per time e.g. step = 0, 3, 6, 9, 12
887            for i in range(len(index_keys)):
888                grib_index_select(iid, index_keys[i], prod[i])
889
890            # get first id from current product
891            gid = grib_new_from_index(iid)
892
893            # if there is data for this product combination
894            # prepare some date and time parameter before reading the data
895            if gid is not None:
896                # Combine all temporary data files into final grib file if
897                # gid is at least one time not None. Therefore set convertFlag
898                # to save information. The fortran program CONVERT2 is also
899                # only done if convertFlag is True
900                convertFlag = True
901                # remove old fort.* files and open new ones
902                # they are just valid for a single product
903                for k, f in fdict.iteritems():
904                    silent_remove(c.inputdir + "/fort." + k)
905                    fdict[k] = open(c.inputdir + '/fort.' + k, 'w')
906
907                cdate = str(grib_get(gid, 'date'))
908                time = grib_get(gid, 'time')
909                step = grib_get(gid, 'step')
910                # create correct timestamp from the three time informations
911                # date, time, step
912                timestamp = datetime.strptime(cdate + '{:0>2}'.format(time/100),
913                                              '%Y%m%d%H')
914                timestamp += timedelta(hours=int(step))
915                cdateH = datetime.strftime(timestamp, '%Y%m%d%H')
916
917                if c.basetime is not None:
918                    slimit = datetime.strptime(c.start_date + '00', '%Y%m%d%H')
919                    bt = '23'
920                    if c.basetime == '00':
921                        bt = '00'
922                        slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
923                            - timedelta(hours=12-int(c.dtime))
924                    if c.basetime == '12':
925                        bt = '12'
926                        slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
927                            - timedelta(hours=12-int(c.dtime))
928
929                    elimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')
930
931                    if timestamp < slimit or timestamp > elimit:
932                        continue
933
934            try:
935                if c.wrf == '1':
936                    if 'olddate' not in locals():
937                        fwrf = open(c.outputdir + '/WRF' + cdate +
938                                    '.{:0>2}'.format(time) + '.000.grb2', 'w')
939                        olddate = cdate[:]
940                    else:
941                        if cdate != olddate:
942                            fwrf = open(c.outputdir + '/WRF' + cdate +
943                                        '.{:0>2}'.format(time) + '.000.grb2',
944                                        'w')
945                            olddate = cdate[:]
946            except AttributeError:
947                pass
948
949            # helper variable to remember which fields are already used.
950            savedfields = []
951            while 1:
952                if gid is None:
953                    break
954                paramId = grib_get(gid, 'paramId')
955                gridtype = grib_get(gid, 'gridType')
956                levtype = grib_get(gid, 'typeOfLevel')
957                if paramId == 133 and gridtype == 'reduced_gg':
958                # Specific humidity (Q.grb) is used as a template only
959                # so we need the first we "meet"
960                    with open(c.inputdir + '/fort.18', 'w') as fout:
961                        grib_write(gid, fout)
962                elif paramId == 131 or paramId == 132:
963                    grib_write(gid, fdict['10'])
964                elif paramId == 130:
965                    grib_write(gid, fdict['11'])
966                elif paramId == 133 and gridtype != 'reduced_gg':
967                    grib_write(gid, fdict['17'])
968                elif paramId == 152:
969                    grib_write(gid, fdict['12'])
970                elif paramId == 155 and gridtype == 'sh':
971                    grib_write(gid, fdict['13'])
972                elif  paramId in [129, 138, 155] and levtype == 'hybrid' \
973                        and c.wrf == '1':
974                    pass
975                elif paramId == 246 or paramId == 247:
976                    # cloud liquid water and ice
977                    if paramId == 246:
978                        clwc = grib_get_values(gid)
979                    else:
980                        clwc += grib_get_values(gid)
981                        grib_set_values(gid, clwc)
982                        grib_set(gid, 'paramId', 201031)
983                        grib_write(gid, fdict['22'])
984                elif paramId == 135:
985                    grib_write(gid, fdict['19'])
986                elif paramId == 77:
987                    grib_write(gid, fdict['21'])
988                else:
989                    if paramId not in savedfields:
990                        grib_write(gid, fdict['16'])
991                        savedfields.append(paramId)
992                    else:
993                        print 'duplicate ' + str(paramId) + ' not written'
994
995                try:
996                    if c.wrf == '1':
997                        if levtype == 'hybrid': # model layer
998                            if paramId in [129, 130, 131, 132, 133, 138, 155]:
999                                grib_write(gid, fwrf)
1000                        else: # sfc layer
1001                            if paramId in wrfpars:
1002                                grib_write(gid, fwrf)
1003                except AttributeError:
1004                    pass
1005
1006                grib_release(gid)
1007                gid = grib_new_from_index(iid)
1008
1009            for f in fdict.values():
1010                f.close()
1011
1012            # call for CONVERT2 if flag is True
1013            if convertFlag:
1014                pwd = os.getcwd()
1015                os.chdir(c.inputdir)
1016                if os.stat('fort.21').st_size == 0 and int(c.eta) == 1:
1017                    print 'Parameter 77 (etadot) is missing, most likely it is \
1018                           not available for this type or date/time\n'
1019                    print 'Check parameters CLASS, TYPE, STREAM, START_DATE\n'
1020                    my_error(c.mailfail, 'fort.21 is empty while parameter eta \
1021                             is set to 1 in CONTROL file')
1022
1023                # create the corresponding output file fort.15
1024                # (generated by CONVERT2) + fort.16 (paramId 167 and 168)
1025                p = subprocess.check_call(
1026                    [os.path.expandvars(os.path.expanduser(c.exedir)) +
1027                     '/CONVERT2'], shell=True)
1028                os.chdir(pwd)
1029
1030                # create final output filename, e.g. EN13040500 (ENYYMMDDHH)
1031                fnout = c.inputdir + '/' + c.prefix
1032                if c.maxstep > 12:
1033                    suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \
1034                             '.{:0>3}'.format(step)
1035                else:
1036                    suffix = cdateH[2:10]
1037                fnout += suffix
1038                print "outputfile = " + fnout
1039                self.outputfilelist.append(fnout) # needed for final processing
1040
1041                # create outputfile and copy all data from intermediate files
1042                # to the outputfile (final GRIB files)
1043                orolsm = os.path.basename(glob.glob(
1044                    c.inputdir + '/OG_OROLSM__SL.*.' + c.ppid + '*')[0])
1045                fluxfile = 'flux' + cdate[0:2] + suffix
1046                if c.cwc != '1':
1047                    flist = ['fort.15', fluxfile, 'fort.16', orolsm]
1048                else:
1049                    flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]
1050
1051                with open(fnout, 'wb') as fout:
1052                    for f in flist:
1053                        shutil.copyfileobj(
1054                            open(c.inputdir + '/' + f, 'rb'), fout)
1055
1056                if c.omega == '1':
1057                    with open(c.outputdir + '/OMEGA', 'wb') as fout:
1058                        shutil.copyfileobj(
1059                            open(c.inputdir + '/fort.25', 'rb'), fout)
1060
1061        if hasattr(c, 'wrf') and c.wrf == '1':
1062            fwrf.close()
1063
1064        grib_index_release(iid)
1065
1066        return
1067
1068    def deacc_fluxes(self, inputfiles, c):
1069        '''
1070        @Description:
1071            Goes through all flux fields in ordered time and de-accumulate
1072            the fields. Afterwards the fields are disaggregated in time.
1073            Different versions of disaggregation is provided for rainfall
1074            data (darain, modified linear) and the surface fluxes and
1075            stress data (dapoly, cubic polynomial).
1076
1077        @Input:
1078            self: instance of EcFlexpart
1079                The current object of the class.
1080
1081            inputfiles: instance of UioFiles
1082                Contains a list of files.
1083
1084            c: instance of class ControlFile
1085                Contains all the parameters of CONTROL file, which are e.g.:
1086                DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
1087                STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
1088                LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
1089                OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
1090                ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
1091                MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
1092                DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
1093
1094                For more information about format and content of the parameter
1095                see documentation.
1096
1097        @Return:
1098            <nothing>
1099        '''
1100
1101        table128 = init128(_config.PATH_GRIBTABLE)
1102        pars = to_param_id(self.params['OG_acc_SL'][0], table128)
1103        index_keys = ["date", "time", "step"]
1104        indexfile = c.inputdir + "/date_time_stepRange.idx"
1105        silent_remove(indexfile)
1106        grib = GribTools(inputfiles.files)
1107        # creates new index file
1108        iid = grib.index(index_keys=index_keys, index_file=indexfile)
1109
1110        # read values of index keys
1111        index_vals = []
1112        for key in index_keys:
1113            key_vals = grib_index_get(iid, key)
1114            print key_vals
1115            # have to sort the steps for disaggregation,
1116            # therefore convert to int first
1117            if key == 'step':
1118                key_vals = [int(k) for k in key_vals]
1119                key_vals.sort()
1120                key_vals = [str(k) for k in key_vals]
1121            index_vals.append(key_vals)
1122            # index_vals looks for example like:
1123            # index_vals[0]: ('20171106', '20171107', '20171108') ; date
1124            # index_vals[1]: ('0', '1200') ; time
1125            # index_vals[2]: (3', '6', '9', '12') ; stepRange
1126
1127        valsdict = {}
1128        svalsdict = {}
1129        stepsdict = {}
1130        for p in pars:
1131            valsdict[str(p)] = []
1132            svalsdict[str(p)] = []
1133            stepsdict[str(p)] = []
1134
1135        print 'maxstep: ', c.maxstep
1136
1137        for prod in product(*index_vals):
1138            # e.g. prod = ('20170505', '0', '12')
1139            #             (  date    ,time, step)
1140            # per date e.g. time = 0, 1200
1141            # per time e.g. step = 3, 6, 9, 12
1142            for i in range(len(index_keys)):
1143                grib_index_select(iid, index_keys[i], prod[i])
1144
1145            gid = grib_new_from_index(iid)
1146            if gid is not None:
1147                cdate = grib_get(gid, 'date')
1148                time = grib_get(gid, 'time')
1149                step = grib_get(gid, 'step')
1150                # date+time+step-2*dtime
1151                # (since interpolated value valid for step-2*dtime)
1152                sdate = datetime(year=cdate/10000,
1153                                 month=(cdate % 10000)/100,
1154                                 day=(cdate % 100),
1155                                 hour=time/100)
1156                fdate = sdate + timedelta(hours=step-2*int(c.dtime))
1157                sdates = sdate + timedelta(hours=step)
1158                elimit = None
1159            else:
1160                break
1161
1162            if c.maxstep > 12:
1163                fnout = c.inputdir + '/flux' + \
1164                    sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
1165                    '.{:0>3}'.format(step-2*int(c.dtime))
1166                gnout = c.inputdir + '/flux' + \
1167                    sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
1168                    '.{:0>3}'.format(step-int(c.dtime))
1169                hnout = c.inputdir + '/flux' + \
1170                    sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
1171                    '.{:0>3}'.format(step)
1172                g = open(gnout, 'w')
1173                h = open(hnout, 'w')
1174            else:
1175                fnout = c.inputdir + '/flux' + fdate.strftime('%Y%m%d%H')
1176                gnout = c.inputdir + '/flux' + (fdate +
1177                                                timedelta(hours=int(c.dtime))
1178                                               ).strftime('%Y%m%d%H')
1179                hnout = c.inputdir + '/flux' + sdates.strftime('%Y%m%d%H')
1180                g = open(gnout, 'w')
1181                h = open(hnout, 'w')
1182
1183            print "outputfile = " + fnout
1184            f = open(fnout, 'w')
1185
1186            # read message for message and store relevant data fields
1187            # data keywords are stored in pars
1188            while 1:
1189                if gid is None:
1190                    break
1191                cparamId = str(grib_get(gid, 'paramId'))
1192                step = grib_get(gid, 'step')
1193                atime = grib_get(gid, 'time')
1194                ni = grib_get(gid, 'Ni')
1195                nj = grib_get(gid, 'Nj')
1196                if cparamId in valsdict.keys():
1197                    values = grib_get_values(gid)
1198                    vdp = valsdict[cparamId]
1199                    svdp = svalsdict[cparamId]
1200                    sd = stepsdict[cparamId]
1201
1202                    if cparamId == '142' or cparamId == '143':
1203                        fak = 1. / 1000.
1204                    else:
1205                        fak = 3600.
1206
1207                    values = (np.reshape(values, (nj, ni))).flatten() / fak
1208                    vdp.append(values[:])  # save the accumulated values
1209                    if step <= int(c.dtime):
1210                        svdp.append(values[:] / int(c.dtime))
1211                    else:  # deaccumulate values
1212                        svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
1213
1214                    print(cparamId, atime, step, len(values),
1215                          values[0], np.std(values))
1216                    # save the 1/3-hourly or specific values
1217                    # svdp.append(values[:])
1218                    sd.append(step)
1219                    # len(svdp) correspond to the time
1220                    if len(svdp) >= 3:
1221                        if len(svdp) > 3:
1222                            if cparamId == '142' or cparamId == '143':
1223                                values = disaggregation.darain(svdp)
1224                            else:
1225                                values = disaggregation.dapoly(svdp)
1226
1227                            if not (step == c.maxstep and c.maxstep > 12 \
1228                                    or sdates == elimit):
1229                                vdp.pop(0)
1230                                svdp.pop(0)
1231                        else:
1232                            if c.maxstep > 12:
1233                                values = svdp[1]
1234                            else:
1235                                values = svdp[0]
1236
1237                        grib_set_values(gid, values)
1238                        if c.maxstep > 12:
1239                            grib_set(gid, 'step', max(0, step-2*int(c.dtime)))
1240                        else:
1241                            grib_set(gid, 'step', 0)
1242                            grib_set(gid, 'time', fdate.hour*100)
1243                            grib_set(gid, 'date', fdate.year*10000 +
1244                                     fdate.month*100+fdate.day)
1245                        grib_write(gid, f)
1246
1247                        if c.basetime is not None:
1248                            elimit = datetime.strptime(c.end_date +
1249                                                       c.basetime, '%Y%m%d%H')
1250                        else:
1251                            elimit = sdate + timedelta(2*int(c.dtime))
1252
1253                        # squeeze out information of last two steps contained
1254                        # in svdp
1255                        # if step+int(c.dtime) == c.maxstep and c.maxstep>12
1256                        # or sdates+timedelta(hours = int(c.dtime))
1257                        # >= elimit:
1258                        # Note that svdp[0] has not been popped in this case
1259
1260                        if step == c.maxstep and c.maxstep > 12 or \
1261                           sdates == elimit:
1262
1263                            values = svdp[3]
1264                            grib_set_values(gid, values)
1265                            grib_set(gid, 'step', 0)
1266                            truedatetime = fdate + timedelta(hours=
1267                                                             2*int(c.dtime))
1268                            grib_set(gid, 'time', truedatetime.hour * 100)
1269                            grib_set(gid, 'date', truedatetime.year * 10000 +
1270                                     truedatetime.month * 100 +
1271                                     truedatetime.day)
1272                            grib_write(gid, h)
1273
1274                            #values = (svdp[1]+svdp[2])/2.
1275                            if cparamId == '142' or cparamId == '143':
1276                                values = disaggregation.darain(list(reversed(svdp)))
1277                            else:
1278                                values = disaggregation.dapoly(list(reversed(svdp)))
1279
1280                            grib_set(gid, 'step', 0)
1281                            truedatetime = fdate + timedelta(hours=int(c.dtime))
1282                            grib_set(gid, 'time', truedatetime.hour * 100)
1283                            grib_set(gid, 'date', truedatetime.year * 10000 +
1284                                     truedatetime.month * 100 +
1285                                     truedatetime.day)
1286                            grib_set_values(gid, values)
1287                            grib_write(gid, g)
1288
1289                    grib_release(gid)
1290
1291                    gid = grib_new_from_index(iid)
1292
1293            f.close()
1294            g.close()
1295            h.close()
1296
1297        grib_index_release(iid)
1298
1299        return
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG