Context Navigation

source: flex_extract.git/python/EcFlexpart.py @ ff99eae

ctbtodev

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

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

Download in other formats: