Context Navigation

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

ctbtodev

Last change on this file since 51f9853 was 51f9853, checked in by anphi <anne.philipp@…>, 6 years ago
added request parameter for writing mars requests into csv file
Property mode set to `100644`
File size: 66.6 KB

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

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

Download in other formats: