source: flex_extract.git/source/python/classes/EcFlexpart.py @ aa275fc

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

substituted grib_api with eccodes

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