Context Navigation

source: flex_extract.git/source/python/classes/EcFlexpart.py @ e5884c9

ctbtodev

Last change on this file since e5884c9 was e5884c9, checked in by Anne Philipp <anne.philipp@…>, 5 years ago
minor fix for purefc check
Property mode set to `100644`
File size: 55.4 KB

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

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

Download in other formats: