Context Navigation

source: flex_extract.git/Source/Python/Classes/EcFlexpart.py @ a676cf7

dev

Last change on this file since a676cf7 was a676cf7, checked in by Anne Tipka <anne.tipka@…>, 18 months ago
modifications and adjustments to dates and times for correctly retrieving analysis and forecast data in operational mode
Property mode set to `100644`
File size: 80.2 KB

Line
1	#!/usr/bin/env python3
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	#
27	# February 2018 - Anne Philipp (University of Vienna):
28	# - applied PEP8 style guide
29	# - added documentation
30	# - removed function getFlexpartTime in class EcFlexpart
31	# - outsourced class ControlFile
32	# - outsourced class MarsRetrieval
33	# - changed class name from EIFlexpart to EcFlexpart
34	# - applied minor code changes (style)
35	# - removed "dead code" , e.g. retrieval of Q since it is not needed
36	# - removed "times" parameter from retrieve-method since it is not used
37	# - seperated function "retrieve" into smaller functions (less code
38	# duplication, easier testing)
39	#
40	# @License:
41	# (C) Copyright 2014-2020.
42	# Anne Philipp, Leopold Haimberger
43	#
44	# SPDX-License-Identifier: CC-BY-4.0
45	#
46	# This work is licensed under the Creative Commons Attribution 4.0
47	# International License. To view a copy of this license, visit
48	# http://creativecommons.org/licenses/by/4.0/ or send a letter to
49	# Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
50	#*******************************************************************************
51	#pylint: disable=unsupported-assignment-operation
52	# this is disabled because for this specific case its an error in pylint
53	#pylint: disable=consider-using-enumerate
54	# this is not useful in this case
55	#pylint: disable=unsubscriptable-object
56	# this error is a bug
57	#pylint: disable=ungrouped-imports
58	# not necessary that we group the imports
59	# ------------------------------------------------------------------------------
60	# MODULES
61	# ------------------------------------------------------------------------------
62	from __future__ import print_function
63
64	import os
65	import sys
66	import glob
67	import shutil
68	from datetime import datetime, timedelta
69
70	# software specific classes and modules from flex_extract
71	#pylint: disable=wrong-import-position
72	sys.path.append('../')
73	import _config
74	from Classes.GribUtil import GribUtil
75	from Mods.tools import (init128, to_param_id, silent_remove, product,
76	my_error, get_informations, get_dimensions,
77	execute_subprocess, to_param_id_with_tablenumber,
78	generate_retrieval_period_boundary)
79	from Classes.MarsRetrieval import MarsRetrieval
80	from Classes.UioFiles import UioFiles
81	import Mods.disaggregation as disaggregation
82	#pylint: enable=wrong-import-position
83	# ------------------------------------------------------------------------------
84	# CLASS
85	# ------------------------------------------------------------------------------
86	class EcFlexpart(object):
87	'''
88	Class to represent FLEXPART specific ECMWF data.
89
90	FLEXPART needs grib files in a specifc format. All necessary data fields
91	for one time step are stored in a single file. The class represents an
92	instance with all the parameter and settings necessary for retrieving
93	MARS data and modifing them so they are fitting FLEXPART needs. The class
94	is able to disaggregate the fluxes and convert grid types to the one needed
95	by FLEXPART, therefore using the FORTRAN program.
96
97	Attributes
98	----------
99	mreq_count : int
100	Counter for the number of generated mars requests.
101
102	inputdir : str
103	Path to the directory where the retrieved data is stored.
104
105	dataset : str
106	For public datasets there is the specific naming and parameter
107	dataset which has to be used to characterize the type of
108	data.
109
110	basetime : int
111	The time for a half day retrieval. The 12 hours upfront are to be
112	retrieved.
113
114	dtime : str
115	Time step in hours.
116
117	acctype : str
118	The field type for the accumulated forecast fields.
119
120	acctime : str
121	The starting time from the accumulated forecasts.
122
123	accmaxstep : str
124	The maximum forecast step for the accumulated forecast fields.
125
126	marsclass : str
127	Characterisation of dataset.
128
129	stream : str
130	Identifies the forecasting system used to generate the data.
131
132	number : str
133	Selects the member in ensemble forecast run.
134
135	resol : str
136	Specifies the desired triangular truncation of retrieved data,
137	before carrying out any other selected post-processing.
138
139	accuracy : str
140	Specifies the number of bits per value to be used in the
141	generated GRIB coded fields.
142
143	addpar : str
144	List of additional parameters to be retrieved.
145
146	level : str
147	Specifies the maximum level.
148
149	expver : str
150	The version of the dataset.
151
152	levelist : str
153	Specifies the required levels.
154
155	glevelist : str
156	Specifies the required levels for gaussian grids.
157
158	gaussian : str
159	This parameter is deprecated and should no longer be used.
160	Specifies the desired type of Gaussian grid for the output.
161
162	grid : str
163	Specifies the output grid which can be either a Gaussian grid
164	or a Latitude/Longitude grid.
165
166	area : str
167	Specifies the desired sub-area of data to be extracted.
168
169	purefc : int
170	Switch for definition of pure forecast mode or not.
171
172	outputfilelist : list of str
173	The final list of FLEXPART ready input files.
174
175	types : dictionary
176	Determines the combination of type of fields, time and forecast step
177	to be retrieved.
178
179	params : dictionary
180	Collection of grid types and their corresponding parameters,
181	levels, level types and the grid definition.
182
183	server : ECMWFService or ECMWFDataServer
184	This is the connection to the ECMWF data servers.
185
186	public : int
187	Decides which Web API Server version is used.
188
189	dates : str
190	Contains start and end date of the retrieval in the format
191	"YYYYMMDD/to/YYYYMMDD"
192	'''
193
194	# --------------------------------------------------------------------------
195	# CLASS FUNCTIONS
196	# --------------------------------------------------------------------------
197	def __init__(self, c, fluxes=False):
198	'''Creates an object/instance of EcFlexpart with the associated
199	settings of its attributes for the retrieval.
200
201	Parameters:
202	-----------
203	c : ControlFile
204	Contains all the parameters of CONTROL file and
205	command line.
206
207	fluxes : boolean, optional
208	Decides if the flux parameter settings are stored or
209	the rest of the parameter list.
210	Default value is False.
211
212	Return
213	------
214
215	'''
216	# set a counter for the number of generated mars requests
217	self.mreq_count = 0
218
219	self.inputdir = c.inputdir
220	self.dataset = c.dataset
221	self.basetime = c.basetime
222	self.dtime = c.dtime
223	self.acctype = c.acctype
224	self.acctime = c.acctime
225	self.accmaxstep = c.accmaxstep
226	self.marsclass = c.marsclass
227	self.stream = c.stream
228	self.number = c.number
229	self.resol = c.resol
230	self.accuracy = c.accuracy
231	self.addpar = c.addpar
232	self.level = c.level
233	self.expver = c.expver
234	self.levelist = c.levelist
235	self.glevelist = '1/to/' + c.level # in case of gaussian grid
236	self.gaussian = c.gaussian
237	self.grid = c.grid
238	self.area = c.area
239	self.purefc = c.purefc
240	self.outputfilelist = []
241
242	# Define the different types of field combinations (type, time, step)
243	self.types = {}
244	# Define the parameters and their level types, level list and
245	# grid resolution for the retrieval job
246	self.params = {}
247
248	if fluxes:
249	self._create_params_fluxes()
250	else:
251	self._create_params(c.gauss, c.eta, c.omega, c.cwc, c.wrf)
252
253	if fluxes:# and not c.purefc:
254	self._create_field_types_fluxes()
255	else:
256	self._create_field_types(c.type, c.time, c.step)
257	return
258
259	def _create_field_types(self, ftype, ftime, fstep):
260	'''Create the combination of field type, time and forecast step.
261
262	Parameters:
263	-----------
264	ftype : list of str
265	List of field types.
266
267	ftime : list of str
268	The time in hours of the field.
269
270	fstep : str
271	Specifies the forecast time step from forecast base time.
272	Valid values are hours (HH) from forecast base time.
273
274	Return
275	------
276
277	'''
278	i = 0
279	for ty, st, ti in zip(ftype, fstep, ftime):
280	btlist = list(range(len(ftime)))
281	if self.basetime == 12:
282	btlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
283	if self.basetime == 0:
284	btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0]
285
286	# if ((ty.upper() == 'AN' and (int(c.time[i]) % int(c.dtime)) == 0) or
287	# (ty.upper() != 'AN' and (int(c.step[i]) % int(c.dtime)) == 0 and
288	# (int(c.step[i]) % int(c.dtime) == 0)) ) and \
289	# (int(c.time[i]) in btlist or c.purefc):
290
291	if (i in btlist) or self.purefc:
292
293	if ((ty.upper() == 'AN' and (int(ti) % int(self.dtime)) == 0) or
294	(ty.upper() != 'AN' and (int(st) % int(self.dtime)) == 0)):
295
296	if ty not in self.types.keys():
297	self.types[ty] = {'times': '', 'steps': ''}
298
299	if ti not in self.types[ty]['times']:
300	if self.types[ty]['times']:
301	self.types[ty]['times'] += '/'
302	self.types[ty]['times'] += ti
303
304	if st not in self.types[ty]['steps']:
305	if self.types[ty]['steps']:
306	self.types[ty]['steps'] += '/'
307	self.types[ty]['steps'] += st
308	i += 1
309
310	return
311
312	def _create_field_types_fluxes(self):
313	'''Create the combination of field type, time and forecast step
314	for the flux data.
315
316	Parameters:
317	-----------
318
319	Return
320	------
321
322	'''
323	if self.purefc:
324	# need to retrieve forecasts for step 000 in case of pure forecast
325	steps = '{}/to/{}/by/{}'.format(0, self.accmaxstep, self.dtime)
326	else:
327	steps = '{}/to/{}/by/{}'.format(self.dtime,
328	self.accmaxstep,
329	self.dtime)
330
331	self.types[str(self.acctype)] = {'times': str(self.acctime),
332	'steps': steps}
333	return
334
335	def _create_params(self, gauss, eta, omega, cwc, wrf):
336	'''Define the specific parameter settings for retrievment.
337
338	The different parameters need specific grid types and level types
339	for retrievement. We might get following combination of types
340	(depending on selection and availability):
341	(These are short cuts for the grib file names (leading sequence)
342	SH__ML, OG__ML, GG__ML, SH__SL, OG__SL, GG__SL, OG_OROLSM_SL
343	where:
344	SH = Spherical Harmonics, GG = Gaussian Grid, OG = Output Grid,
345	ML = Model Level, SL = Surface Level
346
347	For each of this combination there is a list of parameter names,
348	the level type, the level list and the grid resolution.
349
350	There are different scenarios for data extraction from MARS:
351	1) Retrieval of etadot
352	eta=1, gauss=0, omega=0
353	2) Calculation of etadot from divergence
354	eta=0, gauss=1, omega=0
355	3) Calculation of etadot from omega (for makes sense for debugging)
356	eta=0, gauss=0, omega=1
357	4) Retrieval and Calculation of etadot (only for debugging)
358	eta=1, gauss=1, omega=0
359	5) Download also specific model and surface level data for FLEXPART-WRF
360
361	Parameters:
362	-----------
363	gauss : int
364	Gaussian grid is retrieved.
365
366	eta : int
367	Etadot parameter will be directly retrieved.
368
369	omega : int
370	The omega paramterwill be retrieved.
371
372	cwc : int
373	The cloud liquid and ice water content will be retrieved.
374
375	wrf : int
376	Additional model level and surface level data will be retrieved for
377	WRF/FLEXPART-WRF simulations.
378
379	Return
380	------
381
382	'''
383	# SURFACE FIELDS
384	#-----------------------------------------------------------------------
385	self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF']
386	self.params['OG__SL'] = ['SD/MSL/TCC/10U/10V/2T/2D/Z/LSM', \
387	'SFC', '1', self.grid]
388	if self.addpar:
389	self.params['OG__SL'][0] += self.addpar
390
391	if self.marsclass.upper() == 'EA' or self.marsclass.upper() == 'EP':
392	self.params['OG_OROLSM__SL'] = ["SDOR/CVL/CVH/FSR",
393	'SFC', '1', self.grid]
394	else:
395	self.params['OG_OROLSM__SL'] = ["SDOR/CVL/CVH/SR", \
396	'SFC', '1', self.grid]
397
398	# MODEL LEVEL FIELDS
399	#-----------------------------------------------------------------------
400	self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]
401
402	if not gauss and eta:
403	self.params['OG__ML'][0] += '/U/V/ETADOT'
404	elif gauss and not eta:
405	self.params['GG__SL'] = ['Q', 'ML', '1',
406	'{}'.format((int(self.resol) + 1) // 2)]
407	self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF']
408	elif not gauss and not eta:
409	self.params['OG__ML'][0] += '/U/V'
410	else: # GAUSS and ETA
411	print('Warning: Collecting etadot and parameters for gaussian grid '
412	'is a very costly parameter combination, '
413	'use this combination only for debugging!')
414	self.params['GG__SL'] = ['Q', 'ML', '1',
415	'{}'.format((int(self.resol) + 1) // 2)]
416	self.params['GG__ML'] = ['U/V/D/ETADOT', 'ML', self.glevelist,
417	'{}'.format((int(self.resol) + 1) // 2)]
418
419	if omega:
420	self.params['OG__ML'][0] += '/W'
421
422	if cwc:
423	self.params['OG__ML'][0] += '/CLWC/CIWC'
424
425	# ADDITIONAL FIELDS FOR FLEXPART-WRF MODEL (IF QUESTIONED)
426	# ----------------------------------------------------------------------
427	if wrf:
428	# @WRF
429	# THIS IS NOT YET CORRECTLY IMPLEMENTED !!!
430	#
431	# UNDER CONSTRUCTION !!!
432	#
433
434	print('WRF VERSION IS UNDER CONSTRUCTION!') # dummy argument
435
436	#self.params['OG__ML'][0] += '/Z/VO'
437	#if '/D' not in self.params['OG__ML'][0]:
438	# self.params['OG__ML'][0] += '/D'
439
440	#wrf_sfc = ['SP','SKT','SST','CI','STL1','STL2', 'STL3','STL4',
441	# 'SWVL1','SWVL2','SWVL3','SWVL4']
442	#for par in wrf_sfc:
443	# if par not in self.params['OG__SL'][0]:
444	# self.params['OG__SL'][0] += '/' + par
445
446	return
447
448
449	def _create_params_fluxes(self):
450	'''Define the parameter setting for flux data.
451
452	Flux data are accumulated fields in time and are stored on the
453	surface level. The leading short cut name for the grib files is:
454	"OG_acc_SL" with OG for Regular Output Grid, SL for Surface Level, and
455	acc for Accumulated Grid.
456	The params dictionary stores a list of parameter names, the level type,
457	the level list and the grid resolution.
458
459	The flux data are: LSP/CP/SSHF/EWSS/NSSS/SSR
460
461	Parameters:
462	-----------
463
464	Return
465	------
466
467	'''
468	self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR",
469	'SFC', '1', self.grid]
470	return
471
472
473	def _mk_targetname(self, ftype, param, date):
474	'''Creates the filename for the requested grib data to be stored in.
475	This name is passed as the "target" parameter in the request.
476
477	Parameters
478	----------
479	ftype : str
480	Shortcut name of the type of the field. E.g. AN, FC, PF, ...
481
482	param : str
483	Shortcut of the grid type. E.g. SH__ML, SH__SL, GG__ML,
484	GG__SL, OG__ML, OG__SL, OG_OROLSM_SL, OG_acc_SL
485
486	date : str
487	The date period of the grib data to be stored in this file.
488
489	Return
490	------
491	targetname : str
492	The target filename for the grib data.
493	'''
494	targetname = (self.inputdir + '/' + ftype + param + '.' + date + '.' +
495	str(os.getppid()) + '.' + str(os.getpid()) + '.grb')
496
497	return targetname
498
499
500	def _start_retrievement(self, request, par_dict):
501	'''Creates the Mars Retrieval and prints or submits the request
502	depending on the status of the request variable.
503
504	Parameters
505	----------
506	request : int
507	Selects the mode of retrieval.
508	0: Retrieves the data from ECMWF.
509	1: Prints the mars requests to an output file.
510	2: Retrieves the data and prints the mars request.
511
512	par_dict : dictionary
513	Contains all parameter which have to be set for creating the
514	Mars Retrievals. The parameter are:
515	marsclass, dataset, stream, type, levtype, levelist, resol,
516	gaussian, accuracy, grid, target, area, date, time, number,
517	step, expver, param
518
519	Return
520	------
521
522	'''
523	# increase number of mars requests
524	self.mreq_count += 1
525
526	MR = MarsRetrieval(self.server,
527	self.public,
528	marsclass=par_dict['marsclass'],
529	dataset=par_dict['dataset'],
530	stream=par_dict['stream'],
531	type=par_dict['type'],
532	levtype=par_dict['levtype'],
533	levelist=par_dict['levelist'],
534	resol=par_dict['resol'],
535	gaussian=par_dict['gaussian'],
536	accuracy=par_dict['accuracy'],
537	grid=par_dict['grid'],
538	target=par_dict['target'],
539	area=par_dict['area'],
540	date=par_dict['date'],
541	time=par_dict['time'],
542	number=par_dict['number'],
543	step=par_dict['step'],
544	expver=par_dict['expver'],
545	param=par_dict['param'])
546
547	if request == 0:
548	MR.display_info()
549	MR.data_retrieve()
550	elif request == 1:
551	MR.print_infodata_csv(self.inputdir, self.mreq_count)
552	elif request == 2:
553	MR.print_infodata_csv(self.inputdir, self.mreq_count)
554	MR.display_info()
555	MR.data_retrieve()
556	else:
557	print('Failure')
558
559	return
560
561
562	def _mk_index_values(self, inputdir, inputfiles, keys):
563	'''Creates an index file for a set of grib parameter keys.
564	The values from the index keys are returned in a list.
565
566	Parameters
567	----------
568	keys : dictionary
569	List of parameter names which serves as index.
570
571	inputfiles : UioFiles
572	Contains a list of files.
573
574	Return
575	------
576	iid : codes_index
577	This is a grib specific index structure to access
578	messages in a file.
579
580	index_vals : list of list of str
581	Contains the values from the keys used for a distinct selection
582	of grib messages in processing the grib files.
583	Content looks like e.g.:
584	index_vals[0]: ('20171106', '20171107', '20171108') ; date
585	index_vals[1]: ('0', '1200', '1800', '600') ; time
586	index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
587	'''
588	from eccodes import codes_index_get
589
590	iid = None
591	index_keys = keys
592
593	indexfile = os.path.join(inputdir, _config.FILE_GRIB_INDEX)
594	silent_remove(indexfile)
595	grib = GribUtil(inputfiles.files)
596	# creates new index file
597	iid = grib.index(index_keys=index_keys, index_file=indexfile)
598
599	# read the values of index keys
600	index_vals = []
601	for key in index_keys:
602	key_vals = codes_index_get(iid, key)
603	# have to sort the key values for correct order,
604	# therefore convert to int first
605	key_vals = [int(k) for k in key_vals]
606	key_vals.sort()
607	key_vals = [str(k) for k in key_vals]
608	index_vals.append(key_vals)
609	# index_vals looks for example like:
610	# index_vals[0]: ('20171106', '20171107', '20171108') ; date
611	# index_vals[1]: ('0', '1200') ; time
612	# index_vals[2]: (3', '6', '9', '12') ; stepRange
613
614	return iid, index_vals
615
616
617	def retrieve(self, server, dates, public, purefc, request, inputdir='.'):
618	'''Finalizing the retrieval information by setting final details
619	depending on grid type.
620	Prepares MARS retrievals per grid type and submits them.
621
622	Parameters
623	----------
624	server : ECMWFService or ECMWFDataServer
625	The connection to the ECMWF server. This is different
626	for member state users which have full access and non
627	member state users which have only access to the public
628	data sets. The decision is made from command line argument
629	"public"; for public access its True (ECMWFDataServer)
630	for member state users its False (ECMWFService)
631
632	dates : str
633	Contains start and end date of the retrieval in the format
634	"YYYYMMDD/to/YYYYMMDD"
635
636	public : int
637	Switch to select kind of ECMWF Web API access and the possible data sets.
638	Public data sets (1) and Member states data sets (0).
639
640	purefc : int
641	Switch to decide whether the job is a pure forecast retrieval or
642	coupled with analysis data.
643
644	request : int
645	Selects the mode of retrieval.
646	0: Retrieves the data from ECMWF.
647	1: Prints the mars requests to an output file.
648	2: Retrieves the data and prints the mars request.
649
650	inputdir : str, optional
651	Path to the directory where the retrieved data is about
652	to be stored. The default is the current directory ('.').
653
654	Return
655	------
656
657	'''
658	self.dates = dates
659	self.server = server
660	self.public = public
661	self.inputdir = inputdir
662	self.purefc = purefc
663	oro = False
664
665	# define times with datetime module
666	t12h = timedelta(hours=12)
667	t24h = timedelta(hours=24)
668
669	# dictionary which contains all parameter for the mars request,
670	# entries with a "None" will change in different requests and will
671	# therefore be set in each request seperately
672	retr_param_dict = {'marsclass':self.marsclass,
673	'dataset':self.dataset,
674	'stream':None,
675	'type':None,
676	'levtype':None,
677	'levelist':None,
678	'resol':self.resol,
679	'gaussian':None,
680	'accuracy':self.accuracy,
681	'grid':None,
682	'target':None,
683	'area':None,
684	'date':None,
685	'time':None,
686	'number':self.number,
687	'step':None,
688	'expver':self.expver,
689	'param':None}
690
691	for ftype in sorted(self.types):
692	# ftype contains field types such as
693	# [AN, FC, PF, CV]
694	for pk, pv in sorted(self.params.items()):
695	# pk contains one of these keys of params
696	# [SH__ML, SH__SL, GG__ML, GG__SL, OG__ML, OG__SL,
697	# OG_OROLSM_SL, OG_acc_SL]
698	# pv contains all of the items of the belonging key
699	# [param, levtype, levelist, grid]
700	if isinstance(pv, str):
701	continue
702	retr_param_dict['type'] = ftype
703	retr_param_dict['time'] = self.types[ftype]['times']
704	retr_param_dict['step'] = self.types[ftype]['steps']
705	retr_param_dict['date'] = self.dates
706	retr_param_dict['stream'] = self.stream
707	retr_param_dict['target'] = \
708	self._mk_targetname(ftype,
709	pk,
710	retr_param_dict['date'].split('/')[0])
711	table128 = init128(_config.PATH_GRIBTABLE)
712	ids = to_param_id_with_tablenumber(pv[0], table128)
713	retr_param_dict['param'] = ids
714	retr_param_dict['levtype'] = pv[1]
715	retr_param_dict['levelist'] = pv[2]
716	retr_param_dict['grid'] = pv[3]
717	retr_param_dict['area'] = self.area
718	retr_param_dict['gaussian'] = self.gaussian
719
720	if pk == 'OG_OROLSM__SL' and not oro:
721	oro = True
722	# in CERA20C (class EP) there is no stream "OPER"!
723	if self.marsclass.upper() != 'EP':
724	retr_param_dict['stream'] = 'OPER'
725	retr_param_dict['type'] = 'AN'
726	retr_param_dict['time'] = '00'
727	retr_param_dict['step'] = '000'
728	retr_param_dict['date'] = self.dates.split('/')[0]
729	retr_param_dict['target'] = self._mk_targetname('',
730	pk,
731	retr_param_dict['date'])
732	elif pk == 'OG_OROLSM__SL' and oro:
733	continue
734	if pk == 'GG__SL' and pv[0] == 'Q':
735	retr_param_dict['area'] = ""
736	retr_param_dict['gaussian'] = 'reduced'
737	if ftype.upper() == 'FC' and \
738	'acc' not in retr_param_dict['target']:
739	if (int(retr_param_dict['time'][0]) +
740	int(retr_param_dict['step'][0])) > 23:
741	dates = retr_param_dict['date'].split('/')
742	sdate = datetime.strptime(dates[0], '%Y%m%d%H')
743	sdate = sdate - timedelta(days=1)
744	retr_param_dict['date'] = '/'.join(
745	[sdate.strftime("%Y%m%d")] +
746	retr_param_dict['date'][1:])
747
748	print('CHANGED FC start date to ' +
749	sdate.strftime("%Y%m%d") +
750	' to accommodate TIME=' +
751	retr_param_dict['time'][0] +
752	', STEP=' +
753	retr_param_dict['time'][0])
754
755	# ------ on demand path --------------------------------------------------
756	if self.basetime is None:
757	# ******* start retrievement
758	self._start_retrievement(request, retr_param_dict)
759	# ------ operational path ------------------------------------------------
760	else:
761	# check if mars job requests fields beyond basetime.
762	# if yes eliminate those fields since they may not
763	# be accessible with user's credentials
764	dates = retr_param_dict['date'].split('/')
765	enddate = retr_param_dict['date'].split('/')[-1]
766	elimit = datetime.strptime(enddate, '%Y%m%d')
767
768	if self.basetime == 12:
769	# -------------- flux data ----------------------------
770	if 'acc' in pk and not self.purefc:
771
772	retr_param_dict['date'] = dates[0]
773	retr_param_dict['time'] = '12'
774	retr_param_dict['target'] = \
775	self._mk_targetname(ftype, pk,
776	retr_param_dict['date'])
777
778	# ******* start retrievement
779	self._start_retrievement(request, retr_param_dict)
780
781	retr_param_dict['date'] = dates[-1]
782	retr_param_dict['time'] = '00/12'
783	retr_param_dict['target'] = \
784	self._mk_targetname(ftype, pk,
785	retr_param_dict['date'])
786
787	# ******* start retrievement
788	self._start_retrievement(request, retr_param_dict)
789
790	# -------------- non flux data ------------------------
791	else:
792	# ******* start retrievement
793	self._start_retrievement(request, retr_param_dict)
794
795	elif self.basetime == 0:
796	# -------------- flux data ----------------------------
797	if 'acc' in pk and not self.purefc:
798
799	retr_param_dict['date'] = dates[0]
800	retr_param_dict['time'] = '00/12'
801	retr_param_dict['target'] = \
802	self._mk_targetname(ftype, pk,
803	retr_param_dict['date'])
804
805	# ******* start retrievement
806	self._start_retrievement(request, retr_param_dict)
807
808	retr_param_dict['date'] = dates[-1]
809	retr_param_dict['time'] = '00'
810	retr_param_dict['target'] = \
811	self._mk_targetname(ftype, pk,
812	retr_param_dict['date'])
813
814	# ******* start retrievement
815	self._start_retrievement(request, retr_param_dict)
816	elif 'acc' in pk and self.purefc:
817	# ******* start retrievement
818	self._start_retrievement(request, retr_param_dict)
819	# -------------- non flux data ------------------------
820	else: # 'acc' not in pk
821	timesave = ''.join(retr_param_dict['time'])
822
823	if all(['/' in retr_param_dict['time'],
824	pk != 'OG_OROLSM__SL']):
825	times = retr_param_dict['time'].split('/')
826	steps = retr_param_dict['step'].split('/')
827
828	while int(times[0]) + int(steps[0]) <= 12:
829	times = times[1:]
830	if len(times) > 1:
831	retr_param_dict['time'] = '/'.join(times)
832	else:
833	retr_param_dict['time'] = times[0]
834
835	if all([pk != 'OG_OROLSM__SL',
836	int(retr_param_dict['step'].split('/')[0]) == 0,
837	int(timesave.split('/')[0]) == 0]):
838
839	retr_param_dict['date'] = \
840	datetime.strftime(elimit, '%Y%m%d')
841	retr_param_dict['time'] = '00'
842	retr_param_dict['step'] = '000'
843	retr_param_dict['target'] = \
844	self._mk_targetname(ftype, pk,
845	retr_param_dict['date'])
846
847	if ftype.upper() == 'FC' and not self.purefc:
848
849	retr_param_dict['date'] = \
850	datetime.strftime(elimit - t24h, '%Y%m%d')
851
852	# ******* start retrievement
853	self._start_retrievement(request, retr_param_dict)
854	else:
855	raise ValueError('ERROR: Basetime has an invalid value '
856	'-> {}'.format(str(self.basetime)))
857
858	if request == 0 or request == 2:
859	print('MARS retrieve done ... ')
860	elif request == 1:
861	print('MARS request printed ...')
862
863	return
864
865
866	def write_namelist(self, c):
867	'''Creates a namelist file in the temporary directory and writes
868	the following values to it: maxl, maxb, mlevel,
869	mlevelist, mnauf, metapar, rlo0, rlo1, rla0, rla1,
870	momega, momegadiff, mgauss, msmooth, meta, metadiff, mdpdeta
871
872	Parameters
873	----------
874	c : ControlFile
875	Contains all the parameters of CONTROL file and
876	command line.
877
878	filename : str
879	Name of the namelist file.
880
881	Return
882	------
883
884	'''
885
886	from genshi.template.text import NewTextTemplate
887	from genshi.template import TemplateLoader
888	from genshi.template.eval import UndefinedError
889	import numpy as np
890
891	try:
892	loader = TemplateLoader(_config.PATH_TEMPLATES, auto_reload=False)
893	namelist_template = loader.load(_config.TEMPFILE_NAMELIST,
894	cls=NewTextTemplate)
895
896	self.inputdir = c.inputdir
897	area = np.asarray(self.area.split('/')).astype(float)
898	grid = np.asarray(self.grid.split('/')).astype(float)
899
900	if area[1] > area[3]:
901	area[1] -= 360
902	maxl = int(round((area[3] - area[1]) / grid[1])) + 1
903	maxb = int(round((area[0] - area[2]) / grid[0])) + 1
904
905	stream = namelist_template.generate(
906	maxl=str(maxl),
907	maxb=str(maxb),
908	mlevel=str(self.level),
909	mlevelist=str(self.levelist),
910	mnauf=str(self.resol),
911	metapar='77',
912	rlo0=str(area[1]),
913	rlo1=str(area[3]),
914	rla0=str(area[2]),
915	rla1=str(area[0]),
916	momega=str(c.omega),
917	momegadiff=str(c.omegadiff),
918	mgauss=str(c.gauss),
919	msmooth=str(c.smooth),
920	meta=str(c.eta),
921	metadiff=str(c.etadiff),
922	mdpdeta=str(c.dpdeta)
923	)
924	except UndefinedError as e:
925	print('... ERROR ' + str(e))
926
927	sys.exit('\n... error occured while trying to generate namelist ' +
928	_config.TEMPFILE_NAMELIST)
929	except OSError as e:
930	print('... ERROR CODE: ' + str(e.errno))
931	print('... ERROR MESSAGE:\n \t ' + str(e.strerror))
932
933	sys.exit('\n... error occured while trying to generate template ' +
934	_config.TEMPFILE_NAMELIST)
935
936	try:
937	namelistfile = os.path.join(self.inputdir, _config.FILE_NAMELIST)
938
939	with open(namelistfile, 'w') as f:
940	f.write(stream.render('text'))
941	except OSError as e:
942	print('... ERROR CODE: ' + str(e.errno))
943	print('... ERROR MESSAGE:\n \t ' + str(e.strerror))
944
945	sys.exit('\n... error occured while trying to write ' +
946	namelistfile)
947
948	return
949
950
951	def deacc_fluxes(self, inputfiles, c):
952	'''De-accumulate and disaggregate flux data.
953
954	Goes through all flux fields in ordered time and de-accumulate
955	the fields. Afterwards the fields are disaggregated in time.
956	Different versions of disaggregation is provided for rainfall
957	data (darain, modified linear) and the surface fluxes and
958	stress data (dapoly, cubic polynomial).
959
960	Parameters
961	----------
962	inputfiles : UioFiles
963	Contains the list of files that contain flux data.
964
965	c : ControlFile
966	Contains all the parameters of CONTROL file and
967	command line.
968
969	Return
970	------
971
972	'''
973	import numpy as np
974	from eccodes import (codes_index_select, codes_get,
975	codes_get_values, codes_set_values, codes_set,
976	codes_write, codes_release, codes_new_from_index,
977	codes_index_release)
978
979	table128 = init128(_config.PATH_GRIBTABLE)
980	# get ids from the flux parameter names
981	pars = to_param_id(self.params['OG_acc_SL'][0], table128)
982
983	iid = None
984	index_vals = None
985
986	# get the values of the keys which are used for distinct access
987	# of grib messages via product and save the maximum number of
988	# ensemble members if there is more than one
989	if '/' in self.number:
990	# more than one ensemble member is selected
991	index_keys = ["number", "date", "time", "step"]
992	# maximum ensemble number retrieved
993	# + 1 for the control run (ensemble number 0)
994	maxnum = int(self.number.split('/')[-1]) + 1
995	# remember the index of the number values
996	index_number = index_keys.index('number')
997	# empty set to save ensemble numbers which were already processed
998	ens_numbers = set()
999	# index for the ensemble number
1000	inumb = 0
1001	else:
1002	index_keys = ["date", "time", "step"]
1003	# maximum ensemble number
1004	maxnum = None
1005
1006	# get sorted lists of the index values
1007	# this is very important for disaggregating
1008	# the flux data in correct order
1009	iid, index_vals = self._mk_index_values(c.inputdir,
1010	inputfiles,
1011	index_keys)
1012	# index_vals looks like e.g.:
1013	# index_vals[0]: ('20171106', '20171107', '20171108') ; date
1014	# index_vals[1]: ('0', '600', '1200', '1800') ; time
1015	# index_vals[2]: ('0', '3', '6', '9', '12') ; stepRange
1016
1017	if c.rrint:
1018	# set start and end timestamps for retrieval period
1019	if not c.purefc:
1020	start_date = datetime.strptime(c.start_date + '00', '%Y%m%d%H')
1021	end_date = datetime.strptime(c.end_date + '23', '%Y%m%d%H')
1022	else:
1023	sdate_str = c.start_date + '{:0>2}'.format(index_vals[1][0])
1024	start_date = datetime.strptime(sdate_str, '%Y%m%d%H')
1025	edate_str = c.end_date + '{:0>2}'.format(index_vals[1][-1])
1026	end_date = datetime.strptime(edate_str, '%Y%m%d%H')
1027	end_date = end_date + timedelta(hours=c.maxstep)
1028
1029	# get necessary grid dimensions from grib files for storing the
1030	# precipitation fields
1031	info = get_informations(os.path.join(c.inputdir,
1032	inputfiles.files[0]))
1033	dims = get_dimensions(info, c.purefc, c.dtime, index_vals,
1034	start_date, end_date)
1035
1036	# create empty numpy arrays
1037	if not maxnum:
1038	lsp_np = np.zeros((dims[1] * dims[0], dims[2]), dtype=np.float64)
1039	cp_np = np.zeros((dims[1] * dims[0], dims[2]), dtype=np.float64)
1040	else:
1041	lsp_np = np.zeros((maxnum, dims[1] * dims[0], dims[2]), dtype=np.float64)
1042	cp_np = np.zeros((maxnum, dims[1] * dims[0], dims[2]), dtype=np.float64)
1043
1044	# index counter for time line
1045	it_lsp = 0
1046	it_cp = 0
1047
1048	# store the order of date and step
1049	date_list = []
1050	step_list = []
1051
1052	# initialize dictionaries to store flux values per parameter
1053	orig_vals = {}
1054	deac_vals = {}
1055	for p in pars:
1056	orig_vals[p] = []
1057	deac_vals[p] = []
1058
1059	# "product" genereates each possible combination between the
1060	# values of the index keys
1061	for prod in product(*index_vals):
1062	# e.g. prod = ('20170505', '0', '12')
1063	# ( date ,time, step)
1064
1065	print('CURRENT PRODUCT: ', prod)
1066
1067	# the whole process has to be done for each seperate ensemble member
1068	# therefore, for each new ensemble member we delete old flux values
1069	# and start collecting flux data from the beginning time step
1070	if maxnum and prod[index_number] not in ens_numbers:
1071	ens_numbers.add(prod[index_number])
1072	inumb = len(ens_numbers) - 1
1073	# re-initialize dictionaries to store flux values per parameter
1074	# for the next ensemble member
1075	it_lsp = 0
1076	it_cp = 0
1077	orig_vals = {}
1078	deac_vals = {}
1079	for p in pars:
1080	orig_vals[p] = []
1081	deac_vals[p] = []
1082
1083	for i in range(len(index_keys)):
1084	codes_index_select(iid, index_keys[i], prod[i])
1085
1086	# get first id from current product
1087	gid = codes_new_from_index(iid)
1088
1089	# if there is no data for this specific time combination / product
1090	# skip the rest of the for loop and start with next timestep/product
1091	if not gid:
1092	continue
1093
1094	# create correct timestamp from the three time informations
1095	cdate = str(codes_get(gid, 'date'))
1096	time = codes_get(gid, 'time') // 100 # integer
1097	step = codes_get(gid, 'step') # integer
1098	ctime = '{:0>2}'.format(time)
1099
1100	t_date = datetime.strptime(cdate + ctime, '%Y%m%d%H')
1101	t_dt = t_date + timedelta(hours=step)
1102	t_m1dt = t_date + timedelta(hours=step-int(c.dtime))
1103	t_m2dt = t_date + timedelta(hours=step-2*int(c.dtime))
1104	if c.basetime is not None:
1105	t_enddate = datetime.strptime(c.end_date + str(c.basetime),
1106	'%Y%m%d%H')
1107	else:
1108	t_enddate = t_date + timedelta(2*int(c.dtime))
1109
1110	# if necessary, add ensemble member number to filename suffix
1111	# otherwise, add empty string
1112	if maxnum:
1113	numbersuffix = '.N{:0>3}'.format(int(prod[index_number]))
1114	else:
1115	numbersuffix = ''
1116
1117	if c.purefc:
1118	fnout = os.path.join(c.inputdir, 'flux' +
1119	t_date.strftime('%Y%m%d.%H') +
1120	'.{:0>3}'.format(step-2*int(c.dtime)) +
1121	numbersuffix)
1122	gnout = os.path.join(c.inputdir, 'flux' +
1123	t_date.strftime('%Y%m%d.%H') +
1124	'.{:0>3}'.format(step-int(c.dtime)) +
1125	numbersuffix)
1126	hnout = os.path.join(c.inputdir, 'flux' +
1127	t_date.strftime('%Y%m%d.%H') +
1128	'.{:0>3}'.format(step) +
1129	numbersuffix)
1130	else:
1131	fnout = os.path.join(c.inputdir, 'flux' +
1132	t_m2dt.strftime('%Y%m%d%H') + numbersuffix)
1133	gnout = os.path.join(c.inputdir, 'flux' +
1134	t_m1dt.strftime('%Y%m%d%H') + numbersuffix)
1135	hnout = os.path.join(c.inputdir, 'flux' +
1136	t_dt.strftime('%Y%m%d%H') + numbersuffix)
1137
1138	print("outputfile = " + fnout)
1139	f_handle = open(fnout, 'wb')
1140	h_handle = open(hnout, 'wb')
1141	g_handle = open(gnout, 'wb')
1142
1143	# read message for message and store relevant data fields, where
1144	# data keywords are stored in pars
1145	while True:
1146	if not gid:
1147	break
1148	parId = codes_get(gid, 'paramId') # integer
1149	step = codes_get(gid, 'step') # integer
1150	time = codes_get(gid, 'time') # integer
1151	ni = codes_get(gid, 'Ni') # integer
1152	nj = codes_get(gid, 'Nj') # integer
1153	if parId not in orig_vals.keys():
1154	# parameter is not a flux, find next one
1155	continue
1156
1157	# define conversion factor
1158	if parId == 142 or parId == 143:
1159	fak = 1. / 1000.
1160	else:
1161	fak = 3600.
1162
1163	# get parameter values and reshape
1164	values = codes_get_values(gid)
1165	values = (np.reshape(values, (nj, ni))).flatten() / fak
1166
1167	# save the original and accumulated values
1168	orig_vals[parId].append(values[:])
1169
1170	if c.marsclass.upper() == 'EA' or step <= int(c.dtime):
1171	# no de-accumulation needed
1172	deac_vals[parId].append(values[:] / int(c.dtime))
1173	else:
1174	# do de-accumulation
1175	deac_vals[parId].append(
1176	(orig_vals[parId][-1] - orig_vals[parId][-2]) /
1177	int(c.dtime))
1178
1179	# store precipitation if new disaggregation method is selected
1180	# only the exact days are needed
1181	if c.rrint:
1182	if start_date <= t_dt <= end_date:
1183	if not c.purefc:
1184	if t_dt not in date_list:
1185	date_list.append(t_dt)
1186	step_list = [0]
1187	else:
1188	if t_date not in date_list:
1189	date_list.append(t_date)
1190	if step not in step_list:
1191	step_list.append(step)
1192	# store precipitation values
1193	if maxnum and parId == 142:
1194	lsp_np[inumb, :, it_lsp] = deac_vals[parId][-1][:]
1195	it_lsp += 1
1196	elif not maxnum and parId == 142:
1197	lsp_np[:, it_lsp] = deac_vals[parId][-1][:]
1198	it_lsp += 1
1199	elif maxnum and parId == 143:
1200	cp_np[inumb, :, it_cp] = deac_vals[parId][-1][:]
1201	it_cp += 1
1202	elif not maxnum and parId == 143:
1203	cp_np[:, it_cp] = deac_vals[parId][-1][:]
1204	it_cp += 1
1205
1206	# information printout
1207	print(parId, time, step, len(values), values[0], np.std(values))
1208
1209	# length of deac_vals[parId] corresponds to the
1210	# number of time steps, max. 4 are needed for disaggegration
1211	# with the old and original method
1212	# run over all grib messages and perform
1213	# shifting in time
1214	if len(deac_vals[parId]) >= 3:
1215	if len(deac_vals[parId]) > 3:
1216	if not c.rrint and (parId == 142 or parId == 143):
1217	values = disaggregation.darain(deac_vals[parId])
1218	else:
1219	values = disaggregation.dapoly(deac_vals[parId])
1220
1221	if not (step == c.maxstep and c.purefc \
1222	or t_dt == t_enddate):
1223	# remove first time step in list to shift
1224	# time line
1225	orig_vals[parId].pop(0)
1226	deac_vals[parId].pop(0)
1227	else:
1228	# if the third time step is read (per parId),
1229	# write out the first one as a boundary value
1230	if c.purefc:
1231	values = deac_vals[parId][1]
1232	else:
1233	values = deac_vals[parId][0]
1234
1235	if not (c.rrint and (parId == 142 or parId == 143)):
1236	codes_set_values(gid, values)
1237
1238	if c.purefc:
1239	codes_set(gid, 'stepRange', max(0, step-2*int(c.dtime)))
1240	else:
1241	codes_set(gid, 'stepRange', 0)
1242	codes_set(gid, 'time', t_m2dt.hour*100)
1243	codes_set(gid, 'date', int(t_m2dt.strftime('%Y%m%d')))
1244
1245	codes_write(gid, f_handle)
1246
1247	# squeeze out information of last two steps
1248	# contained in deac_vals[parId]
1249	# Note that deac_vals[parId][0] has not been popped
1250	# in this case
1251
1252	if step == c.maxstep and c.purefc or \
1253	t_dt == t_enddate:
1254	# last step
1255	if c.purefc:
1256	values = deac_vals[parId][3]
1257	codes_set_values(gid, values)
1258	codes_set(gid, 'stepRange', step)
1259	#truedatetime = t_m2dt + timedelta(hours=2*int(c.dtime))
1260	codes_write(gid, h_handle)
1261	else:
1262	values = deac_vals[parId][3]
1263	codes_set_values(gid, values)
1264	codes_set(gid, 'stepRange', 0)
1265	truedatetime = t_m2dt + timedelta(hours=2*int(c.dtime))
1266	codes_set(gid, 'time', truedatetime.hour * 100)
1267	codes_set(gid, 'date', int(truedatetime.strftime('%Y%m%d')))
1268	codes_write(gid, h_handle)
1269
1270	if parId == 142 or parId == 143:
1271	values = disaggregation.darain(list(reversed(deac_vals[parId])))
1272	else:
1273	values = disaggregation.dapoly(list(reversed(deac_vals[parId])))
1274
1275	# step before last step
1276	if c.purefc:
1277	codes_set(gid, 'stepRange', step-int(c.dtime))
1278	#truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
1279	codes_set_values(gid, values)
1280	codes_write(gid, g_handle)
1281	else:
1282	codes_set(gid, 'stepRange', 0)
1283	truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
1284	codes_set(gid, 'time', truedatetime.hour * 100)
1285	codes_set(gid, 'date', int(truedatetime.strftime('%Y%m%d')))
1286	codes_set_values(gid, values)
1287	codes_write(gid, g_handle)
1288
1289	codes_release(gid)
1290
1291	gid = codes_new_from_index(iid)
1292
1293	f_handle.close()
1294	g_handle.close()
1295	h_handle.close()
1296
1297	codes_index_release(iid)
1298
1299	if c.rrint:
1300	self._create_rr_grib_dummy(inputfiles.files[0], c.inputdir)
1301
1302	self._prep_new_rrint(dims[0], dims[1], dims[2], lsp_np,
1303	cp_np, maxnum, index_keys, index_vals, c)
1304
1305	return
1306
1307	def _prep_new_rrint(self, ni, nj, nt, lsp_np, cp_np, maxnum, index_keys, index_vals, c):
1308	'''Calculates and writes out the disaggregated precipitation fields.
1309
1310	Disaggregation is done in time and original times are written to the
1311	flux files, while the additional subgrid times are written to
1312	extra files output files. They are named like the original files with
1313	suffixes "_1" and "_2" for the first and second subgrid point.
1314
1315	Parameters
1316	----------
1317	ni : int
1318	Amount of zonal grid points.
1319
1320	nj : int
1321	Amount of meridional grid points.
1322
1323	nt : int
1324	Number of time steps.
1325
1326	lsp_np : numpy array of float
1327	The large scale precipitation fields for each time step.
1328	Shape (ni * nj, nt).
1329
1330	cp_np : numpy array of float
1331	The convective precipitation fields for each time step.
1332	Shape (ni * nj, nt).
1333
1334	maxnum : int
1335	The maximum number of ensemble members. It is None
1336	if there are no or just one ensemble.
1337
1338	index_keys : dictionary
1339	List of parameter names which serves as index.
1340
1341	index_vals : list of list of str
1342	Contains the values from the keys used for a distinct selection
1343	of grib messages in processing the grib files.
1344	Content looks like e.g.:
1345	index_vals[0]: ('20171106', '20171107', '20171108') ; date
1346	index_vals[1]: ('0', '1200', '1800', '600') ; time
1347	index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
1348
1349	c : ControlFile
1350	Contains all the parameters of CONTROL file and
1351	command line.
1352
1353	Return
1354	------
1355
1356	'''
1357	import numpy as np
1358
1359	print('... disaggregation of precipitation with new method.')
1360
1361	tmpfile = os.path.join(c.inputdir, 'rr_grib_dummy.grb')
1362
1363	# initialize new numpy arrays for disaggregated fields
1364	if maxnum:
1365	lsp_new_np = np.zeros((maxnum, ni * nj, nt * 3), dtype=np.float64)
1366	cp_new_np = np.zeros((maxnum, ni * nj, nt * 3), dtype=np.float64)
1367	else:
1368	lsp_new_np = np.zeros((1, ni * nj, nt * 3), dtype=np.float64)
1369	cp_new_np = np.zeros((1, ni * nj, nt * 3), dtype=np.float64)
1370
1371	# do the disaggregation, but neglect the last value of the
1372	# original time series. This one corresponds for example to
1373	# 24 hour, which we don't need. we use 0 - 23 UTC for a day.
1374	if maxnum:
1375	for inum in range(maxnum):
1376	for ix in range(ni*nj):
1377	lsp_new_np[inum, ix, :] = disaggregation.IA3(lsp_np[inum, ix, :])[:-1]
1378	cp_new_np[inum, ix, :] = disaggregation.IA3(cp_np[inum, ix, :])[:-1]
1379	else:
1380	for ix in range(ni*nj):
1381	lsp_new_np[0, ix, :] = disaggregation.IA3(lsp_np[ix, :])[:-1]
1382	cp_new_np[0, ix, :] = disaggregation.IA3(cp_np[ix, :])[:-1]
1383
1384	# write to grib files (full/orig times to flux file and inbetween
1385	# times with step 1 and 2, respectively)
1386	print('... write disaggregated precipitation to files.')
1387
1388	if maxnum:
1389	# remember the index of the number values
1390	index_number = index_keys.index('number')
1391	# empty set to save unique ensemble numbers which were already processed
1392	ens_numbers = set()
1393	# index for the ensemble number
1394	inumb = 0
1395	else:
1396	inumb = 0
1397
1398	# index variable of disaggregated fields
1399	it = 0
1400
1401	# "product" genereates each possible combination between the
1402	# values of the index keys
1403	for prod in product(*index_vals):
1404	# e.g. prod = ('20170505', '0', '12')
1405	# ( date ,time, step)
1406	# or prod = ('0' , '20170505', '0', '12')
1407	# (number, date ,time, step)
1408
1409	cdate = prod[index_keys.index('date')]
1410	ctime = '{:0>2}'.format(int(prod[index_keys.index('time')])//100)
1411	cstep = '{:0>3}'.format(int(prod[index_keys.index('step')]))
1412
1413	date = datetime.strptime(cdate + ctime, '%Y%m%d%H')
1414	date += timedelta(hours=int(cstep))
1415
1416	start_period, end_period = generate_retrieval_period_boundary(c)
1417	# skip all temporary times
1418	# which are outside the retrieval period
1419	if date < start_period or \
1420	date > end_period:
1421	continue
1422
1423	# the whole process has to be done for each seperate ensemble member
1424	# therefore, for each new ensemble member we delete old flux values
1425	# and start collecting flux data from the beginning time step
1426	if maxnum and prod[index_number] not in ens_numbers:
1427	ens_numbers.add(prod[index_number])
1428	inumb = int(prod[index_number])
1429	it = 0
1430
1431	# if necessary, add ensemble member number to filename suffix
1432	# otherwise, add empty string
1433	if maxnum:
1434	numbersuffix = '.N{:0>3}'.format(int(prod[index_number]))
1435	else:
1436	numbersuffix = ''
1437
1438	# per original time stamp: write original time step and
1439	# the two newly generated sub time steps
1440	if c.purefc:
1441	fluxfilename = 'flux' + date.strftime('%Y%m%d.%H') + '.' + cstep
1442	else:
1443	fluxfilename = 'flux' + date.strftime('%Y%m%d%H') + numbersuffix
1444
1445	# write original time step to flux file as usual
1446	fluxfile = GribUtil(os.path.join(c.inputdir, fluxfilename))
1447	fluxfile.set_keys(tmpfile, filemode='ab',
1448	wherekeynames=['paramId'], wherekeyvalues=[142],
1449	keynames=['perturbationNumber', 'date', 'time',
1450	'stepRange', 'values'],
1451	keyvalues=[inumb, int(date.strftime('%Y%m%d')),
1452	date.hour*100, 0, lsp_new_np[inumb, :, it]]
1453	)
1454	fluxfile.set_keys(tmpfile, filemode='ab',
1455	wherekeynames=['paramId'], wherekeyvalues=[143],
1456	keynames=['perturbationNumber', 'date', 'time',
1457	'stepRange', 'values'],
1458	keyvalues=[inumb, int(date.strftime('%Y%m%d')),
1459	date.hour*100, 0, cp_new_np[inumb, :, it]]
1460	)
1461
1462	# rr for first subgrid point is identified by step = 1
1463	fluxfile.set_keys(tmpfile, filemode='ab',
1464	wherekeynames=['paramId'], wherekeyvalues=[142],
1465	keynames=['perturbationNumber', 'date', 'time',
1466	'stepRange', 'values'],
1467	keyvalues=[inumb, int(date.strftime('%Y%m%d')),
1468	date.hour*100, '1', lsp_new_np[inumb, :, it+1]]
1469	)
1470	fluxfile.set_keys(tmpfile, filemode='ab',
1471	wherekeynames=['paramId'], wherekeyvalues=[143],
1472	keynames=['perturbationNumber', 'date', 'time',
1473	'stepRange', 'values'],
1474	keyvalues=[inumb, int(date.strftime('%Y%m%d')),
1475	date.hour*100, '1', cp_new_np[inumb, :, it+1]]
1476	)
1477
1478	# rr for second subgrid point is identified by step = 2
1479	fluxfile.set_keys(tmpfile, filemode='ab',
1480	wherekeynames=['paramId'], wherekeyvalues=[142],
1481	keynames=['perturbationNumber', 'date', 'time',
1482	'stepRange', 'values'],
1483	keyvalues=[inumb, int(date.strftime('%Y%m%d')),
1484	date.hour*100, '2', lsp_new_np[inumb, :, it+2]]
1485	)
1486	fluxfile.set_keys(tmpfile, filemode='ab',
1487	wherekeynames=['paramId'], wherekeyvalues=[143],
1488	keynames=['perturbationNumber', 'date', 'time',
1489	'stepRange', 'values'],
1490	keyvalues=[inumb, int(date.strftime('%Y%m%d')),
1491	date.hour*100, '2', cp_new_np[inumb, :, it+2]]
1492	)
1493
1494	it = it + 3 # jump to next original time step in rr fields
1495	return
1496
1497	def _create_rr_grib_dummy(self, ifile, inputdir):
1498	'''Creates a grib file with a dummy message for the two precipitation
1499	types lsp and cp each.
1500
1501	Parameters
1502	----------
1503	ifile : str
1504	Filename of the input file to read the grib messages from.
1505
1506	inputdir : str, optional
1507	Path to the directory where the retrieved data is stored.
1508
1509	Return
1510	------
1511
1512	'''
1513
1514	gribfile = GribUtil(os.path.join(inputdir, 'rr_grib_dummy.grb'))
1515
1516	gribfile.copy_dummy_msg(ifile, keynames=['paramId','paramId'],
1517	keyvalues=[142,143], filemode='wb')
1518
1519	return
1520
1521	def create(self, inputfiles, c):
1522	'''An index file will be created which depends on the combination
1523	of "date", "time" and "stepRange" values. This is used to iterate
1524	over all messages in each grib file which were passed through the
1525	parameter "inputfiles" to seperate specific parameters into fort.*
1526	files. Afterwards the FORTRAN program is called to convert
1527	the data fields all to the same grid and put them in one file
1528	per unique time step (combination of "date", "time" and
1529	"stepRange").
1530
1531	Note
1532	----
1533	This method is based on the ECMWF example index.py
1534	https://software.ecmwf.int/wiki/display/GRIB/index.py
1535
1536	Parameters
1537	----------
1538	inputfiles : UioFiles
1539	Contains a list of files.
1540
1541	c : ControlFile
1542	Contains all the parameters of CONTROL file and
1543	command line.
1544
1545	Return
1546	------
1547
1548	'''
1549	from eccodes import (codes_index_select, codes_get,
1550	codes_get_values, codes_set_values, codes_set,
1551	codes_write, codes_release, codes_new_from_index,
1552	codes_index_release)
1553
1554	# generate start and end timestamp of the retrieval period
1555	start_period = datetime.strptime(c.start_date + c.time[0], '%Y%m%d%H')
1556	start_period = start_period + timedelta(hours=int(c.step[0]))
1557	end_period = datetime.strptime(c.end_date + c.time[-1], '%Y%m%d%H')
1558	end_period = end_period + timedelta(hours=int(c.step[-1]))
1559
1560	# @WRF
1561	# THIS IS NOT YET CORRECTLY IMPLEMENTED !!!
1562	#
1563	# UNDER CONSTRUCTION !!!
1564	#
1565	#if c.wrf:
1566	# table128 = init128(_config.PATH_GRIBTABLE)
1567	# wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
1568	# stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
1569	# table128)
1570
1571	# these numbers are indices for the temporary files "fort.xx"
1572	# which are used to seperate the grib fields to,
1573	# for the Fortran program input
1574	# 10: U,V \| 11: T \| 12: lnsp \| 13: D \| 16: sfc fields
1575	# 17: Q \| 18: Q, SL, GG\| 19: omega \| 21: etadot \| 22: clwc+ciwc
1576	fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
1577	'17':None, '18':None, '19':None, '21':None, '22':None}
1578
1579	iid = None
1580	index_vals = None
1581
1582	# get the values of the keys which are used for distinct access
1583	# of grib messages via product
1584	if '/' in self.number:
1585	# more than one ensemble member is selected
1586	index_keys = ["number", "date", "time", "step"]
1587	else:
1588	index_keys = ["date", "time", "step"]
1589	iid, index_vals = self._mk_index_values(c.inputdir,
1590	inputfiles,
1591	index_keys)
1592	# index_vals looks like e.g.:
1593	# index_vals[0]: ('20171106', '20171107', '20171108') ; date
1594	# index_vals[1]: ('0', '600', '1200', '1800') ; time
1595	# index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
1596
1597	# "product" genereates each possible combination between the
1598	# values of the index keys
1599	for prod in product(*index_vals):
1600	# e.g. prod = ('20170505', '0', '12')
1601	# ( date ,time, step)
1602
1603	print('current product: ', prod)
1604
1605	for i in range(len(index_keys)):
1606	codes_index_select(iid, index_keys[i], prod[i])
1607
1608	# get first id from current product
1609	gid = codes_new_from_index(iid)
1610
1611	# if there is no data for this specific time combination / product
1612	# skip the rest of the for loop and start with next timestep/product
1613	if not gid:
1614	continue
1615	#============================================================================================
1616	# remove old fort.* files and open new ones
1617	# they are just valid for a single product
1618	for k, f in fdict.items():
1619	fortfile = os.path.join(c.inputdir, 'fort.' + k)
1620	silent_remove(fortfile)
1621	fdict[k] = open(fortfile, 'wb')
1622	#============================================================================================
1623	# create correct timestamp from the three time informations
1624	cdate = str(codes_get(gid, 'date'))
1625	ctime = '{:0>2}'.format(codes_get(gid, 'time') // 100)
1626	cstep = '{:0>3}'.format(codes_get(gid, 'step'))
1627	timestamp = datetime.strptime(cdate + ctime, '%Y%m%d%H')
1628	timestamp += timedelta(hours=int(cstep))
1629	cdate_hour = datetime.strftime(timestamp, '%Y%m%d%H')
1630
1631	# if basetime is used, adapt start/end date period
1632	if c.basetime is not None:
1633	time_delta = timedelta(hours=12-int(c.dtime))
1634	start_period = datetime.strptime(c.end_date + str(c.basetime),
1635	'%Y%m%d%H') - time_delta
1636	end_period = datetime.strptime(c.end_date + str(c.basetime),
1637	'%Y%m%d%H')
1638
1639	# skip all temporary times
1640	# which are outside the retrieval period
1641	if timestamp < start_period or \
1642	timestamp > end_period:
1643	continue
1644
1645
1646	# @WRF
1647	# THIS IS NOT YET CORRECTLY IMPLEMENTED !!!
1648	#
1649	# UNDER CONSTRUCTION !!!
1650	#
1651	#if c.wrf:
1652	# if 'olddate' not in locals() or cdate != olddate:
1653	# fwrf = open(os.path.join(c.outputdir,
1654	# 'WRF' + cdate + '.' + ctime + '.000.grb2'), 'wb')
1655	# olddate = cdate[:]
1656	#============================================================================================
1657	# savedfields remembers which fields were already used.
1658	savedfields = []
1659	# sum of cloud liquid and ice water content
1660	scwc = None
1661	while 1:
1662	if not gid:
1663	break
1664	paramId = codes_get(gid, 'paramId')
1665	gridtype = codes_get(gid, 'gridType')
1666	if paramId == 77: # ETADOT
1667	codes_write(gid, fdict['21'])
1668	elif paramId == 130: # T
1669	codes_write(gid, fdict['11'])
1670	elif paramId == 131 or paramId == 132: # U, V wind component
1671	codes_write(gid, fdict['10'])
1672	elif paramId == 133 and gridtype != 'reduced_gg': # Q
1673	codes_write(gid, fdict['17'])
1674	elif paramId == 133 and gridtype == 'reduced_gg': # Q, gaussian
1675	codes_write(gid, fdict['18'])
1676	elif paramId == 135: # W
1677	codes_write(gid, fdict['19'])
1678	elif paramId == 152: # LNSP
1679	codes_write(gid, fdict['12'])
1680	elif paramId == 155 and gridtype == 'sh': # D
1681	codes_write(gid, fdict['13'])
1682	elif paramId == 246 or paramId == 247: # CLWC, CIWC
1683	# sum cloud liquid water and ice
1684	if scwc is None:
1685	scwc = codes_get_values(gid)
1686	else:
1687	scwc += codes_get_values(gid)
1688	codes_set_values(gid, scwc)
1689	codes_set(gid, 'paramId', 201031)
1690	codes_write(gid, fdict['22'])
1691	scwc = None
1692	# @WRF
1693	# THIS IS NOT YET CORRECTLY IMPLEMENTED !!!
1694	#
1695	# UNDER CONSTRUCTION !!!
1696	#
1697	#elif c.wrf and paramId in [129, 138, 155] and \
1698	# levtype == 'hybrid': # Z, VO, D
1699	# # do not do anything right now
1700	# # these are specific parameter for WRF
1701	# pass
1702	else:
1703	if paramId not in savedfields:
1704	# SD/MSL/TCC/10U/10V/2T/2D/Z/LSM/SDOR/CVL/CVH/SR
1705	# and all ADDPAR parameter
1706	codes_write(gid, fdict['16'])
1707	savedfields.append(paramId)
1708	else:
1709	print('duplicate ' + str(paramId) + ' not written')
1710	# @WRF
1711	# THIS IS NOT YET CORRECTLY IMPLEMENTED !!!
1712	#
1713	# UNDER CONSTRUCTION !!!
1714	#
1715	#try:
1716	# if c.wrf:
1717	# # model layer
1718	# if levtype == 'hybrid' and \
1719	# paramId in [129, 130, 131, 132, 133, 138, 155]:
1720	# codes_write(gid, fwrf)
1721	# # sfc layer
1722	# elif paramId in wrfpars:
1723	# codes_write(gid, fwrf)
1724	#except AttributeError:
1725	# pass
1726
1727	codes_release(gid)
1728	gid = codes_new_from_index(iid)
1729	#============================================================================================
1730	for f in fdict.values():
1731	f.close()
1732	#============================================================================================
1733	# call for Fortran program to convert e.g. reduced_gg grids to
1734	# regular_ll and calculate detadot/dp
1735	pwd = os.getcwd()
1736	os.chdir(c.inputdir)
1737	if os.stat('fort.21').st_size == 0 and c.eta:
1738	print('Parameter 77 (etadot) is missing, most likely it is '
1739	'not available for this type or date / time\n')
1740	print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
1741	my_error('fort.21 is empty while parameter eta '
1742	'is set to 1 in CONTROL file')
1743	# ============================================================================================
1744	# write out all output to log file before starting fortran programm
1745	sys.stdout.flush()
1746
1747	# Fortran program creates file fort.15 (with u,v,etadot,t,sp,q)
1748	execute_subprocess([os.path.join(c.exedir,
1749	_config.FORTRAN_EXECUTABLE)],
1750	error_msg='FORTRAN PROGRAM FAILED!')#shell=True)
1751
1752	os.chdir(pwd)
1753	# ============================================================================================
1754	# create name of final output file, e.g. EN13040500 (ENYYMMDDHH)
1755	# for CERA-20C we need all 4 digits for the year sinc 1900 - 2010
1756	if c.purefc:
1757	if c.marsclass == 'EP':
1758	suffix = cdate[0:8] + '.' + ctime + '.' + cstep
1759	else:
1760	suffix = cdate[2:8] + '.' + ctime + '.' + cstep
1761	else:
1762	if c.marsclass == 'EP':
1763	suffix = cdate_hour[0:10]
1764	else:
1765	suffix = cdate_hour[2:10]
1766
1767	# if necessary, add ensemble member number to filename suffix
1768	if 'number' in index_keys:
1769	index_number = index_keys.index('number')
1770	if len(index_vals[index_number]) > 1:
1771	suffix = suffix + '.N{:0>3}'.format(int(prod[index_number]))
1772
1773	fnout = os.path.join(c.inputdir, c.prefix + suffix)
1774	print("outputfile = " + fnout)
1775	# collect for final processing
1776	self.outputfilelist.append(os.path.basename(fnout))
1777	# # get additional precipitation subgrid data if available
1778	# if c.rrint:
1779	# self.outputfilelist.append(os.path.basename(fnout + '_1'))
1780	# self.outputfilelist.append(os.path.basename(fnout + '_2'))
1781	# ============================================================================================
1782	# create outputfile and copy all data from intermediate files
1783	# to the outputfile (final GRIB input files for FLEXPART)
1784	orolsm = os.path.basename(glob.glob(c.inputdir +
1785	'/OG_OROLSM__SL.*.' +
1786	c.ppid +
1787	'*')[0])
1788	if c.marsclass == 'EP':
1789	fluxfile = 'flux' + suffix
1790	else:
1791	fluxfile = 'flux' + cdate[0:2] + suffix
1792	if not c.cwc:
1793	flist = ['fort.15', fluxfile, 'fort.16', orolsm]
1794	else:
1795	flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]
1796
1797	with open(fnout, 'wb') as fout:
1798	for f in flist:
1799	shutil.copyfileobj(open(os.path.join(c.inputdir, f), 'rb'),
1800	fout)
1801
1802	if c.omega:
1803	with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout:
1804	shutil.copyfileobj(open(os.path.join(c.inputdir, 'fort.25'),
1805	'rb'), fout)
1806	# ============================================================================================
1807
1808	# @WRF
1809	# THIS IS NOT YET CORRECTLY IMPLEMENTED !!!
1810	#
1811	# UNDER CONSTRUCTION !!!
1812	#
1813	#if c.wrf:
1814	# fwrf.close()
1815
1816	codes_index_release(iid)
1817
1818	return
1819
1820
1821	def calc_extra_elda(self, path, prefix):
1822	''' Calculates extra ensemble members for ELDA - Stream.
1823
1824	This is a specific feature which doubles the number of ensemble members
1825	for the ELDA Stream.
1826
1827	Parameters
1828	----------
1829	path : str
1830	Path to the output files.
1831
1832	prefix : str
1833	The prefix of the output filenames as defined in Control file.
1834
1835	Return
1836	------
1837
1838	'''
1839	from eccodes import (codes_grib_new_from_file, codes_get_array,
1840	codes_set_array, codes_release,
1841	codes_set, codes_write)
1842
1843	# max number
1844	maxnum = int(self.number.split('/')[-1])
1845
1846	# get a list of all prepared output files with control forecast (CF)
1847	cf_filelist = UioFiles(path, prefix + '*.N000')
1848	cf_filelist.files = sorted(cf_filelist.files)
1849
1850	for cffile in cf_filelist.files:
1851	with open(cffile, 'rb') as f:
1852	cfvalues = []
1853	while True:
1854	fid = codes_grib_new_from_file(f)
1855	if fid is None:
1856	break
1857	cfvalues.append(codes_get_array(fid, 'values'))
1858	codes_release(fid)
1859
1860	filename = cffile.split('N000')[0]
1861	for i in range(1, maxnum + 1):
1862	# read an ensemble member
1863	g = open(filename + 'N{:0>3}'.format(i), 'rb')
1864	# create file for newly calculated ensemble member
1865	h = open(filename + 'N{:0>3}'.format(i+maxnum), 'wb')
1866	# number of message in grib file
1867	j = 0
1868	while True:
1869	gid = codes_grib_new_from_file(g)
1870	if gid is None:
1871	break
1872	values = codes_get_array(gid, 'values')
1873	# generate a new ensemble member by subtracting
1874	# 2 * ( current time step value - last time step value )
1875	codes_set_array(gid, 'values',
1876	values-2*(values-cfvalues[j]))
1877	codes_set(gid, 'number', i+maxnum)
1878	codes_write(gid, h)
1879	codes_release(gid)
1880	j += 1
1881
1882	g.close()
1883	h.close()
1884	print('wrote ' + filename + 'N{:0>3}'.format(i+maxnum))
1885	self.outputfilelist.append(
1886	os.path.basename(filename + 'N{:0>3}'.format(i+maxnum)))
1887
1888	return
1889
1890
1891	def process_output(self, c):
1892	'''Postprocessing of FLEXPART input files.
1893
1894	The grib files are postprocessed depending on the selection in
1895	CONTROL file. The resulting files are moved to the output
1896	directory if its not equal to the input directory.
1897	The following modifications might be done if
1898	properly switched in CONTROL file:
1899	GRIB2 - Conversion to GRIB2
1900	ECTRANS - Transfer of files to gateway server
1901	ECSTORAGE - Storage at ECMWF server
1902
1903	Parameters
1904	----------
1905	c : ControlFile
1906	Contains all the parameters of CONTROL file and
1907	command line.
1908
1909	Return
1910	------
1911
1912	'''
1913
1914	print('\n\nPostprocessing:\n Format: {}\n'.format(c.format))
1915
1916	print('\n\nGrib compression type:\n packingType: {}\n'.format(c.compression))
1917
1918	if _config.FLAG_ON_ECMWFSERVER:
1919	print('ecstorage: {}\n ecfsdir: {}\n'.
1920	format(c.ecstorage, c.ecfsdir))
1921	print('ectrans: {}\n gateway: {}\n destination: {}\n '
1922	.format(c.ectrans, c.gateway, c.destination))
1923
1924	print('Output filelist: ')
1925	print(sorted(self.outputfilelist))
1926
1927	for ofile in self.outputfilelist:
1928	ofile = os.path.join(self.inputdir, ofile)
1929
1930	if c.format.lower() == 'grib2':
1931	execute_subprocess(['grib_set', '-s', 'edition=2,' +
1932	'productDefinitionTemplateNumber=8',
1933	ofile, ofile + '_2'],
1934	error_msg='GRIB2 CONVERSION FAILED!')
1935
1936	execute_subprocess(['mv', ofile + '_2', ofile],
1937	error_msg='RENAMING FOR NEW GRIB2 FORMAT '
1938	'FILES FAILED!')
1939
1940	if c.compression.lower() != 'grid_simple':
1941	execute_subprocess(['grib_set', '-r', '-s',
1942	'packingType=' + c.compression,
1943	ofile, ofile + '_2'],
1944	error_msg='GRIB COMPRESSION FAILED!')
1945
1946	execute_subprocess(['mv', ofile + '_2', ofile],
1947	error_msg='RENAMING FOR NEW GRIB COMPRESSION '
1948	'FILES FAILED!')
1949
1950	if c.ectrans and _config.FLAG_ON_ECMWFSERVER:
1951	execute_subprocess(['ectrans', '-overwrite', '-gateway',
1952	c.gateway, '-remote', c.destination,
1953	'-source', ofile],
1954	error_msg='TRANSFER TO LOCAL SERVER FAILED!')
1955
1956	if c.ecstorage and _config.FLAG_ON_ECMWFSERVER:
1957	execute_subprocess(['ecp', '-o', ofile,
1958	os.path.expandvars(c.ecfsdir)],
1959	error_msg='COPY OF FILES TO ECSTORAGE '
1960	'AREA FAILED!')
1961
1962	if c.outputdir != c.inputdir:
1963	execute_subprocess(['mv', os.path.join(c.inputdir, ofile),
1964	c.outputdir],
1965	error_msg='RELOCATION OF OUTPUT FILES '
1966	'TO OUTPUTDIR FAILED!')
1967
1968	return

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

Download in other formats: