Context Navigation

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

ctbtodev

Last change on this file since c274d9a was c274d9a, checked in by Anne Philipp <anne.philipp@…>, 5 years ago
added rrint param and updated deacc_fluxes function to store rain fields
Property mode set to `100644`
File size: 55.7 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.maxstep > len(c.type) and 'AN' not in c.type:
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 c.maxstep <= 24:
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.maxstep > 24):
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.maxstep > 12:
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.maxstep > 12 \
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.maxstep > 12:
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.maxstep > 12:
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 contained
920	# in deac_vals[parId]
921	# if step+int(c.dtime) == c.maxstep and c.maxstep>12
922	# or t_dt+timedelta(hours = int(c.dtime))
923	# >= t_enddate:
924	# Note that deac_vals[parId][0] has not been popped in this case
925
926	if step == c.maxstep and c.maxstep > 12 or \
927	t_dt == t_enddate:
928
929	values = deac_vals[parId][3]
930	codes_set_values(gid, values)
931	codes_set(gid, 'stepRange', 0)
932	truedatetime = t_m2dt + timedelta(hours=
933	2*int(c.dtime))
934	codes_set(gid, 'time', truedatetime.hour * 100)
935	codes_set(gid, 'date', truedatetime.year * 10000 +
936	truedatetime.month * 100 +
937	truedatetime.day)
938	codes_write(gid, h_handle)
939
940	#values = (deac_vals[parId][1]+deac_vals[parId][2])/2.
941	if parId == 142 or parId == 143:
942	values = disaggregation.darain(list(reversed(deac_vals[parId])))
943	else:
944	values = disaggregation.dapoly(list(reversed(deac_vals[parId])))
945
946	codes_set(gid, 'stepRange', 0)
947	truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
948	codes_set(gid, 'time', truedatetime.hour * 100)
949	codes_set(gid, 'date', truedatetime.year * 10000 +
950	truedatetime.month * 100 +
951	truedatetime.day)
952	codes_set_values(gid, values)
953	codes_write(gid, g_handle)
954
955	codes_release(gid)
956
957	gid = codes_new_from_index(iid)
958
959	f_handle.close()
960	g_handle.close()
961	h_handle.close()
962
963	codes_index_release(iid)
964
965	# disaggregation.IA3(lsp_np)
966	# pro gitterpunkt die zeitserie disagg, oder geht es auch mit feld
967	# wenn zurück dann über zeitschritte itterieren und rausschreiben.
968	# original in die flux
969	# subgrid wohin? extrafiles die beide oder jeweils einen subgrid beinhalten?
970	# kann time hh und mm fassen?
971
972	print lsp_np.shape
973	print lsp_np
974
975
976	return
977
978
979	def create(self, inputfiles, c):
980	'''An index file will be created which depends on the combination
981	of "date", "time" and "stepRange" values. This is used to iterate
982	over all messages in each grib file which were passed through the
983	parameter "inputfiles" to seperate specific parameters into fort.*
984	files. Afterwards the FORTRAN program is called to convert
985	the data fields all to the same grid and put them in one file
986	per unique time step (combination of "date", "time" and
987	"stepRange").
988
989	Note
990	----
991	This method is based on the ECMWF example index.py
992	https://software.ecmwf.int/wiki/display/GRIB/index.py
993
994	Parameters
995	----------
996	inputfiles : :obj:`UioFiles`
997	Contains a list of files.
998
999	c : :obj:`ControlFile`
1000	Contains all the parameters of CONTROL file and
1001	command line.
1002
1003	Return
1004	------
1005
1006	'''
1007
1008	if c.wrf:
1009	table128 = init128(_config.PATH_GRIBTABLE)
1010	wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
1011	stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
1012	table128)
1013
1014	# these numbers are indices for the temporary files "fort.xx"
1015	# which are used to seperate the grib fields to,
1016	# for the Fortran program input
1017	# 10: U,V \| 11: T \| 12: lnsp \| 13: D \| 16: sfc fields
1018	# 17: Q \| 18: Q , gaussian\| 19: w \| 21: etadot \| 22: clwc+ciwc
1019	fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
1020	'17':None, '18':None, '19':None, '21':None, '22':None}
1021
1022	iid = None
1023	index_vals = None
1024
1025	# get the values of the keys which are used for distinct access
1026	# of grib messages via product
1027	index_keys = ["date", "time", "step"]
1028	iid, index_vals = self._mk_index_values(c.inputdir,
1029	inputfiles,
1030	index_keys)
1031	# index_vals looks like e.g.:
1032	# index_vals[0]: ('20171106', '20171107', '20171108') ; date
1033	# index_vals[1]: ('0', '1200', '1800', '600') ; time
1034	# index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
1035
1036	# "product" genereates each possible combination between the
1037	# values of the index keys
1038	for prod in product(*index_vals):
1039	# e.g. prod = ('20170505', '0', '12')
1040	# ( date ,time, step)
1041
1042	print('current product: ', prod)
1043
1044	for i in range(len(index_keys)):
1045	codes_index_select(iid, index_keys[i], prod[i])
1046
1047	# get first id from current product
1048	gid = codes_new_from_index(iid)
1049
1050	# if there is no data for this specific time combination / product
1051	# skip the rest of the for loop and start with next timestep/product
1052	if not gid:
1053	continue
1054	#============================================================================================
1055	# remove old fort.* files and open new ones
1056	# they are just valid for a single product
1057	for k, f in fdict.iteritems():
1058	fortfile = os.path.join(c.inputdir, 'fort.' + k)
1059	silent_remove(fortfile)
1060	fdict[k] = open(fortfile, 'w')
1061	#============================================================================================
1062	# create correct timestamp from the three time informations
1063	cdate = str(codes_get(gid, 'date'))
1064	ctime = '{:0>2}'.format(codes_get(gid, 'time')/100)
1065	cstep = '{:0>3}'.format(codes_get(gid, 'step'))
1066	timestamp = datetime.strptime(cdate + ctime, '%Y%m%d%H')
1067	timestamp += timedelta(hours=int(cstep))
1068	cdate_hour = datetime.strftime(timestamp, '%Y%m%d%H')
1069
1070	# if the timestamp is out of basetime start/end date period,
1071	# skip this specific product
1072	if c.basetime:
1073	start_time = datetime.strptime(c.end_date + c.basetime,
1074	'%Y%m%d%H') - time_delta
1075	end_time = datetime.strptime(c.end_date + c.basetime,
1076	'%Y%m%d%H')
1077	if timestamp < start_time or timestamp > end_time:
1078	continue
1079
1080	if c.wrf:
1081	if 'olddate' not in locals() or cdate != olddate:
1082	fwrf = open(os.path.join(c.outputdir,
1083	'WRF' + cdate + '.' + ctime + '.000.grb2'), 'w')
1084	olddate = cdate[:]
1085	#============================================================================================
1086	# savedfields remembers which fields were already used.
1087	savedfields = []
1088	# sum of cloud liquid and ice water content
1089	scwc = None
1090	while 1:
1091	if not gid:
1092	break
1093	paramId = codes_get(gid, 'paramId')
1094	gridtype = codes_get(gid, 'gridType')
1095	levtype = codes_get(gid, 'typeOfLevel')
1096	if paramId == 77: # ETADOT
1097	codes_write(gid, fdict['21'])
1098	elif paramId == 130: # T
1099	codes_write(gid, fdict['11'])
1100	elif paramId == 131 or paramId == 132: # U, V wind component
1101	codes_write(gid, fdict['10'])
1102	elif paramId == 133 and gridtype != 'reduced_gg': # Q
1103	codes_write(gid, fdict['17'])
1104	elif paramId == 133 and gridtype == 'reduced_gg': # Q, gaussian
1105	codes_write(gid, fdict['18'])
1106	elif paramId == 135: # W
1107	codes_write(gid, fdict['19'])
1108	elif paramId == 152: # LNSP
1109	codes_write(gid, fdict['12'])
1110	elif paramId == 155 and gridtype == 'sh': # D
1111	codes_write(gid, fdict['13'])
1112	elif paramId == 246 or paramId == 247: # CLWC, CIWC
1113	# sum cloud liquid water and ice
1114	if scwc is None:
1115	scwc = codes_get_values(gid)
1116	else:
1117	scwc += codes_get_values(gid)
1118	codes_set_values(gid, scwc)
1119	codes_set(gid, 'paramId', 201031)
1120	codes_write(gid, fdict['22'])
1121	elif c.wrf and paramId in [129, 138, 155] and \
1122	levtype == 'hybrid': # Z, VO, D
1123	# do not do anything right now
1124	# these are specific parameter for WRF
1125	pass
1126	else:
1127	if paramId not in savedfields:
1128	# SD/MSL/TCC/10U/10V/2T/2D/Z/LSM/SDOR/CVL/CVH/SR
1129	# and all ADDPAR parameter
1130	codes_write(gid, fdict['16'])
1131	savedfields.append(paramId)
1132	else:
1133	print('duplicate ' + str(paramId) + ' not written')
1134
1135	try:
1136	if c.wrf:
1137	# model layer
1138	if levtype == 'hybrid' and \
1139	paramId in [129, 130, 131, 132, 133, 138, 155]:
1140	codes_write(gid, fwrf)
1141	# sfc layer
1142	elif paramId in wrfpars:
1143	codes_write(gid, fwrf)
1144	except AttributeError:
1145	pass
1146
1147	codes_release(gid)
1148	gid = codes_new_from_index(iid)
1149	#============================================================================================
1150	for f in fdict.values():
1151	f.close()
1152	#============================================================================================
1153	# call for Fortran program to convert e.g. reduced_gg grids to
1154	# regular_ll and calculate detadot/dp
1155	pwd = os.getcwd()
1156	os.chdir(c.inputdir)
1157	if os.stat('fort.21').st_size == 0 and c.eta:
1158	print('Parameter 77 (etadot) is missing, most likely it is \
1159	not available for this type or date/time\n')
1160	print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
1161	my_error(c.mailfail, 'fort.21 is empty while parameter eta \
1162	is set to 1 in CONTROL file')
1163	#============================================================================================
1164	# write out all output to log file before starting fortran programm
1165	sys.stdout.flush()
1166
1167	# Fortran program creates file fort.15 (with u,v,etadot,t,sp,q)
1168	p = subprocess.check_call([os.path.join(
1169	c.exedir, _config.FORTRAN_EXECUTABLE)], shell=True)
1170	os.chdir(pwd)
1171	#============================================================================================
1172	# create name of final output file, e.g. EN13040500 (ENYYMMDDHH)
1173	if c.maxstep > 12:
1174	suffix = cdate[2:8] + '.' + ctime + '.' + cstep
1175	else:
1176	suffix = cdate_hour[2:10]
1177	fnout = os.path.join(c.inputdir, c.prefix + suffix)
1178	print("outputfile = " + fnout)
1179	# collect for final processing
1180	self.outputfilelist.append(os.path.basename(fnout))
1181	#============================================================================================
1182	# create outputfile and copy all data from intermediate files
1183	# to the outputfile (final GRIB input files for FLEXPART)
1184	orolsm = os.path.basename(glob.glob(c.inputdir +
1185	'/OG_OROLSM__SL..' + c.ppid + '')[0])
1186	fluxfile = 'flux' + cdate[0:2] + suffix
1187	if not c.cwc:
1188	flist = ['fort.15', fluxfile, 'fort.16', orolsm]
1189	else:
1190	flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]
1191
1192	with open(fnout, 'wb') as fout:
1193	for f in flist:
1194	shutil.copyfileobj(open(os.path.join(c.inputdir, f), 'rb'),
1195	fout)
1196
1197	if c.omega:
1198	with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout:
1199	shutil.copyfileobj(open(os.path.join(c.inputdir, 'fort.25'),
1200	'rb'), fout)
1201	#============================================================================================
1202	if c.wrf:
1203	fwrf.close()
1204
1205	codes_index_release(iid)
1206
1207	return
1208
1209
1210	def process_output(self, c):
1211	'''Postprocessing of FLEXPART input files.
1212
1213	The grib files are postprocessed depending on the selection in
1214	CONTROL file. The resulting files are moved to the output
1215	directory if its not equal to the input directory.
1216	The following modifications might be done if
1217	properly switched in CONTROL file:
1218	GRIB2 - Conversion to GRIB2
1219	ECTRANS - Transfer of files to gateway server
1220	ECSTORAGE - Storage at ECMWF server
1221
1222	Parameters
1223	----------
1224	c : :obj:`ControlFile`
1225	Contains all the parameters of CONTROL file and
1226	command line.
1227
1228	Return
1229	------
1230
1231	'''
1232
1233	print('\n\nPostprocessing:\n Format: {}\n'.format(c.format))
1234
1235	if not c.ecapi:
1236	print('ecstorage: {}\n ecfsdir: {}\n'.
1237	format(c.ecstorage, c.ecfsdir))
1238	print('ectrans: {}\n gateway: {}\n destination: {}\n '
1239	.format(c.ectrans, c.gateway, c.destination))
1240
1241	print('Output filelist: ')
1242	print(self.outputfilelist)
1243
1244	for ofile in self.outputfilelist:
1245	ofile = os.path.join(self.inputdir, ofile)
1246
1247	if c.format.lower() == 'grib2':
1248	p = subprocess.check_call(['grib_set', '-s', 'edition=2,',
1249	'productDefinitionTemplateNumber=8',
1250	ofile, ofile + '_2'])
1251	p = subprocess.check_call(['mv', ofile + '_2', ofile])
1252
1253	if c.ectrans and not c.ecapi:
1254	p = subprocess.check_call(['ectrans', '-overwrite', '-gateway',
1255	c.gateway, '-remote', c.destination,
1256	'-source', ofile])
1257
1258	if c.ecstorage and not c.ecapi:
1259	p = subprocess.check_call(['ecp', '-o', ofile,
1260	os.path.expandvars(c.ecfsdir)])
1261
1262	if c.outputdir != c.inputdir:
1263	p = subprocess.check_call(['mv',
1264	os.path.join(c.inputdir, ofile),
1265	c.outputdir])
1266
1267	return
1268
1269
1270	def prepare_fp_files(self, c):
1271	'''Conversion of GRIB files to FLEXPART binary format.
1272
1273	Parameters
1274	----------
1275	c : :obj:`ControlFile`
1276	Contains all the parameters of CONTROL file and
1277	command line.
1278
1279	Return
1280	------
1281
1282	'''
1283	# generate AVAILABLE file
1284	# Example of AVAILABLE file data:
1285	# 20131107 000000 EN13110700 ON DISC
1286	clist = []
1287	for ofile in self.outputfilelist:
1288	fname = ofile.split('/')
1289	if '.' in fname[-1]:
1290	l = fname[-1].split('.')
1291	timestamp = datetime.strptime(l[0][-6:] + l[1],
1292	'%y%m%d%H')
1293	timestamp += timedelta(hours=int(l[2]))
1294	cdate = datetime.strftime(timestamp, '%Y%m%d')
1295	chms = datetime.strftime(timestamp, '%H%M%S')
1296	else:
1297	cdate = '20' + fname[-1][-8:-2]
1298	chms = fname[-1][-2:] + '0000'
1299	clist.append(cdate + ' ' + chms + ' '*6 +
1300	fname[-1] + ' '*14 + 'ON DISC')
1301	clist.sort()
1302	with open(c.outputdir + '/' + 'AVAILABLE', 'w') as f:
1303	f.write('\n'.join(clist) + '\n')
1304
1305	# generate pathnames file
1306	pwd = os.path.abspath(c.outputdir)
1307	with open(pwd + '/pathnames', 'w') as f:
1308	f.write(pwd + '/Options/\n')
1309	f.write(pwd + '/\n')
1310	f.write(pwd + '/\n')
1311	f.write(pwd + '/AVAILABLE\n')
1312	f.write(' = == = == = == = == = == == = \n')
1313
1314	# create Options dir if necessary
1315	if not os.path.exists(pwd + '/Options'):
1316	make_dir(pwd+'/Options')
1317
1318	# read template COMMAND file
1319	with open(os.path.expandvars(os.path.expanduser(
1320	c.flexpart_root_scripts)) + '/../Options/COMMAND', 'r') as f:
1321	lflist = f.read().split('\n')
1322
1323	# find index of list where to put in the
1324	# date and time information
1325	# usually after the LDIRECT parameter
1326	i = 0
1327	for l in lflist:
1328	if 'LDIRECT' in l.upper():
1329	break
1330	i += 1
1331
1332	# insert the date and time information of run start and end
1333	# into the list of lines of COMMAND file
1334	lflist = lflist[:i+1] + \
1335	[clist[0][:16], clist[-1][:16]] + \
1336	lflist[i+3:]
1337
1338	# write the new COMMAND file
1339	with open(pwd + '/Options/COMMAND', 'w') as g:
1340	g.write('\n'.join(lflist) + '\n')
1341
1342	# change to outputdir and start the grib2flexpart run
1343	# afterwards switch back to the working dir
1344	os.chdir(c.outputdir)
1345	p = subprocess.check_call([
1346	os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))
1347	+ '/../FLEXPART_PROGRAM/grib2flexpart', 'useAvailable', '.'])
1348	os.chdir(pwd)
1349
1350	return

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

Download in other formats: