source: flex_extract.git/python/FlexpartTools.py @ aa1af8c

ctbtodev
Last change on this file since aa1af8c was 68f5e1c, checked in by anphi <anne.philipp@…>, 6 years ago

BugFix? of Reported bug in Ticket #202

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