Context Navigation

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

ctbtodev

Last change on this file since aa275fc was aa275fc, checked in by Anne Philipp <anne.philipp@…>, 6 years ago
substituted grib_api with eccodes
Property mode set to `100644`
File size: 54.1 KB

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

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

Download in other formats: