source: flex_extract.git/source/python/classes/EcFlexpart.py @ 4971f63

dev
Last change on this file since 4971f63 was 4971f63, checked in by Anne Philipp <anne.philipp@…>, 11 months ago

eliminated some redundancy and exchanged CONTROL2 with the config-filename

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