source: flex_extract.git/source/python/classes/EcFlexpart.py @ 0e576fc

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

implemented extraction possibility of EA5 und CERA

  • 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
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, filename):
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        self.inputdir = c.inputdir
686        area = np.asarray(self.area.split('/')).astype(float)
687        grid = np.asarray(self.grid.split('/')).astype(float)
688
689        if area[1] > area[3]:
690            area[1] -= 360
691        maxl = int((area[3] - area[1]) / grid[1]) + 1
692        maxb = int((area[0] - area[2]) / grid[0]) + 1
693
694        with open(self.inputdir + '/' + filename, 'w') as f:
695            f.write('&NAMGEN\n')
696            f.write(',\n  '.join(['maxl = ' + str(maxl), 'maxb = ' + str(maxb),
697                                  'mlevel = ' + str(self.level),
698                                  'mlevelist = ' + '"' + str(self.levelist)
699                                                 + '"',
700                                  'mnauf = ' + str(self.resol),
701                                  'metapar = ' + '77',
702                                  'rlo0 = ' + str(area[1]),
703                                  'rlo1 = ' + str(area[3]),
704                                  'rla0 = ' + str(area[2]),
705                                  'rla1 = ' + str(area[0]),
706                                  'momega = ' + str(c.omega),
707                                  'momegadiff = ' + str(c.omegadiff),
708                                  'mgauss = ' + str(c.gauss),
709                                  'msmooth = ' + str(c.smooth),
710                                  'meta = ' + str(c.eta),
711                                  'metadiff = ' + str(c.etadiff),
712                                  'mdpdeta = ' + str(c.dpdeta)]))
713
714            f.write('\n/\n')
715
716        return
717
718
719    def deacc_fluxes(self, inputfiles, c):
720        '''
721        @Description:
722            Goes through all flux fields in ordered time and de-accumulate
723            the fields. Afterwards the fields are disaggregated in time.
724            Different versions of disaggregation is provided for rainfall
725            data (darain, modified linear) and the surface fluxes and
726            stress data (dapoly, cubic polynomial).
727
728        @Input:
729            self: instance of EcFlexpart
730                The current object of the class.
731
732            inputfiles: instance of UioFiles
733                Contains a list of files.
734
735            c: instance of class ControlFile
736                Contains all the parameters of CONTROL file and
737                command line.
738                For more information about format and content of the parameter
739                see documentation.
740
741        @Return:
742            <nothing>
743        '''
744
745        table128 = init128(_config.PATH_GRIBTABLE)
746        pars = to_param_id(self.params['OG_acc_SL'][0], table128)
747
748        iid = None
749        index_vals = None
750
751        # get the values of the keys which are used for distinct access
752        # of grib messages via product
753        index_keys = ["date", "time", "step"]
754        iid, index_vals = self._mk_index_values(c.inputdir,
755                                                inputfiles,
756                                                index_keys)
757        # index_vals looks like e.g.:
758        # index_vals[0]: ('20171106', '20171107', '20171108') ; date
759        # index_vals[1]: ('0', '1200', '1800', '600') ; time
760        # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
761
762        valsdict = {}
763        svalsdict = {}
764#        stepsdict = {}
765        for p in pars:
766            valsdict[str(p)] = []
767            svalsdict[str(p)] = []
768#            stepsdict[str(p)] = []
769
770        print('maxstep: ', c.maxstep)
771
772        # "product" genereates each possible combination between the
773        # values of the index keys
774        for prod in product(*index_vals):
775            # e.g. prod = ('20170505', '0', '12')
776            #             (  date    ,time, step)
777
778            print('current product: ', prod)
779
780            for i in range(len(index_keys)):
781                grib_index_select(iid, index_keys[i], prod[i])
782
783            # get first id from current product
784            gid = grib_new_from_index(iid)
785
786            # if there is no data for this specific time combination / product
787            # skip the rest of the for loop and start with next timestep/product
788            if not gid:
789                continue
790
791            # create correct timestamp from the three time informations
792            cdate = str(grib_get(gid, 'date'))
793            ctime = '{:0>2}'.format(grib_get(gid, 'time')/100)
794            cstep = '{:0>3}'.format(grib_get(gid, 'step'))
795            t_date = datetime.strptime(cdate + ctime, '%Y%m%d%H')
796            t_dt = t_date + timedelta(hours=int(cstep))
797            t_m1dt = t_date + timedelta(hours=int(cstep)-int(c.dtime))
798            t_m2dt = t_date + timedelta(hours=int(cstep)-2*int(c.dtime))
799            t_enddate = None
800
801            if c.maxstep > 12:
802                fnout = os.path.join(c.inputdir, 'flux' +
803                                     t_date.strftime('%Y%m%d.%H') +
804                                     '.{:0>3}'.format(step-2*int(c.dtime)))
805                gnout = os.path.join(c.inputdir, 'flux' +
806                                     t_date.strftime('%Y%m%d.%H') +
807                                     '.{:0>3}'.format(step-int(c.dtime)))
808                hnout = os.path.join(c.inputdir, 'flux' +
809                                     t_date.strftime('%Y%m%d.%H') +
810                                     '.{:0>3}'.format(step))
811            else:
812                fnout = os.path.join(c.inputdir, 'flux' +
813                                     t_m2dt.strftime('%Y%m%d%H'))
814                gnout = os.path.join(c.inputdir, 'flux' +
815                                     t_m1dt.strftime('%Y%m%d%H'))
816                hnout = os.path.join(c.inputdir, 'flux' +
817                                     t_dt.strftime('%Y%m%d%H'))
818
819            print("outputfile = " + fnout)
820
821            # read message for message and store relevant data fields
822            # data keywords are stored in pars
823            while 1:
824                if not gid:
825                    break
826                cparamId = str(grib_get(gid, 'paramId'))
827                step = grib_get(gid, 'step')
828                time = grib_get(gid, 'time')
829                ni = grib_get(gid, 'Ni')
830                nj = grib_get(gid, 'Nj')
831                if cparamId in valsdict.keys():
832                    values = grib_get_values(gid)
833                    vdp = valsdict[cparamId]
834                    svdp = svalsdict[cparamId]
835 #                   sd = stepsdict[cparamId]
836
837                    if cparamId == '142' or cparamId == '143':
838                        fak = 1. / 1000.
839                    else:
840                        fak = 3600.
841
842                    values = (np.reshape(values, (nj, ni))).flatten() / fak
843                    vdp.append(values[:])  # save the accumulated values
844                    if c.marsclass.upper() == 'EA' or \
845                       step <= int(c.dtime):
846                        svdp.append(values[:] / int(c.dtime))
847                    else:  # deaccumulate values
848                        svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
849
850                    print(cparamId, time, step, len(values),
851                          values[0], np.std(values))
852
853                    # len(svdp) correspond to the time
854                    if len(svdp) >= 3:
855                        if len(svdp) > 3:
856                            if cparamId == '142' or cparamId == '143':
857                                values = disaggregation.darain(svdp)
858                            else:
859                                values = disaggregation.dapoly(svdp)
860
861                            if not (step == c.maxstep and c.maxstep > 12 \
862                                    or t_dt == t_enddate):
863                                vdp.pop(0)
864                                svdp.pop(0)
865                        else:
866                            if c.maxstep > 12:
867                                values = svdp[1]
868                            else:
869                                values = svdp[0]
870
871                        grib_set_values(gid, values)
872
873                        if c.maxstep > 12:
874                            grib_set(gid, 'step', max(0, step-2*int(c.dtime)))
875                        else:
876                            grib_set(gid, 'step', 0)
877                            grib_set(gid, 'time', t_m2dt.hour*100)
878                            grib_set(gid, 'date', int(t_m2dt.strftime('%Y%m%d')))
879
880                        with open(fnout, 'w') as f_handle:
881                            grib_write(gid, f_handle)
882
883                        if c.basetime:
884                            t_enddate = datetime.strptime(c.end_date +
885                                                          c.basetime,
886                                                          '%Y%m%d%H')
887                        else:
888                            t_enddate = t_date + timedelta(2*int(c.dtime))
889
890                        # squeeze out information of last two steps contained
891                        # in svdp
892                        # if step+int(c.dtime) == c.maxstep and c.maxstep>12
893                        # or t_dt+timedelta(hours = int(c.dtime))
894                        # >= t_enddate:
895                        # Note that svdp[0] has not been popped in this case
896
897                        if step == c.maxstep and c.maxstep > 12 or \
898                           t_dt == t_enddate:
899
900                            values = svdp[3]
901                            grib_set_values(gid, values)
902                            grib_set(gid, 'step', 0)
903                            truedatetime = t_m2dt + timedelta(hours=
904                                                             2*int(c.dtime))
905                            grib_set(gid, 'time', truedatetime.hour * 100)
906                            grib_set(gid, 'date', truedatetime.year * 10000 +
907                                     truedatetime.month * 100 +
908                                     truedatetime.day)
909                            with open(hnout, 'w') as h_handle:
910                                grib_write(gid, h_handle)
911
912                            #values = (svdp[1]+svdp[2])/2.
913                            if cparamId == '142' or cparamId == '143':
914                                values = disaggregation.darain(list(reversed(svdp)))
915                            else:
916                                values = disaggregation.dapoly(list(reversed(svdp)))
917
918                            grib_set(gid, 'step', 0)
919                            truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
920                            grib_set(gid, 'time', truedatetime.hour * 100)
921                            grib_set(gid, 'date', truedatetime.year * 10000 +
922                                     truedatetime.month * 100 +
923                                     truedatetime.day)
924                            grib_set_values(gid, values)
925                            with open(gnout, 'w') as g_handle:
926                                grib_write(gid, g_handle)
927
928                grib_release(gid)
929
930                gid = grib_new_from_index(iid)
931
932        grib_index_release(iid)
933
934        return
935
936
937    def create(self, inputfiles, c):
938        '''
939        @Description:
940            This method is based on the ECMWF example index.py
941            https://software.ecmwf.int/wiki/display/GRIB/index.py
942
943            An index file will be created which depends on the combination
944            of "date", "time" and "stepRange" values. This is used to iterate
945            over all messages in each grib file which were passed through the
946            parameter "inputfiles" to seperate specific parameters into fort.*
947            files. Afterwards the FORTRAN program is called to convert
948            the data fields all to the same grid and put them in one file
949            per unique time step (combination of "date", "time" and
950            "stepRange").
951
952        @Input:
953            self: instance of EcFlexpart
954                The current object of the class.
955
956            inputfiles: instance of UioFiles
957                Contains a list of files.
958
959            c: instance of class ControlFile
960                Contains all the parameters of CONTROL file and
961                command line.
962                For more information about format and content of the parameter
963                see documentation.
964
965        @Return:
966            <nothing>
967        '''
968
969        if c.wrf:
970            table128 = init128(_config.PATH_GRIBTABLE)
971            wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
972                                   stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
973                                  table128)
974
975        # these numbers are indices for the temporary files "fort.xx"
976        # which are used to seperate the grib fields to,
977        # for the Fortran program input
978        # 10: U,V | 11: T | 12: lnsp | 13: D | 16: sfc fields
979        # 17: Q | 18: Q , gaussian| 19: w | 21: etadot | 22: clwc+ciwc
980        fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
981                 '17':None, '18':None, '19':None, '21':None, '22':None}
982
983        iid = None
984        index_vals = None
985
986        # get the values of the keys which are used for distinct access
987        # of grib messages via product
988        index_keys = ["date", "time", "step"]
989        iid, index_vals = self._mk_index_values(c.inputdir,
990                                                inputfiles,
991                                                index_keys)
992        # index_vals looks like e.g.:
993        # index_vals[0]: ('20171106', '20171107', '20171108') ; date
994        # index_vals[1]: ('0', '1200', '1800', '600') ; time
995        # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
996
997        # "product" genereates each possible combination between the
998        # values of the index keys
999        for prod in product(*index_vals):
1000            # e.g. prod = ('20170505', '0', '12')
1001            #             (  date    ,time, step)
1002
1003            print('current product: ', prod)
1004
1005            for i in range(len(index_keys)):
1006                grib_index_select(iid, index_keys[i], prod[i])
1007
1008            # get first id from current product
1009            gid = grib_new_from_index(iid)
1010
1011            # if there is no data for this specific time combination / product
1012            # skip the rest of the for loop and start with next timestep/product
1013            if not gid:
1014                continue
1015
1016            # remove old fort.* files and open new ones
1017            # they are just valid for a single product
1018            for k, f in fdict.iteritems():
1019                fortfile = os.path.join(c.inputdir, 'fort.' + k)
1020                silent_remove(fortfile)
1021                fdict[k] = open(fortfile, 'w')
1022
1023            # create correct timestamp from the three time informations
1024            cdate = str(grib_get(gid, 'date'))
1025            ctime = '{:0>2}'.format(grib_get(gid, 'time')/100)
1026            cstep = '{:0>3}'.format(grib_get(gid, 'step'))
1027            timestamp = datetime.strptime(cdate + ctime, '%Y%m%d%H')
1028            timestamp += timedelta(hours=int(cstep))
1029            cdate_hour = datetime.strftime(timestamp, '%Y%m%d%H')
1030
1031            # if the timestamp is out of basetime start/end date period,
1032            # skip this specific product
1033            if c.basetime:
1034                start_time = datetime.strptime(c.end_date + c.basetime,
1035                                                '%Y%m%d%H') - time_delta
1036                end_time = datetime.strptime(c.end_date + c.basetime,
1037                                             '%Y%m%d%H')
1038                if timestamp < start_time or timestamp > end_time:
1039                    continue
1040
1041            if c.wrf:
1042                if 'olddate' not in locals() or cdate != olddate:
1043                    fwrf = open(os.path.join(c.outputdir,
1044                                'WRF' + cdate + '.' + ctime + '.000.grb2'), 'w')
1045                    olddate = cdate[:]
1046
1047            # savedfields remembers which fields were already used.
1048            savedfields = []
1049            # sum of cloud liquid and ice water content
1050            scwc = None
1051            while 1:
1052                if not gid:
1053                    break
1054                paramId = grib_get(gid, 'paramId')
1055                gridtype = grib_get(gid, 'gridType')
1056                levtype = grib_get(gid, 'typeOfLevel')
1057                if paramId == 77: # ETADOT
1058                    grib_write(gid, fdict['21'])
1059                elif paramId == 130: # T
1060                    grib_write(gid, fdict['11'])
1061                elif paramId == 131 or paramId == 132: # U, V wind component
1062                    grib_write(gid, fdict['10'])
1063                elif paramId == 133 and gridtype != 'reduced_gg': # Q
1064                    grib_write(gid, fdict['17'])
1065                elif paramId == 133 and gridtype == 'reduced_gg': # Q, gaussian
1066                    grib_write(gid, fdict['18'])
1067                elif paramId == 135: # W
1068                    grib_write(gid, fdict['19'])
1069                elif paramId == 152: # LNSP
1070                    grib_write(gid, fdict['12'])
1071                elif paramId == 155 and gridtype == 'sh': # D
1072                    grib_write(gid, fdict['13'])
1073                elif paramId == 246 or paramId == 247: # CLWC, CIWC
1074                    # sum cloud liquid water and ice
1075                    if not scwc:
1076                        scwc = grib_get_values(gid)
1077                    else:
1078                        scwc += grib_get_values(gid)
1079                        grib_set_values(gid, scwc)
1080                        grib_set(gid, 'paramId', 201031)
1081                        grib_write(gid, fdict['22'])
1082                elif c.wrf and paramId in [129, 138, 155] and \
1083                      levtype == 'hybrid': # Z, VO, D
1084                    # do not do anything right now
1085                    # these are specific parameter for WRF
1086                    pass
1087                else:
1088                    if paramId not in savedfields:
1089                        # SD/MSL/TCC/10U/10V/2T/2D/Z/LSM/SDOR/CVL/CVH/SR
1090                        # and all ADDPAR parameter
1091                        grib_write(gid, fdict['16'])
1092                        savedfields.append(paramId)
1093                    else:
1094                        print('duplicate ' + str(paramId) + ' not written')
1095
1096                try:
1097                    if c.wrf:
1098                        # model layer
1099                        if levtype == 'hybrid' and \
1100                           paramId in [129, 130, 131, 132, 133, 138, 155]:
1101                            grib_write(gid, fwrf)
1102                        # sfc layer
1103                        elif paramId in wrfpars:
1104                            grib_write(gid, fwrf)
1105                except AttributeError:
1106                    pass
1107
1108                grib_release(gid)
1109                gid = grib_new_from_index(iid)
1110
1111            for f in fdict.values():
1112                f.close()
1113
1114            # call for Fortran program to convert e.g. reduced_gg grids to
1115            # regular_ll and calculate detadot/dp
1116            pwd = os.getcwd()
1117            os.chdir(c.inputdir)
1118            if os.stat('fort.21').st_size == 0 and c.eta:
1119                print('Parameter 77 (etadot) is missing, most likely it is \
1120                       not available for this type or date/time\n')
1121                print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
1122                my_error(c.mailfail, 'fort.21 is empty while parameter eta \
1123                         is set to 1 in CONTROL file')
1124
1125            # Fortran program creates file fort.15 (with u,v,etadot,t,sp,q)
1126            p = subprocess.check_call([os.path.join(
1127                c.exedir, _config.FORTRAN_EXECUTABLE)], shell=True)
1128            os.chdir(pwd)
1129
1130            # create name of final output file, e.g. EN13040500 (ENYYMMDDHH)
1131            if c.maxstep > 12:
1132                suffix = cdate[2:8] + '.' + ctime + '.' + cstep
1133            else:
1134                suffix = cdate_hour[2:10]
1135            fnout = os.path.join(c.inputdir, c.prefix + suffix)
1136            print("outputfile = " + fnout)
1137            # collect for final processing
1138            self.outputfilelist.append(os.path.basename(fnout))
1139
1140            # create outputfile and copy all data from intermediate files
1141            # to the outputfile (final GRIB input files for FLEXPART)
1142            orolsm = os.path.basename(glob.glob(c.inputdir +
1143                                        '/OG_OROLSM__SL.*.' + c.ppid + '*')[0])
1144            fluxfile = 'flux' + cdate[0:2] + suffix
1145            if not c.cwc:
1146                flist = ['fort.15', fluxfile, 'fort.16', orolsm]
1147            else:
1148                flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]
1149
1150            with open(fnout, 'wb') as fout:
1151                for f in flist:
1152                    shutil.copyfileobj(open(os.path.join(c.inputdir, f), 'rb'),
1153                                       fout)
1154
1155            if c.omega:
1156                with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout:
1157                    shutil.copyfileobj(open(os.path.join(c.inputdir, 'fort.25'),
1158                                            'rb'), fout)
1159
1160        if c.wrf:
1161            fwrf.close()
1162
1163        grib_index_release(iid)
1164
1165        return
1166
1167
1168    def process_output(self, c):
1169        '''
1170        @Description:
1171            The grib files are postprocessed depending on the selection in
1172            CONTROL file. The resulting files are moved to the output
1173            directory if its not equal to the input directory.
1174            The following modifications might be done if
1175            properly switched in CONTROL file:
1176            GRIB2 - Conversion to GRIB2
1177            ECTRANS - Transfer of files to gateway server
1178            ECSTORAGE - Storage at ECMWF server
1179
1180        @Input:
1181            self: instance of EcFlexpart
1182                The current object of the class.
1183
1184            c: instance of class ControlFile
1185                Contains all the parameters of CONTROL file and
1186                command line.
1187                For more information about format and content of the parameter
1188                see documentation.
1189
1190        @Return:
1191            <nothing>
1192
1193        '''
1194
1195        print('\n\nPostprocessing:\n Format: {}\n'.format(c.format))
1196
1197        if not c.ecapi:
1198            print('ecstorage: {}\n ecfsdir: {}\n'.
1199                  format(c.ecstorage, c.ecfsdir))
1200            print('ectrans: {}\n gateway: {}\n destination: {}\n '
1201                  .format(c.ectrans, c.gateway, c.destination))
1202
1203        print('Output filelist: ')
1204        print(self.outputfilelist)
1205
1206        if c.format.lower() == 'grib2':
1207            for ofile in self.outputfilelist:
1208                p = subprocess.check_call(['grib_set', '-s', 'edition=2, \
1209                                           productDefinitionTemplateNumber=8',
1210                                           ofile, ofile + '_2'])
1211                p = subprocess.check_call(['mv', ofile + '_2', ofile])
1212
1213        if c.ectrans and not c.ecapi:
1214            for ofile in self.outputfilelist:
1215                p = subprocess.check_call(['ectrans', '-overwrite', '-gateway',
1216                                           c.gateway, '-remote', c.destination,
1217                                           '-source', ofile])
1218
1219        if c.ecstorage and not c.ecapi:
1220            for ofile in self.outputfilelist:
1221                p = subprocess.check_call(['ecp', '-o', ofile,
1222                                           os.path.expandvars(c.ecfsdir)])
1223
1224        if c.outputdir != c.inputdir:
1225            for ofile in self.outputfilelist:
1226                p = subprocess.check_call(['mv',
1227                                           os.path.join(c.inputdir, ofile),
1228                                           c.outputdir])
1229
1230        return
1231
1232
1233    def prepare_fp_files(self, c):
1234        '''
1235        @Description:
1236            Conversion of GRIB files to FLEXPART binary format.
1237
1238        @Input:
1239            c: instance of class ControlFile
1240                Contains all the parameters of CONTROL file and
1241                command line.
1242                For more information about format and content of the parameter
1243                see documentation.
1244
1245
1246        @Return:
1247            <nothing>
1248        '''
1249        # generate AVAILABLE file
1250        # Example of AVAILABLE file data:
1251        # 20131107 000000      EN13110700              ON DISC
1252        clist = []
1253        for ofile in self.outputfilelist:
1254            fname = ofile.split('/')
1255            if '.' in fname[-1]:
1256                l = fname[-1].split('.')
1257                timestamp = datetime.strptime(l[0][-6:] + l[1],
1258                                              '%y%m%d%H')
1259                timestamp += timedelta(hours=int(l[2]))
1260                cdate = datetime.strftime(timestamp, '%Y%m%d')
1261                chms = datetime.strftime(timestamp, '%H%M%S')
1262            else:
1263                cdate = '20' + fname[-1][-8:-2]
1264                chms = fname[-1][-2:] + '0000'
1265            clist.append(cdate + ' ' + chms + ' '*6 +
1266                         fname[-1] + ' '*14 + 'ON DISC')
1267        clist.sort()
1268        with open(c.outputdir + '/' + 'AVAILABLE', 'w') as f:
1269            f.write('\n'.join(clist) + '\n')
1270
1271        # generate pathnames file
1272        pwd = os.path.abspath(c.outputdir)
1273        with open(pwd + '/pathnames', 'w') as f:
1274            f.write(pwd + '/Options/\n')
1275            f.write(pwd + '/\n')
1276            f.write(pwd + '/\n')
1277            f.write(pwd + '/AVAILABLE\n')
1278            f.write(' = == = == = == = == = == ==  = \n')
1279
1280        # create Options dir if necessary
1281        if not os.path.exists(pwd + '/Options'):
1282            make_dir(pwd+'/Options')
1283
1284        # read template COMMAND file
1285        with open(os.path.expandvars(os.path.expanduser(
1286            c.flexpart_root_scripts)) + '/../Options/COMMAND', 'r') as f:
1287            lflist = f.read().split('\n')
1288
1289        # find index of list where to put in the
1290        # date and time information
1291        # usually after the LDIRECT parameter
1292        i = 0
1293        for l in lflist:
1294            if 'LDIRECT' in l.upper():
1295                break
1296            i += 1
1297
1298        # insert the date and time information of run start and end
1299        # into the list of lines of COMMAND file
1300        lflist = lflist[:i+1] + \
1301                 [clist[0][:16], clist[-1][:16]] + \
1302                 lflist[i+3:]
1303
1304        # write the new COMMAND file
1305        with open(pwd + '/Options/COMMAND', 'w') as g:
1306            g.write('\n'.join(lflist) + '\n')
1307
1308        # change to outputdir and start the grib2flexpart run
1309        # afterwards switch back to the working dir
1310        os.chdir(c.outputdir)
1311        p = subprocess.check_call([
1312            os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))
1313            + '/../FLEXPART_PROGRAM/grib2flexpart', 'useAvailable', '.'])
1314        os.chdir(pwd)
1315
1316        return
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG