Context Navigation

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

ctbtodev

Last change on this file since e879023 was e879023, checked in by anphi <anne.philipp@…>, 5 years ago
fixed bug in retrieving with basetime
Property mode set to `100644`
File size: 68.5 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(len(c.time))
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	(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 pk=='OG_OROLSM__SL':
1158	mfdate=self.dates.split('/')[0]
1159	mftarget=self.inputdir+"/"+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
1160
1161	if '/' in mftime:
1162	times=mftime.split('/')
1163	while int(times[0])+int(mfstep.split('/')[0])<=12 and pk!='OG_OROLSM__SL' and 'acc' not in pk:
1164	times=times[1:]
1165	if len(times)>1:
1166	mftime='/'.join(times)
1167	else:
1168	mftime=times[0]
1169
1170	if int(mftimesave.split('/')[0])==0 and int(mfstep.split('/')[0])==0 and pk!='OG_OROLSM__SL':
1171	mfdate=datetime.datetime.strftime(elimit,'%Y%m%d')
1172	mftime='00'
1173	mfstep='000'
1174	mftarget=self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb"
1175
1176	# increase number of mars requests
1177	self.mreq_count += 1
1178	MR= MARSretrieval(self.server, self.public, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream,
1179	type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian,
1180	accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area,
1181	date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0])
1182
1183	if request == 0:
1184	MR.displayInfo()
1185	MR.dataRetrieve()
1186	elif request == 1:
1187	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1188	elif request == 2:
1189	MR.print_infodata_csv(self.inputdir, self.mreq_count)
1190	MR.displayInfo()
1191	MR.dataRetrieve()
1192
1193
1194	if request == 0 or request == 2:
1195	print('MARS retrieve done ... ')
1196	elif request == 1:
1197	print('MARS request printed ...')
1198
1199	def getFlexpartTime(self, type,step, time):
1200	cstep='{:0>3}'.format(step)
1201	ctime='{:0>2}'.format(int(time/100))
1202
1203	ctype = str(type).upper()
1204	myinfo = [ctype,ctime, cstep]
1205	cflextime = None
1206	for t, marsinfo in self.mars.items():
1207	if myinfo == marsinfo:
1208	cflextime=t
1209	return cflextime
1210
1211	def process_output(self, c):
1212
1213	print 'Postprocessing:\n Format: {}\n'.format(c.format)
1214	if c.ecapi==False:
1215	print 'ecstorage: {}\n ecfsdir: {}\n'.format(c.ecstorage,c.ecfsdir)
1216	if not hasattr(c,'gateway'):
1217	c.gateway=os.getenv('GATEWAY')
1218	if not hasattr(c,'destination'):
1219	c.destination=os.getenv('DESTINATION')
1220	print 'ectrans: {}\n gateway: {}\n destination: {}\n '.format(c.ectrans, c.gateway,c.destination)
1221	print 'Output filelist: \n',self.outputfilelist
1222
1223	if c.format.lower()=='grib2':
1224	for ofile in self.outputfilelist:
1225	p=subprocess.check_call(['grib_set','-s','edition=2,productDefinitionTemplateNumber=8',ofile,ofile+'_2'])
1226	p=subprocess.check_call(['mv',ofile+'_2',ofile])
1227	if c.debug==0:
1228	inputfiles=glob.glob('*.grb')
1229	for grb in inputfiles:
1230	try:
1231	os.remove(grb)
1232	except:
1233	pass
1234	if c.stream=='ELDA':
1235	opposite(self.inputdir+'/'+c.prefix)
1236	for i in range(len(self.outputfilelist)):
1237	if self.outputfilelist[i][-4:]!='N000' :
1238	j=int(self.outputfilelist[i][-3:])
1239	self.outputfilelist.append(self.outputfilelist[i][:-3]+'{:0>3}'.format(j+25))
1240
1241	if int(c.ectrans)==1 and c.ecapi==False:
1242	for ofile in self.outputfilelist:
1243	p=subprocess.check_call(['ectrans','-overwrite','-gateway',c.gateway,
1244	'-remote',c.destination,'-source',ofile])
1245	print 'ectrans:',p
1246	if int(c.ecstorage)==1 and c.ecapi==False:
1247	for ofile in self.outputfilelist:
1248	p=subprocess.check_call(['ecp','-o',ofile,os.path.expandvars(c.ecfsdir)])
1249
1250	#20131107 000000 EN13110700 ON DISC
1251	if c.outputdir!=c.inputdir:
1252	for ofile in self.outputfilelist:
1253	p=subprocess.check_call(['mv',ofile,c.outputdir])
1254
1255	if c.grib2flexpart=='1':
1256	f=open(c.outputdir+'/'+'AVAILABLE','w')
1257	clist=[]
1258	for ofile in self.outputfilelist: # generate AVAILABLE file
1259	fname=ofile.split('/')
1260	if '.' in fname[-1]:
1261	l=fname[-1].split('.')
1262	timestamp=datetime.datetime.strptime(l[0][-6:]+l[1],'%y%m%d%H')
1263	timestamp+=datetime.timedelta(hours=int(l[2]))
1264	cdate=datetime.datetime.strftime(timestamp,'%Y%m%d')
1265	chms=datetime.datetime.strftime(timestamp,'%H%M%S')
1266
1267	else:
1268	cdate='20'+fname[-1][-8:-2]
1269	chms=fname[-1][-2:]+'0000'
1270	clist.append(cdate+' '+chms+' '6+fname[-1]+' '14+'ON DISC')
1271	clist.sort()
1272	f.write('\n'.join(clist)+'\n')
1273	f.close()
1274
1275	pwd=os.path.abspath(c.outputdir)
1276	f=open(pwd+'/pathnames','w')
1277	f.write(pwd+'/Options/\n')
1278	f.write(pwd+'/\n')
1279	f.write(pwd+'/\n')
1280	f.write(pwd+'/AVAILABLE\n')
1281	f.write('==================\n')
1282	f.close()
1283
1284	if not os.path.exists(pwd+'/Options'):
1285	os.makedirs(pwd+'/Options')
1286	f=open(os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))+'/../Options/COMMAND','r')
1287	lflist=f.read().split('\n')
1288	i=0
1289	for l in lflist:
1290	if 'LDIRECT' in l.upper():
1291	break
1292	i+=1
1293
1294	clist.sort()
1295	lflist=lflist[:i+1]+[clist[0][:16],clist[-1][:16]]+lflist[i+3:]
1296	g=open(pwd+'/Options/COMMAND','w')
1297	g.write('\n'.join(lflist)+'\n')
1298	g.close()
1299
1300	os.chdir(c.outputdir)
1301	p=subprocess.check_call([os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))+'/../FLEXPART_PROGRAM/grib2flexpart',
1302	'useAvailable','.'])
1303	os.chdir(pwd)
1304
1305	def create(self, inputfiles, c,):
1306
1307	table128=init128(c.flexextractdir+'/grib_templates/ecmwf_grib1_table_128')
1308	wrfpars=toparamId('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',table128)
1309	# index_keys=["date","time","stepRange"]
1310	if '/' in c.number:
1311	index_keys=["number","date","time","step"]
1312	else:
1313	index_keys=["date","time","step"]
1314
1315	indexfile=c.inputdir+"/date_time_stepRange.idx"
1316	silentremove(indexfile)
1317	grib=GribTools(inputfiles.files)
1318	iid=grib.index(index_keys=index_keys, index_file = indexfile)
1319
1320	print 'index done...'
1321	fdict={'10':None,'11':None,'12':None,'13':None,'16':None,'17':None,'19':None,'21':None,'22':None,'20':None}
1322	for f in fdict.keys():
1323	silentremove(c.inputdir+"/fort."+f)
1324
1325	index_vals = []
1326	for key in index_keys:
1327	key_vals = grib_index_get(iid,key)
1328	print key_vals
1329
1330	index_vals.append(key_vals)
1331
1332	for prod in product(*index_vals):
1333	for i in range(len(index_keys)):
1334	grib_index_select(iid,index_keys[i],prod[i])
1335
1336
1337	gid = grib_new_from_index(iid)
1338	hid = gid
1339	cflextime = None
1340	for k,f in fdict.iteritems():
1341	fdict[k] = open(c.inputdir+'/fort.'+k,'w')
1342	if gid is not None:
1343	cdate = str(grib_get(gid, 'date'))
1344	time = grib_get(gid, 'time')
1345	type = grib_get(gid, 'type')
1346	step = grib_get(gid, 'step')
1347	# step = grib_get(gid, 'stepRange')
1348	# cflextime = self.getFlexpartTime(type,step, time)
1349	timestamp=datetime.datetime.strptime(cdate+'{:0>2}'.format(time/100),'%Y%m%d%H')
1350	timestamp+=datetime.timedelta(hours=int(step))
1351	# print gid,index_keys[i],prod[i],cdate,time,step,timestamp
1352
1353	cdateH=datetime.datetime.strftime(timestamp,'%Y%m%d%H')
1354	chms=datetime.datetime.strftime(timestamp,'%H%M%S')
1355
1356	if c.basetime !=None:
1357	slimit=datetime.datetime.strptime(c.start_date+'00','%Y%m%d%H')
1358	bt='23'
1359	if c.basetime=='00':
1360	bt='00'
1361	slimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')-datetime.timedelta(hours=12-int(c.dtime))
1362
1363	if c.basetime=='12':
1364	bt='12'
1365	slimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')-datetime.timedelta(hours=12-int(c.dtime))
1366
1367	elimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')
1368
1369	if timestamp<slimit or timestamp>elimit:
1370	continue
1371	else:
1372	if c.maxstep<24:
1373	if cdateH<c.start_date+'00':
1374	continue
1375	if cdateH>c.end_date+'23':
1376	continue
1377
1378
1379
1380	try:
1381	if c.wrf=='1':
1382	if 'olddate' not in locals():
1383	fwrf = open(c.outputdir+'/WRF'+cdate+'.{:0>2}'.format(time)+'.000.grb2','w')
1384	olddate=cdate[:]
1385	else:
1386	if cdate!=olddate:
1387	fwrf = open(c.outputdir+'/WRF'+cdate+'.{:0>2}'.format(time)+'.000.grb2','w')
1388	olddate=cdate[:]
1389	except AttributeError:
1390	pass
1391
1392	# print 'cyear '+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime
1393
1394	savedfields=[]
1395	while 1:
1396	if gid is None: break
1397	paramId = grib_get(gid, 'paramId')
1398	gridtype = grib_get(gid, 'gridType')
1399	datatype = grib_get(gid, 'dataType')
1400	levtype = grib_get(gid, 'typeOfLevel')
1401	if paramId == 133 and gridtype=='reduced_gg':
1402	# Relative humidity (Q.grb) is used as a template only so we need the first we "meet"
1403	fout=open(c.inputdir+'/fort.18','w')
1404	grib_write(gid,fout)
1405	fout.close()
1406	elif paramId == 131 or paramId == 132:
1407	grib_write(gid, fdict['10'])
1408	elif paramId == 130:
1409	grib_write(gid, fdict['11'])
1410	elif paramId == 133 and gridtype!='reduced_gg':
1411	grib_write(gid, fdict['17'])
1412	elif paramId == 152:
1413	grib_write(gid, fdict['12'])
1414	elif paramId == 155 and gridtype=='sh':
1415	grib_write(gid, fdict['13'])
1416	elif paramId in [129,138,155] and levtype=='hybrid' and c.wrf=='1':
1417	# print paramId,'not written'
1418	pass
1419	elif paramId == 246 or paramId == 247: # cloud liquid water and ice
1420	if paramId==246:
1421	clwc=grib_get_values(gid)
1422	else:
1423	clwc+=grib_get_values(gid)
1424	grib_set_values(gid,clwc)
1425	# grib_set(gid,'shortName','qc')
1426	grib_set(gid,'paramId',201031)
1427	grib_write(gid, fdict['22'])
1428
1429	elif paramId == 135:
1430	grib_write(gid, fdict['19'])
1431	elif paramId == 77:
1432	grib_write(gid, fdict['21'])
1433	else:
1434	if paramId not in savedfields:
1435	grib_write(gid, fdict['16'])
1436	savedfields.append(paramId)
1437	else:
1438	print 'duplicate '+str(paramId)+' not written'
1439
1440	try:
1441	if c.wrf=='1':
1442	if levtype=='hybrid':
1443	if paramId in [129,130,131,132,133,138,155]:
1444	grib_write(gid,fwrf)
1445	else:
1446	if paramId in wrfpars:
1447	grib_write(gid,fwrf)
1448	except AttributeError:
1449	pass
1450
1451
1452
1453	grib_release(gid)
1454	gid = grib_new_from_index(iid)
1455	# call for CONVERT2
1456
1457	for f in fdict.values():
1458	f.close()
1459
1460	if hid is not None:
1461	pwd=os.getcwd()
1462	os.chdir(c.inputdir)
1463	if os.stat('fort.21').st_size==0 and int(c.eta)==1:
1464	print 'Parameter 77 (etadot) is missing, most likely it is not available for this type or date/time\n'
1465	print 'Check parameters CLASS, TYPE, STREAM, START_DATE\n'
1466	myerror(c,'fort.21 is empty while parameter eta is set to 1 in CONTROL file')
1467
1468	p=subprocess.check_call([os.path.expandvars(os.path.expanduser(c.exedir))+'/CONVERT2'],shell=True)
1469	os.chdir(pwd)
1470	# create the corresponding output file fort.15 (generated by CONVERT2) + fort.16 (paramId 167 and paramId 168)
1471	fnout=c.inputdir+'/'+c.prefix
1472	if c.maxstep>12:
1473	suffix=cdate[2:8]+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step)
1474	else:
1475	suffix=cdateH[2:10]
1476	try:
1477	numberindex=index_keys.index('number')
1478	if len(index_vals[numberindex])>1:
1479	suffix=suffix+'.N{:0>3}'.format(int(prod[numberindex]))
1480	except:
1481	pass
1482
1483	fnout+=suffix
1484	print "outputfile = " + fnout
1485	self.outputfilelist.append(fnout) # needed for final processing
1486	fout = open(fnout,'wb')
1487	shutil.copyfileobj(open(c.inputdir+'/fort.15','rb'), fout)
1488	if c.cwc=='1':
1489	shutil.copyfileobj(open(c.inputdir+'/fort.22','rb'), fout)
1490	shutil.copyfileobj(open(c.inputdir+'/flux'+cdate[0:2]+suffix,'rb'), fout)
1491	shutil.copyfileobj(open(c.inputdir+'/fort.16','rb'), fout)
1492	orolsm=glob.glob(c.inputdir+'/OG_OROLSM__SL..'+c.ppid+'')[0]
1493	shutil.copyfileobj(open(orolsm,'rb'), fout)
1494	fout.close()
1495	if c.omega=='1':
1496	fnout=c.outputdir+'/OMEGA'+suffix
1497	fout = open(fnout,'wb')
1498	shutil.copyfileobj(open(c.inputdir+'/fort.25','rb'), fout)
1499
1500
1501
1502	try:
1503	if c.wrf=='1':
1504	fwrf.close()
1505	except:
1506	pass
1507
1508	grib_index_release(iid)
1509
1510	return
1511
1512	def deacc_fluxes(self, inputfiles, c):
1513
1514	table128=init128(c.flexextractdir+'/grib_templates/ecmwf_grib1_table_128')
1515	pars=toparamId(self.params['OG_acc_SL'][0],table128)
1516	if '/' in c.number:
1517	index_keys=["number","date","time","step"]
1518	else:
1519	index_keys=["date","time","step"]
1520	indexfile=c.inputdir+"/date_time_stepRange.idx"
1521	silentremove(indexfile)
1522	grib=GribTools(inputfiles.files)
1523	iid=grib.index(index_keys=index_keys, index_file = indexfile)
1524
1525	print 'index done...'
1526
1527	index_vals = []
1528	for key in index_keys:
1529	key_vals = grib_index_get(iid,key)
1530	print(key_vals)
1531	l=[]
1532	for k in key_vals:
1533	l.append(int(k))
1534	l.sort()
1535	key_vals=[]
1536	for k in l:
1537	key_vals.append(str(k))
1538	print key_vals
1539	index_vals.append(key_vals)
1540
1541
1542	valsdict={}
1543	svalsdict={}
1544	stepsdict={}
1545	for p in pars:
1546	valsdict[str(p)]=[]
1547	svalsdict[str(p)]=[]
1548	stepsdict[str(p)]=[]
1549
1550	for prod in product(*index_vals):
1551	for i in range(len(index_keys)):
1552	grib_index_select(iid,index_keys[i],prod[i])
1553
1554	#for k,f in fdict.iteritems():
1555	#fdict[k] = open('fort.'+k,'w')
1556	gid = grib_new_from_index(iid)
1557	hid = gid
1558	cflextime = None
1559	if gid is not None:
1560	cdate = grib_get(gid, 'date')
1561	#cyear = cdate[:4]
1562	#cmonth = cdate[4:6]
1563	#cday = cdate[6:8]
1564	time = grib_get(gid, 'time')
1565	type = grib_get(gid, 'type')
1566	step = grib_get(gid, 'step')
1567	# date+time+step-2dtime (since interpolated value valid for step-2dtime)
1568	sdate=datetime.datetime(year=cdate/10000,month=mod(cdate,10000)/100,
1569	day=mod(cdate,100),hour=time/100)
1570	fdate=sdate+datetime.timedelta(hours=step-2*int(c.dtime))
1571	sdates=sdate+datetime.timedelta(hours=step)
1572	else:
1573	break
1574
1575	fnout=c.inputdir+'/'
1576	numbersuffix=''
1577	try:
1578	numberindex=index_keys.index('number')
1579	if len(index_vals[numberindex])>1:
1580	numbersuffix='.N{:0>3}'.format(int(prod[numberindex]))
1581	except:
1582	pass
1583
1584	if c.maxstep>12:
1585	fnout+='flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step-2*int(c.dtime))+numbersuffix
1586	gnout=c.inputdir+'/flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step-int(c.dtime))+numbersuffix
1587	hnout=c.inputdir+'/flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step)+numbersuffix
1588	g=open(gnout,'w')
1589	h=open(hnout,'w')
1590	else:
1591	fnout+='flux'+fdate.strftime('%Y%m%d%H')+numbersuffix
1592	gnout=c.inputdir+'/flux'+(fdate+datetime.timedelta(hours=int(c.dtime))).strftime('%Y%m%d%H')+numbersuffix
1593	hnout=c.inputdir+'/flux'+sdates.strftime('%Y%m%d%H')+numbersuffix
1594	g=open(gnout,'w')
1595	h=open(hnout,'w')
1596	print "outputfile = " + fnout
1597	f=open(fnout,'w')
1598
1599	while 1:
1600	if gid is None: break
1601	cparamId = str(grib_get(gid, 'paramId'))
1602	step = grib_get(gid, 'step')
1603	atime = grib_get(gid, 'time')
1604	ni=grib_get(gid, 'Ni')
1605	nj=grib_get(gid, 'Nj')
1606	if cparamId in valsdict.keys():
1607	values = grib_get_values(gid)
1608	vdp=valsdict[cparamId]
1609	svdp=svalsdict[cparamId]
1610	sd=stepsdict[cparamId]
1611
1612	if cparamId=='142' or cparamId=='143':
1613	fak=1./1000.
1614	else:
1615	fak=3600.
1616
1617	values=(reshape(values,(nj,ni))).flatten()/fak
1618	vdp.append(values[:]) # save the accumulated values
1619	if c.marsclass.upper() in ('EA') or \
1620	step<=int(c.dtime):
1621	svdp.append(values[:]/int(c.dtime))
1622	else:
1623	svdp.append((vdp[-1]-vdp[-2])/int(c.dtime))
1624
1625	print cparamId,atime,step,len(values),values[0],std(values)
1626	#svdp.append(values[:]) # save the 1-hourly or 3-hourly specific values
1627	sd.append(step)
1628	if len(svdp)>=3:
1629	if len(svdp)>3:
1630	if cparamId=='142' or cparamId=='143':
1631	values=darain(svdp)
1632	else:
1633	values=dapoly(svdp)
1634
1635	if not (step==c.maxstep and c.maxstep>12 or sdates==elimit):
1636	vdp.pop(0)
1637	svdp.pop(0)
1638	else:
1639	if c.maxstep>12:
1640	values=svdp[1]
1641	else:
1642	values=svdp[0]
1643
1644	grib_set_values(gid, values)
1645	if c.maxstep>12:
1646	grib_set(gid,'stepRange',max(0,step-2*int(c.dtime)))
1647	else:
1648	grib_set(gid,'stepRange',0)
1649	grib_set(gid,'time',fdate.hour*100)
1650	grib_set(gid,'date',fdate.year10000+fdate.month100+fdate.day)
1651	grib_write(gid, f)
1652
1653	if c.basetime is not None:
1654	elimit=datetime.datetime.strptime(c.end_date+c.basetime,'%Y%m%d%H')
1655	else:
1656	elimit=sdate+datetime.timedelta(2*int(c.dtime))
1657
1658	# squeeze out information of last two steps contained in svdp
1659	# if step+int(c.dtime)==c.maxstep and c.maxstep>12 or sdates+datetime.timedelta(hours=int(c.dtime))>=elimit:
1660	# Note that svdp[0] has not been popped in this case
1661
1662
1663	if step==c.maxstep and c.maxstep>12 or sdates==elimit:
1664	values=svdp[3]
1665	grib_set_values(gid, values)
1666	grib_set(gid,'stepRange',0)
1667	truedatetime=fdate+datetime.timedelta(hours=2*int(c.dtime))
1668	grib_set(gid,'time',truedatetime.hour*100)
1669	grib_set(gid,'date',truedatetime.year10000+truedatetime.month100+truedatetime.day)
1670	grib_write(gid, h)
1671
1672	#values=(svdp[1]+svdp[2])/2.
1673	if cparamId=='142' or cparamId=='143':
1674	values=darain(list(reversed(svdp)))
1675	else:
1676	values=dapoly(list(reversed(svdp)))
1677
1678	grib_set(gid,'stepRange',0)
1679	truedatetime=fdate+datetime.timedelta(hours=int(c.dtime))
1680	grib_set(gid,'time',truedatetime.hour*100)
1681	grib_set(gid,'date',truedatetime.year10000+truedatetime.month100+truedatetime.day)
1682	grib_set_values(gid, values)
1683	grib_write(gid, g)
1684
1685	grib_release(gid)
1686
1687	gid = grib_new_from_index(iid)
1688
1689	f.close()
1690	g.close()
1691	h.close()
1692
1693
1694	grib_index_release(iid)

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

Download in other formats: