Context Navigation

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

ctbtodev

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

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

Download in other formats: