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

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

added making of namelist file and jobscript via genshi templates

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