Context Navigation

source: flex_extract.git/python/ECFlexpart.py @ 991df6a

ctbtodev

Last change on this file since 991df6a was 991df6a, checked in by Anne Philipp <anne.philipp@…>, 6 years ago
finished documentation (except plot_retrieved)
Property mode set to `100644`
File size: 54.1 KB

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

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

Download in other formats: