Context Navigation

source: flex_extract.git/python/EcFlexpart.py @ 5d42acd

ctbtodev

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

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

Download in other formats: