Context Navigation

source: flex_extract.git/python/ECFlexpart.py @ efdb01a

ctbtodev

Last change on this file since efdb01a was efdb01a, checked in by Anne Philipp <anne.philipp@…>, 6 years ago
whole bunch of modifications due to new structure of ECMWFDATA, added basics of documentation, minor programming corrections
Property mode set to `100644`
File size: 53.0 KB

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

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

Download in other formats: