source: flex_extract.git/source/python/classes/EcFlexpart.py @ 295ff45

dev
Last change on this file since 295ff45 was 295ff45, checked in by Anne Philipp <anne.philipp@…>, 9 months ago

added request output option also for basetime requests; changed output to csv format

  • Property mode set to 100644
File size: 55.5 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
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 for different requests and therefore
500        # will 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 fields 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'] = self.inputdir + "/" + ftype + \
538                        pk + '.' + self.dates.split('/')[0] + '.' + \
539                        str(os.getppid()) + '.' + str(os.getpid()) + ".grb"
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                    self._start_retrievement(request, retr_param_dict)
567    # ------  operational path  ------------------------------------------------
568                else:
569                    # check if mars job requests fields beyond basetime.
570                    # If yes eliminate those fields since they may not
571                    # be accessible with user's credentials
572                    if 'by' in retr_param_dict['step']:
573                        sm1 = 2
574                    else:
575                        sm1 = -1
576
577                    if 'by' in retr_param_dict['time']:
578                        tm1 = 2
579                    else:
580                        tm1 = -1
581
582                    maxdate = datetime.strptime(retr_param_dict['date'].split('/')[-1] +
583                                                retr_param_dict['time'].split('/')[tm1],
584                                                '%Y%m%d%H')
585                    maxtime = maxdate + \
586                        timedelta(hours=int(retr_param_dict['step'].split('/')[sm1]))
587
588                    elimit = datetime.strptime(retr_param_dict['date'].split('/')[-1] +
589                                               self.basetime, '%Y%m%d%H')
590
591                    if self.basetime == '12':
592                        # --------------  flux data ----------------------------
593                        if 'acc' in pk:
594
595                        # Strategy:
596                        # if maxtime-elimit >= 24h reduce date by 1,
597                        # if 12h <= maxtime-elimit<12h reduce time for last date
598                        # if maxtime-elimit<12h reduce step for last time
599                        # A split of the MARS job into 2 is likely necessary.
600
601
602                            startdate = retr_param_dict['date'].split('/')[0]
603                            enddate = datetime.strftime(elimit - t24h,'%Y%m%d')
604                            retr_param_dict['date'] = '/'.join([startdate, 'to', enddate])
605
606                            self._start_retrievement(request, retr_param_dict)
607
608                            retr_param_dict['date'] = \
609                                datetime.strftime(elimit - t12h, '%Y%m%d')
610                            retr_param_dict['time'] = '00'
611                            retr_param_dict['target'] = \
612                                self._mk_targetname(ftype, pk,
613                                                    retr_param_dict['date'])
614
615                            self._start_retrievement(request, retr_param_dict)
616
617                        # --------------  non flux data ------------------------
618                        else:
619                            self._start_retrievement(request, retr_param_dict)
620
621                    else: # basetime = 0
622                        retr_param_dict['date'] = \
623                            datetime.strftime(elimit - t24h, '%Y%m%d')
624
625                        timesave = ''.join(retr_param_dict['time'])
626
627                        if '/' in retr_param_dict['time']:
628                            times = retr_param_dict['time'].split('/')
629                            steps = retr_param_dict['step'].split('/')
630                            while (pk != 'OG_OROLSM__SL' and
631                                   'acc' not in pk and
632                                   (int(times[0]) + int(steps[0])) <= 12):
633                                times = times[1:]
634
635                            if len(times) > 1:
636                                retr_param_dict['time'] = '/'.join(times)
637                            else:
638                                retr_param_dict['time'] = times[0]
639
640                        self._start_retrievement(request, retr_param_dict)
641
642                        if (pk != 'OG_OROLSM__SL' and
643                            int(retr_param_dict['step'].split('/')[0]) == 0 and
644                            int(timesave.split('/')[0]) == 0):
645
646                            retr_param_dict['date'] = \
647                                datetime.strftime(elimit, '%Y%m%d')
648                            retr_param_dict['time'] = '00'
649                            retr_param_dict['step'] = '000'
650                            retr_param_dict['target'] = \
651                                self._mk_targetname(ftype, pk,
652                                                    retr_param_dict['date'])
653
654                            self._start_retrievement(request, retr_param_dict)
655
656        if request == 0 or request == 2:
657            print('MARS retrieve done ... ')
658        elif request == 1:
659            print('MARS request printed ...')
660
661        return
662
663
664    def process_output(self, c):
665        '''
666        @Description:
667            The grib files are postprocessed depending on the selection in
668            CONTROL file. The resulting files are moved to the output
669            directory if its not equla to the input directory.
670            The following modifications might be done if
671            properly switched in CONTROL file:
672            GRIB2 - Conversion to GRIB2
673            ECTRANS - Transfer of files to gateway server
674            ECSTORAGE - Storage at ECMWF server
675            GRIB2FLEXPART - Conversion of GRIB files to FLEXPART binary format
676
677        @Input:
678            self: instance of EcFlexpart
679                The current object of the class.
680
681            c: instance of class ControlFile
682                Contains all the parameters of CONTROL file, which are e.g.:
683                DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
684                STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
685                LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
686                OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
687                ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
688                MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
689                DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
690
691                For more information about format and content of the parameter
692                see documentation.
693
694        @Return:
695            <nothing>
696
697        '''
698
699        print('\n\nPostprocessing:\n Format: {}\n'.format(c.format))
700
701        if not c.ecapi:
702            print('ecstorage: {}\n ecfsdir: {}\n'.
703                  format(c.ecstorage, c.ecfsdir))
704            #if not hasattr(c, 'gateway'):
705            #    c.gateway = os.getenv('GATEWAY')
706            #if not hasattr(c, 'destination'):
707            #    c.destination = os.getenv('DESTINATION')
708            print('ectrans: {}\n gateway: {}\n destination: {}\n '
709                  .format(c.ectrans, c.gateway, c.destination))
710
711        print('Output filelist: \n')
712        print(self.outputfilelist)
713
714        if c.format.lower() == 'grib2':
715            for ofile in self.outputfilelist:
716                p = subprocess.check_call(['grib_set', '-s', 'edition=2, \
717                                           productDefinitionTemplateNumber=8',
718                                           ofile, ofile + '_2'])
719                p = subprocess.check_call(['mv', ofile + '_2', ofile])
720
721        if c.ectrans and not c.ecapi:
722            for ofile in self.outputfilelist:
723                p = subprocess.check_call(['ectrans', '-overwrite', '-gateway',
724                                           c.gateway, '-remote', c.destination,
725                                           '-source', ofile])
726                #print('ectrans:', p)
727
728        if c.ecstorage and not c.ecapi:
729            for ofile in self.outputfilelist:
730                p = subprocess.check_call(['ecp', '-o', ofile,
731                                           os.path.expandvars(c.ecfsdir)])
732
733        if c.outputdir != c.inputdir:
734            for ofile in self.outputfilelist:
735                p = subprocess.check_call(['mv', ofile, c.outputdir])
736
737        # prepare environment for the grib2flexpart run
738        # to convert grib to flexpart binary
739        if c.grib2flexpart:
740
741            # generate AVAILABLE file
742            # Example of AVAILABLE file data:
743            # 20131107 000000      EN13110700              ON DISC
744            clist = []
745            for ofile in self.outputfilelist:
746                fname = ofile.split('/')
747                if '.' in fname[-1]:
748                    l = fname[-1].split('.')
749                    timestamp = datetime.strptime(l[0][-6:] + l[1],
750                                                  '%y%m%d%H')
751                    timestamp += timedelta(hours=int(l[2]))
752                    cdate = datetime.strftime(timestamp, '%Y%m%d')
753                    chms = datetime.strftime(timestamp, '%H%M%S')
754                else:
755                    cdate = '20' + fname[-1][-8:-2]
756                    chms = fname[-1][-2:] + '0000'
757                clist.append(cdate + ' ' + chms + ' '*6 +
758                             fname[-1] + ' '*14 + 'ON DISC')
759            clist.sort()
760            with open(c.outputdir + '/' + 'AVAILABLE', 'w') as f:
761                f.write('\n'.join(clist) + '\n')
762
763            # generate pathnames file
764            pwd = os.path.abspath(c.outputdir)
765            with open(pwd + '/pathnames', 'w') as f:
766                f.write(pwd + '/Options/\n')
767                f.write(pwd + '/\n')
768                f.write(pwd + '/\n')
769                f.write(pwd + '/AVAILABLE\n')
770                f.write(' = == = == = == = == = == ==  = \n')
771
772            # create Options dir if necessary
773            if not os.path.exists(pwd + '/Options'):
774                os.makedirs(pwd+'/Options')
775
776            # read template COMMAND file
777            with open(os.path.expandvars(os.path.expanduser(
778                c.flexpart_root_scripts)) + '/../Options/COMMAND', 'r') as f:
779                lflist = f.read().split('\n')
780
781            # find index of list where to put in the
782            # date and time information
783            # usually after the LDIRECT parameter
784            i = 0
785            for l in lflist:
786                if 'LDIRECT' in l.upper():
787                    break
788                i += 1
789
790            # insert the date and time information of run start and end
791            # into the list of lines of COMMAND file
792            lflist = lflist[:i+1] + \
793                     [clist[0][:16], clist[-1][:16]] + \
794                     lflist[i+3:]
795
796            # write the new COMMAND file
797            with open(pwd + '/Options/COMMAND', 'w') as g:
798                g.write('\n'.join(lflist) + '\n')
799
800            # change to outputdir and start the grib2flexpart run
801            # afterwards switch back to the working dir
802            os.chdir(c.outputdir)
803            p = subprocess.check_call([
804                os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))
805                + '/../FLEXPART_PROGRAM/grib2flexpart', 'useAvailable', '.'])
806            os.chdir(pwd)
807
808        return
809
810    def create(self, inputfiles, c):
811        '''
812        @Description:
813            This method is based on the ECMWF example index.py
814            https://software.ecmwf.int/wiki/display/GRIB/index.py
815
816            An index file will be created which depends on the combination
817            of "date", "time" and "stepRange" values. This is used to iterate
818            over all messages in each grib file which were passed through the
819            parameter "inputfiles" to seperate specific parameters into fort.*
820            files. Afterwards the FORTRAN program is called to convert
821            the data fields all to the same grid and put them in one file
822            per unique time step (combination of "date", "time" and
823            "stepRange").
824
825        @Input:
826            self: instance of EcFlexpart
827                The current object of the class.
828
829            inputfiles: instance of UioFiles
830                Contains a list of files.
831
832            c: instance of class ControlFile
833                Contains all the parameters of CONTROL files, which are e.g.:
834                DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
835                STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
836                LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
837                OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
838                ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
839                MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
840                DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
841
842                For more information about format and content of the parameter
843                see documentation.
844
845        @Return:
846            <nothing>
847        '''
848
849        table128 = init128(_config.PATH_GRIBTABLE)
850        wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
851                            stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
852                            table128)
853
854        index_keys = ["date", "time", "step"]
855        indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX)
856        silent_remove(indexfile)
857        grib = GribTools(inputfiles.files)
858        # creates new index file
859        iid = grib.index(index_keys=index_keys, index_file=indexfile)
860
861        # read values of index keys
862        index_vals = []
863        for key in index_keys:
864            index_vals.append(grib_index_get(iid, key))
865            print(index_vals[-1])
866            # index_vals looks for example like:
867            # index_vals[0]: ('20171106', '20171107', '20171108') ; date
868            # index_vals[1]: ('0', '1200', '1800', '600') ; time
869            # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
870
871        fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
872                 '17':None, '19':None, '21':None, '22':None, '20':None}
873
874        for prod in product(*index_vals):
875            # e.g. prod = ('20170505', '0', '12')
876            #             (  date    ,time, step)
877            # per date e.g. time = 0, 1200
878            # per time e.g. step = 3, 6, 9, 12
879            # flag for Fortran program and file merging
880            convertFlag = False
881            print('current prod: ', prod)
882            # e.g. prod = ('20170505', '0', '12')
883            #             (  date    ,time, step)
884            # per date e.g. time = 0, 600, 1200, 1800
885            # per time e.g. step = 0, 3, 6, 9, 12
886            for i in range(len(index_keys)):
887                grib_index_select(iid, index_keys[i], prod[i])
888
889            # get first id from current product
890            gid = grib_new_from_index(iid)
891
892            # if there is data for this product combination
893            # prepare some date and time parameter before reading the data
894            if gid is not None:
895                # Combine all temporary data files into final grib file if
896                # gid is at least one time not None. Therefore set convertFlag
897                # to save information. The Fortran program is also
898                # only executed if convertFlag is True
899                convertFlag = True
900                # remove old fort.* files and open new ones
901                # they are just valid for a single product
902                for k, f in fdict.iteritems():
903                    fortfile = os.path.join(c.inputdir, 'fort.' + k)
904                    silent_remove(fortfile)
905                    fdict[k] = open(fortfile, 'w')
906
907                cdate = str(grib_get(gid, 'date'))
908                time = grib_get(gid, 'time')
909                step = grib_get(gid, 'step')
910                # create correct timestamp from the three time informations
911                # date, time, step
912                timestamp = datetime.strptime(cdate + '{:0>2}'.format(time/100),
913                                              '%Y%m%d%H')
914                timestamp += timedelta(hours=int(step))
915                cdateH = datetime.strftime(timestamp, '%Y%m%d%H')
916
917                if c.basetime:
918                    slimit = datetime.strptime(c.start_date + '00', '%Y%m%d%H')
919                    bt = '23'
920                    if c.basetime == '00':
921                        bt = '00'
922                        slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
923                            - timedelta(hours=12-int(c.dtime))
924                    if c.basetime == '12':
925                        bt = '12'
926                        slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
927                            - timedelta(hours=12-int(c.dtime))
928
929                    elimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')
930
931                    if timestamp < slimit or timestamp > elimit:
932                        continue
933            else:
934                pass
935
936            try:
937                if c.wrf:
938                    if 'olddate' not in locals() or cdate != olddate:
939                        fwrf = open(os.path.join(c.outputdir,
940                                    'WRF' + cdate + '.{:0>2}'.format(time) +
941                                    '.000.grb2'), 'w')
942                        olddate = cdate[:]
943            except AttributeError:
944                pass
945
946            # helper variable to remember which fields were already used.
947            savedfields = []
948            while 1:
949                if gid is None:
950                    break
951                paramId = grib_get(gid, 'paramId')
952                gridtype = grib_get(gid, 'gridType')
953                levtype = grib_get(gid, 'typeOfLevel')
954                if paramId == 133 and gridtype == 'reduced_gg':
955                # Specific humidity (Q.grb) is used as a template only
956                # so we need the first we "meet"
957                    with open(os.path.join(c.inputdir, 'fort.18'), 'w') as fout:
958                        grib_write(gid, fout)
959                elif paramId == 131 or paramId == 132:
960                    grib_write(gid, fdict['10'])
961                elif paramId == 130:
962                    grib_write(gid, fdict['11'])
963                elif paramId == 133 and gridtype != 'reduced_gg':
964                    grib_write(gid, fdict['17'])
965                elif paramId == 152:
966                    grib_write(gid, fdict['12'])
967                elif paramId == 155 and gridtype == 'sh':
968                    grib_write(gid, fdict['13'])
969                elif  paramId in [129, 138, 155] and levtype == 'hybrid' \
970                        and c.wrf:
971                    pass
972                elif paramId == 246 or paramId == 247:
973                    # cloud liquid water and ice
974                    if paramId == 246:
975                        clwc = grib_get_values(gid)
976                    else:
977                        clwc += grib_get_values(gid)
978                        grib_set_values(gid, clwc)
979                        grib_set(gid, 'paramId', 201031)
980                        grib_write(gid, fdict['22'])
981                elif paramId == 135:
982                    grib_write(gid, fdict['19'])
983                elif paramId == 77:
984                    grib_write(gid, fdict['21'])
985                else:
986                    if paramId not in savedfields:
987                        grib_write(gid, fdict['16'])
988                        savedfields.append(paramId)
989                    else:
990                        print('duplicate ' + str(paramId) + ' not written')
991
992                try:
993                    if c.wrf:
994                        # model layer
995                        if levtype == 'hybrid' and \
996                           paramId in [129, 130, 131, 132, 133, 138, 155]:
997                            grib_write(gid, fwrf)
998                        # sfc layer
999                        elif paramId in wrfpars:
1000                            grib_write(gid, fwrf)
1001                except AttributeError:
1002                    pass
1003
1004                grib_release(gid)
1005                gid = grib_new_from_index(iid)
1006
1007            for f in fdict.values():
1008                f.close()
1009
1010            # call for Fortran program if flag is True
1011            if convertFlag:
1012                pwd = os.getcwd()
1013                os.chdir(c.inputdir)
1014                if os.stat('fort.21').st_size == 0 and c.eta:
1015                    print('Parameter 77 (etadot) is missing, most likely it is \
1016                           not available for this type or date/time\n')
1017                    print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
1018                    my_error(c.mailfail, 'fort.21 is empty while parameter eta \
1019                             is set to 1 in CONTROL file')
1020
1021                # create the corresponding output file fort.15
1022                # (generated by Fortran program) + fort.16 (paramId 167 and 168)
1023                p = subprocess.check_call([os.path.join(
1024                    c.exedir, _config.FORTRAN_EXECUTABLE)], shell=True)
1025                os.chdir(pwd)
1026
1027                # create final output filename, e.g. EN13040500 (ENYYMMDDHH)
1028                fnout = os.path.join(c.inputdir, c.prefix)
1029                if c.maxstep > 12:
1030                    suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \
1031                             '.{:0>3}'.format(step)
1032                else:
1033                    suffix = cdateH[2:10]
1034                fnout += suffix
1035                print("outputfile = " + fnout)
1036                self.outputfilelist.append(fnout) # needed for final processing
1037
1038                # create outputfile and copy all data from intermediate files
1039                # to the outputfile (final GRIB files)
1040                orolsm = os.path.basename(glob.glob(
1041                    c.inputdir + '/OG_OROLSM__SL.*.' + c.ppid + '*')[0])
1042                fluxfile = 'flux' + cdate[0:2] + suffix
1043                if not c.cwc:
1044                    flist = ['fort.15', fluxfile, 'fort.16', orolsm]
1045                else:
1046                    flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]
1047
1048                with open(fnout, 'wb') as fout:
1049                    for f in flist:
1050                        shutil.copyfileobj(open(os.path.join(c.inputdir, f),
1051                                                'rb'), fout)
1052
1053                if c.omega:
1054                    with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout:
1055                        shutil.copyfileobj(
1056                            open(os.path.join(c.inputdir, 'fort.25'),
1057                                 'rb'), fout)
1058            else:
1059                pass
1060
1061        if c.wrf:
1062            fwrf.close()
1063
1064        grib_index_release(iid)
1065
1066        return
1067
1068    def deacc_fluxes(self, inputfiles, c):
1069        '''
1070        @Description:
1071            Goes through all flux fields in ordered time and de-accumulate
1072            the fields. Afterwards the fields are disaggregated in time.
1073            Different versions of disaggregation is provided for rainfall
1074            data (darain, modified linear) and the surface fluxes and
1075            stress data (dapoly, cubic polynomial).
1076
1077        @Input:
1078            self: instance of EcFlexpart
1079                The current object of the class.
1080
1081            inputfiles: instance of UioFiles
1082                Contains a list of files.
1083
1084            c: instance of class ControlFile
1085                Contains all the parameters of CONTROL file, which are e.g.:
1086                DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
1087                STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
1088                LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
1089                OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
1090                ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
1091                MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
1092                DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
1093
1094                For more information about format and content of the parameter
1095                see documentation.
1096
1097        @Return:
1098            <nothing>
1099        '''
1100
1101        table128 = init128(_config.PATH_GRIBTABLE)
1102        pars = to_param_id(self.params['OG_acc_SL'][0], table128)
1103        index_keys = ["date", "time", "step"]
1104        indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX)
1105        silent_remove(indexfile)
1106        grib = GribTools(inputfiles.files)
1107        # creates new index file
1108        iid = grib.index(index_keys=index_keys, index_file=indexfile)
1109
1110        # read values of index keys
1111        index_vals = []
1112        for key in index_keys:
1113            key_vals = grib_index_get(iid, key)
1114            print(key_vals)
1115            # have to sort the steps for disaggregation,
1116            # therefore convert to int first
1117            if key == 'step':
1118                key_vals = [int(k) for k in key_vals]
1119                key_vals.sort()
1120                key_vals = [str(k) for k in key_vals]
1121            index_vals.append(key_vals)
1122            # index_vals looks for example like:
1123            # index_vals[0]: ('20171106', '20171107', '20171108') ; date
1124            # index_vals[1]: ('0', '1200') ; time
1125            # index_vals[2]: (3', '6', '9', '12') ; stepRange
1126
1127        valsdict = {}
1128        svalsdict = {}
1129        stepsdict = {}
1130        for p in pars:
1131            valsdict[str(p)] = []
1132            svalsdict[str(p)] = []
1133            stepsdict[str(p)] = []
1134
1135        print('maxstep: ', c.maxstep)
1136
1137        for prod in product(*index_vals):
1138            # e.g. prod = ('20170505', '0', '12')
1139            #             (  date    ,time, step)
1140            # per date e.g. time = 0, 1200
1141            # per time e.g. step = 3, 6, 9, 12
1142            print('current prod: ', prod)
1143            for i in range(len(index_keys)):
1144                grib_index_select(iid, index_keys[i], prod[i])
1145
1146            # get first id from current product
1147            gid = grib_new_from_index(iid)
1148
1149            # if there is data for this product combination
1150            # prepare some date and time parameter before reading the data
1151            if gid is not None:
1152                cdate = grib_get(gid, 'date')
1153                time = grib_get(gid, 'time')
1154                step = grib_get(gid, 'step')
1155                # date+time+step-2*dtime
1156                # (since interpolated value valid for step-2*dtime)
1157                sdate = datetime(year=cdate/10000,
1158                                 month=(cdate % 10000)/100,
1159                                 day=(cdate % 100),
1160                                 hour=time/100)
1161                fdate = sdate + timedelta(hours=step-2*int(c.dtime))
1162                sdates = sdate + timedelta(hours=step)
1163                elimit = None
1164            else:
1165                break
1166
1167            if c.maxstep > 12:
1168                fnout = os.path.join(c.inputdir, 'flux' +
1169                                     sdate.strftime('%Y%m%d') +
1170                                     '.{:0>2}'.format(time/100) +
1171                                     '.{:0>3}'.format(step-2*int(c.dtime)))
1172                gnout = os.path.join(c.inputdir, 'flux' +
1173                                     sdate.strftime('%Y%m%d') +
1174                                     '.{:0>2}'.format(time/100) +
1175                                     '.{:0>3}'.format(step-int(c.dtime)))
1176                hnout = os.path.join(c.inputdir, 'flux' +
1177                                     sdate.strftime('%Y%m%d') +
1178                                     '.{:0>2}'.format(time/100) +
1179                                     '.{:0>3}'.format(step))
1180            else:
1181                fnout = os.path.join(c.inputdir, 'flux' +
1182                                     fdate.strftime('%Y%m%d%H'))
1183                gnout = os.path.join(c.inputdir, 'flux' +
1184                                     (fdate + timedelta(hours=int(c.dtime))).
1185                                     strftime('%Y%m%d%H'))
1186                hnout = os.path.join(c.inputdir, 'flux' +
1187                                     sdates.strftime('%Y%m%d%H'))
1188
1189            print("outputfile = " + fnout)
1190            f_handle = open(fnout, 'w')
1191            g_handle = open(gnout, 'w')
1192            h_handle = open(hnout, 'w')
1193
1194            # read message for message and store relevant data fields
1195            # data keywords are stored in pars
1196            while 1:
1197                if gid is None:
1198                    break
1199                cparamId = str(grib_get(gid, 'paramId'))
1200                step = grib_get(gid, 'step')
1201                atime = grib_get(gid, 'time')
1202                ni = grib_get(gid, 'Ni')
1203                nj = grib_get(gid, 'Nj')
1204                if cparamId in valsdict.keys():
1205                    values = grib_get_values(gid)
1206                    vdp = valsdict[cparamId]
1207                    svdp = svalsdict[cparamId]
1208                    sd = stepsdict[cparamId]
1209
1210                    if cparamId == '142' or cparamId == '143':
1211                        fak = 1. / 1000.
1212                    else:
1213                        fak = 3600.
1214
1215                    values = (np.reshape(values, (nj, ni))).flatten() / fak
1216                    vdp.append(values[:])  # save the accumulated values
1217                    if step <= int(c.dtime):
1218                        svdp.append(values[:] / int(c.dtime))
1219                    else:  # deaccumulate values
1220                        svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
1221
1222                    print(cparamId, atime, step, len(values),
1223                          values[0], np.std(values))
1224                    # save the 1/3-hourly or specific values
1225                    # svdp.append(values[:])
1226                    sd.append(step)
1227                    # len(svdp) correspond to the time
1228                    if len(svdp) >= 3:
1229                        if len(svdp) > 3:
1230                            if cparamId == '142' or cparamId == '143':
1231                                values = disaggregation.darain(svdp)
1232                            else:
1233                                values = disaggregation.dapoly(svdp)
1234
1235                            if not (step == c.maxstep and c.maxstep > 12 \
1236                                    or sdates == elimit):
1237                                vdp.pop(0)
1238                                svdp.pop(0)
1239                        else:
1240                            if c.maxstep > 12:
1241                                values = svdp[1]
1242                            else:
1243                                values = svdp[0]
1244
1245                        grib_set_values(gid, values)
1246                        if c.maxstep > 12:
1247                            grib_set(gid, 'step', max(0, step-2*int(c.dtime)))
1248                        else:
1249                            grib_set(gid, 'step', 0)
1250                            grib_set(gid, 'time', fdate.hour*100)
1251                            grib_set(gid, 'date', fdate.year*10000 +
1252                                     fdate.month*100+fdate.day)
1253                        grib_write(gid, f_handle)
1254
1255                        if c.basetime:
1256                            elimit = datetime.strptime(c.end_date +
1257                                                       c.basetime, '%Y%m%d%H')
1258                        else:
1259                            elimit = sdate + timedelta(2*int(c.dtime))
1260
1261                        # squeeze out information of last two steps contained
1262                        # in svdp
1263                        # if step+int(c.dtime) == c.maxstep and c.maxstep>12
1264                        # or sdates+timedelta(hours = int(c.dtime))
1265                        # >= elimit:
1266                        # Note that svdp[0] has not been popped in this case
1267
1268                        if step == c.maxstep and c.maxstep > 12 or \
1269                           sdates == elimit:
1270
1271                            values = svdp[3]
1272                            grib_set_values(gid, values)
1273                            grib_set(gid, 'step', 0)
1274                            truedatetime = fdate + timedelta(hours=
1275                                                             2*int(c.dtime))
1276                            grib_set(gid, 'time', truedatetime.hour * 100)
1277                            grib_set(gid, 'date', truedatetime.year * 10000 +
1278                                     truedatetime.month * 100 +
1279                                     truedatetime.day)
1280                            grib_write(gid, h_handle)
1281
1282                            #values = (svdp[1]+svdp[2])/2.
1283                            if cparamId == '142' or cparamId == '143':
1284                                values = disaggregation.darain(list(reversed(svdp)))
1285                            else:
1286                                values = disaggregation.dapoly(list(reversed(svdp)))
1287
1288                            grib_set(gid, 'step', 0)
1289                            truedatetime = fdate + timedelta(hours=int(c.dtime))
1290                            grib_set(gid, 'time', truedatetime.hour * 100)
1291                            grib_set(gid, 'date', truedatetime.year * 10000 +
1292                                     truedatetime.month * 100 +
1293                                     truedatetime.day)
1294                            grib_set_values(gid, values)
1295                            grib_write(gid, g_handle)
1296
1297                    grib_release(gid)
1298
1299                    gid = grib_new_from_index(iid)
1300
1301            f_handle.close()
1302            g_handle.close()
1303            h_handle.close()
1304
1305        grib_index_release(iid)
1306
1307        return
1308
1309
1310
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG