source: flex_extract.git/python/FlexpartTools.py @ 51f9853

Last change on this file since 51f9853 was 51f9853, checked in by anphi <anne.philipp@…>, 8 months ago

added request parameter for writing mars requests into csv file

  • Property mode set to 100644
File size: 66.6 KB
Line 
1#
2# This software is licensed under the terms of the Apache Licence Version 2.0
3# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
4#
5# Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs
6#
7# Creation: October  2014 - Anne Fouilloux - University of Oslo
8# Extension November 2015 - Leopold Haimberger - University of Vienna for:
9# - using the WebAPI also for general MARS retrievals
10# - job submission on ecgate and cca
11# - job templates suitable for twice daily operational dissemination
12# - dividing retrievals of longer periods into digestable chunks
13# - retrieve also longer term forecasts, not only analyses and short term forecast data
14# - conversion into GRIB2
15# - conversion into .fp format for faster execution of FLEXPART
16#
17# Extension May 2018 - Anne Philipp - University of Vienna
18# - introduced new CONTROL parameter to seperate accumulation and other
19#   data specification
20# - modified code to be able to retrieve ERA-5 data
21# - modified code to be able to retrieve CERA-20C data
22# - added possibility to retrieve data from public data server als as
23#   public user
24# - allow specifying TIME/TYPE/STEP combinations only with needed data instead
25#   of the full 24 hours
26#
27# Requirements:
28# in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed
29# ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/
30# dateutils
31# matplotlib (optional, for debugging)
32#
33#
34import subprocess
35from argparse import ArgumentParser,ArgumentDefaultsHelpFormatter
36import traceback
37import shutil
38import os, errno,sys,inspect,glob
39import datetime
40
41from string import atoi
42from numpy  import *
43
44ecapi=True
45try:
46    import ecmwfapi
47except ImportError:
48    ecapi=False
49from gribapi import *
50from GribTools import GribTools
51
52def interpret_args_and_control(*args,**kwargs):
53
54    parser = ArgumentParser(description='Retrieve FLEXPART input from ECMWF MARS archive',
55                            formatter_class=ArgumentDefaultsHelpFormatter)
56
57# the most important arguments
58    parser.add_argument("--start_date", dest="start_date",
59                        help="start date YYYYMMDD")
60    parser.add_argument( "--end_date", dest="end_date",
61                         help="end_date YYYYMMDD")
62    parser.add_argument( "--date_chunk", dest="date_chunk",default=None,
63                         help="# of days to be retrieved at once")
64
65# some arguments that override the default in the control file
66    parser.add_argument("--basetime", dest="basetime", help="base such as 00/12 (for half day retrievals)")
67    parser.add_argument("--step", dest="step", help="steps such as 00/to/48")
68    parser.add_argument("--levelist", dest="levelist",help="Vertical levels to be retrieved, e.g. 30/to/60")
69    parser.add_argument("--area", dest="area", help="area defined as north/west/south/east")
70
71# set the working directories
72    parser.add_argument("--inputdir", dest="inputdir",default=None,
73                        help="root directory for storing intermediate files")
74    parser.add_argument("--outputdir", dest="outputdir",default=None,
75                        help="root directory for storing output files")
76    parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
77                        help="FLEXPART root directory (to find grib2flexpart and COMMAND file)\n\
78                              Normally flex_extract resides in the scripts directory of the FLEXPART distribution")
79
80# this is only used by prepareFLEXPART.py to rerun a postprocessing step
81    parser.add_argument("--ppid", dest="ppid",
82                        help="Specify parent process id for rerun of prepareFLEXPART")
83
84# arguments for job submission to ECMWF, only needed by submit.py
85    parser.add_argument("--job_template", dest='job_template',default="job.temp",
86                        help="job template file for submission to ECMWF")
87    #parser.add_argument("--remote", dest="remote",
88                        #help="target for submission to ECMWF (ecgate or cca etc.)")
89    parser.add_argument("--queue", dest="queue",
90                        help="queue for submission to ECMWF (e.g. ecgate or cca )")
91
92
93    parser.add_argument("--controlfile", dest="controlfile",default='CONTROL.temp',
94                        help="file with control parameters")
95    parser.add_argument("--debug", dest="debug",default=0,
96                        help="Debug mode - leave temporary files intact")
97    parser.add_argument("--public", dest="public", default=0,
98                        help="Public mode - retrieves the public datasets")
99
100    parser.add_argument("--request", dest="request", default=0,
101                        help="list all mars request in file mars_requests.dat \
102                        and skip submission to mars")
103
104    args = parser.parse_args()
105
106    try:
107        c=Control(args.controlfile)
108    except IOError:
109        try:
110            c=Control(localpythonpath+args.controlfile)
111
112        except:
113            print 'Could not read control file "'+args.controlfile+'"'
114            print 'Either it does not exist or its syntax is wrong.'
115            print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information'
116            exit(1)
117
118    if  args.start_date==None and getattr(c,'start_date')==None:
119        print 'start_date specified neither in command line nor in control file '+args.controlfile
120        print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information'
121        exit(1)
122
123    if args.start_date!=None:
124        c.start_date=args.start_date
125    if args.end_date!=None:
126        c.end_date=args.end_date
127    if c.end_date==None:
128        c.end_date=c.start_date
129    if args.date_chunk!=None:
130        c.date_chunk=args.date_chunk
131
132    if not hasattr(c,'debug'):
133        c.debug=args.debug
134
135    if not hasattr(c, 'public'):
136        setattr(c, 'public', args.public)
137
138    if args.inputdir==None and args.outputdir==None:
139        c.inputdir='../work'
140        c.outputdir='../work'
141    else:
142        if args.inputdir!=None:
143            c.inputdir=args.inputdir
144            if args.outputdir==None:
145                c.outputdir=args.inputdir
146        if args.outputdir!=None:
147            c.outputdir=args.outputdir
148            if args.inputdir==None:
149                c.inputdir=args.outputdir
150
151    if hasattr(c,'outputdir')==False and args.outputdir==None:
152        c.outputdir=c.inputdir
153    else:
154        if args.outputdir!=None:
155            c.outputdir=args.outputdir
156
157    if args.area!=None:
158        afloat='.' in args.area
159        l=args.area.split('/')
160        if afloat:
161            for i in range(len(l)):
162                l[i]=str(int(float(l[i])*1000))
163        c.upper,c.left,c.lower,c.right=l
164
165# basetime aktiviert den ''operational mode''
166    if args.basetime!=None:
167        c.basetime=args.basetime
168    #if int(c.basetime)==0:
169        #c.start_date=datetime.datetime.strftime(datetime.datetime.strptime(c.start_date,'%Y%m%d')-datetime.timedelta(days=1),'%Y%m%d')
170        #c.end_date=c.start_date
171
172    if args.step!=None:
173        l=args.step.split('/')
174        if 'to' in args.step.lower():
175            if 'by' in args.step.lower():
176                ilist=arange(int(l[0]),int(l[2])+1,int(l[4]))
177                c.step=['{:0>3}'.format(i) for i in ilist]
178            else:
179                myerror(None,args.step+':\n'+'please use "by" as well if "to" is used')
180        else:
181            c.step=l
182
183    if args.levelist!=None:
184        c.levelist=args.levelist
185        if 'to' in c.levelist:
186            c.level=c.levelist.split('/')[2]
187        else:
188            c.level=c.levelist.split('/')[-1]
189
190
191    if args.flexpart_root_scripts!=None:
192        c.flexpart_root_scripts=args.flexpart_root_scripts
193
194    # set request attribute to control file
195    if args.request != '0':
196        c.request=args.request
197
198    if c.request != '0':
199        marsfile = os.path.join(c.inputdir, 'mars_request.csv')
200        if os.path.isfile(marsfile):
201            os.remove(marsfile)
202
203    return args,c
204
205def install_args_and_control():
206
207    parser = ArgumentParser(description='Install flex_extract software locally or on ECMWF machines',
208                            formatter_class=ArgumentDefaultsHelpFormatter)
209    parser.add_argument('--target',dest='install_target',
210                        help="Valid targets: local | ecgate | cca , the latter two are at ECMWF")
211    parser.add_argument("--makefile", dest="makefile",help='Name of Makefile to use for compiling CONVERT2')
212    parser.add_argument("--ecuid", dest="ecuid",help='user id at ECMWF')
213    parser.add_argument("--ecgid", dest="ecgid",help='group id at ECMWF')
214    parser.add_argument("--gateway", dest="gateway",help='name of local gateway server')
215    parser.add_argument("--destination", dest="destination",help='ecaccess destination, e.g. leo@genericSftp')
216
217    parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
218                        help="FLEXPART root directory on ECMWF servers (to find grib2flexpart and COMMAND file)\n\
219                              Normally flex_extract resides in the scripts directory of the FLEXPART distribution, thus the:")
220
221# arguments for job submission to ECMWF, only needed by submit.py
222    parser.add_argument("--job_template", dest='job_template',default="job.temp.o",
223                        help="job template file for submission to ECMWF")
224
225    parser.add_argument("--controlfile", dest="controlfile",default='CONTROL.temp',
226                        help="file with control parameters")
227    args = parser.parse_args()
228
229    try:
230#    if True:
231        c=Control(args.controlfile)
232#    else:
233    except:
234        print 'Could not read control file "'+args.controlfile+'"'
235        print 'Either it does not exist or its syntax is wrong.'
236        print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information'
237        exit(1)
238
239    if args.install_target!='local':
240        if args.ecgid==None or args.ecuid==None or args.gateway==None or args.destination==None :
241            print 'Please enter your ECMWF user id and group id as well as the \nname of the local gateway and the ectrans destination '
242            print 'with command line options --ecuid --ecgid --gateway --destination'
243            print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information'
244            print 'Please consult ecaccess documentation or ECMWF user support for further details'
245            sys.exit(1)
246        else:
247            c.ecuid=args.ecuid
248            c.ecgid=args.ecgid
249            c.gateway=args.gateway
250            c.destination=args.destination
251
252    try:
253        c.makefile=args.makefile
254    except:
255        pass
256
257    if args.install_target=='local':
258        if args.flexpart_root_scripts==None:
259            c.flexpart_root_scripts='../'
260        else:
261            c.flexpart_root_scripts=args.flexpart_root_scripts
262
263    if args.install_target!='local':
264        if args.flexpart_root_scripts==None:
265            c.ec_flexpart_root_scripts='${HOME}'
266        else:
267            c.ec_flexpart_root_scripts=args.flexpart_root_scripts
268
269    return args,c
270
271def cleanup(c):
272    print "cleanup"
273    cleanlist=glob.glob(c.inputdir+"/*")
274#    cleanlist[-1:-1]=glob.glob(c.inputdir+"/flux*")
275#    cleanlist[-1:-1]=glob.glob(c.inputdir+"/*.grb")
276#    if c.ecapi==False and (c.ectrans=='1' or c.ecstorage=='1'):
277#       cleanlist[-1:-1]=glob.glob(c.inputdir+"/"+c.prefix+"*")
278
279
280#    silentremove(c.inputdir+"/VERTICAL.EC")
281#    silentremove(c.inputdir+"/date_time_stepRange.idx")
282    for cl in cleanlist:
283        if c.prefix not in cl:
284            silentremove(cl)
285        if c.ecapi==False and (c.ectrans=='1' or c.ecstorage=='1'):
286            silentremove(cl)
287
288
289
290    print "Done"
291
292
293def myerror(c,message='ERROR'):
294
295    try:
296        target=c.mailfail
297    except AttributeError:
298        target=os.getenv('USER')
299
300    if(type(target) is not list):
301        target=[target]
302    print message
303    # uncomment if user wants email notification directly from python
304    #for t in target:
305        #p=subprocess.Popen(['mail','-s flex_extract v7.0.3 ERROR', os.path.expandvars(t)], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,bufsize=1)
306        #tr='\n'.join(traceback.format_stack())
307        #pout=p.communicate(input=message+'\n\n'+tr)[0]
308        #print 'Email sent to '+os.path.expandvars(t) # +' '+pout.decode()
309
310    exit(1)
311
312def normalexit(c,message='Done!'):
313
314    try:
315        target=c.mailops
316
317
318        if(type(target) is not list):
319            target=[target]
320        # Uncomment if user wants notification directly from python
321        #for t in target:
322            #p=subprocess.Popen(['mail','-s flex_extract v7.0.3 normal exit', os.path.expandvars(t)], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,bufsize=1)
323            #pout=p.communicate(input=message+'\n\n')[0]
324            #print pout.decode()
325
326    except:
327        pass
328
329    print message
330
331    return
332
333def product(*args, **kwds):
334    # product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
335    # product(range(2), repeat=3) --> 000 001 010 011 100 101 110 111
336    pools = map(tuple, args) * kwds.get('repeat', 1)
337    result = [[]]
338    for pool in pools:
339        result = [x+[y] for x in result for y in pool]
340    for prod in result:
341        yield tuple(prod)
342
343###############################################################
344# utility to remove a file if it exists
345# it does not fail if the file does not exist
346###############################################################
347def silentremove(filename):
348    try:
349        os.remove(filename)
350    except OSError as e: # this would be "except OSError, e:" before Python 2.6
351        if e.errno != errno.ENOENT: # errno.ENOENT = no such file or directory
352            raise # re-raise exception if a different error occured
353
354
355def init128(fn):
356
357    table128=dict()
358    with open(fn) as f:
359        fdata = f.read().split('\n')
360        for data in fdata:
361            if data[0]!='!':
362                table128[data[0:3]]=data[59:64].strip()
363
364    return table128
365
366def toparamId(pars,table):
367
368    cpar=pars.upper().split('/')
369    ipar=[]
370    for par in cpar:
371        found=False
372        for k,v in table.iteritems():
373            if par==k or par==v:
374                ipar.append(int(k))
375                found=True
376                break
377        if found==False:
378            print 'Warning: par '+par+' not found in table 128'
379    return ipar
380
381def dapoly(alist):
382
383
384    pya=(alist[3]-alist[0]+3.*(alist[1]-alist[2]))/6.
385    pyb=(alist[2]+alist[0])/2.-alist[1]-9.*pya/2.
386    pyc=alist[1]-alist[0]-7.*pya/2.-2.*pyb
387    pyd=alist[0]-pya/4.-pyb/3.-pyc/2.
388    sfeld=8.*pya+4.*pyb+2.*pyc+pyd
389
390    return sfeld
391
392def darain(alist):
393    xa=alist[0]
394    xb=alist[1]
395    xc=alist[2]
396    xd=alist[3]
397    xa[xa<0.]=0.
398    xb[xb<0.]=0.
399    xc[xc<0.]=0.
400    xd[xd<0.]=0.
401
402    xac=0.5*xb
403    mask=xa+xc>0.
404    xac[mask]=xb[mask]*xc[mask]/(xa[mask]+xc[mask])
405    xbd=0.5*xc
406    mask=xb+xd>0.
407    xbd[mask]=xb[mask]*xc[mask]/(xb[mask]+xd[mask])
408    sfeld=xac+xbd
409
410    return sfeld
411
412class Control:
413    'class containing the information of the flex_extract control file'
414
415    def __init__(self,filename):
416
417        #try:
418        with open(filename) as f:
419            fdata = f.read().split('\n')
420            for ldata in fdata:
421                data=ldata.split()
422                if len(data)>1:
423                    if 'm_' in data[0].lower():
424                        data[0]=data[0][2:]
425                    if data[0].lower()=='class':
426                        data[0]='marsclass'
427                    if data[0].lower()=='day1':
428                        data[0]='start_date'
429                    if data[0].lower()=='day2':
430                        data[0]='end_date'
431                    if data[0].lower()=='addpar':
432                        if '/' in data[1]:
433                            if data[1][0]=='/':
434                                data[1]=data[1][1:]
435                            dd=data[1].split('/')
436                            data=[data[0]]
437                            for d in dd:
438                                data.append(d)
439                        pass
440                    if len(data)==2:
441                        if '$' in data[1]:
442#                               print data[1]
443                            setattr(self,data[0].lower(),data[1])
444                            while '$' in data[1]:
445                                i=data[1].index('$')
446                                j=data[1].find('{')
447                                k=data[1].find('}')
448                                var=os.getenv(data[1][j+1:k])
449                                if var is not None:
450                                    data[1]=data[1][:i]+var+data[1][k+1:]
451                                else:
452                                    myerror(None,'Could not find variable '+data[1][j+1:k]+' while reading '+filename)
453
454#                                   print data[0],data[1]
455                            setattr(self,data[0].lower()+'_expanded',data[1])
456                        else:
457                            if data[1].lower()!='none':
458                                setattr(self,data[0].lower(),data[1])
459                            else:
460                                setattr(self,data[0].lower(),None)
461
462                    elif len(data)>2:
463                        setattr(self,data[0].lower(),(data[1:]))
464                else:
465                    pass
466
467            if not hasattr(self,'start_date'):
468                self.start_date=None
469            if not hasattr(self,'end_date'):
470                self.end_date=self.start_date
471            if not hasattr(self,'accuracy'):
472                self.accuracy=24
473            if not hasattr(self,'omega'):
474                self.omega='0'
475            if not hasattr(self,'cwc'):
476                self.cwc='0'
477            if not hasattr(self,'omegadiff'):
478                self.omegadiff='0'
479            if not hasattr(self,'etadiff'):
480                self.etadiff='0'
481
482            if not isinstance(self.type, list):
483                self.type = [self.type]
484            if not isinstance(self.time, list):
485                self.time = [self.time]
486            if not isinstance(self.step, list):
487                self.step = [self.step]
488
489            max_level_list = [16, 19, 31, 40, 50, 60, 62, 91, 137]
490            if not hasattr(self,'levelist'):
491                if not hasattr(self,'level'):
492                    print 'WARNING: neither levelist nor level specified in CONTROL file'
493                else:
494                    if self.level in max_level_list:
495                        self.levelist='1/to/'+self.level
496            else:
497                if not hasattr(self, 'level'):
498                    if 'to' in self.levelist.lower():
499                        max_level = self.levelist.split('/')[2]
500                        if int(max_level) in max_level_list:
501                            self.level = max_level
502                        else:
503                            print 'ERROR: LEVEL-parameter could not be set'
504                            print 'LEVEL must be the maximum level of a specified level list from ECMWF, e.g.'
505                            print '16, 19, 31, 40, 50, 60, 62, 91 or 137'
506                            sys.exit(1)
507                    else:
508                        max_level = self.levelist.split('/')[-1]
509                        if int(max_level) in max_level_list:
510                            self.level = max_level
511                        else:
512                            print 'ERROR: LEVEL-parameter could not be set'
513                            print 'LEVEL must be the maximum level of a specified level list from ECMWF, e.g.'
514                            print '16, 19, 31, 40, 50, 60, 62, 91 or 137'
515                            sys.exit(1)
516                else:
517                    # check if consistent
518                    if int(self.level) in max_level_list:
519                        pass
520                    else:
521                        print 'ERROR: LEVEL-parameter is not properly selected'
522                        print 'LEVEL must be the maximum level of a specified level list from ECMWF, e.g.'
523                        print '16, 19, 31, 40, 50, 60, 62, 91 or 137'
524                        sys.exit(1)
525
526            if not hasattr(self,'maxstep'):
527                self.maxstep=0
528                for s in self.step:
529                    if int(s)>self.maxstep:
530                        self.maxstep=int(s)
531            else:
532                self.maxstep=int(self.maxstep)
533            if not hasattr(self,'prefix'):
534                self.prefix='EN'
535            if not hasattr(self,'makefile'):
536                self.makefile=None
537            if not hasattr(self,'basetime'):
538                self.basetime=None
539            if not hasattr(self,'date_chunk'):
540                self.date_chunk='3'
541
542            self.flexextractdir=os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))+'/../' # script directory
543#               if not hasattr(self,'exedir'):
544            self.exedir=self.flexextractdir+'src/'
545            if not hasattr(self,'grib2flexpart'):
546                self.grib2flexpart='0'
547            if not hasattr(self,'flexpart_root_scripts'):
548                self.flexpart_root_scripts=self.flexextractdir
549
550            if not hasattr(self,'acctype'):
551                print('... Control paramter ACCTYPE was not defined in CONTROL file.')
552                print('Important for (forecast) flux data extraction.')
553                try:
554                    if len(self.type) > 1 and self.type[1] != 'AN':
555                        print('Use old setting by using TYPE[1].')
556                        self.acctype = self.type[1]
557                except:
558                    print('Use default value "FC"!')
559                    self.acctype='FC'
560
561            if not hasattr(self,'acctime'):
562                print('... Control paramter ACCTIME was not defined in CONTROL file.')
563                print('Use default value "00/12" for flux forecast!')
564                print('This may be different for single data sets, see documentation if retrieval fails!')
565                self.acctime='00/12'
566
567            if not hasattr(self, 'accmaxstep'):
568                print('... Control paramter ACCMAXSTEP was not defined in CONTROL file.')
569                print('Use default value "12" for flux forecast!')
570                print('This may be different for single data sets, see documentation if retrieval fails!')
571                self.accmaxstep='12'
572
573            for i in range(len(self.type)):
574                if self.type[i] == 'AN' and int(self.step[i]) != 0:
575                    print('Analysis retrievals must have STEP = 0 (setting to 0)')
576                    self.type[i] = 0
577
578            if not hasattr(self,'request'):
579                self.request='0'
580
581        return
582    def __str__(self):
583
584        attrs = vars(self)
585        # {'kids': 0, 'name': 'Dog', 'color': 'Spotted', 'age': 10, 'legs': 2, 'smell': 'Alot'}
586        # now dump this in some way or another
587        return ', '.join("%s: %s" % item for item in attrs.items())
588
589    def tolist(self):
590        attrs=vars(self)
591        l=list()
592        for item in attrs.items():
593            if '_expanded' in item[0]:
594                pass
595            elif 'exedir' in item[0]:
596                pass
597            elif 'flexpart_root_scripts' in item[0]:
598                pass
599            elif 'flexextractdir' in item[0]:
600                pass
601
602            else:
603                if type(item[1]) is list:
604                    stot=''
605                    for s in item[1]:
606                        stot+=s+' '
607
608                    l.append("%s %s" % (item[0],stot))
609                else:
610                    l.append("%s %s" % item )
611        return sorted(l)
612
613
614
615##############################################################
616# MARSretrieval class
617##############################################################
618class MARSretrieval:
619    'class for MARS retrievals'
620
621    def __init__(self,server, public, marsclass="ei",type="",levtype="",levelist="",
622                 repres="", date="",resol="",stream="",area="",time="",step="",expver="1",
623                 number="",accuracy="", grid="", gaussian="", target="",param="",dataset="",expect=""):
624        self.dataset=dataset
625        self.marsclass=marsclass
626        self.type=type
627        self.levtype=levtype
628        self.levelist=levelist
629        self.repres=repres
630        self.date=date
631        self.resol=resol
632        self.stream=stream
633        self.area=area
634        self.time=time
635        self.step=step
636        self.expver=expver
637        self.target=target
638        self.param=param
639        self.number=number
640        self.accuracy=accuracy
641        self.grid=grid
642        self.gaussian=gaussian
643#       self.expect=expect
644        self.server=server
645        self.public=public
646
647    def displayInfo(self):
648        attrs=vars(self)
649        for item in attrs.items():
650            if item[0] in ('server'):
651                pass
652            else:
653                print item[0]+': '+str(item[1])
654
655        return
656
657    def print_infodata_csv(self, inputdir, request_number):
658        '''
659        @Description:
660            Write all request parameter in alpabetical order into a "csv" file.
661
662        @Input:
663            self: instance of MarsRetrieval
664                For description see class documentation.
665
666            inputdir: string
667                The path where all data from the retrievals are stored.
668
669            request_number: integer
670                Number of mars requests for flux and non-flux data.
671
672        @Return:
673            <nothing>
674        '''
675
676        # Get all class attributes and their values as a dictionary
677        attrs = vars(self)
678        del attrs['server']
679
680        # open a file to store all requests to
681        with open(os.path.join(inputdir, 'mars_request.csv'), 'a') as f:
682            f.write(str(request_number) + ', ')
683            f.write(', '.join(str(attrs[key])
684                        for key in sorted(attrs.iterkeys())))
685            f.write('\n')
686
687        return
688
689    def dataRetrieve(self):
690        attrs=vars(self).copy()
691
692        del attrs['server']
693        del attrs['public']
694        mclass = attrs.get('marsclass')
695        del attrs['marsclass']
696        attrs['class'] = mclass
697        target = attrs.get('target')
698        if not int(self.public):
699            del attrs['target']
700        print 'target: ', target
701
702        delete_keys = []
703        for k, v in attrs.iteritems():
704            if v == '':
705                delete_keys.append(str(k))
706            else:
707                attrs[k] = str(v)
708
709        for k in delete_keys:
710            del attrs[k]
711
712        if self.server != False:
713            try:
714                if int(self.public):
715                    print 'RETRIEVE DATA!'
716                    self.server.retrieve(attrs)
717                else:
718                    print 'EXECUTE RETRIEVAL!'
719                    self.server.execute(attrs, target)
720            except:
721                e = sys.exc_info()[0]
722                print "Error: ", e
723                print 'MARS Request failed'#, have you already registered at apps.ecmwf.int?'
724                #raise IOError
725            if not int(self.public) and os.stat(target).st_size==0:
726                print 'MARS Request returned no data - please check request'
727                raise IOError
728            elif int(self.public) and os.stat(attrs.get('target')).st_size==0:
729                print 'MARS Request returned no data - please check request'
730                raise IOError
731        else:
732            s='ret'
733            for k,v in attrs.iteritems():
734                s=s+','+k+'='+str(v)
735
736            s+=',target="'+target+'"'
737            p=subprocess.Popen(['mars'], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,bufsize=1)
738            pout=p.communicate(input=s)[0]
739            print pout.decode()
740            if 'Some errors reported' in pout.decode():
741                print 'MARS Request failed - please check request'
742                raise IOError
743
744            if os.stat(target).st_size==0:
745                print 'MARS Request returned no data - please check request'
746                raise IOError
747
748        return
749
750##############################################################
751class EIFlexpart:
752    'class to retrieve Era Interim data for running FLEXPART'
753##############################################################
754
755    def __init__(self,c,fluxes=False):
756        # different mars types for retrieving reanalysis data for flexpart
757
758        # set a counter for the number of mars requests generated
759        self.mreq_count = 0
760        self.types=dict()
761        try:
762            if c.maxstep>len(c.type):    # Pure forecast mode
763                c.type=[c.type[0]] # AP changed this index from 1 to 0
764                c.step=['{:0>3}'.format(int(c.step[0]))]
765                c.time=[c.time[0]]
766                for i in range(1,c.maxstep+1):
767                    c.type.append(c.type[0])
768                    c.step.append('{:0>3}'.format(i))
769                    c.time.append(c.time[0])
770        except:
771            pass
772
773        self.inputdir=c.inputdir
774        self.basetime=c.basetime
775        self.dtime=c.dtime
776        self.mars={}
777        i=0
778        j=0
779        if fluxes==True and c.maxstep<=24:
780            self.types[str(c.acctype)]={'times':str(c.acctime),
781                                        'steps':'{}/to/{}/by/{}'.format(
782                                               c.dtime,c.accmaxstep,c.dtime)}
783            i=1
784            for k in [0,12]:
785                for j in range(int(c.dtime),13,int(c.dtime)):
786                    self.mars['{:0>3}'.format(i*int(c.dtime))]=[c.type[1],'{:0>2}'.format(k),'{:0>3}'.format(j)]
787                    i+=1
788        else:
789            for ty,st,ti in zip(c.type,c.step,c.time):
790                btlist=range(24)
791                if c.basetime=='12':
792                    btlist=[1,2,3,4,5,6,7,8,9,10,11,12]
793                if c.basetime=='00':
794                    btlist=[13,14,15,16,17,18,19,20,21,22,23,0]
795
796
797                if ((ty.upper() == 'AN' and mod(int(c.time[i]),int(c.dtime))==0) or \
798                   (ty.upper() != 'AN' and mod(int(c.step[i]),int(c.dtime))==0 and \
799                    mod(int(c.step[i]),int(c.dtime))==0) ) and \
800                   (int(c.time[i]) in btlist or c.maxstep>24):
801                    if ty not in self.types.keys():
802                        self.types[ty]={'times':'','steps':''}
803                    if ti not in self.types[ty]['times']:
804                        if len(self.types[ty]['times'])>0:
805                            self.types[ty]['times']+='/'
806                        self.types[ty]['times']+=ti
807                    if st not in self.types[ty]['steps']:
808                        if len(self.types[ty]['steps'])>0:
809                            self.types[ty]['steps']+='/'
810                        self.types[ty]['steps']+=st
811
812                    self.mars['{:0>3}'.format(j)]=[ty,'{:0>2}'.format(int(ti)),'{:0>3}'.format(int(st))]
813                    j+=int(c.dtime)
814
815                i+=1
816
817# Different grids need different retrievals
818# SH=Spherical Harmonics, GG=Gaussian Grid, OG=Output Grid, ML=MultiLevel, SL=SingleLevel
819        self.params={'SH__ML':'','SH__SL':'',
820                     'GG__ML':'','GG__SL':'',
821                     'OG__ML':'','OG__SL':'',
822                     'OG_OROLSM_SL':'','OG_acc_SL':''}
823
824        self.marsclass=c.marsclass
825        self.stream=c.stream
826        self.number=c.number
827        self.resol=c.resol
828        if 'N' in c.grid:  # Gaussian output grid
829            self.grid=c.grid
830            self.area='G'
831        else:
832            self.grid='{}/{}'.format(int(c.grid)/1000.,int(c.grid)/1000.)
833            self.area='{}/{}/{}/{}'.format(int(c.upper)/1000.,int(c.left)/1000.,int(c.lower)/1000.,int(c.right)/1000.)
834
835        self.accuracy=c.accuracy
836        self.level=c.level
837        try:
838            self.levelist=c.levelist
839        except:
840            self.levelist='1/to/'+c.level
841        self.glevelist='1/to/'+c.level
842        try:
843            self.gaussian=c.gaussian
844        except:
845            self.gaussian=''
846        try:
847            self.dataset=c.dataset
848        except:
849            self.dataset=''
850        try:
851            self.expver=c.expver
852        except:
853            self.expver='1'
854        try:
855            self.number=c.number
856        except:
857            self.number='0'
858
859        self.outputfilelist=[]
860
861
862# Now comes the nasty part that deals with the different scenarios we have:
863# 1) Calculation of etadot on
864#    a) Gaussian grid
865#    b) Output grid
866#    c) Output grid using parameter 77 retrieved from MARS
867# 3) Calculation/Retrieval of omega
868# 4) Download also data for WRF
869
870        if fluxes==False:
871            self.params['SH__SL']=['LNSP','ML','1','OFF']
872#           self.params['OG__SL']=["SD/MSL/TCC/10U/10V/2T/2D/129/172",'SFC','1',self.grid]
873            self.params['OG__SL']=["141/151/164/165/166/167/168/129/172",'SFC','1',self.grid]
874            if c.marsclass.upper() == 'EA' or c.marsclass.upper() == 'EP':
875                self.params['OG_OROLSM__SL']=["160/27/28/244",'SFC','1',self.grid]
876            else:
877                self.params['OG_OROLSM__SL']=["160/27/28/173",'SFC','1',self.grid]
878
879            if len(c.addpar)>0:
880                if c.addpar[0]=='/':
881                    c.addpar=c.addpar[1:]
882                self.params['OG__SL'][0]+='/'+'/'.join(c.addpar)
883            self.params['OG__ML']=['T/Q','ML',self.levelist,self.grid]
884
885            if c.gauss=='0' and c.eta=='1': # the simplest case
886                self.params['OG__ML'][0]+='/U/V/77'
887            elif c.gauss=='0' and c.eta=='0': # this is not recommended (inaccurate)
888                self.params['OG__ML'][0]+='/U/V'
889            elif c.gauss=='1' and c.eta=='0': # this is needed for data before 2008, or for reanalysis data
890                self.params['GG__SL']=['Q','ML','1','{}'.format((int(self.resol)+1)/2)]
891                self.params['SH__ML']=['U/V/D','ML',self.glevelist,'OFF']
892            else:
893                print 'Warning: This is a very costly parameter combination, use only for debugging!'
894                self.params['GG__SL']=['Q','ML','1','{}'.format((int(self.resol)+1)/2)]
895                self.params['GG__ML']=['U/V/D/77','ML',self.glevelist,'{}'.format((int(self.resol)+1)/2)]
896
897            if c.omega=='1':
898                self.params['OG__ML'][0]+='/W'
899
900            try:            # add cloud water content if necessary
901                if c.cwc=='1':
902                    self.params['OG__ML'][0]+='/CLWC/CIWC'
903            except:
904                pass
905
906            try:            # add vorticity and geopotential height for WRF if necessary
907                if c.wrf=='1':
908                    self.params['OG__ML'][0]+='/Z/VO'
909                    if '/D' not in self.params['OG__ML'][0]:
910                        self.params['OG__ML'][0]+='/D'
911#                   wrf_sfc='sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper()
912#                            134/151/235/167/165/166/168/129/172/034/031/141/139/170/183/236/039/040/041/042
913                    wrf_sfc='134/235/167/165/166/168/129/172/34/31/141/139/170/183/236/39/40/41/42'.upper()
914                    lwrt_sfc=wrf_sfc.split('/')
915                    for par in lwrt_sfc:
916                        if par not in self.params['OG__SL'][0]:
917                            self.params['OG__SL'][0]+='/'+par
918
919            except:
920                pass
921        else:
922            self.params['OG_acc_SL']=["LSP/CP/SSHF/EWSS/NSSS/SSR",'SFC','1',self.grid]
923
924
925        return
926        # add additional WRF specific parameters here
927
928    def create_namelist(self,c,filename):
929
930#       if inputdir=="":
931#           self.inputdir='.'
932#       else:
933#           self.inputdir=inputdir
934
935        self.inputdir=c.inputdir
936        area=asarray(self.area.split('/')).astype(float)
937        grid=asarray(self.grid.split('/')).astype(float)
938
939#       if area[3]<0:
940#           area[3]+=360
941        if area[1]>area[3]:
942            area[1]-=360
943        zyk=abs((area[3]-area[1]-360.)+grid[1])<1.e-6
944        maxl=int((area[3]-area[1])/grid[1])+1
945        maxb=int((area[0]-area[2])/grid[0])+1
946
947        f=open(self.inputdir+'/'+filename,'w')
948        f.write('&NAMGEN\n')
949        f.write(',\n  '.join(['maxl='+str(maxl),'maxb='+str(maxb),
950                              'mlevel='+self.level,'mlevelist='+'"'+self.levelist+'"',
951                    'mnauf='+self.resol,'metapar='+'77',
952                    'rlo0='+str(area[1]),'rlo1='+str(area[3]),'rla0='+str(area[2]),'rla1='+str(area[0]),
953                    'momega='+c.omega,'momegadiff='+c.omegadiff,'mgauss='+c.gauss,
954                    'msmooth='+c.smooth,'meta='+c.eta,'metadiff='+c.etadiff,'mdpdeta='+c.dpdeta]))
955
956        f.write('\n/\n')
957        f.close()
958        return
959
960
961    def retrieve(self, server, public, dates, request, times, inputdir=''):
962        self.dates=dates
963        self.server=server
964        self.public=public
965
966        if inputdir=="":
967            self.inputdir='.'
968        else:
969            self.inputdir=inputdir
970
971            # Retrieve Q not for using Q but as a template for a reduced gaussian grid one date and time is enough
972# Take analysis at 00
973#        qdate=self.dates
974#        idx=qdate.find("/")
975#        if (idx >0):
976#            qdate=self.dates[:idx]
977
978        #QG= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream, type="an", levtype="ML", levelist="1",
979                            #gaussian="reduced",grid='{}'.format((int(self.resol)+1)/2), resol=self.resol,accuracy=self.accuracy,target=self.inputdir+"/"+"QG.grb",
980                            #date=qdate, time="00",expver=self.expver, param="133.128")
981        #QG.displayInfo()
982        #QG.dataRetrieve()
983
984        oro=False
985        for ftype in self.types:
986            for pk,pv in self.params.iteritems():
987                if isinstance(pv,str):     # or pk!='GG__SL' :
988                    continue
989                mftype=''+ftype
990                mftime=self.types[ftype]['times']
991                mfstep=self.types[ftype]['steps']
992                mfdate=self.dates
993                mfstream=self.stream
994                mftarget=self.inputdir+"/"+ftype+pk+'.'+self.dates.split('/')[0]+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
995
996                if pk=='OG__SL':
997                    pass
998                if pk=='OG_OROLSM__SL':
999                    if oro==False:
1000                        # in CERA20C there is no stream "OPER"!
1001                        if self.marsclass.upper() != 'EP':
1002                            mfstream='OPER'
1003                        mftype='AN'
1004                        mftime='00'
1005                        mfstep='000'
1006                        mfdate=self.dates.split('/')[0]
1007                        mftarget=self.inputdir+"/"+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
1008                        oro=True
1009                    else:
1010                        continue
1011
1012                if pk=='GG__SL' and pv[0]=='Q':
1013                    area=""
1014                    gaussian='reduced'
1015                else:
1016                    area=self.area
1017                    gaussian=self.gaussian
1018
1019                if self.basetime==None:
1020                    # increase number of mars requests
1021                    self.mreq_count += 1
1022                    MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=mfstream,
1023                                      type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1024                                  accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1025                                  date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1026
1027                    if request == "0":
1028                        MR.display_info()
1029                        MR.data_retrieve()
1030                    elif request == "1":
1031                        MR.print_infodata_csv(self.inputdir, self.mreq_count)
1032                    elif request == "2":
1033                        MR.print_infodata_csv(self.inputdir, self.mreq_count)
1034                        MR.display_info()
1035                        MR.data_retrieve()
1036# The whole else section is only necessary for operational scripts. It could be removed
1037                else: # check if mars job requests fields beyond basetime. If yes eliminate those fields since they may not
1038                        # be accessible with user's credentials
1039                    sm1=-1
1040                    if 'by' in mfstep:
1041                        sm1=2
1042                    tm1=-1
1043                    if 'by' in mftime:
1044                        tm1=2
1045                    maxtime=datetime.datetime.strptime(mfdate.split('/')[-1]+mftime.split('/')[tm1],'%Y%m%d%H')+ \
1046                        datetime.timedelta(hours=int(mfstep.split('/')[sm1]))
1047
1048                    elimit=datetime.datetime.strptime(mfdate.split('/')[-1]+self.basetime,'%Y%m%d%H')
1049
1050                    if self.basetime=='12':
1051                        if 'acc' in pk:
1052
1053                    # Strategy: if maxtime-elimit>=24h reduce date by 1,
1054                    # if 12h<=maxtime-elimit<12h reduce time for last date
1055                    # if maxtime-elimit<12h reduce step for last time
1056                    # A split of the MARS job into 2 is likely necessary.
1057                            maxtime=elimit-datetime.timedelta(hours=24)
1058                            mfdate='/'.join(('/'.join(mfdate.split('/')[:-1]),datetime.datetime.strftime(maxtime,'%Y%m%d')))
1059                            # increase number of mars requests
1060                            self.mreq_count += 1
1061                            MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
1062                                              type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1063                                                accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1064                                                date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1065
1066                            if request == "0":
1067                                MR.display_info()
1068                                MR.data_retrieve()
1069                            elif request == "1":
1070                                MR.print_infodata_csv(self.inputdir, self.mreq_count)
1071                            elif request == "2":
1072                                MR.print_infodata_csv(self.inputdir, self.mreq_count)
1073                                MR.display_info()
1074                                MR.data_retrieve()
1075
1076                            maxtime=elimit-datetime.timedelta(hours=12)
1077                            mfdate=datetime.datetime.strftime(maxtime,'%Y%m%d')
1078                            mftime='00'
1079                            mftarget=self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
1080                            # increase number of mars requests
1081                            self.mreq_count += 1
1082                            MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
1083                                              type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1084                                                accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1085                                                date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1086                            if request == "0":
1087                                MR.display_info()
1088                                MR.data_retrieve()
1089                            elif request == "1":
1090                                MR.print_infodata_csv(self.inputdir, self.mreq_count)
1091                            elif request == "2":
1092                                MR.print_infodata_csv(self.inputdir, self.mreq_count)
1093                                MR.display_info()
1094                                MR.data_retrieve()
1095                        else:
1096                            # increase number of mars requests
1097                            self.mreq_count += 1
1098                            MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
1099                                              type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1100                                                accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1101                                                date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1102
1103                            if request == "0":
1104                                MR.display_info()
1105                                MR.data_retrieve()
1106                            elif request == "1":
1107                                MR.print_infodata_csv(self.inputdir, self.mreq_count)
1108                            elif request == "2":
1109                                MR.print_infodata_csv(self.inputdir, self.mreq_count)
1110                                MR.display_info()
1111                                MR.data_retrieve()
1112                    else:
1113                        maxtime=elimit-datetime.timedelta(hours=24)
1114                        mfdate=datetime.datetime.strftime(maxtime,'%Y%m%d')
1115
1116                        mftimesave=''.join(mftime)
1117
1118                        if '/' in mftime:
1119                            times=mftime.split('/')
1120                            while int(times[0])+int(mfstep.split('/')[0])<=12 and pk!='OG_OROLSM__SL' and 'acc' not in pk:
1121                                times=times[1:]
1122                                if len(times)>1:
1123                                    mftime='/'.join(times)
1124                                else:
1125                                    mftime=times[0]
1126                        # increase number of mars requests
1127                        self.mreq_count += 1
1128                        MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
1129                                          type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1130                                      accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1131                                      date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1132
1133                        if request == "0":
1134                            MR.display_info()
1135                            MR.data_retrieve()
1136                        elif request == "1":
1137                            MR.print_infodata_csv(self.inputdir, self.mreq_count)
1138                        elif request == "2":
1139                            MR.print_infodata_csv(self.inputdir, self.mreq_count)
1140                            MR.display_info()
1141                            MR.data_retrieve()
1142
1143                        if int(mftimesave.split('/')[0])==0 and int(mfstep.split('/')[0])==0 and pk!='OG_OROLSM__SL':
1144                            mfdate=datetime.datetime.strftime(elimit,'%Y%m%d')
1145                            mftime='00'
1146                            mfstep='000'
1147                            mftarget=self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
1148                            # increase number of mars requests
1149                            self.mreq_count += 1
1150                            MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
1151                                              type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1152                                          accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1153                                          date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1154
1155                            if request == "0":
1156                                MR.display_info()
1157                                MR.data_retrieve()
1158                            elif request == "1":
1159                                MR.print_infodata_csv(self.inputdir, self.mreq_count)
1160                            elif request == "2":
1161                                MR.print_infodata_csv(self.inputdir, self.mreq_count)
1162                                MR.display_info()
1163                                MR.data_retrieve()
1164
1165
1166        if request == "0" or request == "2":
1167            print('MARS retrieve done ... ')
1168        elif request == "1":
1169            print('MARS request printed ...')
1170
1171    def getFlexpartTime(self, type,step, time):
1172        cstep='{:0>3}'.format(step)
1173        ctime='{:0>2}'.format(int(time/100))
1174
1175        ctype = str(type).upper()
1176        myinfo = [ctype,ctime, cstep]
1177        cflextime = None
1178        for t, marsinfo in self.mars.items():
1179            if myinfo == marsinfo:
1180                cflextime=t
1181        return cflextime
1182
1183    def process_output(self, c):
1184
1185        print 'Postprocessing:\n Format: {}\n'.format(c.format)
1186        if c.ecapi==False:
1187            print 'ecstorage: {}\n ecfsdir: {}\n'.format(c.ecstorage,c.ecfsdir)
1188            if not hasattr(c,'gateway'):
1189                c.gateway=os.getenv('GATEWAY')
1190            if not hasattr(c,'destination'):
1191                c.destination=os.getenv('DESTINATION')
1192            print 'ectrans: {}\n gateway: {}\n destination: {}\n '.format(c.ectrans, c.gateway,c.destination)
1193        print 'Output filelist: \n',self.outputfilelist
1194
1195        if c.format.lower()=='grib2':
1196            for ofile in self.outputfilelist:
1197                p=subprocess.check_call(['grib_set','-s','edition=2,productDefinitionTemplateNumber=8',ofile,ofile+'_2'])
1198                p=subprocess.check_call(['mv',ofile+'_2',ofile])
1199
1200        if int(c.ectrans)==1 and c.ecapi==False:
1201            for ofile in self.outputfilelist:
1202                p=subprocess.check_call(['ectrans','-overwrite','-gateway',c.gateway,
1203                                         '-remote',c.destination,'-source',ofile])
1204                print 'ectrans:',p
1205        if int(c.ecstorage)==1 and c.ecapi==False:
1206            for ofile in self.outputfilelist:
1207                p=subprocess.check_call(['ecp','-o',ofile,os.path.expandvars(c.ecfsdir)])
1208
1209#20131107 000000      EN13110700              ON DISC
1210        if c.outputdir!=c.inputdir:
1211            for ofile in self.outputfilelist:
1212                p=subprocess.check_call(['mv',ofile,c.outputdir])
1213
1214        if c.grib2flexpart=='1':
1215            f=open(c.outputdir+'/'+'AVAILABLE','w')
1216            clist=[]
1217            for ofile in self.outputfilelist: # generate AVAILABLE file
1218                fname=ofile.split('/')
1219                if '.' in fname[-1]:
1220                    l=fname[-1].split('.')
1221                    timestamp=datetime.datetime.strptime(l[0][-6:]+l[1],'%y%m%d%H')
1222                    timestamp+=datetime.timedelta(hours=int(l[2]))
1223                    cdate=datetime.datetime.strftime(timestamp,'%Y%m%d')
1224                    chms=datetime.datetime.strftime(timestamp,'%H%M%S')
1225
1226                else:
1227                    cdate='20'+fname[-1][-8:-2]
1228                    chms=fname[-1][-2:]+'0000'
1229                clist.append(cdate+' '+chms+' '*6+fname[-1]+' '*14+'ON DISC')
1230            clist.sort()
1231            f.write('\n'.join(clist)+'\n')
1232            f.close()
1233
1234            pwd=os.path.abspath(c.outputdir)
1235            f=open(pwd+'/pathnames','w')
1236            f.write(pwd+'/Options/\n')
1237            f.write(pwd+'/\n')
1238            f.write(pwd+'/\n')
1239            f.write(pwd+'/AVAILABLE\n')
1240            f.write('==================\n')
1241            f.close()
1242
1243            if not os.path.exists(pwd+'/Options'):
1244                os.makedirs(pwd+'/Options')
1245            f=open(os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))+'/../Options/COMMAND','r')
1246            lflist=f.read().split('\n')
1247            i=0
1248            for l in lflist:
1249                if 'LDIRECT' in l.upper():
1250                    break
1251                i+=1
1252
1253            clist.sort()
1254            lflist=lflist[:i+1]+[clist[0][:16],clist[-1][:16]]+lflist[i+3:]
1255            g=open(pwd+'/Options/COMMAND','w')
1256            g.write('\n'.join(lflist)+'\n')
1257            g.close()
1258
1259            os.chdir(c.outputdir)
1260            p=subprocess.check_call([os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))+'/../FLEXPART_PROGRAM/grib2flexpart',
1261                                     'useAvailable','.'])
1262            os.chdir(pwd)
1263
1264    def create(self, inputfiles, c,):
1265
1266        table128=init128(c.flexextractdir+'/grib_templates/ecmwf_grib1_table_128')
1267        wrfpars=toparamId('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',table128)
1268#        index_keys=["date","time","stepRange"]
1269        index_keys=["date","time","step"]
1270        indexfile=c.inputdir+"/date_time_stepRange.idx"
1271        silentremove(indexfile)
1272        grib=GribTools(inputfiles.files)
1273        iid=grib.index(index_keys=index_keys, index_file = indexfile)
1274
1275        print 'index done...'
1276        fdict={'10':None,'11':None,'12':None,'13':None,'16':None,'17':None,'19':None,'21':None,'22':None,'20':None}
1277        for f in fdict.keys():
1278            silentremove(c.inputdir+"/fort."+f)
1279
1280        index_vals = []
1281        for key in index_keys:
1282            key_vals = grib_index_get(iid,key)
1283            print key_vals
1284
1285            index_vals.append(key_vals)
1286
1287        for prod in product(*index_vals):
1288            for i in range(len(index_keys)):
1289                grib_index_select(iid,index_keys[i],prod[i])
1290
1291
1292            gid = grib_new_from_index(iid)
1293            hid = gid
1294            cflextime = None
1295            for k,f in fdict.iteritems():
1296                fdict[k] = open(c.inputdir+'/fort.'+k,'w')
1297            if gid is not None:
1298                cdate = str(grib_get(gid, 'date'))
1299                time = grib_get(gid, 'time')
1300                type = grib_get(gid, 'type')
1301                step = grib_get(gid, 'step')
1302#                step = grib_get(gid, 'stepRange')
1303#               cflextime = self.getFlexpartTime(type,step, time)
1304                timestamp=datetime.datetime.strptime(cdate+'{:0>2}'.format(time/100),'%Y%m%d%H')
1305                timestamp+=datetime.timedelta(hours=int(step))
1306#               print gid,index_keys[i],prod[i],cdate,time,step,timestamp
1307
1308                cdateH=datetime.datetime.strftime(timestamp,'%Y%m%d%H')
1309                chms=datetime.datetime.strftime(timestamp,'%H%M%S')
1310
1311                if c.basetime !=None:
1312                    slimit=datetime.datetime.strptime(c.start_date+'00','%Y%m%d%H')
1313                    bt='23'
1314                    if c.basetime=='00':
1315                        bt='00'
1316                        slimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')-datetime.timedelta(hours=12-int(c.dtime))
1317
1318                    if c.basetime=='12':
1319                        bt='12'
1320                        slimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')-datetime.timedelta(hours=12-int(c.dtime))
1321
1322                    elimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')
1323
1324                    if timestamp<slimit or timestamp>elimit:
1325                        continue
1326
1327
1328
1329            try:
1330                if c.wrf=='1':
1331                    if 'olddate' not in locals():
1332                        fwrf = open(c.outputdir+'/WRF'+cdate+'.{:0>2}'.format(time)+'.000.grb2','w')
1333                        olddate=cdate[:]
1334                    else:
1335                        if cdate!=olddate:
1336                            fwrf = open(c.outputdir+'/WRF'+cdate+'.{:0>2}'.format(time)+'.000.grb2','w')
1337                            olddate=cdate[:]
1338            except AttributeError:
1339                pass
1340
1341#                print 'cyear '+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime
1342
1343            savedfields=[]
1344            while 1:
1345                if gid is None: break
1346                paramId = grib_get(gid, 'paramId')
1347                gridtype = grib_get(gid, 'gridType')
1348                datatype = grib_get(gid, 'dataType')
1349                levtype = grib_get(gid, 'typeOfLevel')
1350                if paramId == 133 and gridtype=='reduced_gg':
1351# Relative humidity (Q.grb) is used as a template only so we need the first we "meet"
1352                    fout=open(c.inputdir+'/fort.18','w')
1353                    grib_write(gid,fout)
1354                    fout.close()
1355                elif paramId == 131 or paramId == 132:
1356                    grib_write(gid, fdict['10'])
1357                elif paramId == 130:
1358                    grib_write(gid, fdict['11'])
1359                elif paramId == 133 and gridtype!='reduced_gg':
1360                    grib_write(gid, fdict['17'])
1361                elif paramId == 152:
1362                    grib_write(gid, fdict['12'])
1363                elif paramId == 155 and gridtype=='sh':
1364                    grib_write(gid, fdict['13'])
1365                elif paramId in [129,138,155] and levtype=='hybrid' and c.wrf=='1':
1366#                   print paramId,'not written'
1367                    pass
1368                elif paramId == 246 or paramId == 247:  # cloud liquid water and ice
1369                    if paramId==246:
1370                        clwc=grib_get_values(gid)
1371                    else:
1372                        clwc+=grib_get_values(gid)
1373                        grib_set_values(gid,clwc)
1374#                       grib_set(gid,'shortName','qc')
1375                        grib_set(gid,'paramId',201031)
1376                        grib_write(gid, fdict['22'])
1377
1378                elif paramId == 135:
1379                    grib_write(gid, fdict['19'])
1380                elif paramId == 77:
1381                    grib_write(gid, fdict['21'])
1382                else:
1383                    if paramId not in savedfields:
1384                        grib_write(gid, fdict['16'])
1385                        savedfields.append(paramId)
1386                    else:
1387                        print 'duplicate '+str(paramId)+' not written'
1388
1389                try:
1390                    if c.wrf=='1':
1391                        if levtype=='hybrid':
1392                            if paramId in [129,130,131,132,133,138,155]:
1393                                grib_write(gid,fwrf)
1394                        else:
1395                            if paramId in wrfpars:
1396                                grib_write(gid,fwrf)
1397                except AttributeError:
1398                    pass
1399
1400
1401
1402                grib_release(gid)
1403                gid = grib_new_from_index(iid)
1404# call for CONVERT2
1405
1406            for f in fdict.values():
1407                f.close()
1408
1409            if hid is not None:
1410                pwd=os.getcwd()
1411                os.chdir(c.inputdir)
1412                if os.stat('fort.21').st_size==0 and int(c.eta)==1:
1413                    print 'Parameter 77 (etadot) is missing, most likely it is not available for this type or date/time\n'
1414                    print 'Check parameters CLASS, TYPE, STREAM, START_DATE\n'
1415                    myerror(c,'fort.21 is empty while parameter eta is set to 1 in CONTROL file')
1416
1417                p=subprocess.check_call([os.path.expandvars(os.path.expanduser(c.exedir))+'/CONVERT2'],shell=True)
1418                os.chdir(pwd)
1419# create the corresponding output file  fort.15 (generated by CONVERT2) + fort.16 (paramId 167 and paramId 168)
1420                fnout=c.inputdir+'/'+c.prefix
1421                if c.maxstep>12:
1422                    suffix=cdate[2:8]+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step)
1423                else:
1424                    suffix=cdateH[2:10]
1425
1426                fnout+=suffix
1427                print "outputfile = " + fnout
1428                self.outputfilelist.append(fnout) # needed for final processing
1429                fout = open(fnout,'wb')
1430                shutil.copyfileobj(open(c.inputdir+'/fort.15','rb'), fout)
1431                if c.cwc=='1':
1432                    shutil.copyfileobj(open(c.inputdir+'/fort.22','rb'), fout)
1433                shutil.copyfileobj(open(c.inputdir+'/flux'+cdate[0:2]+suffix,'rb'), fout)
1434                shutil.copyfileobj(open(c.inputdir+'/fort.16','rb'), fout)
1435                orolsm=glob.glob(c.inputdir+'/OG_OROLSM__SL.*.'+c.ppid+'*')[0]
1436                shutil.copyfileobj(open(orolsm,'rb'), fout)
1437                fout.close()
1438                if c.omega=='1':
1439                    fnout=c.outputdir+'/OMEGA'+suffix
1440                    fout = open(fnout,'wb')
1441                    shutil.copyfileobj(open(c.inputdir+'/fort.25','rb'), fout)
1442
1443
1444
1445        try:
1446            if c.wrf=='1':
1447                fwrf.close()
1448        except:
1449            pass
1450
1451        grib_index_release(iid)
1452
1453    def deacc_fluxes(self, inputfiles, c):
1454
1455        table128=init128(c.flexextractdir+'/grib_templates/ecmwf_grib1_table_128')
1456        pars=toparamId(self.params['OG_acc_SL'][0],table128)
1457        index_keys=["date","time","step"]
1458        indexfile=c.inputdir+"/date_time_stepRange.idx"
1459        silentremove(indexfile)
1460        grib=GribTools(inputfiles.files)
1461        iid=grib.index(index_keys=index_keys, index_file = indexfile)
1462
1463        print 'index done...'
1464
1465        index_vals = []
1466        for key in index_keys:
1467            key_vals = grib_index_get(iid,key)
1468            print key_vals
1469            if key=='step':
1470                l=[]
1471                for k in key_vals:
1472                    l.append(int(k))
1473                l.sort()
1474                key_vals=[]
1475                for k in l:
1476                    key_vals.append(str(k))
1477
1478            index_vals.append(key_vals)
1479
1480
1481        valsdict={}
1482        svalsdict={}
1483        stepsdict={}
1484        for p in pars:
1485            valsdict[str(p)]=[]
1486            svalsdict[str(p)]=[]
1487            stepsdict[str(p)]=[]
1488
1489        for prod in product(*index_vals):
1490            for i in range(len(index_keys)):
1491                grib_index_select(iid,index_keys[i],prod[i])
1492
1493            #for k,f in fdict.iteritems():
1494                #fdict[k] = open('fort.'+k,'w')
1495            gid = grib_new_from_index(iid)
1496            hid = gid
1497            cflextime = None
1498            if gid is not None:
1499                cdate = grib_get(gid, 'date')
1500                #cyear = cdate[:4]
1501                #cmonth = cdate[4:6]
1502                #cday = cdate[6:8]
1503                time = grib_get(gid, 'time')
1504                type = grib_get(gid, 'type')
1505                step = grib_get(gid, 'step')
1506                # date+time+step-2*dtime (since interpolated value valid for step-2*dtime)
1507                sdate=datetime.datetime(year=cdate/10000,month=mod(cdate,10000)/100,
1508                                        day=mod(cdate,100),hour=time/100)
1509                fdate=sdate+datetime.timedelta(hours=step-2*int(c.dtime))
1510                sdates=sdate+datetime.timedelta(hours=step)
1511            else:
1512                break
1513
1514            fnout=c.inputdir+'/'
1515            if c.maxstep>12:
1516                fnout+='flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step-2*int(c.dtime))
1517                gnout=c.inputdir+'/flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step-int(c.dtime))
1518                hnout=c.inputdir+'/flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step)
1519                g=open(gnout,'w')
1520                h=open(hnout,'w')
1521            else:
1522                fnout+='flux'+fdate.strftime('%Y%m%d%H')
1523                gnout=c.inputdir+'/flux'+(fdate+datetime.timedelta(hours=int(c.dtime))).strftime('%Y%m%d%H')
1524                hnout=c.inputdir+'/flux'+sdates.strftime('%Y%m%d%H')
1525                g=open(gnout,'w')
1526                h=open(hnout,'w')
1527            print "outputfile = " + fnout
1528            f=open(fnout,'w')
1529
1530            while 1:
1531                if gid is None: break
1532                cparamId = str(grib_get(gid, 'paramId'))
1533                step = grib_get(gid, 'step')
1534                atime = grib_get(gid, 'time')
1535                ni=grib_get(gid, 'Ni')
1536                nj=grib_get(gid, 'Nj')
1537                if cparamId in valsdict.keys():
1538                    values = grib_get_values(gid)
1539                    vdp=valsdict[cparamId]
1540                    svdp=svalsdict[cparamId]
1541                    sd=stepsdict[cparamId]
1542
1543                    if cparamId=='142' or cparamId=='143':
1544                        fak=1./1000.
1545                    else:
1546                        fak=3600.
1547
1548                    values=(reshape(values,(nj,ni))).flatten()/fak
1549                    vdp.append(values[:]) # save the accumulated values
1550                    if c.marsclass.upper() == 'EA' or \
1551                       step<=int(c.dtime):
1552                        svdp.append(values[:]/int(c.dtime))
1553                    else:
1554                        svdp.append((vdp[-1]-vdp[-2])/int(c.dtime))
1555
1556                    print cparamId,atime,step,len(values),values[0],std(values)
1557                    #svdp.append(values[:]) # save the 1-hourly or 3-hourly specific values
1558                    sd.append(step)
1559                    if len(svdp)>=3:
1560                        if len(svdp)>3:
1561                            if cparamId=='142' or cparamId=='143':
1562                                values=darain(svdp)
1563                            else:
1564                                values=dapoly(svdp)
1565
1566                            if not (step==c.maxstep and c.maxstep>12 or sdates==elimit):
1567                                vdp.pop(0)
1568                                svdp.pop(0)
1569                        else:
1570                            if c.maxstep>12:
1571                                values=svdp[1]
1572                            else:
1573                                values=svdp[0]
1574
1575                        grib_set_values(gid, values)
1576                        if c.maxstep>12:
1577                            grib_set(gid,'step',max(0,step-2*int(c.dtime)))
1578                        else:
1579                            grib_set(gid,'step',0)
1580                            grib_set(gid,'time',fdate.hour*100)
1581                            grib_set(gid,'date',fdate.year*10000+fdate.month*100+fdate.day)
1582                        grib_write(gid, f)
1583
1584                        if c.basetime is not None:
1585                            elimit=datetime.datetime.strptime(c.end_date+c.basetime,'%Y%m%d%H')
1586                        else:
1587                            elimit=sdate+datetime.timedelta(2*int(c.dtime))
1588
1589                        # squeeze out information of last two steps contained in svdp
1590    #                   if step+int(c.dtime)==c.maxstep and c.maxstep>12 or sdates+datetime.timedelta(hours=int(c.dtime))>=elimit:
1591    #                       Note that svdp[0] has not been popped in this case
1592
1593
1594                        if step==c.maxstep and c.maxstep>12 or sdates==elimit:
1595                            values=svdp[3]
1596                            grib_set_values(gid, values)
1597                            grib_set(gid,'step',0)
1598                            truedatetime=fdate+datetime.timedelta(hours=2*int(c.dtime))
1599                            grib_set(gid,'time',truedatetime.hour*100)
1600                            grib_set(gid,'date',truedatetime.year*10000+truedatetime.month*100+truedatetime.day)
1601                            grib_write(gid, h)
1602
1603                            #values=(svdp[1]+svdp[2])/2.
1604                            if cparamId=='142' or cparamId=='143':
1605                                values=darain(list(reversed(svdp)))
1606                            else:
1607                                values=dapoly(list(reversed(svdp)))
1608
1609                            grib_set(gid,'step',0)
1610                            truedatetime=fdate+datetime.timedelta(hours=int(c.dtime))
1611                            grib_set(gid,'time',truedatetime.hour*100)
1612                            grib_set(gid,'date',truedatetime.year*10000+truedatetime.month*100+truedatetime.day)
1613                            grib_set_values(gid, values)
1614                            grib_write(gid, g)
1615
1616                    grib_release(gid)
1617
1618                gid = grib_new_from_index(iid)
1619
1620            f.close()
1621            g.close()
1622            h.close()
1623
1624
1625        grib_index_release(iid)
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG