Context Navigation

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

ctbtodev

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: