Context Navigation

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

ctbtodev

Last change on this file since 2193808 was 2193808, checked in by anphi <anne.philipp@…>, 5 years ago
added ensemble analysis/forecast retrievement option
Property mode set to `100644`
File size: 69.4 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 isinstance(self.type, list):
512	self.type = [self.type]
513	if not isinstance(self.time, list):
514	self.time = [self.time]
515	if not isinstance(self.step, list):
516	self.step = [self.step]
517
518	max_level_list = [16, 19, 31, 40, 50, 60, 62, 91, 137]
519	if not hasattr(self,'levelist'):
520	if not hasattr(self,'level'):
521	print 'WARNING: neither levelist nor level specified in CONTROL file'
522	else:
523	if self.level in max_level_list:
524	self.levelist='1/to/'+self.level
525	else:
526	if not hasattr(self, 'level'):
527	if 'to' in self.levelist.lower():
528	max_level = self.levelist.split('/')[2]
529	if int(max_level) in max_level_list:
530	self.level = max_level
531	else:
532	print 'ERROR: LEVEL-parameter could not be set'
533	print 'LEVEL must be the maximum level of a specified level list from ECMWF, e.g.'
534	print '16, 19, 31, 40, 50, 60, 62, 91 or 137'
535	sys.exit(1)
536	else:
537	max_level = self.levelist.split('/')[-1]
538	if int(max_level) in max_level_list:
539	self.level = max_level
540	else:
541	print 'ERROR: LEVEL-parameter could not be set'
542	print 'LEVEL must be the maximum level of a specified level list from ECMWF, e.g.'
543	print '16, 19, 31, 40, 50, 60, 62, 91 or 137'
544	sys.exit(1)
545	else:
546	# check if consistent
547	if int(self.level) in max_level_list:
548	pass
549	else:
550	print 'ERROR: LEVEL-parameter is not properly selected'
551	print 'LEVEL must be the maximum level of a specified level list from ECMWF, e.g.'
552	print '16, 19, 31, 40, 50, 60, 62, 91 or 137'
553	sys.exit(1)
554
555	if not hasattr(self,'maxstep'):
556	self.maxstep=0
557	for s in self.step:
558	if int(s)>self.maxstep:
559	self.maxstep=int(s)
560	else:
561	self.maxstep=int(self.maxstep)
562	if not hasattr(self,'prefix'):
563	self.prefix='EN'
564	if not hasattr(self,'makefile'):
565	self.makefile=None
566	if not hasattr(self,'basetime'):
567	self.basetime=None
568	if not hasattr(self,'date_chunk'):
569	self.date_chunk='3'
570
571	self.flexextractdir=os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))+'/../' # script directory
572	# if not hasattr(self,'exedir'):
573	self.exedir=self.flexextractdir+'src/'
574	if not hasattr(self,'grib2flexpart'):
575	self.grib2flexpart='0'
576	if not hasattr(self,'flexpart_root_scripts'):
577	self.flexpart_root_scripts=self.flexextractdir
578
579	if not hasattr(self,'acctype'):
580	print('... Control paramter ACCTYPE was not defined in CONTROL file.')
581	print('Important for (forecast) flux data extraction.')
582	try:
583	if len(self.type) > 1 and self.type[1] != 'AN':
584	print('Use old setting by using TYPE[1].')
585	self.acctype = self.type[1]
586	except:
587	print('Use default value "FC"!')
588	self.acctype='FC'
589
590	if not hasattr(self,'acctime'):
591	print('... Control paramter ACCTIME was not defined in CONTROL file.')
592	print('Use default value "00/12" for flux forecast!')
593	print('This may be different for single data sets, see documentation if retrieval fails!')
594	self.acctime='00/12'
595
596	if not hasattr(self, 'accmaxstep'):
597	print('... Control paramter ACCMAXSTEP was not defined in CONTROL file.')
598	print('Use default value "12" for flux forecast!')
599	print('This may be different for single data sets, see documentation if retrieval fails!')
600	self.accmaxstep='12'
601
602	for i in range(len(self.type)):
603	if self.type[i] == 'AN' and int(self.step[i]) != 0:
604	print('Analysis retrievals must have STEP = 0 (setting to 0)')
605	self.type[i] = 0
606
607	if not hasattr(self,'request'):
608	self.request=0
609	elif int(self.request) != 0:
610	self.request=int(self.request)
611	marsfile = os.path.join(self.inputdir,
612	'mars_requests.csv')
613	if os.path.isfile(marsfile):
614	silentremove(marsfile)
615
616	return
617	def __str__(self):
618
619	attrs = vars(self)
620	# {'kids': 0, 'name': 'Dog', 'color': 'Spotted', 'age': 10, 'legs': 2, 'smell': 'Alot'}
621	# now dump this in some way or another
622	return ', '.join("%s: %s" % item for item in attrs.items())
623
624	def tolist(self):
625	attrs=vars(self)
626	l=list()
627	for item in attrs.items():
628	if '_expanded' in item[0]:
629	pass
630	elif 'exedir' in item[0]:
631	pass
632	elif 'flexpart_root_scripts' in item[0]:
633	pass
634	elif 'flexextractdir' in item[0]:
635	pass
636
637	else:
638	if type(item[1]) is list:
639	stot=''
640	for s in item[1]:
641	stot+=s+' '
642
643	l.append("%s %s" % (item[0],stot))
644	else:
645	l.append("%s %s" % item )
646	return sorted(l)
647
648
649
650	##############################################################
651	# MARSretrieval class
652	##############################################################
653	class MARSretrieval:
654	'class for MARS retrievals'
655
656	def __init__(self,server, public, marsclass="ei",type="",levtype="",levelist="",
657	repres="", date="",resol="",stream="",area="",time="",step="",expver="1",
658	number="",accuracy="", grid="", gaussian="", target="",param="",dataset="",expect=""):
659	self.dataset=dataset
660	self.marsclass=marsclass
661	self.type=type
662	self.levtype=levtype
663	self.levelist=levelist
664	self.repres=repres
665	self.date=date
666	self.resol=resol
667	self.stream=stream
668	self.area=area
669	self.time=time
670	self.step=step
671	self.expver=expver
672	self.target=target
673	self.param=param
674	self.number=number
675	self.accuracy=accuracy
676	self.grid=grid
677	self.gaussian=gaussian
678	# self.expect=expect
679	self.server=server
680	self.public=public
681
682	def displayInfo(self):
683	attrs=vars(self)
684	for item in attrs.items():
685	if item[0] in ('server'):
686	pass
687	else:
688	print item[0]+': '+str(item[1])
689
690	return
691
692	def print_infodata_csv(self, inputdir, request_number):
693	'''
694	@Description:
695	Write all request parameter in alpabetical order into a "csv" file.
696
697	@Input:
698	self: instance of MarsRetrieval
699	For description see class documentation.
700
701	inputdir: string
702	The path where all data from the retrievals are stored.
703
704	request_number: integer
705	Number of mars requests for flux and non-flux data.
706
707	@Return:
708	<nothing>
709	'''
710
711	# Get all class attributes and their values as a dictionary
712	attrs = vars(self).copy()
713	del attrs['server']
714	del attrs['public']
715
716	# open a file to store all requests to
717	with open(os.path.join(inputdir, 'mars_requests.csv'), 'a') as f:
718	f.write(str(request_number) + ', ')
719	f.write(', '.join(str(attrs[key])
720	for key in sorted(attrs.iterkeys())))
721	f.write('\n')
722
723	return
724
725	def dataRetrieve(self):
726	attrs=vars(self).copy()
727
728	del attrs['server']
729	del attrs['public']
730	mclass = attrs.get('marsclass')
731	del attrs['marsclass']
732	attrs['class'] = mclass
733
734	add_previousday_ifneeded(attrs)
735
736	target = attrs.get('target')
737	if not int(self.public):
738	del attrs['target']
739	print 'target: ', target
740
741	delete_keys = []
742	for k, v in attrs.iteritems():
743	if v == '':
744	delete_keys.append(str(k))
745	else:
746	attrs[k] = str(v)
747
748	for k in delete_keys:
749	del attrs[k]
750
751	if self.server != False:
752	try:
753	if int(self.public):
754	print 'RETRIEVE DATA!'
755	self.server.retrieve(attrs)
756	else:
757	print 'EXECUTE RETRIEVAL!'
758	self.server.execute(attrs, target)
759	except:
760	e = sys.exc_info()[0]
761	print "Error: ", e
762	print 'MARS Request failed'#, have you already registered at apps.ecmwf.int?'
763	#raise IOError
764	if not int(self.public) and os.stat(target).st_size==0:
765	print 'MARS Request returned no data - please check request'
766	raise IOError
767	elif int(self.public) and os.stat(attrs.get('target')).st_size==0:
768	print 'MARS Request returned no data - please check request'
769	raise IOError
770	else:
771	s='ret'
772	for k,v in attrs.iteritems():
773	s=s+','+k+'='+str(v)
774
775	s+=',target="'+target+'"'
776	p=subprocess.Popen(['mars'], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,bufsize=1)
777	pout=p.communicate(input=s)[0]
778	print pout.decode()
779	if 'Some errors reported' in pout.decode():
780	print 'MARS Request failed - please check request'
781	raise IOError
782
783	if os.stat(target).st_size==0:
784	print 'MARS Request returned no data - please check request'
785	raise IOError
786
787	return
788
789	##############################################################
790	class EIFlexpart:
791	'class to retrieve Era Interim data for running FLEXPART'
792	##############################################################
793
794	def __init__(self,c,fluxes=False):
795	# different mars types for retrieving reanalysis data for flexpart
796
797	# set a counter for the number of mars requests generated
798	self.mreq_count = 0
799	self.types=dict()
800	try:
801	if c.maxstep>24: #len(c.type): # Pure forecast mode
802	c.type=[c.type[0]] # AP changed this index from 1 to 0
803	c.step=['{:0>3}'.format(int(c.step[0]))]
804	c.time=[c.time[0]]
805	for i in range(1,c.maxstep+1):
806	c.type.append(c.type[0])
807	c.step.append('{:0>3}'.format(i))
808	c.time.append(c.time[0])
809	except:
810	pass
811
812	self.inputdir=c.inputdir
813	self.basetime=c.basetime
814	self.dtime=c.dtime
815	self.mars={}
816	i=0
817	j=0
818	if fluxes==True and c.maxstep<=24:
819	self.types[str(c.acctype)]={'times':str(c.acctime),
820	'steps':'{}/to/{}/by/{}'.format(
821	c.dtime,c.accmaxstep,c.dtime)}
822	i=1
823	for k in [0,12]:
824	for j in range(int(c.dtime),13,int(c.dtime)):
825	self.mars['{:0>3}'.format(i*int(c.dtime))]=[c.type[1],'{:0>2}'.format(k),'{:0>3}'.format(j)]
826	i+=1
827	else:
828	for ty,st,ti in zip(c.type,c.step,c.time):
829	btlist=range(24)
830	if c.basetime=='12':
831	btlist=[1,2,3,4,5,6,7,8,9,10,11,12]
832	if c.basetime=='00':
833	btlist=[13,14,15,16,17,18,19,20,21,22,23,0]
834
835
836	if ((ty.upper() == 'AN' and mod(int(c.time[i]),int(c.dtime))==0 and int(c.step[i])==0) or \
837	(ty.upper() != 'AN' and mod(int(c.step[i]),int(c.dtime))==0 and \
838	mod(int(c.step[i]),int(c.dtime))==0) ) and \
839	(int(c.time[i]) in btlist or c.maxstep>24):
840	if ty not in self.types.keys():
841	self.types[ty]={'times':'','steps':''}
842	if ti not in self.types[ty]['times']:
843	if len(self.types[ty]['times'])>0:
844	self.types[ty]['times']+='/'
845	self.types[ty]['times']+=ti
846	if st not in self.types[ty]['steps']:
847	if len(self.types[ty]['steps'])>0:
848	self.types[ty]['steps']+='/'
849	self.types[ty]['steps']+=st
850
851	self.mars['{:0>3}'.format(j)]=[ty,'{:0>2}'.format(int(ti)),'{:0>3}'.format(int(st))]
852	j+=int(c.dtime)
853
854	i+=1
855
856	# Different grids need different retrievals
857	# SH=Spherical Harmonics, GG=Gaussian Grid, OG=Output Grid, ML=MultiLevel, SL=SingleLevel
858	self.params={'SH__ML':'','SH__SL':'',
859	'GG__ML':'','GG__SL':'',
860	'OG__ML':'','OG__SL':'',
861	'OG_OROLSM_SL':'','OG_acc_SL':''}
862
863	self.marsclass=c.marsclass
864	self.stream=c.stream
865	self.number=c.number
866	self.resol=c.resol
867	if 'N' in c.grid: # Gaussian output grid
868	self.grid=c.grid
869	self.area='G'
870	else:
871	self.grid='{}/{}'.format(int(c.grid)/1000.,int(c.grid)/1000.)
872	self.area='{}/{}/{}/{}'.format(int(c.upper)/1000.,int(c.left)/1000.,int(c.lower)/1000.,int(c.right)/1000.)
873
874	self.accuracy=c.accuracy
875	self.level=c.level
876	try:
877	self.levelist=c.levelist
878	except:
879	self.levelist='1/to/'+c.level
880	self.glevelist='1/to/'+c.level
881	try:
882	self.gaussian=c.gaussian
883	except:
884	self.gaussian=''
885	try:
886	self.dataset=c.dataset
887	except:
888	self.dataset=''
889	try:
890	self.expver=c.expver
891	except:
892	self.expver='1'
893	try:
894	self.number=c.number
895	except:
896	self.number='0'
897
898	self.outputfilelist=[]
899
900
901	# Now comes the nasty part that deals with the different scenarios we have:
902	# 1) Calculation of etadot on
903	# a) Gaussian grid
904	# b) Output grid
905	# c) Output grid using parameter 77 retrieved from MARS
906	# 3) Calculation/Retrieval of omega
907	# 4) Download also data for WRF
908
909	if fluxes==False:
910	self.params['SH__SL']=['LNSP','ML','1','OFF']
911	# self.params['OG__SL']=["SD/MSL/TCC/10U/10V/2T/2D/129/172",'SFC','1',self.grid]
912	self.params['OG__SL']=["141/151/164/165/166/167/168/129/172",'SFC','1',self.grid]
913	if c.marsclass.upper() == 'EA' or c.marsclass.upper() == 'EP':
914	self.params['OG_OROLSM__SL']=["160/27/28/244",'SFC','1',self.grid]
915	else:
916	self.params['OG_OROLSM__SL']=["160/27/28/173",'SFC','1',self.grid]
917
918	if len(c.addpar)>0:
919	if c.addpar[0]=='/':
920	c.addpar=c.addpar[1:]
921	self.params['OG__SL'][0]+='/'+'/'.join(c.addpar)
922	self.params['OG__ML']=['T/Q','ML',self.levelist,self.grid]
923
924	if c.gauss=='0' and c.eta=='1': # the simplest case
925	self.params['OG__ML'][0]+='/U/V/77'
926	elif c.gauss=='0' and c.eta=='0': # this is not recommended (inaccurate)
927	self.params['OG__ML'][0]+='/U/V'
928	elif c.gauss=='1' and c.eta=='0': # this is needed for data before 2008, or for reanalysis data
929	self.params['GG__SL']=['Q','ML','1','{}'.format((int(self.resol)+1)/2)]
930	self.params['SH__ML']=['U/V/D','ML',self.glevelist,'OFF']
931	else:
932	print 'Warning: This is a very costly parameter combination, use only for debugging!'
933	self.params['GG__SL']=['Q','ML','1','{}'.format((int(self.resol)+1)/2)]
934	self.params['GG__ML']=['U/V/D/77','ML',self.glevelist,'{}'.format((int(self.resol)+1)/2)]
935
936	if c.omega=='1':
937	self.params['OG__ML'][0]+='/W'
938
939	try: # add cloud water content if necessary
940	if c.cwc=='1':
941	self.params['OG__ML'][0]+='/CLWC/CIWC'
942	except:
943	pass
944
945	try: # add vorticity and geopotential height for WRF if necessary
946	if c.wrf=='1':
947	self.params['OG__ML'][0]+='/Z/VO'
948	if '/D' not in self.params['OG__ML'][0]:
949	self.params['OG__ML'][0]+='/D'
950	# wrf_sfc='sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper()
951	# 134/151/235/167/165/166/168/129/172/034/031/141/139/170/183/236/039/040/041/042
952	wrf_sfc='134/235/167/165/166/168/129/172/34/31/141/139/170/183/236/39/40/41/42'.upper()
953	lwrt_sfc=wrf_sfc.split('/')
954	for par in lwrt_sfc:
955	if par not in self.params['OG__SL'][0]:
956	self.params['OG__SL'][0]+='/'+par
957
958	except:
959	pass
960	else:
961	self.params['OG_acc_SL']=["LSP/CP/SSHF/EWSS/NSSS/SSR",'SFC','1',self.grid]
962
963
964	return
965	# add additional WRF specific parameters here
966
967	def create_namelist(self,c,filename):
968
969	# if inputdir=="":
970	# self.inputdir='.'
971	# else:
972	# self.inputdir=inputdir
973
974	self.inputdir=c.inputdir
975	area=asarray(self.area.split('/')).astype(float)
976	grid=asarray(self.grid.split('/')).astype(float)
977
978	# if area[3]<0:
979	# area[3]+=360
980	if area[1]>area[3]:
981	area[1]-=360
982	zyk=abs((area[3]-area[1]-360.)+grid[1])<1.e-6
983	maxl=int((area[3]-area[1])/grid[1])+1
984	maxb=int((area[0]-area[2])/grid[0])+1
985
986	f=open(self.inputdir+'/'+filename,'w')
987	f.write('&NAMGEN\n')
988	f.write(',\n '.join(['maxl='+str(maxl),'maxb='+str(maxb),
989	'mlevel='+self.level,'mlevelist='+'"'+self.levelist+'"',
990	'mnauf='+self.resol,'metapar='+'77',
991	'rlo0='+str(area[1]),'rlo1='+str(area[3]),'rla0='+str(area[2]),'rla1='+str(area[0]),
992	'momega='+c.omega,'momegadiff='+c.omegadiff,'mgauss='+c.gauss,
993	'msmooth='+c.smooth,'meta='+c.eta,'metadiff='+c.etadiff,'mdpdeta='+c.dpdeta]))
994
995	f.write('\n/\n')
996	f.close()
997	return
998
999
1000	def retrieve(self, server, public, dates, request, times, inputdir=''):
1001	self.dates=dates
1002	self.server=server
1003	self.public=public
1004
1005	if inputdir=="":
1006	self.inputdir='.'
1007	else:
1008	self.inputdir=inputdir
1009
1010	# Retrieve Q not for using Q but as a template for a reduced gaussian grid one date and time is enough
1011	# Take analysis at 00
1012	# qdate=self.dates
1013	# idx=qdate.find("/")
1014	# if (idx >0):
1015	# qdate=self.dates[:idx]
1016
1017	#QG= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream, type="an", levtype="ML", levelist="1",
1018	#gaussian="reduced",grid='{}'.format((int(self.resol)+1)/2), resol=self.resol,accuracy=self.accuracy,target=self.inputdir+"/"+"QG.grb",
1019	#date=qdate, time="00",expver=self.expver, param="133.128")
1020	#QG.displayInfo()
1021	#QG.dataRetrieve()
1022
1023	oro=False
1024	for ftype in self.types:
1025	for pk,pv in self.params.iteritems():
1026	if isinstance(pv,str): # or pk!='GG__SL' :
1027	continue
1028	mftype=''+ftype
1029	mftime=self.types[ftype]['times']
1030	mfstep=self.types[ftype]['steps']
1031	mfdate=self.dates
1032	mfstream=self.stream
1033	mftarget=self.inputdir+"/"+ftype+pk+'.'+self.dates.split('/')[0]+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
1034
1035	if pk=='OG__SL':
1036	pass
1037	if pk=='OG_OROLSM__SL':
1038	if oro==False:
1039	# in CERA20C there is no stream "OPER"!
1040	if self.marsclass.upper() != 'EP':
1041	mfstream='OPER'
1042	mftype='AN'
1043	mftime='00'
1044	mfstep='000'
1045	mfdate=self.dates.split('/')[0]
1046	mftarget=self.inputdir+"/"+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
1047	oro=True
1048	else:
1049	continue
1050
1051	if pk=='GG__SL' and pv[0]=='Q':
1052	area=""
1053	gaussian='reduced'
1054	else:
1055	area=self.area
1056	gaussian=self.gaussian
1057
1058	if self.basetime==None:
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=mfstream,
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.displayInfo()
1068	MR.dataRetrieve()
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.displayInfo()
1074	MR.dataRetrieve()
1075	# The whole else section is only necessary for operational scripts. It could be removed
1076	else: # check if mars job requests fields beyond basetime. If yes eliminate those fields since they may not
1077	# be accessible with user's credentials
1078	sm1=-1
1079	if 'by' in mfstep:
1080	sm1=2
1081	tm1=-1
1082	if 'by' in mftime:
1083	tm1=2
1084	maxtime=datetime.datetime.strptime(mfdate.split('/')[-1]+mftime.split('/')[tm1],'%Y%m%d%H')+ \
1085	datetime.timedelta(hours=int(mfstep.split('/')[sm1]))
1086
1087	elimit=datetime.datetime.strptime(mfdate.split('/')[-1]+self.basetime,'%Y%m%d%H')
1088
1089	if self.basetime=='12':
1090	if 'acc' in pk:
1091
1092	# Strategy: if maxtime-elimit>=24h reduce date by 1,
1093	# if 12h<=maxtime-elimit<12h reduce time for last date
1094	# if maxtime-elimit<12h reduce step for last time
1095	# A split of the MARS job into 2 is likely necessary.
1096	maxtime=elimit-datetime.timedelta(hours=24)
1097	mfdate='/'.join(('/'.join(mfdate.split('/')[:-1]),datetime.datetime.strftime(maxtime,'%Y%m%d')))
1098	# increase number of mars requests
1099	self.mreq_count += 1
1100	MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
1101	type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1102	accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1103	date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1104
1105	if request == 0:
1106	MR.displayInfo()
1107	MR.dataRetrieve()
1108	elif request == 1:
1109	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1110	elif request == 2:
1111	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1112	MR.displayInfo()
1113	MR.dataRetrieve()
1114
1115	maxtime=elimit-datetime.timedelta(hours=12)
1116	mfdate=datetime.datetime.strftime(maxtime,'%Y%m%d')
1117	mftime='00'
1118	mftarget=self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
1119	# increase number of mars requests
1120	self.mreq_count += 1
1121	MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
1122	type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1123	accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1124	date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1125	if request == 0:
1126	MR.displayInfo()
1127	MR.dataRetrieve()
1128	elif request == 1:
1129	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1130	elif request == 2:
1131	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1132	MR.displayInfo()
1133	MR.dataRetrieve()
1134	else:
1135	# increase number of mars requests
1136	self.mreq_count += 1
1137	MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
1138	type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1139	accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1140	date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1141
1142	if request == 0:
1143	MR.displayInfo()
1144	MR.dataRetrieve()
1145	elif request == 1:
1146	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1147	elif request == 2:
1148	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1149	MR.displayInfo()
1150	MR.dataRetrieve()
1151	else:
1152	maxtime=elimit-datetime.timedelta(hours=24)
1153	mfdate=datetime.datetime.strftime(maxtime,'%Y%m%d')
1154
1155	mftimesave=''.join(mftime)
1156
1157	if '/' in mftime:
1158	times=mftime.split('/')
1159	while int(times[0])+int(mfstep.split('/')[0])<=12 and pk!='OG_OROLSM__SL' and 'acc' not in pk:
1160	times=times[1:]
1161	if len(times)>1:
1162	mftime='/'.join(times)
1163	else:
1164	mftime=times[0]
1165	# increase number of mars requests
1166	self.mreq_count += 1
1167	MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
1168	type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1169	accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1170	date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1171
1172	if request == 0:
1173	MR.displayInfo()
1174	MR.dataRetrieve()
1175	elif request == 1:
1176	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1177	elif request == 2:
1178	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1179	MR.displayInfo()
1180	MR.dataRetrieve()
1181
1182	if int(mftimesave.split('/')[0])==0 and int(mfstep.split('/')[0])==0 and pk!='OG_OROLSM__SL':
1183	mfdate=datetime.datetime.strftime(elimit,'%Y%m%d')
1184	mftime='00'
1185	mfstep='000'
1186	mftarget=self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
1187	# increase number of mars requests
1188	self.mreq_count += 1
1189	MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
1190	type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1191	accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1192	date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1193
1194	if request == 0:
1195	MR.displayInfo()
1196	MR.dataRetrieve()
1197	elif request == 1:
1198	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1199	elif request == 2:
1200	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1201	MR.displayInfo()
1202	MR.dataRetrieve()
1203
1204
1205	if request == 0 or request == 2:
1206	print('MARS retrieve done ... ')
1207	elif request == 1:
1208	print('MARS request printed ...')
1209
1210	def getFlexpartTime(self, type,step, time):
1211	cstep='{:0>3}'.format(step)
1212	ctime='{:0>2}'.format(int(time/100))
1213
1214	ctype = str(type).upper()
1215	myinfo = [ctype,ctime, cstep]
1216	cflextime = None
1217	for t, marsinfo in self.mars.items():
1218	if myinfo == marsinfo:
1219	cflextime=t
1220	return cflextime
1221
1222	def process_output(self, c):
1223
1224	print 'Postprocessing:\n Format: {}\n'.format(c.format)
1225	if c.ecapi==False:
1226	print 'ecstorage: {}\n ecfsdir: {}\n'.format(c.ecstorage,c.ecfsdir)
1227	if not hasattr(c,'gateway'):
1228	c.gateway=os.getenv('GATEWAY')
1229	if not hasattr(c,'destination'):
1230	c.destination=os.getenv('DESTINATION')
1231	print 'ectrans: {}\n gateway: {}\n destination: {}\n '.format(c.ectrans, c.gateway,c.destination)
1232	print 'Output filelist: \n',self.outputfilelist
1233
1234	if c.format.lower()=='grib2':
1235	for ofile in self.outputfilelist:
1236	p=subprocess.check_call(['grib_set','-s','edition=2,productDefinitionTemplateNumber=8',ofile,ofile+'_2'])
1237	p=subprocess.check_call(['mv',ofile+'_2',ofile])
1238	if c.debug==0:
1239	inputfiles=glob.glob('*.grb')
1240	for grb in inputfiles:
1241	try:
1242	os.remove(grb)
1243	except:
1244	pass
1245	if c.stream=='ELDA':
1246	opposite(self.inputdir+'/'+c.prefix)
1247	for i in range(len(self.outputfilelist)):
1248	if self.outputfilelist[i][-4:]!='N000' :
1249	j=int(self.outputfilelist[i][-3:])
1250	self.outputfilelist.append(self.outputfilelist[i][:-3]+'{:0>3}'.format(j+25))
1251
1252	if int(c.ectrans)==1 and c.ecapi==False:
1253	for ofile in self.outputfilelist:
1254	p=subprocess.check_call(['ectrans','-overwrite','-gateway',c.gateway,
1255	'-remote',c.destination,'-source',ofile])
1256	print 'ectrans:',p
1257	if int(c.ecstorage)==1 and c.ecapi==False:
1258	for ofile in self.outputfilelist:
1259	p=subprocess.check_call(['ecp','-o',ofile,os.path.expandvars(c.ecfsdir)])
1260
1261	#20131107 000000 EN13110700 ON DISC
1262	if c.outputdir!=c.inputdir:
1263	for ofile in self.outputfilelist:
1264	p=subprocess.check_call(['mv',ofile,c.outputdir])
1265
1266	if c.grib2flexpart=='1':
1267	f=open(c.outputdir+'/'+'AVAILABLE','w')
1268	clist=[]
1269	for ofile in self.outputfilelist: # generate AVAILABLE file
1270	fname=ofile.split('/')
1271	if '.' in fname[-1]:
1272	l=fname[-1].split('.')
1273	timestamp=datetime.datetime.strptime(l[0][-6:]+l[1],'%y%m%d%H')
1274	timestamp+=datetime.timedelta(hours=int(l[2]))
1275	cdate=datetime.datetime.strftime(timestamp,'%Y%m%d')
1276	chms=datetime.datetime.strftime(timestamp,'%H%M%S')
1277
1278	else:
1279	cdate='20'+fname[-1][-8:-2]
1280	chms=fname[-1][-2:]+'0000'
1281	clist.append(cdate+' '+chms+' '6+fname[-1]+' '14+'ON DISC')
1282	clist.sort()
1283	f.write('\n'.join(clist)+'\n')
1284	f.close()
1285
1286	pwd=os.path.abspath(c.outputdir)
1287	f=open(pwd+'/pathnames','w')
1288	f.write(pwd+'/Options/\n')
1289	f.write(pwd+'/\n')
1290	f.write(pwd+'/\n')
1291	f.write(pwd+'/AVAILABLE\n')
1292	f.write('==================\n')
1293	f.close()
1294
1295	if not os.path.exists(pwd+'/Options'):
1296	os.makedirs(pwd+'/Options')
1297	f=open(os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))+'/../Options/COMMAND','r')
1298	lflist=f.read().split('\n')
1299	i=0
1300	for l in lflist:
1301	if 'LDIRECT' in l.upper():
1302	break
1303	i+=1
1304
1305	clist.sort()
1306	lflist=lflist[:i+1]+[clist[0][:16],clist[-1][:16]]+lflist[i+3:]
1307	g=open(pwd+'/Options/COMMAND','w')
1308	g.write('\n'.join(lflist)+'\n')
1309	g.close()
1310
1311	os.chdir(c.outputdir)
1312	p=subprocess.check_call([os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))+'/../FLEXPART_PROGRAM/grib2flexpart',
1313	'useAvailable','.'])
1314	os.chdir(pwd)
1315
1316	def create(self, inputfiles, c,):
1317
1318	table128=init128(c.flexextractdir+'/grib_templates/ecmwf_grib1_table_128')
1319	wrfpars=toparamId('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',table128)
1320	# index_keys=["date","time","stepRange"]
1321	if '/' in c.number:
1322	index_keys=["number","date","time","step"]
1323	else:
1324	index_keys=["date","time","step"]
1325
1326	indexfile=c.inputdir+"/date_time_stepRange.idx"
1327	silentremove(indexfile)
1328	grib=GribTools(inputfiles.files)
1329	iid=grib.index(index_keys=index_keys, index_file = indexfile)
1330
1331	print 'index done...'
1332	fdict={'10':None,'11':None,'12':None,'13':None,'16':None,'17':None,'19':None,'21':None,'22':None,'20':None}
1333	for f in fdict.keys():
1334	silentremove(c.inputdir+"/fort."+f)
1335
1336	index_vals = []
1337	for key in index_keys:
1338	key_vals = grib_index_get(iid,key)
1339	print key_vals
1340
1341	index_vals.append(key_vals)
1342
1343	for prod in product(*index_vals):
1344	for i in range(len(index_keys)):
1345	grib_index_select(iid,index_keys[i],prod[i])
1346
1347
1348	gid = grib_new_from_index(iid)
1349	hid = gid
1350	cflextime = None
1351	for k,f in fdict.iteritems():
1352	fdict[k] = open(c.inputdir+'/fort.'+k,'w')
1353	if gid is not None:
1354	cdate = str(grib_get(gid, 'date'))
1355	time = grib_get(gid, 'time')
1356	type = grib_get(gid, 'type')
1357	step = grib_get(gid, 'step')
1358	# step = grib_get(gid, 'stepRange')
1359	# cflextime = self.getFlexpartTime(type,step, time)
1360	timestamp=datetime.datetime.strptime(cdate+'{:0>2}'.format(time/100),'%Y%m%d%H')
1361	timestamp+=datetime.timedelta(hours=int(step))
1362	# print gid,index_keys[i],prod[i],cdate,time,step,timestamp
1363
1364	cdateH=datetime.datetime.strftime(timestamp,'%Y%m%d%H')
1365	chms=datetime.datetime.strftime(timestamp,'%H%M%S')
1366
1367	if c.basetime !=None:
1368	slimit=datetime.datetime.strptime(c.start_date+'00','%Y%m%d%H')
1369	bt='23'
1370	if c.basetime=='00':
1371	bt='00'
1372	slimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')-datetime.timedelta(hours=12-int(c.dtime))
1373
1374	if c.basetime=='12':
1375	bt='12'
1376	slimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')-datetime.timedelta(hours=12-int(c.dtime))
1377
1378	elimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')
1379
1380	if timestamp<slimit or timestamp>elimit:
1381	continue
1382	else:
1383	if c.maxstep<24:
1384	if cdateH<c.start_date+'00':
1385	continue
1386	if cdateH>c.end_date+'23':
1387	continue
1388
1389
1390
1391	try:
1392	if c.wrf=='1':
1393	if 'olddate' not in locals():
1394	fwrf = open(c.outputdir+'/WRF'+cdate+'.{:0>2}'.format(time)+'.000.grb2','w')
1395	olddate=cdate[:]
1396	else:
1397	if cdate!=olddate:
1398	fwrf = open(c.outputdir+'/WRF'+cdate+'.{:0>2}'.format(time)+'.000.grb2','w')
1399	olddate=cdate[:]
1400	except AttributeError:
1401	pass
1402
1403	# print 'cyear '+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime
1404
1405	savedfields=[]
1406	while 1:
1407	if gid is None: break
1408	paramId = grib_get(gid, 'paramId')
1409	gridtype = grib_get(gid, 'gridType')
1410	datatype = grib_get(gid, 'dataType')
1411	levtype = grib_get(gid, 'typeOfLevel')
1412	if paramId == 133 and gridtype=='reduced_gg':
1413	# Relative humidity (Q.grb) is used as a template only so we need the first we "meet"
1414	fout=open(c.inputdir+'/fort.18','w')
1415	grib_write(gid,fout)
1416	fout.close()
1417	elif paramId == 131 or paramId == 132:
1418	grib_write(gid, fdict['10'])
1419	elif paramId == 130:
1420	grib_write(gid, fdict['11'])
1421	elif paramId == 133 and gridtype!='reduced_gg':
1422	grib_write(gid, fdict['17'])
1423	elif paramId == 152:
1424	grib_write(gid, fdict['12'])
1425	elif paramId == 155 and gridtype=='sh':
1426	grib_write(gid, fdict['13'])
1427	elif paramId in [129,138,155] and levtype=='hybrid' and c.wrf=='1':
1428	# print paramId,'not written'
1429	pass
1430	elif paramId == 246 or paramId == 247: # cloud liquid water and ice
1431	if paramId==246:
1432	clwc=grib_get_values(gid)
1433	else:
1434	clwc+=grib_get_values(gid)
1435	grib_set_values(gid,clwc)
1436	# grib_set(gid,'shortName','qc')
1437	grib_set(gid,'paramId',201031)
1438	grib_write(gid, fdict['22'])
1439
1440	elif paramId == 135:
1441	grib_write(gid, fdict['19'])
1442	elif paramId == 77:
1443	grib_write(gid, fdict['21'])
1444	else:
1445	if paramId not in savedfields:
1446	grib_write(gid, fdict['16'])
1447	savedfields.append(paramId)
1448	else:
1449	print 'duplicate '+str(paramId)+' not written'
1450
1451	try:
1452	if c.wrf=='1':
1453	if levtype=='hybrid':
1454	if paramId in [129,130,131,132,133,138,155]:
1455	grib_write(gid,fwrf)
1456	else:
1457	if paramId in wrfpars:
1458	grib_write(gid,fwrf)
1459	except AttributeError:
1460	pass
1461
1462
1463
1464	grib_release(gid)
1465	gid = grib_new_from_index(iid)
1466	# call for CONVERT2
1467
1468	for f in fdict.values():
1469	f.close()
1470
1471	if hid is not None:
1472	pwd=os.getcwd()
1473	os.chdir(c.inputdir)
1474	if os.stat('fort.21').st_size==0 and int(c.eta)==1:
1475	print 'Parameter 77 (etadot) is missing, most likely it is not available for this type or date/time\n'
1476	print 'Check parameters CLASS, TYPE, STREAM, START_DATE\n'
1477	myerror(c,'fort.21 is empty while parameter eta is set to 1 in CONTROL file')
1478
1479	p=subprocess.check_call([os.path.expandvars(os.path.expanduser(c.exedir))+'/CONVERT2'],shell=True)
1480	os.chdir(pwd)
1481	# create the corresponding output file fort.15 (generated by CONVERT2) + fort.16 (paramId 167 and paramId 168)
1482	fnout=c.inputdir+'/'+c.prefix
1483	if c.maxstep>12:
1484	suffix=cdate[2:8]+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step)
1485	else:
1486	suffix=cdateH[2:10]
1487	try:
1488	numberindex=index_keys.index('number')
1489	if len(index_vals[numberindex])>1:
1490	suffix=suffix+'.N{:0>3}'.format(int(prod[numberindex]))
1491	except:
1492	pass
1493
1494	fnout+=suffix
1495	print "outputfile = " + fnout
1496	self.outputfilelist.append(fnout) # needed for final processing
1497	fout = open(fnout,'wb')
1498	shutil.copyfileobj(open(c.inputdir+'/fort.15','rb'), fout)
1499	if c.cwc=='1':
1500	shutil.copyfileobj(open(c.inputdir+'/fort.22','rb'), fout)
1501	shutil.copyfileobj(open(c.inputdir+'/flux'+cdate[0:2]+suffix,'rb'), fout)
1502	shutil.copyfileobj(open(c.inputdir+'/fort.16','rb'), fout)
1503	orolsm=glob.glob(c.inputdir+'/OG_OROLSM__SL..'+c.ppid+'')[0]
1504	shutil.copyfileobj(open(orolsm,'rb'), fout)
1505	fout.close()
1506	if c.omega=='1':
1507	fnout=c.outputdir+'/OMEGA'+suffix
1508	fout = open(fnout,'wb')
1509	shutil.copyfileobj(open(c.inputdir+'/fort.25','rb'), fout)
1510
1511
1512
1513	try:
1514	if c.wrf=='1':
1515	fwrf.close()
1516	except:
1517	pass
1518
1519	grib_index_release(iid)
1520
1521	return
1522
1523	def deacc_fluxes(self, inputfiles, c):
1524
1525	table128=init128(c.flexextractdir+'/grib_templates/ecmwf_grib1_table_128')
1526	pars=toparamId(self.params['OG_acc_SL'][0],table128)
1527	if '/' in c.number:
1528	index_keys=["number","date","time","step"]
1529	else:
1530	index_keys=["date","time","step"]
1531	indexfile=c.inputdir+"/date_time_stepRange.idx"
1532	silentremove(indexfile)
1533	grib=GribTools(inputfiles.files)
1534	iid=grib.index(index_keys=index_keys, index_file = indexfile)
1535
1536	print 'index done...'
1537
1538	index_vals = []
1539	for key in index_keys:
1540	key_vals = grib_index_get(iid,key)
1541	print(key_vals)
1542	l=[]
1543	for k in key_vals:
1544	l.append(int(k))
1545	l.sort()
1546	key_vals=[]
1547	for k in l:
1548	key_vals.append(str(k))
1549	print key_vals
1550	index_vals.append(key_vals)
1551
1552
1553	valsdict={}
1554	svalsdict={}
1555	stepsdict={}
1556	for p in pars:
1557	valsdict[str(p)]=[]
1558	svalsdict[str(p)]=[]
1559	stepsdict[str(p)]=[]
1560
1561	for prod in product(*index_vals):
1562	for i in range(len(index_keys)):
1563	grib_index_select(iid,index_keys[i],prod[i])
1564
1565	#for k,f in fdict.iteritems():
1566	#fdict[k] = open('fort.'+k,'w')
1567	gid = grib_new_from_index(iid)
1568	hid = gid
1569	cflextime = None
1570	if gid is not None:
1571	cdate = grib_get(gid, 'date')
1572	#cyear = cdate[:4]
1573	#cmonth = cdate[4:6]
1574	#cday = cdate[6:8]
1575	time = grib_get(gid, 'time')
1576	type = grib_get(gid, 'type')
1577	step = grib_get(gid, 'step')
1578	# date+time+step-2dtime (since interpolated value valid for step-2dtime)
1579	sdate=datetime.datetime(year=cdate/10000,month=mod(cdate,10000)/100,
1580	day=mod(cdate,100),hour=time/100)
1581	fdate=sdate+datetime.timedelta(hours=step-2*int(c.dtime))
1582	sdates=sdate+datetime.timedelta(hours=step)
1583	else:
1584	break
1585
1586	fnout=c.inputdir+'/'
1587	numbersuffix=''
1588	try:
1589	numberindex=index_keys.index('number')
1590	if len(index_vals[numberindex])>1:
1591	numbersuffix='.N{:0>3}'.format(int(prod[numberindex]))
1592	except:
1593	pass
1594
1595	if c.maxstep>12:
1596	fnout+='flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step-2*int(c.dtime))+numbersuffix
1597	gnout=c.inputdir+'/flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step-int(c.dtime))+numbersuffix
1598	hnout=c.inputdir+'/flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step)+numbersuffix
1599	g=open(gnout,'w')
1600	h=open(hnout,'w')
1601	else:
1602	fnout+='flux'+fdate.strftime('%Y%m%d%H')+numbersuffix
1603	gnout=c.inputdir+'/flux'+(fdate+datetime.timedelta(hours=int(c.dtime))).strftime('%Y%m%d%H')+numbersuffix
1604	hnout=c.inputdir+'/flux'+sdates.strftime('%Y%m%d%H')+numbersuffix
1605	g=open(gnout,'w')
1606	h=open(hnout,'w')
1607	print "outputfile = " + fnout
1608	f=open(fnout,'w')
1609
1610	while 1:
1611	if gid is None: break
1612	cparamId = str(grib_get(gid, 'paramId'))
1613	step = grib_get(gid, 'step')
1614	atime = grib_get(gid, 'time')
1615	ni=grib_get(gid, 'Ni')
1616	nj=grib_get(gid, 'Nj')
1617	if cparamId in valsdict.keys():
1618	values = grib_get_values(gid)
1619	vdp=valsdict[cparamId]
1620	svdp=svalsdict[cparamId]
1621	sd=stepsdict[cparamId]
1622
1623	if cparamId=='142' or cparamId=='143':
1624	fak=1./1000.
1625	else:
1626	fak=3600.
1627
1628	values=(reshape(values,(nj,ni))).flatten()/fak
1629	vdp.append(values[:]) # save the accumulated values
1630	if c.marsclass.upper() in ('EA') or \
1631	step<=int(c.dtime):
1632	svdp.append(values[:]/int(c.dtime))
1633	else:
1634	svdp.append((vdp[-1]-vdp[-2])/int(c.dtime))
1635
1636	print cparamId,atime,step,len(values),values[0],std(values)
1637	#svdp.append(values[:]) # save the 1-hourly or 3-hourly specific values
1638	sd.append(step)
1639	if len(svdp)>=3:
1640	if len(svdp)>3:
1641	if cparamId=='142' or cparamId=='143':
1642	values=darain(svdp)
1643	else:
1644	values=dapoly(svdp)
1645
1646	if not (step==c.maxstep and c.maxstep>12 or sdates==elimit):
1647	vdp.pop(0)
1648	svdp.pop(0)
1649	else:
1650	if c.maxstep>12:
1651	values=svdp[1]
1652	else:
1653	values=svdp[0]
1654
1655	grib_set_values(gid, values)
1656	if c.maxstep>12:
1657	grib_set(gid,'stepRange',max(0,step-2*int(c.dtime)))
1658	else:
1659	grib_set(gid,'stepRange',0)
1660	grib_set(gid,'time',fdate.hour*100)
1661	grib_set(gid,'date',fdate.year10000+fdate.month100+fdate.day)
1662	grib_write(gid, f)
1663
1664	if c.basetime is not None:
1665	elimit=datetime.datetime.strptime(c.end_date+c.basetime,'%Y%m%d%H')
1666	else:
1667	elimit=sdate+datetime.timedelta(2*int(c.dtime))
1668
1669	# squeeze out information of last two steps contained in svdp
1670	# if step+int(c.dtime)==c.maxstep and c.maxstep>12 or sdates+datetime.timedelta(hours=int(c.dtime))>=elimit:
1671	# Note that svdp[0] has not been popped in this case
1672
1673
1674	if step==c.maxstep and c.maxstep>12 or sdates==elimit:
1675	values=svdp[3]
1676	grib_set_values(gid, values)
1677	grib_set(gid,'stepRange',0)
1678	truedatetime=fdate+datetime.timedelta(hours=2*int(c.dtime))
1679	grib_set(gid,'time',truedatetime.hour*100)
1680	grib_set(gid,'date',truedatetime.year10000+truedatetime.month100+truedatetime.day)
1681	grib_write(gid, h)
1682
1683	#values=(svdp[1]+svdp[2])/2.
1684	if cparamId=='142' or cparamId=='143':
1685	values=darain(list(reversed(svdp)))
1686	else:
1687	values=dapoly(list(reversed(svdp)))
1688
1689	grib_set(gid,'stepRange',0)
1690	truedatetime=fdate+datetime.timedelta(hours=int(c.dtime))
1691	grib_set(gid,'time',truedatetime.hour*100)
1692	grib_set(gid,'date',truedatetime.year10000+truedatetime.month100+truedatetime.day)
1693	grib_set_values(gid, values)
1694	grib_write(gid, g)
1695
1696	grib_release(gid)
1697
1698	gid = grib_new_from_index(iid)
1699
1700	f.close()
1701	g.close()
1702	h.close()
1703
1704
1705	grib_index_release(iid)

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

Download in other formats: