Context Navigation

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

ctbtodev

Last change on this file since e1228f3 was e1228f3, checked in by Anne Philipp <anne.philipp@…>, 6 years ago
resolved loop import between controlfile and tools
Property mode set to `100644`
File size: 56.3 KB

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

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

Download in other formats: