Context Navigation

source: flex_extract.git/python/FlexpartTools.py @ 9e660c4

ctbtodev

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

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

Download in other formats: