Context Navigation

source: flex_extract.git/python/EcFlexpart.py @ 2fb99de

ctbtodev

Last change on this file since 2fb99de was 2fb99de, checked in by Anne Philipp <anne.philipp@…>, 6 years ago
introduced config with path definitions and changed py files accordingly; Installation works; some tests were added for tarball making; Problems in submission to ecgate
Property mode set to `100644`
File size: 56.9 KB

Line
1	#!/usr/bin/env python
2	# -- coding: utf-8 --
3	#*******************************************************************************
4	# @Author: Anne Fouilloux (University of Oslo)
5	#
6	# @Date: October 2014
7	#
8	# @Change History:
9	#
10	# November 2015 - Leopold Haimberger (University of Vienna):
11	# - extended with class Control
12	# - removed functions mkdir_p, daterange, years_between, months_between
13	# - added functions darain, dapoly, to_param_id, init128, normal_exit,
14	# my_error, clean_up, install_args_and_control,
15	# interpret_args_and_control,
16	# - removed function __del__ in class EIFLexpart
17	# - added the following functions in EIFlexpart:
18	# - create_namelist
19	# - process_output
20	# - deacc_fluxes
21	# - modified existing EIFlexpart - functions for the use in
22	# flex_extract
23	# - retrieve also longer term forecasts, not only analyses and
24	# short term forecast data
25	# - added conversion into GRIB2
26	# - added conversion into .fp format for faster execution of FLEXPART
27	# (see https://www.flexpart.eu/wiki/FpCtbtoWo4FpFormat)
28	#
29	# February 2018 - Anne Philipp (University of Vienna):
30	# - applied PEP8 style guide
31	# - added documentation
32	# - removed function getFlexpartTime in class EcFlexpart
33	# - outsourced class ControlFile
34	# - outsourced class MarsRetrieval
35	# - changed class name from EIFlexpart to EcFlexpart
36	# - applied minor code changes (style)
37	# - removed "dead code" , e.g. retrieval of Q since it is not needed
38	# - removed "times" parameter from retrieve-method since it is not used
39	#
40	# @License:
41	# (C) Copyright 2014-2018.
42	#
43	# This software is licensed under the terms of the Apache Licence Version 2.0
44	# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
45	#
46	# @Class Description:
47	# FLEXPART needs grib files in a specifc format. All necessary data fields
48	# for one time step are stored in a single file. The class represents an
49	# instance with all the parameter and settings necessary for retrieving
50	# MARS data and modifing them so they are fitting FLEXPART need. The class
51	# is able to disaggregate the fluxes and convert grid types to the one needed
52	# by FLEXPART, therefore using the FORTRAN program.
53	#
54	# @Class Content:
55	# - __init__
56	# - write_namelist
57	# - retrieve
58	# - process_output
59	# - create
60	# - deacc_fluxes
61	#
62	# @Class Attributes:
63	# - dtime
64	# - basetime
65	# - server
66	# - marsclass
67	# - stream
68	# - resol
69	# - accuracy
70	# - number
71	# - expver
72	# - glevelist
73	# - area
74	# - grid
75	# - level
76	# - levelist
77	# - types
78	# - dates
79	# - area
80	# - gaussian
81	# - params
82	# - inputdir
83	# - outputfilelist
84	#
85	#*******************************************************************************
86	#pylint: disable=unsupported-assignment-operation
87	# this is disabled because its an error in pylint for this specific case
88	#pylint: disable=consider-using-enumerate
89	# this is not useful in this case
90	# ------------------------------------------------------------------------------
91	# MODULES
92	# ------------------------------------------------------------------------------
93	import subprocess
94	import shutil
95	import os
96	import glob
97	from datetime import datetime, timedelta
98	import numpy as np
99	from gribapi import grib_set, grib_index_select, grib_new_from_index, grib_get,\
100	grib_write, grib_get_values, grib_set_values, grib_release,\
101	grib_index_release, grib_index_get
102
103	# software specific classes and modules from flex_extract
104	import _config
105	from GribTools import GribTools
106	from tools import init128, to_param_id, silent_remove, product, my_error
107	from MarsRetrieval import MarsRetrieval
108	import disaggregation
109
110	# ------------------------------------------------------------------------------
111	# CLASS
112	# ------------------------------------------------------------------------------
113	class EcFlexpart(object):
114	'''
115	Class to retrieve FLEXPART specific ECMWF data.
116	'''
117	# --------------------------------------------------------------------------
118	# CLASS FUNCTIONS
119	# --------------------------------------------------------------------------
120	def __init__(self, c, fluxes=False):
121	'''
122	@Description:
123	Creates an object/instance of EcFlexpart with the
124	associated settings of its attributes for the retrieval.
125
126	@Input:
127	self: instance of EcFlexpart
128	The current object of the class.
129
130	c: instance of class ControlFile
131	Contains all the parameters of CONTROL file, which are e.g.:
132	DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
133	STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
134	LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
135	OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
136	ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
137	MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
138	DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
139
140	For more information about format and content of the parameter
141	see documentation.
142
143	fluxes: boolean, optional
144	Decides if the flux parameter settings are stored or
145	the rest of the parameter list.
146	Default value is False.
147
148	@Return:
149	<nothing>
150	'''
151
152	# different mars types for retrieving data for flexpart
153	self.types = dict()
154
155	if c.maxstep > len(c.type): # Pure forecast mode
156	c.type = [c.type[1]]
157	c.step = ['{:0>3}'.format(int(c.step[0]))]
158	c.time = [c.time[0]]
159	for i in range(1, c.maxstep + 1):
160	c.type.append(c.type[0])
161	c.step.append('{:0>3}'.format(i))
162	c.time.append(c.time[0])
163
164	self.inputdir = c.inputdir
165	self.basetime = c.basetime
166	self.dtime = c.dtime
167	i = 0
168	if fluxes and c.maxstep <= 24:
169	# no forecast beyond one day is needed!
170	# Thus, prepare flux data manually as usual
171	# with only forecast fields with start times at 00/12
172	# (but without 00/12 fields since these are
173	# the initialisation times of the flux fields
174	# and therefore are zero all the time)
175	self.types[c.type[1]] = {'times': '00/12', 'steps':
176	'{}/to/12/by/{}'.format(c.dtime, c.dtime)}
177	else:
178	for ty, st, ti in zip(c.type, c.step, c.time):
179	btlist = range(24)
180	if c.basetime == '12':
181	btlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
182	if c.basetime == '00':
183	btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0]
184
185	if i % int(c.dtime) == 0 and (i in btlist or c.maxstep > 24):
186
187	if ty not in self.types.keys():
188	self.types[ty] = {'times': '', 'steps': ''}
189
190	if ti not in self.types[ty]['times']:
191	if self.types[ty]['times']:
192	self.types[ty]['times'] += '/'
193	self.types[ty]['times'] += ti
194
195	if st not in self.types[ty]['steps']:
196	if self.types[ty]['steps']:
197	self.types[ty]['steps'] += '/'
198	self.types[ty]['steps'] += st
199	i += 1
200
201	self.marsclass = c.marsclass
202	self.stream = c.stream
203	self.number = c.number
204	self.resol = c.resol
205	self.accuracy = c.accuracy
206	self.level = c.level
207	self.expver = c.expver
208	# for gaussian grid retrieval
209	self.glevelist = '1/to/' + c.level
210
211	if hasattr(c, 'gaussian') and c.gaussian:
212	self.gaussian = c.gaussian
213	else:
214	self.gaussian = ''
215
216	if 'N' in c.grid: # Gaussian output grid
217	self.grid = c.grid
218	self.area = 'G'
219	else:
220	self.grid = '{}/{}'.format(int(c.grid) / 1000., int(c.grid) / 1000.)
221	self.area = '{}/{}/{}/{}'.format(int(c.upper) / 1000.,
222	int(c.left) / 1000.,
223	int(c.lower) / 1000.,
224	int(c.right) / 1000.)
225
226	self.outputfilelist = []
227
228
229	# Now comes the nasty part that deals with the different
230	# scenarios we have:
231	# 1) Calculation of etadot on
232	# a) Gaussian grid
233	# b) Output grid
234	# c) Output grid using parameter 77 retrieved from MARS
235	# 3) Calculation/Retrieval of omega
236	# 4) Download also data for WRF
237
238	# Different grids need different retrievals
239	# SH = Spherical Harmonics, GG = Gaussian Grid,
240	# OG = Output Grid, ML = MultiLevel, SL = SingleLevel
241	self.params = {'SH__ML': '', 'SH__SL': '',
242	'GG__ML': '', 'GG__SL': '',
243	'OG__ML': '', 'OG__SL': '',
244	'OG_OROLSM_SL': '', 'OG_acc_SL': ''}
245
246	if fluxes is False:
247	self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF']
248	# "SD/MSL/TCC/10U/10V/2T/2D/129/172"
249	self.params['OG__SL'] = ["141/151/164/165/166/167/168/129/172", \
250	'SFC', '1', self.grid]
251	if c.addpar:
252	if c.addpar[0] == '/':
253	c.addpar = c.addpar[1:]
254	self.params['OG__SL'][0] += '/' + '/'.join(c.addpar)
255
256	self.params['OG_OROLSM__SL'] = ["160/27/28/173", \
257	'SFC', '1', self.grid]
258
259	self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]
260
261	#if c.gauss == '0' and c.eta == '1':
262	if not c.gauss and c.eta:
263	# the simplest case
264	self.params['OG__ML'][0] += '/U/V/77'
265	#elif c.gauss == '0' and c.eta == '0':
266	elif not c.gauss and not c.eta:
267	# this is not recommended (inaccurate)
268	self.params['OG__ML'][0] += '/U/V'
269	#elif c.gauss == '1' and c.eta == '0':
270	elif c.gauss and not c.eta:
271	# this is needed for data before 2008, or for reanalysis data
272	self.params['GG__SL'] = ['Q', 'ML', '1', \
273	'{}'.format((int(self.resol) + 1) / 2)]
274	self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF']
275	else:
276	print('Warning: This is a very costly parameter combination, \
277	use only for debugging!')
278	self.params['GG__SL'] = ['Q', 'ML', '1', \
279	'{}'.format((int(self.resol) + 1) / 2)]
280	self.params['GG__ML'] = ['U/V/D/77', 'ML', self.glevelist, \
281	'{}'.format((int(self.resol) + 1) / 2)]
282
283	if c.omega:
284	self.params['OG__ML'][0] += '/W'
285
286	# add cloud water content if necessary
287	if c.cwc:
288	self.params['OG__ML'][0] += '/CLWC/CIWC'
289
290	# add vorticity and geopotential height for WRF if necessary
291	if c.wrf:
292	self.params['OG__ML'][0] += '/Z/VO'
293	if '/D' not in self.params['OG__ML'][0]:
294	self.params['OG__ML'][0] += '/D'
295	#wrf_sfc = 'sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/ /
296	# stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper()
297	wrf_sfc = '134/235/167/165/166/168/129/172/34/31/141/ \
298	139/170/183/236/39/40/41/42'.upper()
299	lwrt_sfc = wrf_sfc.split('/')
300	for par in lwrt_sfc:
301	if par not in self.params['OG__SL'][0]:
302	self.params['OG__SL'][0] += '/' + par
303
304	else:
305	self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR", \
306	'SFC', '1', self.grid]
307
308	# if needed, add additional WRF specific parameters here
309
310	return
311
312
313	def write_namelist(self, c, filename):
314	'''
315	@Description:
316	Creates a namelist file in the temporary directory and writes
317	the following values to it: maxl, maxb, mlevel,
318	mlevelist, mnauf, metapar, rlo0, rlo1, rla0, rla1,
319	momega, momegadiff, mgauss, msmooth, meta, metadiff, mdpdeta
320
321	@Input:
322	self: instance of EcFlexpart
323	The current object of the class.
324
325	c: instance of class ControlFile
326	Contains all the parameters of CONTROL files, which are e.g.:
327	DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
328	STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
329	LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
330	OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
331	ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
332	MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
333	DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
334
335	For more information about format and content of the parameter
336	see documentation.
337
338	filename: string
339	Name of the namelist file.
340
341	@Return:
342	<nothing>
343	'''
344
345	self.inputdir = c.inputdir
346	area = np.asarray(self.area.split('/')).astype(float)
347	grid = np.asarray(self.grid.split('/')).astype(float)
348
349	if area[1] > area[3]:
350	area[1] -= 360
351	maxl = int((area[3] - area[1]) / grid[1]) + 1
352	maxb = int((area[0] - area[2]) / grid[0]) + 1
353
354	with open(self.inputdir + '/' + filename, 'w') as f:
355	f.write('&NAMGEN\n')
356	f.write(',\n '.join(['maxl = ' + str(maxl), 'maxb = ' + str(maxb),
357	'mlevel = ' + str(self.level),
358	'mlevelist = ' + '"' + str(self.levelist)
359	+ '"',
360	'mnauf = ' + str(self.resol),
361	'metapar = ' + '77',
362	'rlo0 = ' + str(area[1]),
363	'rlo1 = ' + str(area[3]),
364	'rla0 = ' + str(area[2]),
365	'rla1 = ' + str(area[0]),
366	'momega = ' + str(c.omega),
367	'momegadiff = ' + str(c.omegadiff),
368	'mgauss = ' + str(c.gauss),
369	'msmooth = ' + str(c.smooth),
370	'meta = ' + str(c.eta),
371	'metadiff = ' + str(c.etadiff),
372	'mdpdeta = ' + str(c.dpdeta)]))
373
374	f.write('\n/\n')
375
376	return
377
378	def retrieve(self, server, dates, request, inputdir='.'):
379	'''
380	@Description:
381	Finalizing the retrieval information by setting final details
382	depending on grid type.
383	Prepares MARS retrievals per grid type and submits them.
384
385	@Input:
386	self: instance of EcFlexpart
387	The current object of the class.
388
389	server: instance of ECMWFService or ECMWFDataServer
390	The connection to the ECMWF server. This is different
391	for member state users which have full access and non
392	member state users which have only access to the public
393	data sets. The decision is made from command line argument
394	"public"; for public access its True (ECMWFDataServer)
395	for member state users its False (ECMWFService)
396
397	dates: string
398	Contains start and end date of the retrieval in the format
399	"YYYYMMDD/to/YYYYMMDD"
400
401	inputdir: string, optional
402	Path to the directory where the retrieved data is about
403	to be stored. The default is the current directory ('.').
404
405	@Return:
406	<nothing>
407	'''
408	self.dates = dates
409	self.server = server
410	self.inputdir = inputdir
411	oro = False
412
413	for ftype in self.types:
414	for pk, pv in self.params.iteritems():
415	if isinstance(pv, str):
416	continue
417	mftype = '' + ftype
418	mftime = self.types[ftype]['times']
419	mfstep = self.types[ftype]['steps']
420	mfdate = self.dates
421	mfstream = self.stream
422	mftarget = self.inputdir + "/" + ftype + pk + '.' + \
423	self.dates.split('/')[0] + '.' + str(os.getppid()) +\
424	'.' + str(os.getpid()) + ".grb"
425	if pk == 'OG__SL':
426	pass
427	if pk == 'OG_OROLSM__SL':
428	if not oro:
429	mfstream = 'OPER'
430	mftype = 'AN'
431	mftime = '00'
432	mfstep = '000'
433	mfdate = self.dates.split('/')[0]
434	mftarget = self.inputdir + "/" + pk + '.' + mfdate + \
435	'.' + str(os.getppid()) + '.' + \
436	str(os.getpid()) + ".grb"
437	oro = True
438	else:
439	continue
440	if pk == 'GG__SL' and pv[0] == 'Q':
441	area = ""
442	gaussian = 'reduced'
443	else:
444	area = self.area
445	gaussian = self.gaussian
446
447	# ------ on demand path --------------------------------------------------
448	if not self.basetime:
449	MR = MarsRetrieval(self.server,
450	marsclass=self.marsclass,
451	stream=mfstream,
452	type=mftype,
453	levtype=pv[1],
454	levelist=pv[2],
455	resol=self.resol,
456	gaussian=gaussian,
457	accuracy=self.accuracy,
458	grid=pv[3],
459	target=mftarget,
460	area=area,
461	date=mfdate,
462	time=mftime,
463	number=self.number,
464	step=mfstep,
465	expver=self.expver,
466	param=pv[0])
467
468	if request == 0:
469	MR.display_info()
470	MR.data_retrieve()
471	elif request == 1:
472	MR.print_info()
473	elif request == 2:
474	MR.print_info()
475	MR.display_info()
476	MR.data_retrieve()
477	else:
478	print('Failure')
479	# ------ operational path ------------------------------------------------
480	else:
481	# check if mars job requests fields beyond basetime.
482	# If yes eliminate those fields since they may not
483	# be accessible with user's credentials
484	if 'by' in mfstep:
485	sm1 = 2
486	else:
487	sm1 = -1
488
489	if 'by' in mftime:
490	tm1 = 2
491	else:
492	tm1 = -1
493
494	maxdate = datetime.strptime(mfdate.split('/')[-1] +
495	mftime.split('/')[tm1],
496	'%Y%m%d%H')
497	istep = int(mfstep.split('/')[sm1])
498	maxtime = maxdate + timedelta(hours=istep)
499
500	elimit = datetime.strptime(mfdate.split('/')[-1] +
501	self.basetime, '%Y%m%d%H')
502
503	if self.basetime == '12':
504	# -------------- flux data ----------------------------
505	if 'acc' in pk:
506
507	# Strategy:
508	# if maxtime-elimit >= 24h reduce date by 1,
509	# if 12h <= maxtime-elimit<12h reduce time for last date
510	# if maxtime-elimit<12h reduce step for last time
511	# A split of the MARS job into 2 is likely necessary.
512	maxtime = elimit - timedelta(hours=24)
513	mfdate = '/'.join(['/'.join(mfdate.split('/')[:-1]),
514	datetime.strftime(maxtime,
515	'%Y%m%d')])
516
517	MR = MarsRetrieval(self.server,
518	marsclass=self.marsclass,
519	stream=self.stream,
520	type=mftype,
521	levtype=pv[1],
522	levelist=pv[2],
523	resol=self.resol,
524	gaussian=gaussian,
525	accuracy=self.accuracy,
526	grid=pv[3],
527	target=mftarget,
528	area=area,
529	date=mfdate,
530	time=mftime,
531	number=self.number,
532	step=mfstep,
533	expver=self.expver,
534	param=pv[0])
535
536	MR.display_info()
537	MR.data_retrieve()
538
539	maxtime = elimit - timedelta(hours=12)
540	mfdate = datetime.strftime(maxtime, '%Y%m%d')
541	mftime = '00'
542	mftarget = self.inputdir + "/" + ftype + pk + \
543	'.' + mfdate + '.' + str(os.getppid()) +\
544	'.' + str(os.getpid()) + ".grb"
545
546	MR = MarsRetrieval(self.server,
547	marsclass=self.marsclass,
548	stream=self.stream,
549	type=mftype,
550	levtype=pv[1],
551	levelist=pv[2],
552	resol=self.resol,
553	gaussian=gaussian,
554	accuracy=self.accuracy,
555	grid=pv[3],
556	target=mftarget,
557	area=area,
558	date=mfdate,
559	time=mftime,
560	number=self.number,
561	step=mfstep,
562	expver=self.expver,
563	param=pv[0])
564
565	MR.display_info()
566
567	MR.data_retrieve()
568	# -------------- non flux data ------------------------
569	else:
570	MR = MarsRetrieval(self.server,
571	marsclass=self.marsclass,
572	stream=self.stream,
573	type=mftype,
574	levtype=pv[1],
575	levelist=pv[2],
576	resol=self.resol,
577	gaussian=gaussian,
578	accuracy=self.accuracy,
579	grid=pv[3],
580	target=mftarget,
581	area=area,
582	date=mfdate,
583	time=mftime,
584	number=self.number,
585	step=mfstep,
586	expver=self.expver,
587	param=pv[0])
588
589	MR.display_info()
590	MR.data_retrieve()
591	else: # basetime == 0 ??? #AP
592
593	maxtime = elimit - timedelta(hours=24)
594	mfdate = datetime.strftime(maxtime, '%Y%m%d')
595	mftimesave = ''.join(mftime)
596
597	if '/' in mftime:
598	times = mftime.split('/')
599	while ((int(times[0]) +
600	int(mfstep.split('/')[0]) <= 12) and
601	(pk != 'OG_OROLSM__SL') and 'acc' not in pk):
602	times = times[1:]
603	if len(times) > 1:
604	mftime = '/'.join(times)
605	else:
606	mftime = times[0]
607
608	MR = MarsRetrieval(self.server,
609	marsclass=self.marsclass,
610	stream=self.stream,
611	type=mftype,
612	levtype=pv[1],
613	levelist=pv[2],
614	resol=self.resol,
615	gaussian=gaussian,
616	accuracy=self.accuracy,
617	grid=pv[3],
618	target=mftarget,
619	area=area,
620	date=mfdate,
621	time=mftime,
622	number=self.number,
623	step=mfstep,
624	expver=self.expver,
625	param=pv[0])
626
627	MR.display_info()
628	MR.data_retrieve()
629
630	if (int(mftimesave.split('/')[0]) == 0 and
631	int(mfstep.split('/')[0]) == 0 and
632	pk != 'OG_OROLSM__SL'):
633
634	mfdate = datetime.strftime(elimit, '%Y%m%d')
635	mftime = '00'
636	mfstep = '000'
637	mftarget = self.inputdir + "/" + ftype + pk + \
638	'.' + mfdate + '.' + str(os.getppid()) +\
639	'.' + str(os.getpid()) + ".grb"
640
641	MR = MarsRetrieval(self.server,
642	marsclass=self.marsclass,
643	stream=self.stream,
644	type=mftype,
645	levtype=pv[1],
646	levelist=pv[2],
647	resol=self.resol,
648	gaussian=gaussian,
649	accuracy=self.accuracy,
650	grid=pv[3],
651	target=mftarget,
652	area=area,
653	date=mfdate,
654	time=mftime,
655	number=self.number,
656	step=mfstep,
657	expver=self.expver,
658	param=pv[0])
659
660	MR.display_info()
661	MR.data_retrieve()
662
663	if request == 0 or request == 2:
664	print('MARS retrieve done ... ')
665	elif request == 1:
666	print('MARS request printed ...')
667
668	return
669
670
671	def process_output(self, c):
672	'''
673	@Description:
674	The grib files are postprocessed depending on the selection in
675	CONTROL file. The resulting files are moved to the output
676	directory if its not equla to the input directory.
677	The following modifications might be done if
678	properly switched in CONTROL file:
679	GRIB2 - Conversion to GRIB2
680	ECTRANS - Transfer of files to gateway server
681	ECSTORAGE - Storage at ECMWF server
682	GRIB2FLEXPART - Conversion of GRIB files to FLEXPART binary format
683
684	@Input:
685	self: instance of EcFlexpart
686	The current object of the class.
687
688	c: instance of class ControlFile
689	Contains all the parameters of CONTROL file, which are e.g.:
690	DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
691	STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
692	LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
693	OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
694	ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
695	MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
696	DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
697
698	For more information about format and content of the parameter
699	see documentation.
700
701	@Return:
702	<nothing>
703
704	'''
705
706	print('\n\nPostprocessing:\n Format: {}\n'.format(c.format))
707
708	if not c.ecapi:
709	print('ecstorage: {}\n ecfsdir: {}\n'.
710	format(c.ecstorage, c.ecfsdir))
711	#if not hasattr(c, 'gateway'):
712	# c.gateway = os.getenv('GATEWAY')
713	#if not hasattr(c, 'destination'):
714	# c.destination = os.getenv('DESTINATION')
715	print('ectrans: {}\n gateway: {}\n destination: {}\n '
716	.format(c.ectrans, c.gateway, c.destination))
717
718	print('Output filelist: \n')
719	print(self.outputfilelist)
720
721	if c.format.lower() == 'grib2':
722	for ofile in self.outputfilelist:
723	p = subprocess.check_call(['grib_set', '-s', 'edition=2, \
724	productDefinitionTemplateNumber=8',
725	ofile, ofile + '_2'])
726	p = subprocess.check_call(['mv', ofile + '_2', ofile])
727
728	if c.ectrans and not c.ecapi:
729	for ofile in self.outputfilelist:
730	p = subprocess.check_call(['ectrans', '-overwrite', '-gateway',
731	c.gateway, '-remote', c.destination,
732	'-source', ofile])
733	#print('ectrans:', p)
734
735	if c.ecstorage and not c.ecapi:
736	for ofile in self.outputfilelist:
737	p = subprocess.check_call(['ecp', '-o', ofile,
738	os.path.expandvars(c.ecfsdir)])
739
740	if c.outputdir != c.inputdir:
741	for ofile in self.outputfilelist:
742	p = subprocess.check_call(['mv', ofile, c.outputdir])
743
744	# prepare environment for the grib2flexpart run
745	# to convert grib to flexpart binary
746	if c.grib2flexpart:
747
748	# generate AVAILABLE file
749	# Example of AVAILABLE file data:
750	# 20131107 000000 EN13110700 ON DISC
751	clist = []
752	for ofile in self.outputfilelist:
753	fname = ofile.split('/')
754	if '.' in fname[-1]:
755	l = fname[-1].split('.')
756	timestamp = datetime.strptime(l[0][-6:] + l[1],
757	'%y%m%d%H')
758	timestamp += timedelta(hours=int(l[2]))
759	cdate = datetime.strftime(timestamp, '%Y%m%d')
760	chms = datetime.strftime(timestamp, '%H%M%S')
761	else:
762	cdate = '20' + fname[-1][-8:-2]
763	chms = fname[-1][-2:] + '0000'
764	clist.append(cdate + ' ' + chms + ' '*6 +
765	fname[-1] + ' '*14 + 'ON DISC')
766	clist.sort()
767	with open(c.outputdir + '/' + 'AVAILABLE', 'w') as f:
768	f.write('\n'.join(clist) + '\n')
769
770	# generate pathnames file
771	pwd = os.path.abspath(c.outputdir)
772	with open(pwd + '/pathnames', 'w') as f:
773	f.write(pwd + '/Options/\n')
774	f.write(pwd + '/\n')
775	f.write(pwd + '/\n')
776	f.write(pwd + '/AVAILABLE\n')
777	f.write(' = == = == = == = == = == == = \n')
778
779	# create Options dir if necessary
780	if not os.path.exists(pwd + '/Options'):
781	os.makedirs(pwd+'/Options')
782
783	# read template COMMAND file
784	with open(os.path.expandvars(os.path.expanduser(
785	c.flexpart_root_scripts)) + '/../Options/COMMAND', 'r') as f:
786	lflist = f.read().split('\n')
787
788	# find index of list where to put in the
789	# date and time information
790	# usually after the LDIRECT parameter
791	i = 0
792	for l in lflist:
793	if 'LDIRECT' in l.upper():
794	break
795	i += 1
796
797	# insert the date and time information of run start and end
798	# into the list of lines of COMMAND file
799	lflist = lflist[:i+1] + \
800	[clist[0][:16], clist[-1][:16]] + \
801	lflist[i+3:]
802
803	# write the new COMMAND file
804	with open(pwd + '/Options/COMMAND', 'w') as g:
805	g.write('\n'.join(lflist) + '\n')
806
807	# change to outputdir and start the grib2flexpart run
808	# afterwards switch back to the working dir
809	os.chdir(c.outputdir)
810	p = subprocess.check_call([
811	os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))
812	+ '/../FLEXPART_PROGRAM/grib2flexpart', 'useAvailable', '.'])
813	os.chdir(pwd)
814
815	return
816
817	def create(self, inputfiles, c):
818	'''
819	@Description:
820	This method is based on the ECMWF example index.py
821	https://software.ecmwf.int/wiki/display/GRIB/index.py
822
823	An index file will be created which depends on the combination
824	of "date", "time" and "stepRange" values. This is used to iterate
825	over all messages in each grib file which were passed through the
826	parameter "inputfiles" to seperate specific parameters into fort.*
827	files. Afterwards the FORTRAN program Convert2 is called to convert
828	the data fields all to the same grid and put them in one file
829	per unique time step (combination of "date", "time" and
830	"stepRange").
831
832	@Input:
833	self: instance of EcFlexpart
834	The current object of the class.
835
836	inputfiles: instance of UioFiles
837	Contains a list of files.
838
839	c: instance of class ControlFile
840	Contains all the parameters of CONTROL files, which are e.g.:
841	DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
842	STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
843	LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
844	OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
845	ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
846	MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
847	DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
848
849	For more information about format and content of the parameter
850	see documentation.
851
852	@Return:
853	<nothing>
854	'''
855
856	table128 = init128(_config.PATH_GRIBTABLE)
857	wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
858	stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
859	table128)
860
861	index_keys = ["date", "time", "step"]
862	indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX)
863	silent_remove(indexfile)
864	grib = GribTools(inputfiles.files)
865	# creates new index file
866	iid = grib.index(index_keys=index_keys, index_file=indexfile)
867
868	# read values of index keys
869	index_vals = []
870	for key in index_keys:
871	index_vals.append(grib_index_get(iid, key))
872	print(index_vals[-1])
873	# index_vals looks for example like:
874	# index_vals[0]: ('20171106', '20171107', '20171108') ; date
875	# index_vals[1]: ('0', '1200', '1800', '600') ; time
876	# index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
877
878	fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
879	'17':None, '19':None, '21':None, '22':None, '20':None}
880
881	for prod in product(*index_vals):
882	# e.g. prod = ('20170505', '0', '12')
883	# ( date ,time, step)
884	# per date e.g. time = 0, 1200
885	# per time e.g. step = 3, 6, 9, 12
886	# flag for Fortran program CONVERT2 and file merging
887	convertFlag = False
888	print('current prod: ', prod)
889	# e.g. prod = ('20170505', '0', '12')
890	# ( date ,time, step)
891	# per date e.g. time = 0, 600, 1200, 1800
892	# per time e.g. step = 0, 3, 6, 9, 12
893	for i in range(len(index_keys)):
894	grib_index_select(iid, index_keys[i], prod[i])
895
896	# get first id from current product
897	gid = grib_new_from_index(iid)
898
899	# if there is data for this product combination
900	# prepare some date and time parameter before reading the data
901	if gid is not None:
902	# Combine all temporary data files into final grib file if
903	# gid is at least one time not None. Therefore set convertFlag
904	# to save information. The fortran program CONVERT2 is also
905	# only done if convertFlag is True
906	convertFlag = True
907	# remove old fort.* files and open new ones
908	# they are just valid for a single product
909	for k, f in fdict.iteritems():
910	fortfile = os.path.join(c.inputdir, 'fort.' + k)
911	silent_remove(fortfile)
912	fdict[k] = open(fortfile, 'w')
913
914	cdate = str(grib_get(gid, 'date'))
915	time = grib_get(gid, 'time')
916	step = grib_get(gid, 'step')
917	# create correct timestamp from the three time informations
918	# date, time, step
919	timestamp = datetime.strptime(cdate + '{:0>2}'.format(time/100),
920	'%Y%m%d%H')
921	timestamp += timedelta(hours=int(step))
922	cdateH = datetime.strftime(timestamp, '%Y%m%d%H')
923
924	if c.basetime is not None:
925	slimit = datetime.strptime(c.start_date + '00', '%Y%m%d%H')
926	bt = '23'
927	if c.basetime == '00':
928	bt = '00'
929	slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
930	- timedelta(hours=12-int(c.dtime))
931	if c.basetime == '12':
932	bt = '12'
933	slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
934	- timedelta(hours=12-int(c.dtime))
935
936	elimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')
937
938	if timestamp < slimit or timestamp > elimit:
939	continue
940	else:
941	pass
942
943	try:
944	if c.wrf:
945	if 'olddate' not in locals() or cdate != olddate:
946	fwrf = open(os.path.join(c.outputdir,
947	'WRF' + cdate + '.{:0>2}'.format(time) +
948	'.000.grb2'), 'w')
949	olddate = cdate[:]
950	except AttributeError:
951	pass
952
953	# helper variable to remember which fields were already used.
954	savedfields = []
955	while 1:
956	if gid is None:
957	break
958	paramId = grib_get(gid, 'paramId')
959	gridtype = grib_get(gid, 'gridType')
960	levtype = grib_get(gid, 'typeOfLevel')
961	if paramId == 133 and gridtype == 'reduced_gg':
962	# Specific humidity (Q.grb) is used as a template only
963	# so we need the first we "meet"
964	with open(os.path.join(c.inputdir, 'fort.18'), 'w') as fout:
965	grib_write(gid, fout)
966	elif paramId == 131 or paramId == 132:
967	grib_write(gid, fdict['10'])
968	elif paramId == 130:
969	grib_write(gid, fdict['11'])
970	elif paramId == 133 and gridtype != 'reduced_gg':
971	grib_write(gid, fdict['17'])
972	elif paramId == 152:
973	grib_write(gid, fdict['12'])
974	elif paramId == 155 and gridtype == 'sh':
975	grib_write(gid, fdict['13'])
976	elif paramId in [129, 138, 155] and levtype == 'hybrid' \
977	and c.wrf:
978	pass
979	elif paramId == 246 or paramId == 247:
980	# cloud liquid water and ice
981	if paramId == 246:
982	clwc = grib_get_values(gid)
983	else:
984	clwc += grib_get_values(gid)
985	grib_set_values(gid, clwc)
986	grib_set(gid, 'paramId', 201031)
987	grib_write(gid, fdict['22'])
988	elif paramId == 135:
989	grib_write(gid, fdict['19'])
990	elif paramId == 77:
991	grib_write(gid, fdict['21'])
992	else:
993	if paramId not in savedfields:
994	grib_write(gid, fdict['16'])
995	savedfields.append(paramId)
996	else:
997	print('duplicate ' + str(paramId) + ' not written')
998
999	try:
1000	if c.wrf:
1001	# model layer
1002	if levtype == 'hybrid' and \
1003	paramId in [129, 130, 131, 132, 133, 138, 155]:
1004	grib_write(gid, fwrf)
1005	# sfc layer
1006	elif paramId in wrfpars:
1007	grib_write(gid, fwrf)
1008	except AttributeError:
1009	pass
1010
1011	grib_release(gid)
1012	gid = grib_new_from_index(iid)
1013
1014	for f in fdict.values():
1015	f.close()
1016
1017	# call for CONVERT2 if flag is True
1018	if convertFlag:
1019	pwd = os.getcwd()
1020	os.chdir(c.inputdir)
1021	if os.stat('fort.21').st_size == 0 and c.eta:
1022	print('Parameter 77 (etadot) is missing, most likely it is \
1023	not available for this type or date/time\n')
1024	print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
1025	my_error(c.mailfail, 'fort.21 is empty while parameter eta \
1026	is set to 1 in CONTROL file')
1027
1028	# create the corresponding output file fort.15
1029	# (generated by CONVERT2) + fort.16 (paramId 167 and 168)
1030	p = subprocess.check_call([os.path.join(c.exedir, 'CONVERT2')],
1031	shell=True)
1032	os.chdir(pwd)
1033
1034	# create final output filename, e.g. EN13040500 (ENYYMMDDHH)
1035	fnout = os.path.join(c.inputdir, c.prefix)
1036	if c.maxstep > 12:
1037	suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \
1038	'.{:0>3}'.format(step)
1039	else:
1040	suffix = cdateH[2:10]
1041	fnout += suffix
1042	print("outputfile = " + fnout)
1043	self.outputfilelist.append(fnout) # needed for final processing
1044
1045	# create outputfile and copy all data from intermediate files
1046	# to the outputfile (final GRIB files)
1047	orolsm = os.path.basename(glob.glob(
1048	c.inputdir + '/OG_OROLSM__SL..' + c.ppid + '')[0])
1049	fluxfile = 'flux' + cdate[0:2] + suffix
1050	if not c.cwc:
1051	flist = ['fort.15', fluxfile, 'fort.16', orolsm]
1052	else:
1053	flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]
1054
1055	with open(fnout, 'wb') as fout:
1056	for f in flist:
1057	shutil.copyfileobj(open(os.path.join(c.inputdir, f),
1058	'rb'), fout)
1059
1060	if c.omega:
1061	with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout:
1062	shutil.copyfileobj(
1063	open(os.path.join(c.inputdir, 'fort.25'),
1064	'rb'), fout)
1065	else:
1066	pass
1067
1068	if c.wrf:
1069	fwrf.close()
1070
1071	grib_index_release(iid)
1072
1073	return
1074
1075	def deacc_fluxes(self, inputfiles, c):
1076	'''
1077	@Description:
1078	Goes through all flux fields in ordered time and de-accumulate
1079	the fields. Afterwards the fields are disaggregated in time.
1080	Different versions of disaggregation is provided for rainfall
1081	data (darain, modified linear) and the surface fluxes and
1082	stress data (dapoly, cubic polynomial).
1083
1084	@Input:
1085	self: instance of EcFlexpart
1086	The current object of the class.
1087
1088	inputfiles: instance of UioFiles
1089	Contains a list of files.
1090
1091	c: instance of class ControlFile
1092	Contains all the parameters of CONTROL file, which are e.g.:
1093	DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
1094	STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
1095	LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
1096	OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
1097	ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
1098	MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
1099	DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
1100
1101	For more information about format and content of the parameter
1102	see documentation.
1103
1104	@Return:
1105	<nothing>
1106	'''
1107
1108	table128 = init128(_config.PATH_GRIBTABLE)
1109	pars = to_param_id(self.params['OG_acc_SL'][0], table128)
1110	index_keys = ["date", "time", "step"]
1111	indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX)
1112	silent_remove(indexfile)
1113	grib = GribTools(inputfiles.files)
1114	# creates new index file
1115	iid = grib.index(index_keys=index_keys, index_file=indexfile)
1116
1117	# read values of index keys
1118	index_vals = []
1119	for key in index_keys:
1120	key_vals = grib_index_get(iid, key)
1121	print(key_vals)
1122	# have to sort the steps for disaggregation,
1123	# therefore convert to int first
1124	if key == 'step':
1125	key_vals = [int(k) for k in key_vals]
1126	key_vals.sort()
1127	key_vals = [str(k) for k in key_vals]
1128	index_vals.append(key_vals)
1129	# index_vals looks for example like:
1130	# index_vals[0]: ('20171106', '20171107', '20171108') ; date
1131	# index_vals[1]: ('0', '1200') ; time
1132	# index_vals[2]: (3', '6', '9', '12') ; stepRange
1133
1134	valsdict = {}
1135	svalsdict = {}
1136	stepsdict = {}
1137	for p in pars:
1138	valsdict[str(p)] = []
1139	svalsdict[str(p)] = []
1140	stepsdict[str(p)] = []
1141
1142	print('maxstep: ', c.maxstep)
1143
1144	for prod in product(*index_vals):
1145	# e.g. prod = ('20170505', '0', '12')
1146	# ( date ,time, step)
1147	# per date e.g. time = 0, 1200
1148	# per time e.g. step = 3, 6, 9, 12
1149	print('current prod: ', prod)
1150	for i in range(len(index_keys)):
1151	grib_index_select(iid, index_keys[i], prod[i])
1152
1153	# get first id from current product
1154	gid = grib_new_from_index(iid)
1155
1156	# if there is data for this product combination
1157	# prepare some date and time parameter before reading the data
1158	if gid is not None:
1159	cdate = grib_get(gid, 'date')
1160	time = grib_get(gid, 'time')
1161	step = grib_get(gid, 'step')
1162	# date+time+step-2*dtime
1163	# (since interpolated value valid for step-2*dtime)
1164	sdate = datetime(year=cdate/10000,
1165	month=(cdate % 10000)/100,
1166	day=(cdate % 100),
1167	hour=time/100)
1168	fdate = sdate + timedelta(hours=step-2*int(c.dtime))
1169	sdates = sdate + timedelta(hours=step)
1170	elimit = None
1171	else:
1172	break
1173
1174	if c.maxstep > 12:
1175	fnout = os.path.join(c.inputdir, 'flux' +
1176	sdate.strftime('%Y%m%d') +
1177	'.{:0>2}'.format(time/100) +
1178	'.{:0>3}'.format(step-2*int(c.dtime)))
1179	gnout = os.path.join(c.inputdir, 'flux' +
1180	sdate.strftime('%Y%m%d') +
1181	'.{:0>2}'.format(time/100) +
1182	'.{:0>3}'.format(step-int(c.dtime)))
1183	hnout = os.path.join(c.inputdir, 'flux' +
1184	sdate.strftime('%Y%m%d') +
1185	'.{:0>2}'.format(time/100) +
1186	'.{:0>3}'.format(step))
1187	else:
1188	fnout = os.path.join(c.inputdir, 'flux' +
1189	fdate.strftime('%Y%m%d%H'))
1190	gnout = os.path.join(c.inputdir, 'flux' +
1191	(fdate + timedelta(hours=int(c.dtime))).
1192	strftime('%Y%m%d%H'))
1193	hnout = os.path.join(c.inputdir, 'flux' +
1194	sdates.strftime('%Y%m%d%H'))
1195
1196	print("outputfile = " + fnout)
1197	f_handle = open(fnout, 'w')
1198	g_handle = open(gnout, 'w')
1199	h_handle = open(hnout, 'w')
1200
1201	# read message for message and store relevant data fields
1202	# data keywords are stored in pars
1203	while 1:
1204	if gid is None:
1205	break
1206	cparamId = str(grib_get(gid, 'paramId'))
1207	step = grib_get(gid, 'step')
1208	atime = grib_get(gid, 'time')
1209	ni = grib_get(gid, 'Ni')
1210	nj = grib_get(gid, 'Nj')
1211	if cparamId in valsdict.keys():
1212	values = grib_get_values(gid)
1213	vdp = valsdict[cparamId]
1214	svdp = svalsdict[cparamId]
1215	sd = stepsdict[cparamId]
1216
1217	if cparamId == '142' or cparamId == '143':
1218	fak = 1. / 1000.
1219	else:
1220	fak = 3600.
1221
1222	values = (np.reshape(values, (nj, ni))).flatten() / fak
1223	vdp.append(values[:]) # save the accumulated values
1224	if step <= int(c.dtime):
1225	svdp.append(values[:] / int(c.dtime))
1226	else: # deaccumulate values
1227	svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
1228
1229	print(cparamId, atime, step, len(values),
1230	values[0], np.std(values))
1231	# save the 1/3-hourly or specific values
1232	# svdp.append(values[:])
1233	sd.append(step)
1234	# len(svdp) correspond to the time
1235	if len(svdp) >= 3:
1236	if len(svdp) > 3:
1237	if cparamId == '142' or cparamId == '143':
1238	values = disaggregation.darain(svdp)
1239	else:
1240	values = disaggregation.dapoly(svdp)
1241
1242	if not (step == c.maxstep and c.maxstep > 12 \
1243	or sdates == elimit):
1244	vdp.pop(0)
1245	svdp.pop(0)
1246	else:
1247	if c.maxstep > 12:
1248	values = svdp[1]
1249	else:
1250	values = svdp[0]
1251
1252	grib_set_values(gid, values)
1253	if c.maxstep > 12:
1254	grib_set(gid, 'step', max(0, step-2*int(c.dtime)))
1255	else:
1256	grib_set(gid, 'step', 0)
1257	grib_set(gid, 'time', fdate.hour*100)
1258	grib_set(gid, 'date', fdate.year*10000 +
1259	fdate.month*100+fdate.day)
1260	grib_write(gid, f_handle)
1261
1262	if c.basetime is not None:
1263	elimit = datetime.strptime(c.end_date +
1264	c.basetime, '%Y%m%d%H')
1265	else:
1266	elimit = sdate + timedelta(2*int(c.dtime))
1267
1268	# squeeze out information of last two steps contained
1269	# in svdp
1270	# if step+int(c.dtime) == c.maxstep and c.maxstep>12
1271	# or sdates+timedelta(hours = int(c.dtime))
1272	# >= elimit:
1273	# Note that svdp[0] has not been popped in this case
1274
1275	if step == c.maxstep and c.maxstep > 12 or \
1276	sdates == elimit:
1277
1278	values = svdp[3]
1279	grib_set_values(gid, values)
1280	grib_set(gid, 'step', 0)
1281	truedatetime = fdate + timedelta(hours=
1282	2*int(c.dtime))
1283	grib_set(gid, 'time', truedatetime.hour * 100)
1284	grib_set(gid, 'date', truedatetime.year * 10000 +
1285	truedatetime.month * 100 +
1286	truedatetime.day)
1287	grib_write(gid, h_handle)
1288
1289	#values = (svdp[1]+svdp[2])/2.
1290	if cparamId == '142' or cparamId == '143':
1291	values = disaggregation.darain(list(reversed(svdp)))
1292	else:
1293	values = disaggregation.dapoly(list(reversed(svdp)))
1294
1295	grib_set(gid, 'step', 0)
1296	truedatetime = fdate + timedelta(hours=int(c.dtime))
1297	grib_set(gid, 'time', truedatetime.hour * 100)
1298	grib_set(gid, 'date', truedatetime.year * 10000 +
1299	truedatetime.month * 100 +
1300	truedatetime.day)
1301	grib_set_values(gid, values)
1302	grib_write(gid, g_handle)
1303
1304	grib_release(gid)
1305
1306	gid = grib_new_from_index(iid)
1307
1308	f_handle.close()
1309	g_handle.close()
1310	h_handle.close()
1311
1312	grib_index_release(iid)
1313
1314	return

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

Download in other formats: