source: flex_extract.git/source/python/classes/EcFlexpart.py @ 65748f4

ctbtodev
Last change on this file since 65748f4 was 65748f4, checked in by Anne Philipp <anne.philipp@…>, 5 years ago

deleted unused variables, code style changes

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